text:
mutable struct NCS_growth_params{T<:Real, I<:Int64, C<:String}
    #Define paramter for the Neoclassical stochastic (Static labour) growth model
    #prefernces
    β :: T
    σ :: T
    #production function
    δ :: T
    α :: T
    #log of TFP process
    rho :: T
    sigma_e :: T
    #grid_point definitions
    lb :: T #lower bound for the capital
    ub :: T
    #method of comuptation, for now, only supproted "Smolyak"
    Meth :: C
    # simulations number
    n_sim :: I
    t_sim :: I
    burnout :: I
end
NCS_growth_params() = NCS_growth_params(0.984, 2.0, 0.025, 0.33, 0.979, 0.0072, 0.001, 50.0, "Smolyak", 1000,1000, 200)

function Smolyak_k_prim(K_guess,lb, ub, k_bord, β,σ, α, δ, grid,ϵ_vals,ρ, P,multi_ind,domain, inv_matrix)
"""
Solve future capital, given the guess for the future capital policy function. Then plug it to fixed point equation. Use the Smolyak grid
INPUTS:
K_guess : guessed policy function
lb : minimal level of capital in the economy
ub: maximal level of capital in the economy
k_bord:
β,σ, α, δ: parameters of the model
grid : Smolyak grid, multi_ind: index of the gird
inv_matrix: inversion of the tnesor product of the polynomial evaluations in the gird points

OUTPUT:
K_prim: New guess for the policy function
"""

    @assert β < 1.0
    @assert β > 0.0
    @assert σ >= 1.0
    @assert α < 1.0
    @assert α > 0.0
    @assert δ <= 1.0
    @assert δ >= 0.0
    @assert ρ < 1.0
    @assert ρ >= 0.0
    @assert ub > lb

    log_lb = log(lb)
    log_ub = log(ub)
    n = size(grid)[1]
    X = smolyak_weights(K_guess, inv_matrix)
    A =  length(ϵ_vals)
    K_prim = zeros(n)
    for i in 1:n

         k_prim_i = log(max(smolyak_evaluate(X,[grid[i,1],grid[i,2]],multi_ind,domain),lb))
         y = exp(grid[i,2])*exp(grid[i,1])^(α)
         c_i = max(y+(1-δ)*exp(grid[i,1]) -exp(k_prim_i),1e-8)
         EC =0.0 #expected future marginal utility
         for kk in 1:A
             a = ρ*grid[i,2] + ϵ_vals[kk]
             r_prim = α*exp(a)*exp(k_prim_i)^(α-1)-δ
             k_prim_ii = log(max(smolyak_evaluate(X,[k_prim_i,a],multi_ind,domain),lb))
             y_prim = exp(a)*exp(k_prim_i)^(α)
             c_prim = max(y_prim+(1-δ)*exp(k_prim_i) -exp(k_prim_ii),1e-9)
             EC = copy(EC)+ (c_prim/c_i)^(-σ)*P[kk]*(1+r_prim)*exp(k_prim_i)
         end
         K_prim[i] = β*EC # the fixed point equation from FOC
     end
     return(K_prim)
end

function Solve_Smolyak(Params::NCS_growth_params)
    """
    Solve Smolyak model for given parameters. For future information for the algorithm, check: Judd, Maliar, Maliar, Valero (2014)

    """

    #Read the parameters of the utility function:
    β = Params.β
    σ = Params.σ
    δ = Params.δ
    α = Params.α
    rho = Params.rho
    sigma_e = Params.sigma_e
    lb = Params.lb
    ub = Params.ub

    @assert β < 1.0
    @assert β > 0.0
    @assert σ >= 1.0
    @assert α < 1.0
    @assert α > 0.0
    @assert δ <= 1.0
    @assert δ >= 0.0
    @assert rho < 1.0
    @assert rho >= 0.0
    @assert ub > lb

    #Now compute the discretization of the probability of the shocks next period (using Gauss quadrature)
    sigma_un = sigma_e/((1.0-rho^2)^(0.5))

    function norm_pdf(x;std = sigma_e) #helpful function, it's for sure more elegant way to o this
        return pdf(Normal(0,std),x)
    end

    #compute Gauss quadrature
    ϵ_vals, P = gauss(norm_pdf, 12, -3.0*sigma_e, 3.0*sigma_e, rtol=1e-9)
    P = P./sum(P) #(normalize quadrature)

    a_bord = [-4.0*sigma_un,3.0*sigma_un] #choose maximial and minimal values for the log of TFP

    #grid for log of capital
    log_lb = log(lb)
    log_ub = log(ub)

    #Now compute the Smolyak grid and interpolants interpolnats matrix:
    d =2
    mu = 5
    domain = [log_lb a_bord[1]; log_ub a_bord[2]]
    grid, multi_ind = smolyak_grid(chebyshev_gauss_lobatto,d,mu,domain)
    inv_matrix = smolyak_inverse_interpolation_matrix_1(grid,multi_ind,domain)

    #define initial guess
    y = exp.(grid[:,2]).*exp.(grid[:,1]).^α

    #define function to find fixed point
    Smolyak_projection_1(x) = Smolyak_k_prim(x,lb, ub, a_bord, β,σ, α, δ, grid,ϵ_vals,rho, P,multi_ind,domain, inv_matrix)

    #solve the model using Andersson method
    solution = fixedpoint(Smolyak_projection_1, y, show_trace= true, xtol = 1e-8, ftol =1e-9, beta=0.05, iterations=1000)

    #stop if model not converge
    println("solution found")
    #define evaluation method
    weights = smolyak_weights(solution.zero,inv_matrix)
    smolyak_ev(point1, point2) = smolyak_evaluate(weights,[point1,point2],multi_ind,domain)

    return smolyak_ev, grid
end

# routines for the simulation of the model
#simulate TFP process
function mc_sample_path(P, sample_T, sample_size)

    # create vector of discrete RVs for each row
    dists = Categorical(P)

    # setup the simulation
    X = fill(0, sample_size) # allocate memory, or zeros(Int64, sample_size)
    SIM = zeros(sample_size,sample_T)
    for i in 1:sample_size
        X = rand.(dists, sample_T)
        SIM[i,:] = X
    end
    SIM_int = convert(Array{Int64},SIM)
    return SIM_int
end

function simulate_smolyak_NGS( Params::NCS_growth_params, Ev)
    """
    SImulate the model,for the given smolyak Evaluation function

    """
    println("simulation")
    #unpack params
    δ = Params.δ
    α = Params.α
    ρ = Params.rho
    lb = Params.lb
    ub = Params.ub
    n_sim = Params.n_sim
    t_sim = Params.t_sim
    burnout = Params.burnout
    sigma_e = Params.sigma_e

    @assert sigma_e >= 0.0
    @assert α < 1.0
    @assert α > 0.0
    @assert δ <= 1.0
    @assert δ >= 0.0
    @assert ρ < 1.0
    @assert ρ >= 0.0
    @assert ub > lb
    @assert n_sim >= 100
    @assert t_sim >= 100

    function norm_pdf(x;std = sigma_e) #helpful function, it's for sure more elegant way to o this
        return pdf(Normal(0,std),x)
    end

    #compute Gauss quadrature
    ϵ_vals, P = gauss(norm_pdf, 12, -3.0*sigma_e, 3.0*sigma_e, rtol=1e-9)
    P = P./sum(P) #(normalize quadrature)
    #allocate memory
    A = ones(n_sim, t_sim)
    W = zeros(n_sim, t_sim)
    C = zeros(n_sim, t_sim)
    Y = zeros(n_sim, t_sim)
    I = zeros(n_sim, t_sim)
    R = zeros(n_sim, t_sim)

    #define starting values

    w_0 = 0.0
    W_0 = w_0*zeros(n_sim)
    A_0 = ones(n_sim)
    W[:,1] = W_0
    A[:,1] = A_0

    #simulate wage process
    SIM = mc_sample_path(P, t_sim, n_sim)
    for t in 2:t_sim
        for i in 1:n_sim
            W[i,t] = ρ*W[i,t-1]+ϵ_vals[SIM[i,t]]
        end
    end

    #compute stats for wage process
    ρ_emp_sum = 0.0
    std_emp_sum = 0

    for i in 1: n_sim
        ρ_emp_sum = ρ_emp_sum+ cor(vec(W[i,1:t_sim-1]), vec(W[i,2:t_sim]))
        std_emp_sum = std_emp_sum+ std(W[i,2:t_sim].- W[i,1:t_sim-1])

    end
    ρ_emp = 1/n_sim*ρ_emp_sum
    std_emp = 1/n_sim*std_emp_sum

    #solve the model for given simulated values
    for t in 2:t_sim
        for i in 1:n_sim
            log_k = log(max(A[i,t-1],lb ))

            A[i,t] = max(Ev(log_k,W[i,t]),lb)
            Y[i,t] = exp(W[i,t])*A[i,t-1]^(α)
            C[i,t] = Y[i,t]+ (1 - δ)*A[i,t-1]-A[i,t]
            I[i,t] = -(1 - δ)*A[i,t-1]+A[i,t]
            R[i,t] = 1.0+α*exp(W[i,t])*A[i,t-1]^(α-1) - δ
        end
    end

    #compute basic stats from the simulation
    mean_R = mean(R[:,burnout:t_sim])

    Std_Y = std(Y[:,burnout:t_sim])
    Std_C = std(C[:,burnout:t_sim])
    Std_I = std(I[:,burnout:t_sim])
    Cov_CY = cor(vec(Y[:,burnout:t_sim]), vec(C[:,burnout:t_sim]))
    Sim_results = vcat(ρ_emp, std_emp, mean_R, Std_Y, Std_C, Std_I,Cov_CY)
    return Sim_results
end

#save for the solution of the Smolyak model
function Save_smolyak(Params::NCS_growth_params, Ev, Sim_results, grid)
    #make some plots and print solution

    #print the results:
    println("Model was solved and simulated using the Smolyak method")
    println("Simualtion results:")
    println("empirical ρ of TFP: ", Sim_results[1], " vs, in the model: ", Params.rho)
    println("empirical std: ", Sim_results[2], " vs, in the model: ", Params.sigma_e)
    println("mean interest rate: " , Sim_results[3])
    println("Std of output: " , Sim_results[4])
    println("Std of consumption: " , Sim_results[5])
    println("Std of investments: " , Sim_results[6])
    println("cor Y and C: " , Sim_results[7])

    #plot grid
    F = Ev.(grid[:,1], grid[:,2])
    scatter(grid[:,1], grid[:,2])
    savefig("capital_policy")

    scatter(exp.(grid[:,1]), F)
    savefig("grid_sol")
    #plot 3d consumption policy function

end

function Find_eq_NCS_Smolyak(Params::NCS_growth_params)

    if(Params.Meth == "Smolyak")
        SolveRAEq(Params, Solve_Smolyak, simulate_smolyak_NGS, Save_smolyak)
    else
        println("method not supported,.. yet ")
    end
end


text:
"""
Re-factoring and re-implementation of MGT
"""
module GapLas

import Printf: @printf
import GraphIdx
import GraphIdx: Graph, EdgeGraph, Edge, Weights, WeightedGraph
import GraphIdx: PrimMstMem, prim_mst_edges
import GraphIdx.Tree: RootedTree

import ..Dual: dual!, gap_vec!
import ..TreeDP: TreeDPMem, tree_dp!
import ..TreeFlow: extract_non_tree!, update_tree!

struct Sol{N}
    x::Array{Float64, N}
    α::Vector{Float64}
end

struct GapMem{N, WL<:Weights{Float64}}
    sol::Sol{N}
    y_tree::Array{Float64,N}
    α_tree::Vector{Float64}
    λ_tree::WL
    gamma::Vector{Float64}
    wgraph::WeightedGraph
    egraph::EdgeGraph
    dp_mem::TreeDPMem
    mst::PrimMstMem
    tree::RootedTree
end

function GapMem(y::Array{Float64,N}, graph::Graph, lambda::Weights{Float64}) where {N}
    @assert length(y) == GraphIdx.num_nodes(graph)
    root_node = 1
    m = GraphIdx.num_edges(graph)
    n = GraphIdx.num_nodes(graph)
    λ_tree = similar(lambda, n)
    α_tree = Vector{Float64}(undef, n)
    egraph = collect(graph)
    mst = PrimMstMem(egraph)
    lam = Float64[lambda[i] for i=1:m]
    tree = RootedTree(root_node, mst.parent)
    return GapMem(
        Sol(copy(y), zeros(Float64, m)), # sol
        similar(y),                     # y_tree
        α_tree,
        λ_tree,
        Vector{Float64}(undef, m),      # gamma
        WeightedGraph(mst.neighbors, lam),
        egraph,
        TreeDPMem(n),
        mst,
        tree,
    )
end

"""
General framework for iterated tree optimizer:
1. Allocate memory (of type `Mem`)
2. Call `gaplas!(mem, y, graph, λ, μ)` in every iteration (`max_iter` many times).

!!! note

    `mem` must contain a `gamma` and an `x` field

"""
function gaplas(
    y::Array{Float64,N},
    graph::GraphT,
    lambda::Weights{Float64},
    mu::W2,
    ::Type{Mem} = GapMem;
    max_iter::Int = 5,
    verbose::Bool = true,
)::Array{Float64,N} where {N, GraphT<:Graph, W2<:Weights{Float64}, Mem}
    mem = Mem(y, graph, lambda)
    for it in 1:max_iter
        gaplas!(mem, y, graph, lambda, mu)
        if verbose
            let gap = -sum(mem.gamma)
                @printf("%4d %12.4f\n", it, gap)
            end
        end
    end
    return mem.sol.x
end

"""

Optional: node weights `μ` are set to `GraphIdx.Ones{Float64}`
"""
function gaplas(
    y::Array{Float64,N},
    graph::Graph,
    lambda::Weights{Float64},
    ::Type{Mem} = GapMem;
    max_iter::Int = 5,
    verbose::Bool = true,
)::Array{Float64,N} where {N, Mem}
    gaplas(
        y, graph, lambda, GraphIdx.Ones{Float64}(), Mem; max_iter=max_iter, verbose=verbose
    )
end

"""
Called in every iteration.
Implements the actual optimization, i.e.
"""
function gaplas!(
    mem::GapMem{N, W1},
    y::Array{Float64,N},
    graph::Graph,
    lambda::Weights{Float64},
    mu::Weights{Float64},
) where {N, W1<:Weights{Float64}}
    find_gap_tree!(mem, y, graph, lambda)
    tree_dp!(mem.sol.x, mem.y_tree, mem.tree, mem.λ_tree, mu, mem.dp_mem)
    mem.α_tree .= vec(mem.sol.x) .- vec(mem.y_tree)
    dual!(mem.α_tree, mem.dp_mem.proc_order, mem.mst.parent)
    update_tree!(
        mem.sol.α,
        mem.α_tree,
        mem.mst.selected,
        mem.egraph,
        mem.mst.parent,
    )
end

"""
Determine a tree by choosing edges with high gap value.
"""
function find_gap_tree!(
    mem::GapMem{N, W1}, y::Array{Float64,N}, graph::Graph, lambda::Weights{Float64}
) where {N, W1<:Weights{Float64}}
    gap_vec!(mem.gamma, mem.sol.x, mem.sol.α, mem.wgraph, -1.0)
    prim_mst_edges(mem.gamma, mem.tree.root, mem.mst)
    mem.y_tree .= y
    extract_non_tree!(
        graph,
        mem.mst.parent,
        mem.y_tree,
        mem.sol.α,
        mem.λ_tree,
        lambda,
    )
end

end

ext: jl

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
include("typerecipes.jl")

@recipe(FieldPlot) do scene
    Attributes(
        symmetric_cbar = false;
        :arrowsize => automatic,
        :arrowcolor => automatic,
        :color => automatic,
        :colormap => :viridis,
        :colorrange => automatic,
        :linecolor => automatic,
        :lengthscale => 1.0f0,
        :linewidth => automatic,
        :lowclip => nothing,
        :highclip => nothing,
        :inspectable => theme(scene, :inspectable),
        :levels => 6,
    )
end

@recipe(ScatterVariable) do scene
    Attributes(;
        lengthscale_factor = 1,
        # linewidth_factor = 1,
        :arrowsize => 1,
        :linewidth => 1,
        :color => automatic,
        :markersize => 8,
        :strokewidth => 1.0,
        :strokecolor => :black,
        :colormap => :viridis,
        :colorrange => automatic,
        :lowclip => nothing,
        :highclip => nothing,
        :inspectable => theme(scene, :inspectable),
    )
end

function MakieCore.plot!(sc::ScatterVariable{<:Tuple{ScalarVariable{T}}}) where {T}
    v = sc[1]
    grid, data = unwrap(v)

    scattergrid = Observable(empty(grid[]))
    scattercolor = Observable(Float32[])
    cl = @lift ustrip(max(abs.(extrema($v))...))
    valuerange = @lift (-$cl, $cl)

    function update_plot(grid, color)
        empty!(scattergrid[])
        empty!(scattercolor[])

        append!(scattergrid[], grid)
        append!(scattercolor[], color)

        scattergrid[] = scattergrid[]
        scattercolor[] = scattercolor[]
    end

    replace_automatic!(sc, :colorrange) do
        valuerange
    end

    onany(update_plot, grid, data)

    update_plot(grid[], data[])

    scatter!(sc, scattergrid;
        sc.markersize,
        sc.strokewidth,
        sc.strokecolor,
        color=scattercolor,
        sc.colormap,
        sc.colorrange,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end

function MakieCore.plot!(sc::ScatterVariable{<:Tuple{VectorVariable{T}}}) where {T}
    v = sc[1]

    if hasunit(v)
        _v = @lift ustrip($v)
    else
        _v = v
    end

    arrow_norm = @lift Float32.(vec(norm.($_v)))
    maxarrow = @lift maximum(norm.($_v))

    # arrowsize = @lift $(sc.arrowsize_factor)*$arrow_norm
    lengthscale = @lift 1/$maxarrow
    # linewidth = @lift $(sc.linewidth_factor)*$arrow_norm
    valuerange = @lift extrema($arrow_norm)
    replace_automatic!(sc, :colorrange) do
        valuerange
    end

    replace_automatic!(sc, :colorrange) do
        valuerange
    end

    arrows!(sc, v;
        arrowcolor=arrow_norm,
        sc.arrowsize,
        linecolor=arrow_norm,
        lengthscale,
        sc.linewidth,
        sc.color,
        sc.colormap,
        sc.colorrange,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end

function symmetric_colorrange!(sc, f)
    onany(sc[:symmetric_cbar], f) do symmetric, f
        extr = map(vals->ustrip.(vals), extrema(f))
        if symmetric
            cl = map(e->max(abs.(e)...), extr)
            valuerange = map(c->(-c, c), cl)
        else
            valuerange = extr
        end
        replace_automatic!(sc, :colorrange) do
            valuerange
        end
    end
end

function MakieCore.plot!(sc::FieldPlot{<:Tuple{ScalarField{1}}})
    f = sc[1]

    grid, data = unwrap(f)

    symmetric_colorrange!(sc, f)
    notify(sc[:symmetric_cbar])

    lines!(sc, grid..., data;
        sc.colorrange,
        sc.colormap,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end

function MakieCore.plot!(sc::FieldPlot{<:Tuple{ScalarField{2}}})
    f = sc[1]

    grid, data = unwrap(f)

    symmetric_colorrange!(sc, f)
    notify(sc[:symmetric_cbar])

    heatmap!(sc, grid..., data;
        sc.colorrange,
        sc.colormap,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end

function MakieCore.plot!(sc::FieldPlot{<:Tuple{ScalarField{3}}})
    f = sc[1]

    symmetric_colorrange!(sc, f)
    notify(sc[:symmetric_cbar])
    @lift @debug "Contour plot for 3D ScalarField with " * string($(sc.levels)) * " levels"

    contour!(sc, f;
        sc.colorrange,
        sc.colormap,
        sc.color,
        sc.levels,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end

function MakieCore.plot!(sc::FieldPlot{<:Tuple{VectorField{3}}})
    f = sc[1]

    if hasunit(f)
        _f = @lift ustrip($f)
    else
        _f = f
    end

    f_norm = @lift vec(norm.($_f))
    maxarrow = @lift maximum(norm.($_f))
    normed_f = @lift $_f ./ $maxarrow

    # linewidth = @lift $(sc.linewidth_factor)*$arrow_norm
    valuerange = @lift extrema($f_norm)
    replace_automatic!(sc, :colorrange) do
        valuerange
    end

    arrows!(sc, normed_f;
        # arrowscale = @lift(norm($_f)/$maxarrow),
        # lengthscale = @lift(norm($_f)/$maxarrow),
        # sc.linewidth,
        # markerspace = SceneSpace,
        color = f_norm,
        sc.colormap,
        sc.colorrange,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end

function MakieCore.plot!(sc::FieldPlot{<:Tuple{VectorField{2}}})
    f = sc[1]

    if hasunit(f)
        _f = @lift ustrip($f)
    else
        _f = f
    end

    f_norm = @lift vec(norm.($_f))
    replace_automatic!(sc, :linecolor) do
        f_norm
    end
    replace_automatic!(sc, :arrowcolor) do
        sc.linecolor
    end
    maxarrow = @lift maximum(norm.($_f))
    # normed_f = @lift $_f ./ $maxarrow

    # linewidth = @lift $(sc.linewidth_factor)*$arrow_norm
    valuerange = @lift extrema($f_norm)
    replace_automatic!(sc, :colorrange) do
        valuerange
    end

    arrows!(sc, _f;
        arrowsize = lift(v->norm(v)/maxarrow[]*1.4, _f),
        sc.arrowcolor,
        lengthscale = 3,
        linewidth = lift(v->norm(v)/maxarrow[]*0.7, _f),
        sc.linecolor,
        sc.colormap,
        sc.colorrange,
        sc.highclip,
        sc.lowclip,
        sc.inspectable
    )

    return sc
end


text:
using Underscores
using Utils

const Num = Int8
const Deck = Array{Num,1}

struct Problem
    player::Array{Deck}
end

function parse_file(f)
    a = Deck()
    b = Deck()
    cur = a

    for line in eachline(f)
        if isempty(line)
            cur = b
        elseif isdigit(line[1])
            push!(cur, parse(Num, line))
        end
    end

    Problem([a, b])
end

function load()
    open(joinpath(@__DIR__, "22.txt"), "r") do f
        parse_file(f)
    end
end

function score(pone, ptwo)
    winner = isempty(pone) ? ptwo : pone
    reverse!(winner)
    @_ mapreduce(_[1] * _[2], +, enumerate(winner))
end

score((pone, ptwo)) = score(pone, ptwo)

function ex1(problem::Problem)
    pone = copy(problem.player[1])
    ptwo = copy(problem.player[2])
    while !isempty(pone) && !isempty(ptwo)
        one = popfirst!(pone)
        two = popfirst!(ptwo)
        if one < two
            push!(ptwo, two)
            push!(ptwo, one)
        else
            push!(pone, one)
            push!(pone, two)
        end
    end
    score(pone, ptwo)
end

const Cache = Set

combine(pone::Deck, ptwo::Deck) = hash((pone, ptwo))

function ex2(cache::Cache, pone::Deck, ptwo::Deck, game=1, round=1)
    isempty(pone) || isempty(ptwo) && return (pone, ptwo)

    while !isempty(pone) && !isempty(ptwo)
        marker = combine(pone, ptwo)
        if marker in cache
            throw(:one)
        else
            push!(cache, marker)
        end

        one = popfirst!(pone)
        two = popfirst!(ptwo)

        if length(pone) >= one && length(ptwo) >= two
            try
                sub = ex2(Cache(), pone[1:one], ptwo[1:two], game + 1)
                if isempty(sub[1])
                    winner = ptwo
                    card = two
                else
                    winner = pone
                    card = one
                end
            catch
                winner = pone
                card = one
            end
        else
            if one < two
                winner = ptwo
                card = two
            else
                winner = pone
                card = one
            end
        end

        push!(winner, card)
        push!(winner, one ⊻ two ⊻ card)

        round += 1
    end
    pone, ptwo
end

function ex2(problem::Problem)
    score(ex2(Cache(), copy(problem.player[1]), copy(problem.player[2])))
end

problem = load()

println("Winning player's score: ", ex1(problem))
println("Recursive winner's score: ", ex2(problem))

using Test

@testset "Crab Combat" begin
    @testset "example 1" begin
        input1 = """
        Player 1:
        9
        2
        6
        3
        1

        Player 2:
        5
        8
        4
        7
        10
        """

        example1 = parse_file(IOBuffer(input1))

        @test ex1(example1) == 306
        @test ex2(example1) == 291
    end

    @testset "example 2" begin
        input2 = """
        Player 1:
        43
        19

        Player 2:
        2
        29
        14
        """

        example2 = parse_file(IOBuffer(input2))

        @test_throws :one ex2(example2)
    end

    @testset "results" begin
        @test ex1(problem) == 32677
        @test ex2(problem) == 33661
    end
end

max_stars_repo_name: CmdQ/AoC


text:
### Schur Eigenspace solvers

struct BaseSolver end
struct ArnoldiSolver end
struct KrylovSolver end

import ArnoldiMethod
function schurvectors(T, pi, n, israte, ::ArnoldiSolver)
	D = Diagonal(sqrt.(pi))
	T̃ = D * T * D^-1
	s, hist = ArnoldiMethod.partialschur(T̃; nev=n, which=israte ? ArnoldiMethod.LR() : ArnoldiMethod.LM())
	X̃ = collect(s.Q)
	X = D^-1 * X̃
	X ./= X[1,1]
	X, s.eigenvalues
end

import KrylovKit
function schurvectors(T, pi, n, israte, ::KrylovSolver)
	D = Diagonal(sqrt.(pi))
	T̃ = D * T * D^-1
	R, Q, v, info = KrylovKit.schursolve(T̃, pi, n, israte ? :LR : :LM, KrylovKit.Arnoldi())
	X̃ = reduce(hcat, Q)
	X = D^-1 * X̃
	X = X ./ X[1,1]
	X, v
end

function schurvectors(T, pi, n, israte, ::BaseSolver)
    Tw = Diagonal(sqrt.(pi))*T*Diagonal(1 ./ sqrt.(pi)) # rescale to keep markov property
    Sw = schur!(Tw)                       # returns orthonormal vecs by def
    Xw, λ = selclusters!(Sw, n, israte)
    X  = Diagonal(1 ./sqrt.(pi)) * Xw              # scale back
    X  = X[1,1]>0 ? X : -X
    X, λ
end

# select the schurvectors corresponding to the n abs-largest eigenvalues
# if reverse==true select highest abs value, otherwise select lowest (for rate matrices)
function selclusters!(S, n, ratematrix)
    ind = sortperm(abs.(S.values), rev=!ratematrix) # get indices for dominant eigenvalues
    select = zeros(Bool, size(ind))           # create selection vector
    select[ind[1:n]] .= true
    S = ordschur!(S, select)                  # reorder selected vectors to the left
    if !isapprox(S.T[n+1, n], 0)                       # check if we are cutting along a schur block
        @error("conjugated eigenvector missing")
        display(S.T)
    end
    S.vectors[:,1:n], S.values[1:n]       # select first n vectors
end


text:
@testset "maximum_dims" begin
    ind1 = [1,2,3,4,5,6]
    @test NNlib.maximum_dims(ind1) == (6,)
    ind2 = [(3,4,5), (1,2,3), (2,3,9)]
    @test NNlib.maximum_dims(ind2) == (3,4,9)
    ind3 = [(3,4,5) (1,2,3) (2,3,9);
            (4,6,2) (5,3,2) (4,4,4)]
    @test NNlib.maximum_dims(ind3) == (5,6,9)
    ind4 = CartesianIndex.(
           [(3,4,5) (1,2,3) (2,3,9);
            (4,6,2) (5,3,2) (4,4,4)])
    @test NNlib.maximum_dims(ind4) == (5,6,9)
end

@testset "reverse_indices" begin
    res = [
        CartesianIndex.([(1,1), (2,3)]),
        CartesianIndex.([(1,2), (2,2)]),
        CartesianIndex.([(3,1), (1,3), (2,4), (3,4)]),
        CartesianIndex.([(2,1), (1,4)]),
        CartesianIndex.([(3,2), (3,3)])
    ]
    idx = [1 2 3 4;
           4 2 1 3;
           3 5 5 3]
    @test NNlib.reverse_indices(idx) == res
    @test NNlib.reverse_indices(idx) isa typeof(res)
    idx = [(1,) (2,) (3,) (4,);
           (4,) (2,) (1,) (3,);
           (3,) (5,) (5,) (3,)]
    @test NNlib.reverse_indices(idx) == res
    @test NNlib.reverse_indices(idx) isa typeof(res)
    idx = CartesianIndex.(
        [(1,) (2,) (3,) (4,);
        (4,) (2,) (1,) (3,);
        (3,) (5,) (5,) (3,)])
    @test NNlib.reverse_indices(idx) == res
    @test NNlib.reverse_indices(idx) isa typeof(res)
end


text:
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
function MOI.set(
    opt::Optimizer,
    ::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    opt.obj_sense = sense
    return nothing
end
MOI.get(opt::Optimizer, ::MOI.ObjectiveSense) = opt.obj_sense

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F<:_O_FUNCS} = true
function MOI.set(
    opt::Optimizer,
    ::MOI.ObjectiveFunction{F},
    obj::F,
) where {F<:_O_FUNCS}
    opt.form.obj = convert(CSAF, obj)
    return nothing
end


text:
using SafeTestsets
using Test

@testset "Common Utilities   " begin
    @safetestset "States" begin include("test_states.jl") end
end

@testset "Echo State Networks" begin
    @safetestset "ESN Input Layers" begin include("esn/test_input_layers.jl") end
    @safetestset "ESN Reservoirs" begin include("esn/test_reservoirs.jl") end
    @safetestset "ESN States" begin include("esn/test_states.jl") end
    @safetestset "ESN Train and Predict" begin include("esn/test_train.jl") end
    @safetestset "ESN Drivers" begin include("esn/test_drivers.jl") end
    @safetestset "ESN Non Linear Algos" begin include("esn/test_nla.jl") end
    @safetestset "Hybrid ESN" begin include("esn/test_hybrid.jl") end
end

@testset "CA based Reservoirs" begin
    @safetestset "RECA" begin include("reca/test_predictive.jl") end
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
export SingleObservable, of

import Base: ==
import Base: show

"""
    SingleObservable{D, H}(value::D, scheduler::H)

SingleObservable wraps single value of type `D` into a observable.

# Constructor arguments
- `value`: a single value to emit
- `scheduler`: scheduler-like object

# See also: [`of`](@ref), [`Subscribable`](@ref)
"""
struct SingleObservable{D, H} <: ScheduledSubscribable{D}
    value     :: D
    scheduler :: H
end

getscheduler(observable::SingleObservable) = observable.scheduler

function on_subscribe!(observable::SingleObservable, actor, scheduler)
    next!(actor, observable.value, scheduler)
    complete!(actor, scheduler)
    return voidTeardown
end

"""
    of(x; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

Creation operator for the `SingleObservable` that emits a single value x and then completes.

# Arguments
- `x`: value to be emmited before completion

# Examples

```jldoctest
using Rocket

source = of(1)
subscribe!(source, logger())
;

# output

[LogActor] Data: 1
[LogActor] Completed

```

See also: [`SingleObservable`](@ref), [`subscribe!`](@ref), [`logger`](@ref)
"""
of(x::T; scheduler::H = AsapScheduler()) where { T, H <: AbstractScheduler } = SingleObservable{T, H}(x, scheduler)

Base.:(==)(left::SingleObservable{D, H},  right::SingleObservable{D, H}) where { D, H } = left.value == right.value
Base.:(==)(left::SingleObservable,        right::SingleObservable)                      = false

Base.show(io::IO, ::SingleObservable{D, H}) where { D, H } = print(io, "SingleObservable($D, $H)")


text:
#=
To improve:
-Only one possible best alignment calculated... Revise so that all possible alignments are calculated
-Improve speed
-More user input options?
=#
using ArgParse

cli = ArgParseSettings()
@add_arg_table cli begin
	"--match"
		help = "match score"
		arg_type = Int
		default = 1
	"--mismatch"
		help = "mismatch score"
		arg_type = Int
		default = -1
	"--gap"
		help = "gap score"
		arg_type = Int
		default = -1
end
arg = parse_args(ARGS, cli)

seq1 = "AGGTCC"
seq2 = "AGTTC"

match = arg["match"]
mismatch = arg["mismatch"]
gap = arg["gap"]

function solvenw(seq_1, seq_2, match, mismatch, gap)
	#Info

	println("""seq1: $(seq1)
		seq2: $(seq2)
		""")

	println("""
		match score = $(match)
		mismatch score = $(mismatch)
		gap score = $(gap)
		""")

	#Initializing Matrix
	seq1_length = sizeof(seq1)
	seq2_length = sizeof(seq2)

	workspace = fill(0,(seq1_length+1,seq2_length+1))

	#Initialization
	for i in 2:seq1_length+1, j in 2:seq2_length+1
		workspace[i,1] = workspace[i-1,1]+gap
		workspace[1,j] = workspace[1,j-1]+gap
	end

	#Fill
	for i in 2:seq1_length+1
		for j in 2:seq2_length+1
			left_gap = workspace[i,j-1]+gap
			top_gap = workspace[i-1,j]+gap
			if seq1[i-1] == seq2[j-1]
				match_mismatch = workspace[i-1,j-1]+match
			else
				match_mismatch = workspace[i-1,j-1]+mismatch
			end
			workspace[i,j]=max(left_gap,top_gap,match_mismatch)
		end
	end

	#Display
	println("Matrix: ")
	display(workspace)

	#Trace Back
	x = seq1_length+1
	y = seq2_length+1

	traceback = String[]
	score = 0
	while x!=1
		trace_gap = workspace[x,y] - gap
		if seq1[x-1] == seq2[y-1]
			trace_diagnal = workspace[x,y] - match
			local_mm = match
		else
			trace_diagnal = workspace[x,y] - mismatch
			local_mm = mismatch
		end

		if trace_diagnal == workspace[x-1,y-1]
			push!(traceback, "match/mismatch")
			x = x-1
			y = y-1
			score += local_mm
		elseif trace_gap == workspace[x,y-1]
			push!(traceback,"left_gap")
			y = y-1
			score += gap
		elseif trace_gap == workspace[x-1,y]
			push!(traceback,"top_gap")
			x = x-1
			score += gap
		end
	end

	println("\r", traceback)

	aligned_seq1 = ""
	aligned_seq2 = ""
	counter_seq1 = 1
	counter_seq2 = 1

	for i in length(traceback):-1:1
		if traceback[i] == "match/mismatch"
			aligned_seq1 *= seq1[counter_seq1]
			aligned_seq2 *= seq2[counter_seq2]
			counter_seq1 += 1
			counter_seq2 += 1
		elseif traceback[i] == "left_gap"
			aligned_seq1 *= '-'
			aligned_seq2 *= seq2[counter_seq2]
			counter_seq2 += 1
		elseif traceback[i] == "top_gap"
			aligned_seq2 *= '-'
			aligned_seq1 *= seq1[counter_seq1]
			counter_seq1 +=1
		end
	end

	println(aligned_seq1)
	println(aligned_seq2)

	#Score
	println("score = $(score)")

end

solvenw(seq1, seq2, match, mismatch, gap)

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
using FactCheck
using MarketData
using DataFrames
using Base.Dates
using Compat

include("signals.jl")
include("orders.jl")
include("exchange_backtest.jl")
include("perfmetrics.jl")
include("backtest.jl")

# testing built-in trading strategies
include("teststrategy_goldencross.jl")
include("teststrategy_luxor.jl")

#exitstatus()


text:
module Animations

import Colors

# easings.jl
export Easing,
    eased,
    expin,
    expout,
    funcease,
    linear,
    mixed,
    multiplied,
    noease,
    opposite,
    polyin,
    polyio,
    polyout,
    saccadic,
    sineio,
    stepease

# keyframes.jl
export Keyframe

# animation.jl
export Animation,
    at,
    easings,
    keyframes,
    keyvalues,
    rel,
    timestamps

# loop.jl
export Loop

# sequence.jl
export Sequence

# animate_async.jl
export animate_async, stop

include("easings.jl")
include("keyframes.jl")
include("animation.jl")
include("interpolation.jl")
include("loop.jl")
include("sequence.jl")
include("animate_async.jl")

end


text:
@testset "utils tests" begin
    @testset "combine_fuzzy_sets tests" begin
        A = Float32[1.0 0.1; 0.4 1.0]
        union_res = Float32[1.0 0.46; 0.46 1.0]
        res = combine_fuzzy_sets(A, 1f0)
        @test isapprox(res, union_res)
        inter_res = [1.0 0.04; 0.04 1.0]
        res = combine_fuzzy_sets(A, 0f0)
        @test isapprox(res, inter_res)

    end

    @testset "fit_ab tests" begin
        @test all((1, 2) .== fit_ab(-1, 0, 1, 2))
    end
end


text:
# Copyright (c) 2013: Iain Dunning, Miles Lubin, and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestCWrapper

using Ipopt
using Test

function test_hs071()
    # hs071
    # min x1 * x4 * (x1 + x2 + x3) + x3
    # st  x1 * x2 * x3 * x4 >= 25
    #     x1^2 + x2^2 + x3^2 + x4^2 = 40
    #     1 <= x1, x2, x3, x4 <= 5
    # Start at (1,5,5,1)
    # End at (1.000..., 4.743..., 3.821..., 1.379...)

    function eval_f(x::Vector{Float64})
        return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]
    end

    function eval_g(x::Vector{Float64}, g::Vector{Float64})
        # Bad: g    = zeros(2)  # Allocates new array
        # OK:  g[:] = zeros(2)  # Modifies 'in place'
        g[1] = x[1] * x[2] * x[3] * x[4]
        return g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
    end

    function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
        # Bad: grad_f    = zeros(4)  # Allocates new array
        # OK:  grad_f[:] = zeros(4)  # Modifies 'in place'
        grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
        grad_f[2] = x[1] * x[4]
        grad_f[3] = x[1] * x[4] + 1
        return grad_f[4] = x[1] * (x[1] + x[2] + x[3])
    end

    function eval_jac_g(
        x::Vector{Float64},
        rows::Vector{Int32},
        cols::Vector{Int32},
        values::Union{Nothing,Vector{Float64}},
    )
        if values === nothing
            # Constraint (row) 1
            rows[1] = 1
            cols[1] = 1
            rows[2] = 1
            cols[2] = 2
            rows[3] = 1
            cols[3] = 3
            rows[4] = 1
            cols[4] = 4
            # Constraint (row) 2
            rows[5] = 2
            cols[5] = 1
            rows[6] = 2
            cols[6] = 2
            rows[7] = 2
            cols[7] = 3
            rows[8] = 2
            cols[8] = 4
        else
            # Constraint (row) 1
            values[1] = x[2] * x[3] * x[4]  # 1,1
            values[2] = x[1] * x[3] * x[4]  # 1,2
            values[3] = x[1] * x[2] * x[4]  # 1,3
            values[4] = x[1] * x[2] * x[3]  # 1,4
            # Constraint (row) 2
            values[5] = 2 * x[1]  # 2,1
            values[6] = 2 * x[2]  # 2,2
            values[7] = 2 * x[3]  # 2,3
            values[8] = 2 * x[4]  # 2,4
        end
        return
    end

    function eval_h(
        x::Vector{Float64},
        rows::Vector{Int32},
        cols::Vector{Int32},
        obj_factor::Float64,
        lambda::Vector{Float64},
        values::Union{Nothing,Vector{Float64}},
    )
        if values === nothing
            # Symmetric matrix, fill the lower left triangle only
            idx = 1
            for row in 1:4
                for col in 1:row
                    rows[idx] = row
                    cols[idx] = col
                    idx += 1
                end
            end
        else
            # Again, only lower left triangle
            # Objective
            values[1] = obj_factor * (2 * x[4])  # 1,1
            values[2] = obj_factor * (x[4])  # 2,1
            values[3] = 0                      # 2,2
            values[4] = obj_factor * (x[4])  # 3,1
            values[5] = 0                      # 3,2
            values[6] = 0                      # 3,3
            values[7] = obj_factor * (2 * x[1] + x[2] + x[3])  # 4,1
            values[8] = obj_factor * (x[1])  # 4,2
            values[9] = obj_factor * (x[1])  # 4,3
            values[10] = 0                     # 4,4

            # First constraint
            values[2] += lambda[1] * (x[3] * x[4])  # 2,1
            values[4] += lambda[1] * (x[2] * x[4])  # 3,1
            values[5] += lambda[1] * (x[1] * x[4])  # 3,2
            values[7] += lambda[1] * (x[2] * x[3])  # 4,1
            values[8] += lambda[1] * (x[1] * x[3])  # 4,2
            values[9] += lambda[1] * (x[1] * x[2])  # 4,3

            # Second constraint
            values[1] += lambda[2] * 2  # 1,1
            values[3] += lambda[2] * 2  # 2,2
            values[6] += lambda[2] * 2  # 3,3
            values[10] += lambda[2] * 2  # 4,4
        end
        return
    end

    n = 4
    x_L = [1.0, 1.0, 1.0, 1.0]
    x_U = [5.0, 5.0, 5.0, 5.0]

    m = 2
    g_L = [25.0, 40.0]
    g_U = [2.0e19, 40.0]

    prob = Ipopt.CreateIpoptProblem(
        n,
        x_L,
        x_U,
        m,
        g_L,
        g_U,
        8,
        10,
        eval_f,
        eval_g,
        eval_grad_f,
        eval_jac_g,
        eval_h,
    )

    prob.x = [1.0, 5.0, 5.0, 1.0]
    solvestat = Ipopt.IpoptSolve(prob)

    @test solvestat == 0
    @test prob.x[1] ≈ 1.0000000000000000 atol = 1e-5
    @test prob.x[2] ≈ 4.7429996418092970 atol = 1e-5
    @test prob.x[3] ≈ 3.8211499817883077 atol = 1e-5
    @test prob.x[4] ≈ 1.3794082897556983 atol = 1e-5
    @test prob.obj_val ≈ 17.014017145179164 atol = 1e-5

    # This tests callbacks.
    function intermediate(
        alg_mod::Cint,
        iter_count::Cint,
        obj_value::Float64,
        inf_pr::Float64,
        inf_du::Float64,
        mu::Float64,
        d_norm::Float64,
        regularization_size::Float64,
        alpha_du::Float64,
        alpha_pr::Float64,
        ls_trials::Cint,
    )
        return iter_count < 1  # Interrupts after one iteration.
    end

    Ipopt.SetIntermediateCallback(prob, intermediate)

    solvestat = Ipopt.IpoptSolve(prob)

    @test solvestat == 5

    # Test setting some options
    # String option
    println("\nString option")
    Ipopt.AddIpoptStrOption(prob, "hessian_approximation", "exact")
    @test_throws(
        ErrorException,
        Ipopt.AddIpoptStrOption(prob, "hessian_approximation", "badoption"),
    )
    println("\nInt option")
    # Int option
    @test Ipopt.AddIpoptIntOption(prob, "file_print_level", 3) === nothing
    @test_throws(
        ErrorException,
        Ipopt.AddIpoptIntOption(prob, "file_print_level", -1),
    )
    # Double option
    println("\nFloat option")
    Ipopt.AddIpoptNumOption(prob, "derivative_test_tol", 0.5)
    @test_throws(
        ErrorException,
        Ipopt.AddIpoptNumOption(prob, "derivative_test_tol", -1.0),
    )

    # Test opening an output file
    Ipopt.OpenIpoptOutputFile(prob, "blah.txt", 5)
    finalize(prob) # Needed before the `rm` on Windows.
    # unlink the output file
    rm("blah.txt")
    return
end

end  # TestCWrapper

runtests(TestCWrapper)


text:
##########################################################################
# Channels
abstract Channel

immutable SinglecarrierChannel <: Channel
    H::Array{Matrix{Complex128},2}
    Ns::Vector{Int}
    Ms::Vector{Int}
    K::Int
    I::Int

    large_scale_fading_factor::Matrix{Float64}
end
SinglecarrierChannel(H, Ns, Ms, K, I) =
    SinglecarrierChannel(H, Ns, Ms, K, I, ones(K, I))

immutable MulticarrierChannel <: Channel
    H::Array{Array{Complex128,3},2}
    Ns::Vector{Int}
    Ms::Vector{Int}
    K::Int
    I::Int
    Lc::Int

    large_scale_fading_factor::Matrix{Float64}
end
MulticarrierChannel(H, Ns, Ms, K, I, Lc) =
    MulticarrierChannel(H, Ns, Ms, K, I, Lc, ones(K, I))

get_channel_gains_dB(channel::SinglecarrierChannel) =
    Float64[ 20*log10(vecnorm(channel.H[k,i])) for k = 1:channel.K, i = 1:channel.I ]

function get_average_channel_gains_dB(channel::MulticarrierChannel)
    average_channel_gains = zeros(Float64, channel.K, channel.I)

    for k = 1:channel.K
        for i = 1:channel.I
            for l = 1:channel.Lc
                average_channel_gains[k,i] += (1/channel.Lc)*vecnorm(channel.H[k,i,l])^2
            end
        end
    end

    return 10*log10(average_channel_gains)
end


text:
using Test
using StaticArrays
using RootedTrees

trees_array = (rootedtree([1,2,3]),
               rootedtree([1,2,3]),
               rootedtree([1,2,2]),
               rootedtree([1,2,3,3]))

for (t1,t2,t3,t4) in (trees_array,)
  @test t1 == t1
  @test t1 == t2
  @test !(t1 == t3)
  @test !(t1 == t4)

  @test t3 < t2    && t2 > t3
  @test !(t2 < t3) && !(t3 > t2)
  @test t1 < t4    && t4 > t1
  @test !(t4 < t1) && !(t1 > t4)
  @test t1 <= t2   && t2 >= t1
  @test t2 <= t2   && t2 >= t2

  println(devnull, t1)
  println(devnull, t2)
  println(devnull, t3)
  println(devnull, t4)
end

t = rootedtree([1])
@test order(t) == 1
@test σ(t) == 1
@test γ(t) == 1
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2])
@test order(t) == 2
@test σ(t) == 1
@test γ(t) == 2
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 2])
@test order(t) == 3
@test σ(t) == 2
@test γ(t) == 3
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3])
@test order(t) == 3
@test σ(t) == 1
@test γ(t) == 6
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 2, 2])
@test order(t) == 4
@test σ(t) == 6
@test γ(t) == 4
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 2, 3])
@inferred RootedTrees.subtrees(t)
@test order(t) == 4
@test σ(t) == 1
@test γ(t) == 8
@test α(t) == 3
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 3])
@test order(t) == 4
@test σ(t) == 2
@test γ(t) == 12
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 4])
@test order(t) == 4
@test σ(t) == 1
@test γ(t) == 24
@test α(t) == 1

t = rootedtree([1, 2, 2, 2, 2])
@test order(t) == 5
@test σ(t) == 24
@test γ(t) == 5
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 2, 2, 3])
@test order(t) == 5
@test σ(t) == 2
@test γ(t) == 10
@test α(t) == 6
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 2, 3, 3])
@test order(t) == 5
@test σ(t) == 2
@test γ(t) == 15
@test α(t) == 4

t = rootedtree([1, 2, 2, 3, 4])
@test order(t) == 5
@test σ(t) == 1
@test γ(t) == 30
@test α(t) == 4
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 2, 3])
@test order(t) == 5
@test σ(t) == 2
@test γ(t) == 20
@test α(t) == 3
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 3, 3])
@test order(t) == 5
@test σ(t) == 6
@test γ(t) == 20
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 3, 4])
@test order(t) == 5
@test σ(t) == 1
@test γ(t) == 40
@test α(t) == 3
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 4, 4])
@test order(t) == 5
@test σ(t) == 2
@test γ(t) == 60
@test α(t) == 1
@test β(t) == α(t)*γ(t)

t = rootedtree([1, 2, 3, 4, 5])
@test order(t) == 5
@test σ(t) == 1
@test γ(t) == 120
@test α(t) == 1
@test β(t) == α(t)*γ(t)

# see butcher2008numerical, Table 302(I)
number_of_rooted_trees = [1, 1, 2, 4, 9, 20, 48, 115, 286, 719]
for order in 1:10
  num = 0
  for t in RootedTreeIterator(order)
    num += 1
  end
  @test num == number_of_rooted_trees[order] == count_trees(order)
end

# Runge-Kutta method SSPRK33 of order 3
A = [0 0 0; 1 0 0; 1/4 1/4 0]
b = [1/6, 1/6, 2/3]
c = A * fill(1, length(b))
for order in 1:3
  for t in RootedTreeIterator(order)
    @test residual_order_condition(t, A, b, c) ≈ 0 atol=eps()
  end
end
let order=4
  res = 0.
  for t in RootedTreeIterator(order)
    res += abs(residual_order_condition(t, A, b, c))
  end
  @test res > 10*eps()
end

A = @SArray [0 0 0; 1 0 0; 1/4 1/4 0]
b = @SArray [1/6, 1/6, 2/3]
c = A * SVector(1, 1, 1)
for order in 1:3
  for t in RootedTreeIterator(order)
    @test residual_order_condition(t, A, b, c) ≈ 0 atol=eps()
  end
end
let order=4
  res = 0.
  for t in RootedTreeIterator(order)
    res += abs(residual_order_condition(t, A, b, c))
  end
  @test res > 10*eps()
end


text:
"""
    setStopCondition(; genj::GeneticJuliaStructure = GenJ, maxEvaluations::Integer = -1,
                     maxIterations::Integer = -1, maxIterNotImproving::Integer = -1,
                     maxTime::AbstractFloat = Inf)

Sets the stop conditions for `genj`, however if none are introduced, `maxIterations`
will be set to 100 by default. This function constructs a structure of type
`StopCondition` and adds it to the main structure.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.
- `maxEvaluations::Integer = -1`: number of calls to compute fitness before the
    algorithm stops. If it is set to a `negative value`, the verification over
    this parameter will be `useless`, as default.
- `maxIterations::Integer = -1`: number of iterations before the algorithm stops.
    If it is set to a `negative value`, the verification over this parameter
    will be `useless`, as default.
- `maxIterNotImproving::Integer = -1`: number of iterations in which the best
    individual doesn't improve before the algorithm stops. If it is set to a
    `negative value`, the verification over this parameter will be `useless`, as default.
- `maxTime::AbstractFloat = Inf`: time, in seconds, passed since the beginning of the
    experiment before the algorithm stops. If it is set to `inf`, the verification
    over this parameter will be `useless`, as default.

See also: [`StopCondition`](@ref)
"""
function setStopCondition(; genj::GeneticJuliaStructure = GenJ, maxEvaluations::Integer = -1,
                          maxIterations::Integer = -1, maxIterNotImproving::Integer = -1,
                          maxTime::AbstractFloat = Inf)

    if maxTime < 0
        error("maxTime must be a positive number")
    end

    if (maxEvaluations < 1) && (maxIterations < 1) && (maxIterNotImproving < 1) && (maxTime == Inf)
        error("At least one of the stop conditions must be set to a value that allows checking")
    end

    stopCondition = StopCondition(maxEvaluations, maxIterations, maxIterNotImproving, maxTime, 0, 0, 0, 0.0)
    genj._stopCondition = stopCondition
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
#=
TODO: Each graph representation is assumed to be static... so
=#

abstract type AbstractDirectedGraph end
abstract type AbstractDirectedAcyclicGraph <: AbstractDirectedGraph end

const AbstractDG = AbstractDirectedGraph
const AbstractDAG = AbstractDirectedAcyclicGraph

@enum Linearity LIN_CONSTANT LIN_LINEAR LIN_PIECEWISE_LINEAR LIN_NONLINEAR
@enum VariableType VT_BIN VT_INT VT_CONT

function _variable_count(g::AbstractDG)::Int
    error("Variable count not defined for graph type = $(typeof(g))")
end

# added id field from JuMP UserOperatorRegistry, expect more extensive changes in future.
struct OperatorRegistry
    multivariate_id::Vector{Symbol}
    multivariate_operator_to_id::Dict{Symbol,Int}
    multivariate_operator_evaluator::Vector{MOI.AbstractNLPEvaluator}
    univariate_operator_id::Vector{Symbol}
    univariate_operator_to_id::Dict{Symbol,Int}
    univariate_operator_f::Vector{Any}
    univariate_operator_fprime::Vector{Any}
    univariate_operator_fprimeprime::Vector{Any}
end
function OperatorRegistry()
    return OperatorRegistry(
        Symbol[],
        Dict{Symbol,Int}(),
        MOI.AbstractNLPEvaluator[],
        Symbol[],
        Dict{Symbol,Int}(),
        [],
        [],
        [],
    )
end

include(joinpath(@__DIR__, "expressions.jl"))
include(joinpath(@__DIR__, "node.jl"))
include(joinpath(@__DIR__, "utilities.jl"))
include(joinpath(@__DIR__, "graphs", "directed_tree.jl"))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "Methods" begin
    @testset "Correctness" begin
        # Finite difference methods to test.
        methods = [forward_fdm, backward_fdm, central_fdm]

        # The different floating point types to try.
        types = [Float32, Float64]

        # The different functions to evaluate (`.f`), their first derivative at 1 (`.d1`),
        # and their second derivative at 1 (`.d2`).
        fs = [
            (f=sin, d1=cos(1), d2=-sin(1), only_forward=false),
            (f=exp, d1=exp(1), d2=exp(1), only_forward=false),
            (f=abs2, d1=2, d2=2, only_forward=false),
            (f=x -> sqrt(x + 1), d1=0.5 / sqrt(2), d2=-0.25 / 2^(3/2), only_forward=true),
        ]

        # Test all combinations of the above settings, i.e. differentiate all functions
        # using all methods and data types.
        @testset "f=$(f.f), method=$m, type=$(T)" for f in fs, m in methods, T in types
            # Check if only `forward_fdm` is allowed.
            f.only_forward && m != forward_fdm && continue

            @testset "method-order=$order" for order in [1, 2, 3, 4, 5]
                @test m(order, 0; adapt=2)(f.f, T(1)) isa T
                @test m(order, 0; adapt=2)(f.f, T(1)) == T(f.f(1))
            end

            @test m(10, 1)(f.f, T(1)) isa T
            @test m(10, 1)(f.f, T(1)) ≈ T(f.d1)

            @test m(10, 2)(f.f, T(1)) isa T
            T == Float64 && @test m(10, 2)(f.f, T(1)) ≈ T(f.d2)
        end
    end

    @testset "Accuracy" begin
        # Finite difference methods to test.
        methods = [forward_fdm, backward_fdm, central_fdm]

        # `f`s, `x`s, the derivatives of `f` at `x`, and a factor that loosens tolerances.
        fs = [
            # See #124.
            (f=x -> 0, x=0, d=0,           atol=0,     atol_central=0),
            (f=x -> x, x=0, d=1,           atol=5e-15, atol_central=5e-16),
            (f=exp, x=0, d=1,              atol=5e-13, atol_central=5e-14),
            (f=sin, x=0, d=1,              atol=1e-14, atol_central=5e-16),
            (f=cos, x=0, d=0,              atol=5e-14, atol_central=5e-14),
            (f=exp, x=1, d=exp(1),         atol=5e-12, atol_central=5e-13),
            (f=sin, x=1, d=cos(1),         atol=5e-13, atol_central=5e-14),
            (f=cos, x=1, d=-sin(1),        atol=5e-13, atol_central=5e-14),
            (f=sinc, x=0, d=0,             atol=5e-12, atol_central=5e-14),
            # See #125. The tolerances are lower because `cosc` is hard.
            (f=cosc, x=0, d=-(pi ^ 2) / 3, atol=5e-9,  atol_central=5e-10)
        ]
        @testset "f=$(f.f), method=$m" for f in fs, m in methods
            atol = m == central_fdm ? f.atol_central : f.atol
            @test m(6, 1)(f.f, f.x) ≈ f.d rtol=0 atol=atol
            @test m(7, 1)(f.f, f.x) ≈ f.d rtol=0 atol=atol
            if m == central_fdm
                if f.f == cosc
                    # `cosc` is hard.
                    @test m(14, 1)(f.f, f.x) ≈ f.d atol=1e-13
                else
                    @test m(14, 1)(f.f, f.x) ≈ f.d atol=5e-15
                end
            end
        end
    end

    @testset "Test custom grid" begin
        @test FiniteDifferenceMethod([-2, 0, 5], 1)(sin, 1) ≈ cos(1)
    end

    @testset "Test allocations" begin
        m = central_fdm(5, 2, adapt=2)
        @test @ballocated($m(sin, 1)) == 0
    end

    # Integration test to ensure that Integer-output functions can be tested.
    @testset "Integer output" begin
        @test isapprox(central_fdm(5, 1)(x -> 5, 0), 0; rtol=1e-12, atol=1e-12)
    end

    @testset "Adaptation improves estimate" begin
        @test abs(forward_fdm(5, 1, adapt=0)(log, 0.001) - 1000) > 1
        @test forward_fdm(5, 1, adapt=1)(log, 0.001) ≈ 1000
    end

    @testset "Limiting step size" begin
        @test !isfinite(central_fdm(5, 1, max_range=0)(abs, 0.001))
        @test central_fdm(10, 1, max_range=9e-4)(log, 1e-3) ≈ 1000
        @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0
    end

    @testset "Accuracy at high orders and high adaptation (issue #64)" begin
        # Regression test against issues with precision during computation of coefficients.
        @test central_fdm(9, 5, adapt=4)(exp, 1.0) ≈ exp(1) atol=2e-7
        @test central_fdm(15, 5, adapt=2)(exp, 1.0) ≈ exp(1) atol=1e-10
        poly(x) = 4x^3 + 3x^2 + 2x + 1
        @test central_fdm(9, 3, adapt=4)(poly, 1.0) ≈ 24 atol=1e-11
    end

    @testset "Printing FiniteDifferenceMethods" begin
        @test sprint(show, "text/plain", central_fdm(2, 1)) == """
            FiniteDifferenceMethod:
              order of method:       2
              order of derivative:   1
              grid:                  [-1, 1]
              coefficients:          [-0.5, 0.5]
            """
    end

    @testset "_is_symmetric" begin
        # Test odd grids:
        @test FiniteDifferences._is_symmetric(SVector{5}(2, 1, 0, 1, 2))
        @test !FiniteDifferences._is_symmetric(SVector{5}(2, 1, 0, 3, 2))
        @test !FiniteDifferences._is_symmetric(SVector{5}(4, 1, 0, 1, 2))

        # Test even grids:
        @test FiniteDifferences._is_symmetric(SVector{4}(2, 1, 1, 2))
        @test !FiniteDifferences._is_symmetric(SVector{4}(2, 1, 3, 2))
        @test !FiniteDifferences._is_symmetric(SVector{4}(4, 1, 1, 2))

        # Test zero at centre:
        @test FiniteDifferences._is_symmetric(SVector{5}(2, 1, 4, 1, 2))
        @test !FiniteDifferences._is_symmetric(SVector{5}(2, 1, 4, 1, 2), centre_zero=true)
        @test FiniteDifferences._is_symmetric(SVector{4}(2, 1, 1, 2), centre_zero=true)

        # Test negation of a half:
        @test !FiniteDifferences._is_symmetric(SVector{4}(2, 1, -1, -2))
        @test FiniteDifferences._is_symmetric(SVector{4}(2, 1, -1, -2), negate_half=true)
        @test FiniteDifferences._is_symmetric(SVector{5}(2, 1, 0, -1, -2), negate_half=true)
        @test FiniteDifferences._is_symmetric(SVector{5}(2, 1, 4, -1, -2), negate_half=true)
        @test !FiniteDifferences._is_symmetric(
            SVector{5}(2, 1, 4, -1, -2);
            negate_half=true,
            centre_zero=true
        )

        # Test symmetry of `central_fdm`.
        for p in 2:10
            m = central_fdm(p, 1)
            @test FiniteDifferences._is_symmetric(m)
        end

        # Test asymmetry of `forward_fdm` and `backward_fdm`.
        for p in 2:10
            for f in [forward_fdm, backward_fdm]
                m = f(p, 1)
                @test !FiniteDifferences._is_symmetric(m)
            end
        end
    end

    @testset "extrapolate_fdm" begin
        # Also test an `Integer` argument as input.
        for x in [1, 1.0]
            for f in [forward_fdm, central_fdm, backward_fdm]
                estimate, _ = extrapolate_fdm(f(4, 3), exp, x, contract=0.8)
                @test estimate ≈ exp(1.0) atol=1e-7
            end
        end
    end
end


text: Dict{String, Any}("" => Dict{String, Any}("value" => "blank", "type" => "string"))


text:
# Batch simulations

"
$SIGNATURES

Settings for batch data generation
"
Base.@kwdef struct BatchDataSettings
  "Master Random Number Generator"
  rng::MersenneTwister
  "Experiment name - used for results subdirectory"
  name::String
  "Process labels"
  proclabels::Vector{Symbol}
  "Average sampling time"
  Tav::Float64
  "Minimal time difference between samples"
  Tmin::Float64
  "Number of observations"
  Ns::Vector{Int}
  "Number of simulation runs"
  Nrun::Int
end

function Base.show(io::IO,x::BatchDataSettings)
  for k in fieldnames(BatchDataSettings)
    v = getfield(x,k)
    vp = k==:rng ? "$(typeof(v))(...)" : v
    print(io,@sprintf("% 10s : %s\n",k,vp))
  end
end

settings_resultsdir(name) = resultsdir("ICMLA",name)

"""
$SIGNATURES

Batch data generation
"""
function batch_data(S::BatchDataSettings)
  config = (proclabel=S.proclabels,N=S.Ns,i=1:S.Nrun)
  partial_spec = vec([(;(keys(config).=>vs)...) for vs in Iterators.product(config...)])

  # generate datas
  dfrows = Vector{Any}(undef,length(partial_spec))
  prog = Progress(10)
  @threads for index in 1:length(partial_spec)
    ps = partial_spec[index]
    spec = (rng=copy(S.rng),N=ps.N,proclabel=ps.proclabel,Tav=S.Tav,Tmin=S.Tmin)
    proc = testprocess(spec.proclabel)
    data = rand(S.rng,proc,spec.N,spec.Tav,spec.Tmin)
    dfrows[index] = merge((index=index,),spec,(data=data,))
    next!(prog)
  end
  df = DataFrame(dfrows)

  # Save to file
  f = joinpath(settings_resultsdir(S.name),"data.jld")
  R = (df=df,settings=S,status=systemstatus())
  serialize(f,R); printstyled("\nSaved to file:\n$f\n",color=:blue,bold=true)

  return R
end

"
Settings for batch maximum likelihood estimation

$FIELDS
"
Base.@kwdef struct BatchMlEstimation
  "Model order"
  p::Int
  "Time scale of interest (Kullback-Leiber Discrepancy)"
  Tkld::Float64
  "Kalman engines"
  engines::Vector
  "Initiate at actual process?"
  init_actuals::Vector{Bool}
end

function Base.show(io::IO,x::BatchMlEstimation)
  for k in fieldnames(BatchMlEstimation)
    v = getfield(x,k)
    print(io,@sprintf("% 20s : %s\n",k,v))
  end
end

"
$SIGNATURES

Single maximum likelihood estimate for batch estimation. Returns
`nothing` if estimation trace cannot be computed.
"
function batch_mle_trace_single(S::BatchMlEstimation,spec,row_data;throw_errors=false)
  # Initial estimators
  tor = MlEstimator()
  Ts1 = [0.5,1,2]*S.Tkld
  Ts2 = [0.5,1,2]*row_data.Tav
  Ts = sort(unique([Ts1;Ts2]))
  init_tors = default_init_tors(row_data.rng,Ts)

  if spec.init_actual
    proc = testprocess(row_data.proclabel)
    if !(proc isa SdeProcess)
      # cannot initiate from a non-SDE process such as the white noise or squared exponential
      return nothing
    end
    push!(init_tors,InitialModel(proc))
  end

  problem = ProblemLike(row_data.data;p=S.p,engine=spec.engine)

  _,trace = estimate(tor,problem;init_tors=init_tors,ret_trace=true,throw_errors=throw_errors)

  return trace
end

function batch_mle_trace_body(S::BatchMlEstimation,df_data,specs)
  Ndata = nrow(df_data)
  df_arr = Array{Any,2}(undef,Ndata,length(specs))
  df_arr[:] .= missing
  @threads for idata in 1:Ndata
    row_data = df_data[idata,:]
    @threads for i in 1:length(specs)
      spec = specs[i]
      mle_trace = batch_mle_trace_single(S,spec,row_data)
      if !isnothing(mle_trace)
        df_arr[idata,i] = merge((index=row_data.index,),spec,(mle_trace=mle_trace,))
      end
    end
    print("[$(threadid())] $idata/$Ndata, ")
  end
  df = DataFrame(skipmissing(vec(df_arr)))
  return df
end

"
$SIGNATURES

Batch Maximum likelihood trace
"
function batch_mle_trace(S::BatchMlEstimation,df_data,results_dir)
  config = (init_actual=S.init_actuals,engine=S.engines)
  specs = vec([(;(keys(config).=>vs)...) for vs in Iterators.product(config...)])
  printstyled("$(length(specs)) Maximum likelihood configurations:\n",color=:blue)
  println.(specs)
  Ndata = nrow(df_data); Nspec = length(specs)
  t = Ndata*Nspec
  println("ML configurations: $Nspec Dataset count : $Ndata ⟹  Total estimates : $t")

  df = batch_mle_trace_body(S,df_data,specs)

  # Save to file
  f = joinpath(results_dir,"mle-trace.jld")
  output_mle_trace = (df=df,settings=S,status=systemstatus())
  serialize(f,output_mle_trace)
  printstyled("\nSaved to file:\n$f\n",color=:blue,bold=true)

  return output_mle_trace
end

function batchsettings(f)
  bs = OrderedDict{Symbol,Any}()
  for (k,v) in TOML.parsefile(f)
    k = Symbol(k)
    if k==:proclabel v=Symbol(v)  end
    bs[k] = v
  end
  batchset = (;bs...)
  return batchset
end

function batch_posterior_samples(S,output_mle)
  estimates = output_mle.estimates
  # Tav = output_mle.settings.settings_data.Tav
  ## Body
  Nest = length(estimates)
  if !isnothing(S.Nest)
    @info "Sampling is done for $(S.Nest) out of $Nest maximum likelihood estimates"
    Nest = S.Nest
  end
  samples = Vector{Any}(undef,Nest)
  @threads for iest in 1:Nest
    estml = estimates[iest]
    data_c = estml.data_c
    tor = PmeanEstimator(;rng=data_c.rng,S.tor_kwargs...)
    proc_ml = estml.rmle.est
    samples[iest] = if proc_ml isa WhiteNoise
      (estml=estml,procs=nothing,samtrace=nothing)
    else
      p = EstimateSDE.order(proc_ml)
      problem = ProblemProb(data_c.data;p=p,engine=SdeRootsEngine,prior=ReferencePrior(Tav))
      _,procs,samtrace = estimate(tor,problem,proc_ml)``
      (estml=estml,procs=procs,samtrace=samtrace)
    end
    @info "Sampling for $iest completed"
  end

  # Save to file
  # TODO Make resdir an input argument to the function
  resdir = settings_resultsdir(output_mle.settings.settings_data)
  f = joinpath(resdir,"posterior-samples.jld")
  R = (samples=samples,settings=S,status=systemstatus())
  serialize(f,R)
  printstyled("\nSaved to file:\n$f\n",color=:blue,bold=true)

  return R
end

colors = [:blue, :light_red, :cyan, :magenta, :green, :light_blue, :light_magenta,
    :light_cyan,:light_green, :light_yellow,  :red, :yellow]
function tree(t::NamedTuple,depth=1;ind="————",maxdepth=Inf)
  color = colors[depth]
  identation = repeat(ind,depth)*" "
  for (k,v) in pairs(t)
    printstyled("$identation$k = ";color=color)
    if v isa NamedTuple
      if depth==maxdepth printstyled("…\n";color=color); return end
      printstyled("(\n";color=color)
      tree(v,depth+1;ind=ind,maxdepth=maxdepth)
      printstyled("$identation)\n";color=color)
    else
      println(repr(v)[1:min(end,80)])
    end
  end
  return nothing
end

"
$SIGNATURES

Estimates derived from maximum likelihood traces using `trace_estimates`.
"
function derived_estimates(df_mle_trace,results_dir)
  deriv_row = []
  for row_mle_trace in eachrow(df_mle_trace)
    table_trace_estimate = trace_estimates(row_mle_trace.mle_trace)
    for r in table_trace_estimate
      ei = r.estinfo
      tinfo = (p=ei.p,init_tor=ei.init_tor,est=ei.est)
      row_deriv = merge((;row_mle_trace...),(method=r.method,pmethod=r.pmethod),tinfo)
      push!(deriv_row,row_deriv)
    end
  end
  df = DataFrame(deriv_row)

  ## Save to file
  f = joinpath(results_dir,"derived-estimates.jld")
  R = (settings=nothing,df=df,status=systemstatus())
  serialize(f,R)
  printstyled("\nSaved to file:\n$f\n",color=:blue,bold=true)

  return R
end

"
$SIGNATURES

Batch model evaluation
"
function batch_model_evaluation(df_data,Smle_trace,df_derived,results_dir)
  ## KLD per estimate
  Nderived = nrow(df_derived)
  Nmetrics = 3 # Number of metrics
  rows_eval = Array{Any,2}(undef,Nmetrics,Nderived)
  @threads for i_derived in 1:Nderived
    row_derived = df_derived[i_derived,:]
    row_data = df_data[row_derived.index,:]
    est = row_derived.est
    proc = testprocess(row_data.proclabel)

    for i_metric in 1:Nmetrics
      label,T,t_kld = if i_metric==3
        t = getfield.(row_data.data,:t)
        "same",missing,t
      else
        T = i_metric==1 ? Smle_trace.Tkld : row_data.Tav
        duration = row_data.Tav*row_data.N
        t = (0:T:duration)
        "regular",T,t
      end
      kld = safe_kldiv(proc,est,t_kld,missing)
      row_eval = merge(row_derived,(label=label,T=T,kld=kld))
      rows_eval[i_metric,i_derived] = row_eval
    end
    print("[$(threadid())] $i_derived/$Nderived, ")
  end
  df = DataFrame(vec(rows_eval))

  ## Save to file
  f = joinpath(results_dir,"model-evaluation.jld")
  R = (settings=nothing,df=df,status=systemstatus())
  serialize(f,R); printstyled("\nSaved to file\n$f\n: ",color=:blue,bold=true)

  return R
end

const icmla_paper_dir = abspath(joinpath(@__DIR__,"..","ICMLA2020"))

"Copy to ICMLA paper figures directory"
function copy_icmla_figures(figure_dir,choice=nothing)
  if isnothing(choice)
    choice = request("Copy figures to ICMLA paper figures directory?",RadioMenu(["yes","no"]))
  end
  if choice==1
    dest_dir = joinpath(icmla_paper_dir,"figures")
    for fn in readdir(figure_dir)
      if !(endswith(fn,"pdf") || endswith(fn,"tex")) continue end
      src = joinpath(figure_dir,fn)
      if !isfile(src) continue end
      dst = joinpath(dest_dir,fn)
      cp(src,dst,force=true)
      println("Copied to ",dst)
    end
    println("\nDone.")
  end
  return choice
end


text:
@testset "group" begin
    @test group(identity, 11:20)::Dictionary == Dictionary(11:20, (x->[x]).(11:20))
    @test group(iseven, 1:10)::Dictionary == dictionary([false => [1,3,5,7,9], true => [2,4,6,8,10]])

    @test group(iseven, x -> x*2, 1:10)::Dictionary == dictionary([false => [2,6,10,14,18], true => [4,8,12,16,20]])

    @test group((x,y) -> iseven(x+y), (x,y) -> x, 1:10, [1,3,4,2,5,6,4,2,3,9])::Dictionary == dictionary([true => [1,4,5,6,8,9], false => [2,3,7,10]])
end

@testset "grouponly" begin
    @test grouponly(identity, 11:20)::Dictionary == Dictionary(11:20, 11:20)
    @test_throws IndexError grouponly(iseven, 11:20)
end

@testset "groupunique" begin
    @test groupunique(identity, 11:20)::Dictionary{Int, Indices{Int}} == Dictionary(11:20, map(x -> Indices([x]), 11:20))
    @test groupunique(iseven, 11:20)::Dictionary{Bool, Indices{Int}} == dictionary([false => Indices([11,13,15,17,19]), true => Indices([12,14,16,18,20])])
end

@testset "groupfind" begin
    @test groupfind(identity, 11:20) == Dictionary(11:20, (x->[x]).(1:10))
    @test groupfind(iseven, 11:20) == dictionary([false => [1,3,5,7,9], true => [2,4,6,8,10]])
end

@testset "groupview" begin
    @test groupview(identity, 11:20)::GroupDictionary == group(identity, 11:20)::Dictionary
    @test groupview(iseven, 11:20)::GroupDictionary == group(iseven, 11:20)::Dictionary
end

@testset "groupreduce" begin
    @test groupreduce(identity, +, 1:10) == Dictionary(1:10, 1:10)
    @test groupreduce(iseven, +, 1:10) == dictionary([false => 25, true => 30])

    @test groupreduce(iseven, x -> x*2, +, 1:10) == dictionary([false => 50, true => 60])

    @test groupreduce(iseven, x -> x*2, +, 1:10; init=10) == dictionary([false => 60, true => 70])

    @test groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, 1:10; init=10) == dictionary([true => 120])
    @test groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, [1,3,4,2,5,6,4,2,3,9]; init=10) == dictionary([true => 62, false => 52])

    @test groupcount(iseven, 1:10) == dictionary([false => 5, true => 5])
    @test groupsum(iseven, 1:10) == dictionary([false => 25, true => 30])
    @test groupprod(iseven, 1:10) == dictionary([false => 1*3*5*7*9, true => 2*4*6*8*10])
    @test groupfirst(iseven, 1:10) == dictionary([false => 1, true => 2])
    @test grouplast(iseven, 1:10) == dictionary([false => 9, true => 10])
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# add new event to list, maintain sorting by time
function addEvent!(eventList::Vector{Event};
	parentEvent::Union{Event,Nothing} = nothing, form::EventForm = nullEvent, time::Float = nullTime,
	ambulance::Union{Ambulance,Nothing} = nothing, call::Union{Call,Nothing} = nothing, station::Union{Station,Nothing} = nothing,
	addEventToAmb::Bool = true)

	event = Event()
	event.parentIndex = parentEvent == nothing ? nullIndex : parentEvent.index
	event.form = form
	event.time = time
	event.ambIndex = ambulance == nothing ? nullIndex : ambulance.index
	event.callIndex = call == nothing ? nullIndex : call.index
	event.stationIndex = station == nothing ? nullIndex : station.index

	# set next event of ambulance
	if addEventToAmb && ambulance != nothing
		ambulance.event = event
	end

	# find where to insert event into list
	# maintain sorting by time, events nearer to end of list are sooner
	i = something(findlast(e -> e.time >= event.time, eventList), 0) + 1
	insert!(eventList, i, event)

	if checkMode
		# check time ordering of events
		for i = 1:length(eventList)-1
			@assert(eventList[i].time >= eventList[i+1].time)
		end
		# @assert(issorted(eventList, by = e -> e.time, rev = true)) # slow
	end

	return event
end

# for adding call arrival to event list
function addEvent!(eventList::Vector{Event}, call::Call)
	addEvent!(eventList; form = callArrives, time = call.arrivalTime, call = call)
end

# get next event from list
function getNextEvent!(eventList::Vector{Event})
	return length(eventList) > 0 ? pop!(eventList) : Event()
end

function printEvent(event::Event)
	println("Event:")
	println("form = ", string(event.form))
	println("time = ", event.time)
	println("ambIndex = ", event.ambIndex)
	println("callIndex = ", event.callIndex)
end

# delete event in eventList
function deleteEvent!(eventList::Vector{Event}, event::Event)
	i = something(findfirst(isequal(event), eventList))
	@assert(i != nothing)
	@assert(findnext(e -> e == event, eventList, i + 1) == nothing)
	deleteat!(eventList, i)
end


text:
"""
    hs_boltz_kn(mu_ref, alpha)

Calculate effective Knudsen number for fast spectral method with hard sphere (HS) model

"""
hs_boltz_kn(mu_ref, alpha) =
    64 * sqrt(2.0)^alpha / 5.0 * gamma((alpha + 3) / 2) * gamma(2.0) * sqrt(pi) * mu_ref

"""
    kernel_mode(
        M::TI,
        umax::TR,
        vmax::TR,
        du::TR,
        dv::TR,
        unum::TI,
        vnum::TI;
        quad_num = 64::TI,
    ) where {TI<:Integer,TR<:Real}

    kernel_mode(
        M::I,
        umax::R,
        vmax::R,
        wmax::R,
        du::R,
        dv::R,
        dw::R,
        unum::I,
        vnum::I,
        wnum::I,
        alpha::R;
        quad_num = 64::I,
    ) where {I<:Integer,R<:Real}

Calculate collision kernel for fast spectral method

"""
function kernel_mode(
    M::I,
    umax::R,
    vmax::R,
    wmax::R,
    du::R,
    dv::R,
    dw::R,
    unum::I,
    vnum::I,
    wnum::I,
    alpha::R;
    quad_num = 64::I,
) where {I<:Integer,R<:Real}

    supp = sqrt(2.0) * 2.0 * max(umax, vmax, wmax) / (3.0 + sqrt(2.0))

    fre_vx = range(-π / du, (unum ÷ 2 - 1) * 2.0 * π / unum / du, length = unum)
    fre_vy = range(-π / dv, (vnum ÷ 2 - 1) * 2.0 * π / vnum / dv, length = vnum)
    fre_vz = range(-π / dw, (wnum ÷ 2 - 1) * 2.0 * π / wnum / dw, length = wnum)

    # abscissa, gweight = gausslegendre(quad_num)
    # @. abscissa = (0. * (1. - abscissa) + supp * (1. + abscissa)) / 2
    # @. gweight *= (supp - 0.) / 2

    abscissa, gweight = lgwt(quad_num, 0, supp)

    phi = zeros(unum, vnum, wnum, M * (M - 1))
    psi = zeros(unum, vnum, wnum, M * (M - 1))
    phipsi = zeros(unum, vnum, wnum)
    for loop = 1:M-1
        theta = π / M * loop
        for loop2 = 1:M
            theta2 = π / M * loop2
            idx = (loop - 1) * M + loop2
            for k = 1:wnum, j = 1:vnum, i = 1:unum
                s =
                    fre_vx[i] * sin(theta) * cos(theta2) +
                    fre_vy[j] * sin(theta) * sin(theta2) +
                    fre_vz[k] * cos(theta)
                # phi
                int_temp = 0.0
                for id = 1:quad_num
                    int_temp +=
                        2.0 * gweight[id] * cos(s * abscissa[id]) * (abscissa[id]^alpha)
                end
                phi[i, j, k, idx] = int_temp * sin(theta)
                # psi
                s = fre_vx[i]^2 + fre_vy[j]^2 + fre_vz[k]^2 - s^2
                if s <= 0.0
                    psi[i, j, k, idx] = π * supp^2
                else
                    s = sqrt(s)
                    bel = supp * s
                    bessel = besselj(1, bel)
                    psi[i, j, k, idx] = 2.0 * π * supp * bessel / s
                end
                # phipsi
                phipsi[i, j, k] += phi[i, j, k, idx] * psi[i, j, k, idx]
            end
        end
    end

    return phi, psi, phipsi

end

function kernel_mode(
    M::I,
    umax::R,
    vmax::R,
    wmax::R,
    unum::I,
    vnum::I,
    wnum::I,
    alpha::R;
    quad_num = 64::I,
) where {I<:Integer,R<:Real}

    supp = sqrt(2.0) * 2.0 * max(umax, vmax, wmax) / (3.0 + sqrt(2.0))

    fre_vx = [i * π / umax for i = -unum÷2:unum÷2-1]
    fre_vy = [i * π / vmax for i = -vnum÷2:vnum÷2-1]
    fre_vz = [i * π / wmax for i = -wnum÷2:wnum÷2-1]

    abscissa, gweight = KitBase.lgwt(quad_num, 0, supp)

    phi = zeros(unum, vnum, wnum, M * (M - 1))
    psi = zeros(unum, vnum, wnum, M * (M - 1))
    phipsi = zeros(unum, vnum, wnum)
    for loop = 1:M-1
        theta = π / M * loop
        for loop2 = 1:M
            theta2 = π / M * loop2
            idx = (loop - 1) * M + loop2
            for k = 1:wnum, j = 1:vnum, i = 1:unum
                s =
                    fre_vx[i] * sin(theta) * cos(theta2) +
                    fre_vy[j] * sin(theta) * sin(theta2) +
                    fre_vz[k] * cos(theta)
                # phi
                int_temp = 0.0
                for id = 1:quad_num
                    int_temp +=
                        2.0 * gweight[id] * cos(s * abscissa[id]) * (abscissa[id]^alpha)
                end
                phi[i, j, k, idx] = int_temp * sin(theta)
                # psi
                s = fre_vx[i]^2 + fre_vy[j]^2 + fre_vz[k]^2 - s^2
                if s <= 0.0
                    psi[i, j, k, idx] = π * supp^2
                else
                    s = sqrt(s)
                    bel = supp * s
                    bessel = besselj(1, bel)
                    psi[i, j, k, idx] = 2.0 * π * supp * bessel / s
                end
                # phipsi
                phipsi[i, j, k] += phi[i, j, k, idx] * psi[i, j, k, idx]
            end
        end
    end

    return phi, psi, phipsi

end

function kernel_mode(
    M::TI,
    umax::TR,
    vmax::TR,
    du::TR,
    dv::TR,
    unum::TI,
    vnum::TI;
    quad_num = 64::TI,
) where {TI<:Integer,TR<:Real}

    supp = sqrt(2.0) * 2.0 * max(umax, vmax) / (3.0 + sqrt(2.0))

    fre_vx = range(-π / du, (unum ÷ 2 - 1) * 2.0 * π / unum / du, length = unum)
    fre_vy = range(-π / dv, (vnum ÷ 2 - 1) * 2.0 * π / vnum / dv, length = vnum)
    Fre_vx, Fre_vy = ndgrid(fre_vx, fre_vy)

    αp = zeros(unum, vnum, M)
    αq = zeros(unum, vnum, M)
    for lp = 1:M
        θ = lp * π / M

        s = Fre_vx * cos(θ + π / 2) + Fre_vy * sin(θ + π / 2)
        αq[:, :, lp] .= 2.0 * supp .* sin.(supp .* s) ./ (supp .* s) * (π / M)

        s = Fre_vx * cos(θ) + Fre_vy * sin(θ)
        αp[:, :, lp] .= 2.0 * supp .* sin.(supp .* s) ./ (supp .* s)
    end
    αp ./= π

    pq = αp .* αq
    α0 = zeros(unum, vnum)
    for j in axes(α0, 2), i in axes(α0, 1)
        α0[i, j] = sum(pq[i, j, :])
    end

    return αp, αq, α0

end

"""
    boltzmann_fft(
        f::X,
        Kn,
        M::I,
        ϕ::Y,
        ψ::Y,
        phipsi::Z,
    ) where {
        X<:AbstractArray{<:Real,3},
        Y<:AbstractArray{<:Real,4},
        Z<:AbstractArray{<:Real,3},
        I<:Integer,
    }

Calculate collision operator with FFT-based fast spectral method

"""
function boltzmann_fft(
    f::X,
    Kn,
    M::I,
    ϕ::Y,
    ψ::Y,
    phipsi::Z,
) where {
    X<:AbstractArray{<:Real,3},
    Y<:AbstractArray{<:Real,4},
    Z<:AbstractArray{<:Real,3},
    I<:Integer,
}

    f_spec = f .+ 0im
    bfft!(f_spec)
    f_spec ./= size(f, 1) * size(f, 2) * size(f, 3)
    f_spec .= fftshift(f_spec)

    #--- gain term ---#
    f_temp = zeros(axes(f_spec)) .+ 0im
    for i = 1:M*(M-1)
        fg1 = f_spec .* ϕ[:, :, :, i]
        fg2 = f_spec .* ψ[:, :, :, i]
        fg11 = fft(fg1)
        fg22 = fft(fg2)
        f_temp .+= fg11 .* fg22
    end

    #--- loss term ---#
    fl1 = f_spec .* phipsi
    fl2 = f_spec
    fl11 = fft(fl1)
    fl22 = fft(fl2)
    f_temp .-= fl11 .* fl22

    Q = @. 4.0 * π^2 / Kn / M^2 * real(f_temp)

    return Q

end

"""
    boltzmann_fft!(
        Q::T1,
        f::T2,
        Kn::TR,
        M::TI,
        ϕ::TY,
        ψ::TY,
        phipsi::TZ,
    ) where {
        T1<:AbstractArray{<:Real,3},
        T2<:AbstractArray{<:Real,3},
        TR<:Real,
        TI<:Integer,
        TY<:AbstractArray{<:Real,4},
        TZ<:AbstractArray{<:Real,3},
    }

Calculate collision operator with FFT-based fast spectral method

"""
function boltzmann_fft!(
    Q::T1,
    f::T2,
    Kn::TR,
    M::TI,
    ϕ::TY,
    ψ::TY,
    phipsi::TZ,
) where {
    T1<:AbstractArray{<:Real,3},
    T2<:AbstractArray{<:Real,3},
    TR<:Real,
    TI<:Integer,
    TY<:AbstractArray{<:Real,4},
    TZ<:AbstractArray{<:Real,3},
}

    f_spec = f .+ 0im
    bfft!(f_spec)
    f_spec ./= size(f, 1) * size(f, 2) * size(f, 3)
    f_spec .= fftshift(f_spec)

    # --- gain term ---#
    f_temp = zeros(axes(f_spec)) .+ 0im
    for i = 1:M*(M-1)
        fg1 = f_spec .* ϕ[:, :, :, i]
        fg2 = f_spec .* ψ[:, :, :, i]
        fg11 = fft(fg1)
        fg22 = fft(fg2)
        f_temp .+= fg11 .* fg22
    end

    # --- loss term ---#
    fl1 = f_spec .* phipsi
    fl2 = f_spec
    fl11 = fft(fl1)
    fl22 = fft(fl2)
    f_temp .-= fl11 .* fl22

    @. Q = 4.0 * π^2 / Kn / M^2 * real(f_temp)

end


text:
module GPCRAnalysis

using Downloads
using Statistics
using LinearAlgebra

using MIToS
using MIToS.Pfam
using MIToS.MSA
using MIToS.PDB

using MultivariateStats
using Distances
using Hungarian
using StaticArrays

export @res_str

export SequenceMapping
export species, uniprotX
export try_download_alphafold, download_alphafolds, getchain, findall_subseq
export filter_species!, filter_long!, chimerax_script
export project_sequences, columnwise_entropy, align, residue_centroid, residue_centroid_matrix, mapclosest, chargelocations, positive_locations, negative_locations
export StructAlign, residueindex, ismapped
export BWScheme, lookupbw
export aa_properties, aa_properties_zscored

include("naming_conventions.jl")
include("utils.jl")
include("msa.jl")
include("analyze.jl")
include("tmalign.jl")
include("ballesteros_weinstein.jl")
include("properties.jl")

end


text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#

### Data structures that are additionally available for media <: PartialPureSubstance -----------------------

"""
    state = setState_pT(medium, p,T)

Generate a state object for medium `medium::PureSubstance` for
pressure `p` [Pa] and temperature `T` [K].
"""
setState_pT( m::PureSubstance,p,T) = setState_pTX(m,p,T,fill(0.0,0))

"""
    state = setState_ph(medium, p,h)

Generate a state object for medium `medium::PureSubstance` for
pressure `p` [Pa] and specific enthalpy `h` [J/kg].
"""
setState_ph( m::PureSubstance,p,h) = setState_phX(m,p,h,fill(0.0,0))

"""
    state = setState_ps(medium, p,s)

Generate a state object for medium `medium::PureSubstance` for
pressure `p` [Pa] and specific entropy `s` [J/(kg*K)].
"""
setState_ps( m::PureSubstance,p,s) = setState_psX(m,p,s,fill(0.0,0))

"""
    state = setState_dT(medium, d,T)

Generate a state object for medium `medium::PureSubstance` for
density `d` [kg/m^3] and temperature `T` [K].
"""
setState_dT( m::PureSubstance,d,T) = setState_dTX(m,d,T,fill(0.0,0))

### Medium functions that are additionally available for media <: PartialPureSubstance -----------------------

"d = density_ph(medium,p,h) - return density in [kg/m^3] for `medium::PureSubstance` from pressure `p` [Pa] and specific enthalpy `h` [J/kg]."
density_ph(m::PureSubstance, p,h) = density(m, setState_ph(m,p,h))

"T = temperature_ph(medium,p,h) - return temperature in [K] for `medium::PureSubstance` from pressure `p` [Pa] and specific enthalpy `h` [J/kg]."
temperature_ph(m::PureSubstance, p,h) = temperature(m, setState_ph(m,p,h))

"p = pressure_dT(medium,d,T) - return pressure in [Pa] for `medium::PureSubstance` from density `c` in [kg/m^3] and temperature `T` [K]."
pressure_dT(m::PureSubstance, d,T) = pressure(m, setState_dT(m,d,T))

"h = specificEnthalpy_dT(medium,d,T) - return specific enthalpy in [J/kg] for `medium` from density `d` [kg/m^3] and and temperature `T` [K]."
specificEnthalpy_dT(m::PureSubstance, d,T) = specificEnthalpy(m, setState_dT(m,d,T))

### Derivatives needed for mass and energy balance ------------------------------------------------------

density_ph_der_1(m::PureSubstance, p, h) = 0.0
density_ph_der_2(m::PureSubstance, p, h) = undefinedFunction("density_ph_der_2")
density_ph_der_3(m::PureSubstance, p, h) = undefinedFunction("density_ph_der_3")


text:
list_test_str = split(
"NNCB

CH -> B
HH -> N
CB -> H
NH -> C
HB -> C
HC -> B
HN -> C
NN -> C
BH -> H
NC -> B
NB -> B
BN -> B
BB -> N
BC -> B
CC -> N
CN -> C",
'\n'
)

list_str = let input = "day14_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list, n)
    polymer = list[1]
    rules = Dict{String, String}()
    for l in list[3:end]
        a, b = split(l, " -> ")
        push!(rules, a => first(a) * b)
    end
    while n > 0
        oldpolymer = polymer
        polymer = ""
        for i in 1:length(oldpolymer)-1
            polymer *= rules[oldpolymer[i:i+1]]
        end
        polymer *= last(oldpolymer)
        n -= 1
    end
    res = extrema([count(c -> c == d, polymer) for d in unique(polymer)]) |> u -> last(u) - first(u)
    return polymer, res
end

testit("Test Part 1", part1(list_test_str, 1)[1], ==("NCNBCHB"))
testit("Test Part 1", part1(list_test_str, 2)[1], ==("NBCCNBBBCBHCB"))
testit("Test Part 1", part1(list_test_str, 3)[1], ==("NBBBCNCCNBBNBNBBCHBHHBCHB"))
testit("Test Part 1", part1(list_test_str, 4)[1], ==("NBBNBNBBCCNBCNCCNBBNBBNBBBNBBNBBCBHCBHHNHCBBCBHCB"))
testit("Test Part 1", part1(list_test_str, 10)[2], ==(1588))
testit("Challenge Part 1", part1(list_str, 10)[2], ==(3342))

# Part 2

function part2(list, n)
    polymer = list[1]
    rules = Dict{String, Tuple{String, String}}()
    cs = Dict{String, Int}()

    for l in list[3:end]
        a, b = split(l, " -> ")
        push!(rules, a => (first(a) * b, b * last(a)))
        push!(cs, a => length(collect(eachmatch(Regex(a), polymer))))
    end

    while n > 0
        for (k, v) in filter(kv -> last(kv) > 0, cs)
            cs[k] -= v
            cs[first(rules[k])] += v
            cs[last(rules[k])] += v
        end
        n -= 1
    end
    ks = Dict{Char, Int}()
    for (k, v) in filter(kv -> last(kv) > 0, cs)
        if first(k) in keys(ks)
            ks[first(k)] += v
        else
            push!(ks, first(k) => v)
        end
    end

    if last(polymer) in keys(ks)
        ks[last(polymer)] += 1
    else
        push!(ks, last(polymer) => 1)
    end
    m, n = extrema(values(ks))
    res = n - m
    return rules, cs, ks, res
end

testit("Test Part 2", part2(list_test_str, 10)[4], ==(1588))
testit("Test Part 2", part2(list_str, 10)[4], ==(3342))
testit("Test Part 2", part2(list_test_str, 40)[4], ==(2188189693529))
testit("Challenge Part 2", part2(list_str, 40)[4], ==(3776553567525))


text:
function move_to_position(init_position; init_velocity = [1.0, -0.01],
        goal_position=[2.0,0.0,0.0], T=5.0, stride_length_x=0.5,
        stride_length_y=0.4, a=1.0, b=1.0, delta_t=0.02)
    COM_pos_x = []
    COM_pos_y = []
    left_foot_pos_x = []
    left_foot_pos_y = []
    left_foot_pos_z = []
    right_foot_pos_x = []
    right_foot_pos_y = []
    right_foot_pos_z = []
    s_x = stride_length_x
    s_y = stride_length_y
    theta = goal_position[3]

    # Initialize the COM position and velocity
    COM_pos_0 =  init_position
    COM_v0 = init_velocity

    # Initialize the foot positions
    left_foot_pos = [COM_pos_0[1]+0.2, COM_pos_0[2]+0.1, 0 ]
    right_foot_pos = [COM_pos_0[1]+0.6, COM_pos_0[2]-0.5, 0]

    LIPM_model = LIP(dtt=delta_t, Ts=0.5)
    initializeModel!(LIPM_model,COM_pos_0, left_foot_pos, right_foot_pos)

    LIPM_model.support_leg = "left_leg"
    if LIPM_model.support_leg == "left_leg"
        support_foot_pos = LIPM_model.left_foot_pos
        LIPM_model.p_x = LIPM_model.left_foot_pos[1]
        LIPM_model.p_y = LIPM_model.left_foot_pos[2]
    else
        support_foot_pos = LIPM_model.right_foot_pos
        LIPM_model.p_x = LIPM_model.right_foot_pos[1]
        LIPM_model.p_y = LIPM_model.right_foot_pos[2]
    end

    LIPM_model.x_0 = LIPM_model.COM_pos[1] - support_foot_pos[1]
    LIPM_model.y_0 = LIPM_model.COM_pos[2] - support_foot_pos[2]
    LIPM_model.vx_0 = COM_v0[1]
    LIPM_model.vy_0 = COM_v0[2]

    step_num = 0
    total_time = T # seconds
    global_time = 0

    swing_data_len = Int(LIPM_model.T_sup/delta_t)
    swing_foot_pos = zeros((swing_data_len, 3))
    j = 1

    switch_index = swing_data_len

    for i = 0:(Int(total_time/delta_t)-1)
        global_time += delta_t
        step!(LIPM_model)
        if step_num >= 1
            if LIPM_model.support_leg == "left_leg"
                LIPM_model.right_foot_pos = [swing_foot_pos[j,1],
                             swing_foot_pos[j,2], swing_foot_pos[j,3]]
            else
                LIPM_model.left_foot_pos = [swing_foot_pos[j,1],
                            swing_foot_pos[j,2], swing_foot_pos[j,3]]

            end
            j += 1
        end

        # record data
        push!(COM_pos_x, LIPM_model.x_t + support_foot_pos[1])
        push!(COM_pos_y, LIPM_model.y_t + support_foot_pos[2])
        push!(left_foot_pos_x, LIPM_model.left_foot_pos[1])
        push!(left_foot_pos_y, LIPM_model.left_foot_pos[2])
        push!(left_foot_pos_z, LIPM_model.left_foot_pos[3])
        push!(right_foot_pos_x, LIPM_model.right_foot_pos[1])
        push!(right_foot_pos_y, LIPM_model.right_foot_pos[2])
        push!(right_foot_pos_z, LIPM_model.right_foot_pos[3])

        # switch the support leg
        if (i > 0) && (i % switch_index == 0)
            j = 1
            switchSupportLeg!(LIPM_model) # switch the support leg
            step_num += 1
            if at_goal([COM_pos_x[end], COM_pos_y[end]], goal_position)
                break
            end
            if LIPM_model.support_leg == "left_leg"
                support_foot_pos = LIPM_model.left_foot_pos
                LIPM_model.p_x = LIPM_model.left_foot_pos[1]
                LIPM_model.p_y = LIPM_model.left_foot_pos[2]
            else
                support_foot_pos = LIPM_model.right_foot_pos
                LIPM_model.p_x = LIPM_model.right_foot_pos[1]
                LIPM_model.p_y = LIPM_model.right_foot_pos[2]
            end
            x_0, vx_0, y_0, vy_0 = calculateXtVt!(LIPM_model, LIPM_model.T_sup)
            if LIPM_model.support_leg == "left_leg"
                x_0 = x_0 + LIPM_model.left_foot_pos[1]
                y_0 = y_0 + LIPM_model.left_foot_pos[2]
            else
                x_0 = x_0 + LIPM_model.right_foot_pos[1]
                y_0 = y_0 + LIPM_model.right_foot_pos[2]
            end

            calculateFootLocationForNextStep!(LIPM_model,s_x, s_y, a, b, theta,
                                        x_0, vx_0, y_0, vy_0)

            if LIPM_model.support_leg == "left_leg"
                right_foot_target_pos = [LIPM_model.p_x_star,
                                        LIPM_model.p_y_star, 0]
                swing_foot_pos[:,1] = collect(range(LIPM_model.right_foot_pos[1]
                        , stop=right_foot_target_pos[1], length=swing_data_len))
                swing_foot_pos[:,2] = collect(range(LIPM_model.right_foot_pos[2]
                        , stop=right_foot_target_pos[2], length=swing_data_len))
                swing_foot_pos[2:swing_data_len, 3] .= 0.1
            else
                left_foot_target_pos = [LIPM_model.p_x_star,
                                LIPM_model.p_y_star, 0]
                swing_foot_pos[:,1] = collect(range(LIPM_model.left_foot_pos[1],
                         stop=left_foot_target_pos[1], length=swing_data_len))
                swing_foot_pos[:,2] = collect(range(LIPM_model.left_foot_pos[2],
                 stop=left_foot_target_pos[2], length=swing_data_len))
                swing_foot_pos[2:swing_data_len, 3] .= 0.1
            end
        end
    end

    X = [COM_pos_x,
            COM_pos_y,
            left_foot_pos_x,
            left_foot_pos_y,
            left_foot_pos_z,
            right_foot_pos_x,
            right_foot_pos_y,
            right_foot_pos_z]
end


text:
module RecurrentTests

using Photon, Test
using Photon.Layers
using Photon.Losses
using Photon.Callbacks: SilentMeter

gethistoric(steps=100, features=10) = KorA(randn(features,steps,4))

function basic_layer(T)

    d = gethistoric(100)

    l = T(20)
    pred = l(d)
    @test size(pred) == (20,4)

    l = T(20,2)
    pred = l(d)
    @test size(pred) == (20,4)

    l = T(20;dropout=0.5, bidirectional=true)
    pred = l(d)
    @test size(pred) == (20,4)

    l = T(20;dropout=0.5, last_only=false)
    pred = l(d)
    @test size(pred) == (20,100,4)

    l = T(20; initw=zeros, initb=zeros, initf=zeros)
    pred = l(d)
    @test size(pred) == (20,4)
    @test sum(pred) == 0.0

end

function train_lstm_model()

    model = Sequential(
        LSTM(50, last_only=false),
        Dense(100, :relu),
        Dense(10)
    )

    w = Workout(model, MSELoss())

    X = [randn(10,20,8) for i in 1:10]
    Y = [randn(10,8)    for i in 1:10]
    data = zip(X,Y)
    for _ in 1:100
        train!(w, data;cb=SilentMeter())
        # GC.gc() # this line fixed some issue
    end
end

function lstm_model()

    model = Sequential(
        LSTM(20),
        Dense(100, :relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)

    pred = model(gethistoric(200))
    @test size(pred) == (10,4)

end

function lstm_model2()

    model = Sequential(
        LSTM(20, last_only=false),
        Dense(100, :relu),
        Dense(10)
    )

    pred = model(gethistoric(100,10))
    @test size(pred) == (10,4)

    # @test_throws DimensionMismatch model(gethistoric(200,20))
end

function gru_model()

    model = Sequential(
        GRU(20,2;dropout=0.5),
        Dense(100, :relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)
end

function rnn_model()

    model = Sequential(
        RNN(20,2;activation=:relu, dropout=0.5),
        Dense(100, :relu),
        Dense(10)
    )

    pred = model(gethistoric(100))
    @test size(pred) == (10,4)
end

@testset "Recurrent" begin
    train_lstm_model()
    basic_layer(LSTM)
    basic_layer(GRU)
    basic_layer(RNN)
    lstm_model()
    lstm_model2()
    gru_model()
    rnn_model()
end

end


text:
module zip

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh