text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ ae22b97a-9a72-11ec-17d5-5be07c4373a9
begin
	using Plots
	using Distributions
	using PlutoUI
end

# ╔═╡ a9658aa9-9f2b-4242-866a-22fb587e8a13
md"""
# Gamma distribution
$(@bind a Slider(1:1000))
$(@bind t Slider(1:1000))
``\alpha \qquad \qquad\qquad\qquad\Theta``
"""

# ╔═╡ 2ee86f44-144a-451f-801a-85c8d0e05ae1
begin
	alpha = a * 0.01
	theta = t * 0.01
end

# ╔═╡ 7ffbaaf8-78d1-4061-99b6-90c055142091
g = Gamma(alpha, theta)

# ╔═╡ bc0166c7-8b3a-4616-8c33-6d9e18d88c65
begin
	x = 0:0.1:10
	plot(x,pdf.(g,x))
end

# ╔═╡ 17e131ad-8382-4b42-ad6c-1896027bad60
md"""
mean:   $(mean(g))
std:   $(std(g))
"""

# ╔═╡ Cell order:
# ╠═ae22b97a-9a72-11ec-17d5-5be07c4373a9
# ╟─a9658aa9-9f2b-4242-866a-22fb587e8a13
# ╟─2ee86f44-144a-451f-801a-85c8d0e05ae1
# ╟─7ffbaaf8-78d1-4061-99b6-90c055142091
# ╟─bc0166c7-8b3a-4616-8c33-6d9e18d88c65
# ╠═17e131ad-8382-4b42-ad6c-1896027bad60


text:
using ApproxFun, SingularIntegralEquations, Test
    import ApproxFun: testbandedoperator, testraggedbelowoperator, testblockbandedoperator

include("runtests.jl")

@testset "Full" begin
    ## Memory CurveTests

    d = exp(im*Interval(0.1,0.2))
    x = Fun(d)
    w = 1/(sqrt(abs(first(d)-x))*sqrt(abs(last(d)-x)))
    testbandedoperator(Hilbert(space(w)))

    ## 3 domain ideal fluid flow

    Γ = Segment(-im,1.0-im) ∪ Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6))) ∪ Circle(2.0,0.2)
        z = Fun(Γ)

    S = PiecewiseSpace(map(d->isa(d,Circle) ? Fourier(d) : JacobiWeight(0.5,0.5,Ultraspherical(1,d)),components(Γ)))
    H = Hilbert(S)
    testblockbandedoperator(H)

    #  TODO: fix testraggedbelowoperator(H)

    Ai = [Operator(Fun(ones(component(Γ,1)),Γ)) Fun(ones(component(Γ,2)),Γ) Fun(ones(component(Γ,3)),Γ) real(H)]

    @test ApproxFun.israggedbelow(Ai)
    @test ApproxFun.israggedbelow(Ai.ops[4])
    @test ApproxFun.israggedbelow(Ai.ops[4].op)

    B=ApproxFun.SpaceOperator(ApproxFun.BasisFunctional(3),S,ApproxFun.ConstantSpace(Float64))

    Ai=[Operator(0)                 0                 0                 B;
        Fun(ones(component(Γ,1)),Γ) Fun(ones(component(Γ,2)),Γ) Fun(ones(component(Γ,3)),Γ) real(H)]

    @time testblockbandedoperator(Ai)

    k=114;
        α=exp(k/50*im)
        @time a,b,c,ui= Ai \ [0; imag(α*z)]

    u =(x,y)->α*(x+im*y)+2cauchy(ui,x+im*y)
    @test u(1.1,0.2) ≈ (-0.8290718508107162+0.511097153754im)
end

println("Example Tests")
include("ExamplesTest.jl")

include("WienerHopfTest.jl")

include("IdealFluidFlowTest.jl")

ext: jl


text:
######################################################################
# pf.jl
#
# starting to migrate over to ParticleFilters.jl
######################################################################

# generate new target position
#ParticleFilters.generate_s(sm::Vehicle, s, a, rng::AbstractRNG) = s

# To shorten some of these names
const LVR = LowVarianceResampler
const MerTwist = MersenneTwister

struct Model{V <: Vehicle, S <: Sensor, M <: MotionModel}
    x::V
    sensor::S
    motion_model::M
end
Model(x::Vehicle) = Model(x, x.sensor)
Model(x::Vehicle, s::Sensor) = Model(x, s, NoMotion())

function ParticleFilters.generate_s(m::Model, s, a, rng::AbstractRNG)
    return move_target(m.motion_model, s, 0.0)
end
# recall that "a" here is the state of the drone
ParticleFilters.obs_weight(m::Model, a::Pose, sp, o) = O(m.sensor, sp, a, o)

# TODO: should this type of thing go in model subtype?
export CircularRange
struct CircularRange
    h_min::Float64      # minimum heading, in degrees
    h_max::Float64      # maximum heading, in degrees
    v_max::Float64      # maximum velocity
end
function Base.rand(c::CircularRange)
    heading = rand()*(c.h_max - c.h_min) + c.h_min
    velocity = rand()*c.v_max

    vx = velocity * sind(heading)
    vy = velocity * cosd(heading)

    return vx, vy
end

mutable struct PF{M <: Model, OL} <: AbstractFilter
    model::M            # contains vehicle, sensor, target motion models
    n::Int              # number of particles
    obs_list::OL        # list of possible observations
    L::Float64          # length of square domain
    vr::CircularRange   # range of possible velocities

    # User should never see this
    _pf::SimpleParticleFilter{TargetTuple, M, LVR, MerTwist}
    _b::ParticleCollection{TargetTuple}       # particle set
end

function PF(model::Model, n::Int, obs_list=0:0;
            L::Real=100.0,
            vr::CircularRange=CircularRange(0,360,3)
           )
    pf = SimpleParticleFilter(model, LVR(n), Base.GLOBAL_RNG)
    b = initialize_particles(n, L, vr)

    return PF(model, n, obs_list, L, vr, pf, b)
end

ParticleFilters.state_type(m::Model) = TargetTuple

# TODO: allow user to dictate original spread of velocities
function initialize_particles(n::Int, L::Real, vr::CircularRange)
    bv = TargetTuple[]
    for i = 1:n
        # it was done like the below first
        #push!(bv, (L*rand(), L*rand(), 2.0*rand(), 2.0*rand()))
        vx, vy = rand(vr)
        push!(bv, (L*rand(), L*rand(), vx, vy) )
    end
    return ParticleCollection(bv)
end

# convenience functions
ParticleFilters.particles(pf::PF) = pf._b.particles
ParticleFilters.particle(pf::PF, i::Int) = particle(pf._b, i)
ParticleFilters.weight(pf::PF, i::Int) = weight(pf._b, i)
ParticleFilters.weight_sum(pf::PF) = weight_sum(pf._b)

function predict(pf::PF)
    m = pf.model.motion_model
    return collect( move_target(m , s, 0.0) for s in particles(pf) )
end
function predict(pf::PF, b)
    m = pf.model.motion_model
    return collect( move_target(m , s, 0.0) for s in b.particles )
end

function predict!(pf::PF)
    m = pf.model.motion_model
    particle_vec = collect( move_target(m , s, 0.0) for s in particles(pf) )
    pf._b = ParticleCollection(particle_vec)
end

# vehicle is included in model
# normally update(filter, belief, a, o), but action not needed
# so I pass in the UAV pose instead
function update!(pf::PF, p::Pose, o)
    pf._b = update(pf._pf, pf._b, p, o)
end

function update_b(pf::PF, b::ParticleCollection, p, o)
    return update(pf._pf, b, p, o)
end
export update_b

function update_bp{S}(up::SimpleParticleFilter{S}, b::ParticleCollection, a, o)
    ps = particles(b)
    pm = up._particle_memory
    wm = up._weight_memory
    resize!(pm, 0)
    resize!(wm, 0)
    sizehint!(pm, n_particles(b))
    sizehint!(wm, n_particles(b))
    for i in 1:n_particles(b)
        s = ps[i]
        push!(pm, s)
        push!(wm, obs_weight(up.model, s, a, s, o))
    end

    return resample(up.resample, WeightedParticleBelief{S}(pm, wm, sum(wm), nothing), up.rng)
end
export update_bp

# TODO: should I pass in something other than Base.GLOBAL_RNG) ?
function sample_state(pf::PF)
    return rand(Base.GLOBAL_RNG, pf._b)
end

# n is the number of states you want to sample
function sample_states(pf::PF, n::Int)
    N = n_particles(pf._b)
    bv = Array{TargetTuple}(n)
    for i = 1:n
        bv[i] = particle(pf, rand(1:N))
    end
    return ParticleCollection(bv)
end

function reset!(f::PF)
    f._b = initialize_particles(f.n, f.L, f.vr)
end

function centroid(pf::PF)
    mx = 0.0
    my = 0.0
    for i = 1:pf.n
        tp = particle(pf, i)
        mx += tp[1]
        my += tp[2]
    end

    return mx/pf.n, my/pf.n
end
function centroid(b::ParticleCollection)
    mx = 0.0
    my = 0.0
    n = length(b.particles)
    for i = 1:n
        tp = particle(b, i)
        mx += tp[1]
        my += tp[2]
    end

    return mx/n, my/n
end

function covariance(pf::PF)
    covariance(pf._b)
end
function covariance(b::ParticleCollection)
    xx = 0.0
    xy = 0.0
    yy = 0.0

    mx = 0.0
    my = 0.0
    n = n_particles(b)
    for i = 1:n
        tp = particle(b, i)

        xx += tp[1] * tp[1]
        xy += tp[1] * tp[2]
        yy += tp[2] * tp[2]

        mx += tp[1]
        my += tp[2]
    end
    xx /= n
    xy /= n
    yy /= n

    mx /= n
    my /= n

    a = xx - mx * mx
    bc = xy - mx * my
    d = yy - my * my

    return [a bc; bc d]
end

export bin_pf
function bin_pf(b::ParticleCollection, discrete_b::Matrix{Float64}, L::Float64)

    # determine the size of each cell
    n_cells = size(discrete_b, 1)        # assuming symmetric matrix
    cell_size = L / n_cells

    # ensure the matrix is filled with zeros
    fill!(discrete_b, 0.0)

    # loop over particles and determine where they belong
    n = n_particles(b)
    w = 1.0 / n
    for i = 1:n
        tp = particle(b, i)

        xi = round(Int, tp[1] / cell_size, RoundUp)
        yi = round(Int, tp[2] / cell_size, RoundUp)

        # ensure that cells are legal
        xi = min(max(xi, 1), n_cells)
        yi = min(max(yi, 1), n_cells)

        discrete_b[xi, yi] += w
    end
end

function cheap_entropy(b::ParticleCollection,m::Matrix{Float64},L::Float64)

    # bin the particle collection into a matrix
    bin_pf(b, m, L)

    # now compute entropy from the matrix
    ent = 0.0
    for p in m
        if p > 0.0
            ent -= p * log(p)
        end
    end

    return ent
end

export cheap_entropy

function cheap_entropy(pf::PF, L::Float64, n_cells::Int)
    cheap_entropy(pf._b, L, n_cells)
end
function cheap_entropy(b::ParticleCollection, L::Float64, n_cells::Int)
    discrete_b = zeros(n_cells, n_cells)
    return cheap_entropy(b, discrete_b, L)
end
export cheap_entropy2
function cheap_entropy2(pf::PF, L::Float64, n_cells::Int)
    cheap_entropy2(pf._b, L, n_cells::Int)
end
function cheap_entropy2(b::ParticleCollection, L::Float64, n_cells::Int)
    d = Dict{Tuple{Int64,Int64}, Float64}()
    cell_size = L / n_cells

    # loop over particles and determine where they belong
    n = n_particles(b)
    w = 1.0 / n
    for i = 1:n
        tp = particle(b, i)

        xi = round(Int, tp[1] / cell_size, RoundUp)
        yi = round(Int, tp[2] / cell_size, RoundUp)

        if !haskey(d, (xi,yi))
            d[(xi,yi)] = 0.0
        end
        d[(xi,yi)] += w
    end

    # loop over entries in dictionary
    ent = 0.0
    for k in keys(d)
        p = d[k]
        ent -= p * log(p)
    end
    return ent
end


text:
using Test, YaoBlocks, YaoArrayRegister

@testset "test constructor" for T in [Float16, Float32, Float64]
    # NOTE: type should follow the axis
    @test RotationGate(X, 0.1) isa PrimitiveBlock{1}
    @test_throws TypeError RotationGate{1, Complex{T}, XGate} # will not accept non-real type

    @test Rx(T(0.1)) isa RotationGate{1, T, XGate}
    @test Ry(T(0.1)) isa RotationGate{1, T, YGate}
    @test Rz(T(0.1)) isa RotationGate{1, T, ZGate}
end

@testset "test matrix" begin
    theta = 2.0
    for (DIRECTION, MAT) in [
        (X, [cos(theta/2) -im*sin(theta/2); -im*sin(theta/2) cos(theta/2)]),
        (Y, [cos(theta/2) -sin(theta/2); sin(theta/2) cos(theta/2)]),
        (Z, [exp(-im*theta/2) 0;0 exp(im*theta/2)]),
        (CNOT, exp(-mat(CNOT)/2*theta*im |> Matrix)),
        (control(2, (1,), 2=>X), exp(-mat(CNOT)/2*theta*im |> Matrix))]

        @test mat(RotationGate(DIRECTION, theta)) ≈ MAT
    end
end

@testset "test apply" begin
    r = rand_state(1)
    @test state(apply!(copy(r), Rx(0.1))) ≈ mat(Rx(0.1)) * state(r)
end

@testset "test dispatch" begin
    @test dispatch!(Rx(0.1), 0.3) == Rx(0.3)
    @test nparameters(Rx(0.1)) == 1

    @testset "test $op" for op in [+, -, *, /]
        @test dispatch!(op, Rx(0.1), π) == Rx(op(0.1, π))
    end

    @test_throws AssertionError dispatch!(Rx(0.1), (0.2, 0.3))
end

@testset "adjoints" begin
    @test Rx(0.1)' == Rx(-0.1)
    @test Rx(0.2)' == Rx(-0.2)
    @test copy(Rx(0.1)) == Rx(0.1)

    g = Rx(0.1) # creates a new one
    @test copy(g) !== g
end


text:
#module TestSystem
# this module implements a mechanism for running tests, including
# subsets of tests defined by tags
module TestSystem
using Input

global const TAG_DEFAULT = "tag_default"  # tag that every function must have
include("tags.jl")  # list of tags

global const EmptyDict = Dict{Any, Any}()

export TAG_DEFAULT, TestList, add_func1!, add_func2!, add_func3!, run_testlist,
       runTestSystem, not_isapprox

"""
  This type stores all the data that describes a set of tests and their
  associated tags.

  ** Fields**

   * funcs: list of functions that contain the tests.  The order in which
           the tests run is the same as the order of insertion

   * func_tags: Vector, same length as funcs, where each element is a vector
               of tags associated with each func.

   * func_type:  Vector, same length as funcs, describing how to run each
                function. 1 = no arguments, 2 = args = (mesh, sbp, eqn, opts),
                3 = (mesh, sbp, eqn, opts), using a modified version of an
                existing input file
   * input_name: Vector, same length as funcs, containing the name of the input
                file to use with each function.  For func_type == 1 this value
                is unused, for func type == 2 it is used directly, for
                func_type == 3 it is modified and written to a new name
                according to mod_input

   * mod_input: Vector, same length as funcs, where each element is a dictionary.
               For func_type == 3 the input file specified by input_name is
               modified with these keys.  There must also be a key called
               "new_fname" that specifies name to write the modified file
               to (not including file extension.
               This modified file will be loaded before the function
               is called.

   * tag_list: collection of all known tags
"""
mutable struct TestList
  funcs::Vector{Function}  # functions containing the tests
  func_tags::Vector{Vector{String}}  # tags associated with each function
  func_type::Vector{Int}  # enum defining how to run the functino
  input_name::Vector{String} # the name of the input file for each function
  input_mod::Vector{Dict}  # dictionary of modifications to input_name
  tag_list::Set{String}  # all known tags

  function TestList()
    sizehint = 0  # initial length of vectors
    funcs = Array{Function}(sizehint)
    func_tags = Array{Vector{String}}(sizehint)
    func_type = Array{Int}(sizehint)
    input_name = Array{String}(sizehint)
    input_mod = Array{Dict}(sizehint)
    tag_list = Set{String}()
    push!(tag_list, TAG_DEFAULT)

    return new(funcs, func_tags, func_type, input_name, input_mod, tag_list)
  end  # end constructor
end  # end type definition

#------------------------------------------------------------------------------
# TestList functions

"""
  This function adds a new test function of func_type == 1 to the list
  (a function that takes no arguments).

  **Inputs**

   * test_list: the list of tests to append the function to
   * func: the function
   * tags: an array of tags associated with this function.  The user is free
          to modify this array afterwards, but not to mutate the strings within
          the array. (optional)

"""
function add_func1!(testlist::TestList, func::Function, tags::Array{String}=String[])

  check_tags(func, tags)

  # function
  push!(testlist.funcs, func)

  # tags
  tag_list_full = Array{String}(length(tags)+1)
  tag_list_full[1] = TAG_DEFAULT
  tag_list_full[2:end] = tags
  push!(testlist.func_tags, tag_list_full)

  # function type
  push!(testlist.func_type, 1)

  # input name
  push!(testlist.input_name,  "")

  # input_mod
  push!(testlist.input_mod, EmptyDict)

  for i=1:length(tags)
    push!(testlist.tag_list, tags[i])
  end

end

"""
  This function adds a new test function of `func_type == 2` to the list
  (a function that takes the arguments (mesh, sbp, eqn, opts), as
  obtained from the input file specified).

  **Inputs**

   * test_list: the list of tests to append the function to
   * func: the function
   * input_name: name of input file to be used with this function
   * tags: an array of tags associated with this function.  The user is free
          to modify this array afterwards, but not to mutate the strings within
          the array. (optional)

"""
function add_func2!(testlist::TestList, func::Function, input_name::String,
                    tags::Array{String}=String[])

  check_tags(func, tags)

  # function
  push!(testlist.funcs, func)

  # tags
  tag_list_full = Array{String}(length(tags)+1)
  tag_list_full[1] = TAG_DEFAULT
  tag_list_full[2:end] = tags
  push!(testlist.func_tags, tag_list_full)

  # function type
  push!(testlist.func_type, 2)

  # input name
  push!(testlist.input_name,  input_name)

  # input_mod
  push!(testlist.input_mod, EmptyDict)

  for i=1:length(tags)
    push!(testlist.tag_list, tags[i])
  end

end

"""
  This function adds a new test function of `func_type == 3` to the list
  (a function that takes the arguments `(mesh, sbp, eqn, opts)`, as
  obtained from modifying the input file specified).

  **Inputs**

   * test_list: the list of tests to append the function to
   * func: the function
   * input_name: name of input file to be modified for use with this function
   * mod_dict: dictionary of keys to be added (or replaced) in the input file
   * tags: an array of tags associated with this function.  The user is free
           to modify this array afterwards, but not to mutate the strings within
           the array. (optional)

"""
function add_func3!(testlist::TestList, func::Function, input_name::String,
                    mod_dict::Dict, tags::Array{String}=String[])

  check_tags(func, tags)

  # function
  push!(testlist.funcs, func)

  # tags
  tag_list_full = Array{String}(length(tags)+1)
  tag_list_full[1] = TAG_DEFAULT
  tag_list_full[2:end] = tags
  push!(testlist.func_tags, tag_list_full)

  # function type
  push!(testlist.func_type, 3)

  # input name
  push!(testlist.input_name,  input_name)

  # input_mod
  push!(testlist.input_mod, mod_dict)

  for i=1:length(tags)
    push!(testlist.tag_list, tags[i])
  end

end

"""
  This function checks that all necessary tags are present and throws
  an exception is they are not.

  Currently this only verifies that one of LengthTags is present

  **Inputs**

   * func: the function (only needed to produce error message)
   * tags: the tags for this function

  Outputs: none
"""
function check_tags(func::Function, tags::Array{String})

  n = length(intersect(tags, LengthTags))
  if n > 1
    throw(ErrorException("function $func has $n > 1 LengthTags"))
  elseif n < 1
    throw(ErrorException("function $func has $n < 1 LengthTags"))
  end

  return nothing
end

"""
  This function runs a test list.  Tests are run in the order they were
  loaded into the TestList object.  This implementation handles the case of
  several tests sets sharing an input file efficiencly, ie. if several test
  functions use the same input file and are placed in the test list
  consecutively, the input file will be loaded only once.

  A list of tags can be optionally supplied.  In this case, only the tests
  that have the specified tags will be run.  If no list is supplied, all tags
  are run.

  **Inputs**

   * testlist: a TestList loaded with functions
   * prep_func = function used with test functions of type 2 or 3.  It must
                 have signature prep_func(fname::String). fname is the
                 name of hte input file associated with the test function
   * tags: an array of tags (optional)

  **Outputs**

   * ntests: number of tests that were run.
"""
function run_testlist(testlist::TestList, prep_func::Function, tags::Vector{String}=String[TAG_DEFAULT])

  println("Running tests with tags matching = ", tags)
  tags_set = Set(tags)
  ntests = length(testlist.funcs)

  # need to declare these here so they have the right scope
  # these variables are type unstable anways, the default the type of the
  # default value doesn't matter
  mesh = 0
  sbp = 0
  eqn = 0
  opts = 0

  ftiming = open("timing.dat", "w")
  input_name_prev = ""  # holds the previous input name
  n_run_tests = 0
  for i=1:ntests
    func_i = testlist.funcs[i]
    func_tags_i = testlist.func_tags[i]
    functype_i = testlist.func_type[i]
    input_name_i = testlist.input_name[i]
    input_mod_i = testlist.input_mod[i]

    for j=1:length(func_tags_i)
      tag_j = func_tags_i[j]
      if tag_j in tags_set
        n_run_tests += 1
        println("running function ", func_i)
        println("function type = ", functype_i)
        # run this test
        t_test = @elapsed if functype_i == 1  # function with no arguments
          func_i()

        elseif functype_i == 2  # function with all 4 arguments
#          if input_name_i != input_name_prev
            println("about to run prep function ", prep_func, " on input ", input_name_i)
            mesh, sbp, eqn, opts = prep_func(input_name_i)
#          end
          func_i(mesh, sbp, eqn, opts)
          input_name_prev = input_name_i

        elseif functype_i == 3  # modify input before running
          new_fname = input_mod_i["new_fname"]*".jl"

          if input_name_prev != new_fname  # don't create a file if it was done already
            arg_dict = evalfile(joinpath(pwd(), input_name_i))

            for (key, val) in input_mod_i
              if key != "new_fname"
                arg_dict[key] = val
              end
            end

            make_input(arg_dict, new_fname)

            mesh, sbp, eqn, opts = prep_func(new_fname)
            input_name_prev = new_fname
          end  # end if file already created

          func_i(mesh, sbp, eqn, opts)
        else
          throw(ErrorException("unsupported function type"))
        end

        println(ftiming, func_i, ": ", t_test, " seconds")
        println("time for test function ", func_i, " = ", t_test)

        break  # don't run this test more than once even if it matches
                  # multiple tags

      end  # end if tag matches
    end  # end loop over tags

  end  # end loop over test functions

  close(ftiming)

  return n_run_tests
end

"""
  This function runs a test list, according to the tags supplied.  Unlike
  [`run_tests`](@ref), this function supplies default behavior for the `tags`
  list.  If the list is empty, then all tests are run, otherwise only the
  specified tag are run.  See [`run_tests`](@ref) for argument descriptions.

  **Inputs**

   * testlist
   * prep_func
   * tags
"""
function runTestSystem(testlist::TestList, prep_func::Function, tags::Vector{String})

  nargs = length(tags)
  if nargs == 0
    tags2 = String[TAG_DEFAULT]
  else
    tags2 = Array{String}(nargs)
    copy!(tags2, tags)
  end

  n = run_testlist(testlist, prep_func, tags2)

  if n == 0
    error("No tests run matching tags: ", tags2)
  end

  return nothing
end

"""
  Helper function to compute !isapprox(args..., kwargs) in @test macros because

  @test !isapprox(args..., kwargs...)

  doesn't work.
"""
function not_isapprox(args...; kwargs...)
  return !isapprox(args...; kwargs...)
end

end  # end module


text:
"""
    fdichkspec(sysf::FDIModel, SFDI::BitMatrix; sdeg, FDtol, FDGainTol, FDfreq,
                 atol, atol1, atol2, atol3, rtol, fast = true, minimal = false) -> (rdims, orders, leastorders)

Check for a given synthesis model `sysf::FDIModel` the feasibility of a set of fault detection and isolation specifications `SFDI`.
If `SFDI` has `N` rows (i.e., contains `N` specifications), then the `N`-dimensional integer vectors `rdims`, `orders`, `leastorders`
are returned and contain information related to the synthesis of FDI filters to achieve the feasible specifications.
For the `i`-th specification contained in `SFDI[i,:]`, `rdims[i]` contains the number of residual
outputs of a minimal nullspace basis based FDI filter which can be used to achieve this specification.
If `rdims[i] = 0`, then the `i`-th specification is not feasible. For a feasible `i`-th specification, `orders[i]`
contains the order of the minimal nullspace basis based FDI filter which can be used to achieve this specification.
If the `i`-th specification is not feasible, then `orders[i]` is set to `-1`.
If `minimal = true`, `leastorders[i]` contains the least achievable order for a scalar output FDI filter which can be used
to achieve the `i`-th specification. If `minimal = false` or if the `i`-th specification is not
feasible, then `leastorders[i]` is set to `-1`.

`FDFreq = freq` specifies a vector of real frequency values or a scalar real frquency value
for strong detectability checks (default: `FDFreq = missing`).

`FDtol = tol1` specifies the threshold `tol1` for assessing weak specifications
                      (see also function [`fditspec`](@ref)) (default: `tol1 = 0.0001`).

`FDGainTol = tol2` specifies the threshold `tol2` for assessing strong specifications,
i.e., the threshold for nonzero frequency responce gains for all frequency values specified in
`freq` (see also function [`fdisspec`](@ref)) (default: `tol2 = 0.01`).

The keyword argument `sdeg = β` specifies a prescribed stability degree `β` for the poles of the internally
generated candidate filters, such that the real parts of
filters poles must be less than or equal to `β`, in the continuous-time case, and
the magnitudes of filter poles must be less than or
equal to `β`, in the discrete-time case. If `sdeg = missing` then no then no stabilization is performed if and `FDFreq = missing`.
If `sdeg = missing` and `FDFreq = freq`, then the fllowing default values are employed : `β = -0.05`, in continuous-time case, and  `β = 0.95`,
in discrete-time case.

The rank determinations in the performed reductions
are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C`, `D`,
the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the order of the system `sysf.sys`.
The keyword argument `atol3` is an absolute tolerance for observability tests
(default: internally determined value).
The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`.

_Method:_ The nullspace method of [1] is successively employed to
determine FDI filters as minimal left nullspace bases which solve
suitably formulated fault detection problems.

_References:_

[1] A. Varga, On computing nullspace bases – a fault detection perspective.
      Proc. IFAC 2008 World Congress, Seoul, Korea, pages 6295–6300, 2008.
"""
function fdichkspec(sysf::FDIModel{T}, SFDI::Union{BitMatrix,Matrix{Bool}}; sdeg::Union{Real,Missing} = missing, minimal::Bool = false,
                    FDtol::Real = 0.0001, FDGainTol::Real = 0.01, FDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing,
                    atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol, atol3::Real = atol,
                    rtol::Real = ((size(sysf.sys.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
                    fast::Bool = true) where T

   mu = length(sysf.controls)
   md = length(sysf.disturbances)
   mud = mu+md
   mf = length(sysf.faults)
   m = mu+md+mf
   strongFD = !ismissing(FDfreq)
   strongFD && !isa(FDfreq,Vector) && (FDfreq = [FDfreq])
   strongFD && ismissing(sdeg) && (sdeg = sysf.sys.Ts != 0 ? 0.95 : -0.05)

   size(SFDI,2) == mf || error("Structure matrix incompatible with the number of faults")

   N = size(SFDI,1)
   p = size(sysf.sys,1)
   strongFD && (lfreq = length(FDfreq))     # number of frequency values
   rdims = fill(0,N)
   orders = fill(-1,N)
   leastorders = fill(-1,N)

   N == 0 && (return rdims, orders, leastorders)

   for i = 1:N
       # set input groups for first design
       inddf = Vector(1:mf)[SFDI[i,:] .== false]
       indff = Vector(1:mf)[SFDI[i,:] .== true]
       isempty(indff) && (rdims[i] = 1; orders[i] = 0; leastorders[i] = minimal ? 0 : -1; continue)
       syse = [sysf.sys[:,[1:mud; mud .+ inddf; mud .+ indff]]; eye(mu,m)]
       m2f = length(indff)
       #
       # compute a left nullspace basis Q of G = [Gu Gd Gdf; I 0 0] = 0
       # obtain QR = [ Q Rf ], where Rf = Q*[Gf;0]
       QR, info = glnull(syse, m2f; sdeg, atol1, atol2, rtol, fast)
       nvec = size(QR,1);          # number of basis vectors
       # check solvability conditions
       nvec == 0 && continue

       nq = order(QR)               # order of minimal basis
       degs = reverse(info.degs)    # degrees of a minimal polynomial basis
       m2f == 0 && (rdims[i] = nvec; orders[i] = nq; leastorders[i] = minimum(degs); continue)

       indf = (p+mu) .+ Vector(1:m2f)
       feasible = true;
       if strongFD
          t = gir(QR[:,indf]; atol = atol1)
          S = fdisspec_(t, FDfreq; stabilize = true, block = !minimal, FDGainTol, atol1, atol2, rtol = 0, fast)[1]
          # check strong detectability conditions
          for ii = 1:lfreq
             all(maximum(S[:,:,ii],dims=1)) || (feasible = false; break)
          end
       else
          # check weak detectability conditions
          S = fditspec_(QR[:,indf]; block = !minimal, atol1, atol2, rtol, FDtol)
          all(maximum(S,dims=1)) || (feasible = false)
       end
       if feasible
          rdims[i] = nvec
          orders[i] = nq
          if minimal
             if nvec > 1
                 finish = false    # set termination flag
                 nout = 1          # initialize number of selected basis vectors
                 #QR = xperm(QR,nq:-1:1);  # permute states to speedup glmcover1
                 while !finish
                    # choose nout basis vectors, which potentially lead to a least order
                    # filter with rdim outputs:
                    # basesel[i,:] contains the indices of candidate basis vectors;
                    # ordsel[i]    contains the presumably achievable least orders
                    basesel, ordsel = efdbasesel(S, degs, 1, nout, false)
                    #
                    # update the synthesis using the selections of candidate vector(s),
                    # starting with the least (potentially) achievable order
                    for ibas = 1:size(basesel,1)
                        baseind = basesel[ibas] # indices of current basis selection
                        ip = [baseind; setdiff(1:nvec,baseind)][:]
                        if nout == 1
                           QRfwtest = glmcover1(QR[ip,:], 1; atol1, atol2, rtol)[1]
                           if !isempty(ordsel) && (order(QRfwtest) != ordsel[ibas])
                              @warn "fdichkspec: expected reduced order not achieved"
                           end
                        else
                           # build a linear combination of the first nout vectors
                           nrest = nvec - nout
                           h = [ rand(1,nout) zeros(T, 1, nrest) ]
                           QRfwtest = glmcover1([h; zeros(T,nrest,nout) eye(nrest)]*QR[ip,:], 1; atol1, atol2, rtol)[1]
                        end
                        # check complete fault detectability of the current design;
                        # dismiss design if check fails
                        t = gminreal(QRfwtest[:,indf]; atol1, atol2, rtol, fast)
                        if strongFD
                           Stest = fdisspec_(t, FDfreq; stabilize = true, block = true,
                                             FDGainTol, atol1, atol2, atol3, rtol = 0, fast)[1]
                        else
                           Stest = fditspec_(t; atol1, atol2, rtol, FDtol, block = true)
                        end
                        if all(Stest)
                           finish = true
                           leastorders[i] = order(QRfwtest)
                           break
                        end
                     end
                     nout += 1
                     nout > nvec && (finish = true; leastorders[i] = nq)  # least order is
                  end
             else
                  leastorders[i] = nq  #least order for nvec = 1
             end
          end
       end
   end
   return rdims, orders, leastorders
   # end FDICHKSPEC
end
"""
    S = fdigenspec(sysf::FDIModel; sdeg, FDtol, FDGainTol, FDfreq, atol, atol1, atol2, atol3, rtol, fast = true)

Generate all achievable specifications `S` for a given synthesis model `sysf` with additive faults.
Each row of the resulting binary matrix `S` contains a nonzero specification (or fault signature) which can
be achieved using a linear fault detection filter (e.g., as obtainable with the help of function [`efdisyn`](@ref)).

`FDFreq = freq` specifies a vector of real frequency values or a scalar real frquency value
for strong detectability checks (default: `FDFreq = missing`).

`FDtol = tol1` specifies the threshold `tol1` for assessing weak specifications
                      (see also function [`fditspec`](@ref)) (default: `tol1 = 0.0001`).

`FDGainTol = tol2` specifies the threshold `tol2` for assessing strong specifications,
i.e., the threshold for nonzero frequency responce gains for all frequency values specified in
`freq` (see also function [`fdisspec`](@ref)) (default: `tol2 = 0.01`).

The keyword argument `sdeg = β` specifies a prescribed stability degree `β` for the poles of the internally
generated candidate filters, such that the real parts of
filters poles must be less than or equal to `β`, in the continuous-time case, and
the magnitudes of filter poles must be less than or
equal to `β`, in the discrete-time case. If `sdeg = missing` then no then no stabilization is performed if and `FDFreq = missing`.
If `sdeg = missing` and `FDFreq = freq`, then the fllowing default values are employed : `β = -0.05`, in continuous-time case, and  `β = 0.95`,
in discrete-time case.

The rank determinations in the performed reductions
are based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C`, `D`,
the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `B`, `C`, `D` and `E`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the order of the system `sysf.sys`.
The keyword argument `atol3` is an absolute tolerance for observability tests
(default: internally determined value).
The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`.

_Method:_ The Procedure GENSPEC from [1] is implemented.
The nullspace method [2] is recursively employed to generate candidate fault detection and isolation filters,
whose internal forms provide the structure matrices corresponding to the achievable weak specifications, if `freq = missing`,
or strong specifications for the frequencies conatined in `freq`. The generation method is also described in [3].

_References:_

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques.
     Springer Verlag, 2017; sec. 5.4.

[2] A. Varga, On computing nullspace bases – a fault detection perspective.
      Proc. IFAC 2008 World Congress, Seoul, Korea, pages 6295–6300, 2008.

[3] A. Varga, On computing achievable fault signatures. Proc. SAFEPROCESS'2009, Barcelona, Spain.
"""
function fdigenspec(sysf::FDIModel{T}; sdeg::Union{Real,Missing} = missing,
                    FDtol::Real = 0.0001, FDGainTol::Real = 0.01, FDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing,
                    atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol, atol3::Real = atol,
                    rtol::Real = ((size(sysf.sys.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
                    fast::Bool = true) where T

   mu = length(sysf.controls)
   md = length(sysf.disturbances)
   mf = length(sysf.faults)
   m = mu+md+mf
   strongFD = !ismissing(FDfreq)
   strongFD && !isa(FDfreq,Vector) && (FDfreq = [FDfreq])
   strongFD && ismissing(sdeg) && (sdeg = sysf.sys.Ts != 0 ? 0.95 : -0.05)

   return  genspec_in([sysf.sys[:,1:m]; eye(T,mu,m)]; m1 = mu+md, sdeg, FDtol, FDGainTol, FDfreq,
                       atol1, atol2, atol3, rtol, fast)

   # end FDIGENSPEC
end
function genspec_in(sys::DescriptorStateSpace{T}; inp1::Vector{Int} = Int[], inp2::Vector{Int} = Int[], m1::Int = 0,
    sdeg::Union{Real,Missing} = missing, FDtol::Real = 0.0001, FDGainTol::Real = 0.01, FDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing,
    atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol, atol3::Real = atol,
    rtol::Real = ((size(sysf.sys.A,1)+1)*eps(real(float(one(T)))))*iszero(max(atol1,atol2)),
    fast::Bool = true) where T
#  Recursive function used to generate the achievable specifications.
#
#     S = genspec_in(sys::DescriptorStateSpace; inp1, inp2, md, FDfreq, FDtol, FDGainTol, atol1, atol2, atol3, rtol, fast)
#
#  determines, in the rows of the binary matrix S, all achievable fault detection specifications for the
#  partitioned state-space system sys = [sys[:,inp1] sys[:,inp2]], with the
#  corresponding input-output representation of the form,
#
#                y = Gd*d + Gf*f ,
#
#  where y is the output, d and f are the disturbance and fault
#  inputs, respectively, and Gd and Gf are the transfer-function
#  matrices of sys[:,inp1] and sys[:,inp2], respectively.
#  If inp1 and inp2 are both empty (e.g., on the first call), then
#  the value of m1 specifies the number of disturbance inputs and
#  inp1 and inp2 are set to inp1 = 1:m1 and inp2 = m1+1:m,
#  where m is the number of inputs of sys.
#  Each row of the resulting S contains the an achievable specification,
#  obtainable by using a certain scalar output fault detection and isolation filter
#  with the input-output form
#
#                r = Q*y ,
#
#  which ensures that all disturbance inputs d are decoupled
#  from the residual r (i.e., Q*Gd = 0) and the residual r is
#  sensitive to a subset of fault inputs (i.e., Q*Gf != 0).
#  If FDfreq = missing, then the i-th row S[i,:] is the weak structure matrix of Q*Gf,
#  whose element S[i,j] = true if Q*Gf[:,j] is nonzero and S[i,j] = false if Q*Gf[:,j] = 0.
#  The check for nonzero Q*Gf[:,j] is performed by using the function
#  fditspec to evaluate the corresponding weak structure matrix.
#  If FDfreq = freq specifies the frequency values in freq, then
#  S[i,j] = true if |Q*Gf[:,j]| is larger than or equal to the threshold FDGainTol
#  for all specified frequency values and S[i,j] = false otherwise.
#  For checking purposes, the function fdisspec is used to
#  evaluate the corresponding strong structure matrix.
#  For the description of keyword parameters atol1, atol2, atol3, rtol and fast
#  see the documentation of the function fdispecgen.

   # set persistent variables used for recursive calls: level, ksave, nf
   global level
   global ksave
   global nf

   # determine input information to define SYS1 and SYS2 (relevant for recursive usage)
   md = isempty(inp1) ? 0 : length(inp1)
   mf = isempty(inp2) ? 0 : length(inp2)

   # initialize recursion level
   md == 0 && mf == 0 && (level = 0)   # iteration level is set to zero on first call

   if level == 0
      # at first call the number of disturbance inputs are specified via m1
      # and the rest are fault inputs
      md = m1
      mf = size(sys,2) - md  # the number of fault inputs

      # finish if there are no faults
      mf == 0 && (return falses(0,0))

      # set nf and ksave
      nf = ismissing(FDfreq) ? 0 : length(FDfreq);
      ksave = 1;  # parameter used to reduce the total number of iterations
   end

   # compute SYSN*SYS2, where SYSN is a left nullspace of SYS1
   md > 0 && (sys = glnull(sys, mf; sdeg, atol1, atol2, rtol, fast)[1][:,end-mf+1:end])
   # compute weak specification for the current level
   S = fditspec_(sys; atol1, atol2, rtol, FDtol)

   if nf > 0
      # determine strong structure matrix
      # enforce that strong specifications are a subset of weak specifications
      #S2 = fdisspec_(sys, FDfreq; stabilize = true, atol1, atol2, rtol, FDGainTol)[1];
      t = gminreal(sys; atol1, atol2, rtol)
      Sc = maximum(S,dims=1) # include the cumulated specification
      S2 = fdisspec_(t, FDfreq; stabilize = true, atol1, atol2, rtol, FDGainTol)[1];
      ns = size(S,1)
      tv = trues(ns)
      tvc = true
      for i = 1:nf
          for j = 1:ns
              tv[j] && (tv[j] = isequal(view(S,j,:),view(S2,j,:,i)))
          end
          tvc && (tvc = isequal(Sc,maximum(view(S2,:,:,i),dims=1)))
      end
      S = S[tv,:]
      tvc && (S = [S; Sc] )
   end
   # add cummulated specification
   size(S,1) > 1 && (S = [S; maximum(S,dims=1)])

   # exit level if SYS has only one output, or if S = [] or S = 0
   if size(sys,1) == 1 || isempty(S) || all(.!S)
      level > 0 && (level -= 1)
      # eliminate duplicate or zero specifications
      size(S,1) > 1 && (S = unique(S,dims=1); S = S[any(S,dims=2)[:],:])
      return S
   end
   #   loop over all fault inputs
   for k1 = 1:mf
      k1 < ksave && continue
      level += 1
      # redefine disturbance and fault inputs
      # inp1 = 1;
      # sys.InputGroup.inp2 = 1:mf-1;
      inp1 = collect(1:1)
      inp2 = collect(1:mf-1)
      ksave = k1
      #s1 = genspec_in(sys(:,[k1, 1:k1-1,k1+1:mf]));
      s1 = genspec_in(sys[:,[[k1]; collect(1:k1-1); collect(k1+1:mf)]];
                      inp1, inp2, sdeg, FDtol, FDGainTol, FDfreq, atol1, atol2, atol3, rtol, fast)
      ksave = k1
      #  add specifications from the current level
      S = [S; s1[:,1:k1-1] falses(size(s1,1),1) s1[:,k1:mf-1]]
   end

   # eliminate duplicate or zero specifications
   size(S,1) > 1 && (S = [S; maximum(S,dims=1)])
   S = unique(S,dims=1)
   S = S[any(S,dims=2)[:],:]
   level -= 1
   return S

   # end GENSPEC_IN
end
fdichkspec(sysf::DescriptorStateSpace{T},SFDI::Union{BitMatrix,Matrix{Bool}}; kwargs...) where T = fdichkspec(fdimodset(sysf,faults = 1:size(sysf,2)),SFDI; kwargs...)
fdigenspec(sysf::DescriptorStateSpace{T}; kwargs...) where T = fdigenspec(fdimodset(sysf,faults = 1:size(sysf,2)); kwargs...)

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
"""
	GlobalNetwork <: AbstractNetwork

A global (complete) network.

# Fields
- N              :: Integer -> Number of nodes
"""
struct GlobalNetwork <: AbstractNetwork
	N :: Integer
	_props :: NetworkProperties
end

"""
	GlobalNetwork(N::Integer; directed::Bool=false)

Create a global network with `N` nodes
"""
function GlobalNetwork(N::Integer; directed::Bool=false)
	(N <= 0) && throw(ArgumentError("Number of nodes must be a positive integer!"))
	GlobalNetwork(N, NetworkProperties(N, directed))
end

function show(network::GlobalNetwork)
	print(isdirected(network) ? "Directed " : "Undirected ")
	println("Global Network")
	println("- N = $(network.N)")
end


text:
# Low-level C API for libCEED

module C

using CEnum, Libdl

# Get the path to the libCEED dynamic library, configured during the build step
# of the LibCEED.jl package.
const depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsfile)
    error("LibCEED.jl not properly installed. Please run Pkg.build(\"LibCEED\")")
end
include(depsfile)

include(joinpath(@__DIR__, "generated", "libceed_common.jl"))
include(joinpath(@__DIR__, "generated", "libceed_api.jl"))

const CEED_STRIDES_BACKEND = Ref{Ptr{CeedInt}}()
const CEED_BASIS_COLLOCATED = Ref{CeedBasis}()
const CEED_VECTOR_ACTIVE = Ref{CeedVector}()
const CEED_VECTOR_NONE = Ref{CeedVector}()
const CEED_ELEMRESTRICTION_NONE = Ref{CeedElemRestriction}()
const CEED_QFUNCTION_NONE = Ref{CeedQFunction}()
const CEED_REQUEST_IMMEDIATE = Ref{CeedRequest}()
const CEED_REQUEST_ORDERED = Ref{CeedRequest}()

function __init__()
    global libceed_handle = dlopen(libceed)
    # some global variables
    CEED_STRIDES_BACKEND[] = cglobal((:CEED_STRIDES_BACKEND, libceed))
    CEED_BASIS_COLLOCATED[] =
        unsafe_load(cglobal((:CEED_BASIS_COLLOCATED, libceed), CeedBasis))
    CEED_VECTOR_ACTIVE[] = unsafe_load(cglobal((:CEED_VECTOR_ACTIVE, libceed), CeedVector))
    CEED_VECTOR_NONE[] = unsafe_load(cglobal((:CEED_VECTOR_NONE, libceed), CeedVector))
    CEED_ELEMRESTRICTION_NONE[] =
        unsafe_load(cglobal((:CEED_ELEMRESTRICTION_NONE, libceed), CeedElemRestriction))
    CEED_QFUNCTION_NONE[] =
        unsafe_load(cglobal((:CEED_QFUNCTION_NONE, libceed), CeedQFunction))
    CEED_REQUEST_IMMEDIATE[] =
        unsafe_load(cglobal((:CEED_REQUEST_IMMEDIATE, libceed), Ptr{CeedRequest}))
    CEED_REQUEST_ORDERED[] =
        unsafe_load(cglobal((:CEED_REQUEST_ORDERED, libceed), Ptr{CeedRequest}))
end

end # module

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


text:
abstract Bootstrap{I <: Input} <: OnlineStat{I}
nobs(b::Bootstrap) = b.n
value(b::Bootstrap) = b.replicates
updatecounter!(b::Bootstrap, n2::Int) = (b.n += n2)
weight(b::Bootstrap, n2::Int) = 0.0

#-----------------------------------------------------------# BernoulliBootstrap
"""
`BernoulliBootstrap(o::OnlineStat, f::Function, r::Int = 1000)`

Create a double-or-nothing bootstrap using `r` replicates of `o` for estimate `f(o)`

Example:
```julia
BernoulliBootstrap(Mean(), mean, 1000)
```
"""
type BernoulliBootstrap{S <: OnlineStat{ScalarInput}} <: Bootstrap{ScalarInput}
    replicates::Vector{S}            # replicates of base stat
    cached_state::Vector{Float64}    # cache of replicate states
    f::Function                      # function to generate state. Ex: mean, var, std
    n::Int                           # number of observations
    cache_is_dirty::Bool
end

function BernoulliBootstrap{T <: ScalarInput}(o::OnlineStat{T}, f::Function, r::Int = 1_000)
    replicates = OnlineStat{T}[copy(o) for i in 1:r]
    cached_state = Array(Float64, r)
    return BernoulliBootstrap(replicates, cached_state, f, 0, true)
end

function _fit!(b::BernoulliBootstrap, x::Real, γ::Float64)
    for replicate in b.replicates
        if rand() > 0.5
            fit!(replicate, x)
            fit!(replicate, x)
        end
    end
    b.cache_is_dirty = true
    b
end

#-------------------------------------------------------------# PoissonBootstrap
"""
`PoissonBootstrap(o::OnlineStat, f::Function, r::Int = 1000)`

Create a poisson bootstrap using `r` replicates of `o` for estimate `f(o)`

Example:
```julia
PoissonBootstrap(Mean(), mean, 1000)
```
"""
type ScalarPoissonBootstrap{S <: OnlineStat{ScalarInput}} <: Bootstrap{ScalarInput}
    replicates::Vector{S}           # replicates of base stat
    cached_state::Vector{Float64}  # cache of replicate states
    f::Function
    n::Int                          # number of observations
    cache_is_dirty::Bool
end
function PoissonBootstrap{T <: ScalarInput}(o::OnlineStat{T}, f::Function, r::Int = 1_000)
    replicates = OnlineStat{T}[copy(o) for i in 1:r]
    cached_state = Array(Float64, r)
    ScalarPoissonBootstrap(replicates, cached_state, f, 0, true)
end

type VectorPoissonBootstrap{S <: OnlineStat{VectorInput}} <: Bootstrap{VectorInput}
    replicates::Vector{S}           # replicates of base stat
    cached_state::Matrix{Float64}  # cache of replicate states
    f::Function
    n::Int                          # number of observations
    cache_is_dirty::Bool
end
function PoissonBootstrap{T <: VectorInput}(o::OnlineStat{T}, f::Function, r::Int = 1_000)
    replicates = OnlineStat{T}[copy(o) for i in 1:r]
    cached_state = Array(Float64, length(value(o)), r)
    VectorPoissonBootstrap(replicates, cached_state, f, 0, true)
end

const unitPoissonDist = Ds.Poisson(1)
function _fit!(b::Union{ScalarPoissonBootstrap, VectorPoissonBootstrap},
        x::Union{Real,Vector}, γ::Float64)
    for replicate in b.replicates
        for repetition in 1:rand(unitPoissonDist)
            fit!(replicate, x)
        end
    end
    b.cache_is_dirty = true
    b
end

#--------------------------------------------------------------# FrozenBootstrap
# "Frozen bootstraps object are generated when two bootstrap distributions are combined
#  e.g., if they are differenced."
immutable FrozenBootstrap <: Bootstrap{ScalarInput}
    cached_state::Vector{Float64}  # cache of replicate states
    n::Int                          # number of observations
end

# "Return the value of interest for each of the `OnlineStat` replicates"
cached_state(b::FrozenBootstrap) = copy(b.cached_state)

#-----------------------------------------------------------------------# Common
function Base.show(io::IO, b::Bootstrap)
    printheader(io, string(typeof(b)))
    print_item(io, "Boostrap of", typeof(b.replicates[1]))
    print_item(io, "function", b.f)
    print_item(io, "nreplicates", length(b.replicates))
    print_item(io, "nobs", nobs(b))
end

# update cached_state' states if necessary and return their values
function cached_state(b::Bootstrap{ScalarInput})
    if b.cache_is_dirty
        for i in 1:length(b.replicates)
            b.cached_state[i] = b.f(b.replicates[i])
        end
        b.cache_is_dirty = false
    end
    return b.cached_state
end
function cached_state(b::Bootstrap{VectorInput})
    if b.cache_is_dirty
        for i in 1:length(b.replicates)
            b.cached_state[:,i] = b.f(b.replicates[i])
        end
        b.cache_is_dirty = false
    end
    return b.cached_state
end

Base.mean(b::Bootstrap{ScalarInput}) = mean(cached_state(b))
Base.std(b::Bootstrap{ScalarInput}) = std(cached_state(b))
Base.var(b::Bootstrap{ScalarInput}) = var(cached_state(b))

Base.mean(b::Bootstrap{VectorInput}) = vec(mean(cached_state(b),2))
Base.std(b::Bootstrap{VectorInput}) = vec(std(cached_state(b),2))
Base.var(b::Bootstrap{VectorInput}) = vec(var(cached_state(b),2))

replicates(b::Bootstrap) = copy(b.replicates)

# Assumes a and b are independent.
function Base.:-(a::Bootstrap{ScalarInput}, b::Bootstrap{ScalarInput})
    return FrozenBootstrap(cached_state(a) - cached_state(b), nobs(a) + nobs(b))
end

function StatsBase.confint(b::Bootstrap{ScalarInput}, coverageprob = 0.95, method=:quantile)
    states = cached_state(b)
    # If any NaN, return NaN, NaN
    if any(isnan, states)
        return (NaN, NaN)
    else
        α = 1 - coverageprob
        if method == :quantile
            return (quantile(states, α / 2), quantile(states, 1 - α / 2))
        elseif method == :normal
            norm_approx = Ds.Normal(mean(states), std(states))
            return (quantile(norm_approx, α / 2), quantile(norm_approx, 1 - α / 2))
        else
            error("method $method not recognized.  use :quantile or :normal")
        end
    end
end


text:
@testset "transformedkernel" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    s = rand(rng)
    v = rand(rng, 3)
    P = rand(rng, 3, 2)
    k = SqExponentialKernel()
    @test k ∘ IdentityTransform() === k

    kt = TransformedKernel(k, ScaleTransform(s))
    ktard = TransformedKernel(k, ARDTransform(v))
    @test kt ∘ IdentityTransform() === kt
    @test ktard ∘ IdentityTransform() === ktard
    @test kt(v1, v2) == (k ∘ ScaleTransform(s))(v1, v2)
    @test kt(v1, v2) ≈ k(s * v1, s * v2) atol = 1e-5
    @test ktard(v1, v2) == (k ∘ ARDTransform(v))(v1, v2)
    @test ktard(v1, v2) == k(v .* v1, v .* v2)
    @test (k ∘ LinearTransform(P') ∘ ScaleTransform(s))(v1, v2) ==
          ((k ∘ LinearTransform(P')) ∘ ScaleTransform(s))(v1, v2) ==
          (k ∘ (LinearTransform(P') ∘ ScaleTransform(s)))(v1, v2)

    @test repr(kt) == repr(k) * "\n\t- " * repr(ScaleTransform(s))

    TestUtils.test_interface(k, Float64)
    test_ADs(x -> SqExponentialKernel() ∘ ScaleTransform(x[1]), rand(1))

    # Test implicit gradients
    @testset "Implicit gradients" begin
        k = SqExponentialKernel() ∘ ScaleTransform(2.0)
        ps = Flux.params(k)
        X = rand(10, 1)
        x = vec(X)
        A = rand(10, 10)
        # Implicit
        g1 = Flux.gradient(ps) do
            tr(kernelmatrix(k, X; obsdim=1) * A)
        end
        # Explicit
        g2 = Flux.gradient(k) do k
            tr(kernelmatrix(k, X; obsdim=1) * A)
        end

        # Implicit for a vector
        g3 = Flux.gradient(ps) do
            tr(kernelmatrix(k, x) * A)
        end
        @test g1[first(ps)] ≈ first(g2).transform.s
        @test g1[first(ps)] ≈ g3[first(ps)]
    end

    @testset "Parameters" begin
        k = ConstantKernel(; c=rand(rng))
        c = Chain(Dense(3, 2))

        test_params(k ∘ ScaleTransform(s), (k, [s]))
        test_params(k ∘ ARDTransform(v), (k, v))
        test_params(k ∘ LinearTransform(P), (k, P))
        test_params(k ∘ LinearTransform(P) ∘ ScaleTransform(s), (k, [s], P))
        test_params(k ∘ FunctionTransform(c), (k, c))
    end
end


text:
using BinaryBuilder, Pkg

name = "xfoil_light"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/byuflowlab/xfoil_light.git", "7e1ac9cea9de6941b293a66b47198e7ab4ec4e7f")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd xfoil_light/
install_license LICENSE
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libxfoil_light", :libxfoil_light),
    LibraryProduct("libxfoil_light_cs", :libxfoil_light_cs)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ------------ Generate comilation context ------------ #

mutable struct GenerateContext{J <: CompilationOptions,
                               T <: AddressMap,
                               K <: AddressMap,
                               P <: AddressMap} <: ExecutionContext{J, K, GenerateContext}
    tr::T
    target::K
    weight::Float64
    score::Float64
    visited::Visitor
    params::P
    GenerateContext{J}(tr::T, target::K, weight::Float64, score::Float64, visited::Visitor, params::P) where {J, T, K, P} = new{J, T, K, P}(tr, target, weight, score, visited, params)
end

function Generate(::J, tr::AddressMap, target::AddressMap, params::AddressMap) where J <: CompilationOptions
    GenerateContext{J}(tr,
                       target,
                       0.0,
                       0.0,
                       Visitor(),
                       params)
end

# ------------ Propose compilation context ------------ #

mutable struct ProposeContext{J <: CompilationOptions,
                              T <: AddressMap,
                              P <: AddressMap} <: ExecutionContext{J, Nothing, ProposeContext}
    map::T
    score::Float64
    visited::Visitor
    params::P
    ProposeContext{J}(tr::T, score::Float64, visited::Visitor, params::P) where {J, T, P} = new{J, T, P}(tr, score, visited, params)
end

function Propose(opt::J, tr, params) where J
    ProposeContext{J}(tr,
                      0.0,
                      Visitor(),
                      params)
end

# ------------ Regenerate compilation context ------------ #

mutable struct RegenerateContext{J <: CompilationOptions,
                                 C <: AddressMap,
                                 T <: AddressMap,
                                 K <: AddressMap,
                                 D <: AddressMap,
                                 P <: AddressMap} <: ExecutionContext{J, K, RegenerateContext}
    prev::C
    tr::T
    target::K
    weight::Float64
    score::Float64
    discard::D
    visited::Visitor
    params::P
    RegenerateContext{J}(cl::C, tr::T, target::K, weight::Float64, score::Float64, discard::D, visited::Visitor, params::P) where {J, C, T, K, D, P} = new{J, C, T, K, D, P}(cl, tr, target, weight, score, discard, visited, params)
end

function Regenerate(::J, target::K, ps, cl::C, tr, discard) where {J <: CompilationOptions, K <: AddressMap, C <: CallSite}
    RegenerateContext{J}(cl,
                         tr,
                         target,
                         0.0,
                         0.0,
                         discard,
                         Visitor(),
                         Empty())
end

# ------------ Simulate compilation context ------------ #

mutable struct SimulateContext{J <: CompilationOptions,
                               T <: AddressMap,
                               P <: AddressMap} <: ExecutionContext{J, Nothing, SimulateContext}
    tr::T
    score::Float64
    visited::Visitor
    params::P
    SimulateContext{J}(tr::T, score::Float64, visited::Visitor, params::P) where {J, T, P} = new{J, T, P}(tr, score, visited, params)
end

function Simulate(opt::J, tr::T, params) where {J, T}
    SimulateContext{J}(tr,
                       0.0,
                       Visitor(),
                       params)
end

# ------------ Update compilation context ------------ #

mutable struct UpdateContext{J <: CompilationOptions,
                             C <: CallSite,
                             T <: AddressMap,
                             K <: AddressMap,
                             D <: AddressMap,
                             P <: AddressMap} <: ExecutionContext{J, K, UpdateContext}
    prev::C
    tr::T
    target::K
    weight::Float64
    score::Float64
    discard::D
    visited::Visitor
    params::P
    UpdateContext{J}(cl::C, tr::T, target::K, weight, score, discard::D, vs::Visitor, params::P) where {J, C, T, K, D, P} = new{J, C, T, K, D, P}(cl, tr, target, weight, score, discard, vs, params)
end

function Update(::J, select::K, ps::P, cl::CL, tr, discard) where {J <: CompilationOptions, K <: AddressMap, P <: AddressMap, CL <: CallSite}
    UpdateContext{J}(cl,
                     tr,
                     select,
                     0.0,
                     0.0,
                     discard,
                     Visitor(),
                     ps)
end

# ------------ Forward mode AD compilation context ------------ #

# Support for forward mode automatic differentiation.
mutable struct ForwardModeContext{J <: CompilationOptions,
                                  T <: Tuple,
                                  C <: AddressMap,
                                  D,
                                  P <: AddressMap} <: ExecutionContext{J, Nothing, ForwardModeContext}
    target::T
    map::C
    weight::D
    visited::Visitor
    params::P
    ForwardModeContext{J}(target::T, map::C, weight::D, visited::Visitor, params::P) where {J, T, C, D, P} = new{J, T, C, D, P}(target, map, weight, visited, params)
end

function ForwardMode(addr, params, cl, weight)
    ForwardModeContext{DefaultPipeline}(addr,
                                        cl,
                                        weight,
                                        Visitor(),
                                        params)
end

# ------------ Assess compilation context ------------ #

mutable struct AssessContext{J <: CompilationOptions,
                             M <: AddressMap,
                             P <: AddressMap} <: ExecutionContext{J, Nothing, AssessContext}
    target::M
    weight::Float64
    visited::Visitor
    params::P
    AssessContext{J}(target::M, weight::Float64, visited::Visitor, params::P) where {J, M, P} = new{J, M, P}(target, weight, visited, params)
end

function Assess(opt::J, obs::AddressMap, params) where J
    AssessContext{J}(obs,
                     0.0,
                     Visitor(),
                     params)
end

# ------------ Backpropagation compilation contexts ------------ #

abstract type BackpropagationContext{J} <: ExecutionContext{J, Nothing, BackpropagationContext} end

# Learnable parameters
mutable struct ParameterBackpropagateContext{J <: CompilationOptions,
                                             T <: CallSite,
                                             S <: AddressMap,
                                             P <: AddressMap} <: BackpropagationContext{J}
    call::T
    weight::Float64
    fillables::S
    initial_params::P
    params::Store
    param_grads::Gradients
    ParameterBackpropagateContext{J}(call::T, weight::Float64, fillables::S, initial_params::P, params::Store, param_grads::Gradients) where {J, T, S, P} = new{J, T, S, P}(call, weight, fillables, initial_params, params, param_grads)
end

function ParameterBackpropagate(call::T, sel::S, init, params, param_grads::Gradients) where {T <: CallSite, S <: AddressMap, K <: Target}
    ParameterBackpropagateContext{DefaultPipeline}(call,
                                                  0.0,
                                                  sel,
                                                  init,
                                                  params,
                                                  param_grads)
end

# Choice sites
mutable struct ChoiceBackpropagateContext{J <: CompilationOptions,
                                          T <: CallSite,
                                          S <: AddressMap,
                                          P <: AddressMap,
                                          K <: Target} <: BackpropagationContext{J}
    call::T
    weight::Float64
    fillables::S
    initial_params::P
    choices::Store
    choice_grads::Gradients
    target::K
    ChoiceBackpropagateContext{J}(call::T, weight::Float64, fillables::S, initial_params::P, choices::Store, choice_grads::Gradients, target::K) where {J, T, S, P, K} = new{J, T, S, P, K}(call, weight, fillables, initial_params, choices, choice_grads, target)
end

function ChoiceBackpropagate(call::T, fillables::S, init, choice_store, choice_grads, sel::K) where {T <: CallSite, S <: AddressMap, K <: Target}
    ChoiceBackpropagateContext{DefaultPipeline}(call,
                                                0.0,
                                                fillables,
                                                init,
                                                choice_store,
                                                choice_grads,
                                                sel)
end


text:
using ADSeismic
using ADCME
using PyPlot
using MAT
using DelimitedFiles
matplotlib.use("Agg")
close("all")

################### Generate synthetic data #####################
model_name = "../../../examples/nn_fwi/models/marmousi2-model-true.mat"
# model_name = "models/BP-model-true.mat"

## load model setting
PropagatorKernel = 2
params = load_params(model_name, vp_ref=1000, PropagatorKernel=PropagatorKernel)
src = load_acoustic_source(model_name)
rcv = load_acoustic_receiver(model_name)
vp = constant(matread(model_name)["vp"])
## For debug, only run the first source
# src = [src[1]]
# rcv = [rcv[1]]

## assemble acoustic propagator model
model = x->AcousticPropagatorSolver(params, x, vp)

## simulated wavefield
[SimulatedObservation!(model(src[i]), rcv[i]) for i = 1:length(src)]
Rs_ = [rcv[i].rcvv for i = 1:length(rcv)]
sess = Session(); init(sess)

# ts = zeros(11)
# for i = 1:11
#     d = @timed Rs = run(sess, Rs_)
#     ts[i] = d[2]
#     @info "Test $i, time = $(ts[i])"
# end
# @info mean(ts[2:end])

## visualize_wavefield if needed
ii = 1 # source number
u = run(sess, model(src[ii]).u)
p = visualize_wavefield(u, params)
saveanim(p,"test$PropagatorKernel.gif")


text:
using WhereTheWaterAlsoFlows
using Test

@testset "WhereTheWaterAlsoFlows.jl" begin
    # Write your tests here.
end


text:
#
# bounds.jl --
#
# Implement operations involving simple bound constraints.
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2015-2018, Éric Thiébaut.
# <https://github.com/emmt/OptimPackNextGen.jl>.
#

module SimpleBounds

export
    fastclamp,
    fastmax,
    fastmin,
    get_free_variables!,
    get_free_variables,
    project_direction!,
    project_variables!,
    step_limits

using Base: @propagate_inbounds

"""
    fastmin(x, y)

yields the least of `x` and `y` if neither `x` nor `y` are NaNs, or `x`
otherwise.

Calling `fastmin(x,y)` is faster than `min(x,y)` but the latter propagates
NaNs (i.e., `min(x,y)` yields NaN if any of `x` or `y` is a NaN).

"""
fastmin(x::T, y::T) where {T<:Real} = (x > y ? y : x)

"""
    fastmax(x, y)

yields the greatest of `x` and `y` if neither `x` nor `y` are NaNs, or `x`
otherwise.

Calling `fastmax(x,y)` is faster than `max(x,y)` but the latter propagates
NaNs (i.e., `max(x,y)` yields NaN if any of `x` or `y` is a NaN).

"""
fastmax(x::T, y::T) where {T<:Real} = (x < y ? y : x)

"""
    fastclamp(x, lo, hi)

yields `x` subject to simple bound constraints.  That is, returns `x` if
`lo ≤ x ≤ hi`, `lo` if 'x < lo' and `hi` if `x > hi`.  This method is
similar to `clamp` except that bounds may be `nothing` to indicate that
there is no limit for the corresponding bound and that NaNs are not
treated specially.

    fastclamp(x, lo, hi, i)

yields `i`-th variable `x[i]` subject to simple bound constraints as
specified by `lo` and `hi`.

"""
fastclamp(x::T, lo::T, hi::T) where {T<:Real} = fastmax(fastmin(x, hi), lo)
fastclamp(x::T, ::Nothing, hi::T) where {T<:Real} = fastmin(x, hi)
fastclamp(x::T, lo::T, ::Nothing) where {T<:Real} = fastmax(x, lo)

@inline @propagate_inbounds fastclamp(x::AbstractArray, lo, hi, i) =
    fastclamp(x[i], boundvalue(lo, i), boundvalue(hi, i))

"""
    boundvalue(b, i)

yields the bound value at index `i` for bounds `b`.

"""
boundvalue(b::Nothing, i) = b
boundvalue(b::Real, i) = b
@inline @propagate_inbounds boundvalue(b::AbstractArray, i) = b[i]

#-------------------------------------------------------------------------------
# PROJECTING VARIABLES

const Bound{T<:Real,N} = Union{AbstractArray{T,N},T,Nothing}

"""
    project_variables!(dst, src, lo, hi) -> dst

overwrites destination `dst` the projection of the source variables `src` in the
box whose lower bound is `lo` and upper bound is `hi`.  The destination `dst`
is returned.

This is the same as `dst = clamp.(src, lo, hi)` except that the result is
preallocated and that the operation is *much* faster (by a factor of 2-3).

"""
function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::Real,
                            hi::Real) where {T<:AbstractFloat,N}
    project_variables!(dst, src, T(lo), T(hi))
end

function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::T,
                            hi::T) where {T<:AbstractFloat,N}
    @assert size(dst) == size(src)
    @assert lo ≤ hi # this also check for NaN
    bounded_below = (lo > T(-Inf))
    bounded_above = (hi < T(+Inf))
    if bounded_below && bounded_above
        @inbounds @simd for i in eachindex(dst, src)
            dst[i] = fastclamp(src[i], lo, hi)
        end
    elseif bounded_below
        @inbounds @simd for i in eachindex(dst, src)
            dst[i] = fastclamp(src[i], lo, nothing)
        end
    elseif bounded_above
        @inbounds @simd for i in eachindex(dst, src)
            dst[i] = fastclamp(src[i], nothing, hi)
        end
    elseif dst !== src
        vcopy!(dst, src)
    end
    return dst
end

function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::DenseArray{T,N},
                            hi::Real) where {T<:AbstractFloat,N}
    project_variables!(dst, src, lo, T(hi))
end

function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::DenseArray{T,N},
                            hi::T) where {T<:AbstractFloat,N}
    @assert size(dst) == size(src) == size(lo)
    if hi < T(+Inf)
        @inbounds @simd for i in eachindex(dst, src, lo)
            dst[i] = fastclamp(src[i], lo[i], hi)
        end
    else
        @inbounds @simd for i in eachindex(dst, src, lo)
            dst[i] = fastclamp(src[i], lo[i], nothing)
        end
    end
    return dst
end

function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::Real,
                            hi::DenseArray{T,N}) where {T<:AbstractFloat,N}
    project_variables!(dst, src, T(lo), hi)
end

function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::T,
                            hi::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(dst) == size(src) == size(hi)
    if lo > T(-Inf)
        @inbounds @simd for i in eachindex(dst, src, hi)
            dst[i] = fastclamp(src[i], lo, hi[i])
        end
    else
        @inbounds @simd for i in eachindex(dst, src, hi)
            dst[i] = fastclamp(src[i], nothing, hi[i])
        end
    end
    return dst
end

function project_variables!(dst::DenseArray{T,N},
                            src::DenseArray{T,N},
                            lo::DenseArray{T,N},
                            hi::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(dst) == size(src) == size(lo) == size(hi)
    @inbounds  @simd for i in eachindex(dst, src, lo, hi)
        dst[i] = fastclamp(src[i], lo[i], hi[i])
    end
    return dst
end

#------------------------------------------------------------------------------
# PROJECTING DIRECTION

# Orientation is indicated by a singleton.
abstract type Orientation end
struct Forward  <: Orientation; end
struct Backward <: Orientation; end

const FORWARD = Forward()
const BACKWARD = Backward()

Base.convert(::Type{Orientation}, ::Union{Forward,Type{Forward}}) = FORWARD
Base.convert(::Type{Orientation}, ::Union{Backward,Type{Backward}}) = BACKWARD
Base.convert(::Type{Orientation}, s::Real) = (s > 0 ? FORWARD :
                                              s < 0 ? BACKWARD :
                                              error("invalid orientation"))
Base.sign(::Union{Forward,Type{Forward}}) = +1
Base.sign(::Union{Backward,Type{Backward}}) = -1

Orientation(x) = convert(Orientation, x)
orientation(T::DataType, x) = Orientation(x) === FORWARD ? +one(T) : -one(T)

@inline projdir(x::T, lo::T, hi::T, ::Forward, d::T) where {T<:AbstractFloat} =
    (d > zero(T) ? x < hi : x > lo) ? d : zero(T)

@inline projdir(x::T, lo::T, hi::T, ::Backward, d::T) where {T<:AbstractFloat} =
    (d < zero(T) ? x < hi : x > lo) ? d : zero(T)

@inline projdir(x::T, ::Nothing, hi::T, ::Forward, d::T) where {T<:AbstractFloat} =
    (d < zero(T) || x < hi) ? d : zero(T)

@inline projdir(x::T, ::Nothing, hi::T, ::Backward, d::T) where {T<:AbstractFloat} =
    (d > zero(T) || x < hi) ? d : zero(T)

@inline projdir(x::T, lo::T, ::Nothing, ::Forward, d::T) where {T<:AbstractFloat} =
    (d > zero(T) || x > lo) ? d : zero(T)

@inline projdir(x::T, lo::T, ::Nothing, ::Backward, d::T) where {T<:AbstractFloat} =
    (d < zero(T) || x > lo) ? d : zero(T)

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::T,
                            hi::T,
                            o::Orientation,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(dst) == size(x) == size(d)
    @assert lo ≤ hi # this also check for NaN
    bounded_above = (hi < T(+Inf))
    bounded_below = (lo > T(-Inf))
    if bounded_below && bounded_above
        @inbounds @simd for i in eachindex(dst, x, d)
            dst[i] = projdir(x[i], lo, hi, o, d[i])
        end
    elseif bounded_below
        @inbounds @simd for i in eachindex(dst, x, d)
            dst[i] = projdir(x[i], lo, nothing, o, d[i])
        end
    elseif bounded_above
        @inbounds @simd for i in eachindex(dst, x, d)
            dst[i] = projdir(x[i], nothing, hi, o, d[i])
        end
    elseif !is(dst, d)
        vcopy!(dst, d)
    end
    return dst
end

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::DenseArray{T,N},
                            hi::T,
                            o::Orientation,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(dst) == size(x) == size(d) == size(lo)
    if hi < T(+Inf)
        @inbounds @simd for i in eachindex(dst, x, d, lo)
            dst[i] = projdir(x[i], lo[i], hi, o, d[i])
        end
    else
        @inbounds @simd for i in eachindex(dst, x, d, lo)
            dst[i] = projdir(x[i], lo[i], nothing, o, d[i])
        end
    end
    return dst
end

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::T,
                            hi::DenseArray{T,N},
                            o::Orientation,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(dst) == size(x) == size(d) == size(hi)
    if lo > T(-Inf)
        @inbounds @simd for i in eachindex(dst, x, d, hi)
            dst[i] = projdir(x[i], lo, hi[i], o, d[i])
        end
    else
        @inbounds @simd for i in eachindex(dst, x, d, hi)
            dst[i] = projdir(x[i], nothing, hi[i], o, d[i])
        end
    end
    return dst
end

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::DenseArray{T,N},
                            hi::DenseArray{T,N},
                            o::Orientation,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(dst) == size(x) == size(d) == size(lo) == size(hi)
    @inbounds @simd for i in eachindex(dst, x, d, lo, hi)
        dst[i] = projdir(x[i], lo[i], hi[i], o, d[i])
    end
    return dst
end

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::Real,
                            hi::Real,
                            orient,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    project_direction!(dst, x, T(lo), T(hi), Orientation(orient), d)
end

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::DenseArray{T,N},
                            hi::Real,
                            orient,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    project_direction!(dst, x, lo, T(hi), Orientation(orient), d)
end

function project_direction!(dst::DenseArray{T,N},
                            x::DenseArray{T,N},
                            lo::Real,
                            hi::DenseArray{T,N},
                            orient,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    project_direction!(dst, x, T(lo), hi, Orientation(orient), d)
end

function project_gradient!(dst::DenseArray{T,N},
                           x::DenseArray{T,N},
                           lo::Union{Real,DenseArray{T,N}},
                           hi::Union{Real,DenseArray{T,N}},
                           d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    project_direction!(dst, x, lo, hi, BACKWARD, d)
end

#------------------------------------------------------------------------------
# COMPUTING STEP LIMITS

"""
### Compute step limits for line search

When there are separable bound constraints on the variables, the step `smin` to
the closest not yet reached bound and the step `smax` to the farthest bound are
computed by the call:
```
    (smin, smax) = step_limits(x, lo, hi, s, d)
```
where `lo` is the lower bound, `hi` is the upper bound, `x` are the current
variables and `sign(s)*d` is the search direction.

In orther words, `smin` is the smallest step which will bring at least one more
variable "out of bounds" and `smax` is the smallest step which will bring all
variables "out of bounds".  As a consequence, `0 < smin` and `0 ≤ smax`.
"""
function step_limits(x::AbstractArray{T,N},
                     lo::T,
                     hi::T,
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(d)
    @assert lo ≤ hi # this also check for NaN
    ZERO = zero(T)
    INFINITY = T(Inf)
    s = orientation(T, orient)
    bounded_below = lo > -INFINITY
    bounded_above = hi < +INFINITY
    smin = INFINITY
    @inbounds begin
        if bounded_below && bounded_above
            smax = ZERO
            @simd for i in eachindex(x, d)
                p = s*d[i]
                if p != ZERO
                    a = (p > ZERO ? hi - x[i] : lo - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                end
            end
        elseif bounded_below
            smax = ZERO
            @simd for i in eachindex(x, d)
                p = s*d[i]
                if p > ZERO
                    smax = INFINITY
                elseif p < ZERO
                    a = (lo - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                end
            end
        elseif bounded_above
            smax = ZERO
            @simd for i in eachindex(x, d)
                p = s*d[i]
                if p > ZERO
                    a = (hi - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                elseif p < ZERO
                    smax = INFINITY
                end
            end
        else
            smax = INFINITY
        end
    end
    return (smin, smax)
end

function step_limits(x::AbstractArray{T,N},
                     lo::AbstractArray{T,N},
                     hi::T,
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(lo) == size(d)
    ZERO = zero(T)
    INFINITY = T(Inf)
    s = orientation(T, orient)
    bounded_above = hi < +INFINITY
    smin = INFINITY
    smax = ZERO
    @inbounds begin
        if bounded_above
            @simd for i in eachindex(x, lo, d)
                p = s*d[i]
                if p != ZERO
                    a = (p > ZERO ? hi - x[i] : lo[i] - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                end
            end
        else
            @simd for i in eachindex(x, lo, d)
                p = s*d[i]
                if p > ZERO
                    smax = INFINITY
                elseif p < ZERO
                    a = (lo[i] - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                end
            end
        end
    end
    return (smin, smax)
end

function step_limits(x::AbstractArray{T,N},
                     lo::T,
                     hi::AbstractArray{T,N},
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(hi) == size(d)
    ZERO = zero(T)
    INFINITY = T(Inf)
    s = orientation(T, orient)
    bounded_below = lo > -INFINITY
    smin = INFINITY
    smax = ZERO
    @inbounds begin
        if bounded_below
            @simd for i in eachindex(x, hi, d)
                p = s*d[i]
                if p != ZERO
                    a = (p > ZERO ? hi[i] - x[i] : lo - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                end
            end
        else
            @simd for i in eachindex(x, hi, d)
                p = s*d[i]
                if p > ZERO
                    a = (hi[i] - x[i])/p
                    if ZERO < a < smin
                        smin = a
                    end
                    if a > smax
                        smax = a
                    end
                elseif p < ZERO
                    smax = INFINITY
                end
            end
        end
    end
    return (smin, smax)
end

function step_limits(x::AbstractArray{T,N},
                     lo::AbstractArray{T,N},
                     hi::AbstractArray{T,N},
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(lo) == size(hi) == size(d)
    ZERO = zero(T)
    s = orientation(T, orient)
    smin = T(Inf)
    smax = ZERO
    @inbounds begin
        @simd for i in eachindex(x, lo, hi, d)
            p = s*d[i]
            if p != ZERO
                # Step length to reach the upper/lower bound:
                a = (p > ZERO ? hi[i] - x[i] : lo[i] - x[i])/p
                if ZERO < a < smin
                    smin = a
                end
                if a > smax
                    smax = a
                end
            end
        end
    end
    return (smin, smax)
end

function step_limits(x::AbstractArray{T,N},
                     lo::Real,
                     hi::Real,
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    step_limits(x, T(lo), T(hi), orient, d)
end

function step_limits(x::AbstractArray{T,N},
                     lo::AbstractArray{T,N},
                     hi::Real,
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    step_limits(x, lo, T(hi), orient, d)
end

function step_limits(x::AbstractArray{T,N},
                     lo::Real,
                     hi::AbstractArray{T,N},
                     orient,
                     d::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    step_limits(x, T(lo), hi, orient, d)
end

#------------------------------------------------------------------------------
# GETTING FREE VARIABLES

@inline may_move(x::T, lo::T, ::Nothing, ::Forward, d::T) where {T<:AbstractFloat} =
    d > zero(T) || (d != zero(T) && x > lo)

@inline may_move(x::T, lo::T, ::Nothing, ::Backward, d::T) where {T<:AbstractFloat} =
    d < zero(T) || (d != zero(T) && x > lo)

@inline may_move(x::T, ::Nothing, hi::T, ::Forward, d::T) where {T<:AbstractFloat} =
    d < zero(T) || (d != zero(T) && x < hi)

@inline may_move(x::T, ::Nothing, hi::T, ::Backward, d::T) where {T<:AbstractFloat} =
    d > zero(T) || (d != zero(T) && x < hi)

@inline may_move(x::T, lo::T, hi::T, ::Forward, d::T) where {T<:AbstractFloat} =
    d != zero(T) && (d < zero(T) ? x > lo : x < hi)

@inline may_move(x::T, lo::T, hi::T, ::Backward, d::T) where {T<:AbstractFloat} =
    d != zero(T) && (d > zero(T) ? x > lo : x < hi)

"""
## Get free variables when following a direction

    sel = get_free_variables(x, lo, hi, orient, d)

yields the list of components of the variables `x` which are allowed to vary
along the search direction `sign(orient)*d` under box constraints with `lo` and
`hi` the lower and upper bounds.

If the projected gradient `gp` of the objective function is available, the free
variables can be obtained by:

    sel = get_free_variables(gp)

where the projected gradient `gp` has been computed as:

    project_direction!(gp, x, lo, hi, -1, g)

with `g` the gradient of the objective function at `x`.

There are in-place versions:

    get_free_variables!(sel, x, lo, hi, orient, d) -> sel
    get_free_variables!(sel, gp)

"""
function get_free_variables!(sel::Vector{Int},
                             gp::DenseArray{T,N}) where {T<:AbstractFloat,N}
    ZERO = zero(T)
    n = length(gp)
    resize!(sel, n)
    j = 0
    @inbounds @simd for i in 1:n
        if gp[i] != ZERO
            j += 1
            sel[j] = i
        end
    end
    return (j == n ? sel : resize!(sel, j))
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::T,
                             hi::T,
                             o::Orientation,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(d)
    @assert lo ≤ hi # this also check for NaN
    bounded_below = (lo > T(-Inf))
    bounded_above = (hi < T(+Inf))
    n = length(x)
    resize!(sel, n)
    j = 0
    if bounded_below && bounded_above
        @inbounds @simd for i in 1:n
            if may_move(x[i], lo, hi, o, d[i])
                j += 1
                sel[j] = i
            end
        end
    elseif bounded_below
        @inbounds @simd for i in 1:n
            if may_move(x[i], lo, nothing, o, d[i])
                j += 1
                sel[j] = i
            end
        end
    elseif bounded_above
        @inbounds @simd for i in 1:n
            if may_move(x[i], nothing, hi, o, d[i])
                j += 1
                sel[j] = i
            end
        end
    else
        @inbounds @simd for i in 1:n
            sel[i] = i
        end
        j = n
    end
    return (j == n ? sel : resize!(sel, j))
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::DenseArray{T,N},
                             hi::T,
                             o::Orientation,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(d) == size(lo)
    n = length(x)
    resize!(sel, n)
    j = 0
    if hi < T(+Inf)
        @inbounds @simd for i in 1:n
            if may_move(x[i], lo[i], hi, o, d[i])
                j += 1
                sel[j] = i
            end
        end
    else
        @inbounds @simd for i in 1:n
            if may_move(x[i], lo[i], nothing, o, d[i])
                j += 1
                sel[j] = i
            end
        end
    end
    return (j == n ? sel : resize!(sel, j))
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::T,
                             hi::DenseArray{T,N},
                             o::Orientation,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(d) == size(hi)
    n = length(x)
    resize!(sel, n)
    j = 0
    if lo > T(-Inf)
        @inbounds @simd for i in 1:n
            if may_move(x[i], lo, hi[i], o, d[i])
                j += 1
                    sel[j] = i
            end
        end
    else
        @inbounds @simd for i in 1:n
            if may_move(x[i], nothing, hi[i], o, d[i])
                j += 1
                sel[j] = i
            end
        end
    end
    return (j == n ? sel : resize!(sel, j))
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::DenseArray{T,N},
                             hi::DenseArray{T,N},
                             o::Orientation,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    @assert size(x) == size(d) == size(lo) == size(hi)
    n = length(x)
    resize!(sel, n)
    j = 0
    @inbounds @simd for i in 1:n
        if may_move(x[i], lo[i], hi[i], o, d[i])
            j += 1
            sel[j] = i
        end
    end
    return (j == n ? sel : resize!(sel, j))
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::Real,
                             hi::Real,
                             orient,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    get_free_variables!(sel, x, T(lo), T(hi), Orientation(orient), d)
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::DenseArray{T,N},
                             hi::Real,
                             orient,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    get_free_variables!(sel, x, lo, T(hi), Orientation(orient), d)
end

function get_free_variables!(sel::Vector{Int},
                             x::DenseArray{T,N},
                             lo::Real,
                             hi::DenseArray{T,N},
                             orient,
                             d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    get_free_variables!(sel, x, T(lo), hi, Orientation(orient), d)
end

function get_free_variables(gp::DenseArray{T,N}) where {T<:AbstractFloat,N}
    get_free_variables!(newvariablelengthvector(Int, length(gp)), gp)
end

function get_free_variables(x::DenseArray{T,N},
                            lo::Union{Real,DenseArray{T,N}},
                            hi::Union{Real,DenseArray{T,N}},
                            orient,
                            d::DenseArray{T,N}) where {T<:AbstractFloat,N}
    get_free_variables!(newvariablelengthvector(Int, length(x)),
                        x, lo, hi, orient, d)
end

@doc @doc(get_free_variables!) get_free_variables

function newvariablelengthvector(::Type{T}, n::Integer) where {T}
    vec = Vector{T}(undef, n)
    sizehint!(vec, n)
    return vec
end

#------------------------------------------------------------------------------

end # module


text:
include("module.jl")

let rulesets = @.("data/markovrules0" * string(1:5) * ".txt"),
    ruletest = @.("data/markovtest0" * string(1:5) * ".txt")
    for i in eachindex(rulesets, ruletest)
        rules = MarkovAlgos.ruleset(rulesets[i])
        println("# Example n.$i")
        println("Original:\n", readstring(ruletest[i]))
        println("Transformed:\n", MarkovAlgos.apply(ruletest[i], rules))
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export RandomPatch

"""
    RandomPatch()

A random patch. The vector has norm 1.
"""
struct RandomPatch <: AbstractAffinePatch end

struct RandomPatchState{T, N} <: AbstractAffinePatchState{N}
    v̄::PVector{T, N}
end

function state(::RandomPatch, x::PVector)
    v = copy(x)
    Random.randn!(v)
    LinearAlgebra.normalize!(v)
    conj!(v)
    RandomPatchState(v)
end

onpatch!(x::AbstractVector, state::RandomPatchState) = onpatch!(x, state.v̄)
evaluate!(u, state::RandomPatchState, x::PVector) = evaluate_patch!(u, state.v̄, x)
jacobian!(U, state::RandomPatchState, x::PVector) = jacobian_patch!(U, state.v̄, x)


text:
module ISODATA

using StatsBase
using Distances

# source files
include("base.jl")
include("centers.jl")
include("clusters.jl")
include("evaluation.jl")
include("kmeans.jl")
include("params.jl")

export base, base_conv
export randomcenters, kmppcenters, isodata_s
export clusterfit
export kmeans
export summarize, iso_param

end


text:
import Convex.MultiplyAtom, Convex.Vexity

export multivexity

abstract type MultiAffineVexity <: Convex.Vexity             end
abstract type MultiConvexVexity <: Convex.Vexity             end
abstract type MultiConcaveVexity <: Convex.Vexity            end

multi(v::Convex.Vexity) = MultiConvexVexity()

# notes:
# * may add conflicts twice since (x,y) is different (as a tuple) from (y,x)

# multicurvature computes the curvature of the function,
# supposing ...
# records id to variable mapping and
# adds conflicts to `conflicts`
function multicurvature(x::Convex.AbstractExpr,
                        conflicts = Set{@compat Tuple{UInt64, UInt64}}(),
                        id2var = Dict{UInt64, Variable}())
  return curvature(x)
end

function multicurvature(x::MultiplyAtom,
                        conflicts = Set{@compat Tuple{UInt64, UInt64}}(),
                        id2var = Dict{UInt64, Variable}())
  c1, c2 = x.children[1], x.children[2]
  m1 = multivexity(c1, conflicts, id2var)
  m2 = multivexity(c2, conflicts, id2var)
  if m1 != ConstVexity() && m2 != ConstVexity()
    # one or the other of c1 and c2 must be constant for
    # the expression to be DCP
    # so we add the complete bipartite graph on the variables in c1 and c2 to the conflict edges E
    union!(conflicts, outerproduct(variablesin(c1, id2var), variablesin(c2, id2var)))
  end
  return m1 + m2
end

multivexity(x::Variable, conflicts, id2var) = vexity(x)
multivexity(x::Constant, conflicts, id2var) = vexity(x)

function multivexity(x::Convex.AbstractExpr,
                     conflicts = Set{@compat Tuple{UInt64, UInt64}}(),
                     id2var = Dict{UInt64, Variable}())
  monotonicities = monotonicity(x)
  vex = multicurvature(x, conflicts, id2var)
  for i = 1:length(x.children)
    vex += monotonicities[i] * multivexity(x.children[i], conflicts, id2var)
  end
  return vex
end

function multivexity(x::Convex.AbstractExpr)
  conflicts = Set{@compat Tuple{UInt64, UInt64}}()
  id2var = Dict{UInt64, Variable}()
  return multivexity(x, conflicts, id2var), conflicts, id2var
end

# constraints are multiconvex iff they are convex
multivexity(x::Constraint) = vexity(x)

# a problem is multiconvex if there exists a partition of the variables
# in the problem such that
  # 1. the objective is convex in each element of the partition holding the others constant
  # 2. all constraints are convex
  # 3. all variables participating in any constraint are in the same element of the partition
# this function checks multiconvexity of a problem
# and returns the partition that certifies multiconvexity
function multivexity(problem::Problem)
  bad_vex = [ConcaveVexity, NotDcp]

  conflicts = Set{@compat Tuple{UInt64, UInt64}}()
  id2var = Dict{UInt64, Variable}()
  obj_vex = multivexity(problem.objective, conflicts, id2var)
  if problem.head == :maximize
    obj_vex = -obj_vex
  end
  # the objective must be convex
  if typeof(obj_vex) in bad_vex
    warn("Expression not DMCP compliant: objective not multiconvex")
  end

  # the constraints must be convex and must not violate the conflicts graph
  m = length(problem.constraints)
  idpartition = Array(Set{UInt64}, m)
  activeparts = fill(true, m)
  for i=1:m
    constraint = problem.constraints[i]
    if typeof(vexity(constraint)) in bad_vex
      warn("Expression not DMCP compliant: constraint $i not multiconvex")
    end
    vars = variablesin(constraint, id2var)
    # if the variables in two constraints intersect,
    # the the union of their variables must be in the same partition
    for j=(1:i-1)[activeparts[1:i-1]] # for every subset known to be in the partition
      if length(intersect(idpartition[j], vars)) > 0
        union!(vars, idpartition[j])
        activeparts[j] = false
      end
    end
    idpartition[i] = vars
    activeparts[i] = true
  end
  idpartition = idpartition[activeparts]

  # check that the variables that must be simultaneously active don't violate the conflicts graph
  for s in idpartition
    if !isstableset(s, conflicts)
      warn("Expression not DMCP compliant:
        variables appearing together in a constraint
        must be in different subsets of the partition
        to ensure the objective is convex.")
    end
  end

  # ok, now partition the remaining variables
  # and add them to the appropriate subsets in idpartition
  stablepartition!(conflicts, idpartition)

  # make list of vars and stablesets by Variable rather than by id
  vars = Variable[id2var[id] for id in keys(id2var)]
  stablesetvars = Array{Variable,1}[Variable[id2var[id] for id in s] for s in idpartition]

  return obj_vex, stablesetvars, vars
end