text:
import Dolo

path = Dolo.pkg_path

@testset "testing perfect_foresight" begin

    fn = joinpath(path, "examples", "models", "rbc.yaml")

    model = Dolo.yaml_import(fn)

    # we must define the series for exogenous shocks
    n_e = length(model.symbols[:exogenous])

    T_e = 5
    exo = zeros(T_e, n_e)
    exo[1,:] = [0.00]  # this is used to determine initial steady-state of the model
    exo[2,:] = [0.01]
    exo[3,:] = [0.02]
    exo[4,:] = [0.03]
    exo[5,:] = [0.04]  # this is used to determine final steady-state

    @time df = Dolo.perfect_foresight(model, exo, T=20)

    @test true

end

# Second method: Define a dictionary with shocks
# Compare the results with the previous method

path = Dolo.pkg_path
fn = joinpath(path, "examples", "models", "rbc.yaml");

model = Dolo.yaml_import(fn);

n_e = length(model.symbols[:exogenous])
T_e = 5
exo = zeros(T_e, n_e)
exo[1,:] = [0.00]  # this is used to determine initial steady-state of the model
exo[2,:] = [0.01]
exo[3,:] = [0.02]
exo[4,:] = [0.03]
exo[5,:] = [0.04]
@time df_1matrix = Dolo.perfect_foresight(model, exo, T=20)

exo = Dict(:z =>  [0, 0.01, 0.02, 0.03, 0.04])
@time df_1dict = Dolo.perfect_foresight(model, exo, T=20)

if df_1dict == df_1matrix
  println("df is the same")
end

# Two shocks
fn = joinpath(path, "examples", "models", "rbc_catastrophe.yaml");

model = Dolo.yaml_import(fn)

n_e = length(model.symbols[:exogenous])
T_e = 5
exo = zeros(T_e, n_e)
exo[1,:] = [0.00, 0.00]  # this is used to determine initial steady-state of the model
exo[2,:] = [0.01, 0.00]
exo[3,:] = [0.02, 0.00]
exo[4,:] = [0.03, 0.00]
exo[5,:] = [0.04, 0.00]
@time  df_2matrix = Dolo.perfect_foresight(model, exo, T=20)

exo = Dict(:z =>  [0.00, 0.01, 0.02, 0.03, 0.04], :xi =>  [0.00, 0.00, 0.00, 0.00, 0.00])
@time  df_2dict_a = Dolo.perfect_foresight(model, exo, T=20)

if df_2matrix == df_2dict_a
  println("df is the same")
end

# or

exo = Dict(:z =>  [0.00, 0.01, 0.02, 0.03, 0.04])
@time df_2dict_b = Dolo.perfect_foresight(model, exo, T=20)

if df_2matrix == df_2dict_b
  println("df is the same")
end

#or order changed
exo = Dict( :xi =>  [0.00, 0.00, 0.00, 0.00, 0.00], :z =>  [0.00, 0.01, 0.02, 0.03, 0.04])
@time df_2dict_c = Dolo.perfect_foresight(model, exo, T=20)

if df_2matrix == df_2dict_c
  println("df is the same")
end

#or not defined as Float64
exo = Dict( :xi =>  [0, 0, 0, 0, 0], :z =>  [0.00, 0.01, 0.02, 0.03, 0.04])

@time df_2dict_d = Dolo.perfect_foresight(model, exo, T=20)
if df_2matrix == df_2dict_d
  println("df is the same")
end


text:
typecheckArray()="none"
typecheckArray(x::Array{Void,1})="void"
typecheckArray(x::Array{Bool,1})="boolean"
typecheckArray(x::Array{Int64,1})="int"
typecheckArray(x::Array{Float64,1})="float"
typecheckArray(x::Array{UTF8String,1})="string"
concat(x::Array{UTF8String,1})=reduce((x,y)->x*y,x)
true

ext: jl

max_stars_repo_name: Alfheim/node-julia

max_issues_repo_name: Alfheim/node-julia

max_forks_repo_name: Alfheim/node-julia


text:
using Bosonic
using Test
@testset "Bosonic.jl" begin
    @test basis_m(4,2)[1][1,:][4] == 2
    @test basis_m(4,3)[2][4] == 12
    @test bdb(4,2,1,2)[7,4] == 1
    @test bbd(4,3,2,2)[10,10] == 3
    @test bdb(basis_m(4,2)[1],basis_m(4,2)[2],2,1,2)[9,6] == 1
    @test bbd(basis_m(4,2)[1],basis_m(4,2)[2],2,1,2)[6,9] == 1
    @test bdb(Op_fixed(4,2),1,2)[9,6] == 1
    @test bbd(Op_fixed(4,2),1,2)[9,10] == 1
    @test dim(Op_fixed(5,1)) == 5
    @test nume(Op_fixed(4,3)) == 3
    @test le(Op_fixed(4,3)) == binomial(4+3-1,3)
    @test basis(Op_fixed(4,2))[2,4] == 1
    @test bd(4,1,1)[8,2] == 1
    @test b(4,2,1)[1,7] == 1
    @test bd(4,1,[1,2])[15,2] == 1
    @test b(4,3,[1,2])[1,14] == 1
    @test state_form(bd(4,0,[1,3,2,3]))[22] == 1
end


text:
using Dash
using DashVtk
using DashHtmlComponents
using DashCoreComponents
using PyCall
using VTKDataIO
using VTKDataTypes

dash_vtk = pyimport("dash_vtk")
vtk = pyimport("vtk")
dash_vtk_utils = pyimport("dash_vtk.utils")
pv = pyimport("pyvista")
np = pyimport("numpy")
examples = pyimport("pyvista.examples")

rand(42)
dataset = examples.download_lidar()
fp = _VTKDataTypes(dataset)
#pydataset = PyVTK(dataset)
subset = 0.2
selection = np.random.randint(
    low=1, high=dataset.n_points, size=round(Int,dataset.n_points * subset)
)
#selection = [1 , 3, 5]
n, m = size(dataset.points)
#jlpoints = copy(reshape(vec(dataset.points), m, n)')

points = dataset.points[selection, :]
xyz = points'[:]
elevation = points[:, end]
min_elevation = minimum(elevation)
max_elevation = maximum(elevation)
app = dash()

content = vtk_view([
    vtk_pointcloudrepresentation(
        xyz=xyz,
        scalars=elevation,
        colorDataRange=[min_elevation, max_elevation],
        property=Dict(
            "pointSize" => 2
        ),
    )
])

app.layout = html_div() do
        html_div(
            style=Dict(
                "height" => "calc(100vh - 16px)"
            ),
            children=[(
                html_div(
                    content,
                    style=Dict(
                        "height" => "100%",
                        "width" => "100%",
                    )
                )
            )],
        )
end

run_server(app, "0.0.0.0", debug = true)


text:
using DEBmicroTrait
using CSV, DataFrames, Statistics
using JLD
using HypothesisTests

########################################
# I/O
df_isolates     = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/isolates2traits.csv", DataFrame, missingstring="N/A")
########################################

########################################
Genome_size     = convert(Array{Float64,1}, df_isolates.Genome_size)
V_cell          = DEBmicroTrait.genome_size_to_cell_volume(Genome_size)
rrn_copies      = convert(Array{Float64,1}, df_isolates.rRNA_genes)
Min_gen_time    = df_isolates.Min_gen_time
Gram_stain      = convert(Array{String,1}, df_isolates.gram_stain)
########################################

########################################
gmax            = log(2)./Min_gen_time
V_p             = DEBmicroTrait.cell_volume_to_protein_volume(V_cell)
V_r             = DEBmicroTrait.cell_volume_to_ribosome_volume(V_cell, gmax)
k_E             = DEBmicroTrait.translation_power(V_p, V_r, Min_gen_time)
########################################

########################################
y_EV            = DEBmicroTrait.relative_translation_efficiency_regression(rrn_copies)
########################################

########################################
# I/O
df_out              = DataFrame()
df_out.k_E          = k_E
df_out.y_EV         = y_EV
df_out.response     = df_isolates.Rhizosphere_response

save("/Users/glmarschmann/.julia/dev/DEBmicroTrait/final_manuscript/files/isolates_protein_synthesis.jld", "kE", k_E, "yEV", y_EV, "mingt", Min_gen_time)
CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/isolates_protein_synthesis.csv", df_out)

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

########################################
# statistics
df_p       = filter(x->(x.response.=="positive"), df_out)
df_n       = filter(x->(x.response.=="negative"), df_out)
df_u       = filter(x->(x.response.=="undefined"), df_out)

k_E_kw     = KruskalWallisTest(df_p.k_E, df_n.k_E)
k_E_p      = median(df_p.k_E)
k_E_n      = median(df_n.k_E)

y_EV_kw     = KruskalWallisTest(df_p.y_EV, df_n.y_EV)
y_EV_p      = median(df_p.y_EV)
y_EV_n      = median(df_n.y_EV)
########################################

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
abstract Body

"""
ground is the root body that defines the global coordinate frame.

ground has infinite mass an infinite inertia.

There are 2D & 3D types (ground2D & ground3D)
"""
immutable ground2D <: Body
  position::SymFloat
  velocity::SymFloat
  mass::SymFloat
  inertia::SymFloat
end

immutable ground3D <: Body
  position::SymFloat
  velocity::SymFloat
  mass::SymFloat
  inertia::SymFloat
end

function ground2D()
  g = ground2D(zeros(SymFloat,2),
  zeros(SymFloat,2),
  SymFloat(Sym(Inf),Inf),
  SymFloat(Sym(Inf),Inf))

  return g
end

function ground3D()
  I = Sym("Iground")
  return ground3D(zeros(SymFloat,3),zeros(SymFloat,3),SymFloat(Sym(Inf),Inf),SymFloat(diagM([Inf,Inf,Inf]),diagM([Inf,Inf,Inf])))
end

==(b1::ground2D,b2::ground2D) = true
==(b1::ground3D,b2::ground3D) = true


text:
using Revise, FourierGPE
using Test

## Initialize simulation
L = (12.0,12.0)
N = (256,256)
sim = Sim(L,N)
@unpack_Sim sim

# set simulation parameters
μ = 12.0

# Time dependent potential function (here trivial t dep)
import FourierGPE.V
V(x,y,t) = 0.5*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

# make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim
##

## evolve
@time sol = runsim(sim)

# ground state

ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)
##

##
using VortexDistributions
R(w) = sqrt(2*μ/w^2)
R(1)
rv = 2.
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
ξ0 = healinglength.(0.,0.,μ,g)
ξ = healinglength(rv,0.,μ,g)

pv = PointVortex(rv,0.,1)
vi = ScalarVortex(ξ,pv)

psi = Torus(copy(ψg),x,y)
vortex!(psi,vi)
showpsi(x,y,psi.ψ)

ψi .= psi.ψ
ϕi = kspace(ψi,sim) |> fftshift
kx,ky = K .|> fftshift
x,y = X

## Transform onto polar coordinates in k space
# grids

# make transform to polar grid

function cart2pol(ψi,x,y,kx,ky)
    kmax = maximum(kx)
    nx = length(x)
    nk = nx
    Nk = nk/2 |> Int
    Nθ = 2*nk
    k = LinRange(0,kmax,Nk)
    θ = LinRange(0,2*pi,Nθ+1)[1:Nθ]
    xi = reshape(x,1,1,nx)
    ki = reshape(k,1,1,Nk)
    θi = reshape(θ,1,1,Nθ)
    dθ = diff(θ)[1]
    dk = diff(k)[1]
    dx = diff(x)[1]
    dy = diff(y)[1]

    T = zeros(eltype(ψi),Nk,Nθ,nx)
    S = zeros(eltype(ψi),nx,Nk,Nθ)
    @. T = exp(-im*xi*k*cos(θ'))*dx/sqrt(2*pi)
    @. S = exp(-im*y*k'sin(θi))*dy/sqrt(2*pi)

    ϕp = zeros(eltype(ψi),Nk,Nθ)
    for i in eachindex(x), j in eachindex(y)
        t = @view T[:,:,i]
        s = @view S[j,:,:]
        @. ϕp += t*ψi[i,j]*s
    end
    return ϕp
end

##
@time ϕp = cart2pol(ψi,x,y,kx,ky)
sum(abs2.(ϕp).*k)*dθ*dk
sum(abs2.(ψi))*dx*dy

##
heatmap(log.(eps.() .+abs2.(ϕp)))
Ek = sum(abs2.(ϕp).*k.^3*dθ,dims = 2)
plot(k,log.(Ek))
plot!(k,-k .+5)
vline!([2*pi/ξ0])


text:
# max action level length == 5
import Base: ==
==(f1::Funnel, f2::Funnel) = f1.name == f2.name && f1.end_region[:yr] == f2.end_region[:yr] && f1.end_region[:xr] == f2.end_region[:xr]
function create_funnel_graph(domain_name::String, problem_name::String; max_levels=10)
    graph = Graph()
    domain = load_domain(domain_name)
    problem = load_problem(problem_name)
    init_propositions = get_init_propositions(domain, problem)
    ex_constraints = get_external_constraints(domain, problem)
    graph.props[1] = init_propositions
    graph.initprops = init_propositions
    graph.goalprops = get_goal_propositions(domain, problem)
    graph.num_levels = 1

    for _=1:max_levels
        expand!(graph, domain, problem, ex_constraints)
        if goal_reached(graph, domain, problem) break end
    end
    return graph
end

function expand!(graph::Graph, domain::HPD.Parser.GenericDomain,
    problem::HPD.Parser.GenericProblem, constraints)
    k = graph.num_levels ; print("level $k indexes $(graph.indexes) ")
    if k ≥ 2
        graph.μprops[k] = get_proposition_mutexes(graph, k)
    end
    action_list = []
    graph.props[k+1] = Dict()
    graph.props[k+1][:discrete] = []
    graph.props[k+1][:continuous] = []
    graph.acts[k] = []
    actions = get_applicable_actions(graph, k, constraints, domain, problem)
    repped = false
    for act in actions
        propositions = get_action_propositions(act)
        push!(graph.acts[k], act)
        if !isempty(propositions[:discrete])
            push!(graph.props[k+1][:discrete], propositions[:discrete]...)
        end
        if !isempty(propositions[:continuous])
            if k>1 && !action_in_prev(graph, act, k) && (act.name == :pick || act.name == :place)
                println("Level $(k+1) should have 1 cont ",propositions[:continuous])
                graph.props[k+1][:continuous]=propositions[:continuous]
                repped = true
            elseif !repped && !(propositions[:continuous][1] in graph.props[k][:continuous])
                # println("not here")
                push!(graph.props[k+1][:continuous], propositions[:continuous]...)
            end
        end
        push!(action_list, act)
    end
    graph.acts[k] = collect(Set(graph.acts[k]))
    println(" Num actions: $(length(graph.acts[k]))")
    for prop in graph.props[k][:discrete]
        noop = get_noop_action(prop)
        push!(graph.acts[k], noop)
        push!(graph.props[k+1][:discrete], prop)
    end
    # for cont in graph.props[k][:continuous]
    #     if !(cont in graph.props[k+1][:continuous])
    #         push!(graph.props[k+1][:continuous], cont)
    #     end
    # end
    #graph.acts[k] = collect(Set(graph.acts[k]))
    graph.props[k+1][:discrete] = collect(Set(graph.props[k+1][:discrete]))
    graph.props[k+1][:continuous] = collect(Set(graph.props[k+1][:continuous]))
    graph.μacts[k] = get_action_mutexes(graph, k)
    graph.num_levels +=1
    return graph
end

function action_in_prev(graph::Graph, action::Funnel, level::Int)
    for act in graph.acts[level-1]
        if act.name == action.name
            return true
        end
    end
    return false
end

function terminal_props(domain, problem; init=true)
    termprops=[]
    if init props = collect(initstate(domain, problem).facts) else props = collect(goalstate(domain, problem).facts) end
    for prop in props
        objs = prop.args
        push!(termprops, fill_proposition(prop, objs))
    end
    return termprops

end

function get_init_propositions(domain::HPD.Parser.GenericDomain,
                              problem::HPD.Parser.GenericProblem)
    init_props = Dict()
    init = initstate(domain, problem)
    init_props[:discrete] = terminal_props(domain, problem)
    init_props[:continuous] = [init.values]
    return init_props
end

function get_goal_propositions(domain::HPD.Parser.GenericDomain,
                              problem::HPD.Parser.GenericProblem)
    goal_props = Dict()
    goal = goalstate(domain, problem)
    goal_props[:discrete] = terminal_props(domain, problem;init=false)
    goal_props[:continuous] = goal.values
    return goal_props
end

function goal_reached(graph::Graph, domain::HPD.Parser.GenericDomain,
                                    problem::HPD.Parser.GenericProblem)
    goals = get_goal_propositions(domain, problem)
    index = graph.num_levels-1
    props = graph.props[index]
    μprops = graph.μprops[index]
    goal_found = false
    if issubset(goals[:discrete], props[:discrete])
        goal_found = true
        # for goal_pair in collect(permutations(goals[:discrete], 2))
        #     if goal_pair in μprops
        #         goal_found = false
        #         break
        #     end
        # end
    elseif index > 1 && graph.props[index-1] == props
        graph.leveled = true
    end
    if goal_found
        goal_found = false
        for cont in props[:continuous]
            if in_region(goals[:continuous], cont)
                goal_found = true
                break
            end
        end
    end
    return goal_found
end

function in_region(subregion , region )
    for k in keys(subregion)
        if isa(subregion[k], Float64) && isa(region[k], Float64)
            if subregion[k] != region[k] return false end
        elseif isa(subregion[k], Float64) && isa(region[k], Array) #println("2inhere");
            if !(region[k][1] < subregion[k] < region[k][2]) return false; end
        elseif isa(subregion[k], Array) && isa(region[k], Array)
            if subregion[k][2] < region[k][1] || region[k][2] < subregion[k][1]
                return false
            end
        elseif isa(subregion[k], Array) && isa(region[k], Float64)
            if !(subregion[k][1] < region[k] < subregion[k][2]) return false end
        end
    end
    return true
end

function get_external_constraints(domain::HPD.Parser.GenericDomain,
    problem::HPD.Parser.GenericProblem)
    # ext_consts = [Meta.parse(replace(c.name,"&"=>";")) for c in problem.external_constraints.args]
    ext_consts = [Meta.parse(c.name) for c in problem.external_constraints.args]
    return ext_consts
end

function fill_proposition(proposition, objs)
    if isempty(objs) objs=Term[] end
    prop = Compound(Symbol(proposition.name), objs)
    return prop
end

function get_preconditions(domain, action, vars)
    vardict = Dict()
    [vardict[action.args[i]] = vars[i] for i=1:length(action.args)]
    pos, neg = [], []
    for pre in action.precond.args
        if pre.name == :not
            vs = []
            [push!(vs, vardict[arg]) for arg in pre.args[1].args]
            prop = fill_proposition(pre.args[1], vs)
            push!(neg, prop)
        else
            vs = []
            [push!(vs, vardict[arg]) for arg in pre.args]
            prop = fill_proposition(pre, vs)
            push!(pos, prop)
        end
    end
    return pos, neg
end

function get_effects(domain, action, vars)
    vardict = Dict()
    [vardict[action.args[i]] = vars[i] for i=1:length(action.args)]
    pos, neg = [], []
    for eff in action.effect.args
        if eff.name == :not
            vs = []
            [push!(vs, vardict[arg]) for arg in eff.args[1].args]
            prop = fill_proposition(eff.args[1], vs)
            push!(neg, prop)
        else
            vs = []
            [push!(vs, vardict[arg]) for arg in eff.args]
            prop = fill_proposition(eff, vs)
            push!(pos, prop)
        end
    end
    return pos, neg
end

function get_action_causal_fluents(act::HPD.Parser.GenericAction,
    graph::Graph, domain::HPD.Parser.GenericDomain,
    problem::HPD.Parser.GenericProblem)
    obs = collect(HPD.Parser.get_objects(problem))
    ob = obs[graph.indexes[1]]
    pos_prec, neg_prec = get_preconditions(domain, act, [ob])
    pos_eff, neg_eff = get_effects(domain, act, [ob])
    return pos_prec, neg_prec, pos_eff, neg_eff, [ob]
end

# make more general
# only works with one object
function satisfies_precondition(action, vars, cont_prop)
    if action.name == :move || action.name == :move_holding  return true end
    var = vars[1]
    string_exp = action.cont_precond.args[1].name
    # string_exp = replace(string_exp, "xo"=>string(var.name))
    precond_exp = Meta.parse(string_exp)
    # τ = precond_exp.args[3]
    if action.name == :pick
        xr, yr = cont_prop[:xr], cont_prop[:yr]
        box_vars = ["x"*string(var.name), "y"*string(var.name)]
        xo, yo = cont_prop[Symbol(box_vars[1])], cont_prop[Symbol(box_vars[2])]
        # println(xr," ",yr)
        if isa(xr, Array) && isa(xo, Float64)
            return xr[1]<=xo<=xr[2] && yr[1]<=yo<=yr[2]
        elseif isa(xr, Float64) && isa(xo, Array)
            return xo[1]<=xr<=xo[2] && yo[1]<=yr<=yo[2]
        elseif isa(xr, Array) && isa(xo, Array)
            # println("shouldn't happen")
            return xr[1]<=xo[1] && xo[2]<=xr[2] && yr[1]<=yo[1] && yo[2]<=yr[2]
        else
            fun = @eval (xr, yr, xo, yo) -> $precond_exp
        end
        if Base.invokelatest(fun, xr, yr, xo, yo) return true end
    elseif action.name == :place
        xr, yr = cont_prop[:xr], cont_prop[:yr]
        xg, yg = cont_prop[:xg], cont_prop[:yg]
        if isa(xr, Array) func(xr,yr,xg, yg) = xr[1]<=xg<=xr[2] && yr[1]<=yg<=yr[2]
        else func = @eval (xr, yr, xg, yg) -> $precond_exp end
        if Base.invokelatest(func, xr, yr, xg, yg) return true end
    end
    return false
end

function get_satisfying_index(action, var, cont_props)
    if action.name == :move || action.name == :move_holding return 1 end
    string_exp = action.cont_precond.args[1].name
    precond_exp = Meta.parse(string_exp)
    for (i,cont_prop) in enumerate(cont_props)
        xr, yr = cont_prop[:xr], cont_prop[:yr]
        box_vars = ["x"*string(var.name), "y"*string(var.name)]
        xo, yo = cont_prop[Symbol(box_vars[1])], cont_prop[Symbol(box_vars[2])]
        if isa(xr, Array) && isa(xo, Float64)
            if xr[1]<=xo<=xr[2] && yr[1]<=yo<=yr[2] return i end
        elseif isa(xr, Array) && isa(xo, Array)
            if xr[1]<=xo[1] && xo[2]<=xr[2] && yr[1]<=yo[1] && yo[2]<=yr[2] return i end
        else
            fun = @eval (xr, yr, xg, yg) -> $precond_exp
            if Base.invokelatest(fun, xr, yr, xo, yo) return i end
        end
    end
    return -1
end

function instantiate_action!(action::HPD.Parser.GenericAction, graph::Graph,
            domain::HPD.Parser.GenericDomain, problem::HPD.Parser.GenericProblem, ind::Int, contprop)
    funnel = Funnel(action.name)
    d = 1.0
    funnel.pos_prec, funnel.neg_prec, funnel.pos_eff, funnel.neg_eff, var =
                    get_action_causal_fluents(action, graph, domain, problem)
    if funnel.name == :pick || funnel.name == :place funnel.is_continuous = false end
    if funnel.name == :move funnel.params = [] else funnel.params = var end
    if length(action.cont_precond.args)>0
        push!(funnel.continuous_prec, Meta.parse(action.cont_precond.args[1].name))
    end
    cont_props = graph.props[graph.num_levels][:continuous]
    funnel.cont_ind = ind
    if action.dynamics != true
        xr = contprop[:xr]; yr = contprop[:yr]
        vxmin = graph.props[1][:continuous][1][:vxmin]
        vxmax = graph.props[1][:continuous][1][:vxmax]
        vymin = graph.props[1][:continuous][1][:vymin]
        vymax = graph.props[1][:continuous][1][:vymax]
        if isa(xr, Array)
            xrmin = xr[1] + d*vxmin
            xrmax = xr[2] + d*vxmax
            yrmin = yr[1] + d*vymin
            yrmax = yr[2] + d*vymax
        else
            xrmin = xr + d*vxmin
            xrmax = xr + d*vxmax
            yrmin = yr + d*vymin
            yrmax = yr + d*vymax
        end
        funnel.end_region[:xr] = [xrmin, xrmax]
        funnel.end_region[:yr] = [yrmin, yrmax]
    else
        funnel.end_region[:xr] = contprop[:xr]
        funnel.end_region[:yr] = contprop[:yr]
    end
    #propagate other continuous values
    for key in keys(contprop)
        if !(key in keys(funnel.end_region))
            funnel.end_region[key] = contprop[key]
        end
    end
    if action.cont_effect != true
        if action.name == :pick
            funnel.end_region[:xr] = contprop[Symbol("x"*string(var[1]))]
            funnel.end_region[:yr] = contprop[Symbol("y"*string(var[1]))]
            # graph.props[graph.num_levels][:continuous] = [contprop]
        elseif action.name == :place
            funnel.end_region[:xr] = contprop[:xg]
            funnel.end_region[:yr] = contprop[:yg]
            funnel.end_region[Symbol("x"*string(var[1]))] = contprop[:xg]
            funnel.end_region[Symbol("y"*string(var[1]))] = contprop[:yg]
            # graph.props[graph.num_levels][:continuous] = [contprop]
        end
        if action.name == :move_holding
            funnel.end_region[Symbol("x"*string(var[1]))] = contprop[:xr]
            funnel.end_region[Symbol("y"*string(var[1]))] = contprop[:yr]
        end
    end
    return funnel

end

function working_instantiate_action!(action::HPD.Parser.GenericAction, graph::Graph,
    domain::HPD.Parser.GenericDomain, problem::HPD.Parser.GenericProblem, ind::Int, contprop)
    funnel = Funnel(action.name)
    d = 1.0
    funnel.pos_prec, funnel.neg_prec, funnel.pos_eff, funnel.neg_eff, var =
                get_action_causal_fluents(action, graph, domain, problem)
    if funnel.name == :pick || funnel.name == :place funnel.is_continuous = false end
    if funnel.name == :move funnel.params = [] else funnel.params = var end
    if length(action.cont_precond.args)>0
        push!(funnel.continuous_prec, Meta.parse(action.cont_precond.args[1].name))
    end
    cont_props = graph.props[graph.num_levels][:continuous]
    # ind = get_satisfying_index(action, var[1], cont_props)
    # ind = ind ? ind > 0 : 1
    if action.dynamics != true
        xr = contprop[:xr]; yr = contprop[:yr]
        vxmin = graph.props[1][:continuous][1][:vxmin]
        vxmax = graph.props[1][:continuous][1][:vxmax]
        vymin = graph.props[1][:continuous][1][:vymin]
        vymax = graph.props[1][:continuous][1][:vymax]
        if isa(xr, Array)
            xrmin = xr[1] + d*vxmin
            xrmax = xr[2] + d*vxmax
            yrmin = yr[1] + d*vymin
            yrmax = yr[2] + d*vymax
        else
            xrmin = xr + d*vxmin
            xrmax = xr + d*vxmax
            yrmin = yr + d*vymin
            yrmax = yr + d*vymax
        end
        funnel.end_region[:xr] = [xrmin, xrmax]
        funnel.end_region[:yr] = [yrmin, yrmax]
    else
        funnel.end_region[:xr] = contprop[:xr]
        funnel.end_region[:yr] = contprop[:yr]
    end
    #propagate other continuous values
    for key in keys(contprop)
        if !(key in keys(funnel.end_region))
            funnel.end_region[key] = contprop[key]
        end
    end
    if action.cont_effect != true
    # for ce in action.cont_effect.args
    # cexp = Meta.parse(ce.name)
    # funnel.end_region[Symbol("x"*string(var[1]))] = contprop[:xr]
    # funnel.end_region[Symbol("y"*string(var[1]))] = contprop[:yr]
    # if action.name == :pick
    #     # if graph.num_levels == 2 println(" var name x"*string(var[1]), " ",contprop[Symbol("x"*string(var[1]))] , " ",contprop) end
    #     cont_props[ind][:xr] = contprop[Symbol("x"*string(var[1]))]
    #     cont_props[ind][:yr] = contprop[Symbol("y"*string(var[1]))]
    # elseif action.name == :place
    #     cont_props[ind][:xr] = contprop[:xg]
    #     cont_props[ind][:yr] = contprop[:yg]
    #     cont_props[ind][Symbol("x"*string(var[1]))] = contprop[:xg]
    #     cont_props[ind][Symbol("y"*string(var[1]))] = contprop[:yg]
    # end
        if action.name == :move_holding
            funnel.end_region[Symbol("x"*string(var[1]))] = contprop[:xr]
            funnel.end_region[Symbol("y"*string(var[1]))] = contprop[:yr]
        end
    end
    return funnel

end

function get_applicable_actions(graph::Graph, level::Int, constraints::Vector{Expr}, domain::HPD.Parser.GenericDomain,  problem::HPD.Parser.GenericProblem)
    applicable_actions = []
    actions = collect(values(HPD.Parser.get_actions(domain)))
    props = graph.props[level]
    num_obs = length(problem.objects)
    update_indexes!(graph, level, num_obs)
    for act in actions
        fluents = get_action_causal_fluents(act, graph, domain, problem)
        pos_prec, neg_prec, pos_eff, neg_eff, vs = fluents
        if issubset(pos_prec, props[:discrete])
            for (ind,contprop) in enumerate(props[:continuous])
                pose = get_placement_pose(graph, level, constraints)
                contprop[:xg] = pose[1]; contprop[:yg] = pose[2]
                if satisfies_precondition(act, vs, contprop) #same action can satisfy multiple contprobs
                    instantiated_action = instantiate_action!(act, graph, domain, problem,ind,contprop)
                    if !(instantiated_action in applicable_actions)
                        # if level <=2
                        #     push!(applicable_actions, instantiated_action)
                        # elseif level > 2 && !(instantiated_action in graph.acts[level-2])
                            push!(applicable_actions, instantiated_action)
                        # end
                    end
                    # num_obs = length(problem.objects)
                    # if act.name == :pick && graph.indexes[2]<num_obs graph.indexes[2]+=1 end #++ placepose index
                    # if act.name == :place && graph.indexes[1]<num_obs graph.indexes[1]+=1 end #++ object index
                end
            end
        end
    end
    # println("before ",applicable_actions)
    applicable_actions = prune_funnels(applicable_actions)
    # println("after ",applicable_actions)

    return applicable_actions
end

function int_range(funnel::Union{Funnel, Nothing})
    if isnothing(funnel) return -99 end
    array = funnel.end_region[:yr]
    return array[2]-array[1]
end

function prune_funnels(applicable_actions::Vector{Any})
    pruned_actions = []
    fun_dicts = Dict(:move=>[], :move_holding=>[], :pick=>[], :place=>[])
    for action in applicable_actions push!(fun_dicts[action.name], action) end
    for key in keys(fun_dicts)
        if !isempty(fun_dicts[key])
            max_fun = nothing
            dim1 = []
            for act in fun_dicts[key]
                if isa(act.end_region[:yr], Array)
                    if int_range(act) > int_range(max_fun) max_fun = act  end
                else
                    push!(dim1, act)
                end
            end
            fun_dicts[key] = [max_fun, dim1...]
        end
    end
    pruned_actions = reduce(vcat, collect(values(fun_dicts)))
    filter!(x->x!=nothing, pruned_actions)
    return pruned_actions
end

function update_indexes!(graph::Graph, level::Int, num_obs::Int)
    props = graph.props[level]
    laid_fluent = Compound(:laid, [Const(Symbol("b"*string(graph.indexes[1])))])
    if laid_fluent in props[:discrete] && graph.indexes[1] < num_obs
        graph.indexes[1]+=1
    end
end

function get_proposition_mutexes(graph::Graph, level::Int)
    disc_props = graph.props[level][:discrete]
    μprops=[]
    action_list = graph.acts[level-1]
    for a in action_list
        for (p,q) in collect(permutations(disc_props, 2))
            if p in a.pos_eff && q in a.neg_eff
                push!(μprops, [p,q])
            end
        end

    end
    return μprops
end

function action_pairs_are_mutex(a, b, μprops)
    #1
    if !isempty(intersect(a.neg_eff, union(b.pos_prec, b.pos_eff))) ||
        !isempty(intersect(b.neg_eff, union(a.pos_prec, a.pos_eff)))
         return true
    end
    if !isempty(μprops)
        for μ in μprops
            p = μ[1]
            q = μ[2]
            if (p in a.pos_prec && q in b.pos_prec)
                return true
            end
        end
    end
    return false
end

function get_action_mutexes(graph::Graph, level::Int)
    actions = graph.acts[level]
    μacts = []
    μprops = graph.μprops[level]
    for (a, b) in collect(permutations(actions, 2))
        if action_pairs_are_mutex(a, b, μprops)
            push!(μacts, [a, b])
        end
    end
    return μacts
end

function get_action_propositions(action::Funnel)
    propositions = Dict()
    propositions[:discrete] = action.pos_eff
    propositions[:continuous] = [action.end_region]
    return propositions
end

function get_noop_action(proposition)
    noop = Funnel(:noop)
    noop.pos_prec = [proposition]
    noop.pos_eff = [proposition]
    noop.is_continuous = false
    return noop
end

function get_external_constraint(graph::Graph, level::Int, constraints)
    return constraints #[graph.indexes[1]]
end

#problem specific
function get_placement_pose(graph::Graph, level::Int, constraints)
    ls = get_poses_from_constraints(constraints)
    return ls[graph.indexes[1]]
end

function get_pose_from_fxn!(exp::Expr, ls::Vector{Any})
    xs = 0.0:0.1:5
    ys = 0.0:0.1:5
    fun = @eval (xg, yg) -> $exp
    for x in xs
        for y in ys
            if Base.invokelatest(fun, x, y)
                push!(ls, (x, y))
            end
        end
    end
end

function get_poses_from_constraints(exps::Vector{Expr})
    ls = []
    for exp in exps
        get_pose_from_fxn!(exp, ls)
    end
    return collect(Set(ls))
end

# function get_pose_from_fxn(exp)
#     xs = 0.0:0.1:5
#     ys = 0.0:0.1:5
#     ls = []
#     for xg in xs
#         for yg in ys
#             @eval f(xg, yg)=$exp
#             if f(xg, yg)
#                 push!(ls, (xg, yg))
#             end
#         end
#     end
#     return ls
# end

# box_vars = ["x"*string(var.name), "y"*string(var.name)]
#         for xb = box_vars
#             @eval $xb = cont_prob[Symbol(xb)]
#         end


text:
"""
    write_dicts(deid_dicts)

Writes DeIdDicts structure to file. The dictionaries are written to josn. The
files are written to the  `output_path` specified in the configuration YAML.
"""
function write_dicts(deid::DeIdDicts, logger, outdir)

    currentdate = getcurrentdate()

    idfile = joinpath(outdir, "id_dicts_" * currentdate * ".json")
    Memento.info(logger, "$(Dates.now()) Writing ID dictionaries to $(idfile)")
    write(idfile, JSON.Writer.json(deid.id, 4))

    dateshift_file = joinpath(outdir, "dateshifts_" * currentdate * ".json")
    Memento.info(logger, "$(Dates.now()) Writing dateshift values to $(dateshift_file)")
    write(dateshift_file, JSON.Writer.json(deid.dateshift, 4))

    saltfile = joinpath(outdir, "salts_" * currentdate * ".json")
    Memento.info(logger, "$(Dates.now()) Writing salt values to $(saltfile)")
    write(saltfile, JSON.Writer.json(deid.salt, 4))

    return nothing
end


text:
struct Guardian <: Agent
    id::ID
end

struct Registrator <: Agent
    id::ID
end

struct Member <: Agent
    id::ID
end

struct Braider <: Agent
    id::ID
end

### I could have S argument but I would not need to dispatch on it
### The signature could be a just a String if necessary
struct Transaction{T} # where T <: Union{Agent,Action,Braid}
    document::T
    signature
end

Transaction(x, s::Signature) = Transaction(x, Dict(s))
Transaction(x, s::Signer) = Transaction(x, sign(x, s))

# The next step would be to define expanded Transaction with signers unpacked

### For each pseudonym I could also keep track of anonymity set size
### That could be something to consider in the future

struct State
    guardian::ID
    members::Set{ID}
    pseudonyms::Set{ID}
    registrators::Set{ID}
    braiders::Set{ID}
end

State(guardian::Guardian) = State(guardian.id,Set{ID}(),Set{ID}(),Set{ID}(),Set{ID}())

function Base.push!(state::State, registrator::Registrator)
    push!(state.registrators, registrator.id)
end

function Base.push!(state::State, member::Member)
    push!(state.members, member.id)
    push!(state.pseudonyms, member.id)
end

function Base.push!(state::State, braider::Braider)
    push!(state.braiders, braider.id)
end

function Base.push!(state::State, braid::Braid, signers::Vector{ID})
    for signer in signers
        pop!(state.pseudonyms, signer)
    end

    for id in pseudonyms(braid)
        push!(state.pseudonyms, id)
    end
end

# Validation of power of the signer to add an element of the sort. Also validates that something like dublicate is not beeing added.

function validate(state::State, registrator::Registrator, signer::ID)
    if registrator.id in state.registrators return false end
    return signer == state.guardian
end

function validate(state::State, member::Member, signer::ID)
    if member.id in state.members return false end
    return signer == state.guardian || signer in state.registrators
end

function validate(state::State, braider::Braider, signer::ID)
    if braider.id in state.braiders return false end
    return signer == state.guardian
end

# struct Transactions
#     transactions::Vector{Transaction}
# end

# Tables.is

### I could add a Trigger as a layer for an array

abstract type TransactionLog end

Base.lastindex(t::TransactionLog) = length(t)

### This one is necessary to implement a validation validation check of existing transactions which could happen in the BraidChain constructor. Incidetally that would also is necessary to obtain the state for the BraidChain.
mutable struct TransactionHole <: TransactionLog
    n::Int
end

### Need to think about what to do with hash!
TransactionHole() = TransactionHole(0)

Base.length(t::TransactionHole) = t.n

function Base.push!(tlog::TransactionHole, t::Transaction)
    tlog.n += 1
    return nothing
end

### A simple implementation for a list of transactions stored in memory
struct TransactionVector <: TransactionLog
    transactions::Vector{Transaction}
end

TransactionVector() = TransactionVector(Transaction[])

Base.length(t::TransactionVector) = length(t.transactions)
Base.push!(tlog::TransactionVector, t::Transaction) = push!(tlog.transactions, t)
Base.getindex(t::TransactionVector, i::Int) = t.transactions[i]
Base.getindex(t::TransactionVector, r::UnitRange{Int}) = t.transactions[r]
Base.iterate(t::TransactionVector) = iterate(t.transactions)
Base.iterate(t::TransactionVector, i::Int) = iterate(t.transactions, i)

struct Trigger{T <: TransactionLog} <: TransactionLog
    transactions::T
    trigger::Function
end

Base.length(t::Trigger) = length(t.transactions)

function Base.push!(tlog::Trigger, t::Transaction)
    tlog.trigger(t)
    push!(tlog.transactions, t)
end

Base.getindex(t::Trigger, i::Int) = getindex(t.transactions, i)
Base.getindex(t::Trigger, r::UnitRange{Int}) = getindex(t.transactions, r)

Base.iterate(t::Trigger) = iterate(t.transactions)
Base.iterate(t::Trigger, i::Int) = iterate(t.transactions, i)

# Before one pushes the state it is necessary to validate. On the other hand if validation had already bben done then push! as a raw thing makes sense for recovery of intermidiate states.

struct BraidChain
    transactions::TransactionLog #Vector{Transaction} ### Table
    notary::Notary
    state::State
end

function BraidChain(t::TransactionLog, guardian::Guardian, notary::Notary) # ;validate=true

    tvector = TransactionVector() ### for TransactionHole we need to think how to keep hash.
    bc = BraidChain(tvector, notary, State(guardian))

    for ti in t
        push!(bc, ti)
    end

    @assert length(tvector)==length(t)

    return BraidChain(t, notary, bc.state)
end

state(bc::BraidChain) = bc.state

Tables.istable(::Type{BraidChain}) = true

Tables.isrowtable(::Type{BraidChain}) = true

Tables.rows(log::BraidChain) = log.transactions

Base.eltype(log::BraidChain) = Transaction
Base.length(log::BraidChain) = length(log.transactions)

### I could perhaps add another field for the row!

# struct TransactionRow{T} <: Tables.AbstractRow
#     row::Int
#     source::BraidChain
# end

# Base.iterate(log::BraidChain, st=1) = st > length(log) ? nothing : (TransactionRow(st, log), st + 1)

# getcolumn(m::TransactionRow, ::Type, col::Int, nm::Symbol) =
#     getfield(getfield(m, :source), :matrix)[getfield(m, :row), col]
# getcolumn(m::MatrixRow, i::Int) =
#     getfield(getfield(m, :source), :matrix)[getfield(m, :row), i]
# getcolumn(m::MatrixRow, nm::Symbol) =
#     getfield(getfield(m, :source), :matrix)[getfield(m, :row), getfield(getfield(m, :source), :lookup)[nm]]
# columnnames(m::MatrixRow) = names(getfield(m, :source))

function Base.push!(log::BraidChain, row::Transaction{T}) where T
    signerid = id(row.signature, log.notary)

    @assert validate(log.state, row.document, signerid)
    @assert verify(row.document, row.signature, log.notary)

    push!(log.state, row.document)
    push!(log.transactions, row)
end

function Base.append!(log::BraidChain, rows::Vector{Transaction})
    for row in rows
        push!(log, row)
    end
end

pseudonyms(tlog::BraidChain) = tlog.state.pseudonyms

function pseudonyms(tlog::BraidChain, n::Int)
    @assert length(tlog.transactions) >= n

    state = State(Guardian(tlog.state.guardian))

    for transaction in tlog.transactions[1:n]
        if transaction.document isa Braid
            signers = [id(s, tlog.crypto) for s in row.signature]
            push!(state, transaction.document, signers)
        else
            push!(state, transaction.document)
        end
    end

    return state.pseudonyms
end

function PeaceCypher.hash(tlog::BraidChain, N::Int)
    str = "$(tlog.transactions[1:N])"
    return hash(str, tlog.notary)
end

struct Hash
    N::Int
    hash
end

Hash(tlog::BraidChain, N::Int) = Hash(N,hash(tlog, N))
Hash(tlog::BraidChain) = Hash(tlog, length(tlog))

validate(tlog::BraidChain, th::Hash) = hash(tlog, th.N) == th.hash


text:
# NOTE:
# this is almost just the simple version as QASM 2.0 toplevel
# we should think about how to merge these two

# TODO: improve this to actually validate:
# 1. reset
# 2. intrinsic gate set from given target
# 3. validate the code has no loop
function validate(::QobjQASMTarget, ci::CodeInfo)
    for (v, stmt) in enumerate(ci.code)
        @switch stmt begin
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, _, QuoteNode(locs::Locations))
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), _, _, QuoteNode(locs::Locations), QuoteNode(ctrl::CtrlLocations))
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), _, QuoteNode(locs::Locations))
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :barrier), _, QuoteNode(locs::Locations))
            @case Expr(:invoke, _, GlobalRef(YaoCompiler, :measure_cmp), ::SSAValue, ::Int)
            @case Expr(:invoke, _, GlobalRef(YaoCompiler, :measure_cmp), ::Int, ::SSAValue)
            @case Expr(:new, ::Type{<:NamedTuple}, _...)
            @case Expr(:new, ::Type{<:IntrinsicRoutine}, _...)
            @case ::ReturnNode
            @case ::GotoIfNot
            # return args
            @case Expr(:call, GlobalRef(&Core, :tuple), _...)
            @case _
                error("unsupported operation: $stmt")
        end
    end
end


text:
#!/usr/bin/env julia

using Pkg
Pkg.add("JuliaFormatter")

using JuliaFormatter: format_file
for file in ARGS
    format_file(file, verbose = true)
end

max_stars_repo_name: clemenswasser/benchmarks

max_issues_repo_name: clemenswasser/benchmarks

max_forks_repo_name: clemenswasser/benchmarks


text:
include("../tcm_code/tcm_estimateForest.jl")

function tcm_estimateRanking_L1_ECG(N, assortments, choice_probs, time_limit_overall)

	M = size(assortments,1)

	dim_V = sum(assortments)

	m_dual = Model(solver = GurobiSolver(OutputFlag = 0))
	@variable(m_dual, -1 <= p[1:M, 1:N] <= +1)
	@variable(m_dual, q)

	@objective(m_dual, Max, sum( p[m,i]* choice_probs[m,i] for m in 1:M, i in 1:N ) + q)

    orderings = zeros(Int64,0,N)
	v_by_class = Array{Int64,2}[];

	tol_RC = -1e-4

    # numInitializations_LS = 5

    m_sub = Model(solver = GurobiSolver(OutputFlag = 0))
    @variable(m_sub, x[1:N,1:N], Bin)
    @variable(m_sub, a[1:M,1:N], Bin)

    # VVM - constraint below is not in van Ryzin and Vulcano, but maybe its needed?
    for m in 1:M
        @constraint(m_sub,sum( a[m,p1] for p1 in 1:N) == 1)
    end

    for m in 1:M
        for p1 in 1:N
            @constraint(m_sub,a[m,p1] <= assortments[m,p1])
            for p2 in 1:N
                if ( (p1 != p2) & (assortments[m,p1] == 1) & (assortments[m,p2] == 1) )
                    @constraint(m_sub, a[m,p1] <= x[p1,p2])
                end
           end
        end
    end

    for p1 in 1:N
        for p2 in 1:N
            if p1 != p2
                @constraint(m_sub, x[p1,p2] + x[p2,p1] == 1)
            end
        end
    end

    for p1 in 1:N
        for p2 in 1:N
            for p3 in 1:N
                if ( (p1 != p2) & (p2 != p3) & (p1 != p3) )
                    @constraint(m_sub, x[p1,p2] + x[p2,p3] + x[p3,p1] <= 2)
                end
            end
        end
    end

    function solve_subproblem(u)
        @objective(m_sub, Max, sum(u[m,p] * a[m,p] for m in 1:M, p in 1:N))
        solve(m_sub)

        x_val = getvalue(x)
        objval = getobjectivevalue(m_sub)

        single_ranking = sum(x_val, 2)
        single_ranking = single_ranking[:]
        single_ordering = sortperm(single_ranking, rev = true)

        return single_ordering, objval
    end

	start_time = time()

	solve(m_dual)

	dual_objval = getobjectivevalue(m_dual)

	iter = 0;

    println("Ranking CG iter ", iter, " -- time elapsed = ", 0.0, " objective = ",dual_objval)

	while (time() - start_time < time_limit_overall)
        iter += 1

        p_val = getvalue(p)
        q_val = getvalue(q)

        # single_ordering, sub_objval = tcm_solveRankingSubproblem_L1(N, p_val, assortments, numInitializations_LS)
        single_ordering, sub_objval = solve_subproblem(p_val)

        # @show single_ordering
        # @show typeof(single_ordering)

        RC = -sub_objval - q_val

        println("Ranking CG iter ", iter, " -- time elapsed = ", time() - start_time, " objective = ", dual_objval, " -- RC = ", RC)

        if (RC >= tol_RC)
        	# Stop here -- we are at optimality
        	println("\t RC = ", RC, " -- terminating ... ")
        	break;
        else
        	# Add the ranking to the dual problem, and re-solve the master.
        	orderings = [orderings; single_ordering']

        	v_one_class = tcm_orderingToA( single_ordering, assortments)
        	push!(v_by_class, v_one_class)
        	@constraint(m_dual, sum( v_one_class[m,i] * p[m,i] for m in 1:M, i in 1:N) + q <= 0)

        	solve(m_dual)
        	dual_objval = getobjectivevalue(m_dual)
        end

        # break;
	end

	end_time = time()
	elapsed_time = end_time - start_time

	# Have terminated; solve primal to get what the lambda's are.
	K = size(orderings,1)

	m_primal = Model(solver = GurobiSolver())
	@variable(m_primal, lambda_var[1:K] >= 0)
	@variable(m_primal, eps_plus[1:M,1:N] >= 0)
	@variable(m_primal, eps_minus[1:M,1:N] >= 0)

	for m in 1:M
		for i in 1:N
			@constraint(m_primal, sum( v_by_class[k][m,i] * lambda_var[k] for k in 1:K) + eps_plus[m,i] - eps_minus[m,i] == choice_probs[m,i])
		end
	end
	@constraint(m_primal, sum(lambda_var[k] for k in 1:K) == 1)

	@objective(m_primal, Min, sum( eps_minus[m,i] for m in 1:M, i in 1:N) + sum(eps_plus[m,i] for m in 1:M, i in 1:N))

	solve(m_primal)

	lambda_val = getvalue(lambda_var)

	primal_objval = getobjectivevalue(m_primal)

	return lambda_val, orderings, primal_objval, elapsed_time

end

function tcm_orderingToA(single_ordering, assortments)
	M, N = size(assortments)

	v_one_class = zeros(Float64,size(assortments))
        for m = 1 : M
            for j in 1:N
                if (assortments[m, single_ordering[j]] > 0)
                    v_one_class[m, single_ordering[j]] = 1
                    break;
                end
            end
        end

    return v_one_class
end

# function tcm_solveRankingSubproblem_L1(N, u, assortments, numInitializations_LS)

#     orderings_by_rep = zeros(Int64,numInitializations_LS, N)
#     objval_by_rep = zeros(numInitializations_LS)

#     for rep in 1:numInitializations_LS
#         current_solution = randperm(N)
#         v = tcm_orderingToA(current_solution, assortments)
#         current_obj = dot(u, v)

#         while (true)
#             best_solution = copy(current_solution)
#             best_obj = current_obj

#             for i in 1:(N-1)
#                 for j in (i+1):N
#                     candidate_solution = copy(current_solution)
#                     temp = candidate_solution[i]
#                     candidate_solution[i] = candidate_solution[j]
#                     candidate_solution[j] = temp

#                     v = tcm_orderingToA(candidate_solution, assortments)
#                     candidate_obj = dot(u, v)

#                     if (candidate_obj > best_obj)
#                         best_solution = candidate_solution
#                         best_obj = candidate_obj
#                     end
#                 end
#             end

#             if ( best_obj > current_obj)
#                 current_solution = best_solution
#                 current_obj = best_obj
#             else
#                 break;
#             end
#         end

#         orderings_by_rep[rep,:] = current_solution
#         objval_by_rep[rep] = current_obj
#     end

#     best_rep = indmax(objval_by_rep)
#     best_ordering = orderings_by_rep[best_rep, :]
#     best_objval = objval_by_rep[best_rep]

#     return best_ordering, best_objval
# end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
module PreSumm

include("data/Data.jl")
using PreSumm.Data
export SummaryPair, xsum, cnn_dm, out_dir, snapshot_file, snapshot_files
include("search/Search.jl")
using PreSumm.Search
include("model/Model.jl")
using PreSumm.Model
export Encoder, Decoder, WordPositionEmbed, Generator, TransformersModel, BeamSearch, Translator, BertAbs, TransformerAbs, TransformerAbsTiny, params_count
include("evaluation/Evaluation.jl")
using PreSumm.Evaluation
export mean_rouge_n, mean_rouge_l
include("training/Training.jl")
using PreSumm.Training
export preprocess, loss, logtranslationloss, Warmup, WarmupADAM

end # module


text:
using LinearAlgebra
using GEMPIC
using ParticleInCell
using Test

include("test_poisson_1d.jl")
include("test_landau_damping.jl")
include("test_push.jl")
include("test_fdtd.jl")
include("test_deposition.jl")
include("test_particles.jl")
include("test_particle_group.jl")
include("test_splitting.jl")
include("test_poisson_2d.jl")
include("test_particle_mesh_coupling_spline_2d.jl")
include("test_spline_pp.jl")
include("test_particle_mesh_coupling_spline_1d.jl")
include("test_maxwell_1d_fem.jl")


text:
#
# Here, we aim to slightly more efficiently calculate gradients that may be
# scaled to behave better under optimization algorithms.
#

function scale_gradient!(result::GradResult{V,P,L}, f::F, x::AbstractArray, scale_target,
        cfg::ForwardDiff.GradientConfig{T} = ForwardDiff.GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK, F, V,P,L}
    CHK && ForwardDiff.checktag(T, f, x)
    if ForwardDiff.chunksize(cfg) == length(x)
        ForwardDiff.vector_mode_gradient!(result, f, x, cfg)
    else
        ForwardDiff.chunk_mode_gradient!(result, f, x, cfg)
    end
    ∇ = DiffResults.gradient(result)
    # Need to track down how the buffer becomes NaN in the first place
    # Insert check before multiplication??
    @inbounds for i ∈ P+1:L
        ∇[i] = zero(V)
    end
    # Do not scale up the gradient.
    scale = min(one(V), scale_target / norm(∇))
    result.value *= scale
    @inbounds @simd for l ∈ 1:L
        ∇[l] *= scale
    end
    return scale
end
function scaled_gradient!(result::GradResult, f::F, x::AbstractArray, scale,
        cfg::ForwardDiff.GradientConfig{T} = ForwardDiff.GradientConfig(f, x), ::Val{CHK}=Val{true}()) where {T, CHK, F}
    CHK && ForwardDiff.checktag(T, f, x)
    if ForwardDiff.chunksize(cfg) == length(x)
        scaled_vector_mode_gradient!(result, f, x, scale, cfg)
    else
        scaled_chunk_mode_gradient!(result, f, x, scale, cfg)
    end
    return result
end

function scaled_vector_mode_gradient!(result, f::F, x, scale, cfg::ForwardDiff.GradientConfig{T}) where {T, F}
    ydual = ForwardDiff.vector_mode_dual_eval(f, x, cfg) * scale
    result = ForwardDiff.extract_gradient!(T, result, ydual)
    return result
end

# function extract_gradient!(::Type{T}, result::DiffResult, dual::Dual) where {T}
#     result = DiffResults.value!(result, value(T, dual))
#     result = DiffResults.gradient!(result, partials(T, dual))
#     return result
# end

function chunk_mode_gradient_expr(N, xlen, ydual_expr)

    remainder = xlen % N
    lastchunksize = ifelse(remainder == 0, N, remainder)
    # precalculate loop bounds
    lastchunkindex = xlen - lastchunksize + 1
    middlechunks = 2:div(xlen - lastchunksize, N)
    quote

        # seed work vectors
        xdual = cfg.duals
        seeds = cfg.seeds
        ForwardDiff.seed!(xdual, x)

        # do first chunk manually to calculate output type
        ForwardDiff.seed!(xdual, x, 1, seeds)
        $(ydual_expr)

        ForwardDiff.extract_gradient_chunk!(T, result, ydual, 1, $N)
        ForwardDiff.seed!(xdual, x, 1)

        # do middle chunks
        for c in $middlechunks
            i = ((c - 1) * N + 1)
            ForwardDiff.seed!(xdual, x, i, seeds)
            $(ydual_expr)
            ForwardDiff.extract_gradient_chunk!(T, result, ydual, i, $N)
            ForwardDiff.seed!(xdual, x, i)
        end

        # do final chunk
        ForwardDiff.seed!(xdual, x, $lastchunkindex, seeds, $lastchunksize)
        $(ydual_expr)
        ForwardDiff.extract_gradient_chunk!(T, result, ydual, $lastchunkindex, $lastchunksize)

        # get the value, this is a no-op unless result is a DiffResult
        ForwardDiff.extract_value!(T, result, ydual)

        return result
    end
end

## method defined for better performance of regular chunk mode gradients.
@generated function ForwardDiff.chunk_mode_gradient!(result::GradResult{V,P}, f::F, x,
                                            cfg::ForwardDiff.GradientConfig{T,V,N}) where {F,V,T,N,P}
    chunk_mode_gradient_expr(N, P, :(ydual = f(xdual)))
end
@generated function scaled_chunk_mode_gradient!(result::GradResult{V,P}, f::F, x,
                                            scale, cfg::ForwardDiff.GradientConfig{T,V,N}) where {F,T,V,N,P}
    chunk_mode_gradient_expr(N, P, :(ydual = f(xdual) * scale))
end

#
# function DiffResults.derivative!(r::GradResult, x::AbstractArray)
#     copyto!(r.grad, x)
#     return r
# end
# function DiffResults.derivative!(r::GradResult, x::AbstractArray, ::Type{Val{1}})
#     copyto!(r.grad, x)
#     return r
# end
# function DiffResults.derivative!(r::HessianDiffResult, x::AbstractArray, ::Type{Val{2}})
#     copyto!(r.hess, x)
#     return r
# end
#
# @inline function hessian!(c::HessianConfiguration, x::AbstractArray)
#     ForwardDiff.jacobian!(DiffResults.hessian(c.result), c, DiffResults.gradient(c.result), x, c.jacobian_config, Val{false}())
#     # DiffResults.hessian(c.result)
#     c.result.hess
# end
#
# function jacobian!(result::Union{AbstractArray,DiffResult}, f!, y::AbstractArray, x::AbstractArray, cfg::JacobianConfig{T} = JacobianConfig(f!, y, x), ::Val{CHK}=Val{true}()) where {T,CHK}
#     CHK && checktag(T, f!, x)
#     if chunksize(cfg) == length(x)
#         vector_mode_jacobian!(result, f!, y, x, cfg)
#     else
#         chunk_mode_jacobian!(result, f!, y, x, cfg)
#     end
#     return result
# end
# function vector_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
#     ydual = vector_mode_dual_eval(f!, y, x, cfg)
#     map!(d -> value(T,d), y, ydual)
#     extract_jacobian!(T, result, ydual, N)
#     extract_value!(T, result, y, ydual)
#     return result
# end
#
#
# @generated function chunk_mode_jacobian!(result, f!::F, y, x, cfg::JacobianConfig{T,V,N}) where {F,T,V,N}
#     jacobian_chunk_mode_expr(quote
#                                    ydual, xdual = cfg.duals
#                                    seed!(xdual, x)
#                                end,
#                                :(f!(seed!(ydual, y), xdual)),
#                                :(),
#                                :(extract_value!(T, result, y, ydual))))
# end


text:
# -----------------------------------------------------------------
# Utilities for summarizing a dataset

# TODO: would be useful for Cas to move this in here since he has done a lot
# on that in `IPFitting.jl`. File an issue once this is ready to be
# brought in.


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

module FEMQuad

# Gaussian-Legendre quadratures
include("quaddata.jl")
include("glquad.jl")
# Gaussian-Legendre quadratures in 2d triangles
include("gltri.jl")
# Gaussian-Legendre quadratures in 3d hexahedrons
include("gltet.jl")
# Gaussian-Legendre quadratures in 3d wedges
include("glwed.jl")
# Gaussian-Legendre quadratures in 3d pyramid
include("glpyr.jl")

function get_rule(order::Int, rules::Vararg{Symbol})
    for rule in rules
        if get_order(Val{rule}) >= order
            return rule
        end
    end
    @warn("No accurate rule enough found, picking last.", order, rules)
    return rules[end]
end

function integrate_1d(f::Function, rule::Symbol)
    points = get_quadrature_points(Val{rule})
    result = sum(w*f(ip) for (w, ip) in points)
    return result
end

function integrate_2d(f::Function, rule::Symbol)
    points = get_quadrature_points(Val{rule})
    result = sum(w*f(ip) for (w, ip) in points)
    return result
end

function integrate_3d(f::Function, rule::Symbol)
    points = get_quadrature_points(Val{rule})
    result = sum(w*f(ip) for (w, ip) in points)
    return result
end

export get_quadrature_points

end


text:
using Publish
using BitSAD
using Artifacts, LazyArtifacts

# override default theme
Publish.Themes.default() = artifact"flux-theme"

p = Publish.Project(BitSAD)

# serve documentation
serve(BitSAD)


text:
import LinearAlgebra.cross

struct precipMesh
  sMesh::ScalarMesh
  fMesh::FluidMesh
end

struct precipSolution
  ψA::ScalarSolution
  ψB::ScalarSolution
  ψC::ScalarSolution
  θs::ScalarSolution
  θf::ScalarSolution
  fluid::FluidSolution      # velocity and pressure solution
end

struct precipProblem
  ψA::Problem
  ψB::Problem
  ψC::Problem
  fluid::Problem
end

function initializeSolution(p::precipParam)
  sMesh,fMesh = loadMesh(p)
  N = length(sMesh.xy); Np = length(fMesh.xyp)
  zs = zeros(N); zq = zeros(N); zp = zeros(Np)
  ψA = ScalarSolution(copy(zs))
  ψB = ScalarSolution(copy(zs))
  ψC = ScalarSolution(copy(zs))
  θs = ScalarSolution(copy(zs))
  θf = ScalarSolution(ones(N))
  fluidqp = FluidSolution(zq,copy(zq),zp)

  if p.runID == 6
    for i=1:N
      x = sMesh.xy[i].x; y = sMesh.xy[i].y
      ψA.u[i] = 5e-4*diffusionAndReactionIC(x,y,:ψA)
      ψB.u[i] = 5e-4*diffusionAndReactionIC(x,y,:ψB)
    end
  end

  return precipSolution(ψA,ψB,ψC,θs,θf,fluidqp)
end

function diffusionAndReactionIC(x,y,symbol)
  if symbol == :ψA
    #if inLeftArea(x,y)
    #  return 1.0
    #else
    #  return 0.0
    #end
    return abs(x)
  elseif symbol == :ψB
    #if inRightArea(x,y)
    #  return 1.0
    #else
    #  return 0.0
    #end
    return abs(y)
  end
  error("symbol not recognized")
end
function inLeftArea(x,y)
  if x < -0.5 && x > -1.5 && y < 0.5 && y > -0.5
    return true
  else
    return false
  end
end
function inRightArea(x,y)
  if x < 1.5 && x > 0.5 && y < 0.5 && y > -0.5
    return true
  else
    return false
  end
end

function randomizeScalarSolutions!(soln::precipSolution)
  N = length(soln.ψA.u)
  soln.ψA.u = 0.5*1e-2*rand(N)
  soln.ψB.u = 0.5*1e-2*rand(N)
  soln.ψC.u = zeros(N)
  soln.θs.u = zeros(N)
  soln.θf.u = 1 .- soln.θs.u
  nothing
end

function defineGMSHChemicalA(mesh::ScalarMesh)
  # define boundary conditions
  dNodes = Dirichlet(:inflow)
  nNodes = Neumann(:outflow,:wall)
  dBCf = Dirichlet((x,y) -> scalarBCf_ydom(x,y,:lower))
  nBCf = Neumann((x,y) -> 0.0)
  Nodes = [dNodes,nNodes]; bcfun = [dBCf,nBCf]

  # define problems
  chemAProb = Problem(mesh,Nodes,bcfun,:AdvDiff2D)
end

function defineGMSHChemicalB(mesh::ScalarMesh)
  # define boundary conditions
  dNodes = Dirichlet(:inflow)
  nNodes = Neumann(:outflow,:wall)
  dBCf = Dirichlet((x,y) -> scalarBCf_ydom(x,y,:upper))
  nBCf = Neumann((x,y) -> 0.0)
  Nodes = [dNodes,nNodes]; bcfun = [dBCf,nBCf]

  # define problems
  chemBProb = Problem(mesh,Nodes,bcfun,:AdvDiff2D)
end

function defineGMSHChemicalC(mesh::ScalarMesh)
  # define boundary conditions
  nNodes = Neumann(:inflow,:outflow,:wall)
  nBCf = Neumann((x,y) -> 0.0)
  Nodes = [nNodes]; bcfun = [nBCf]

  # define problems
  chemCProb = Problem(mesh,Nodes,bcfun,:AdvDiff2D)
end

function defineGMSHMPBrinkman(mesh::FluidMesh)
  dNodes = Dirichlet(:inflow,:wall)
  nNodes = Neumann(:outflow)
  dUBCf = Dirichlet(uBCf_ydom)
  dVBCf = Dirichlet(vBCf_ydom)
  nBCf = Neumann((x,y) -> 0.0)
  Nodes = [dNodes,nNodes]; bcfun = [dUBCf,dVBCf,nBCf]

  # define problems
  fluidProb = Problem(mesh,Nodes,bcfun,:BrinkmanMP2D)
end

function defineSquareChemical(mesh::ScalarMesh)
  # define boundary conditions
  nNodes = Neumann(:left,:right,:top,:bottom)
  nBCf = Neumann((x,y) -> 0.0)
  Nodes = [nNodes]; bcfun = [nBCf]

  # define problems
  chemAProb = Problem(mesh,Nodes,bcfun,:Poisson2D)
end

function defineSquareMPBrinkman(mesh::FluidMesh)
  dNodes = Dirichlet(:left,:right,:bottom,:top)
  dUBCf = Dirichlet((x,y) -> 0.0)
  dVBCf = Dirichlet((x,y) -> 0.0)
  Nodes = [dNodes]; bcfun = [dUBCf,dVBCf]

  # define problems
  fluidProb = Problem(mesh,Nodes,bcfun,:BrinkmanMP2D)
end

function defineProblems(p::precipParam,mesh::precipMesh)
  if p.runID == 1
    ψAprob = defineGMSHChemicalA(mesh.sMesh)
    ψBprob = defineGMSHChemicalB(mesh.sMesh)
    ψCprob = defineGMSHChemicalC(mesh.sMesh)
    fluidprob = defineGMSHMPBrinkman(mesh.fMesh)
  elseif p.runID == 2
    ψAprob = defineGMSHChemicalA(mesh.sMesh)
    ψBprob = defineGMSHChemicalB(mesh.sMesh)
    ψCprob = defineGMSHChemicalC(mesh.sMesh)
    fluidprob = defineGMSHMPBrinkman(mesh.fMesh)
  elseif p.runID == 3
    error("undefine problem. see problem_structs.jl")
  elseif p.runID == 4
    error("undefine problem. see problem_structs.jl")
  elseif p.runID == 5
    error("undefine problem. see problem_structs.jl")
  elseif p.runID == 6
    ψAprob = defineSquareChemical(mesh.sMesh)
    ψBprob = defineSquareChemical(mesh.sMesh)
    ψCprob = defineSquareChemical(mesh.sMesh)
    fluidprob = defineSquareMPBrinkman(mesh.fMesh)
  else
    error("you shouldn't get here")
  end

  return precipProblem(ψAprob,ψBprob,ψCprob,fluidprob)
end

function uBCf_ydom(x,y)
  # define relavant parameters for mesh
  θ = 0.7227342478134156
  α = pi/2 - θ
  sα = sin(α)
  sθ = sin(θ)
  cα = cos(α)
  cθ = cos(θ)
  IL = 1.0
  IW = 0.5
  OW = 1.5*IW
  mp = sqrt(IW^2-(0.5*OW)^2)

  # define (x,y) nodes defining inflow boundary
  x1  = -mp-IL*sα;  x2  = -IL*sα
  y1L = -IL*cα;  y2L = -IL*cα - OW/2
  y1U =  IL*cα;  y2U =  IL*cα + OW/2

  # determine if (x,y) sits on interface
  pTest = [x,y,0.0]
  p1L   = [x1,y1L,0.0]; p2L = [x2,y2L,0.0]
  p1U   = [x1,y1U,0.0]; p2U = [x2,y2U,0.0]

  v1L = pTest - p1L;  v2L = p2L - p1L
  v1U = pTest - p1U;  v2U = p2U - p1U

  cL = cross(v1L,v2L)[3]; cU = cross(v1U,v2U)[3]

  sectionName = :temp

  if abs(cL)<1e-8
    sectionName = :lower
  elseif abs(cU)<1e-8
    sectionName = :upper
  else
    sectionName = :wall
  end

  # if on inflow transformation = (x,y) -> (c*(-W+1)-SL, c*(-W+1))
  if sectionName == :lower
    xstart = x1; xend   = x2
    s      = linearInterp(x,xstart,xend,0,1)

    return 4*s*(1-s)*sα

  elseif sectionName == :upper  # top inflow
    xstart = x1; xend   = x2
    s      = linearInterp(x,xstart,xend,0,1)

    return 4*s*(1-s)*sα

  elseif sectionName == :wall
    return 0.0
  end
end

function vBCf_ydom(x,y)
  # define relavant parameters for mesh
  θ = 0.7227342478134156
  α = pi/2 - θ
  sα = sin(α)
  sθ = sin(θ)
  cα = cos(α)
  cθ = cos(θ)
  IL = 1.0
  IW = 0.5
  OW = 1.5*IW
  mp = sqrt(IW^2-(0.5*OW)^2)

  # define (x,y) nodes defining inflow boundary
  x1  = -mp-IL*sα;  x2  = -IL*sα
  y1L = -IL*cα;  y2L = -IL*cα - OW/2
  y1U =  IL*cα;  y2U =  IL*cα + OW/2

  # determine if (x,y) sits on interface
  pTest = [x,y,0.0]
  p1L   = [x1,y1L,0.0]; p2L = [x2,y2L,0.0]
  p1U   = [x1,y1U,0.0]; p2U = [x2,y2U,0.0]

  v1L = pTest - p1L;  v2L = p2L - p1L
  v1U = pTest - p1U;  v2U = p2U - p1U

  cL = cross(v1L,v2L)[3]; cU = cross(v1U,v2U)[3]

  sectionName = :temp

  if abs(cL)<1e-8
    sectionName = :lower
  elseif abs(cU)<1e-8
    sectionName = :upper
  else
    sectionName = :wall
  end

  # if on inflow transformation = (x,y) -> (c*(-W+1)-SL, c*(-W+1))
  if sectionName == :lower          # bottom inflow
    xstart = x1; xend = x2
    s      = linearInterp(x,xstart,xend,0,1)

    return 4*s*(1-s)*cα
  elseif sectionName == :upper      # top inflow
    xstart = x1; xend = x2
    s      = linearInterp(x,xstart,xend,0,1)

    return -4*s*(1-s)*cα
  elseif sectionName == :wall
    return 0.0
  end
end

function scalarBCf_ydom(x,y,ENTRANCE::Symbol)
  # define relavant parameters for mesh
  θ = 0.7227342478134156
  α = pi/2 - θ
  sα = sin(α)
  sθ = sin(θ)
  cα = cos(α)
  cθ = cos(θ)
  IL = 1.0
  IW = 0.5
  OW = 1.5*IW
  mp = sqrt(IW^2-(0.5*OW)^2)

  # define (x,y) nodes defining inflow boundary
  x1  = -mp-IL*sα;  x2  = -IL*sα
  y1L = -IL*cα;  y2L = -IL*cα - OW/2
  y1U =  IL*cα;  y2U =  IL*cα + OW/2

  # determine if (x,y) sits on interface
  pTest = [x,y,0.0]
  p1L   = [x1,y1L,0.0]; p2L = [x2,y2L,0.0]
  p1U   = [x1,y1U,0.0]; p2U = [x2,y2U,0.0]

  v1L = pTest - p1L;  v2L = p2L - p1L
  v1U = pTest - p1U;  v2U = p2U - p1U

  cL = cross(v1L,v2L)[3]; cU = cross(v1U,v2U)[3]

  sectionName = :temp

  if abs(cL)<1e-8
    sectionName = :lower
  elseif abs(cU)<1e-8
    sectionName = :upper
  else
    sectionName = :wall
  end

  # if on inflow transformation = (x,y) -> (c*(-W+1)-SL, c*(-W+1))
  if sectionName == ENTRANCE
    return 0.005
  else
    return 0.0
  end
end

"""
  linearInterp()

linear interpolates the point x0 which lies in (x1,x2) onto the domain (s1,s2)
"""
function linearInterp(x0,x1,x2,s1,s2)
  A = [x1 1; x2 1]; rhs=[s1,s2]
  mbtemp = A\rhs
  m = mbtemp[1]
  b = mbtemp[2]

  s0 = m*x0 + b

  return s0
end


text:
using MultivariateStats
using Distances
using Catch22
using ManifoldLearning
using LinearAlgebra
using LowRankModels
using SparseArrays

# Function wrapping projectors
function project(F::AbstractArray, projectionFunc, args...; kwargs...)
    M = projectionFunc(F, args...; kwargs...)
end
function project(F::AbstractFeatureArray, projectionFunc, args...; kwargs...)
    M = projectionFunc(Array(F), args...; kwargs...)
end
export project

"""
PCA
"""
function principalcomponents(F::AbstractArray; pratio=1.0, kwargs...)
    M = MultivariateStats.fit(MultivariateStats.PCA, F; pratio=pratio, kwargs...)
end
principalcomponents(F::AbstractFeatureArray, args...; kwargs...) = principalcomponents(Array(F), args...; kwargs...)
export principalcomponents

function embed(M::MultivariateStats.PCA, F::AbstractArray, PCs::Union{Int, Vector{Int64}, UnitRange}=1:length(M.prinvars))
    P = MultivariateStats.projection(M)
    P = P[:, PCs]
    D = (P'*(F .- mean(M)))
    if size(D, 1) == 1
        D = D[:]
    end
    return D
end

export PCA

function embed(M::MultivariateStats.PCA, F::AbstractFeatureArray, PCs::Union{Int, Vector{Int64}, UnitRange}=1:length(M.prinvars); kwargs...)
    D = embed(M, Array(F), PCs; kwargs...)
    Catch22.featureMatrix(D, [Symbol("PC$x") for x ∈ 1:size(D, 1)])
end
export embed

"""
Explained Variance
"""
function explainedVariance(M::MultivariateStats.PCA)
    ev = principalvars(M)
    ev = cumsum(unitL1(ev))
end
export explainedVariance

# add method to explainedVariance for new data

function explainedVariance(F::AbstractArray)
# The varainces of each feature as a proportion of the total variance
    Σ² = StatsBase.cov(F, dims=2)
    ev = diag(Σ²)./sum(diag(Σ²))
end

export principalvars

"""
Weights on features
"""
function featureweights(model::MultivariateStats.PCA, pc=1:outdim(model))
    P = projection(model)[:, pc]
end
function featureweights(I::Inference, pc::Int=1)
    P = projection(I.model)
    Catch22.featureVector(P[:, pc], val(dims(I.F̂, :feature)))
end
export featureweights

"""
Isomap
"""
function isomap(F::AbstractArray; maxoutdim=max(size(F, 1)÷2, 10), kwargs...)
    # Strange bug here; cant have the maxoutdim too high otherwise:
    #ERROR: DomainError with -8.109628524266554e-10:
    #sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
    # Watch out for this in the future
    M = MultivariateStats.fit(ManifoldLearning.Isomap, Array(F); maxoutdim=maxoutdim, kwargs...)
end
export isomap

function embed(M::ManifoldLearning.Isomap, F::AbstractArray, PCs::Union{Int, Vector{Int64}, UnitRange}=1:ManifoldLearning.outdim(M))
    D = ManifoldLearning.transform(M, F)
    D = D[PCs, :]
    if size(D, 1) == 1
        D = D[:]
    end
    return D
end
export Isomap
# If you just want the original data transformed, this is slightly faster:
function embed(M::ManifoldLearning.Isomap, PCs::Union{Int, Vector{Int64}, UnitRange}=1:ManifoldLearning.outdim(M))
    D = ManifoldLearning.transform(M)
    D = D[PCs, :]
    if size(D, 1) == 1
        D = D[:]
    end
    return D
end
export embed

function residualVariance(M, F::AbstractArray, nPCs::Int64)
    D = embed(M, Array(F), 1:nPCs) # Model self-embedding can be run more quickly...
    if D isa Vector
        D = D'
    end
    dD = pairwise(Euclidean(), D, D, dims=2)
    dF = pairwise(Euclidean(), F, F, dims=2)
    residualVariance(dF, dD)
end
function residualVariance(M, F::AbstractArray)
    𝛔² = residualVariance.((M,), (F,), 1:outdim(M))
end
function residualVariance(dF::Array, dD::Array{Float64, 2})
    # Calculate a distance matrix for F and D, and correlate them
    σ² = 1 - (cor(dF[:], dD[:]))^2
    # Sqrt of 1 minus the correlation between distances in feature space and the low dimensional space, squared
end
export residualVariance

function robustprincipalcomponents(F::AbstractArray, npcs=size(F, 1)::Int, scale::Float64=1.0; abs_tol=1e-4^3, rel_tol=1e-4^3, max_iter=10000, kwargs...)
	loss = HuberLoss()
	r = ZeroReg()
    # Start with PCA
    m = GLRM(F, QuadLoss(), ZeroReg(), ZeroReg(), npcs; kwargs...)
    fit!(m)
    # Then do the other thing
    m = GLRM(A, loss, r, r, k; kwargs...)
    fit!(m, ProxGradParams(;abs_tol, rel_tol, max_iter))
	return m
end
MultivariateStats.projection(M::LowRankModels.GLRM) = inv(M.X)
function embed(M::LowRankModels.GLRM, F::AbstractArray, PCs::Union{Int, Vector{Int64}, UnitRange}=1:M.k)
    @assert M.k == size(F,  1) # For now, to invert X it must be square whcih means keeping all PCs
    P = projection(M)
    P = P[:, PCs]
    D = P'*F
    if size(D, 1) == 1
        D = D[:]
    end
    return D
end
export robustprincipalcomponents

nanmedian = x -> all(isnan.(x)) ? NaN : median(x[.!isnan.(x)])
naniqr = x -> all(isnan.(x)) ? NaN : iqr(x[.!isnan.(x)])

function nanprincipalcomponents(F::AbstractArray, args...; kwargs...)
    # Convert F to a sparse array, sending NaNs to sparse elements
    init = deepcopy(F) |> Array # Chunky?
    for i ∈ 1:size(init, 1)
        repnans = all(isnan.(init[i, :])) ? randn(sum(isnan.(init[i, :]))).*eps() : nanmedian(init[i, :][.!isnan.(init[i, :])]) # If all nans, just use some tiny noise.
        init[i, isnan.(init[i, :])] .= repnans
    end
    A = SparseMatrixCSC(F)
    SparseArrays.fkeep!(A, (i,j,x) -> !isnan(x))
    m = _nanprincipalcomponents(A, init, args...; kwargs...)
end

function _nanprincipalcomponents(A::AbstractArray, init, npcs=size(A, 1)::Int, scale::Float64=1.0; abs_tol=1e-4^3, rel_tol=1e-4^3, max_iter=20000, kwargs...)
    m = GLRM(A, QuadLoss(), ZeroReg(), ZeroReg(), npcs; kwargs...)
    U, D, V = svd(init)
    D = D |> Diagonal
    m.X, m.Y = collect((U*sqrt(D))'), sqrt(D)*V';
    fit!(m, ProxGradParams(;abs_tol, rel_tol, max_iter))
    return m
end
export nanprincipalcomponents

function outlierprincipalcomponents(F::AbstractArray, args...; kwargs...)
    # * The goal is the remove any outliers (say, greater than 10 iqrs from the median) and then do PCA
    A = deepcopy(F)
    centres = mapslices(nanmedian, F, dims=2)
    scales = mapslices(naniqr, F, dims=2)
    for i ∈ 1:size(A, 1)
        A[i, abs.(A[i, :] .- centres[i]) .>= 5*scales[i]] .= NaN
    end
    return nanprincipalcomponents(A, args...; kwargs...)
end
export outlierprincipalcomponents

"""
Only keep the PC's explaining more than `thresh` proportion of variance.
"""
function significantprincipalcomponents(F::AbstractArray; thresh=1/size(F, 1), kwargs...)
    M = principalcomponents(F; kwargs...)
    v = principalvars(M)./sum(principalvars(M))
    npcs = sum(v .> thresh)
    M = principalcomponents(F; maxoutdim=npcs, kwargs...)
end
significantprincipalcomponents(F::AbstractFeatureArray, args...; kwargs...) = significantprincipalcomponents(Array(F), args...; kwargs...)
export significantprincipalcomponents


text:
using Distributions
using Parameters
using LightGraphs
using TransformVariables
using LogDensityProblems
using DynamicHMC
using ForwardDiff
using BenchmarkTools
using Plots

include("model.jl")
include("generate.jl")

S = 10

nicheweb = niche_network(S, 0.1)

heatmap(nicheweb.fw, legend = false)

m = model(nicheweb.fw)

heatmap(m.foodweb, legend = false)

θ = (n = nicheweb.n[m.reorder],
     r = nicheweb.r[m.reorder[(m.basal + 1):end]],
     c = nicheweb.c[m.reorder[(m.basal + 1):end]])

m(θ)

toθ = as((n = as(Array, as𝕀, m.S),
          r = as(Array, as𝕀, m.nonbasal),
          c = as(Array, as𝕀, m.nonbasal)))

lp_t = TransformedLogDensity(toθ, m)
lp∇ = ForwardDiffLogDensity(lp_t, chunk = ForwardDiff.Chunk{24}())

inits = inverse(toθ, θ)
logdensity(LogDensityProblems.ValueGradient, lp∇, inits)

chain, smplr = NUTS_init_tune_mcmc(lp∇, 500, max_depth = 10)

NUTS_statistics(chain)

posterior = transform.(Ref(lp∇.transformation), get_position.(chain))

n = mean(first, posterior)

foo = zip(posterior...) |> collect

scatter(θ.n, n)

max_stars_repo_name: clint-leach/niche-model

max_issues_repo_name: clint-leach/niche-model

max_forks_repo_name: clint-leach/niche-model


text:
#=
This is a modified version of MS_Import by Saer Samanipour: https://bitbucket.org/SSamanipour/ms_import.jl/src/master/

Copyright (c) 2020 Saer Samanipour, PhD, Computational Mass Spec Lab (CMSL)

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
=#

using LightXML
using Codecs

#############################################################################
# import_files_light

function import_files_light(pathin,filenames,mz_thresh=[0,0],Int_thresh=0)

    m=split(filenames[1],".")

    if m[end] == "mzxml" || m[end] == "mzXML" || m[end] == "MZXML"

        if isa(filenames,Array)==1
            path_in=joinpath(pathin,filenames[])
        else
            path_in=joinpath(pathin,filenames)
        end

        chrom = mzxml_read_light(path_in,mz_thresh,Int_thresh)
    end

    return (chrom)

end

#############################################################################
# mzxml_read_light

function mzxml_read_light(path2mzxml,mz_thresh,Int_thresh)
        xdoc = parse_file(path2mzxml);
        xroot = root(xdoc);

        if LightXML.name(xroot) != "mzXML"
            error("Not an mzXML file")
        end
        # Find the msRun node
        msRun = find_element(xroot, "msRun");
        mslevels = ["1"]
        row_length, col_length = determine_matrix_size(msRun,mslevels)
        retentionTime,mz,mz_int=read_scan_light(msRun,mz_thresh,Int_thresh, row_length, col_length, mslevels)

    chrom=Dict("MS1" => Dict{String, Any}(
        "Rt" => retentionTime,
        "Mz_values" => mz,
        "Mz_intensity" => mz_int

    ))
    free(xdoc)
    return chrom

end # function

###############################################################################
# Reading the scans

function read_scan_light(msRun,mz_thresh,Int_thresh, row_length, col_length, mslevels)
    retentionTime = zeros(Float32, row_length)
    mz, mz_int = zeros(Float32, row_length, col_length), zeros(Int64, row_length, col_length)
    i = 0
    for line in child_elements(msRun)
        if LightXML.name(line) == "scan" && attribute(line, "msLevel") in mslevels
            i += 1
            retentionTime1,mz1,I=read_scan_info_light(line,mz_thresh,Int_thresh)
            # if retentionTime1 == false #TEMP
            #     continue
            # end
            retentionTime[i] = round(retentionTime1 / 60, digits=2)
            mz[i, 1:length(mz1)] = mz1
            mz_int[i, 1:length(I)] = I
            # retentionTime=vcat(retentionTime,retentionTime1) #TEMP
            # mz=vcat(mz,[mz1])
            # mz_int=vcat(mz_int,[I])
        end

    end
    return(retentionTime,mz,mz_int)

end # function

###############################################################################
# Extracting info from the scans

function read_scan_info_light(line,mz_thresh,Int_thresh)

    precisiondict = Dict("32" => (Float32, Float32, Float32), "64" => (Float64, Float64, Float64))
    # line = find_element(msRun,"scan")
    # println(line)

    tstr=attribute(line, "retentionTime")
    retentionTime = parse(Float64, tstr[3:end-1])

    # # Only read ms level 1 #TEMP
    # if attribute(line, "msLevel") != "1"
    #     return false, 0, 0
    # end

    peak = find_element(line, "peaks")
    # 1840-element Array{UInt8,1}:
    if attribute(peak, "compressionType") ==  "zlib"
        data = decode(Zlib, decode(Base64, content(peak)))
    else
        data = decode(Base64, content(peak))
    end

    # data = decode(Zlib, decode(Base64, content(peak)))
    # println(bytestring(data))
    TI, T, nochildren = precisiondict[attribute(peak, "precision")]
    A = reinterpret(TI, data)

    bo = attribute(peak, "byteOrder")
    if bo == "network"
        ntoh!(A)
    else
        error("Don't know what to do with byteOrder $bo")
    end
    I = A[2:2:end]

    mz = reinterpret(T, A[1:2:end])

    if Int_thresh > 0
        mz1=mz[findall(x -> Int_thresh<= x, I)]
        I1=I[findall(x -> Int_thresh<= x, I)]
    else
        mz1=mz
        I1=I
    end

    if mz_thresh[2]>0
        mz2=mz1[findall(x -> mz_thresh[1]<= x <= mz_thresh[2], mz1)]
        I2=I1[findall(x -> mz_thresh[1]<= x <= mz_thresh[2], mz1)]
    elseif mz_thresh[1] == 0
        mz2=mz1[mz1 .>= mz_thresh[1]]
        I2=I1[mz1 .>= mz_thresh[1]]

    end

    return(retentionTime,mz2,I2)
end # function

###############################################################################
# Determines scan with highest peak count and number of scans
# This determines spectrum mz and intensity matrix size

function determine_matrix_size(msRun,mslevels)

    row_length, col_length = 0, 0
    for line in child_elements(msRun)
        if LightXML.name(line) == "scan" && attribute(line, "msLevel") in mslevels
            row_length += 1
            peakscount = parse(Int, attribute(line, "peaksCount"))

            # If peakscount is longer than current length, change to peakscount
            col_length = peakscount > col_length ? peakscount : col_length
        end
    end

    return row_length, col_length

end # function

##############################################################################
# XML encoding

function ntoh!(A)
    for i = 1:length(A)
        A[i] = ntoh(A[i])
    end
end

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


text:
using RobotDynamics
using StaticArrays
using Test
using BenchmarkTools
using ForwardDiff
using FiniteDiff
using LinearAlgebra
using Random

using RobotDynamics: @autodiff, state, control, KnotPoint, getdata, getstate, getcontrol, getparams,
                     FiniteDifference, ForwardAD, StaticReturn, InPlace

using RobotDynamics: evaluate, evaluate!, jacobian!
function test_allocs(fun)
    n,m,p = size(fun)
    z_ = @SVector randn(n+m)
    z = KnotPoint{n,m}(z_,1.0,0.1)
    y = zeros(p)
    J = zeros(p,n+n)
    allocs = 0
    allocs += @allocated evaluate(fun, z)
    allocs += @allocated evaluate!(fun, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.ForwardAD(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.FiniteDifference(), fun, J, y, z)

    # check for zero allocations with normal array inputs for in place methods
    z = KnotPoint{n,m}(Vector(z_),1.0,0.1)
    allocs += @allocated evaluate(fun, z)
    allocs += @allocated evaluate!(fun, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.ForwardAD(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.InPlace(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    allocs += @allocated jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.FiniteDifference(), fun, J, y, z)
end

function test_fun(fun)
    Random.seed!(1)
    @test size(fun) == (2,2,3)

    z_ = @SVector randn(4)
    z = KnotPoint{2,2}(z_,1.0,0.1)
    y = zeros(3)
    J = zeros(3,4)
    J0 = zeros(3,4)

    evaluate!(fun, y, z)
    @test y ≈ evaluate(fun, z)

    jacobian!(fun, J0, y, z)

    jacobian!(RobotDynamics.InPlace(), RobotDynamics.UserDefined(), fun, J, y, z)
    @test J ≈ J0
    jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.UserDefined(), fun, J, y, z)
    @test J ≈ J0

    jacobian!(RobotDynamics.InPlace(), RobotDynamics.ForwardAD(), fun, J, y, z)
    @test J ≈ J0
    jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
    @test J ≈ J0

    jacobian!(RobotDynamics.InPlace(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    @test J ≈ J0 atol = 1e-6
    jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.FiniteDifference(), fun, J, y, z)
    @test J ≈ J0 atol = 1e-6

    run_alloc_tests && (@test test_allocs(fun) == 0)
end

##
struct TestFun0{CH} <: RobotDynamics.AbstractFunction
    cfg::ForwardDiff.JacobianConfig{Nothing, Float64, CH, Tuple{Vector{ForwardDiff.Dual{Nothing, Float64, CH}}, Vector{ForwardDiff.Dual{Nothing, Float64, CH}}}}
    cache::FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}
    function TestFun0()
        n,m,p = 2,2,3
        cfg = ForwardDiff.JacobianConfig(nothing, zeros(p), zeros(n+m))
        cache = FiniteDiff.JacobianCache(zeros(n+m), zeros(p))
        new{length(cfg.seeds)}(cfg, cache)
    end
end
function RobotDynamics.evaluate(::TestFun0, x, u, p)
    return SA[cos(x[1]) * u[1], sin(x[2]^2 * x[1]) * u[2], exp(x[2] + x[1]/10)] * p[1]
end
function RobotDynamics.evaluate!(::TestFun0, y, x, u, p)
    y[1] = cos(x[1])  * u[1]
    y[2] = sin(x[2]^2 * x[1]) * u[2]
    y[3] = exp(x[2] + x[1]/10)
    y .*= p[1]
    return nothing
end
RobotDynamics.state_dim(::TestFun0) = 2
RobotDynamics.control_dim(::TestFun0) = 2
RobotDynamics.output_dim(::TestFun0) = 3

function RobotDynamics.jacobian!(::TestFun0, J, y, x, u, p)
    J .= 0
    J[1,1] = -sin(x[1]) * u[1]
    J[1,3] = cos(x[1])
    J[2,1] = x[2]^2 * cos(x[2]^2 * x[1]) * u[2]
    J[2,2] = 2 * x[1] * x[2] * cos(x[2]^2 * x[1]) * u[2]
    J[2,4] = sin(x[2]^2 * x[1])
    J[3,1] = exp(x[2] + x[1]/10) / 10
    J[3,2] = exp(x[2] + x[1]/10)
    J .*= p[1]
    return nothing
end

function RobotDynamics.jacobian!(::StaticReturn, ::ForwardAD, fun::TestFun0, J, y, z)
    f(_z) = RobotDynamics.evaluate(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
    J .= ForwardDiff.jacobian(f, getdata(z))
    return nothing
end

function RobotDynamics.jacobian!(::InPlace, ::ForwardAD, fun::TestFun0, J, y, z)
    f!(_y,_z) = RobotDynamics.evaluate!(fun, _y, getstate(z, _z), getcontrol(z, _z), getparams(z))
    ForwardDiff.jacobian!(J, f!, y, getdata(z), fun.cfg)
    return nothing
end

function RobotDynamics.jacobian!(::StaticReturn, ::FiniteDifference, fun::TestFun0, J, y, z)
    f!(_y,_z) = _y .= RobotDynamics.evaluate(fun, getstate(z, _z), getcontrol(z, _z), getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, getdata(z), fun.cache)
    return nothing
end

function RobotDynamics.jacobian!(::InPlace, ::FiniteDifference, fun::TestFun0, J, y, z)
    f!(_y,_z) = RobotDynamics.evaluate!(fun, _y, getstate(z, _z), getcontrol(z, _z), getparams(z))
    FiniteDiff.finite_difference_jacobian!(J, f!, getdata(z), fun.cache)
    return nothing
end

fun = TestFun0()
n,m,p = size(fun)
x = @SVector randn(n)
u = @SVector randn(m)
t = 1.2
dt = 0.1

zs = KnotPoint{n,m}([x;u],t,dt)
z = KnotPoint{n,m}(Vector([x;u]),t,dt)
z_ = copy(z.z)
@test getstate(zs, z_) isa SVector{n}
@test getcontrol(zs, z_) isa SVector{m}
@test getstate(z, z_) isa SubArray
@test getcontrol(z, z_) isa SubArray

test_fun(fun)

##############################
# Autogen With Inner Constructor
##############################
@autodiff struct TestFun <: RobotDynamics.AbstractFunction
    a::Int
    function TestFun()
        new(1.0)
    end
end

RobotDynamics.state_dim(::TestFun) = 2
RobotDynamics.control_dim(::TestFun) = 2
RobotDynamics.output_dim(::TestFun) = 3

function RobotDynamics.evaluate(::TestFun, x, u)
    return SA[sin(x[1]), cos(x[2]), u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFun, y, x, u)
    y[1] = sin(x[1])
    y[2] = cos(x[2])
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFun, J, y, x, u)
    J .= 0
    J[1,1] = cos(x[1])
    J[2,2] = -sin(x[2])
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end

fun = TestFun()
@test :cfg ∈ fieldnames(typeof(fun))
@test :cache ∈ fieldnames(typeof(fun))
@test fun.a == 1.0
test_fun(fun)

##############################
# With Input Parameters
##############################
@autodiff struct TestFunTime <: RobotDynamics.AbstractFunction end

function RobotDynamics.evaluate(::TestFunTime, x, u, t)
    return SA[sin(x[1]) * t.t, cos(x[2]) * t.t, u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFunTime, y, x, u, t)
    y[1] = sin(x[1]) * t.t
    y[2] = cos(x[2]) * t.t
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFunTime, J, y, x, u, t)
    J .= 0
    J[1,1] = cos(x[1]) * t.t
    J[2,2] = -sin(x[2]) * t.t
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end
RobotDynamics.state_dim(::TestFunTime) = 2
RobotDynamics.control_dim(::TestFunTime) = 2
RobotDynamics.output_dim(::TestFunTime) = 3

fun = TestFunTime()
z_ = @SVector zeros(4)
z = KnotPoint{2,2}(z_, 1.2, 0.1)
y = zeros(3)
J = zeros(3,4)
out = RobotDynamics.evaluate(fun, z)
RobotDynamics.evaluate!(fun, y, z)
@test y ≈ out
@test out[2] == z.t

jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
@test J[1,1] == z.t
z.t = 2.4
jacobian!(RobotDynamics.StaticReturn(), RobotDynamics.ForwardAD(), fun, J, y, z)
@test J[1,1] == z.t

test_fun(fun)

##############################
# With Type Parameters
##############################
@autodiff struct TestFunParam{T} <: RobotDynamics.AbstractFunction
    a::T
    function TestFunParam(a::T) where T
        new{T}(a)
    end
end

RobotDynamics.state_dim(::TestFunParam) = 2
RobotDynamics.control_dim(::TestFunParam) = 2
RobotDynamics.output_dim(::TestFunParam) = 3

function RobotDynamics.evaluate(::TestFunParam, x, u)
    return SA[sin(x[1]), cos(x[2]), u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFunParam, y, x, u)
    y[1] = sin(x[1])
    y[2] = cos(x[2])
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFunParam, J, y, x, u)
    J .= 0
    J[1,1] = cos(x[1])
    J[2,2] = -sin(x[2])
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end

fun = TestFunParam(1.0)
@test fun.a == 1.0
test_fun(fun)

##############################
# Without inner constructor and inherited
##############################
abstract type TestFunBase <: RobotDynamics.AbstractFunction end

@autodiff struct TestFunInherited <: TestFunBase end
RobotDynamics.state_dim(::TestFunBase) = 2
RobotDynamics.control_dim(::TestFunBase) = 2
RobotDynamics.output_dim(::TestFunBase) = 3

function RobotDynamics.evaluate(::TestFunBase, x, u)
    return SA[sin(x[1]), cos(x[2]), u[1] * exp(u[2])]
end

function RobotDynamics.evaluate!(::TestFunBase, y, x, u)
    y[1] = sin(x[1])
    y[2] = cos(x[2])
    y[3] = u[1] * exp(u[2])
end

function RobotDynamics.jacobian!(::TestFunBase, J, y, x, u)
    J .= 0
    J[1,1] = cos(x[1])
    J[2,2] = -sin(x[2])
    J[3,3] = exp(u[2])
    J[3,4] = u[1] * exp(u[2])
    return nothing
end

fun = TestFunInherited()
test_fun(fun)


text:
abstract type ArchimedeanCopula <: AbstractCopula end

function rosenblatt(M::AbstractMatrix, c::ArchimedeanCopula)
    u1 = M[:, 1]
    u2 = M[:, 2]

    return [u1 (φ².(φ⁻¹.(u1, c) + φ⁻¹.(u2, c), c) ./ φ².(φ⁻¹.(u1, c), c))]
end

function cdf(c::ArchimedeanCopula, u1::Real, u2::Real)
    return φ(φ⁻¹(u1, c) + φ⁻¹(u2, c), c)
end

function density(c::ArchimedeanCopula, u1::Real, u2::Real)
    return φ²(φ⁻¹(u1, c) + φ⁻¹(u2, c), c) * ∇φ⁻¹(u1, c) * ∇φ⁻¹(u2, c)
end

function sample(c::ArchimedeanCopula, n::Int64)
    @assert n > 0
    return inverse_rosenblatt(rand(n, 2), c)
end

function (c::ArchimedeanCopula)(
    X::ContinuousUnivariateDistribution, Y::ContinuousUnivariateDistribution
)
    # Construct a joint distribution using the copula
    return Joint(X, Y, c)
end

include("independence.jl")
include("clayton.jl")


text:
using DataStructures
using Base.Meta
using AutoHashEquals

import Base.ImmutableDict

# @auto_hash_equals struct EClass
#     id::Int64
# end

@auto_hash_equals struct ENode{X}
    head::Any
    args::Vector{Int64}
    metadata::Union{Nothing, NamedTuple}
end

ariety(n::ENode) = length(n.args)

function ENode(e, c_ids::AbstractVector{Int64})
    # @assert length(getargs(e)) == length(c_ids)
    # static_args = MVector{length(c_ids), Int64}(c_ids...)
    ENode{typeof(e)}(gethead(e), c_ids, getmetadata(e))
end

ENode(e) = ENode(e, Int64[])

ENode(a::ENode) =
    error("constructor of ENode called on enode. This should never happen")

# string representation of the rule
function Base.show(io::IO, x::ENode)
    print(io, "(", x.head)
    n = ariety(x)
    if n == 0
        print(io, ")")
        return
    else
        print(io, " ")
    end
    for i ∈ 1:n
        if i < n
            print(io, x.args[i], " ")
        else
            print(io, x.args[i], ")")
        end
    end
end


text:
type Pair
    class1 :: Class
    class2 :: Class
    pclass :: Class

    pair :: DFunction
    first :: DFunction
    second :: DFunction
    first! :: DFunction
    second! :: DFunction

    axioms :: Array{Axiom, 1}
end

function Base.show(io :: IO, m :: Pair)
    print(io, string("Module.Pair(", m.class1.class_name, ",", m.class2.class_name, "->", m.pclass.class_name, ")"))
end

function Pair(class1 :: Class, class2 :: Class, pclass :: Class, symbols = [:xy, :x, :y])
    pair = DFunction("pair", [class1, class2], pclass)
    first = DFunction("first", [pclass], class1)
    second = DFunction("second", [pclass], class2)
    first! = DFunction("first!", [pclass, class1], pclass)
    second! = DFunction("second!", [pclass, class2], pclass)

    tsym = symbols[1]
    fsym = symbols[2]
    ssym = symbols[3]
    axiom_first = Axiom([first, [pair, fsym, ssym]], fsym)
    axiom_second = Axiom([second, [pair, fsym, ssym]], ssym)
    axiom_first! = Axiom([first!, tsym, fsym], [pair, fsym, [second, tsym]])
    axiom_second! = Axiom([second!, tsym, ssym], [pair, [first, tsym], ssym])
    axioms = [axiom_first, axiom_second, axiom_first!, axiom_second!]

    Pair(class1, class2, pclass, pair, first, second, first!, second!, axioms)
end

function Pair(class1 :: Class, class2 :: Class, type_len :: Int64, symbols = [:xy, :x, :y])
    post_name = string("_", class1.class_name, "_", class2.class_name)
    pair_name = string("Pair", post_name)
    pclass = Class(pair_name, type_len)
    Pair(class1, class2, pclass, symbols)
end

type List
    class :: Class
    lclass :: Class

    empty :: Object

    cons :: DFunction
    first :: DFunction
    rest :: DFunction
    first! :: DFunction
    rest! :: DFunction

    axioms :: Array{Axiom, 1}
end

function Base.show(io :: IO, m :: List)
    print(io, string("Module.List(", m.class.class_name, "->", m.lclass.class_name, ")"))
end

function List(class :: Class, lclass :: Class, symbols = [:xxs, :x, :xs])
    empty = Object("empty", lclass)
    empty.value .= -1.

    cons = DFunction("cons", [class, lclass], lclass)
    first = DFunction("first", [lclass], class)
    rest = DFunction("rest", [lclass], lclass)
    first! = DFunction("first!", [lclass, class], lclass)
    rest! = DFunction("rest!", [lclass, lclass], lclass)

    tsym = symbols[1]
    fsym = symbols[2]
    rsym = symbols[3]
    axiom_first = Axiom([first, [cons, fsym, rsym]], fsym)
    axiom_rest = Axiom([rest, [cons, fsym, rsym]], rsym)
    axiom_first! = Axiom([first!, tsym, fsym], [cons, fsym, [rest, tsym]])
    axiom_rest! = Axiom([rest!, tsym, rsym], [cons, [first, tsym], rsym])
    axioms = [axiom_first, axiom_rest, axiom_first!, axiom_rest!]

    List(class, lclass, empty, cons, first, rest, first!, rest!, axioms)
end

function List(class :: Class, type_len :: Int64, symbols = [:xxs, :x, :xs])
    list_name = string("List", "_", class.class_name)
    lclass = Class(list_name, type_len)
    List(class, lclass, symbols)
end

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
export read_SDPPrimal, set_constraints!, set_vartypes!, set_blocks!, set_linvars!, set_matrices!, set_linear!, set_const!, print

function read_SDPPrimal(path::String)
  BLOCKS = readdlm(joinpath(path, "matrix.sdp"), String)
  if isfile(joinpath(path, "lin.sdp")) && length(matchall(r"\n", readstring(joinpath(path, "lin.sdp")))) > 7
    LINEAR = readdlm(joinpath(path, "lin.sdp"), String)
  else
    LINEAR = []
  end
  CONST = readdlm(joinpath(path, "const.sdp"), String)
  VAR_TYPES = readdlm(joinpath(path, "types.sdp"), String)

  SDP_Instance(VAR_TYPES,
               BLOCKS,
               LINEAR,
               CONST)
end

"""
set_constraints!(sdp::SDP_Problem, instance::SDP_Instance)

Build `name_to_ctr` with explicit constraint parameters from instance with ==0 as default.
"""
function set_constraints!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  # Collect constraints names
  # ctr_names = SortedSet{Tuple{String, String, String}}([(instance.BLOCKS[i, 1], instance.BLOCKS[i, 2], instance.BLOCKS[i, 3]) for i=1:size(instance.BLOCKS, 1)])
  # union!(ctr_names, [(instance.LINEAR[i, 1], instance.LINEAR[i, 2], instance.LINEAR[i, 3]) for i=1:size(instance.LINEAR, 1)])
  # union!(ctr_names, [(instance.CONST[i, 1], instance.CONST[i, 2], instance.CONST[i, 3]) for i=1:size(instance.CONST, 1)])

ctr_names = SortedSet{SDP_Moment}([(instance.CONST[i, 1], instance.CONST[i, 2], instance.CONST[i, 3]) for i=1:size(instance.CONST, 1)])

  obj_keys = SortedSet{SDP_Moment}()
  for ctr_name in ctr_names
    if (ctr_name[1:2] == ("1", "1"))
      push!(obj_keys, ctr_name)
      delete!(ctr_names, ctr_name)
    end
  end

  length(obj_keys) == 0 && error(LOGGER, "No ctrkey matching objective key (\"1\", \"1\", .)")
  sdp.obj_keys = obj_keys

  ctr_id = 1
  for ctr_name in ctr_names
    sdp.name_to_ctr[ctr_name] = (ctr_id, "EQ", 0, 0)
    sdp.id_to_ctr[ctr_id] = ctr_name
    ctr_id += 1
  end

  if debug
  end
end

"""
set_vartypes!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)

Input all matrix varibales structural information regarding name and kind in the appropriate attributes of `SDP_Problem`.
"""
function set_vartypes!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  n_sdp = 0
  for i=1:size(instance.VAR_TYPES, 1)
    (block_name, block_type) = instance.VAR_TYPES[i, :]

    if block_type == "SDP"
      n_sdp += 1
      block = SDP_Block(n_sdp, block_name)
      sdp.name_to_sdpblock[block_name] = block
      sdp.id_to_sdpblock[n_sdp] = block

    else
      error(LOGGER, "set_vartypes()!: Unknown blockvar type $block_type for variable $block_name")
    end
  end

  if debug
  end
end

"""
set_blocks!(sdp::SDP_Problem, instance::SDP_Instance)

Fill all variables blocks with their elementary variables previously declared by `set_vartypes!`.
"""
function set_blocks!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  for i in 1:size(instance.BLOCKS, 1)
    block_name, var1, var2 = instance.BLOCKS[i, 4:6]

    if haskey(sdp.name_to_sdpblock, block_name)
      cur_blockvar = sdp.name_to_sdpblock[block_name]

      # Adding vars and ids to SDP block
      if !haskey(cur_blockvar.var_to_id, var1)
        cur_blockvar.var_to_id[var1] = length(cur_blockvar.var_to_id) + 1
      end
      if !haskey(cur_blockvar.var_to_id, var2)
        cur_blockvar.var_to_id[var2] = length(cur_blockvar.var_to_id) + 1
      end
    else
      error(LOGGER, "set_blocks!(): Unknown block_kind $(sdp.block_to_kind) for i=$i")
    end
  end

  if debug
  end
end

"""
set_linvars!(sdp::SDP_Problem, sdp_instance::SDP_Instance)

Set the scalar variable name to id dict of the `sdp` structure.
"""
function set_linvars!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  for i in 1:size(instance.LINEAR, 1)
    var = instance.LINEAR[i, 4]
    if !haskey(sdp.scalvar_to_id, var)
      sdp.scalvar_to_id[var] = length(sdp.scalvar_to_id) + 1
    end
  end
end

function set_matrices!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  for i=1:size(instance.BLOCKS, 1)
    ctr_name = (instance.BLOCKS[i, 1], instance.BLOCKS[i, 2], instance.BLOCKS[i, 3])
    (block_name, var1, var2, coeff) = instance.BLOCKS[i, 4:7]

    # Sort variables for triangular matrix storage
    var1, var2 = min(var1, var2), max(var1, var2)

    if haskey(sdp.name_to_sdpblock, block_name)
      if !haskey(sdp.matrices, (ctr_name, block_name, var1, var2))
        sdp.matrices[(ctr_name, block_name, var1, var2)] = parse(coeff)
      else
        error(LOGGER, "set_matrices!(): sdp.matrices already has key ($ctr_name, $block_name, $var1, $var2) with val $(sdp.matrices[(ctr_name, block_name, var1, var2)]), $(parse(coeff))")
      end

    else
      error(LOGGER, "set_matrices!(): Unhandled matrix var $block_name")
    end
  end

  if debug
  end
end

function set_linear!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  if length(instance.LINEAR) != 0
    for i=1:size(instance.LINEAR, 1)
      ctr_name = (instance.LINEAR[i, 1], instance.LINEAR[i, 2], instance.LINEAR[i, 3])
      (var, coeff) = instance.LINEAR[i, 4:5]

      if !haskey(sdp.linear, (ctr_name, var))
        sdp.linear[(ctr_name, var)] = parse(coeff)
      else
        error(LOGGER, "set_linear!(): sdp.linear already has key ($ctr_name, $var) with val $(sdp.linear[(ctr_name, var)]). New val is $(parse(coeff))")
      end
    end
  end

  if debug
  end
end

function set_const!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  for i=1:size(instance.CONST, 1)
    ctr_name = (instance.CONST[i, 1], instance.CONST[i, 2], instance.CONST[i, 3])
    coeff = instance.CONST[i, 4]

    @assert !haskey(sdp.cst_ctr, ctr_name)
    if coeff != 0
      sdp.cst_ctr[ctr_name] = parse(coeff)
    end
  end

  if debug
  end
end

function print(io::IO, sdp::SDP_Problem)
  for cstr in sort(keys(sdp.name_to_sdpblock))
    block = sdp.name_to_sdpblock[cstr]
    println(io, "  sdp   : $cstr -> $block")
  end

  println(io, "  objk  : $(sdp.obj_keys)")

  for name in sort(keys(sdp.name_to_ctr))
    ctr = sdp.name_to_ctr[name]
    println(io, "  ctr   : $name \t $ctr")
  end

  for name in sort(keys(sdp.matrices))
    mat = sdp.matrices[name]
    println(io, "  matrix: $name \t $mat")
  end

  for name in sort(keys(sdp.linear))
    lin = sdp.linear[name]
    println(io, "  lin   : $name \t $lin")
  end

  for name in sort(keys(sdp.cst_ctr))
    cst = sdp.cst_ctr[name]
    println(io, "  cst   : $name \t $cst")
  end
end