text:
macro test_throws_unwrap(terr, args...)
    _test_throws_unwrap(terr, args...)
end
function _test_throws_unwrap(terr, ex; to_match=[])
    @gensym rerr
    match_expr = Expr(:block)
    for m in to_match
        @assert m.head == :(=)
        push!(match_expr.args, :(@test $rerr.$(m.args[1]) == $(m.args[2])))
    end
    quote
        $rerr = try
            $(esc(ex))
        catch err
            Dagger.Sch.unwrap_nested_exception(err)
        end
        @test $rerr isa $terr
        $match_expr
    end
end
function _test_throws_unwrap(terr, args...)
    ex = last(args)
    to_match = args[1:end-1]
    _test_throws_unwrap(terr, ex; to_match=to_match)
end


text:
using Klimakoffer
using Printf

NT = 48 # Number of time-steps per year

year_start = 1950
year_end = 1951
year_delta = year_end - year_start + 1

co2_concentration_start = 315.0 # in [ppm] for year 1950
co2_concentration_yearly_delta = 1.39

co2_concentration_yearly = zeros(Float64,year_delta)
mean_temperature_yearly = zeros(Float64, year_delta)
# co2 scenario: simple linear growth per year
for years in year_start:year_end
  iter = years+1-year_start
  co2_concentration_yearly[iter] = co2_concentration_start + co2_concentration_yearly_delta * (iter - 1)
end

mesh = Mesh()
model = Model(mesh, NT; co2_concentration = co2_concentration_start)
discretization = Discretization(mesh, model, NT)

for years in year_start:year_end
  iter = years+1-year_start
  set_co2_concentration!(model, co2_concentration_yearly[iter])
  mean_temperature_yearly[iter] = compute_equilibrium!(discretization; verbose = false)
  @printf "Global Temerature in year %4i with CO2 concentration %.2f [ppm] is T=%.3f °C" years co2_concentration_yearly[iter] mean_temperature_yearly[iter]
  println("")
end

ext: jl


text:
#using .StandardRules
using Test, Printf, InteractiveFiction
using TestSetExtensions

function ex2()
    #toggle_dbg!()
    function run_checks(w, rb)
        for thing in w.objects_by_type[:thing]
            if description(thing) == ""
                sayln!(w, @sprintf "%s has no description." name(thing))
            end
        end
    end
    w = make_blank_world("Bic")
    when_play_begins!(w, :run_through_property_checks_at_the_start_of_play, run_checks)
    staffroom = add_room!(w, "Staff Break Room", :dark)
    orange = add_thing!(w, "orange"; description="It's a small hard pinch-skinned thing from the lunch room,
            probably with lots of pips and no juice.")
    pen = add_thing!(w, "Bic pen")
    napkin = add_thing!(w, "napkin"; description="Slightly crumpled")
    #world.now_player_carries(orange, pen, napkin)
    out = run_world!(w)
    out
end

function ex4()
    w = make_blank_world("Slightly Wrong")
    awning = add_room!(w, "Awning"; description="A tan awning is stretched on tent poles over the dig-site, providing a little shade to " *
    "the workers here; you are at the bottom of a square twenty feet on a side, marked out with pegs and lines of string. " *
    "Uncovered in the south face of this square is an awkward opening into the earth.")
    slightlywrong = add_room!(w, "Slightly Wrong Chamber"; description=DynamicString(w,
    "{main_description}A mural on the far wall depicts a woman with a staff,
    tipped with a pine-cone. She appears to be watching you.",
    Dict("main_description" => _ ->
    slightlywrong |> unvisited ? "When you first step into the room, you are bothered by the sense that
    something is not quite right: perhaps the lighting, perhaps the angle of the
    walls." : "")), map_connections=Dict(:south_of => awning))
    run_world!(w)
    #test_with_actions(world, ["looking", ["going", "south"], "looking"])
    #finish_up_world()
    #world
end

@testset ExtendedTestSet "Chapter 3: Things" begin
    @test ex2() ==
    "-----------------
-------Bic-------
-----------------

Staff Break Room

Bic pen has no description."
    #@test ex3()
end


text:
module DANCE

using Base: @kwdef
using Julog

export Domain, Problem, Action, Event, State
export Term, Compound, Var, Const
export parse_domain, parse_problem, parse_pddl, @pddl, @pddl_str
export write_domain, write_problem, write_pddl
export load_domain, load_problem, preprocess
export save_domain, save_problem
export get_static_predicates, get_static_functions
export satisfy, evaluate, find_matches
export init_state, goal_state, initialize, transition, simulate
export get_preconditions, get_effect
export effect_diff, precond_diff, effect_dist, update!, update
export available, relevant, execute, execpar, execseq, regress, trigger
export use_available_action_cache!, use_relevant_action_cache!
export clear_available_action_cache!, clear_relevant_action_cache!

max_stars_repo_name: timo-cmd/DANCE

max_issues_repo_name: timo-cmd/Twister

max_forks_repo_name: timo-cmd/DANCE


text:
#=

amss.jl

@Author: Shunsuke Hori

=#

include("lucas_stokey.jl")
using LS

using QuantEcon
using NLopt
using NLsolve
using Dierckx

Para=LS.Para
Planners_Allocation_Sequential = LS.Planners_Allocation_Sequential
"""
Bellman equation for the continuation of the Lucas-Stokey Problem
"""
type BellmanEquation{TI<:Integer,TF<:AbstractFloat}
    para::Para
    S::TI
    xbar::Array{TF}
    time_0::Bool
    z0::Array{Array}
    cFB::Vector{TF}
    nFB::Vector{TF}
    xFB::Vector{TF}
    zFB::Vector{Array}
end

"""
Compute the planner's allocation by solving Bellman
equation.
"""
type Planners_Allocation_Bellman
    para::Para
    mc::MarkovChain
    S::Integer
    T::BellmanEquation
    mugrid::Vector
    xgrid::Vector
    Vf::Array
    policies
end

"""
Initializes the type from the calibration Para
"""
function Planners_Allocation_Bellman(para::Para,mugrid::AbstractArray)
    G = para.G
    S = size(para.Pi,1) # number of states
    mc=MarkovChain(para.Pi)
    mugrid = collect(mugrid)
    #now find the first best allocation
    Vf, policies, T, xgrid = solve_time1_bellman(para, mugrid)
    T.time_0 = true #Bellman equation now solves time 0 problem
    return Planners_Allocation_Bellman(para,mc,S,T,mugrid,xgrid, Vf, policies)
end

"""
Solve the time  1 Bellman equation for calibration Para and initial grid mugrid
"""
function solve_time1_bellman{TR<:Real}(para::Para,mugrid::AbstractArray{TR})
    Pi = para.Pi
    S = size(para.Pi,1)

    #First get initial fit from lucas stockey solution.
    #Need to change things to be ex_ante
    PP = LS.Planners_Allocation_Sequential(para)

    function incomplete_allocation(PP::LS.Planners_Allocation_Sequential,
                                    mu_::AbstractFloat, s_::Integer)
        c,n,x,V = LS.time1_value(PP,mu_)
        return c,n,dot(Pi[s_,:], x), dot(Pi[s_,:], V)
    end
    cf = Array{Function}(S,S)
    nf = Array{Function}(S,S)
    xprimef = Array{Function}(S,S)
    Vf = Array{Function}(S)
    # TR=Float64
    xgrid = Array{TR}(S,length(mugrid))
    for s_ in 1:S
        c = Array{TR}(length(mugrid),S)
        n = Array{TR}(length(mugrid),S)
        x = Array{TR}(length(mugrid))
        V = Array{TR}(length(mugrid))
        for (i_mu, mu) in enumerate(mugrid)
            c[i_mu,:], n[i_mu,:], x[i_mu], V[i_mu] =
                incomplete_allocation(PP, mu, s_)
        end
        xprimes = repmat(x,1,S)
        xgrid[s_,:] = x
        for sprime = 1:S
            splc = Spline1D(x[end:-1:1],c[:,sprime][end:-1:1],k=3)
            spln = Spline1D(x[end:-1:1],n[:,sprime][end:-1:1],k=3)
            splx = Spline1D(x[end:-1:1],xprimes[:,sprime][end:-1:1],k=3)
            cf[s_,sprime] = y -> splc(y)
            nf[s_,sprime] = y -> spln(y)
            xprimef[s_,sprime] = y -> splx(y)
        end
        splV = Spline1D(x[end:-1:1],V[end:-1:1],k=3)
        Vf[s_] = y -> splV(y)
    end

    policies = [cf,nf,xprimef]

    #create xgrid
    xbar = [maximum(minimum(xgrid)), minimum(maximum(xgrid))]
    xgrid = collect(linspace(xbar[1],xbar[2],length(mugrid)))

    #Now iterate on Bellman equation
    T = BellmanEquation(para,xgrid,policies)
    diff = 1.0
    while diff > 1e-6
        PF = (i_x, x,s) -> get_policies_time1(T, i_x, x, s, Vf, xbar)
        Vfnew, policies = fit_policy_function(T, PF, xgrid)

        diff = 0.0
        for s=1:S
            diff = max(diff,maxabs((Vf[s].(xgrid)-Vfnew[s].(xgrid))./Vf[s].(xgrid)))
        end

        println("diff = $diff")
        Vf = copy(Vfnew)
    end

    return Vf, policies, T, xgrid
end

"""
Fits the policy functions
"""
function fit_policy_function{TF<:AbstractFloat}(
                    T::BellmanEquation,
                    PF::Function, xgrid::AbstractArray{TF})
    S = T.S
    # preallocation
    # TF = Float64
    PFvec = Array{TF}(4S+1,length(xgrid))
    cf = Array{Function}(S,S)
    nf = Array{Function}(S,S)
    xprimef = Array{Function}(S,S)
    TTf = Array{Function}(S,S)
    Vf = Array{Function}(S)
    # fit policy fuctions
    for s_ in 1:S
        for (i_x, x) in enumerate(xgrid)
            PFvec[:,i_x] = PF(i_x,x,s_)
        end
        splV = Spline1D(xgrid, PFvec[1,:], k=3)
        Vf[s_] = y -> splV(y)
        for sprime=1:S
            splc = Spline1D(xgrid, PFvec[1+sprime,:], k=3)
            spln = Spline1D(xgrid, PFvec[1+S+sprime,:], k=3)
            splxprime = Spline1D(xgrid, PFvec[1+2S+sprime,:], k=3)
            splTT = Spline1D(xgrid, PFvec[1+3S+sprime,:], k=3)
            cf[s_,sprime] = y -> splc(y)
            nf[s_,sprime] = y -> spln(y)
            xprimef[s_,sprime] = y -> splxprime(y)
            TTf[s_,sprime] = y -> splTT(y)
        end
    end
    policies = (cf, nf, xprimef, TTf)
    return Vf, policies
end

"""
Computes Tau given c,n
"""
function Tau(pab::Planners_Allocation_Bellman,
            c::Array,n::Array)
    para = pab.para
    Uc, Un = para.Uc(c,n), para.Un(c,n)
    return 1+Un./(para.Theta .* Uc)
end

Tau(pab::Planners_Allocation_Bellman, c::AbstractFloat,n::AbstractFloat) =
    Tau(pab, [c], [n])

"""
Finds the optimal allocation given initial government debt B_ and state s_0
"""
function time0_allocation(pab::Planners_Allocation_Bellman,
                           B_::Real,s0::Integer)
    T, Vf = pab.T, pab.Vf
    xbar = T.xbar
    z0 = get_policies_time0(T, B_, s0, Vf, xbar)

    c0, n0, xprime0, T0 = z0[2], z0[3], z0[4], z0[5]
    return c0, n0, xprime0, T0
end

"""
Simulates planners policies for T periods
"""
function simulate{TI<:Integer,TF<:AbstractFloat}(pab::Planners_Allocation_Bellman,
                   B_::TF,s_0::TI,T::TI,sHist=nothing)
    para, mc, Vf, S = pab.para, pab.mc, pab.Vf, pab.S
    Pi, Uc = para.Pi, para.Uc
    cf,nf,xprimef,TTf = pab.policies

    if sHist == nothing
        sHist = QuantEcon.simulate(mc,T,init=s_0)
    end

    cHist=Array{TF}(T)
    nHist=Array{TF}(T)
    Bhist=Array{TF}(T)
    xHist=Array{TF}(T)
    TauHist=Array{TF}(T)
    THist=Array{TF}(T)
    muHist=Array{TF}(T)

    #time0
    cHist[1],nHist[1],xHist[1],THist[1]  =
        time0_allocation(pab, B_,s_0)
    TauHist[1] = Tau(pab,cHist[1],nHist[1])[s_0]
    Bhist[1] = B_
    muHist[1] = Vf[s_0](xHist[1])

    #time 1 onward
    for t in 2:T
        s_, x, s = sHist[t-1], xHist[t-1], sHist[t]
        c = Array{TF}(S)
        n = Array{TF}(S)
        xprime = Array{TF}(S)
        TT = Array{TF}(S)
        for sprime=1:S
            c[sprime], n[sprime], xprime[sprime], TT[sprime] =
                cf[s_,sprime](x), nf[s_,sprime](x),
                xprimef[s_,sprime](x), TTf[s_,sprime](x)
        end

        Tau_val = Tau(pab, c, n)[s]
        u_c = Uc(c,n)
        Eu_c = dot(Pi[s_,:], u_c)

        muHist[t] = Vf[s](xprime[s])

        cHist[t], nHist[t], Bhist[t], TauHist[t] = c[s], n[s], x/Eu_c, Tau_val
        xHist[t],THist[t] = xprime[s],TT[s]
    end
    return cHist,nHist,Bhist,xHist,TauHist,THist,muHist,sHist
end

"""
Initializes the class from the calibration Para
"""
function BellmanEquation{TF<:AbstractFloat}(para::Para,
        xgrid::AbstractVector{TF}, policies0::Array)
    S = size(para.Pi,1) # number of states
    xbar = [minimum(xgrid),maximum(xgrid)]
    time_0 = false
    z0 = Array{Array}(length(xgrid),S)
    cf, nf, xprimef =
        policies0[1], policies0[2], policies0[3]
    for s in 1:S
        for (i_x, x) in enumerate(xgrid)
            cs=Array{TF}(S)
            ns=Array{TF}(S)
            xprimes=Array{TF}(S)
            for j=1:S
                cs[j], ns[j], xprimes[j] =
                    cf[s,j](x), nf[s,j](x), xprimef[s,j](x)
            end
            z0[i_x,s] =
                vcat(cs,ns,xprimes,zeros(S))
        end
    end
    cFB, nFB, IFB, xFB, zFB = LS.find_first_best(para, S, 2)
    return BellmanEquation(para,S,xbar,time_0,z0,cFB,nFB,xFB,zFB)
end

"""
Finds the optimal policies
"""
function get_policies_time1{TF<:AbstractFloat}(T::BellmanEquation,
        i_x::Integer, x::TF, s_::Integer, Vf::Array{Function}, xbar::Vector)
    para, S = T.para, T.S
    beta, Theta, G, Pi = para.beta, para.Theta, para.G, para.Pi
    U,Uc,Un = para.U, para.Uc, para.Un

    S_possible = sum(Pi[s_,:].>0)
    sprimei_possible = find(Pi[s_,:].>0)

    function objf(z,grad)
        c, xprime = z[1:S_possible], z[S_possible+1:2S_possible]
        n=(c+G[sprimei_possible])./Theta[sprimei_possible]
        Vprime = Vector{TF}(S_possible)
        for (si, s) in enumerate(sprimei_possible)
            Vprime[si] = Vf[s](xprime[si])
        end
        return -dot(Pi[s_,sprimei_possible], U(c,n)+beta*Vprime)
    end

    function cons(out,z,grad)
        c, xprime, TT =
            z[1:S_possible], z[S_possible+1:2S_possible], z[2S_possible+1:3S_possible]
        n=(c+G[sprimei_possible])./Theta[sprimei_possible]
        u_c = Uc(c,n)
        Eu_c = dot(Pi[s_,sprimei_possible], u_c)
        out .= x*u_c/Eu_c - u_c.*(c-TT)-Un(c,n).*n - beta*xprime
    end
    function cons_no_trans(out,z,grad)
        c, xprime =
            z[1:S_possible], z[S_possible+1:2S_possible]
        n=(c+G[sprimei_possible])./Theta[sprimei_possible]
        u_c = Uc(c,n)
        Eu_c = dot(Pi[s_,sprimei_possible], u_c)
        out .= x*u_c/Eu_c - u_c.*c-Un(c,n).*n - beta*xprime
    end

    if para.transfers == true
        lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1], zeros(S_possible))
        if para.n_less_than_one == true
            ub = vcat(100*ones(S_possible)-G[sprimei_possible], ones(S_possible)*xbar[2], 100*ones(S_possible))
        else
            ub = vcat(100*ones(S_possible), ones(S_possible)*xbar[2], 100*ones(S_possible))
        end
            init = vcat(T.z0[i_x,s_][sprimei_possible],
                T.z0[i_x,s_][2S+sprimei_possible], T.z0[i_x,s_][3S+sprimei_possible])
        opt = Opt(:LN_COBYLA, 3S_possible)
        equality_constraint!(opt, cons, zeros(S_possible))
    else
        lb = vcat(zeros(S_possible), ones(S_possible)*xbar[1])
        if para.n_less_than_one == true
            ub = vcat(ones(S_possible)-G[sprimei_possible], ones(S_possible)*xbar[2])
        else
            ub = vcat(ones(S_possible), ones(S_possible)*xbar[2])
        end
        init = vcat(T.z0[i_x,s_][sprimei_possible],
                    T.z0[i_x,s_][2S+sprimei_possible])
        opt = Opt(:LN_COBYLA, 2S_possible)
        equality_constraint!(opt, cons_no_trans, zeros(S_possible))
    end
    init[init.> ub] = ub[init.> ub]
    init[init.< lb] = lb[init.< lb]

    min_objective!(opt, objf)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 10000000)
    maxtime!(opt, 10)
    ftol_rel!(opt, 1e-8)
    ftol_abs!(opt, 1e-8)

    (minf,minx,ret) = optimize(opt, init)

    if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED &&
         ret != :FTOL_REACHED && ret != :MAXTIME_REACHED
        error("optimization failed: ret = $ret")
    end

    T.z0[i_x,s_][sprimei_possible] = minx[1:S_possible]
    T.z0[i_x,s_][S+sprimei_possible] = minx[1:S_possible]+G[sprimei_possible]
    T.z0[i_x,s_][2S+sprimei_possible] = minx[S_possible+1:2S_possible]
    if para.transfers == true
        T.z0[i_x,s_][3S+sprimei_possible] = minx[2S_possible+1:3S_possible]
    else
        T.z0[i_x,s_][3S+sprimei_possible] = zeros(S)
    end

    return vcat(-minf, T.z0[i_x,s_])
end

"""
Finds the optimal policies
"""
function get_policies_time0(T::BellmanEquation,
        B_::Real,s0::Integer,Vf::Array{Function},xbar::Vector)
    para = T.para
    beta, Theta, G = para.beta, para.Theta, para.G
    U, Uc, Un = para.U, para.Uc, para.Un

    function objf(z,grad)
        c, xprime = z[1], z[2]
        n=(c+G[s0])/Theta[s0]
        return -(U(c,n)+beta*Vf[s0](xprime))
    end

    function cons(z,grad)
        c, xprime, TT = z[1], z[2], z[3]
        n=(c+G[s0])/Theta[s0]
        return -Uc(c,n)*(c-B_-TT)-Un(c,n)*n - beta*xprime
    end
    cons_no_trans(z,grad) = cons(vcat(z,0.0), grad)

    if para.transfers == true
        lb = [0.0, xbar[1], 0.0]
        if para.n_less_than_one == true
            ub = [1-G[s0], xbar[2], 100]
        else
            ub = [100.0, xbar[2], 100.0]
        end
        init = vcat(T.zFB[s0][1], T.zFB[s0][3], T.zFB[s0][4])
        init = [0.95124922, -1.15926816,  0.0]
        opt = Opt(:LN_COBYLA, 3)
        equality_constraint!(opt, cons)
    else
        lb = [0.0, xbar[1]]
        if para.n_less_than_one == true
            ub = [1-G[s0], xbar[2]]
        else
            ub = [100, xbar[2]]
        end
        init = vcat(T.zFB[s0][1], T.zFB[s0][3])
        init = [0.95124922, -1.15926816]
        opt = Opt(:LN_COBYLA, 2)
        equality_constraint!(opt, cons_no_trans)
    end
    init[init.> ub] = ub[init.> ub]
    init[init.< lb] = lb[init.< lb]

    min_objective!(opt, objf)
    lower_bounds!(opt, lb)
    upper_bounds!(opt, ub)
    maxeval!(opt, 100000000)
    maxtime!(opt, 30)
    # ftol_rel!(opt, 1e-16)
    # ftol_abs!(opt, 1e-16)

    (minf,minx,ret) = optimize(opt, init)

    if ret != :SUCCESS && ret != :ROUNDOFF_LIMITED && ret != :MAXEVAL_REACHED && ret != :FTOL_REACHED
        error("optimization failed: ret = $ret")
    end

    if para.transfers == true
        return -minf, minx[1], minx[1]+G[s0], minx[2], minx[3]
    else
        return -minf, minx[1], minx[1]+G[s0], minx[2], 0.0
    end
end


text:
module EnKFs

# package code goes here
include("Localizations/Locals.jl")
include("EnKFs/EnKFs.jl")

end # module


text:
## Acceptance probability calculations

# Jumptypes are:
# bg, move, add,remove, split, merge

function checkbounds(rjs::RJStruct, teststate::BAMFState,idx::Int)
    outofbounds=false;
    if (teststate.x[idx] < -rjs.bndpixels)||(teststate.y[idx] < -rjs.bndpixels)||
        (teststate.x[idx] > rjs.data.sz+rjs.bndpixels)||(teststate.y[idx] > rjs.data.sz+rjs.bndpixels)
        outofbounds=true;
    end
    return outofbounds
end

function accept_move(rjs::RJStruct, currentstate::BAMFState, teststate::BAMFState,idx::Int)

    if idx<1 return 0 end
    genmodel!(currentstate, rjs, rjs.modeldata)
    genmodel!(teststate, rjs, rjs.testdata)
    if (teststate.n>0)&& (minimum(teststate.photons)<0)
        println(currentstate)
        println(teststate)
    end

    LR=likelihoodratio(rjs.modeldata, rjs.testdata, rjs.data)
    PR=pdf(rjs.prior_photons,teststate.photons[idx])/
        pdf(rjs.prior_photons,currentstate.photons[idx])

    α = PR*LR
    return α
end

function accept_bg(rjs::RJStruct, currentstate::BAMFState, teststate::BAMFState,idx::Int)
    if teststate.bg<0 return 0 end

    genmodel!(currentstate, rjs, rjs.modeldata)
    genmodel!(teststate, rjs, rjs.testdata)

    LR=likelihoodratio(rjs.modeldata, rjs.testdata, rjs.data)
    PR=pdf(rjs.prior_background,teststate.bg)/
        pdf(rjs.prior_background,currentstate.bg)

    α = PR*LR
    return α
end

function accept_add(rjs::RJStruct, currentstate::BAMFState, teststate::BAMFState,idx::Int)
    genmodel!(currentstate, rjs, rjs.modeldata)
    genmodel!(teststate, rjs, rjs.testdata)
    LLR=likelihoodratio(rjs.modeldata, rjs.testdata, rjs.data)
    #proposal probability
    residuum=calcresiduum(rjs.testdata,rjs.data)
    jj=Int(min(rjs.data.sz,max(1,round(teststate.x[idx]))))
    ii=Int(min(rjs.data.sz,max(1,round(teststate.y[idx]))))
    p=arraypdf(residuum,ii,jj)
    α = LLR*(rjs.data.sz+rjs.bndpixels)^2/p
    #println(("add: ",α,ii,jj,teststate.photons[end]))
    return α
end

function accept_remove(rjs::RJStruct, currentstate::BAMFState, teststate::BAMFState,idx::Int)
    #use the inverse function
    if idx<1 return 0 end
    α=accept_add(rjs,teststate,currentstate,idx)
    return 1/α
end

function accept_split(rjs::RJStruct, currentstate::BAMFState, teststate::BAMFState,vararg=(Int,Int,Float32,Float32,Float32))

    idx1,idx2,u1,u2,u3=vararg

    if idx1==0
        return 0
    end

    # reject if split gives value outside x,y prior
    if checkbounds(rjs,teststate,idx1)
        return 0
    end

    if checkbounds(rjs,teststate,idx2)
        return 0
    end

    genmodel!(currentstate, rjs, rjs.modeldata)
    genmodel!(teststate, rjs, rjs.testdata)
    LLR=likelihoodratio(rjs.modeldata, rjs.testdata, rjs.data)

    #proposal probability
    split_std = rjs.σ_split
    N=Normal(0,split_std)
    p=pdf(N,u2)*pdf(N,u3)

    #ratio on intensity priors
    IR=pdf(rjs.prior_photons,teststate.photons[idx1])*pdf(rjs.prior_photons,teststate.photons[idx2])/
        pdf(rjs.prior_photons,currentstate.photons[idx1])

    #ratio of XY priors
    XYPR=1/(rjs.data.sz+rjs.bndpixels)^2
    #Jacobian
    J=currentstate.photons[idx1]/(1-u1)^2
    PR=LLR*IR*XYPR/p*J
    α = PR/teststate.n
    #println(("add: ",α,ii,jj,teststate.photons[end]))
    return α

end

function accept_merge(rjs::RJStruct, currentstate::BAMFState, teststate::BAMFState,vararg=(Int,Int,Float32,Float32,Float32))
    idx1,idx2,u1,u2,u3=vararg
    if idx1<1 return 0 end
    α=accept_split(rjs,teststate,currentstate,(idx1,idx2,u1,u2,u3))
    return 1/α
end


text:
module Ceres

using CeresCApi_jll
import Libdl

include("libceres.jl")
include("loss_function.jl")

using ForwardDiff

export Problem, AddResidualBlock!, solve, CreateQuaternionParameterization,
        SetLocalParameterization!, SetParameterBlockConstant!, # CostFunction,

        AngleAxisRotatePoint, two_pose_point3d_dist_error, SE3_edge_error,
        reprojection_error,

        R2quat, quat2R, quaternion_parameterization,

        # loss
        LossFunction, EmptyLoss, HuberLoss

quaternion_parameterization = nothing

function __init__()

    h = CeresCApi_jll.libceres_c_api_handle
    ccall(Libdl.dlsym(h, :ceres_init), Cvoid, ())

    quat_param_ptr = ccall(
        Libdl.dlsym(h, :ceres_create_quaternion_parameterization),
        Ptr{Cvoid},
        ()
    )

    global quaternion_parameterization = QuaternionParameterization(quat_param_ptr)

end

```
This is for that ForwardDiff can only take vector of primitive types
```
function recover_params_shape(params_vec::Vector, params_sizes::Vector{Int64})
    parameters::Vector{SubArray} = Vector{SubArray}(undef, 0)
    sidx = 1
    for i in 1:length(params_sizes)
        push!(parameters, view(params_vec, sidx:sidx + params_sizes[i] - 1))
        sidx = sidx + params_sizes[i]
    end
    return parameters
end

function cost_function_autodiff(func,
                                user_data::Vector{Float64},
                                parameters::Vector{Vector{Float64}},
                                residuals::Vector{Float64},
                                jacobians::Vector{Vector{Float64}})::Int32

    resres = func(user_data, parameters)
    for idx = 1:length(residuals)
        residuals[idx] = resres[idx]
    end

    params_sizes = [length(v) for v in parameters]
    params_vec = reduce(vcat, parameters)
    f_params = (params) -> func(user_data, recover_params_shape(params, params_sizes))

    J = ForwardDiff.jacobian(f_params, params_vec)
    # , ForwardDiff.JacobianConfig(f_params, parameters)

    if (length(jacobians) == 0)
        return true;
    end

    sidx = 1
    for i = 1:length(jacobians)
      # we need a row major jacobian
        jvec = vec(view(J, :, sidx:sidx + params_sizes[i] - 1)')
        if !((length(jacobians[i]) == 0) | (length(jacobians[i]) == length(jvec)))
            println("jacobian vector size wrong")
        end
        for j = 1:length(jacobians[i])
            jacobians[i][j] = jvec[j]
        end
        sidx = sidx + params_sizes[i]
    end
    return true;
end

# // (double* user_data,
# //  double** parameters,
# //  double* residuals,
# //  double** jacobians,
# //  int user_data_size,
# //  int num_residuals,
# //  int num_parameter_block,
# //  int* parameter_block_sizes)
"""
typedef int (*ceres_cost_function_t)(void* user_data,
                                     double** parameters,
                                     double* residuals,
                                     double** jacobians);
"""
function costfunctionwrap(cost_function_jl,
                        user_data_ptr::Ptr{Cdouble},
                        params_pptr::Ptr{Ptr{Cdouble}},
                        residuals_ptr::Ptr{Cdouble},
                        jacobians_pptr::Ptr{Ptr{Cdouble}},
                        user_data_size::Cint,
                        num_residuals::Cint,
                        num_parameter_block::Cint,
                        parameter_block_sizes_ptr::Ptr{Cint})::Cint

    user_data = unsafe_wrap(Array, user_data_ptr, user_data_size, own=false)

    pbsizes = unsafe_wrap(Array, parameter_block_sizes_ptr, num_parameter_block, own=false)
    params_ptr = unsafe_wrap(Array, params_pptr, num_parameter_block, own=false)
    params = Vector{Vector{Float64}}(undef, 0)

    jacobians_ptr = Vector{Ptr{Cdouble}}(undef, num_parameter_block)
    fill!(jacobians_ptr, C_NULL)
    if (jacobians_pptr != C_NULL)
        jacobians_ptr = unsafe_wrap(Array, jacobians_pptr, num_parameter_block, own=false)
    end

    jacobians = Vector{Vector{Float64}}(undef, num_parameter_block)
    fill!(jacobians, Vector{Float64}(undef, 0))

    for i = 1:num_parameter_block
        push!(params, unsafe_wrap(Array, params_ptr[i], pbsizes[i], own=false))
        if jacobians_ptr[i] != C_NULL
            jacobians[i] = unsafe_wrap(Array, jacobians_ptr[i],
                                     pbsizes[i] * num_residuals, own=false)
        end
    end

    residuals = unsafe_wrap(Array, residuals_ptr, num_residuals, own=false)

    return cost_function_autodiff(cost_function_jl, user_data, params, residuals, jacobians)
end

struct Problem
    c_ptr::Ptr{Cvoid}
    # costfunctions::Vector{CostFunction}

    # prevent GC collect back params ptr vector
    params_ptr_vec::Vector{Vector{Ptr{Cdouble}}}
    cost_function_data::Vector{Vector{Float64}}
    params_sizes_int32::Vector{Vector{Cint}}
end

function Problem()
    problemc = createproblem_c()

    Problem(problemc,
      Vector{Vector{Ptr{Cdouble}}}(undef, 0),
      Vector{Vector{Float64}}(undef, 0),
      Vector{Vector{Cint}}(undef, 0))
end

function AddResidualBlock!(problem::Problem,
                           cost_function_jl,
                           cost_function_data::Vector{Float64},
                           params::Vector{Vector{Float64}},
                           residuals_num::Int64;
                           loss_function::TL=EmptyLoss()
) where TL <: LossFunction
    params_sizes_int32 = [Cint(length(p)) for p in params]
    pbsizes_ptr = Base.unsafe_convert(Ptr{Cint}, params_sizes_int32)
    # pb_sizes_ptr = Base.unsafe_convert(Ptr{Cint}, params_sizes)
    push!(problem.params_sizes_int32, params_sizes_int32)

    params_ptr = Vector{Ptr{Cdouble}}(undef, 0)
    for i = 1:length(params)
       # User must protect params from GC
        push!(params_ptr, Base.unsafe_convert(Ptr{Cdouble}, params[i]))
    end

    params_pptr = Base.unsafe_convert(Ptr{Ptr{Cdouble}}, params_ptr)
    push!(problem.params_ptr_vec, params_ptr)

    dataptr = Base.unsafe_convert(Ptr{Cdouble}, cost_function_data)
    push!(problem.cost_function_data, cost_function_data)
    # push!(problem.costfunctions, costfunction)
    ############################################################################
    # params_sizes_int32 = convert(Vector{Int32}, params_sizes)

    lambdaf = (user_data_ptr::Ptr{Cdouble},
               params_pptr::Ptr{Ptr{Cdouble}},
               residuals_ptr::Ptr{Cdouble},
               jacobians_pptr::Ptr{Ptr{Cdouble}},
               user_data_size::Cint,
               num_residuals::Cint,
               num_parameter_block::Cint,
               parameter_block_sizes_ptr::Ptr{Cint}) -> costfunctionwrap(
                              cost_function_jl,
                              user_data_ptr,
                              params_pptr,
                              residuals_ptr,
                              jacobians_pptr,
                              user_data_size,
                              num_residuals,
                              num_parameter_block,
                              parameter_block_sizes_ptr)
    costcb_c = @cfunction($lambdaf,
            Cint,
            (Ptr{Cdouble}, Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Ptr{Cdouble}}, Cint, Cint, Cint, Ptr{Cint}))

    loss_function_ptr = if TL == EmptyLoss
        C_NULL
    else
        @assert loss_function.c_ptr != C_NULL
        @show loss_function.c_ptr
        loss_function.c_ptr
    end

    val = addresidualblock_c(problem.c_ptr,
                            costcb_c.ptr, # cfunction
                            dataptr,
                            convert(Cint, length(cost_function_data)),
                            loss_function_ptr,  # loss function c ptr
                            convert(Ptr{Cdouble}, C_NULL),
                            convert(Cint, residuals_num),
                            convert(Cint, length(params)),
                            pbsizes_ptr,
                            params_pptr)
    return val

    ############################################################################

end

################# local paramterization ####################

abstract type LocalParameterization end

struct QuaternionParameterization <: LocalParameterization
    c_ptr::Ptr{Cvoid}
end

function SetLocalParameterization!(problem::Problem,
                                   params::Vector{Float64},
                                   local_parameterization::QuaternionParameterization)# where T<:

    param_ptr = Base.unsafe_convert(Ptr{Cdouble}, params)
    set_local_parameterization_c(problem.c_ptr, param_ptr, local_parameterization.c_ptr)
end

function CreateQuaternionParameterization()
    quater_c_ptr = create_quaternion_parameterization_c()
    return QuaternionParameterization(quater_c_ptr)
end

######################### solve ###############################

function solve(problem::Problem; max_iter_num::Int64=100, solver_type::Int64=1)
    solve_c(problem.c_ptr, max_iter_num, solver_type)
end

include("parameter_block.jl")

end # module


text:
using RICoverAlg

# Queueing problems
m = 2;  n = 40 ; lmd = 1.0
g = ones(n); g /=sum(g)
U = zeros(m, n)
for i in (1:n)
   U[2,i] = 28 - i
end
Z = exp.(U/lmd)

@test cancel_zero_alternatives(U,g) == (U,g,findall(a->(!isapprox(a,0; atol = 0)), g))
@test do_cover_iterations(Z,g,[.5,.5]) ≈ [0.2901315939203102, 0.7098684060796899]
@test run_cover_algorithm(Z,g,
                  initialIterations  = 10,
                  aTolSetChoiceProbToZero = 10e-6,
                  aTolerance4Stopping = 0.0,
                  maxIterations = 9999)[2] ≈ 8.880930086724845
@test cancel_zero_alternatives(Z,g) == (Z,g,collect(1:1:40))
@test cancel_zero_alternatives(Z,[0.,1.]) == (Z[[2],:],[1.],[2])

# Matejka & McKay 2015, RBT
U = [0 1 0 1;
    0 0 1 1;
    1/2 1/2 1/2 1/2]
lmd = 0.4
g = [1+1; 1+1]./4
Z,shiftterm=  exp.((U[:,[1,4]].-.5)/lmd)   ,  .5
m = size(Z)[1]
p = ones(m)./m
@test do_cover_iterations(Z,g,p) == [0.2522629174954305
 0.2522629174954305
 0.495474165009139]


text:
@testset "113.path-sum-ii.jl" begin
    root = TreeNode{Int}([5, 4, 8, 11, nothing, 13, 4, 7, 2, nothing, nothing, 5, 1])
    @test path_sum_113(root, 22) == [[5, 4, 11, 2], [5, 8, 4, 5]]

    @test path_sum_113(TreeNode{Int}([1, 2, 3]), 5) == Vector{Vector{Int}}[]

    @test path_sum_113(TreeNode{Int}([1, 2]), 0) == Vector{Vector{Int}}[]

    root = TreeNode{Int}([1, 2, nothing, 3, nothing, 4, nothing, 5])
    @test path_sum_113(root, 15) == [[1, 2, 3, 4, 5]]
end


text:
module JuliaCon

juliacon2021() = println("Welcome to JuliaCon 2021! Find more information on https://juliacon.org/2021/.")

export juliacon2021

end


text:
# acc satisfiability script
using OVERTVerify
using LazySets
using Dates
using JLD2

# ACC CONTROLLER
# 3x20 controller here: https://github.com/souradeep-111/sherlock/tree/master/systems_with_networks/ARCH_2019/ACC
controller = "acc_controller"
controller_filepath = "nnet_files/jmlr/"*controller*".nnet"
println("Controller is: ", controller)
query = OvertQuery(
    ACC,  # problem
    controller_filepath,    # network file
    Id(),      	# last layer activation layer Id()=linear, or ReLU()=relu
    "MIP",     	# query solver, "MIP" or "ReluPlex"
    55,        	# ntime
    0.1,       	# dt
    -1,        	# N_overt
    )

# x1,x2,x3 are lead vehicle variables
# x4,x5,x6 are ego vehicle variables
x_lead = [90,110]
v_lead = [32, 32.2]
gamma_lead = [0,0]
x_ego = [10,11]
v_ego = [30, 30.2]
gamma_ego = [0, 0]
var_list = [x_lead, v_lead, gamma_lead, x_ego, v_ego, gamma_ego]
input_set = Hyperrectangle(
    low=[v_range[1] for v_range in var_list],
    high=[v_range[2] for v_range in var_list]
    )
# test using Constraint on output instead of hyperrectangle!
# desired property to be proven is:
# D_rel = x_lead - x_ego >= D_safe = D_default + T_gap * v_ego
# x_lead + -1*x_ego - T_gap*v_ego >= D_default = 10
# BUT WE MUST NEGATE THIS (and so we flip >= to <=) (and then want it to always be unsat)
# THe measurement model gives x_lead + -1*x_ego - T_gap*v_ego  so we just have to pass it through with coeff=1.
output_constraint = Constraint([1.], :(<=), 10)
# NOTE: this constraint is written for use with the measurement, so we want to apply the measurement model here, and we set apply_meas to true in the call to symbolic_satisfiability below

t1 = Dates.time()
SATus, vals, stats = symbolic_satisfiability(query, input_set, output_constraint, apply_meas=true)
t2 = Dates.time()
dt = t2 - t1

JLD2.@save "src/examples/jmlr/data/acc_satisfiability_data_55.jld2" query input_set output_constraint SATus vals stats dt controller


text:
const GLOBAL_DEBUG = 1

import Debugger

# macro if_debug(expr)
#     return quote
#         if GLOBAL_DEBUG > 0 # && isinteractive()
#             $(esc(expr))
#         end
#     end
# end

macro if_debug(expr) end


text:
module IncompressibleModels

using KernelAbstractions: @index, @kernel, Event, MultiEvent
using KernelAbstractions.Extras.LoopInfo: @unroll

using Oceananigans.Utils: launch!

import Oceananigans: fields

#####
##### IncompressibleModel definition
#####

include("incompressible_model.jl")
include("non_dimensional_model.jl")
include("show_incompressible_model.jl")
include("set_incompressible_model.jl")

#####
##### Time-stepping IncompressibleModels
#####

"""
    fields(model::IncompressibleModel)

Returns a flattened `NamedTuple` of the fields in `model.velocities` and `model.tracers`.
"""
fields(model::IncompressibleModel) = merge(model.velocities, model.tracers)

include("update_hydrostatic_pressure.jl")
include("update_incompressible_model_state.jl")
include("pressure_correction.jl")
include("velocity_and_tracer_tendencies.jl")
include("calculate_tendencies.jl")
include("correct_incompressible_immersed_tendencies.jl")

end # module


text:
using TropicalNumbers
using Test

@testset "tropical" begin
    include("tropical.jl")
end

@testset "counting_tropical" begin
    include("counting_tropical.jl")
end


text:
struct Line
    password::String
    min::Int
    max::Int
    letter::Char
end

function parse_line(l)
    r = split(l, ':')
    rule = r[1]
    password = r[2]

    # Eep! Arrays are 1-indexed
    rule_parts = split(rule, " ")

    range = rule_parts[1]
    range_parts = split(range, "-")
    min = range_parts[1]
    max = range_parts[2]

    letter = rule_parts[2]

    return Line(strip(password), parse(Int64, min), parse(Int64, max), letter[1])
end

function is_valid(item)
    println(item)
    count = 0
    for c in item.password
        if c == item.letter
            count += 1
        end
    end
    return count >= item.min && count <= item.max
end

function is_valid_v2(item)
    match1 = item.password[item.min] == item.letter
    match2 = item.password[item.max] == item.letter
    return xor(match1, match2)
end

function run(fname)
    f = open(fname, "r")
    lines = readlines(f)
    valid = map((x) -> is_valid_v2(parse_line(x)), lines)
    result = count(x -> x == true, valid)
    close(f)
    println("Ran ", fname, " got ", result)
    return result
end

@assert run("2ex.txt") == 1

run("2.txt")


text:
using OrdinaryDiffEq, Test, LinearAlgebra

T = 1000.0; Ttr = 0.0;
d0=1e-9; threshold=10^4*d0;
dt = 0.1; diff_eq_kwargs = Dict()

@inbounds function eom_lorenz!(du, u, p, t)
  σ = 10.0; ρ = 28.0; β = 8/3
  du[1] = σ*(u[2]-u[1])
  du[2] = u[1]*(ρ-u[3]) - u[2]
  du[3] = u[1]*u[2] - β*u[3]
end

prob = OrdinaryDiffEq.ODEProblem(eom_lorenz!, [0., 10.0, 0], (zero(T), T))
integ1 = init(prob, Tsit5(); diff_eq_kwargs...)
@. prob.u0 = prob.u0 + d0/sqrt(3)
integ2 = init(prob, Tsit5(); diff_eq_kwargs...)
integ1.opts.advance_to_tstop=true
integ2.opts.advance_to_tstop=true

dist = d0
λ = zero(eltype(integ1.u))
i = 0
tvector = dt:dt:T

for τ in tvector
  # evolve until rescaling:
  push!(integ1.opts.tstops, τ); step!(integ1)
  push!(integ2.opts.tstops, τ); step!(integ2)
  dist = norm(integ1.u .- integ2.u)
  # Rescale:
  if dist ≥ threshold
    # Rescale and reset everything:
    integ2.u .= integ1.u #.+ (integ2.u .- integ1.u)./a
    u_modified!(integ2, true)
    set_proposed_dt!(integ2,integ1)
    break
  end
end

τ = tvector[end]
push!(integ1.opts.tstops, τ); step!(integ1)
push!(integ2.opts.tstops, τ); step!(integ2)
@test dist = norm(integ1.u .- integ2.u) == 0


text:
export zoo_cnf,
    zoo_dnf,
    zoo_sdd,
    load_logic_circuit,
    load_smooth_logic_circuit,
    load_struct_smooth_logic_circuit,
    load_cnf,
    load_dnf,
    zoo_cnf_file,
    zoo_dnf_file,
    zoo_lc_file,
    zoo_psdd_file,
    zoo_sdd_file

#####################
# loaders from model zoo
#####################

const zoo_version = "/Circuit-Model-Zoo-0.1.2"

zoo_cnf(name) =
    load_cnf(zoo_cnf_file(name))

zoo_dnf(name) =
    load_dnf(zoo_dnf_file(name))

# zoo file identifiers

zoo_cnf_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/cnfs/$name"

zoo_dnf_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/dnfs/$name"

zoo_lc_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/lcs/$name"

zoo_psdd_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/psdds/$name"

zoo_sdd_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/sdds/$name"

zoo_sdd(name) =
    load_logic_circuit(zoo_sdd_file(name))

#####################
# general parser infrastructure for circuits
#####################

# The `ParserCombinator` library works correctly but is orders of magnitude too slow.
# Instead here we hardcode some simpler parsers to speed things up

"""
    load_logic_circuit(file::String)::PlainLogicCircuit

Load a logical circuit from file.
Support file formats:
 * ".sdd" for SDD files
 * ".psdd" for PSDD files
 * ".circuit" for Logistic Circuit files
"""
function load_logic_circuit(file::String)::PlainLogicCircuit
    compile_logical(parse_circuit_file(file))
end

"""
    load_smooth_logic_circuit(file::String)::PlainLogicCircuit

Load a smooth logic circuit from file.
Support file formats:
 * ".psdd" for PSDD files
 * ".circuit" for Logistic Circuit files
"""
function load_smooth_logic_circuit(file::String)::PlainLogicCircuit
    compile_smooth_logical(parse_circuit_file(file))
end

"""
    load_struct_smooth_logic_circuit(circuit_file::String, vtree_file::String)::Tuple{StructLogicCircuit,PlainVtree}

Load a smooth structured logic circuit and its vtree from file.
Supported circuit file formats:
 * ".psdd" for PSDD files
 * ".circuit" for Logistic Circuit files
Supported vtree file formats:
 * ".vtree" for Vtree files
"""
function load_struct_smooth_logic_circuit(circuit_file::String, vtree_file::String)::Tuple{StructLogicCircuit,PlainVtree}
    circuit_lines = parse_circuit_file(circuit_file)
    vtree_lines = parse_vtree_file(vtree_file)
    compile_smooth_struct_logical(circuit_lines, vtree_lines)
end

#####################
# parse based on file extension
#####################

function parse_circuit_file(file::String)::CircuitFormatLines
    if endswith(file,".circuit")
        parse_lc_file(file)
    elseif endswith(file,".psdd")
        parse_psdd_file(file)
    elseif endswith(file,".sdd")
        parse_sdd_file(file)
    else
        throw("Cannot parse this file type as a circuit: $file")
    end
end

#####################
# parser of logistic circuit file format
#####################

const parens = r"\(([^\)]+)\)"

function parse_lc_decision_line(ln::String)::DecisionLine{LCElement}
    @assert startswith(ln, "D")
    head::SubString, tail::SubString = split(ln,'(',limit=2)
    head_tokens = split(head)
    head_ints::Vector{UInt32} = map(x->parse(UInt32,x),head_tokens[2:4])
    elems_str::String = "("*tail
    elems = Vector{LCElement}()
    for x in eachmatch(parens::Regex, elems_str)
        tokens = split(x[1], limit=3)
        weights::Vector{Float64} = map(x->parse(Float64,x), split(tokens[3]))
        elem = LCElement(parse(UInt32,tokens[1]), parse(UInt32,tokens[2]), weights)
        push!(elems,elem)
    end
    DecisionLine(head_ints[1],head_ints[2],head_ints[3],elems)
end

function parse_lc_literal_line(ln::String)::WeightedLiteralLine
    @assert startswith(ln, "T") || startswith(ln, "F")
    tokens = split(ln)
    head_ints = map(x->parse(UInt32,x),tokens[2:4])
    weights = map(x->parse(Float64,x), tokens[5:end])
    lit = var2lit(head_ints[3])
    if startswith(ln, "F")
       lit = -lit # negative literal
    end
    WeightedLiteralLine(head_ints[1],head_ints[2],lit,true,weights)
end

function parse_comment_line(ln::String)
    @assert startswith(ln, "c")
    CircuitCommentLine(lstrip(chop(ln, head = 1, tail = 0)))
end

function parse_lc_header_line(ln::String)
    @assert (ln == "Logistic Circuit") || (ln == "Logisitic Circuit")
    LcHeaderLine()
end

function parse_bias_line(ln::String)::BiasLine
    @assert startswith(ln, "B")
    tokens = split(ln)
    weights = map(x->parse(Float64,x), tokens[2:end])
    BiasLine(weights)
end

parse_lc_file(file::String)::CircuitFormatLines = open(parse_lc_file, file)

function parse_lc_file(file::Core.IO)::CircuitFormatLines
    q = Vector{CircuitFormatLine}()
    for ln in eachline(file)
        @assert !isempty(ln)
        if ln[1] == 'D'
            push!(q, parse_lc_decision_line(ln))
        elseif ln[1] == 'T' || ln[1] == 'F'
            push!(q, parse_lc_literal_line(ln))
        elseif ln[1] == 'c'
            push!(q, parse_comment_line(ln))
        elseif ln[1] == 'L'
            push!(q, parse_lc_header_line(ln))
        elseif ln[1] == 'B'
            push!(q, parse_bias_line(ln))
        else
            error("Cannot parse logistic circuit file format line '$ln'")
        end
    end
    q
end

#####################
# parser for PSDD circuit file format
#####################

function parse_psdd_decision_line(ln::String)::DecisionLine{PSDDElement}
    @assert startswith(ln, "D")
    tokens = split(ln)
    head_ints::Vector{UInt32} = map(x->parse(UInt32,x),tokens[2:4])
    elems = Vector{PSDDElement}()
    for (p,s,w) in Iterators.partition(tokens[5:end],3)
        prime = parse(UInt32,p)
        sub = parse(UInt32,s)
        weight = parse(Float64,w)
        elem = PSDDElement(prime, sub, weight)
        push!(elems,elem)
    end
    DecisionLine(head_ints[1],head_ints[2],head_ints[3],elems)
end

function parse_psdd_true_leaf_line(ln::String)::WeightedNamedConstantLine
    @assert startswith(ln, "T")
    tokens = split(ln)
    @assert length(tokens)==5
    head_ints = map(x->parse(UInt32,x),tokens[2:4])
    weight = parse(Float64,tokens[5])
    WeightedNamedConstantLine(head_ints[1],head_ints[2],head_ints[3],weight)
end

function parse_literal_line(ln::String, normalized::Bool)::UnweightedLiteralLine
    @assert startswith(ln, "L")
    tokens = split(ln)
    @assert length(tokens)==4 "line has too many tokens: $ln"
    head_ints = map(x->parse(UInt32,x),tokens[2:3])
    lit = parse(Int32,tokens[4])
    UnweightedLiteralLine(head_ints[1],head_ints[2],lit,normalized)
end

function parse_psdd_header_line(ln::String)
    @assert startswith(ln, "psdd")
    PsddHeaderLine(parse(Int,split(ln)[2]))
end

parse_psdd_file(file::String)::CircuitFormatLines = open(parse_psdd_file, file)

function parse_psdd_file(file::Core.IO)::CircuitFormatLines
    q = Vector{CircuitFormatLine}()
    for ln in eachline(file)
        @assert !isempty(ln)
        if ln[1] == 'D'
            push!(q, parse_psdd_decision_line(ln))
        elseif ln[1] == 'T'
            push!(q, parse_psdd_true_leaf_line(ln))
        elseif ln[1] == 'L'
            push!(q, parse_literal_line(ln, true))
        elseif ln[1] == 'c'
            push!(q, parse_comment_line(ln))
        elseif startswith(ln,"psdd")
            push!(q, parse_psdd_header_line(ln))
        else
            error("Cannot parse PSDD file format line '$ln'")
        end
    end
    q
end

#####################
# parser for SDD circuit file format
#####################

function parse_sdd_decision_line(ln::String)::DecisionLine{SDDElement}
    @assert startswith(ln, "D")
    tokens = split(ln)
    head_ints::Vector{UInt32} = map(x->parse(UInt32,x),tokens[2:4])
    elems = Vector{SDDElement}()
    for (p,s) in Iterators.partition(tokens[5:end],2)
        prime = parse(UInt32,p)
        sub = parse(UInt32,s)
        elem = SDDElement(prime, sub)
        push!(elems,elem)
    end
    DecisionLine(head_ints[1],head_ints[2],head_ints[3],elems)
end

function parse_sdd_constant_leaf_line(ln::String)::AnonymousConstantLine
    @assert startswith(ln, "T") || startswith(ln, "F")
    tokens = split(ln)
    @assert length(tokens)==2
    AnonymousConstantLine(parse(UInt32,tokens[2]), startswith(ln, "T"), false)
end

function parse_sdd_header_line(ln::String)
    @assert startswith(ln, "sdd")
    SddHeaderLine(parse(Int,split(ln)[2]))
end

parse_sdd_file(file::String)::CircuitFormatLines = open(parse_sdd_file, file)

function parse_sdd_file(file::Core.IO)::CircuitFormatLines
    q = Vector{CircuitFormatLine}()
    for ln in eachline(file)
        @assert !isempty(ln)
        if ln[1] == 'D'
            push!(q, parse_sdd_decision_line(ln))
        elseif ln[1] == 'T' || ln[1] == 'F'
            push!(q, parse_sdd_constant_leaf_line(ln))
        elseif ln[1] == 'L'
            push!(q, parse_literal_line(ln, false))
        elseif ln[1] == 'c'
            push!(q, parse_comment_line(ln))
        elseif startswith(ln,"sdd")
            push!(q, parse_sdd_header_line(ln))
        else
            error("Cannot parse SDD file format line '$ln'")
        end
    end
    q
end

#####################
# loader for CNF/DNF file format
#####################

"""
    load_cnf(file::String; dual=false)::PlainLogicCircuit

Load a CNF/DNF as a logic circuit from file.
Suppported file formats:
* ".cnf" for CNF files  (`dual=False`)
* ".dnf" for DNF files  (`dual=True`)
"""
function load_cnf(file::String; dual=false)::PlainLogicCircuit
    @assert !dual && endswith(file, ".cnf") || dual && endswith(file, ".dnf")

    Clause = dual ? Plain⋀Node : Plain⋁Node
    # record the current clause
    clause = Clause([])
    # linearized clauses (disjunctions)
    clauses = Vector{Clause}()

    # literal cache is responsible for making leaf literals nodes unique and adding them to `circuit`
    lit_cache = Dict{Lit,PlainLiteralNode}()

    open(file) do file

        for ln in eachline(file)
            @assert !isempty(ln)
            if ln[1] == 'c' ||
                !dual && startswith(ln, "p cnf") ||
                dual && startswith(ln, "p dnf")
                # skip comment and header lines
                continue
            elseif ln[1] == '%'
                # some CNF files appear to have trailing lines after '%' which should not be parsed
                break
            else
                tokens = split(ln)
                for token in tokens
                    occursin(r"^\s*[-]?[0-9]+\s*$", token) || error("Cannot parse CNF file format line '$ln'")
                    literal = parse(Lit, token)
                    if literal == 0
                        push!(clauses, clause)
                        clause = Clause([])
                    else
                        push!(clause.children, get!(lit_cache, literal) do
                            PlainLiteralNode(literal)
                        end)
                    end
                end
            end
        end

    end

    @assert isempty(clause.children)

    dual ? disjoin(clauses) : conjoin(clauses)
end

"""
    load_dnf(file::String)::PlainLogicCircuit

Load a DNF as a logical circuit from file.
"""
load_dnf(file::String)::PlainLogicCircuit =
    load_cnf(file; dual=true)


text:
"""
    internal_M(time, mu, timesteps)

Return `(numberOfSteps, kout)` where `numberOfSteps` is the required number of
internal simulation steps so that the probability of having a noise jump in a
timestep is `jumpProb`. `kout` tells how many internal time steps are there
between timesteps with output.

`timesteps` = number of output steps.
"""
function internal_M(time, mu, timesteps)
    temp=convert(Int64,cld(time,jumpProb/mu )+1) # minimum number of internal steps

    if time/(temp-1)>0.1
        temp=ceil(Int,time/0.1+1)
    end
    if temp<=timesteps                            # if temp<timesteps then I simply use the number of output steps
        kout=1
        return timesteps, kout
    else

    # otherwise we choose the smallest multiple of (timesteps-1) bigger than the required number of internal steps -1 :
    # numberOfSteps=kout*(timesteps-1)+1

        kout=cld((temp-1),(timesteps-1))
        # in finding kout we have to use (temp-1) and (timesteps-1) because the time step=1 is fixed for both
        # internal and output steps.

        return ((timesteps-1)*kout+1),kout

    end
end

internal_dt(time,mu,timesteps)=time/(internal_M(time,mu,timesteps)[1]-1) # unit of time given the number of internal steps

rescaledMu(time,mu,timesteps)=dt(time,mu,timesteps)*mu

"""
    psi = superposition_state(N::Integer)

Return the superposition state for the `N`-dimensional Hilbert space. `psi` is
a `Array{Complex{Float64}}`.

```jldoctest
julia> superposition_state(4)
4-element Array{Complex{Float64},1}:
 0.5 + 0.0im
 0.5 + 0.0im
 0.5 + 0.0im
 0.5 + 0.0im
```
"""
function superposition_state(N::Integer)
    return ones(ComplexF64, N) ./ sqrt(N)
end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
# test ModelSetup/Medium/mesh.jl
include("typedefine.jl")
include("mesh.jl")

#=== Case 1 layered medium ===#
vp = [5000.0,5300,5800,6000]
vs = [3500.0,3200,4000,4500]
rho = 2.8
depth = [700,200,400,600]
Horizon = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
nT = nothing

Model = model(vp,vs,rho,depth,Horizon,pkf,T,ext,iflag,dx,dz,dt,nT)

#=== Case 2 layered medium ===#
vp = 5000:100:6000
vs = 3500:100:4500
rho = 2.8
depth = 100:100:1100
Horizon = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
Tn = nothing

Model = model(vp,vs,rho,depth,Horizon,pkf,T,ext,iflag,dx,dz,dt,Tn)

#=== Case 3 homogeneous medium ===#
vp = 5000
vs = 3500
rho = 2.8
depth = 1000
Horizon = 1000
pkf = 10
T = 2
ext = 10
iflag = 2
dx = nothing
dz = nothing
dt = nothing
Tn = nothing

Model = model(vp,vs,rho,depth,Horizon,pkf,T,ext,iflag,dx,dz,dt,Tn)


text:
try
    @eval using Revise
  catch e
    @warn "Revise init" exception=(e, catch_backtrace())
  end


text:
# Tables.jl interface

Tables.isrowtable(::Type{<:TableRowIterator}) = true
Tables.columnnames(tr::TableRow) = Tuple(tr.index.column_labels)
Tables.getcolumn(tr::TableRow, nm::Symbol) = getdata(tr, nm)

function writetable(filename::AbstractString, x; kw...)

	_as_vector(y::AbstractVector) = y
	_as_vector(y) = collect(y)

    writetable(filename, Any[_as_vector(c) for c in Tables.Columns(x)], collect(Symbol, Tables.columnnames(x)); kw...)
end


text:
using JuMP, Ipopt
using DelimitedFiles
using PyPlot
import ForwardDiff

### JuMP with Ipopt

#---FUNCTIONS FOR CONSTRAINTS---

# set up boundary constraint wrapper function at specific index
# it takes in x as separate elements instead of as a vector (needed for Ipopt)
function boundary_wrapper_atindex(ind, x...)

    # combine the inputs into a single vector
    x = collect(x)

    # include globals for optimization with Ipopt
    global xlast
    global boundary_constraint_values

    # check if the constraint calculation is necessary
    if !isequal(x, xlast)

        global boundary_constraint_gradients
        global spacing_constraint_values
        global spacing_constraint_gradients

        # perform all constraint calculations for this set of design variables
        boundary_constraint_values, boundary_constraint_gradients, spacing_constraint_values, spacing_constraint_gradients = con(x)

        # reset xlast
        xlast = x

    end

    # return the desired value
    return boundary_constraint_values[Int(ind)]

end

# set up boundary constraint gradient wrapper function at specific index
# it takes in x as separate elements instead of as a vector (needed for Ipopt)
function boundary_gradient_wrapper_atindex(g, ind, x...)

    # combine the inputs into a single vector
    x = collect(x)

    # include globals for optimization with Ipopt
    global xlast
    global boundary_constraint_gradients

    # check if the constraint calculation is necessary
    if !isequal(x, xlast)

        global boundary_constraint_values
        global spacing_constraint_values
        global spacing_constraint_gradients

        # perform all constraint calculations for this set of design variables
        boundary_constraint_values, boundary_constraint_gradients, spacing_constraint_values, spacing_constraint_gradients = con(x)

        # reset xlast
        xlast = x

    end

    # fill in the provided gradient vector
    g[1] = 0.0 # the first spot is for the index argument, which isn't an actual design variable, so we set the partial with respect to it as zero
    for i = 2:length(g)
        g[i] = boundary_constraint_gradients[Int(ind), i-1]
    end

    # return the gradient vector
    return g

end

# set up spacing constraint wrapper function at specific index
# it takes in x as separate elements instead of as a vector (needed for Ipopt)
function spacing_wrapper_atindex(ind, x...)

    # combine the inputs into a single vector
    x = collect(x)

    # include globals for optimization with Ipopt
    global xlast
    global spacing_constraint_values

    # check if the constraint calculation is necessary
    if !isequal(x, xlast)

        global boundary_constraint_values
        global boundary_constraint_gradients
        global spacing_constraint_gradients

        # perform all constraint calculations for this set of design variables
        boundary_constraint_values, boundary_constraint_gradients, spacing_constraint_values, spacing_constraint_gradients = con(x)

        # reset xlast
        xlast = x

    end

    # return the desired value
    return spacing_constraint_values[Int(ind)]

end

# set up spacing constraint gradient wrapper function at specific index
# it takes in x as separate elements instead of as a vector (needed for Ipopt)
function spacing_gradient_wrapper_atindex(g, ind, x...)

    # combine the inputs into a single vector
    x = collect(x)

    # include globals for optimization with Ipopt
    global xlast
    global spacing_constraint_gradients

    # check if the constraint calculation is necessary
    if !isequal(x, xlast)

        global boundary_constraint_values
        global boundary_constraint_gradients
        global spacing_constraint_values

        # perform all constraint calculations for this set of design variables
        boundary_constraint_values, boundary_constraint_gradients, spacing_constraint_values, spacing_constraint_gradients = con(x)

        # reset xlast
        xlast = x

    end

    # fill in the provided gradient vector
    g[1] = 0.0 # the first spot is for the index argument, which isn't an actual design variable, so we set the partial with respect to it as zero
    for i = 2:length(g)
        g[i] = spacing_constraint_gradients[Int(ind), i-1]
    end

    # return the gradient vector
    return g

end

# set up constraint function wrapper
# when this function is called, all constaints values and gradients are calculated
function con(x)

    # get boundary constraint values and gradients
    boundary_con = boundary_wrapper(x)
    db_dx = ForwardDiff.jacobian(boundary_wrapper, x)

    # get spacing constraint values and gradients
    spacing_con = spacing_wrapper(x)
    ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)

    return boundary_con, db_dx, spacing_con, ds_dx

end

# set up boundary constraint wrapper function
# this function actually calculates the constraint value
function boundary_wrapper(x)

    # include relevant globals
    global boundary_vertices
    global boundary_normals

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get boundary constraint values
    return ff.ray_trace_boundary(boundary_vertices, boundary_normals, turbine_x, turbine_y)

end

# set up spacing constraint wrapper function
# this function actually calculates the constraint value
function spacing_wrapper(x)

    # include relevant globals
    global rotor_diameter

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get spacing constraint values
    return 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x,turbine_y)

end

#---FUNCTIONS FOR OBJECTIVE---

# set up objective wrapper function
# it takes in x as separate elements instead of as a vector (needed for Ipopt)
function aep_wrapper_atindex(x...)

    # combine inputs into a single vector
    x = collect(x)

    # get the objective
    AEP = aep_wrapper(x)[1]

    # return the objective
    return AEP
end

# set up objective gradient wrapper function
# it takes in x as separate elements instead of as a vector (needed for Ipopt)
function aep_gradient_wrapper(g, x...)

    x = collect(x)
    gradient = ForwardDiff.jacobian(aep_wrapper, x)
    for i = 1:size(g)[1]
        g[i] = gradient[i]
    end

    return g
end

# set up objective wrapper function
# this functino actually calculates the objective value
function aep_wrapper(x)

    # include relevant globals
    global turbine_z
    global rotor_diameter
    global hub_height
    global turbine_yaw
    global ct_model
    global generator_efficiency
    global cut_in_speed
    global cut_out_speed
    global rated_speed
    global rated_power
    global windresource
    global power_models
    global model_set
    global rotor_points_y
    global rotor_points_z
    global obj_scale

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set,
                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z)

    # return the objective as a vector
    return [AEP]
end

# import model set with wind farm and related details
include("./model_sets/model_set_6.jl")

# scale objective to be between 0 and 1
obj_scale = 1E-11

# set wind farm boundary parameters
boundary_vertices = ([0 0; 1 0; 1 .75; .75 .75; .75 1; 0 1] .- .5).*500 # Utah-shape boundary
boundary_normals = [0 1.0; -1 0; 0 -1; -1 0; 0 -1; 1 0]

# get the number of vertices of the boundary
n_vertices = length(boundary_vertices)

# initialize variables for constraint wrapper functions
xlast = 0
boundary_constraint_values = 0
boundary_constraint_gradients = 0
spacing_constraint_values = 0
spacing_constraint_gradients = 0
global xlast
global boundary_constraint_values
global boundary_constraint_gradients
global spacing_constraint_values
global spacing_constraint_gradients

# set globals for use in wrapper functions
global model_set
global rotor_points_y
global rotor_points_z
global turbine_z
global ambient_ti
global rotor_diameter
global boundary_vertices
global boundary_normals
global obj_scale

# initialize design variable array
x_initial = [copy(turbine_x);copy(turbine_y)]
n_designvariables = length(x_initial)

# add initial turbine location to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end

# get number of spacing constraints
function numberofspacingconstraints(x)
    n = 0
    for i = 1:x; n += i; end
    return n
end
n_spacingconstraints = numberofspacingconstraints(length(turbine_x)-1)

# set general lower and upper bounds for design variables
# lb = minimum(boundary_vertices)
# ub = maximum(boundary_vertices)

# initialize the model for optimization
model = Model(Ipopt.Optimizer)

# create the design variables
@variable(model, x[i = 1:n_designvariables], start = x_initial[i])

# register the required functions
register(model, :boundary, n_designvariables + 1, boundary_wrapper_atindex, boundary_gradient_wrapper_atindex)
register(model, :spacing, n_designvariables + 1, spacing_wrapper_atindex, spacing_gradient_wrapper_atindex)
register(model, :AEP, n_designvariables, aep_wrapper_atindex, aep_gradient_wrapper)

# set the objective
@NLobjective(model, Max, AEP(x...))

# set the constraints
@NLconstraint(model, boundary_con[i = 1:nturbines], boundary(i, x...) <= 0)
@NLconstraint(model, spacing_con[i = 1:n_spacingconstraints], spacing(i, x...) <= 0)

# run and time optimization
t1 = time()
optimize!(model)
t2 = time()
clkt = t2-t1
xopt = value.(x)
fopt = objective_value(model)
info = (termination_status(model), primal_status(model), dual_status(model))

# print optimization results
println()
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("start objective value: ", aep_wrapper(x_initial)[1])
println("end objective value: ", aep_wrapper(xopt)[1])

# extract final turbine locations
turbine_x = copy(xopt[1:nturbines])
turbine_y = copy(xopt[nturbines+1:end])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--"))
end

# add wind farm boundary to plot
plt.gcf().gca().plot([boundary_vertices[:,1];boundary_vertices[1,1]],[boundary_vertices[:,2];boundary_vertices[1,2]], color="C2")

# set up and show plot
axis("square")
xlim(minimum(boundary_vertices) - (maximum(boundary_vertices)-minimum(boundary_vertices))/5, maximum(boundary_vertices) + (maximum(boundary_vertices)-minimum(boundary_vertices))/5)
ylim(minimum(boundary_vertices) - (maximum(boundary_vertices)-minimum(boundary_vertices))/5, maximum(boundary_vertices) + (maximum(boundary_vertices)-minimum(boundary_vertices))/5)
plt.show()


text:
# CheckerboardFalse
const DQMC_CBFalse = DQMC{M, CheckerboardFalse} where M

"""
	slice_matrix(mc::DQMC_CBFalse, m::Model, slice::Int, power::Float64=1.)

Direct calculation of effective slice matrix, i.e. no checkerboard.
Calculates `Beff(slice) = exp(−1/2∆tauT)exp(−1/2∆tauT)exp(−∆tauV(slice))`.
"""
function slice_matrix(mc::DQMC_CBFalse, m::Model, slice::Int,
					power::Float64=1.)
	eT = mc.s.hopping_matrix_exp
	eTinv = mc.s.hopping_matrix_exp_inv
	eV = mc.s.eV

	interaction_matrix_exp!(mc, m, eV, mc.conf, slice, power)

	if power > 0
		return eT * eT * eV
	else
		return eV * eTinv * eTinv
	end
end
function slice_matrix!(mc::DQMC_CBFalse, m::Model, slice::Int,
					power::Float64=1., result = mc.s.U)
	eT = mc.s.hopping_matrix_exp
	eTinv = mc.s.hopping_matrix_exp_inv
	eV = mc.s.eV

	interaction_matrix_exp!(mc, m, eV, mc.conf, slice, power)

	if power > 0
		# eT * (eT * eV)
		mul!(mc.s.tmp, eT, eV)
		mul!(result, eT, mc.s.tmp)
	else
		# ev * (eTinv * eTinv)
		mul!(mc.s.tmp, eTinv, eTinv)
		mul!(result, eV, mc.s.tmp)
	end
	return result
end

function multiply_slice_matrix_left!(mc::DQMC_CBFalse, m::Model,
								slice::Int, M::AbstractMatrix)
	slice_matrix!(mc, m, slice, 1.0, mc.s.U)
	mul!(mc.s.tmp, mc.s.U, M)
	M .= mc.s.tmp
	nothing
end
function multiply_slice_matrix_right!(mc::DQMC_CBFalse, m::Model,
								slice::Int, M::AbstractMatrix)
	slice_matrix!(mc, m, slice, 1.0, mc.s.U)
	mul!(mc.s.tmp, M, mc.s.U)
	M .= mc.s.tmp
	nothing
end
function multiply_slice_matrix_inv_right!(mc::DQMC_CBFalse, m::Model,
								slice::Int, M::AbstractMatrix)
	slice_matrix!(mc, m, slice, -1.0, mc.s.U)
	mul!(mc.s.tmp, M, mc.s.U)
	M .= mc.s.tmp
	nothing
end
function multiply_slice_matrix_inv_left!(mc::DQMC_CBFalse, m::Model,
								slice::Int, M::AbstractMatrix)
	slice_matrix!(mc, m, slice, -1.0, mc.s.U)
	mul!(mc.s.tmp, mc.s.U, M)
	M .= mc.s.tmp
	nothing
end
function multiply_daggered_slice_matrix_left!(mc::DQMC_CBFalse, m::Model,
								slice::Int, M::AbstractMatrix)
	slice_matrix!(mc, m, slice, 1.0, mc.s.U)
	mul!(mc.s.tmp, adjoint(mc.s.U), M)
	M .= mc.s.tmp
	nothing
end

# CheckerboardTrue
const DQMC_CBTrue = DQMC{M, CheckerboardTrue} where M

function slice_matrix(mc::DQMC_CBTrue, m::Model, slice::Int,
					power::Float64=1.)
    M = Matrix{heltype(mc)}(I, m.flv*m.l.sites, m.flv*m.l.sites)
  	if power > 0
    	multiply_slice_matrix_left!(mc, m, slice, M)
  	else
    	multiply_slice_matrix_inv_left!(mc, m, slice, M)
  	end
  	return M
end
function slice_matrix!(mc::DQMC_CBTrue, m::Model, slice::Int,
					power::Float64=1., M = mc.s.U)
	copyto!(M, I)
  	if power > 0
    	multiply_slice_matrix_left!(mc, m, slice, M)
  	else
    	multiply_slice_matrix_inv_left!(mc, m, slice, M)
  	end
  	return M
end

function multiply_slice_matrix_left!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.s
    interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, 1.)

	mul!(s.tmp, s.eV, M)
	M .= s.tmp
	mul!(s.tmp, s.chkr_mu, M)
  	M .= s.tmp

    @inbounds begin
    	for i in reverse(2:s.n_groups)
			mul!(s.tmp, s.chkr_hop_half[i], M)
			M .= s.tmp
    	end
		mul!(s.tmp, s.chkr_hop[1], M)
		M .= s.tmp
    	for i in 2:s.n_groups
			mul!(s.tmp, s.chkr_hop_half[i], M)
			M .= s.tmp
    	end
  	end
  	nothing
end
function multiply_slice_matrix_right!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.s
  	@inbounds begin
    	for i in reverse(2:s.n_groups)
			mul!(s.tmp, M, s.chkr_hop_half[i])
			M .= s.tmp
    	end
		mul!(s.tmp, M, s.chkr_hop[1])
		M .= s.tmp
    	for i in 2:s.n_groups
			mul!(s.tmp, M, s.chkr_hop_half[i])
			M .= s.tmp
    	end
  	end

  	interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, 1.)
	mul!(s.tmp, M, s.chkr_mu)
	M .= s.tmp
	mul!(s.tmp, M, s.eV)
	M .= s.tmp
  	nothing
end
function multiply_slice_matrix_inv_left!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.s
  	@inbounds begin
    	for i in reverse(2:s.n_groups)
			mul!(s.tmp, s.chkr_hop_half_inv[i], M)
			M .= s.tmp
    	end
		mul!(s.tmp, s.chkr_hop_inv[1], M)
		M .= s.tmp
    	for i in 2:s.n_groups
			mul!(s.tmp, s.chkr_hop_half_inv[i], M)
			M .= s.tmp
    	end
  	end

  	interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, -1.)
	mul!(s.tmp, s.chkr_mu_inv, M)
	M .= s.tmp
	mul!(s.tmp, s.eV, M)
	M .= s.tmp
  	nothing
end
function multiply_slice_matrix_inv_right!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.s
  	interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, -1.)
	mul!(s.tmp, M, s.eV)
	M .= s.tmp
	mul!(s.tmp, M, s.chkr_mu_inv)
	M .= s.tmp

  	@inbounds begin
    	for i in reverse(2:s.n_groups)
			mul!(s.tmp, M, s.chkr_hop_half_inv[i])
			M .= s.tmp
    	end
		mul!(s.tmp, M, s.chkr_hop_inv[1])
		M .= s.tmp
    	for i in 2:s.n_groups
			mul!(s.tmp, M, s.chkr_hop_half_inv[i])
			M .= s.tmp
    	end
  	end
  	nothing
end
function multiply_daggered_slice_matrix_left!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.s

  	@inbounds begin
    	for i in reverse(2:s.n_groups)
			mul!(s.tmp, s.chkr_hop_half_dagger[i], M)
			M .= s.tmp
    	end
		mul!(s.tmp, s.chkr_hop_dagger[1], M)
		M .= s.tmp
    	for i in 2:s.n_groups
			mul!(s.tmp, s.chkr_hop_half_dagger[i], M)
			M .= s.tmp
    	end
  	end

  	interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, 1.)
  	# s.eV == adjoint(s.eV) and s.chkr_mu == adjoint(s.chkr_mu)
	mul!(s.tmp, s.chkr_mu, M)
	M .= s.tmp
	mul!(s.tmp, s.eV, M)
	M .= s.tmp
  	nothing
end

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
using ChainRulesCore
using ChainRulesTestUtils
using FiniteDifferences

const FDM = FiniteDifferences

# Load AD backends
if GROUP == "All" || GROUP == "ForwardDiff"
    @eval using ForwardDiff
end
if GROUP == "All" || GROUP == "Zygote"
    @eval using Zygote
end
if GROUP == "All" || GROUP == "ReverseDiff"
    @eval using ReverseDiff
end
if GROUP == "All" || GROUP == "Tracker"
    @eval using Tracker
end

function test_reverse_mode_ad(f, ȳ, x...; rtol=1e-6, atol=1e-6)
    # Perform a regular forwards-pass.
    y = f(x...)

    # Use finite differencing to compute reverse-mode sensitivities.
    x̄s_fdm = FDM.j′vp(central_fdm(5, 1), f, ȳ, x...)

    if GROUP == "All" || GROUP == "Zygote"
        # Use Zygote to compute reverse-mode sensitivities.
        y_zygote, back_zygote = Zygote.pullback(f, x...)
        x̄s_zygote = back_zygote(ȳ)

        # Check that Zygpte forwards-pass produces the correct answer.
        @test y ≈ y_zygote atol=atol rtol=rtol

        # Check that Zygote reverse-mode sensitivities are correct.
        @test all(zip(x̄s_zygote, x̄s_fdm)) do (x̄_zygote, x̄_fdm)
            return isapprox(x̄_zygote, x̄_fdm; atol=atol, rtol=rtol)
        end
    end

    if GROUP == "All" || GROUP == "ReverseDiff"
        test_rd = length(x) == 1 && y isa Number
        if test_rd
            # Use ReverseDiff to compute reverse-mode sensitivities.
            if x[1] isa Array
                x̄s_rd = similar(x[1])
                tp = ReverseDiff.GradientTape(x -> f(x), x[1])
                ReverseDiff.gradient!(x̄s_rd, tp, x[1])
                x̄s_rd .*= ȳ
                y_rd = ReverseDiff.value(tp.output)
                @assert y_rd isa Number
            else
                x̄s_rd = [x[1]]
                tp = ReverseDiff.GradientTape(x -> f(x[1]), [x[1]])
                ReverseDiff.gradient!(x̄s_rd, tp, [x[1]])
                y_rd = ReverseDiff.value(tp.output)[1]
                x̄s_rd = x̄s_rd[1] * ȳ
                @assert y_rd isa Number
            end

            # Check that ReverseDiff forwards-pass produces the correct answer.
            @test y ≈ y_rd atol=atol rtol=rtol

            # Check that ReverseDiff reverse-mode sensitivities are correct.
            @test x̄s_rd ≈ x̄s_fdm[1] atol=atol rtol=rtol
        end
    end

    if GROUP == "All" || GROUP == "Tracker"
        # Use Tracker to compute reverse-mode sensitivities.
        y_tracker, back_tracker = Tracker.forward(f, x...)
        x̄s_tracker = back_tracker(ȳ)

        # Check that Tracker forwards-pass produces the correct answer.
        @test y ≈ Tracker.data(y_tracker) atol=atol rtol=rtol

        # Check that Tracker reverse-mode sensitivities are correct.
        @test all(zip(x̄s_tracker, x̄s_fdm)) do (x̄_tracker, x̄_fdm)
            return isapprox(Tracker.data(x̄_tracker), x̄_fdm; atol=atol, rtol=rtol)
        end
    end
end

# Define pullback for `to_simplex`
function to_simplex_pullback(ȳ::AbstractArray, y::AbstractArray)
    x̄ = ȳ .* y
    x̄ .= x̄ .- y .* sum(x̄; dims=1)
    return x̄
end
function ChainRulesCore.rrule(::typeof(to_simplex), x::AbstractArray{<:Real})
    y = to_simplex(x)
    pullback(ȳ) = (NoTangent(), to_simplex_pullback(ȳ, y))
    return y, pullback
end

# Define adjoints for ReverseDiff
if GROUP == "All" || GROUP == "ReverseDiff"
    @eval begin
        function to_simplex(x::AbstractArray{<:ReverseDiff.TrackedReal})
            return ReverseDiff.track(to_simplex, x)
        end
        ReverseDiff.@grad function to_simplex(x)
            _x = ReverseDiff.value(x)
            y = to_simplex(_x)
            pullback(ȳ) = (to_simplex_pullback(ȳ, y),)
            return y, pullback
        end
    end
end

# Define adjoints for Tracker
if GROUP == "All" || GROUP == "Tracker"
    @eval begin
        to_posdef(A::Tracker.TrackedMatrix) = Tracker.track(to_posdef, A)
        Tracker.@grad function to_posdef(A::Tracker.TrackedMatrix)
            data_A = Tracker.data(A)
            S = data_A * data_A' + I
            function pullback(∇)
                return ((∇ + ∇') * data_A,)
            end
            return S, pullback
        end

        to_simplex(x::Tracker.TrackedArray) = Tracker.track(to_simplex, x)
        Tracker.@grad function to_simplex(x::Tracker.TrackedArray)
            data_x = Tracker.data(x)
            y = to_simplex(data_x)
            pullback(ȳ) = (to_simplex_pullback(ȳ, y),)
            return y, pullback
        end
    end
end

# Struct of distribution, corresponding parameters, and a sample.
struct DistSpec{D,F,T,X,G,B<:Tuple}
    f::F
    "Distribution parameters."
    θ::T
    "Sample."
    x::X
    "Transformation of sample `x`."
    xtrans::G
    "Broken backends"
    broken::B

    function DistSpec{D}(f::F, θ, x, xtrans::T, broken) where {D,F,T}
        return new{D,F,typeof(θ),typeof(x),T,typeof(broken)}(f, θ, x, xtrans, broken)
    end
end

function DistSpec(d::Distribution, θ, x, xtrans=nothing; broken=())
    return DistSpec{typeof(d)}(d, θ, x, xtrans, broken)
end

function DistSpec(f::F, θ, x, xtrans=nothing; broken=()) where {F}
    D = typeof(f(θ...))
    return DistSpec{D}(f, θ, x, xtrans, broken)
end

dist_type(::DistSpec{D}) where {D} = D

# Auxiliary methods for vectorizing parameters and samples and unflattening them
# similar to `FDM.to_vec`
# However, some implementations in FDM don't work with overload AD such as Tracker,
# ForwardDiff, and ReverseDiff
# Therefore we add a `_to_vec` function

function _to_vec(x::Real)
    function Real_from_vec(v)
        length(v) == 1 || error("vector has incorrect number of elements")
        return first(v)
    end
    return [x], Real_from_vec
end

function _to_vec(x::AbstractArray{<:Real})
    sz = size(x)
    Array_from_vec(v) = reshape(v, sz)
    return vec(x), Array_from_vec
end

function _to_vec(x::Union{Tuple,AbstractVector{<:AbstractArray}})
    x_vecs_and_backs = map(_to_vec, x)
    x_vecs, x_backs = map(first, x_vecs_and_backs), map(last, x_vecs_and_backs)
    lengths = map(length, x_vecs)
    sz = typeof(lengths)(cumsum(collect(lengths)))
    function Tuple_or_Array_of_Array_from_vec(v)
        map(x_backs, lengths, sz) do x_back, l, s
            return x_back(v[(s - l + 1):s])
        end
    end
    return reduce(vcat, x_vecs), Tuple_or_Array_of_Array_from_vec
end

# Functor that fixes non-differentiable location `x` for discrete distributions
struct FixedLocation{X}
    x::X
end
(f::FixedLocation)(args...) = f.x, args

# Functor that transforms differentiable location `x` for continuous distributions
# from unconstrained to constrained space
struct TransformedLocation{F}
    trans::F
end
(f::TransformedLocation)(x, args...) = f.trans(x), args
(f::TransformedLocation{Nothing})(x, args...) = x, args

# Convenience function that returns the correct functor for
# discrete and continuous distributions
make_unpack_x_θ(_, x, ::Type{<:DiscreteDistribution}) = FixedLocation(x)
make_unpack_x_θ(trans, _, ::Type{<:ContinuousDistribution}) = TransformedLocation(trans)

# "Unignore" arguments, i.e., add default arguments if they were ignored
struct Unignore{A}
    args::A
    ignores::BitVector
end

function Unignore(args, ignores::BitVector)
    n = length(args)
    @assert length(ignores) == n
    return Unignore{typeof(args)}(args, ignores)
end

function (f::Unignore)(x...)
    j = Ref(0)
    newx = map(f.args, f.ignores) do argsi, ignoresi
        return if ignoresi
            argsi
        else
            x[(j[] += 1)]
        end
    end

    @assert length(x) == j[] || error("wrong number of arguments")

    return newx
end

# we define the following two functions to be able to tell Zygote that it should not
# compute derivatives for the fields of the functors `unpack_x_θ`
"""
    loglikelihood_parameterized(unpack_x_θ, dist, args...)

Compute the log-likelihood of distribution `dist(θ...)` for `x` where
`x, θ = unpack_x_θ(args...)` are extracted from the arguments `args` with `unpack_x_θ`.

Internally, computations are performed with `loglikelihood`.

See also: [`sum_logpdf_parameterized`](@ref)
"""
function loglikelihood_parameterized(unpack_x_θ, f, args...)
    x, θ = ignore_derivatives(unpack_x_θ)(args...)
    return loglikelihood(f(θ...), x)
end

"""
    sum_logpdf_parameterized(unpack_x_θ, dist, args...)

Compute the log-likelihood of distribution `dist(θ...)` for `x` where
`x, θ = unpack_x_θ(args...)` are extracted from the arguments `args` with `unpack_x_θ`.

Internally, the log pdf of individual data points is computed with `logpdf` which are then
summed up.

See also: [`loglikelihood_parameterized`](@ref)
"""
function sum_logpdf_parameterized(unpack_x_θ, f, args...)
    x, θ = ignore_derivatives(unpack_x_θ)(args...)
    # we use `_logpdf` to be able to handle univariate distributions correctly (see below)
    return sum(_logpdf(f(θ...), x))
end

# Function that computes arrays of `logpdf` values
# `logpdf` does not handle arrays of samples for univariate distributions
_logpdf(d::Distribution, x) = logpdf(d, x)
_logpdf(d::UnivariateDistribution, x::AbstractArray) = logpdf.((d,), x)

# Run AD tests
function test_ad(dist::DistSpec{D}; kwargs...) where {D}
    f = dist.f
    θ = dist.θ
    x = dist.x
    broken = dist.broken

    # combine all arguments
    # point `x` is not differentiable if the distribution is discrete
    args = D <: ContinuousDistribution ? (x, θ...) : θ

    # Create function that splits arguments and transforms location x if needed
    unpack_x_θ = make_unpack_x_θ(dist.xtrans, x, D)

    # short cut: since Zygote does not use special number types with
    # different dispatches etc., it is suffiient to just test derivatives of
    # all differentiable arguments at once
    if GROUP === "All" || GROUP === "Zygote"
        # is Zygote broken?
        zygote_broken = :Zygote in broken

        if zygote_broken
            testset_zygote_broken(dist, unpack_x_θ, args...; kwargs...)
        else
            testset_zygote(dist, unpack_x_θ, args...; kwargs...)
        end
    end

    # Early exit
    GROUP !== "Zygote" || return

    # Define functions for computing the log-likelihood that ignore some arguments
    # (i.e., set them to their default values)
    # This is used to check if we can differentiate with respect to a subset of arguments
    # with ForwardDiff, Tracker, and ReverseDiff
    n = length(args)
    ignores = falses(n)
    unignore = Unignore(args, ignores)
    function loglikelihood_test(x...)
        return sum_logpdf_parameterized(unpack_x_θ, f, unignore(x...)...)
    end
    sum_logpdf_test(x...) = sum_logpdf_parameterized(unpack_x_θ, f, unignore(x...)...)

    # Quick sanity check
    @test loglikelihood_test(args...) ≈ sum_logpdf_test(args...)

    # For all combinations of arguments
    for inds in powerset(1:n, 1, n)
        # Update boolean vector of ignored arguments
        fill!(ignores, true)
        for i in inds
            @inbounds ignores[i] = false
        end

        # Vectorize to-be-differentiated arguments for ForwardDiff, Tracker, and ReverseDiff
        args_vec, args_unflatten = _to_vec(args[inds])
        loglik_test(x) = loglikelihood_test(args_unflatten(x)...)
        logpdf_test(x) = sum_logpdf_test(args_unflatten(x)...)

        @test loglik_test(args_vec) ≈ logpdf_test(args_vec)

        test_ad(loglik_test, args_vec, broken; kwargs...)
        test_ad(logpdf_test, args_vec, broken; kwargs...)
    end

    return
end

function test_ad(f, x, broken = (); rtol = 1e-6, atol = 1e-6)
    finitediff = FDM.grad(central_fdm(5, 1), f, x)[1]

    if GROUP == "All" || GROUP == "Tracker"
        if :Tracker in broken
            @test_broken Tracker.data(Tracker.gradient(f, x)[1]) ≈ finitediff rtol=rtol atol=atol
        else
            @test Tracker.data(Tracker.gradient(f, x)[1]) ≈ finitediff rtol=rtol atol=atol
        end
    end

    if GROUP == "All" || GROUP == "ForwardDiff"
        if :ForwardDiff in broken
            @test_broken ForwardDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        else
            @test ForwardDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        end
    end

    if GROUP == "All" || GROUP == "ReverseDiff"
        if :ReverseDiff in broken
            @test_broken ReverseDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        else
            @test ReverseDiff.gradient(f, x) ≈ finitediff rtol=rtol atol=atol
        end
    end

    return
end

function testset_zygote(distspec, unpack_x_θ, args...; kwargs...)
    f = distspec.f
    θ = distspec.θ
    x = distspec.x

    @testset "Zygote: $(f(θ...)) at x=$x" begin
        @test loglikelihood_parameterized(unpack_x_θ, f, args...) ≈
            sum_logpdf_parameterized(unpack_x_θ, f, args...)

        for l in (loglikelihood_parameterized, sum_logpdf_parameterized)
            # Zygote has type inference problems so we don't check it
            test_rrule(
                Zygote.ZygoteRuleConfig(), l, unpack_x_θ, f, args...;
                rrule_f=rrule_via_ad, check_inferred=false, kwargs...
            )
        end
    end
end

function testset_zygote_broken(args...; kwargs...)
    # don't show test errors - tests are known to be broken :)
    testset = suppress_stdout() do
        testset_zygote(args...; kwargs...)
    end

    # change errors and fails to broken results, and count number of errors and fails
    efs = errors_to_broken!(testset)

    # ensure that passing tests are not marked as broken
    if iszero(efs)
        error("Zygote tests of $(f(θ...)) at x=$x passed unexpectedly, please mark not as broken")
    end

    return testset
end

# `redirect_stdout(f, devnull)` is only available in Julia >= 1.6
function suppress_stdout(f)
    @static if VERSION < v"1.6"
        open((@static Sys.iswindows() ? "NUL" : "/dev/null"), "w") do devnull
            redirect_stdout(f, devnull)
        end
    else
        redirect_stdout(f, devnull)
    end
end

# change test errors and failures to broken results
function errors_to_broken!(ts::Test.DefaultTestSet)
    results = ts.results
    efs = 0
    for i in eachindex(results)
        @inbounds t = results[i]
        if t isa Test.DefaultTestSet
            efs += errors_to_broken!(t)
        elseif t isa Union{Test.Fail, Test.Error}
            efs += 1
            results[i] = Test.Broken(t.test_type, t.orig_expr)
        end
    end
    return efs
end


text:
@testset "Build directly through Dict" begin
    @swagger """
    /test:
        post:
            description: Testing swagger markdown test!
            responses:
                '200':
                    description: Returns a mysterious string test.
    """

    spec = Dict{String, Any}([
        "info" => Dict{String, Any}([
            "version" => "2.0",
            "title" => "Swagger Petstore"
        ]),
        "swagger" => "2.0"
    ])
    spec = build(spec)
    @test haskey(spec, "paths")
end

@testset "Build from json/yml files" begin
    valid_specs = ["spec1_valid_v3.yml", "spec1_valid_v2.yml", "spec1_valid_v3.json", "spec1_valid_v2.json"]
    for filename in valid_specs
        spec = build(joinpath(SPECS_PATH, filename))
        @test haskey(spec, "info")
        @test haskey(spec, "paths")
        @test (haskey(spec, "swagger") || haskey(spec, "openapi"))
    end
end

@testset "Build OpenAPI correctly" begin
    @swagger """
    /test:
        post:
            description: Testing swagger markdown test!
            responses:
                '200':
                    description: Returns a mysterious string test.
    """
    info = Dict{String, Any}()
    info["title"] = "Swagger Petstore"
    info["version"] = "1.0.5"
    openApi = OpenAPI("2.0", info)

    @test openApi.version == "2.0"
    @test openApi.info["title"] == "Swagger Petstore"
    @test openApi.info["version"] == "1.0.5"
end

@testset "Build final spec correctly" begin
    @swagger """
    /test:
        post:
            description: Testing swagger markdown test!
            responses:
                '200':
                    description: Returns a mysterious string test.
    """
    info = Dict{String, Any}()
    info["title"] = "Swagger Petstore"
    info["version"] = "1.0.5"
    openApi = OpenAPI("2.0", info)

    spec = build(openApi)
    cleanup()

    @assert haskey(spec, "swagger")
    @assert haskey(spec, "paths")
    @assert haskey(spec, "info")

    @test spec["swagger"] == "2.0"
    @test spec["info"]["title"] == "Swagger Petstore"
    @test spec["info"]["version"] == "1.0.5"

    paths = spec["paths"]

    TEST = "/test"
    @test haskey(paths, TEST)
    r = paths[TEST]
    @test haskey(r, POST)
    @assert haskey(r[POST], "description")
    @test r[POST]["description"] == "Testing swagger markdown test!"
    @assert haskey(r[POST], "responses")
    @assert haskey(r[POST]["responses"], "200")
    @assert haskey(r[POST]["responses"]["200"], "description")
    @test r[POST]["responses"]["200"]["description"] == "Returns a mysterious string test."
end


text:
module OdeWrapper

using JavaCall

const AbstractFieldIntegrator = @jimport org.hipparchus.ode.AbstractFieldIntegrator
const AbstractIntegrator = @jimport org.hipparchus.ode.AbstractIntegrator
const AbstractParameterizable = @jimport org.hipparchus.ode.AbstractParameterizable
const AdaptiveStepsizeFieldIntegrator = @jimport org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator
const AdaptiveStepsizeIntegrator = @jimport org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator
const BracketedRealFieldUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedRealFieldUnivariateSolver
const BracketedUnivariateSolver = @jimport org.hipparchus.analysis.solvers.BracketedUnivariateSolver
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const ComplexODEConverter = @jimport org.hipparchus.ode.ComplexODEConverter
const ComplexODEState = @jimport org.hipparchus.ode.ComplexODEState
const ComplexODEStateAndDerivative = @jimport org.hipparchus.ode.ComplexODEStateAndDerivative
const ComplexOrdinaryDifferentialEquation = @jimport org.hipparchus.ode.ComplexOrdinaryDifferentialEquation
const ComplexSecondaryODE = @jimport org.hipparchus.ode.ComplexSecondaryODE
const DenseOutputModel = @jimport org.hipparchus.ode.DenseOutputModel
const Enum = @jimport java.lang.Enum
const EquationsMapper = @jimport org.hipparchus.ode.EquationsMapper
const ExpandableODE = @jimport org.hipparchus.ode.ExpandableODE
const Field = @jimport org.hipparchus.Field
const FieldDenseOutputModel = @jimport org.hipparchus.ode.FieldDenseOutputModel
const FieldEquationsMapper = @jimport org.hipparchus.ode.FieldEquationsMapper
const FieldExpandableODE = @jimport org.hipparchus.ode.FieldExpandableODE
const FieldODEEventHandler = @jimport org.hipparchus.ode.events.FieldODEEventHandler
const FieldODEIntegrator = @jimport org.hipparchus.ode.FieldODEIntegrator
const FieldODEState = @jimport org.hipparchus.ode.FieldODEState
const FieldODEStateAndDerivative = @jimport org.hipparchus.ode.FieldODEStateAndDerivative
const FieldODEStateInterpolator = @jimport org.hipparchus.ode.sampling.FieldODEStateInterpolator
const FieldODEStepHandler = @jimport org.hipparchus.ode.sampling.FieldODEStepHandler
const FieldOrdinaryDifferentialEquation = @jimport org.hipparchus.ode.FieldOrdinaryDifferentialEquation
const FieldSecondaryODE = @jimport org.hipparchus.ode.FieldSecondaryODE
const FirstOrderConverter = @jimport org.hipparchus.ode.FirstOrderConverter
const JComplex = @jimport org.hipparchus.complex.Complex
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Locale = @jimport java.util.Locale
const LocalizedODEFormats = @jimport org.hipparchus.ode.LocalizedODEFormats
const MultistepFieldIntegrator = @jimport org.hipparchus.ode.MultistepFieldIntegrator
const MultistepIntegrator = @jimport org.hipparchus.ode.MultistepIntegrator
const NamedParameterJacobianProvider = @jimport org.hipparchus.ode.NamedParameterJacobianProvider
const ODEEventHandler = @jimport org.hipparchus.ode.events.ODEEventHandler
const ODEIntegrator = @jimport org.hipparchus.ode.ODEIntegrator
const ODEJacobiansProvider = @jimport org.hipparchus.ode.ODEJacobiansProvider
const ODEState = @jimport org.hipparchus.ode.ODEState
const ODEStateAndDerivative = @jimport org.hipparchus.ode.ODEStateAndDerivative
const ODEStateInterpolator = @jimport org.hipparchus.ode.sampling.ODEStateInterpolator
const ODEStepHandler = @jimport org.hipparchus.ode.sampling.ODEStepHandler
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const OrdinaryDifferentialEquation = @jimport org.hipparchus.ode.OrdinaryDifferentialEquation
const ParameterConfiguration = @jimport org.hipparchus.ode.ParameterConfiguration
const ParameterJacobianWrapper = @jimport org.hipparchus.ode.ParameterJacobianWrapper
const Parameterizable = @jimport org.hipparchus.ode.Parameterizable
const ParametersController = @jimport org.hipparchus.ode.ParametersController
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SecondOrderODE = @jimport org.hipparchus.ode.SecondOrderODE
const SecondaryODE = @jimport org.hipparchus.ode.SecondaryODE
const VariationalEquation = @jimport org.hipparchus.ode.VariationalEquation

export AbstractFieldIntegrator
export AbstractIntegrator
export AbstractParameterizable
export ComplexODEConverter
export ComplexODEState
export ComplexODEStateAndDerivative
export ComplexOrdinaryDifferentialEquation
export ComplexSecondaryODE
export DenseOutputModel
export EquationsMapper
export ExpandableODE
export FieldDenseOutputModel
export FieldEquationsMapper
export FieldExpandableODE
export FieldODEIntegrator
export FieldODEState
export FieldODEStateAndDerivative
export FieldOrdinaryDifferentialEquation
export FieldSecondaryODE
export FirstOrderConverter
export LocalizedODEFormats
export MultistepFieldIntegrator
export MultistepIntegrator
export NamedParameterJacobianProvider
export ODEIntegrator
export ODEJacobiansProvider
export ODEState
export ODEStateAndDerivative
export OrdinaryDifferentialEquation
export ParameterConfiguration
export ParameterJacobianWrapper
export Parameterizable
export ParametersController
export SecondOrderODE
export SecondaryODE
export VariationalEquation
export add_event_handler
export add_secondary_equations
export add_step_handler
export append
export clear_event_handlers
export clear_step_handlers
export compare_to
export complain_if_not_supported
export compute_derivatives
export compute_main_state_jacobian
export compute_parameter_jacobian
export compute_second_derivatives
export convert_equations
export convert_secondary_equations
export convert_state
export describe_constable
export equals
export extract_equation_data
export extract_main_set_jacobian
export extract_parameter_jacobian
export get_class
export get_complete_derivative
export get_complete_state
export get_complete_state_dimension
export get_current_signed_stepsize
export get_declaring_class
export get_dimension
export get_evaluations
export get_event_handlers
export get_field
export get_final_time
export get_hp
export get_initial_time
export get_interpolated_state
export get_localized_string
export get_mapper
export get_max_evaluations
export get_max_growth
export get_max_step
export get_min_reduction
export get_min_step
export get_n_steps
export get_name
export get_number_of_equations
export get_number_of_secondary_states
export get_ode
export get_parameter
export get_parameter_name
export get_parameters_names
export get_primary
export get_primary_derivative
export get_primary_state
export get_primary_state_dimension
export get_safety
export get_secondary_derivative
export get_secondary_state
export get_secondary_state_dimension
export get_source_string
export get_starter_integrator
export get_step_handlers
export get_step_start
export get_time
export get_total_dimension
export handle_step
export hash_code
export init
export initialize_step
export insert_equation_data
export integrate
export is_supported
export map_state_and_derivative
export name
export notify
export notify_all
export ordinal
export set_initial_main_state_jacobian
export set_initial_parameter_jacobian
export set_initial_step_size
export set_max_evaluations
export set_max_growth
export set_min_reduction
export set_parameter
export set_safety
export set_starter_integrator
export set_step_size_control
export set_up_initial_state
export to_string
export value_of
export values
export wait

include("EventsWrapper/EventsWrapper.jl")
include("NonstiffWrapper/NonstiffWrapper.jl")
include("SamplingWrapper/SamplingWrapper.jl")
include("abstract_field_integrator.jl")
include("abstract_integrator.jl")
include("abstract_parameterizable.jl")
include("complex_ode_converter.jl")
include("complex_ode_state.jl")
include("complex_ode_state_and_derivative.jl")
include("complex_ordinary_differential_equation.jl")
include("complex_secondary_od_e.jl")
include("dense_output_model.jl")
include("equations_mapper.jl")
include("expandable_od_e.jl")
include("field_dense_output_model.jl")
include("field_equations_mapper.jl")
include("field_expandable_od_e.jl")
include("field_ode_integrator.jl")
include("field_ode_state.jl")
include("field_ode_state_and_derivative.jl")
include("field_ordinary_differential_equation.jl")
include("field_secondary_od_e.jl")
include("first_order_converter.jl")
include("localized_ode_formats.jl")
include("multistep_field_integrator.jl")
include("multistep_integrator.jl")
include("named_parameter_jacobian_provider.jl")
include("ode_integrator.jl")
include("ode_jacobians_provider.jl")
include("ode_state.jl")
include("ode_state_and_derivative.jl")
include("ordinary_differential_equation.jl")
include("parameter_configuration.jl")
include("parameter_jacobian_wrapper.jl")
include("parameterizable.jl")
include("parameters_controller.jl")
include("second_order_od_e.jl")
include("secondary_od_e.jl")
include("variational_equation.jl")

end


text:
module TransportBasedInference

using ColorTypes
using DocStringExtensions
using Distributions
using ForwardDiff
using IterativeSolvers
using JLD
using LaTeXStrings
using LinearAlgebra
using LineSearches
using LoopVectorization
using MLDataUtils
using NLsolve
using Optim
using OrdinaryDiffEq
using OrthoMatchingPursuit
using PlotUtils
using Polynomials
using ProgressMeter
using QuadGK
using RecipesBase
using Roots
using SpecialFunctions
using Statistics
using TensorOperations
using ThreadPools

include("tools/get.jl")
include("tools/parallel.jl")
include("tools/ADtools.jl")
include("tools/metric.jl")
include("tools/product.jl")
include("tools/tools.jl")
include("tools/normal.jl")
include("tools/clenshaw_curtis.jl")
# Tools to apply a linear transformation
include("tools/transform.jl")
# Tools for Banana distribution
include("tools/banana.jl")
# Tools for mixture of Gaussian distributions
# include("tools/mixture.jl")
include("tools/view.jl")

# Tools for state-space model
include("statespace/system.jl")

# Tools for data assimilation
include("DA/ensemble.jl")
include("DA/inflation.jl")
include("DA/model.jl")
include("DA/gaspari.jl")
include("DA/seqfilter.jl")
include("DA/seqassim.jl")
include("DA/postprocess.jl")

# Tools for EnKF
include("enkf/senkf.jl")
include("enkf/etkf.jl")

# Tools for lorenz63
include("lorenz63/lorenz63.jl")
include("lorenz63/metric_lorenz.jl")

# Tools for lorenz96
include("lorenz96/lorenz96.jl")

# Rectifier tools
include("hermitemap/rectifier.jl")

# Hermite Polynomials
include("hermitefunction/hermite.jl")
include("hermitefunction/polyhermite.jl")
include("hermitefunction/phypolyhermite.jl")
include("hermitefunction/propolyhermite.jl")

# Hermite Functions
include("hermitefunction/phyhermite.jl")
include("hermitefunction/prohermite.jl")

# Uni and Multi Basis function
include("hermitemap/basis.jl")
include("hermitemap/multibasis.jl")
include("hermitemap/multifunction.jl")
include("hermitemap/expandedfunction.jl")
include("hermitemap/expandedfunction2.jl")
include("hermitemap/storage.jl")

# Integrated positive function
include("hermitemap/integratedfunction.jl")

# ReducedMargin
include("margin/reducedmargin.jl")
include("margin/totalorder.jl")

# KR-rearrangement and TransportMap structure
include("hermitemap/hermitemapcomponent.jl")
include("hermitemap/linhermitemapcomponent.jl")
# Tools for fast inversion
include("hermitemap/inverse.jl")
include("hermitemap/hybridinverse.jl")
# Tools for greedyfit and optimization
include("hermitemap/qr.jl")
include("hermitemap/qraccelerated.jl")
include("hermitemap/precond.jl")
include("hermitemap/greedyfit.jl")
include("hermitemap/optimize.jl")

include("hermitemap/hermitemap.jl")
include("hermitemap/totalordermap.jl")
include("hermitemap/stochmapfilter.jl")

# Tools for radial maps
include("radialmap/function.jl")
include("radialmap/separablecomponent.jl")
include("radialmap/radialmapcomponent.jl")
include("radialmap/sparseradialmapcomponent.jl")
include("radialmap/radialmap.jl")
include("radialmap/weights.jl")
include("radialmap/quantile.jl")
include("radialmap/cost.jl")
include("radialmap/solver.jl")
include("radialmap/optimize.jl")
include("radialmap/iterativeoptimize.jl")
include("radialmap/greedyfit.jl")
include("radialmap/sparseoptimize.jl")
include("radialmap/inverse.jl")
include("radialmap/stochmapfilter.jl")

include("tools/plot_recipes.jl")

end # module