text:
CRCBS.build_env(solver,pc_mapf::AbstractPC_MAPF,env::SearchEnv,node::ConstraintTreeNode,i::Int) = build_env(solver,pc_mapf,env,node,AgentID(i))
function construct_heuristic_model(solver::PIBTPlanner,env_graph;
        ph = PerfectHeuristic(get_dist_matrix(env_graph)),
        kwargs...)
    construct_composite_heuristic(ph,ph,NullHeuristic())
end
function construct_cost_model(solver::PIBTPlanner,
        schedule, cache, problem_spec, env_graph, primary_objective=SumOfMakeSpans(); kwargs...)
    cost_model = construct_composite_cost_model(
        typeof(primary_objective)(schedule,cache),
        FullCostModel(maximum,TravelTime()),
        FullCostModel(maximum,TravelDistance())
        )
    heuristic_model = construct_heuristic_model(solver,env_graph;kwargs...)
    cost_model, heuristic_model
end

function CRCBS.is_consistent(cache::PIBTCache,pc_mapf::PC_MAPF)
    search_env = CRCBS.get_solution(cache)
    n_jobless = 0
    # for e in CRCBS.get_envs(cache)
    #     if !CRCBS.is_valid(get_goal(e))
    #         n_jobless += 1
    #     end
    # end
    if length(get_cache(search_env).closed_set) + n_jobless >= nv(get_schedule(search_env))
        return true
    elseif !any(CRCBS.is_valid, map(e->CRCBS.get_goal(e),CRCBS.get_envs(cache)))
        # all robots are done with their jobs
        return true
    end
    return false
end

function CRCBS.pibt_priority_law(solver,pc_mapf::PC_MAPF,cache,i)
    env = CRCBS.get_envs(cache)[i]
    return (
        ~CRCBS.is_valid(env,get_goal(env)),
        ~isa(get_node(env),DEPOSIT),
        ~isa(get_node(env),CARRY),
        ~isa(get_node(env),COLLECT),
        # minimum(cache.solution.cache.slack[get_vtx(
        #     cache.solution.schedule,env.node_id)]),
        minimum(get_slack(cache.solution,env.node_id)),
        -CRCBS.get_timers(cache)[i],
        i
    )
end

function CRCBS.pibt_preprocess!(solver,pc_mapf::PC_MAPF,cache)
    update_planning_cache!(solver,CRCBS.get_solution(cache))
    CRCBS.pibt_update_envs!(solver,pc_mapf,cache)
end

"""
    CRCBS.pibt_update_solution!(solver,pc_mapf::PC_MAPF,solution::SearchEnv,cache)

Overridden to update the SearchEnv
"""
function CRCBS.pibt_update_solution!(solver,solution::SearchEnv,cache)
    # update route plan
    CRCBS.pibt_update_solution!(solver,get_route_plan(solution),cache)
    # NOTE Update planning cache so that all goal times can be updated
    update_planning_cache!(solver,solution)
    solution
end
function pibt_info_string(cache,i,
        env=CRCBS.get_envs(cache)[i],
        s=CRCBS.get_states(cache)[i],
        a=CRCBS.get_actions(cache)[i]
    )
    string("\t",i,": ", get_vtx(cache.solution.schedule,env.node_id)," ",
        string(get_node(env))," s=",string(s),", a=",string(a),
        ", goal=",string(env.goal),"\n")
end
function pibt_info_strings(cache)
    info_strings = String[]
    for i in CRCBS.get_ordering(cache)
        str = pibt_info_string(cache,i)
        push!(info_strings,str)
    end
    return info_strings
end
function CRCBS.pibt_update_envs!(solver,pc_mapf::PC_MAPF,cache)
    @log_info(2,verbosity(solver),"PIBT iteration $(iterations(solver)) PRE UPDATE:\n",
        pibt_info_strings(cache)...)
    solution = CRCBS.get_solution(cache)
    node = initialize_root_node(solver,pc_mapf)
    done = true
    sweeps = 0
    while true
        done = true
        for (i,p) in enumerate(get_paths(solution))
            v_next = get_next_vtx_matching_agent_id(solution,RobotID(i))
            if has_vertex(solution.schedule,v_next)
                # rebuild env to ensure that the goal time is up to date
                CRCBS.get_envs(cache)[i] = build_env(solver,pc_mapf,solution,node,AgentID(i))
            end
            env = CRCBS.get_envs(cache)[i]
            v = get_vtx(solution.schedule,env.node_id)
            sp = get_final_state(p)
            if (v in solution.cache.active_set) && is_goal(env,sp) && CRCBS.is_valid(env,get_goal(env))
                done = false
                update_env!(solver,solution,v,p)
                update_planning_cache!(solver,solution)
                v_next = get_next_vtx_matching_agent_id(solution,RobotID(i))
                @assert v_next != v
                if has_vertex(solution.schedule,v_next)
                    CRCBS.get_timers(cache)[i] = 0
                end
            end
        end
        sweeps += 1
        done ? break : nothing
    end
    @log_info(2,verbosity(solver),"PIBT iteration $(iterations(solver)) POST UPDATE:\n",
        pibt_info_strings(cache)...)
    return cache
end


text:
# Support for ZeroMQ, a network and interprocess communication library

module ZMQ
using ZeroMQ_jll

using Base.Libc: EAGAIN
using FileWatching: UV_READABLE, uv_pollcb, _FDWatcher
import Sockets
using Sockets: connect, bind, send, recv
import Base.GC: @preserve

export
    #Types
    StateError,Context,Socket,Message,
    #functions
    set, subscribe, unsubscribe,
    #Constants
    IO_THREADS,MAX_SOCKETS,PAIR,PUB,SUB,REQ,REP,ROUTER,DEALER,PULL,PUSH,XPUB,XSUB,XREQ,XREP,UPSTREAM,DOWNSTREAM,MORE,POLLIN,POLLOUT,POLLERR,STREAMER,FORWARDER,QUEUE,SNDMORE,
    #Sockets
    connect, bind, send, recv

include("constants.jl")
include("optutil.jl")
include("error.jl")
include("context.jl")
include("socket.jl")
include("sockopts.jl")
include("message.jl")
include("comm.jl")

function __init__()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    patch = Ref{Cint}()
    ccall((:zmq_version, libzmq), Cvoid, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), major, minor, patch)
    global version = VersionNumber(major[], minor[], patch[])
    if version < v"3"
        error("ZMQ version $version < 3 is not supported")
    end
    atexit() do
        close(_context)
    end
end

end

ext: jl


text:
# This file was generated, do not modify it. # hide
no_miss = .!ismissing.(y)
y = collect(skipmissing(y))
X = X[no_miss, :]
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=424);


text:
struct WithGenericT{T}
    a::T
end

struct WithNestedGenericT{T}
    a::WithGenericT{T}
end

struct WithSetGeneric
    a::WithGenericT{Int64}
end

struct WithValueType{N}
    a::Int64
end

struct WithGenericUnionAll
    a::WithGenericT
end

struct WithGenericTuple{T}
    a::Tuple{T}
end

struct WithSetGenericTuple
    a::Tuple{WithGenericT{Int64}}
end

struct WithPropagatedLifetime
    a::WithGenericT{Module}
end

struct WithPropagatedLifetimes
    a::WithGenericT{Tuple{Int32, WithGenericT{Array{Int32, 2}}}}
end

@testset "Structs with generic fields" begin
    @test begin
        b = JlrsReflect.reflect([WithGenericT])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(WithGenericT)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithGenericT")]
        pub struct WithGenericT<T>
        where
            T: ::jlrs::layout::valid_layout::ValidLayout + Clone,
        {
            pub a: T,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithNestedGenericT])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(WithNestedGenericT)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithNestedGenericT")]
        pub struct WithNestedGenericT<T>
        where
            T: ::jlrs::layout::valid_layout::ValidLayout + Clone,
        {
            pub a: WithGenericT<T>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithSetGeneric])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithSetGeneric] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.WithSetGeneric")]
        pub struct WithSetGeneric {
            pub a: WithGenericT<i64>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithValueType])
        sb = JlrsReflect.StringWrappers(b)

        sb[JlrsReflect.basetype(WithValueType)] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithValueType")]
        pub struct WithValueType {
            pub a: i64,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithGenericUnionAll])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithGenericUnionAll] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithGenericUnionAll")]
        pub struct WithGenericUnionAll<'frame, 'data> {
            pub a: ::jlrs::wrappers::ptr::ValueRef<'frame, 'data>,
        }"""
    end

    @test_throws ErrorException begin
        JlrsReflect.reflect([WithGenericTuple])
    end

    @test begin
        b = JlrsReflect.reflect([WithSetGenericTuple])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithSetGenericTuple] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.WithSetGenericTuple")]
        pub struct WithSetGenericTuple {
            pub a: ::jlrs::wrappers::inline::tuple::Tuple1<WithGenericT<i64>>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithPropagatedLifetime])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithPropagatedLifetime] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithPropagatedLifetime")]
        pub struct WithPropagatedLifetime<'frame> {
            pub a: WithGenericT<::jlrs::wrappers::ptr::ModuleRef<'frame>>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithPropagatedLifetimes])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithPropagatedLifetimes] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithPropagatedLifetimes")]
        pub struct WithPropagatedLifetimes<'frame, 'data> {
            pub a: WithGenericT<::jlrs::wrappers::inline::tuple::Tuple2<i32, WithGenericT<::jlrs::wrappers::ptr::ArrayRef<'frame, 'data>>>>,
        }"""
    end
end


text:
using Plots
pgfplotsx()
plotly()
# gr()
col = [RGB(i) for i in [0x00ab9f,0x7409b7,0xd2cf00,0xff6f00] ] # mah precious colors
using DelimitedFiles

path = "/home/freshstart/DiplomaThesisData/"

median = (data) -> (s=0;for i in data;s+=i;end;return s/length(data))

function dataloader(data)
    a = data[2:end,1]
    b = data[2:end,2]
    c = data[2:end,3]
    d = data[2:end,4]
    time = data[2:end,5:end]
    return a,b,c,d, time
end

function Σum(data)
    sum = 0
    for j in data
        sum += j
    end
    return sum
end

function PLOTS_COMB(data,bs; save = true)

    acc_up,  acc_down, vld_up, vld_down, time = dataloader(data)

    plot(   acc_up,  label="training loss top face",        c=col[1])
    plot!(  vld_up,  label="verification loss top face",    c=col[4])
    plot!(  acc_down,label="training loss bottom face",     c=col[2])
    plot!(  vld_down,label="verification loss bottom face", c=col[3])
    plot!( yscale=:log10,
        minorgrid =true,
        minortick=true,
        xlabel= "Epochs",
        ylabel= "Mean Square Error", yguidefontrotation = -90,
        title="Training and verification loss\n Batch Size = $bs",
        legend= :topright,size=(1000,9/16*1000))
    ylims!((0.9*minimum((minimum(acc_up),minimum(acc_down),minimum(vld_up),minimum(vld_down))),1e-2))
    save ? (display(plot!());savefig(path*"images/Convergence_NN_bs_$bs.pdf")) : display(plot!())

    labels=["training top side","validation top side","training bottom side","validation bottom side"]
    plot(title="Neural Networks\' timelog\n Batch Size = $bs",xlabel = "Epochs",ylabel="Time [ms]"
    , yguidefontrotation = -90,minorticks = true,)
    sum = 0
    for i in 1:4
        plot!(time[2:end,i]*1000,label = labels[i], c =col[i])
        for j in time[:,i]
            sum+=j
        end
    end

    sum /= 60
    sum = round(sum,digits = 1)
    plot!(time[2:3],c = RGB(0xFFFFFF), label = "Total training and")
    plot!(time[2:3],c = RGB(0xFFFFFF), label = "validation time = $sum mins")
    save ? (display(plot!());savefig(path*"images/Time_NN_bs_$bs.pdf")) : display(plot!())
end

function c_avg(data;sp=1 ::Int ,ep = 1 ::Int)
    (sp > ep)&(ep != 1) ? error("Start point cannot be greater than end point.") :
    ep == 1 ? a = data[sp:end] : a = data[sp:ep]

    sum = 0
    for i in a
        sum += i
    end
    return sum/length(a)
end

function AveragePerStep(data,bs,step)
    acc_up, acc_down, vld_up, vld_down, time = dataloader(data)

    avg_acc_up   = [c_avg(acc_up  ,sp = i-step+1, ep = i )  for i in range(step,length(acc_up  ),step = step)]
    avg_vld_up   = [c_avg(vld_up  ,sp = i-step+1, ep = i )  for i in range(step,length(vld_up  ),step = step)]
    avg_acc_down = [c_avg(acc_down,sp = i-step+1, ep = i )  for i in range(step,length(acc_down),step = step)]
    avg_vld_down = [c_avg(vld_down,sp = i-step+1, ep = i )  for i in range(step,length(vld_down),step = step)]
    x = [i*step for i in 1:length(avg_acc_up)]

    plot( x, avg_acc_up,  label="averaged training loss top face"       ,c=col[1])
    plot!(x, avg_vld_up,  label="averaged verification loss top face"   ,c=col[4])
    plot!(x, avg_acc_down,label="averaged training loss bottom face"    ,c=col[2])
    plot!(x, avg_vld_down,label="averaged verification loss bottom face",c=col[3])
    plot!(  yscale=:log10,
            minorgrid =true,
            minortick=true,
            xlabel= "Epochs",
            ylabel= "Mean Square Error", yguidefontrotation = -90,
            title="Training and verification loss averaged per $step Epochs \n Batch Size = $bs ",
            legend= :topright,size=(1000,9/16*1000))
    ylims!((0.9*minimum((minimum(acc_up),minimum(acc_down),minimum(vld_up),minimum(vld_down))),1e-2))
    display(plot!())
return avg_acc_up, avg_vld_up, avg_acc_down, avg_vld_down
end

data_025 = readdlm(path*"NeuralNetSaves/MSEsTime.dat")
data_050 = readdlm(path*"NeuralNetSaves_050/MSEsTime.dat")
data_100 = readdlm(path*"NeuralNetSaves_100/MSEsTime.dat")
data_150 = readdlm(path*"NeuralNetSaves_150/MSEsTime.dat")
data_200 = readdlm(path*"NeuralNetSaves_200/MSEsTime.dat")

PLOTS_COMB(data_025,025,save=false)
# PLOTS_COMB(data_050[1:end,:],50,save=false)
# PLOTS_COMB(data_100,100)
# PLOTS_COMB(data_150,150)
# PLOTS_COMB(data_200,200)

step = 40

AveragePerStep(data_025,025,step)
a = AveragePerStep(data_050,050,step)
b = AveragePerStep(data_100,100,step)
c = AveragePerStep(data_150,150,step)
d = AveragePerStep(data_200,200,step)

x = [i*step for i in 1:length(a[1])]
labels=["training top side","validation top side","training bottom side","validation bottom side"]

for i in [1,3]
plot(title = "Comparison between the different Batch Sizes",yaxis = :log10)
plot!(x,a[i], c = col[1], label = labels[i]*" Batch Size = 50")
plot!(x,b[i], c = col[4], label = labels[i]*" Batch Size = 100")
plot!(x,c[i], c = col[3], label = labels[i]*" Batch Size = 150")
plot!(x,d[i], c = col[2], label = labels[i]*" Batch Size = 200")
display(plot!())
# savefig(path*"images/Conv_Comp_BS"*labels[i]*".pdf")
end
plot(title = "Comparison between the different Batch Sizes",)
plot!([Σum(data_050[i,5:end]) for i in 3:length(data_050[2:end,5])], c = col[1], label = "Total Epoch time Batch Size = 50")
plot!([Σum(data_100[i,5:end]) for i in 3:length(data_050[2:end,5])], c = col[4], label = "Total Epoch time Batch Size = 100")
plot!([Σum(data_150[i,5:end]) for i in 3:length(data_050[2:end,5])], c = col[3], label = "Total Epoch time Batch Size = 150")
plot!([Σum(data_200[i,5:end]) for i in 3:length(data_050[2:end,5])], c = col[2], label = "Total Epoch time Batch Size = 200")
display(plot!())
savefig(path*"images/Time_Comp_BS.pdf")

println("Epoch @ Batch Size 50  = ", c_avg([Σum(data_025[i,5:end]) for i in 3:length(data_050[2:end,5])]),"\n",
        "Epoch @ Batch Size 100 = ", c_avg([Σum(data_100[i,5:end]) for i in 3:length(data_050[2:end,5])]),"\n",
        "Epoch @ Batch Size 150 = ", c_avg([Σum(data_150[i,5:end]) for i in 3:length(data_050[2:end,5])]),"\n",
        "Epoch @ Batch Size 200 = ", c_avg([Σum(data_200[i,5:end]) for i in 3:length(data_050[2:end,5])]),"\n",
)

println("Total time @ Batch Size 50  = ", Σum([Σum(data_050[i,5:end]) for i in 3:length(data_050[2:end,5])])/60,"\n",
        "Total time @ Batch Size 100 = ", Σum([Σum(data_100[i,5:end]) for i in 3:length(data_050[2:end,5])])/60,"\n",
        "Total time @ Batch Size 150 = ", Σum([Σum(data_150[i,5:end]) for i in 3:length(data_050[2:end,5])])/60,"\n",
        "Total time @ Batch Size 200 = ", Σum([Σum(data_200[i,5:end]) for i in 3:length(data_050[2:end,5])])/60,"\n",
)


text:
using Test
using LinearAlgebra
import Adapt

@testset "constructors" begin
  xs = VEArray{Int}(undef, 2, 3)
  @test collect(VEArray([1 2; 3 4])) == [1 2; 3 4]
  @test testf(vec, rand(5,3))
  @test Base.elsize(xs) == sizeof(Int)
  @test VEArray{Int, 2}(xs) === xs

  @test_throws ArgumentError Base.unsafe_convert(Ptr{Int}, xs)
  @test_throws ArgumentError Base.unsafe_convert(Ptr{Float32}, xs)

  @test collect(VectorEngine.zeros(2, 2)) == zeros(2, 2)
  @test collect(VectorEngine.ones(2, 2)) == ones(2, 2)

  @test collect(VectorEngine.fill(0, 2, 2)) == zeros(2, 2)
  @test collect(VectorEngine.fill(1, 2, 2)) == ones(2, 2)
end

@testset "adapt" begin
  A = rand(Float32, 3, 3)
  dA = VEArray(A)
  @test Adapt.adapt(Array, dA) == A
  @test Adapt.adapt(VEArray, A) isa VEArray
  @test Array(Adapt.adapt(VEArray, A)) == A
end

@testset "reshape" begin
  A = [1 2 3 4
       5 6 7 8]
  gA = reshape(VEArray(A),1,8)
  _A = reshape(A,1,8)
  _gA = Array(gA)
  @test all(_A .== _gA)
  A = [1,2,3,4]
  gA = reshape(VEArray(A),4)
end

@testset "fill(::SubArray)" begin
  xs = VectorEngine.zeros(3)
  fill!(view(xs, 2:2), 1)
  @test Array(xs) == [0,1,0]
end


text:
function persistence(num)
    count = 0
    while num > 9
        num = reduce(*,digits(num))
        count += 1
    end
    count
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
function is_inside(cur_x::Number, cur_a::Number, cur_b::Number; inclusive::Bool=true)
  if inclusive
    ( cur_a <= cur_x <= cur_b ) && return true
    ( cur_a >= cur_x >= cur_b ) && return true
  else
    ( cur_a < cur_x < cur_b ) && return true
    ( cur_a > cur_x > cur_b ) && return true
  end

  false
end


text:
module SlicingExample
    using PicoQuant

    n = 4
    circ = create_simple_preparation_circuit(n, 2).compose(create_qft_circuit(n))

    # contract the regular way
    tn = convert_qiskit_circ_to_network(circ, decompose=true, transpile=true)
    add_input!(tn, "0"^n)
    wf = full_wavefunction_contraction!(tn, "vector")

    # now decompose to 4 partitions and contract each partition separately
    N = 4
    wfs = []
    for p in 1:N
        tn = convert_qiskit_circ_to_network(circ, decompose=true, transpile=true)
        add_input!(tn, "0"^n)
        bond_labels, bond_values = partition_network_on_virtual_bonds(tn, N, p)
        slice_tensor_network(tn, bond_labels, bond_values)
        push!(wfs, full_wavefunction_contraction!(tn, "vector"))
    end
    if sum(wfs) ≈ wf
        println("Results match with and without partitioning")
        print("Overlap: ")
        println(sum(wf .* conj(sum(wfs))))
    else
        println("Discrepency between results with and without partitioning")
        println("$wf")
        println("$(sum(wfs))")
    end

end


text:
using KitBase, Plots
using KitBase.JLD2
using Base.Threads: @threads
using KitBase.ProgressMeter: @showprogress
cd(@__DIR__)

begin
    set = Setup(case = "layer", space = "1d0f0v", maxTime = 0.2, boundary = ["fix", "fix"], cfl = 0.3)
    ps = PSpace1D(-0.5, 0.5, 500, 1)
    vs = nothing
    gas = Gas(Kn = 5e-3, K = 0.0)
    fw = function(x)
        prim = zeros(5)
        if x <= 0
            prim .= [1.0, 0.0, 1.0, 0.0, 1.0]
        else
            prim .= [1.0, 0.0, -1.0, 0.0, 2.0]
        end

        return prim_conserve(prim, ks.gas.γ)
    end
    ib = IB(fw, gas)
    ks = SolverSet(set, ps, vs, gas, ib)
    ctr, face = init_fvm(ks, ks.ps)
end

τ0 = vhs_collision_time(ctr[1].prim, ks.gas.μᵣ, ks.gas.ω)
tmax = 50 * τ0
t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(tmax ÷ dt)
res = zero(ctr[1].w)

@time @showprogress for iter = 1:nt
    #reconstruct!(ks, ctr)

    @inbounds @threads for i = 1:ks.ps.nx+1
        flux_gks!(
            face[i].fw,
            ctr[i-1].w .+ 0.5 .* ctr[i-1].sw .* ks.ps.dx[i-1],
            ctr[i].w .- 0.5 .* ctr[i].sw .* ks.ps.dx[i],
            ks.gas.K,
            ks.gas.γ,
            ks.gas.μᵣ,
            ks.gas.ω,
            dt,
            ks.ps.dx[i-1] / 2,
            ks.ps.dx[i] / 2,
            1.0,
            ctr[i-1].sw,
            ctr[i].sw,
        )
    end

    update!(ks, ctr, face, dt, res)

    #=global t += dt
    if abs(t - τ0) < dt
        @save "solns_t.jld2" ctr face
    elseif abs(t - 10 * τ0) < dt
        @save "solns_10t.jld2" ctr face
    end=#
end
#@save "solns_50t.jld2" ctr face

#=sol = zeros(ks.ps.nx, 5)
for i in axes(sol, 1)
    sol[i, :] .= ctr[i].prim
    sol[i, end] = 1 / sol[i, end]
end

plot(ks.ps.x[1:ks.ps.nx], sol)
=#

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
#-------------------------------------------------------------------------------
# LOWDER types
#-------------------------------------------------------------------------------

export AbstractModel, LOWDEROutput

abstract type AbstractModel end

"""

Define the type for LOWDER output.

    - 'status': solution status.
    - 'true_val': indicates whether the value 'f' is the true value of the
                    function 'fmin' for the point 'solution'.
    - 'iter': number of iterations.
    - 'nf': number of function evaluations.
    - 'index': index of the function in 'func_list'.
    - 'sample_radius': sample set radius 'δ'.
    - 'tr_radius': trust-region radius 'Δ'.
    - 'stationarity': stationarity measure 'π'.
    - 'f': objective function value.
    - 'solution': solution found by LOWDER.

"""
struct LOWDEROutput

    status        :: Symbol
    true_val      :: Bool
    iter          :: Int64
    nf            :: Int64
    index         :: Int64
    sample_radius :: Float64
    tr_radius     :: Float64
    stationarity  :: Float64
    f             :: Float64
    solution      :: Vector{Float64}

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
# @variables a, b, c

# ∨(a, b, c)

# struct Boolean
#     name
#     variable
# end

# function Boolean(x)
#     Boolean(x, Variable(Symbol(x)))
# end

# x = Boolean(:x)


text:
# Everything related to transformation on sparse matrices

export sparse_trafo_scale, sparse_trafo_swap, sparse_trafo_add_scaled,
       sparse_trafo_para_add_scaled, sparse_trafo_partial_dense,
       sparse_trafo_move_row

################################################################################
#
#  Constructors
#
################################################################################

function sparse_trafo_scale(i::Int, c::T) where {T}
  z = SparseTrafoElem{T, dense_matrix_type(T)}()
  z.type = 1
  z.i = i
  z.a = c
  return z
end

function sparse_trafo_swap(::Type{T}, i::Int, j::Int) where {T}
  z = SparseTrafoElem{T, dense_matrix_type(T)}()
  z.type = 2
  z.i = i
  z.j = j
  return z
end

function sparse_trafo_add_scaled(i::Int, j::Int, s::T) where {T}
  z = SparseTrafoElem{T, dense_matrix_type(T)}()
  z.type = 3
  z.i = i
  z.j = j
  z.a = s
  return z
end

function sparse_trafo_para_add_scaled(i::Int, j::Int, a::T, b::T, c::T, d::T) where {T}
  z = SparseTrafoElem{T, dense_matrix_type(T)}()
  z.type = 4
  z.i = i
  z.j = j
  z.a = a
  z.b = b
  z.c = c
  z.d = d
  return z
end

function sparse_trafo_partial_dense(i::Int, rows::UnitRange{Int}, cols::UnitRange{Int}, U::S) where {S}
  z = SparseTrafoElem{coefficient_type(S), S}()
  z.type = 5
  z.i = i
  z.rows = rows
  z.cols = cols
  z.U = U
  return z
end

# this is shorthand for the permutation matrix corresponding to
# (i i+1)(i+1 i+2)...(rows-1 rows)
function sparse_trafo_id(::Type{T}) where {T}
  z = SparseTrafoElem{T, dense_matrix_type(T)}()
  z.type = 7
  return z
end

function sparse_trafo_move_row(::Type{T}, i::Int, j::Int) where {T}
  z = SparseTrafoElem{T, dense_matrix_type(T)}()
  z.type = 6
  z.i = i
  z.j = j
  return z
end

function change_indices!(T::Array{SparseTrafoElem{S, SS}, 1}, st::Int, off::Int) where {S, SS}
  for t in T
    if t.type == 7
        continue
    end
    if t.i > st
      t.i += off
    end
    if t.type == 1 || t.type == 5
      continue
    end
    if t.j > st
      t.j += off
    end
  end
end

function max_index(T::Array{SparseTrafoElem{S, SS}, 1}) where {S, SS}
  i = 0
  for t in T
    if t.type == 7
        continue
    end
    i = max(i, t.i)
    if t.type == 1 || t.type == 5
      continue
    end
    i = max(i, t.j)
  end
  return i
end

################################################################################
#
#  Materialize
#
################################################################################

function matrix(R::Ring, T::SparseTrafoElem, n::Int)
  M = identity_matrix(R, n)
  t = T.type
  if t == 1
    i = T.i
    M[i, i] = T.c
  elseif t == 2
    i = T.i
    j = T.j
    M[i, i] = zero(R)
    M[j, j] = zero(R)
    M[i, j] = one(R)
    M[j, i] = one(R)
  elseif t == 3
    i = T.i
    j = T.j
    M[i, j] = T.a
  elseif t == 4
    i = T.i
    j = T.j
    M[i, i] = T.a
    M[i, j] = T.b
    M[j, j] = T.d
    M[j, i] = T.c
    return M
  elseif t == 6
    error("Not Implemented")
  end
end

################################################################################
#
#  Row scaling
#
################################################################################

@doc Markdown.doc"""
    scale_row!(A::SMat{T}, i::Int, c::T)

Multiply the $i$-th row of $A$ by $c$ inplace.
"""
function scale_row!(A::SMat{T}, i::Int, c::T) where T
  scale_row!(A[i], c)
  return A
end

################################################################################
#
#  Row swapping
#
################################################################################

@doc Markdown.doc"""
    swap_rows!(A::SMat{T}, i::Int, j::Int)

Swap the $i$-th and $j$-th row of $A$ inplace.
"""
function swap_rows!(A::SMat{T}, i::Int, j::Int) where T
  A[i], A[j] = A[j], A[i]
  return A
end

################################################################################
#
#  Row inversion
#
################################################################################

@doc Markdown.doc"""
    reverse_rows!(A::SMat)

Inplace inversion of the rows of $A$.
"""
function reverse_rows!(A::SMat{T}) where T
  for i = 1:div(A.r, 2)
    A[i], A[A.r+1-i] = A[A.r+1-i], A[i]
  end
  return A
end

################################################################################
#
#  Column swapping
#
################################################################################

@doc Markdown.doc"""
    swap_cols!(A::SMat, i::Int, j::Int)

Swap the $i$-th and $j$-th column of $A$ inplace.
"""
function swap_cols!(A::SMat{T}, i::Int, j::Int) where T
  @assert 1 <= i <= ncols(A) && 1 <= j <= ncols(A)

  if i == j
    return A
  end

  for r in A.rows
    if i in r.pos
      i_i = findfirst(r.pos, i)
      val_i = r.values[i_i]
      if j in r.pos
        i_j = findfirst(r.pos, j)
        val_j = r.values[i_j]

        r.values[i_i], r.values[i_j] = r.values[i_j], r.values[i_i]
      else
        t = r.values[i_i]
        deleteat!(r.pos, i_i)
        deleteat!(r.values, i_i)
        k = searchsortedfirst(r.pos, j)
        insert!(r.pos, k, j)
        insert!(r.values, k, t)
      end
    else
      if j in r.pos
        i_j = findfirst(r.pos, j)
        val_j = r.values[i_j]

        t = r.values[i_j]
        deleteat!(r.pos, i_j)
        deleteat!(r.values, i_j)
        k = searchsortedfirst(r.pos, i)
        insert!(r.pos, k, i)
        insert!(r.values, k, t)
      end
    end
  end
  return A
end

################################################################################
#
#  Addition of scaled rows
#
################################################################################
# rows j -> row i*c + row j
@doc Markdown.doc"""
    add_scaled_row!(A::SMat{T}, i::Int, j::Int, c::T)

Add $c$ times the $i$-th row to the $j$-th row of $A$ inplace, that is,
$A_j \rightarrow A_j + c \cdot A_i$, where $(A_i)_i$ denote the rows of $A$.
"""
function add_scaled_row!(A::SMat{T}, i::Int, j::Int, c::T) where T
  A.nnz = A.nnz - length(A[j])
  A.rows[j] = add_scaled_row(A[i], A[j], c)
  A.nnz = A.nnz + length(A[j])
  return A
end

################################################################################
#
#  Addition of scaled cols
#
################################################################################

@doc Markdown.doc"""
    add_scaled_col!(A::SMat{T}, i::Int, j::Int, c::T)

Add $c$ times the $i$-th column to the $j$-th column of $A$ inplace, that is,
$A_j \rightarrow A_j + c \cdot A_i$, where $(A_i)_i$ denote the columns of $A$.
"""
function add_scaled_col!(A::SMat{T}, i::Int, j::Int, c::T) where T
  @assert c != 0

  @assert 1 <= i <= ncols(A) && 1 <= j <= ncols(A)

  for r in A.rows
    if i in r.pos
      i_i = findfirst(r.pos, i)
      val_i = r.values[i_i]
      if j in r.pos
        i_j = findfirst(r.pos, j)
        val_j = r.values[i_j]

        r.values[i_j] += c*r.values[i_i]
      else
        k = searchsortedfirst(r.pos, j)
        insert!(r.pos, k, j)
        insert!(r.values, k, c*r.values[i_i])
      end
    end
  end
  return A
end

################################################################################
#
#  Elementary row transformation
#
################################################################################

@doc Markdown.doc"""
    transform_row!(A::SMat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)

Applies the transformation $(A_i, A_j) \rightarrow (aA_i + bA_j, cA_i + dA_j)$
to $A$, where $(A_i)_i$ are the rows of $A$.
"""
function transform_row!(A::SMat{T}, i::Int, j::Int, a::T, b::T, c::T, d::T) where T
  A.nnz = A.nnz - length(A[i]) - length(A[j])
  A.rows[i], A.rows[j] = transform_row(A[i], A[j], a, b, c, d)
  A.nnz = A.nnz + length(A[i]) + length(A[j])
  return A
end

@doc Markdown.doc"""
    transform_row(A::SRow{T}, B::SRow{T}, i::Int, j::Int, a::T, b::T, c::T, d::T)

Returns the tuple $(aA + bB, cA + dB)$.
"""
function transform_row(Ai::SRow{T}, Aj::SRow{T}, a::T, b::T, c::T, d::T) where T
  sr = SRow{T}()
  tr = SRow{T}()
  pi = 1
  pj = 1
  while pi <= length(Ai) && pj <= length(Aj)
    if Ai.pos[pi] < Aj.pos[pj]
      if a != 0
        push!(sr.pos, Ai.pos[pi])
        push!(sr.values, a*Ai.values[pi])
      end
      if c != 0
        push!(tr.pos, Ai.pos[pi])
        push!(tr.values, c*Ai.values[pi])
      end
      pi += 1
    elseif Ai.pos[pi] > Aj.pos[pj]
      if b != 0
        push!(sr.pos, Aj.pos[pj])
        push!(sr.values, b*Aj.values[pj])
      end
      if d != 0
        push!(tr.pos, Aj.pos[pj])
        push!(tr.values, d*Aj.values[pj])
      end
      pj += 1
    else
      m = a*Ai.values[pi] + b*Aj.values[pj]
      n = c*Ai.values[pi] + d*Aj.values[pj]
      if m != 0
        push!(sr.pos, Ai.pos[pi])
        push!(sr.values, m)
      end
      if n != 0
        push!(tr.pos, Ai.pos[pi])
        push!(tr.values, n)
      end
      pi += 1
      pj += 1
    end
  end
  while pi <= length(Ai.pos)
    if a != 0
      push!(sr.pos, Ai.pos[pi])
      push!(sr.values, a*Ai.values[pi])
    end
    if c != 0
      push!(tr.pos, Ai.pos[pi])
      push!(tr.values, c*Ai.values[pi])
    end
    pi += 1
  end
  while pj <= length(Aj.pos)
    if b != 0
      push!(sr.pos, Aj.pos[pj])
      push!(sr.values, b*Aj.values[pj])
    end
    if d != 0
      push!(tr.pos, Aj.pos[pj])
      push!(tr.values, d*Aj.values[pj])
    end
    pj += 1
  end

  return sr, tr
end

################################################################################
#
#  Inverses of elementary matrices
#
################################################################################

function inv(t::SparseTrafoElem{T, TT}) where {T, TT}
  i = t.type
  if i == 1
    fl, d = divides(one(parent(t.a)), t.a)
    @assert fl
    return sparse_trafo_scale(t.i, d)
  elseif i == 2
    return t
  elseif i == 3
    return sparse_trafo_add_scaled(t.i, t.j, -t.a)
  elseif i == 4
    e = t.a * t.d - t.b * t.c
    fl, e = divides(one(parent(e)), e)
    return sparse_trafo_para_add_scaled(t.i, t.j, e * t.d, e * -t.b, e * -t.c, e * t.a)
  elseif i == 5
    return sparse_trafo_partial_dense(t.i, t.rows, t.cols, inv(t.U))
  elseif i == 6
    return sparse_trafo_move_row(T, t.j, t.i)
  elseif i == 7
    return t
  end
end

################################################################################
#
#  String I/O
#
################################################################################

function Base.show(io::IO, t::SparseTrafoElem)
  print(io, "Sparse transformation: ")
  i = t.type
  if i == 1
    print(io, "Scale ", t.i, " by ", t.a)
  elseif i == 2
    print(io, "Swap ", t.i, " and ", t.j)
  elseif i == 3
    print(io, "Scale ", t.i, " by ", t.a, " and add to ", t.j)
  elseif i == 4
    print(io, "Transform", t.i, ", ", t.j, " by [", t.a, " ", t.b, " ", t.c, " ", t.d, "]")
  elseif i == 5
    print(io, "Dense ", nrows(t.U), "x", nrows(t.U), " at ", t.i)
  elseif i == 6
    print(io, "Move ", t.i, " to ", t.j)
  end
end

################################################################################
#
#  Application of a transformation on the left side of a sparse matrix
#
################################################################################

# The following function do not update the number of nonzero entries
# properly

function apply_left!(A::SMat{T}, t::SparseTrafoElem{T, S}) where {T, S}
  i = t.type
  if i == 1
    scale_row!(A, t.i, t.a)
  elseif i == 2
    swap_rows!(A, t.i, t.j)
  elseif i == 3
    add_scaled_row!(A, t.i, t.j, t.a)
  elseif i == 4
    transform_row!(A, t.i, t.j, t.a, t.b, t.c, t.d)
  elseif i == 5
    R = base_ring(A)
    k = t.i
    h = sub(A, t.rows, t.cols)
    A.nnz -= sum(length(A[i].pos) for i in t.rows)
    deleteat!(A.rows, t.rows)
    A.r -= length(t.rows)
    @assert length(A.rows) == A.r
    # make dense matrix
    hdense = matrix(h)
    hdense = t.U * hdense

    htransformed = sparse_matrix(hdense)
    for r in htransformed
      for j in 1:length(r.pos)
        r.pos[j] = r.pos[j] + (k - 1)
      end
      push!(A, r)
    end
  elseif i == 6
    x = A.rows
    if t.i <= t.j
      for j in t.i:t.j-1
        r = x[j]
        x[j] = x[j + 1]
        x[j + 1] = r
      end
    else
      for j in t.i:-1:t.j+1
        r = x[j]
        x[j] = x[j - 1]
        x[j - 1] = r
      end
    end
    #deleteat!(A.rows, t.i)
    #push!(A.rows, sparse_row(base_ring(A)))
  else
    error("Wrong type")
  end
  return A
end

################################################################################
#
#  Application of a transformation on the right hand side of a row
#
################################################################################

function apply_right!(x::Vector{T}, t::SparseTrafoElem{T, S}) where {T, S}
  i = t.type
  if i == 1
    x[t.i] = x[t.i] * t.a
  elseif i == 2
    r = x[t.i]
    x[t.i] = x[t.j]
    x[t.j] = r
  elseif i == 3
    x[t.i] = x[t.i] + x[t.j]*t.a
  elseif i == 4
    r = t.a * x[t.i] + t.c * x[t.j]
    s = t.b * x[t.i] + t.d * x[t.j]
    x[t.i] = r
    x[t.j] = s
  elseif i == 5
    sm = matrix(parent(x[1]), 1, nrows(t.U), x[t.rows])
    #println("s :$s")
    sm = sm*t.U
    for (i, j) in zip(t.rows,1:nrows(t.U))
      x[i] = sm[1, j]
    end
  elseif i == 6
    if t.j > t.i
      for j in t.j:-1:t.i+1
        r = x[j]
        x[j] = x[j - 1]
        x[j - 1] = r
      end
    else
      for j in t.j:t.i-1
        r = x[j]
        x[j] = x[j + 1]
        x[j + 1] = r
      end
    end
  end
  return x
end

################################################################################
#
#  Application to the left of an Array of ideals
#
################################################################################

#function apply_left!(x::Vector{NfOrdFracIdl}, y::TrafoSwap)
#  x[y.i], x[y.j] = x[y.j], x[y.i]
#end
#
#function apply_left!(x::Vector{NfOrdFracIdl}, y::TrafoAddScaled)
#  x[y.j] = x[y.j] * x[y.i]^Int(y.s)
#end
#
#function apply_left!(x::Vector{NfOrdFracIdl}, y::TrafoPartialDense)
#  z = view(deepcopy(x), y.cols)
#  xx = view(x, y.cols)
#  for i in 1:nrows(y.U)  ## use power product instead
#    xx[i] = z[1]^Int(y.U[i, 1])
#    for j in 2:ncols(y.U)
#      xx[i] *= z[j]^Int(y.U[i, j])
#    end
#  end
#end


text:
struct PieceWiseNormal <: Gen.Distribution{Float64} end
const piecewise_normal = PieceWiseNormal()

# piece 1 has interval (-Inf, boundaries[1])
# piece 2 has interval (boundaries[1], boundaries[2])
# ..
# piece n has interval (boundaries[n-1], boundaries[n])
# piece n+1 has interval (boundaries[n], Inf)

function Gen.random(::PieceWiseNormal, probabilities::Vector{Float64},
                    mus::Vector{Float64}, stds::Vector{Float64},
                    boundaries::Vector{Float64})
    n_segments = length(probabilities)
    @assert n_segments == length(mus)
    @assert n_segments == length(stds)
    @assert n_segments == length(boundaries)+1
    idx = categorical(probabilities)
    dist = Distributions.TruncatedNormal(
        mus[idx],
        stds[idx],
        idx == 1 ? -Inf : boundaries[idx-1], idx == n_segments ? Inf : boundaries[idx])
    rand(dist)
end

function Gen.logpdf(::PieceWiseNormal, x::Float64, probabilities::Vector{Float64},
                    mus::Vector{Float64}, stds::Vector{Float64},
                    boundaries::Vector{Float64})
    n_segments = length(probabilities)
    @assert n_segments == length(mus)
    @assert n_segments == length(stds)
    @assert n_segments == length(boundaries)+1
    log_probs = Vector{Float64}(undef, n_segments)
    for idx=1:length(probabilities)
        dist = Distributions.TruncatedNormal(
            mus[idx],
            stds[idx],
            idx == 1 ? -Inf : boundaries[idx-1], idx == n_segments ? Inf : boundaries[idx])
        log_probs[idx] = log(probabilities[idx]) + Distributions.logpdf(dist, x)
    end
    logsumexp(log_probs)
end

Gen.has_output_grad(::PieceWiseNormal) = false
Gen.has_argument_grads(::PieceWiseNormal) = (false, false, false, false)


text:
module CellFieldsMocksTests

using Gridap
using ..CellFieldsMocks
using ..CellValuesMocks

l = 10
cf = IterCellFieldMock(2,Int,l)

p1 = Point(1,2)
p2 = Point(0,3)
p3 = Point(7,4)
p4 = Point(8,1)
p = [p1,p2,p3,p4]
cp = TestIterCellValue(p,l)

f,_ = iterate(cf)
fg = gradient(f)
r = evaluate(f,p)
rg = evaluate(fg,p)
v = [r for i in 1:l]
g = [rg for i in 1:l]

test_iter_cell_field(cf,cp,v,g)
test_iter_cell_field_without_grad(cf,cp,v)

cf = IndexCellFieldMock(2,Int,l)
cp = TestIndexCellValue(p,l)
test_index_cell_field(cf,cp,v,g)
test_index_cell_field_without_grad(cf,cp,v)

cf = IterCellBasisMock(2,Int,l)
cp = TestIterCellValue(p,l)

f,_ = iterate(cf)
fg = gradient(f)
r = evaluate(f,p)
rg = evaluate(fg,p)
v = [r for i in 1:l]
g = [rg for i in 1:l]

test_iter_cell_basis(cf,cp,v,g)
test_iter_cell_basis_without_grad(cf,cp,v)

cf = IndexCellBasisMock(2,Int,l)
cp = TestIndexCellValue(p,l)
test_index_cell_basis(cf,cp,v,g)
test_index_cell_basis_without_grad(cf,cp,v)

cf = IterCellGeomapMock(2,Int,l)
cp = TestIndexCellValue(p,l)

f,_ = iterate(cf)
fg = gradient(f)
r = evaluate(f,p)
rg = evaluate(fg,p)
v = [r for i in 1:l]
g = [rg for i in 1:l]

test_iter_cell_field(cf,cp,v,g)

cf = IndexCellGeomapMock(2,Int,l)
test_index_cell_field(cf,cp,v,g)

end # module


text:
"""
    Kucoin.list_accounts(api_data[; currency[, type]]) -> JSON3.Array{JSON3.Object}

Get a list of accounts. See Kucoin official [docs](https://docs.kucoin.com/#list-accounts)
for more details.

# Arguments
- `api_data::UserApiData`: holds user API data susch as api key, secret, passphrase, etc.

# Keywords
- `currency::String`: [Optional] a unique currency kucoin code
- `type::String`: [Optional] an account type, the valid values are: `"main"`, `"trade"`,
`"margin"` or `"pool"`

# Returns
- `JSON3.Array{JSON3.Object}`
```json
[
  {
    "id": "5bd6e9286d99522a52e458de",  //accountId
    "currency": "BTC",  //Currency
    "type": "main",  //Account type, including main and trade
    "balance": "237582.04299",  //Total assets of a currency
    "available": "237582.032",  //Available assets of a currency
    "holds": "0.01099"  //Hold assets of a currency
  },
  {
    "id": "5bd6e9216d99522a52e458d6",
    "currency": "BTC",
    "type": "trade",
    "balance": "1234356",
    "available": "1234356",
    "holds": "0"
}]
```
"""
function list_accounts(
    api_data::ApiData;
    currency::Union{String,Nothing}=nothing,
    type::Union{String,Nothing}=nothing,
)::JSON3.Array{JSON3.Object}
    return _handle(
        _kucoin_request(
            api_data, _HTTP_METHOD_GET, "/api/v1/accounts"; currency=currency, type=type
        ),
    )
end

"""
    Kucoin.get_trade_fees(api_data; symbols) -> JSON3.Array{JSON3.Object}

Get the actual fee rate of the trading pair. See Kucoin official
[docs](https://docs.kucoin.com/#actual-fee-rate-of-the-trading-pair) for more details.

# Arguments
- `api_data::UserApiData`: holds user API data susch as api key, secret, passphrase, etc.
- `symbols::AbstractVector{String}`: Trading pairs (optional, you can inquire fee rates of
10 trading pairs each time at most)

# Returns
- `JSON3.Array{JSON3.Object}`
```json
[
  {
    "symbol": "BTC-USDT",
    "takerFeeRate": "0.001",
    "makerFeeRate": "0.001"
  },
  {
    "symbol": "KCS-USDT",
    "takerFeeRate": "0.002",
    "makerFeeRate": "0.0005"
}
```
"""
function get_trade_fees(
    api_data::ApiData; symbols::AbstractVector{String}
)::JSON3.Array{JSON3.Object}
    return _handle(
        _kucoin_request(
            api_data, _HTTP_METHOD_GET, "/api/v1/trade-fees"; symbols=join(symbols, ",")
        ),
    )
end


text:
using CImGui
using CImGui.LibCImGui
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using GLFW, ModernGL
using CSyntax
using Printf
@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end
# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)
# first window and its context
window1 = GLFW.CreateWindow(1280, 720, "Demo")
@assert window1 != C_NULL
GLFW.MakeContextCurrent(window1)
ctx1 = CImGui.CreateContext()
CImGui.StyleColorsDark()
ImGui_ImplGlfw_InitForOpenGL(window1, true)
ImGui_ImplOpenGL3_Init(glsl_version)
# second window and its context
window2 = GLFW.CreateWindow(800, 600, "Demo2")
@assert window2 != C_NULL
GLFW.MakeContextCurrent(window2)
ctx2 = CImGui.CreateContext()
CImGui.SetCurrentContext(ctx2)  # swap context
CImGui.StyleColorsLight()
ImGui_ImplGlfw_InitForOpenGL(window2, true)
ImGui_ImplOpenGL3_Init(glsl_version)
show_demo_window = true
show_another_window = false
counter = 0
f = Cfloat(0.0)
clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
window1_open = true
window2_open = true
while window1_open || window2_open
    # oh my global scope
    global show_demo_window; global show_another_window
    global counter; global f; global clear_color;
    global window1_open; global window2_open;
    GLFW.PollEvents()
    GLFW.WindowShouldClose(window1) && (window1_open = false;)
    GLFW.WindowShouldClose(window2) && (window2_open = false;)
    if window1_open
        GLFW.MakeContextCurrent(window1)
        CImGui.SetCurrentContext(ctx1)
        GLFW.PollEvents()
        # append ctx1's frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame(window1)
        CImGui.NewFrame()
        # show the big demo window
        show_demo_window && @c igShowDemoWindow(&show_demo_window)
        # render window1
        CImGui.Render()
        GLFW.MakeContextCurrent(window1)
        display_w, display_h = GLFW.GetFramebufferSize(window1)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(igGetDrawData())
        GLFW.SwapBuffers(window1)
    else
        GLFW.HideWindow(window1)
    end
    # window 2
    if window2_open
        GLFW.MakeContextCurrent(window2)
        CImGui.SetCurrentContext(ctx2)
        GLFW.PollEvents()
        # append ctx2's frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame(window2)
        CImGui.NewFrame()
        # show a simple window that we create ourselves.
        # we use a Begin/End pair to created a named window.
        CImGui.Begin("Hello, world!")  # create a window called "Hello, world!" and append into it.
        CImGui.Text("This is some useful text.")  # display some text
        @c CImGui.Checkbox("Demo Window", &show_demo_window)  # edit bools storing our window open/close state
        @c CImGui.Checkbox("Another Window", &show_another_window)
        @c CImGui.SliderFloat("float", &f, 0, 1)  # edit 1 float using a slider from 0 to 1
        CImGui.ColorEdit3("clear color", clear_color)  # edit 3 floats representing a color
        CImGui.Button("Button") && (counter += 1)
        CImGui.SameLine()
        CImGui.Text("counter = $counter")
        CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / ImGuiIO_Get_Framerate(CImGui.GetIO()), ImGuiIO_Get_Framerate(CImGui.GetIO())))
        CImGui.End()
        CImGui.Render()
        GLFW.MakeContextCurrent(window2)
        display_w, display_h = GLFW.GetFramebufferSize(window2)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(igGetDrawData())
        GLFW.MakeContextCurrent(window2)
        GLFW.SwapBuffers(window2)
    else
        GLFW.HideWindow(window2)
    end
end
# cleanup
ImGui_ImplOpenGL3_Shutdown()
ImGui_ImplGlfw_Shutdown()
CImGui.DestroyContext(ctx1)
CImGui.DestroyContext(ctx2)
GLFW.DestroyWindow(window1)
GLFW.DestroyWindow(window2)


text:
using LinearAlgebra

# Define function and gradient
f(x) = x[1] + exp(x[2]-x[1]);
grad_f(z) = [1-exp(z[2]-z[1]), exp(z[2]-z[1])];
z = [1, 2];
f_hat(x) = f(z) + grad_f(z)'*(x-z);

# Compare f and f_hat for some specific x’s
f([1,2]),  f_hat([1,2])

f([0.96,1.98]),  f_hat([0.96,1.98])

f([0.96,1.98]),  f_hat([0.96,1.98])

max_stars_repo_name: lynnlangit/julia-linear-algebra

max_issues_repo_name: lynnlangit/julia-linear-algebra

max_forks_repo_name: lynnlangit/julia-linear-algebra


text:
module MvDecisionTrees

export
    LossFunction,
    LossFunction_MSE,
    LossFunction_LOGL,
    LossFunction_LOGL_MEAN_SUBTRACTED,

    Leaf,
    MeanVecLeaf,
    CovLeaf,
    MvNormLeaf,

    InternalNode,
    Ensemble,
    BuildTreeParameters,

    depth,
    print_tree,
    apply_tree,
    apply_forest!,

    build_leaf,
    build_tree,
    build_forest,
    build_tree_parameters,

    get_split_index_set

import Distributions: MvNormal, fit_mle
import StatsBase: cov

if VERSION.minor ≥ 4
    typealias Range1{Int} Range{Int}
    _int(x) = round(Int, x, RoundNearestTiesUp)
else
    _int(x) = int(x)
end

##############################################################################
# TreeData

type TreeData{T<:AbstractFloat, U<:Real}
    # learning data
    Y::Matrix{T} # [n_observations × n_rows]
    X::Matrix{U} # [n_rows × n_predictors]
    assignment::Vector{Int} # [n_rows]

    # preallocated learning memory
    μ_l::Vector{T} # n
    μ_r::Vector{T} # n
    Σ_l::Matrix{T} # n×n
    Σ_r::Matrix{T} # n×n
    subfeature_indeces::Vector{Int}
end

Base.size(d::TreeData) = (size(d.Y,2), size(d.Y,1), size(d.X,2)) # (n, o, p)  n_rows, n_obs, n_predictors
get_num_samples(d::TreeData) = size(d.Y,2)
get_num_observations(d::TreeData) = size(d.Y,1)
function get_id_count(data::TreeData, id::Int)
    id_count = 0
    for i = 1 : size(data.X, 1)
        if data.assignment[i] == id
            id_count += 1
        end
    end
    id_count
end
function get_nth_valid_id(data::TreeData, id::Int, n::Int)
    id_count = 0
    for i = 1 : size(data.X, 1)
        if data.assignment[i] == id
            id_count += 1
            if id_count == n
                return i
            end
        end
    end
    return 0
end
function get_left_right_counts{T,U}(data::TreeData{T,U}, predictor_index::Int, id::Int, threshold::T)

    num_left = 0
    num_right = 0

    for i = 1 : get_num_samples(data)
        if data.assignment[i] == id
            if data.X[i,predictor_index] < threshold
                num_left += 1
            else
                num_right += 1
            end
        end
    end

    (num_left, num_right)
end

##############################################################################
# LossFunction

abstract LossFunction
type LossFunction_MSE <: LossFunction end # mean squared error loss function
type LossFunction_LOGL <: LossFunction end # log likelihood loss function
type LossFunction_LOGL_MEAN_SUBTRACTED <: LossFunction end # log likelihood loss function where labels already has predicted mean subtracted

function loss{T<:AbstractFloat, U<:Real}(::Type{LossFunction_MSE},
    data::TreeData{T,U},
    predictor_index::Int,
    assignment_id::Int,
    threshold::T, # split point
    )

    # mse loss for a vector is tr{E{(y-p)(y-p)ᵀ}}
    # where y is the label
    # and p is the prediction - the average of labels on that side of the threshold

    Σ_l = data.Σ_l
    Σ_r = data.Σ_r

    n, o, p = size(data) # n_rows, n_obs, n_predictors
    mse = 0.0
    nl, nr = get_left_right_counts(data, predictor_index, assignment_id, threshold)

    for i = 1 : n
        if data.assignment[i] == assignment_id
            s_l, s_r = 0.0, 0.0
            for j = 1 : o
                l = data.Y[j,i]
                mse += l*l
                if data.X[j, predictor_index] < threshold
                    s_l += l
                else
                    s_r += l
                end
            end
            mse -= s_l*s_l / nl
            mse -= s_r*s_r / nr
        end
    end

    return (-mse, nl, nr)
end
function loss{T<:AbstractFloat, U<:Real}(::Type{LossFunction_MSE},
    data::TreeData{T,U},
    assignment_id::Int,
    )

    # mse loss for a vector is tr{E{(y-p)(y-p)ᵀ}}
    # where y is the label
    # and p is the prediction - the average of labels on that side of the threshold

    Σ = data.Σ_l
    n, o, p = size(data) # n_samples, n_outcomes, n_predictors
    mse = 0.0

    id_count = get_id_count(data, assignment_id)

    for i = 1 : n
        if data.assignment[i] == assignment_id
            s = 0.0
            for j = 1 : o
                l = data.Y[j,i]
                mse += l*l
                s += l
            end
            mse -= s*s / id_count
        end
    end

    return -mse
end

function loss{T<:AbstractFloat, S<:Real}(::Type{LossFunction_LOGL},
    data::TreeData{T,S},
    predictor_index::Int,
    assignment_id::Int,
    threshold::T, # split point
    )

    μ_l = data.μ_l
    μ_r = data.μ_r
    Σ_l = data.Σ_l
    Σ_r = data.Σ_r

    n, o, p = size(data)
    nl, nr = 0, 0

    # compute the mean
    for i = 1 : o
        μ_l[i] = 0.0
        μ_r[i] = 0.0
    end
    for i = 1 : n
        if data.assignment[i] == assignment_id
            if data.X[i,predictor_index] < threshold
                nl += 1
                for a = 1 : o
                    μ_l[a] += data.Y[a,i]
                end
            else
                nr += 1
                for a = 1 : o
                    μ_r[a] += data.Y[a,i]
                end
            end
        end
    end
    if nl < 2 || nr < 2
        # NOTE(tim): force at least two samples on each side to ensure proper covariance matrix
        return (-Inf, nl, nr)
    end
    for i = 1 : o
        μ_l[i] /= nl
        μ_r[i] /= nr
    end

    # solve for upper triangle of symmetric covariance matrix
    for i = 1 : o*o
        Σ_l[i] = 0.0
        Σ_r[i] = 0.0
    end
    for i = 1 : n
        if data.assignment[i] == assignment_id
            if data.X[i,predictor_index] < threshold
                for a = 1 : o
                    l = data.Y[a,i] - μ_l[a]
                    for b = a : o
                        Σ_l[a,b] += l*(data.Y[b,i] - μ_l[b])
                    end
                end
            else
                for a = 1 : o
                    l = data.Y[a,i] - μ_r[a]
                    for b = a : o
                        Σ_r[a,b] += l*(data.Y[b,i] - μ_r[b])
                    end
                end
            end
        end
    end
    for i = 1 : o*o
        Σ_l[i] /= (nl-1)
        Σ_r[i] /= (nr-1)
    end
    for a = 2:o
        # copy over symmetric component
        for b = 1:a-1
            Σ_l[a,b] = Σ_l[b,a]
            Σ_r[a,b] = Σ_r[b,a]
        end
    end

    logdetΣ_l = log(max(_fast_det(Σ_l), 1e-20))
    logdetΣ_r = log(max(_fast_det(Σ_r), 1e-20))

    # compute log likelihood score under gaussian distribution
    # NOTE: does not compute the true logl likelihood
    #       computation ignores constant components
    logl = 0.0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            if data.X[i] < threshold
                logl -= logdetΣ_l
                for a = 1 : o
                    subval = 0.0
                    for b = 1 : o
                        subval += Σ_l[a,b]*(data.Y[b,i] - μ_l[b])
                    end
                    logl -= subval*(data.Y[a,i] - μ_l[a])
                end
            else
                logl -= logdetΣ_r
                for a = 1 : o
                    subval = 0.0
                    for b = 1 : o
                        subval += Σ_r[a,b]*(data.Y[b,i] - μ_r[b])
                    end
                    logl -= subval*(data.Y[a,i] - μ_r[a])
                end
            end
        end
    end

    return (logl, nl, nr)
end
function loss{T<:AbstractFloat, S<:Real}(::Type{LossFunction_LOGL},
    data::TreeData{T,S},
    assignment_id::Int,
    )

    μ = data.μ_l
    Σ = data.Σ_l

    n, o, p = size(data)
    # n, m = size(data.Y)
    m_assigned = 0

    # compute the sample mean
    for i = 1 : o
        μ[i] = 0.0
    end
    for i = 1 : n
        if data.assignment[i] == assignment_id
            m_assigned += 1
            for a = 1 : o
                μ[a] += data.Y[a,i]
            end
        end
    end
    if m_assigned ≤ 2
        # insufficient entries to split
        return -Inf
    end
    for i = 1 : o
        μ[i] /= m_assigned
    end

    # compute covariance matrix

    for i = 1 : o*o
        Σ[i] = 0.0
    end
    for i = 1 : n
        if data.assignment[i] == assignment_id
            for a = 1 : o
                l = data.Y[a,i] - μ[a]
                for b = a : o
                    Σ[a,b] += l*(data.Y[b,i] - μ[b])
                end
            end
        end
    end
    for i = 1 : o*o
        Σ[i] /= (m_assigned-1)
    end
    for a = 2:o
        # copy over symmetric component
        for b = 1:a-1
            Σ[a,b] = Σ[b,a]
        end
    end

    logdetΣ = log(max(_fast_det(Σ), 1e-20))

    # compute log likelihood score under gaussian distribution
    # NOTE: does not compute the true logl likelihood
    #       computation ignores constant components
    logl = 0.0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            logl -= logdetΣ
            for a = 1 : o
                subval = 0.0
                for b = 1 : o
                    subval += Σ[a,b]*(data.Y[b,i] - μ[b])
                end
                logl -= subval*(data.Y[a,i] - μ[a])
            end
        end
    end
    logl
end

function loss{T<:AbstractFloat, S<:Real}(::Type{LossFunction_LOGL_MEAN_SUBTRACTED},
    data::TreeData{T,S},
    predictor_index::Int,
    assignment_id::Int,
    threshold::T, # split point
    )

    # data.Y is n_samples × n_outcomes and contains the distance from the mean, x-μ
    # we want to predict the covariance matrix Σ, which is determined by the covariance
    # of the samples in a leaf

    Σ_l = data.Σ_l
    Σ_r = data.Σ_r

    n, o, p = size(data)

    nl, nr = 0, 0

    # solve for upper triangle of symmetric matrix
    for i = 1 : o*o
        Σ_l[i] = 0.0
        Σ_r[i] = 0.0
    end

    for i = 1 : n
        if data.assignment[i] == assignment_id
            if data.X[i,predictor_index] < threshold
                for a = 1 : o
                    l = data.Y[a,i]
                    for b = a : o
                        Σ_l[a,b] += l*data.Y[b,i]
                    end
                end
                nl += 1
            else
                for a = 1 : o
                    l = data.Y[a,i]
                    for b = a : o
                        Σ_r[a,b] += l*data.Y[b,i]
                    end
                end
                nr += 1
            end
        end
    end
    for i = 1 : o*o
        Σ_l[i] /= (nl-1)
        Σ_r[i] /= (nr-1)
    end

    # NOTE(tim): force at least two samples on each side to ensure proper covariance matrix
    if nl < 2 || nr < 2
        return (-Inf, nl, nr)
    end

    # copy over symmetric component
    for a = 2:o
        for b = 1:a-1
            Σ_l[a,b] = Σ_l[b,a]
            Σ_r[a,b] = Σ_r[b,a]
        end
    end

    logdetΣ_l = log(max(_fast_det(Σ_l), 1e-20))
    logdetΣ_r = log(max(_fast_det(Σ_r), 1e-20))

    # compute log likelihood (we seek to maximize it)
    # NOTE: ignores constant components
    logl = 0.0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            if data.X[i] < threshold
                logl -= logdetΣ_l
                for a = 1 : o
                    subval = 0.0
                    for b = 1 : o
                        subval += Σ_l[a,b]*data.Y[b,i]
                    end
                    logl -= subval*data.Y[a,i]
                end
            else
                logl -= logdetΣ_r
                for a = 1 : o
                    subval = 0.0
                    for b = 1 : o
                        subval += Σ_r[a,b]*data.Y[b,i]
                    end
                    logl -= subval*data.Y[a,i]
                end
            end
        end
    end

    return (logl, nl, nr)
end
function loss{T<:AbstractFloat, S<:Real}(::Type{LossFunction_LOGL_MEAN_SUBTRACTED},
    data::TreeData{T,S},
    assignment_id::Int,
    )

    # data.Y is n_samples × n_outcomes and contains the distance from the mean, x-μ
    # we want to predict the covariance matrix Σ, which is determined by the covariance
    # of the samples in a leaf

    Σ = data.Σ_l
    n, o, p = size(data)

    for i = 1 : o*o
        Σ[i] = 0
    end

    # solve for upper triangle of symmetric matrix

    m_assigned = 0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            m_assigned += 1
            for a = 1 : o
                l = data.Y[a,i]
                for b = a : o
                    Σ[a,b] += l*data.Y[b,i]
                end
            end
        end
    end
    if m_assigned ≤ 2
        return -Inf
    end

    for i = 1 : o*o
        Σ[i] /= (m_assigned-1)
    end

    # copy over symmetric component
    for a = 2:o
        for b = 1:a-1
            Σ[a,b] = Σ[b,a]
        end
    end

    logdetΣ = log(max(_fast_det(Σ), 1e-20))

    # compute log likelihood score under gaussian distribution
    # NOTE: does not compute the true logl likelihood,
    #       computation ignores constant components
    logl = 0.0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            logl -= logdetΣ
            for a = 1 : o
                subval = 0.0
                for b = 1 : o
                    subval += Σ[a,b]*data.Y[b,i]
                end
                logl -= subval*data.Y[a,i]
            end
        end
    end
    logl
end

##############################################################################
# Leaf

abstract Leaf
immutable MeanVecLeaf{T<:AbstractFloat} <: Leaf
    μ::Vector{T}
end
immutable CovLeaf{T<:AbstractFloat} <: Leaf
    Σ::Matrix{T}
end
immutable MvNormLeaf <: Leaf
    m::MvNormal
end

function build_leaf{T<:AbstractFloat, S<:Real}(::Type{MeanVecLeaf},
    data::TreeData{T,S},
    assignment_id::Int,
    )

    #=
    Construct a leaf whose payload is the mean vector of samples within it
    =#

    n, o, p = size(data) # n_rows, n_obs, n_predictors
    leaf_mean = zeros(T, o)

    n_for_id = 0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            n_for_id += 1
            for j = 1 : o
                leaf_mean[j] += data.Y[j,i]
            end
        end
    end
    for j = 1 : o
        leaf_mean[j] /= n_for_id
    end

    MeanVecLeaf(leaf_mean)
end
function build_leaf{T<:AbstractFloat, S<:Real}(::Type{CovLeaf},
    data::TreeData{T,S},
    assignment_id::Int,
    )

    #=
    Construct a leaf whose payload is the covariance matrix for all samples within it
    =#
    o = get_num_observations(data)
    CovLeaf(_calc_covariance!(zeros(T, o, o), data, assignment_id))
end
function build_leaf{T<:AbstractFloat, S<:Real}(::Type{MvNormLeaf},
    data::TreeData{T,S},
    assignment_id::Int,
    )

    #=
    Construct a leaf whose payload is the mvnormal for all samples within it
    =#

    n, o, p = size(data)
    leaf_mean = zeros(T, o)

    n_for_id = 0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            n_for_id += 1
            for j = 1 : o
                leaf_mean[j] += data.Y[j,i]
            end
        end
    end
    for j = 1 : o
        leaf_mean[j] /= n_for_id
    end

    # calc covariance

    Σ = zeros(T, o, o)

    for i = 1 : n
        if data.assignment[i] == assignment_id
            for a = 1 : o
                l = data.Y[a,i] - leaf_mean[a]
                for b = a : o
                    Σ[a,b] += l*(data.Y[b,i]-leaf_mean[b])
                end
            end
        end
    end

    if n_for_id < 2
        warn("MvDecisionTrees.build_leaf{MvNormLeaf}: Not enough samples ($n_for_id) to build a covariance matrix")
        for i = 1 : o
            Σ[i,i] = 1.0
        end
    else
        for i = 1 : o*o
            Σ[i] /= (n_for_id-1)
        end
        _diagonal_shrinkage!(Σ)
        for a = 2:o
            for b = 1: a-1
                Σ[a,b] = Σ[b,a]
            end
        end
    end

    for a = 2:o
        for b = 1:a-1
            Σ[a,b] = Σ[b,a]
        end
    end

    MvNormLeaf(MvNormal(leaf_mean, Σ))
end

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

immutable InternalNode{T<:Any, LeafType<:Leaf}
    featid::Int # id of the feature we use to split
    featval::T # value over which we split
    left::Union{LeafType, InternalNode}
    right::Union{LeafType, InternalNode}
end

immutable Ensemble{LeafType<:Leaf}
    trees::Vector{Union{InternalNode, LeafType}}
end

type BuildTreeParameters{F<:LossFunction, LeafType<:Any}
    nsubfeatures::Int # number of candidate features to select for split; 0 will use all
    max_depth::Int # will only split if above this depth (depth of just a leaf is 0)
    min_samples_split::Int # minimum number of samples requried in a node to split
    min_samples_leaves::Int # minimum number of samples in newly created leaves
    min_split_improvement::Float64 # will only split if loss(split(leaf)) > loss(leaf) + min_split_improvement
    n_split_tries::Int # number of random splits to try (sampled uniformally from available predictors and data points)
    loss_function::Type{F}
    leaf_type::Type{LeafType} # the type used for the leaf
end
function build_tree_parameters{F<:LossFunction, LeafType<:Leaf}(;
    nsubfeatures::Int=0,
    max_depth::Int=typemax(Int),
    min_samples_split::Int=5,
    min_samples_leaves::Int=0,
    min_split_improvement::Float64=-Inf,
    loss_function::Type{F}=LossFunction_MSE,
    leaf_type::Type{LeafType}=MeanVecLeaf,
    )

    @assert(nsubfeatures ≥ 0)
    @assert(max_depth ≥ 0)
    @assert(min_samples_split ≥ 0)
    @assert(min_samples_leaves ≤ min_samples_split/2)

    BuildTreeParameters{F, LeafType}(nsubfeatures, max_depth, min_samples_split,
               min_samples_leaves, min_split_improvement, loss_function, leaf_type)
end

function Base.print(io::IO, params::BuildTreeParameters)
    println(io, "BuildTreeParameters:")
    @printf(io, "\tnsubfeatures:          %8d\n", params.nsubfeatures)
    @printf(io, "\tmax_depth:             %8d\n", params.max_depth)
    @printf(io, "\tmin_samples_split:     %8d\n", params.min_samples_split)
    @printf(io, "\tmin_samples_leaves:    %8d\n", params.min_samples_leaves)
    @printf(io, "\tmin_split_improvement: %8.3f\n", params.min_split_improvement)
    @printf(io, "\tloss_function:         %8s\n", string(params.loss_function))
end

Base.length(::Leaf) = 1
Base.length(tree::InternalNode) = length(tree.left) + length(tree.right)
Base.length(ensemble::Ensemble) = length(ensemble.trees)

depth(::Leaf) = 0
depth(tree::InternalNode) = 1 + max(depth(tree.left), depth(tree.right))

function print_tree(leaf::Leaf, depth::Integer=-1, indent::Integer=0)
    # NOTE(tim): I think this only works for classification
    print("    " ^ indent * "R-> ")
    println(leaf)
end
function print_tree(tree::InternalNode, depth::Integer=-1, indent::Integer=0)
    if depth == indent
        println()
        return
    end
    println("Feature $(tree.featid), Threshold $(tree.featval)")
    print("    " ^ indent * "L-> ")
    print_tree(tree.left, depth, indent + 1)
    print("    " ^ indent * "R-> ")
    print_tree(tree.right, depth, indent + 1)
end

apply_tree(leaf::Leaf, feature::Vector) = leaf
function apply_tree(tree::InternalNode, features::Vector)
    if features[tree.featid] < tree.featval
        apply_tree(tree.left, features)
    else
        apply_tree(tree.right, features)
    end
end

function apply_forest!{T<:Any}(μ::Vector{T}, forest::Ensemble{MeanVecLeaf}, features::Vector)

    # return the mean among tree's observations

    ntrees = length(forest)
    o = length(μ)

    for i in 1 : o
        μ[i] = 0.0
    end
    for i in 1 : ntrees
        leaf = apply_tree(forest.trees[i], features)::MeanVecLeaf
        for j in 1 : o
            μ[j] += leaf.μ[j]
        end
    end
    for i in 1 : o
        μ[i] /= ntrees
    end
    μ
end
function apply_forest!{T<:Any}(Σ::Matrix{T}, forest::Ensemble{CovLeaf}, features::Vector)

    # return the mean covariance matrix among the tree's covariances

    ntrees = length(forest)
    len = length(Σ)

    for i in 1 : len
        Σ[i] = 0.0
    end
    for i in 1 : ntrees
        leaf = apply_tree(forest.trees[i], features)::CovLeaf
        for j in 1 : len
            Σ[j] += leaf.Σ[j]
        end
    end
    for i in 1 : len
        Σ[i] /= ntrees
    end

    Σ
end

function _split{F<:LossFunction, T<:AbstractFloat, U<:Real}(
    data::TreeData{T,U},
    assignment_id::Int,

    min_samples_leaves::Int,
    min_split_improvement::Float64,
    loss_function::Type{F},
    n_thresholds::Int,
    )

    # returns a tuple: (index_of_feature_we_split_over, threshold)

    nrow, nfeatures = size(data.X)
    best = (0,0)
    best_loss = loss(loss_function, data, assignment_id) + min_split_improvement
    n_thresholds = min(nrow-1, n_thresholds)

    _reservoir_sample!(data.subfeature_indeces, nfeatures)

    for i in data.subfeature_indeces

        f_lo = Inf
        f_hi = -Inf
        for j in 1 : nrow
            if data.assignment[j] == assignment_id
                v = data.X[j,i]
                if v < f_lo
                    f_lo = v
                elseif v > f_hi
                    f_hi = v
                end
            end
        end

        for j in 1 : n_thresholds

            threshold = f_lo + (f_hi - f_lo)*(j/(n_thresholds+1))

            value, num_left, num_right = loss(loss_function, data, i, assignment_id, threshold)

            if value > best_loss && num_left ≥ min_samples_leaves && num_right ≥ min_samples_leaves
                best_loss = value
                best = (i, threshold)
            end
        end
    end

    best
end
function _split_unif_rand{F<:LossFunction, T<:AbstractFloat, U<:Real}(
    data::TreeData{T,U},
    assignment_id::Int,

    min_samples_leaves::Int,
    min_split_improvement::Float64,
    loss_function::Type{F},
    n_split_tries::Int,
    )

    # returns a tuple: (index_of_feature_we_split_over, threshold)

    id_count = get_id_count(data, assignment_id)
    nrows, nfeatures = size(data.X)
    id_range = 1:id_count
    predictor_range = 1:nfeatures

    @assert(id_count > 0)
    @assert(nfeatures > 0)

    best_id = 0
    best_threshold = 0.0
    best_loss = loss(loss_function, data, assignment_id) + min_split_improvement

    for _ in 1 : n_split_tries
        predictor_id = rand(predictor_range)
        nth_id = rand(id_range)
        row_index = get_nth_valid_id(data, assignment_id, nth_id)
        threshold = data.X[row_index, predictor_id]
        value, num_left, num_right = loss(loss_function, data, predictor_id, assignment_id, threshold)

        if value > best_loss && num_left ≥ min_samples_leaves && num_right ≥ min_samples_leaves
            best_loss = value
            best_id = predictor_id
            best_threshold = threshold
        end
    end

    (best_id, best_threshold)
end

function _build_tree{T<:AbstractFloat, U<:Real}(
    data::TreeData{T,U},
    assignment_id::Int,
    params::BuildTreeParameters,
    _depth::Int=0,
    )

    n = get_num_samples(data)
    n_for_id = 0
    for i = 1 : n
        n_for_id += convert(Int, data.assignment[i] == assignment_id)
    end

    # println("size(assignment): ", size(assignment))
    # println(assignment[1:20])
    # println("assignment_id:    ", assignment_id)
    # println("n_for_id: ", n_for_id)
    # println("_depth ≥ params.max_depth   ", _depth ≥ params.max_depth)
    # println("n_for_id < params.min_samples_split:   ", n_for_id < params.min_samples_split)

    if _depth ≥ params.max_depth || n_for_id < params.min_samples_split
        return build_leaf(params.leaf_type, data, assignment_id)
    end

    id, threshold = _split(data, assignment_id,
                        params.min_samples_leaves,
                        params.min_split_improvement,
                        params.loss_function,
                        params.n_split_tries,
                        )

    # println("split id:      ", id)
    # println("split thresh:  ", threshold)

    if id == 0
        return build_leaf(params.leaf_type, data, assignment_id)
    end

    next_id = data.assignment[end]+1
    data.assignment[end] = next_id

    for i in 1 : n
        if data.assignment[i] == assignment_id && data.X[i,id] < threshold
            data.assignment[i] = next_id
        end
    end

    return InternalNode{T, params.leaf_type}(id, threshold,
                _build_tree(data, next_id,       params, _depth+1),
                _build_tree(data, assignment_id, params, _depth+1))
end
function build_tree{T<:AbstractFloat, U<:Real}(labels::Matrix{T}, features::Matrix{U}, params::BuildTreeParameters)

    o, n = size(labels)
    μ_l = Array(T, o)
    μ_r = Array(T, o)
    Σ_l = Array(T, o,o)
    Σ_r = Array(T, o,o)

    nfeatures = size(features,2)
    nsubfeatures = params.nsubfeatures > 0 ? min(params.nsubfeatures, nfeatures) : nfeatures
    subfeature_indeces = Array(Int, nsubfeatures)
    assignment = ones(Int, n+1)

    data = TreeData{T,U}(labels, features, assignment, μ_l, μ_r, Σ_l, Σ_r, subfeature_indeces)

    _build_tree(data, 1, params)
end

function _build_forest_parallel{T<:AbstractFloat, U<:Real}(labels::Matrix{T}, features::Matrix{U}, ntrees::Integer, params::BuildTreeParameters, partialsampling::Float64)
    nlabels = _nsamples_labels(labels)
    nsamples = _int(partialsampling * nlabels)

    forest = @parallel (vcat) for i in 1:ntrees
        inds = rand(1:nlabels, nsamples)
        build_tree(labels[:,inds], features[inds,:], params)
    end
    Ensemble{params.leaf_type}(convert(Vector{Union{InternalNode, params.leaf_type}}, forest))
end
function _build_forest_serial{T<:AbstractFloat, U<:Real}(labels::Matrix{T}, features::Matrix{U}, ntrees::Integer, params::BuildTreeParameters, partialsampling::Float64)

    o, nlabels = size(labels)
    μ_l = Array(T, o)
    μ_r = Array(T, o)
    Σ_l = Array(T, o,o)
    Σ_r = Array(T, o,o)

    nfeatures = size(features,2)
    nsubfeatures = params.nsubfeatures > 0 ? min(params.nsubfeatures, nfeatures) : nfeatures
    subfeature_indeces = Array(Int, nsubfeatures)
    assignment = ones(Int, nlabels+1)
    nsamples = _int(partialsampling * nlabels)

    chosen_samples = Array(Int, nsamples)

    data = TreeData{T,U}(labels, features, assignment, μ_l, μ_r, Σ_l, Σ_r, subfeature_indeces)

    forest = Array(Union{InternalNode, params.leaf_type}, ntrees)
    for i = 1 : ntrees

        # select nsamples
        _reservoir_sample!(chosen_samples, nlabels)
        for j = 1:nlabels
            assignment[j] = 0 # remove labels that we are not considering
        end
        for j in chosen_samples
            assignment[j] = 1
        end
        assignment[end] = 1

        # build the tree
        forest[i] = _build_tree(data, 1, params)
    end

    Ensemble{params.leaf_type}(forest)
end
function build_forest{T<:AbstractFloat, U<:Real}(labels::Matrix{T}, features::Matrix{U}, ntrees::Integer, params::BuildTreeParameters, partialsampling::Float64=0.7)

    @assert(0.0 < partialsampling ≤ 1.0)
    if nprocs() == 1
        _build_forest_serial(labels, features, ntrees, params, partialsampling)
    else
        _build_forest_parallel(labels, features, ntrees, params, partialsampling)
    end
end

#####################################
# analysis functions

get_split_index_set(leaf::Leaf) = Set{Int}()
function get_split_index_set(tree::InternalNode)
    # returns a Set{Int} of indeces used as splits in the tree
    retval = Set{Int}()
    push!(retval, tree.featid)
    for featid in get_split_index_set(tree.left)
        push!(retval, featid)
    end
    for featid in get_split_index_set(tree.right)
        push!(retval, featid)
    end
    retval
end
function get_split_index_set{L<:Leaf}(forest::Ensemble{L})
    # returns a Set{Int} of indeces used as splits in the forest
    retval = Set{Int}()
    for tree in forest.trees
        for featid in get_split_index_set(tree)
            push!(retval, featid)
        end
    end
    retval
end

#####################################
# helper functions

function _calc_covariance!{T, U}(
    Σ::Matrix{T},
    data::TreeData{T,U},
    assignment_id::Int,
    )

    # println(fit_mle(MvNormal, data.Y))

    n, o, p = size(data)
    leaf_mean = zeros(T, o)

    n_for_id = 0
    for i = 1 : n
        if data.assignment[i] == assignment_id
            n_for_id += 1
            for j = 1 : o
                leaf_mean[j] += data.Y[j,i]
            end
        end
    end
    for j = 1 : o
        leaf_mean[j] /= n_for_id
    end

    # calc covariance

    Σ = zeros(T, o, o)

    for i = 1 : n
        if data.assignment[i] == assignment_id
            for a = 1 : o
                l = data.Y[a,i] - leaf_mean[a]
                for b = a : o
                    Σ[a,b] += l*(data.Y[b,i]-leaf_mean[b])
                end
            end
        end
    end

    for i = 1 : o*o
        Σ[i] /= (n_for_id-1) # we assume that we estimated the mean, so we use (n-1)
    end

    _diagonal_shrinkage!(Σ)

    # copy over symmetric part of matrix
    for a = 2:o
        for b = 1:a-1
            Σ[a,b] = Σ[b,a]
        end
    end

    Σ
end
function _diagonal_shrinkage!{T<:Real}(
    Σ::Matrix{T}, # 2×2 matrix with upper diagonal filled in
                  # this method ONLY affects the upper diagonal
    ε::T = convert(T, 1e-6),
    )

    # now we want to use Shrinkage Estimation - https://en.wikipedia.org/wiki/Estimation_of_covariance_matrices#Shrinkage_estimation
    # Since with a small number of samples we might get a Σ that isn't positive semidefinite
    # A matrix is pos. semidefinite it all of its eigenvalues are > 0
    # We choose γ such that A = γΣ + (1-γI) with all of eigenvalues are > ε
    # This method will work in general, but what is coded below assumes Σ is 2×2
    #      Σ = [a b; b d]
    @assert(size(Σ,1)==2 && size(Σ,2)==2)

    a, b, d = Σ[1,1], Σ[1,2], Σ[2,2]

    γ = 2.0*(1.0-ε) / (2.0 - a - d + sqrt(a*a -2*a*d + 4*b*b + d*d))

    if 0.0 < γ < 1.0
        # @assert(γ ≥ 0.0)

        Σ[1,1] = γ*a + (1.0-γ)
        Σ[1,2] = γ*b
        Σ[2,2] = γ*d + (1.0-γ)
    end

    Σ
end
function _fast_det{T<:AbstractFloat}(M::AbstractMatrix{T})
    n = size(M,1)
    if n == 1
        M[1]
    elseif n == 2
        M[1]*M[4] - M[2]*M[3]
    elseif n == 3
        M[1]*M[5]*M[9] + M[4]*M[8]*M[3] + M[7]*M[2]*M[6] -
            M[7]*M[5]*M[3] - M[4]*M[2]*M[9] - M[1]*M[8]*M[6]
    else
        det(M)
    end
end
function _reservoir_sample!(indeces::Vector{Int}, n::Int, k::Int=length(indeces))

    # performs Reservoir Sampling to obtain k elements from [1:n]
    #  indeces = input array
    #  n       = length of set of indeces to choose from (ie, we select from [1:n])
    #  k       = number of samples to choose and place in first k elements of indeces
    # runtime: O(n)
    # https://en.wikipedia.org/wiki/Reservoir_sampling

    if n == 0
        n = k
    end

    for i = 1 : k
        indeces[i] = i
    end

    for i = k+1 : n
        j = rand(1:i)
        if j ≤ k
            indeces[j] = i
        end
    end

    indeces
end

end # end module


text:
module Spid

using Base.Threads
using BioSequences
using CodecZlib
using DataFrames
using CSV
using FASTX

const valid_alleles = "ACGT";
const skip_chars = "N-";

# Only used by script entry-points, but include them here to
# pre-compile
using ArgParse

function read_fasta!(fasta, names, seqs)
    for record in FASTA.Reader(fasta)
        name = FASTA.identifier(record)
        seq = String(sequence(LongCharSeq, record))
        push!(names, name)
        push!(seqs, seq)
    end
end

function concat_fasta_contigs(fasta_path)
    open(fasta_path) do f
        char_seq = Char[]
        stream = f
        if endswith(fasta_path, ".gz")
            stream = GzipDecompressorStream(f)
        end
        for record in FASTA.Reader(stream)
            append!(char_seq, sequence(LongCharSeq, record))
        end
        return LongCharSeq(join(char_seq))
    end
end

function combine_sample_fastas(samples, in_fastas, out_fasta)
    # TODO check that fastas have same contig names and lengths
    open(out_fasta, "w") do f
        w = FASTA.Writer(f)
        for (sample_name, sample_fasta) in zip(samples, in_fastas)
            write(w, FASTA.Record(sample_name,
                                  concat_fasta_contigs(sample_fasta)))
        end
    end
end

"""
    merge_alignments(out_prefix, fastas, samples)

Merge several samples into a combined FASTA and compute distance matrix.

Outputs the following FASTA files:
- `{out_prefix}.fa`: all sites
- `{out_prefix}.core.fa`: core sites only (no missing alleles)
- `{out_prefix}.variantsOnly.fa`: variant sites only
- `{out_prefix}.core.variantsOnly.fa`: core variants only

Also outputs pairwise distances in the following CSV files:
- `{out_prefix}.fa.pairwise_diffs.csv`: all sites
- `{out_prefix}.core.fa.pairwise_diffs.csv`: core sites only

`fastas` should be a list of FASTA filenames aligned against the same
reference, e.g. via [`align_short_reads`](@ref) or
[`align_assembly`](@ref). They should have the same number of contigs,
with the same lengths, in the same order. The contig names are
ignored. The contigs for each sample are concatenated into a
single sequence in the output FASTA.

`samples` should be either a list of the sample names, or a `Regex`
for extracting a sample name from the FASTA filename; in particular,
the first capturing group of the regex is taken as the sample name.
"""
function merge_alignments(
    out_prefix, fastas, samples)

    samples_dict = Dict{String, String}()
    if samples == nothing
        samples = fastas
    elseif typeof(samples) == Regex
        samples = [match(samples, fa)[1] for fa in fastas]
    end

    merged_fa = "$out_prefix.fa"
    combine_sample_fastas(samples, fastas, merged_fa)
    CSV.write("$out_prefix.fa.pairwise_diffs.csv",
              fasta_edit_distances(merged_fa))

    variants_fa = "$out_prefix.variantsOnly.fa"
    open(merged_fa) do in_f
        open(variants_fa, "w") do out_f
            w = FASTA.Writer(out_f)
            for record in get_variant_only_seqs(in_f)
                write(w, record)
            end
        end
    end

    core_fa = "$out_prefix.core.fa"
    core_seqs = Vector{FASTA.Record}()
    open(merged_fa) do in_f
        core_seqs = get_core_only_seqs(in_f)
    end

    if (length(core_seqs) == 0 ||
        length(sequence(LongCharSeq, core_seqs[1])) == 0)
        @warn "Core genome appears to be empty; skipping core genome outputs."
        return
    end

    open(core_fa, "w") do out_f
        w = FASTA.Writer(out_f)
        for record in core_seqs
            write(w, record)
        end
    end
    CSV.write("$out_prefix.core.fa.pairwise_diffs.csv",
              fasta_edit_distances(core_fa))

    core_variants_fa = "$out_prefix.core.variantsOnly.fa"
    open(core_fa) do in_f
        open(core_variants_fa, "w") do out_f
            w = FASTA.Writer(out_f)
            for record in get_variant_only_seqs(in_f)
                write(w, record)
            end
        end
    end
end

export merge_alignments

include("FastaEditDistances.jl")
include("FilterFastaSites.jl")
include("Pileup2Consensus.jl")
include("Align.jl")

end # module


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using RoboDojo
using BenchmarkTools
using Folds

@show Threads.nthreads()

# ## Initial conditions
q1 = nominal_configuration(hopper)
v1 = zeros(hopper.nq)

# ## Time
T = 1
h = 0.01

# ## Simulator
s = Simulator(hopper, T, h=h)

# ## Simulate
@benchmark simulate!(s, q1, v1) setup=(s=s, q1=q1, v1=v1)

# ## Visualizer
vis = Visualizer()
render(vis)

# ## Visualize
visualize!(vis, s)

# ## for loop
function sim_forloop!(S, Q, V, N)
    for i = 1:N
        simulate!(S[i], Q[i], V[i])
    end
    return
end

function sim_threads!(S, Q, V, N)
    Threads.@threads for i = 1:N
        simulate!(S[i], Q[i], V[i])
    end
    return
end

N = 16
S = [deepcopy(s) for i = 1:N]
Q = [q1 for i = 1:N]
V = [v1 for i = 1:N]
X = [(S[i], Q[i], V[i]) for i = 1:N]

sim_forloop!(S, Q, V, N)
@benchmark sim_forloop!(S, Q, V, N) setup=(S=S, Q=Q, V=V, N=N)

sim_threads!(S, Q, V, N)
@benchmark sim_threads!(S, Q, V, N) setup=(S=S, Q=Q, V=V, N=N)
visualize!(vis, S[1])

Folds.map(simulate!, S, Q, V, Folds.ThreadedEx())
@benchmark Folds.map(simulate!, S, Q, V) setup=(simulate! =simulate!, S=S, Q=Q, V=V)
@benchmark Folds.map(simulate!, S, Q, V, Folds.ThreadedEx()) setup=(simulate! =simulate!, S=S, Q=Q, V=V)


text:
function get_FE_spaces(
  strategy::WeakForms.MeshStrategy,
  coupling::WeakForms.Coupling{:strong},
  models,
  order,
  bconds;
  constraint=nothing)

  # Reference FE
  reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
  reffeₚ = ReferenceFE(lagrangian,Float64,order-1)

  # Test FE Spaces
  Vu_FSI = TestFESpace(models[:Ω], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_tags])
  Vv_FSI = TestFESpace(models[:Ω], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_tags])
  Vu_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags])
  Vv_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vv_tags])
  Q = TestFESpace(models[:Ωf], reffeₚ, constraint=constraint, conformity=:C0)

  # Trial FE Spaces
  Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values])
  Uu_FSI = TransientTrialFESpace(Vu_FSI,bconds[:FSI_Vu_values])
  Uv_FSI = TransientTrialFESpace(Vv_FSI,bconds[:FSI_Vv_values])
  P = TrialFESpace(Q)

  # Multifield FE Spaces
  fe_spaces = (
    Y_ST = MultiFieldFESpace([Vu_ST,Vv_ST,Q]),
    X_ST = MultiFieldFESpace([Uu_ST,Uv_ST,P]),
    Y_FSI = MultiFieldFESpace([Vu_FSI,Vv_FSI,Q]),
    X_FSI = TransientMultiFieldFESpace([Uu_FSI,Uv_FSI,P])
  )
end

function get_FE_spaces(
  strategy::WeakForms.MeshStrategy{:biharmonic},
  coupling::WeakForms.Coupling{:strong},
  models,
  order,
  bconds;
  constraint=nothing)

  # Reference FE
  reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
  reffeₚ = ReferenceFE(lagrangian,Float64,order-1)

  # Test FE Spaces
  Vw_FSI = TestFESpace(models[:Ω], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_tags])
  Vu_FSI = TestFESpace(models[:Ω], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_tags])
  Vv_FSI = TestFESpace(models[:Ω], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_tags])
  Vw_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags])
  Vu_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags])
  Vv_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vv_tags])
  Q = TestFESpace(models[:Ωf], reffeₚ, constraint=constraint, conformity=:C0)

  # Trial FE Spaces
  Uw_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values])
  Uw_FSI = TrialFESpace(Vu_FSI,bconds[:ST_Vu_values])
  Uu_FSI = TransientTrialFESpace(Vu_FSI,bconds[:FSI_Vu_values])
  Uv_FSI = TransientTrialFESpace(Vv_FSI,bconds[:FSI_Vv_values])
  P = TrialFESpace(Q)

  # Multifield FE Spaces
  fe_spaces = (
    Y_ST = MultiFieldFESpace([Vw_ST,Vu_ST,Vv_ST,Q]),
    X_ST = MultiFieldFESpace([Uw_ST,Uu_ST,Uv_ST,P]),
    Y_FSI = MultiFieldFESpace([Vw_FSI,Vu_FSI,Vv_FSI,Q]),
    X_FSI = TransientMultiFieldFESpace([Uw_FSI,Uu_FSI,Uv_FSI,P])
  )
end

function get_FE_spaces(
  strategy::WeakForms.MeshStrategy,
  coupling::WeakForms.Coupling{:weak},
  models,
  order,
  bconds;
  constraint=nothing)

  # Reference FE
  reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
  reffeₚ = ReferenceFE(lagrangian,Float64,order-1)

  # Test FE Spaces
  Vu_FSI_f = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_f_tags])
  Vv_FSI_f = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_f_tags])
  Vu_FSI_s = TestFESpace(models[:Ωs], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_s_tags])
  Vv_FSI_s = TestFESpace(models[:Ωs], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_s_tags])
  Vu_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags])
  Vv_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vv_tags])
  Q = TestFESpace(models[:Ωf], reffeₚ, constraint=constraint, conformity=:C0)

  # Trial FE Spaces
  Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values])
  Uu_FSI_f = TransientTrialFESpace(Vu_FSI_f,bconds[:FSI_Vu_f_values])
  Uv_FSI_f = TransientTrialFESpace(Vv_FSI_f,bconds[:FSI_Vv_f_values])
  Uu_FSI_s = TransientTrialFESpace(Vu_FSI_s,bconds[:FSI_Vu_s_values])
  Uv_FSI_s = TransientTrialFESpace(Vv_FSI_s,bconds[:FSI_Vv_s_values])
  P = TrialFESpace(Q)

  # Multifield FE Spaces
  fe_spaces = (
    Y_ST = MultiFieldFESpace([Vu_ST,Vv_ST,Q]),
    X_ST = MultiFieldFESpace([Uu_ST,Uv_ST,P]),
    Y_FSI = MultiFieldFESpace([Vu_FSI_f,Vv_FSI_f,Q,Vu_FSI_s,Vv_FSI_s]),
    X_FSI = TransientMultiFieldFESpace([Uu_FSI_f,Uv_FSI_f,P,Uu_FSI_s,Uv_FSI_s])
  )
end

function get_FE_spaces(
  strategy::WeakForms.MeshStrategy{:biharmonic},
  coupling::WeakForms.Coupling{:weak},
  models,
  order,
  bconds;
  constraint=nothing)

  # Reference FE
  reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
  reffeₚ = ReferenceFE(lagrangian,Float64,order-1)

  # Test FE Spaces
  Vw_FSI_f = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vw_f_tags])
  Vu_FSI_f = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_f_tags])
  Vv_FSI_f = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_f_tags])
  Vu_FSI_s = TestFESpace(models[:Ωs], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vu_s_tags])
  Vv_FSI_s = TestFESpace(models[:Ωs], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:FSI_Vv_s_tags])
  Vw_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags])
  Vu_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vu_tags])
  Vv_ST = TestFESpace(models[:Ωf], reffeᵤ, conformity=:H1, dirichlet_tags=bconds[:ST_Vv_tags])
  Q = TestFESpace(models[:Ωf], reffeₚ, constraint=constraint, conformity=:C0)

  # Trial FE Spaces
  Uw_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uu_ST = TrialFESpace(Vu_ST,bconds[:ST_Vu_values])
  Uv_ST = TrialFESpace(Vv_ST,bconds[:ST_Vv_values])
  Uw_FSI_f = TrialFESpace(Vw_FSI_f,bconds[:FSI_Vw_f_values])
  Uu_FSI_f = TransientTrialFESpace(Vu_FSI_f,bconds[:FSI_Vu_f_values])
  Uv_FSI_f = TransientTrialFESpace(Vv_FSI_f,bconds[:FSI_Vv_f_values])
  Uu_FSI_s = TransientTrialFESpace(Vu_FSI_s,bconds[:FSI_Vu_s_values])
  Uv_FSI_s = TransientTrialFESpace(Vv_FSI_s,bconds[:FSI_Vv_s_values])
  P = TrialFESpace(Q)

  # Multifield FE Spaces
  fe_spaces = (
    Y_ST = MultiFieldFESpace([Vw_ST,Vu_ST,Vv_ST,Q]),
    X_ST = MultiFieldFESpace([Uw_ST,Uu_ST,Uv_ST,P]),
    Y_FSI = MultiFieldFESpace([Vw_FSI_f,Vu_FSI_f,Vv_FSI_f,Q,Vu_FSI_s,Vv_FSI_s]),
    X_FSI = TransientMultiFieldFESpace([Uw_FSI_f,Uu_FSI_f,Uv_FSI_f,P,Uu_FSI_s,Uv_FSI_s])
  )
end


text:
mutable struct RobbinsMonro
    κ::Float64
    τ::Float64
    state::IdDict
end

function RobbinsMonro(κ::Real=0.51, τ::Real=1)
    @assert 0.5 < κ <= 1 "κ should be in the interval (0.5,1]"
    @assert τ > 0 "τ should be positive"
    return RobbinsMonro(κ, τ, IdDict())
end

function Optimise.apply!(o::RobbinsMonro, x, Δ)
    κ = o.κ
    τ = o.τ
    n = get!(o.state, x, 1)
    Δ .*= 1 / (τ + n)^κ
    o.state[x] = n + 1
    return Δ
end

mutable struct ALRSVI
    τ::Int64
    state::IdDict
end

struct ALRSVIBase{T<:Real}
    g::Array{T}
    h::Base.RefValue{T}
    ρ::Base.RefValue{T}
    τ::Base.RefValue{T}
    t::Base.RefValue{Int64}
end

""" Construct Adaptive Learning Rate for Stochastic Variational Inference"""
function ALRSVI(ρ::Real=0.1, τ::Int=100)
    return ALRSVI(ρ, τ, IdDict())
end

function init!(model::AbstractGP{T,L,<:AnalyticVI}) where {T,L}
    for n_s in 1:(model.vi_opt[1].opt.τ)
        model.inference.MBIndices .= StatsBase.sample(
            1:(model.inference.nSamples), model.inference.nMinibatch; replace=false
        )
        model.inference.xview = view(model.X, model.inference.MBIndices, :)
        model.inference.yview = view_y(model.likelihood, model.y, model.inference.MBIndices)
        compute_kernel_matrices!(model)
        local_updates!(likelihood(model), yview(model), mean_f(model), var_f(model))
        natural_gradient!.(
            ∇E_μ(model.likelihood, model.inference.vi_opt[1], yview(model)),
            ∇E_Σ(model.likelihood, model.inference.vi_opt[1], yview(model)),
            model.inference,
            model.inference.vi_opt,
            Zviews(model),
            model.f,
        )
        init_ALRSVI!.(model.inference.vi_opt, model.f, τ)
    end
    return finalize_init_ALRSVI!.(model.inference.vi_opt, model.f)
end

function init_ALRSVI!(vi_opt::AVIOptimizer, gp::AbstractLatent{T}, τ) where {T}
    objη₁ = get!(
        vi_opt.opt.state,
        gp.η₁,
        ALRSVIBase(zeros(T, gp.dims), Ref(zero(T)), Ref(1.0), Ref(τ), Ref(1)),
    )
    objη₁.g .+= inf.∇η₁
    objη₁.h[] += norm(inf.∇η₁)

    objη₂ = get!(
        vi_opt.opt.state,
        gp.η₂,
        ALRSVIBase(zeros(T, gp.dims, gp.dims), Ref(zero(T)), Ref(1.0), Ref(τ), Ref(1)),
    )
    objη₂.g .+= inf.∇η₂
    return objη₂.h[] += sum(abs2, inf.∇η₂)
end

function finalize_init_ALRSVI!(vi_opt::AVIOptimizer, gp::AbstractLatent{T}) where {T}
    objη₁ = get(vi_opt.opt.state, gp.η₁)
    objη₁.g ./= obj.τ
    objη₁.h[] /= obj.τ
    objη₁.ρ[] = dot(objη₁.g, objη₁.g) / objη₁.h
    objη₂ = get(vi_opt.opt.state, gp.η₂)
    objη₂.g ./= obj.τ
    objη₂.h[] /= obj.τ
    return objη₂.ρ[] = sum(abs2, objη₂.g) / objη₂.h
end

function Optimise.apply!(opt::ALRSVI, x, Δ)
    # update timestep
    obj = get(opt.state, x)
    obj.t += 1
    obj.g .= (1.0 - 1.0 / obj.τ) * obj.g .+ 1.0 / obj.τ * Δ
    obj.h[] = (1.0 - 1.0 / obj.τ) * obj.h[] + 1.0 / obj.τ * sum(abs2, Δ)
    obj.ρ[] = dot(obj.g, obj.g) / obj.h[]
    obj.τ[] = obj.τ[] * (1 - obj.ρ[]) + 1.0
    return obj.ρ * Δ
end


text:
proclist = workers()         # Get list of workers

# test that can use tasks on all workers
for proc in proclist
  println("Julia says welcome from proc ",proc," w/ hostname ",@fetchfrom(proc,gethostname() ))
end
#Libc.flush_cstdio()

# Functions for demo calculation
function calc_pi(N::Integer)
  piapprox = @sync @parallel (+) for i in 1:N
    x, y = rand(), rand()
    (x*x+y*y <= 1.0) ? 1 : 0
  end
  piapprox *= 4.0/N
end

# Perform demo calculation
srand(42)        # make sure different processors have different seeds
calc_pi(10)
@time piapprox = calc_pi(1000000)
@printf("pi is estimated to be %14.12f\n", piapprox)


text:
"""
    getParameters(
        N::Int64,
        C::Float64,
        σ::Float64
        )

Get matrix of interactions `A` and instrisic growth rates `r`.

# Arguments:
- `N::Int64` : number of species
- `C::Float64` : connectivity
- `σ::Float64` : interaction strength
"""
function getParameters(
            N::Int64,
            C::Float64,
            σ::Float64
            )
    A = [ i==j ? -1. : rand()<C ? rand(Normal(0.,σ^2)) : 0. for i in 1:N, j in 1:N ]
    r = rand(N,1)
    return A, r
end

"""
    getGLV(
        x,
        pars::Array,
        t
    )

Get GLV dynamics.

# Arguments:
- `pars::Array` : concatenation of matrix `A` and vector `r`
- `x` and `t` : placeholders for ODE solver
"""
function getGLV(
            x,
            pars::Array,
            t
            )
    A = pars[:, 1:end-1]
    r = pars[:, end]
    inx = x .< 0
    x[inx] .= 0.0
    dx = diagm(vec(x)) * ( A * x + r )
    dx[inx] .= 0.0
    return dx
end

"""
    getIC(
        N::Int64,
        M::Int64
        )

Get collection of species `z` and numerical initial condition `x₀`.

# Arguments:
- `N::Int64` : number of species
- `M::Int64` : number of communities
"""
function getIC(
            N::Int64,
            M::Int64
            )
    z = zeros(Int64,N,M)
    inx = Any[]
    while size(inx,1) < M
        push!(inx, sort(unique( rand(1:N,rand(1:N)) )))
        unique!(inx)
    end
    for i in 1:M
        z[inx[i],i] .= 1
    end
    x₀ = z .* rand(N,M)
    return z, x₀
end

"""
    getSteadyState(
        x₀::Vector,
        τ::Float64,
        A::Array{Float64},
        r::Array{Float64}
        )

Get steady-state of dynamics by integrating ODE.

# Arguments:
- `x₀::Vector` : numerical initial condition
- `τ::Float64` : integration time
- `A::Array{Float64}` : interaction matrix
- `r::Array{Float64}` : intrinsic growth rates
"""
function getSteadyState(
            x₀::Vector,
            τ::Float64,
            A::Array{Float64},
            r::Array{Float64}
            )
    prob = ODEProblem(getGLV, x₀, (0.0,τ), hcat(A, r))
    sol = solve(prob,force_dtmin=true)
    xf = sol.u[end]
    return xf
end

"""
    getRewiredNetwork(
        A::Array{Float64},
        p::Float64
        )

Get rewired interaction matrix.

# Arguments:
- `A::Array{Float64}` : interaction matrix
- `p::Float64` : rewiring probability
"""
function getRewiredNetwork(
            A::Array{Float64},
            p::Float64
            )
    N = size(A,1)
    sA = sparse(A + diagm(ones(N)))
    row, col, s = findnz(sA)
    if p > 0.0
        for i in findall(rand(nnz(sA)) .< p)
            alreadyconnected = row[col .== col[i]]
            possibletargets = setdiff(1:N, alreadyconnected)
            length(possibletargets) >= 1 && (row[i] = rand(possibletargets))
        end
    end
    return sparse(row,col,s,N,N) - diagm(ones(N))
end

"""
        generate_data(
            N::Int64,
            M::Int64,
            repetitions::Int64,
            values::Int64,
            params::Array
        )

Generate synthetic data for parameters in `params`.
The function generates a number of `repetitions` for every value of `values`.

# Arguments:
- `N::Int64` : number of species
- `M::Int64` : number of communitites
- `repetitions::Int64` : number of ecological networks
- `values::Int64` : number of values per parameter
- `params::Array` : sweeping parameters
"""
function generate_data(N,M,repetitions,values,params)
    reps = 0:repetitions-1
    num = 1:values
    for par in params
        "Generating data: $par" |> println
        inx = [ (r,p) for r in reps, p in num ]
        @sync @distributed for i in inx
            # Changing parameters
            rep,p = i
            print(p-1,"_",rep)
            if par == "Strength"
                σ,c,η,ϵ,ρ = S[p],C[1],0,0,0.
            elseif par == "Connectivity"
                σ,c,η,ϵ,ρ = S[1],C[p],0,0,0.
            elseif par == "Universality"
                σ,c,η,ϵ,ρ = S[1],C[1],U[p],0,0.
            elseif par == "Noise"
                σ,c,η,ϵ,ρ = S[1],C[1],0,E[p],0.
            elseif par == "Rewiring"
                σ,c,η,ϵ,ρ = S[1],C[1],0,0,R[p]
            end
            print("  $σ,  $c, $η,  $ϵ,  $ρ\t")
            # Allocate
            X = SharedArray{Float64,2}(N,M)
            A = SharedArray{Float64,2}(N,N)
            r = SharedArray{Float64,2}(N,1)
            # Get initial values
            A₀,r = getParameters(N,c,σ)
            Z,x₀ = getIC(N,M)
            conv = 0
            j = 1
            while j <= M
                # Rewiring
                A = getRewiredNetwork(A₀,ρ)
                # Non-Universality: noise in non-zero entries
                if η > 0
                    A += η*( rand(Normal(0,σ^2),N,N) .* (x-> x != 0 ? 1 : 0 ).(A) )
                end
                # Solve
                x = getSteadyState(x₀[:,j], 1000., A, r)
                # Convergence
                conv += sum(diagm(vec(x)) * ( A * x + r))
                # Noise in non-zero entries
                x += ϵ*( rand(Uniform(0,mean(x[x.>0])),N) .* rand([-1,1],N) .* Z[:,j])
                x[x.<0].=0
                X[:,j] = x
                j += 1
                j%100==0&&print("··")
            end
            print("\t== $(conv/M) \n")

            path = "./data/synthetic/$N/$par/"
            !ispath(path) && mkpath(path)

            writedlm(path*"A$(p-1)_$rep.csv", hcat(A₀,r), ',')
            writedlm(path*"Z$(p-1)_$rep.csv", Z, ',')
            writedlm(path*"P$(p-1)_$rep.csv", X, ',')
        end
    end
end


text: @testset "Agent" begin end