text:
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using Evolutionary, Test, Random

function get_trains(p)
    fire = p.records[:fire]
    x, y = Float32[], Float32[]
    for time in eachindex(fire)
        for neuron_id in findall(fire[time])
            push!(x, time)
            push!(y, neuron_id)
        end
    end
    cellsa = Array{Union{Missing,Any}}(undef, 1, Int(findmax(y)[1]))
    nac = Int(findmax(y)[1])
    for (inx, cell_id) in enumerate(1:nac)
        cellsa[inx] = []
    end
    @inbounds for cell_id in unique(y)
        @inbounds for (time, cell) in collect(zip(x, y))
            if Int(cell_id) == cell
                append!(cellsa[Int(cell_id)], time)
            end
        end
    end
    cellsa

end

function Evolutionary.trace!(
    record::Dict{String,Any},
    objfun,
    state,
    population,
    method::GA,
    options,
)
    idx = sortperm(state.fitpop)
    record["fitpop"] = state.fitpop[:]
    record["pop"] = population[:]
end

function get_ranges(ranges)
    ###
    # Code appropriated from:
    # https://github.com/JuliaML/MLPlots.jl/blob/master/src/optional/onlineai.jl
    ###

    lower = Float32[]
    upper = Float32[]
    for (k, v) in ranges
        append!(lower, v[1])
        append!(upper, v[2])
    end
    lower, upper
end

function init_b(lower, upper)
    gene = [] # TODO: should be a typed list one day
    for (i, (l, u)) in enumerate(zip(lower, upper))
        p1 = rand(l:u, 1)
        append!(gene, p1)
    end
    gene
end

function initf(n)
    genesb = []  # TODO: should be a typed list one day
    for i = 1:n
        genes = init_b(lower, upper)
        append!(genesb, [genes])
    end
    genesb
end

function custom_raster(P1::Array, P2::Array)
    y0 = Int32[0]
    X = Float32[]
    Y = Float32[]
    N = 1
    y = 0
    for p in [P1, P2]
        y += 1
        append!(X, P1)
        append!(Y, 1 .+ sum(y0))
        push!(y0, 1)
        append!(X, P2)
        append!(Y, 1 .+ sum(y0))
        push!(y0, 1)
        N += 1
    end
    plt = scatter(
        X,
        Y,
        w = 50,
        m = (1, :black),
        leg = :none,
        marker = :vline,
        xaxis = ("t", (0, Inf)),
        yaxis = ("neuron",),
    )
    y0 = y0[2:end-1]
    !isempty(y0) && hline!(plt, cumsum(y0), linecolor = :red)
    return plt
end

function get_vm(p, simulation_duration)
    vm = p.records[:v]
    vm = [i[1] for i in vm]
    vm = signal(vm, length(vm) / simulation_duration)
    vm
end

function checkmodel(param, cell_type, ngt_spikes)
    if cell_type == "IZHI"
        pp = SNN.IZParameter(; a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(; N = 1, param = pp)
    end
    if cell_type == "ADEXP"
        adparam = SNN.ADEXParameter(;
            a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10],
        )
        E = SNN.AD(; N = 1, param = adparam)
    end

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    current = current_search(cell_type, param, ngt_spikes)
    E.I = [current * nA]

    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])

    ###
    # TODO buid t into neuron models
    # SNN.monitor(E, [:t])
    ###
    simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 343ms
    SNN.sim!(
        [E];
        dt = 1 * ms,
        delay = ALLEN_DELAY,
        stimulus_duration = ALLEN_DURATION,
        simulation_duration = simulation_duration,
    )
    spikes = get_spikes(E)
    spikes = [s / 1000.0 for s in spikes]
    vm = get_vm(E, simulation_duration)
    (vm, spikes)
end

#=
function Evolutionary.value!(
    ::Val{:multiproc},
    fitness,
    objfun,
    population::AbstractVector{IT},
) where {IT}
    fitness = SharedArrays.SharedArray{Float32}(fitness)
    @time @sync @distributed for i = 1:length(population)
        fitness[i] = value(objfun, population[i])
    end
    fitness
end
=#
function Evolutionary.value!(
    ::Val{:serial},
    fitness,
    objfun,
    population::AbstractVector{IT},
) where {IT}
    Threads.@threads for i = 1:length(population)
        fitness[i] = value(objfun, population[i])
    end
end
function get_data()
    file = "../JLD/ground_truth.jld"
    if file
        vmgtv = load(file, "vmgtv")
        ngt_spikes = load(file, "ngt_spikes")
        gt_spikes = load(file, "gt_spikes")

        ground_spikes = gt_spikes
        ngt_spikes = size(gt_spikes)[1]
        vmgtt = load(file, "vmgtt")

        plot(plot(vmgtv, vmgtt, w = 1))

    else
        # if the file doesn't exist complete the expensive operation of making the file.
        # using Python
        py"""
        from neo import AnalogSignal
        from neuronunit.allenapi import make_allen_tests_from_id

        """

        py"""
        specimen_id = (
            325479788,
            324257146,
            476053392,
            623893177,
            623960880,
            482493761,
            471819401
        )
        specimen_id = specimen_id[1]
        target_num_spikes=7
        sweep_numbers, data_set, sweeps = make_allen_tests_from_id.allen_id_to_sweeps(specimen_id)
        (vmm,stimulus,sn,spike_times) = make_allen_tests_from_id.get_model_parts_sweep_from_spk_cnt(
            target_num_spikes, data_set, sweep_numbers, specimen_id
        )
        """
        gt_spikes = py"spike_times"
        ground_spikes = gt_spikes
        ngt_spikes = size(gt_spikes)[1]
        vmgtv = py"vmm.magnitude"
        vmgtt = py"vmm.times"
        julia_version_vm = signal(vmgtv, length(vmgtt) / last(vmgtt))
        plot(plot(vmgtv, vmgtt, w = 1))

        save(
            "ground_truth.jld",
            "vmgtv",
            vmgtv,
            "vmgtt",
            vmgtt,
            "ngt_spikes",
            ngt_spikes,
            "gt_spikes",
            gt_spikes,
            "julia_version_vm",
            julia_version_vm,
        )
    end

    vmgtv = load("ground_truth.jld", "vmgtv")
    ngt_spikes = load("ground_truth.jld", "ngt_spikes")
    ngt_spikes = size(gt_spikes)[1]
    ground_spikes = load("ground_truth.jld", "gt_spikes")

    vmgtt = load("ground_truth.jld", "vmgtt")
    julia_version_vm = load("ground_truth.jld", "julia_version_vm")
    return (vmgtv, vmgtt, ngt_spikes, ground_spikes, julia_version_vm)
end
function get_izhi_ranges()
    ranges_izhi = DataStructures.OrderedDict{Char,Float32}()
    ranges_izhi =
        ("a" => (0.002, 0.3), "b" => (0.02, 0.36), "c" => (-75, -35), "d" => (0.005, 16))#,"I"=>[100,9000])
    lower, upper = get_ranges(ranges_izhi)
    return lower, upper
end

function get_adexp_ranges()
    ranges_adexp = DataStructures.OrderedDict{String,Tuple{Float32,Float32}}()
    ranges_adexp[:"a"] = (2.0, 10)
    ranges_adexp[:"b"] = (5.0, 10)
    ranges_adexp[:"cm"] = (700.0, 983.5)
    ranges_adexp[:"v0"] = (-70, -55)
    ranges_adexp[:"τ_m"] = (10.0, 42.78345)
    ranges_adexp[:"τ_w"] = (300.0, 454.0)  # Tau_w 0, means very low adaption
    ranges_adexp[:"θ"] = (-45.0, -10)
    ranges_adexp[:"delta_T"] = (1.0, 5.0)
    ranges_adexp[:"v_reset"] = (-70.0, -15.0)
    ranges_adexp[:"spike_delta"] = (1.25, 20.0)
    lower, upper = get_ranges(ranges_adexp)
    return lower, upper
end

function vecplot(p, sym)
    v = SNN.getrecord(p, sym)
    y = hcat(v...)'
    x = 1:length(v)
    plot(x, y, leg = :none, xaxis = ("t", extrema(x)), yaxis = (string(sym), extrema(y)))
end

function vecplot(P::Array, sym)
    plts = [vecplot(p, sym) for p in P]
    N = length(plts)
    plot(plts..., size = (600, 400N), layout = (N, 1))
end

function spike_train_difference(spkd_ground, spkd_found)
    if length(spkd_found) == 0
        return sum(ones(Int(length(spkd_ground))))
    end
    if length(spkd_found) > 0
        maxi0 = findmax(spkd_ground[:, :])[1]
        maxi1 = findmax(spkd_found[:, :])[1]
        maxi1 = maxi1[1]
        maxi0 = maxi0[1]
        mini = findmin([maxi0, maxi1])[1]
        l0 = length(spkd_ground)
        l1 = length(spkd_found)
        maxi = findmin([l0, l1])[1]
        spkd = zeros(Int(maxi))
        @inbounds for (_, i) in zip(spkd, eachindex(spkd))
            if !isempty(spkd_ground[i]) && !isempty(spkd_found[i])

                maxt1 = findmax(spkd_ground[i])[1]
                maxt2 = findmax(spkd_found[i])[1]
                maxt = findmax([maxt1, maxt2])[1]

                if maxt1 > 0.0 && maxt2 > 0.0
                    t, S = SpikeSynchrony.SPIKE_distance_profile(
                        unique(sort(spkd_ground[i])),
                        unique(sort(spkd_found[i]));
                        t0 = 0.0,
                        tf = maxt,
                    )
                    spkd[i] = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1]) # == SPIKE_distance(y1, y2)

                end
            else
                spkd[i] = 0.0
            end
        end
        scatter([i for i = 1:mini], spkd) |> display
        error = sum(spkd)
    end
end

function initd()
    population = initd(10)
    garray = zeros((length(population)[1], length(population[1])))
    for (i, p) in enumerate(population)
        garray[i, :] = p
    end
    garray[1, :]
end

function initd(n)
    genesb = []
    lower = Float32[0.0 0.0 0.0 0.0]
    upper = Float32[1.0 1.0 1.0 1.0]
    lower = vec(lower)
    upper = vec(upper)

    for i = 1:n
        genes = initd(lower, upper)
        append!(genesb, [genes])
    end
    genesb
end

function initd(lower, upper)
    gene = []

    for (i, (l, u)) in enumerate(zip(lower, upper))
        p1 = rand(l:u, 1)
        append!(gene, p1)
    end
    gene
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const APIVersion = (;[
    Symbol("APIVERSION_UNSPECIFIED") => Int32(0),
    Symbol("V1") => Int32(1),
]...)

export APIVersion

ext: jl


text:
using WebIO
using Test

@testset "Build JavaScript" begin
    if haskey(ENV, "WEBIO_SKIP_BUNDLEJS")
        @info "Skipping building JavaScript."
        return
    end
    WebIO.bundlejs()
end


text:
using SurfaceGeometry

isdefined(:N) || (N=100)
isdefined(:PasiveStabilisation) || (PasiveStabilisation=false)
isdefined(:NormalVelocity) || (NormalVelocity=false)
isdefined(:ActiveStabilisation) || (ActiveStabilisation=false)
isdefined(:ViewerActive) || (ViewerActive=false)

using JLD
data = load("sphere.jld")
points = data["points"]
faces = data["faces"]

function velocity(t,pos)
    x,y,z = pos

    x = x*0.15 + 0.35
    y = y*0.15 + 0.35
    z = z*0.15 + 0.35

    u = 2*sin(pi*x)^2 * sin(2*pi*y) * sin(2*pi*z) * sin(2/3*pi*t)
    v = - sin(2*pi*x) * sin(pi*y)^2 * sin(2*pi*z) * sin(2/3*pi*t)
    w = - sin(2*pi*x) * sin(2*pi*y) * sin(pi*z)^2 * sin(2/3*pi*t)

    [u,v,w] / 0.5 # / 0.15
end

### In reality this is very complicated function
function velocity{T<:Real}(t::Real,points::Array{T,2},faces::Array{Int,2})

    ### Some magic here on how to calculate the velocity
    v = Array(Float64,size(points)...)
    for xkey in 1:size(points,2)
        v[:,xkey] = velocity(t,points[:,xkey])
    end

    return v
end

memory = []

h = 3/N

t1 = 0
p1 = points
push!(memory,(t1,p1,faces))

### MidPoint method
v1 = velocity(t1,p1,faces)
v1 = velocity(t1 + h/2,p1 + h/2*v1,faces)
t2,p2 = t1 + h, p1 + h*v1

push!(memory,(t2,p2,faces))

for i in 1:50
    println("It is step $i")

    v2 = velocity(t2,p2,faces)

    ### Some stabilisation magic here
    # n = zeros(size(points)...)
    # NormalVectors!(n,points,faces,i->FaceVRing(i,faces))

    # for xkey in 1:size(points,2)
    #     #P = enye(3) - n[:,xkey]*n[:,xkey]'
    #     v2[:,xkey] = n[:,xkey]*dot(n[:,xkey],v2[:,xkey])
    # end

    if PasiveStabilisation==true
        res = stabilise(p2,faces,v2;vp=v1)
        println("Finit $(res.Finit) \t Fres $(res.Fres)")
        v2 = res.vres
    end

    ### Adams Bachforth oneliner
    (t1,p1,v1), (t2,p2) = (t2,p2,v2), (t2+h,p2 + h/2/(t2-t1)*( (2*(t2-t1) + h)*v2 - h*v1))

    push!(memory,(t2,p2,faces))
end

### For viewing results

using Escher
if !ViewerActive
    include(Pkg.dir("Escher", "src", "cli", "serve.jl"))
    @spawn escher_serve(5555,Pkg.dir("SurfaceGeometry","examples","viewers"))
    ViewerActive = true
end


text:
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This example demonstrates how to use Knitro to solve the following
# simple mathematical program with equilibrium/complementarity
# constraints(MPEC/MPCC).
#
#  min  (x0 - 5)^2 +(2 x1 + 1)^2
#  s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2
#        x2 complements(3 x0 - x1 - 3)
#        x3 complements(-x0 + 0.5 x1 + 4)
#        x4 complements(-x0 - x1 + 7)
#        x0, x1, x2, x3, x4 >= 0
#
# The complementarity constraints must be converted so that one
# nonnegative variable complements another nonnegative variable.
#
#  min  (x0 - 5)^2 +(2 x1 + 1)^2
#  s.t.  -1.5 x0 + 2 x1 + x2 - 0.5 x3 + x4 = 2  (c0)
#        3 x0 - x1 - 3 - x5 = 0                 (c1)
#        -x0 + 0.5 x1 + 4 - x6 = 0              (c2)
#        -x0 - x1 + 7 - x7 = 0                  (c3)
#        x2 complements x5
#        x3 complements x6
#        x4 complements x7
#        x0, x1, x2, x3, x4, x5, x6, x7 >= 0
#
# The solution is(1, 0, 3.5, 0, 0, 0, 3, 6), with objective value 17.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

using JuMP, MathOptInterface
using KNITRO
using Test

const MOI = MathOptInterface
b = [-3.0, 4.0, 7.0]
M = [3.0   -1.0   0.0   0.0   0.0;
     -1.0   0.5   0.0   0.0   0.0;
     -1.0  -1.0   0.0   0.0   0.0]

model = JuMP.Model(with_optimizer(KNITRO.Optimizer))

@variable(model, x[1:5] >= 0)
@constraint(model, -1.5 * x[1] + 2.0 * x[2] + x[3] - 0.5 * x[4] + x[5] == 2.0)

@constraint(model, [M*x + b; x[3:5]] in MOI.Complements(3))
@objective(model, Min, (x[1] - 5.0)^2 + (2.0 * x[2] + 1.0)^2)

JuMP.optimize!(model)

@test JuMP.has_values(model)
@test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
@test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
@test JuMP.value.(x) ≈ [1.0, 0.0, 3.5, 0.0, 0.0] atol=1e-3


text:
#!/bin/sh
#= La ligne shell suivante est en commentaire multiligne pour julia
exec julia --color=yes --startup-file=no --depwarn=no -- "$0" "$@"
=#

# IDÉE TODO :
# - faire macros @p (pour print intelligent) et @s (pour to_string générique)

# Le realpath est nécessaire si l'un des répertoire parents est un lien
# symbolique :
global const APPDIR = dirname(dirname(realpath(@__FILE__())))
@show APPDIR
# utile pour les includes mais fonctionnement douteux !?
push!(LOAD_PATH, "$APPDIR/src")
println("LOAD_PATH=\n    ", join(LOAD_PATH, "\n    "))

include("$APPDIR/src/utils/time.jl") # pour la macro @ms()

using Printf
using Random
using Dates

if basename(@__FILE__) == basename(PROGRAM_FILE)
    # Mode d'appel normal : on exécute le programme "bin/xxx.jl"
    @ms include("$APPDIR/src/main.jl")
    main()
else
    # Mode interactif de l'exécutable seqata.lj
    #
    # Exécution par :
    #    julia -iL ./bin/seqata.jl
    #    julia -i --color=yes -L  ./bin/seqata.jl
    #    include("test/test-21-mutation_move2x_solution.jl")
    #
    # - le main n'est pas appelé
    # - on inclue tous les packages et les includes possibles
    # - l'analyse des arguments est faite
    # - on charge le fichier d'utilitatire dédié au mode interactif
    #

    # @ms using Revise # pour rechargement dynamique des fichiers/fonctions modifiés
    # MAIS NE FONCTIONNE QUE POUR DES FONCTIONS INTÉGRÉES DANS UN PACKAGE
    # (donc ne fonctionne par avec seqata pour l'instant)
    # using Revise # pour rechargement dynamique des fichiers/fonctions modifiés

    @ms include("$APPDIR/src/args.jl");
    @ms include("$APPDIR/src/utils/log.jl")
    @ms include("$APPDIR/src/utils/console.jl")
    @ms include("$APPDIR/src/plane.jl")
    @ms include("$APPDIR/src/instance.jl")
    @ms include("$APPDIR/src/processing/instance_generators.jl")
    @ms include("$APPDIR/src/processing/instance_read_alp.jl")
    @ms include("$APPDIR/src/utils/array.jl")
    @ms include("$APPDIR/src/utils/file.jl")
    @ms include("$APPDIR/src/solution.jl")
    @ms include("$APPDIR/src/processing/solution_read.jl")

    @ms using JuMP
    @ms using CPLEX
    @ms using Glob
    @ms using Dates
    @ms using Test

    @ms include("$APPDIR/src/utils/model.jl")
    @ms include("$APPDIR/src/solvers/earliest_timing.jl")
    @ms include("$APPDIR/src/solvers/lp_timing.jl")
    @ms include("$APPDIR/src/solvers/mip.jl")

    @ms include("$APPDIR/src/solvers/annealing.jl")
    @ms include("$APPDIR/src/solvers/descent.jl")
    @ms include("$APPDIR/src/solvers/explore.jl")

    # On impose des arguments par défaut
    # ARGS = ["test"]
    # @ms args = Args.parse_commandline(["test"])
    @ms args = Args.parse_commandline()
    lg1() && Args.show_args()

    @ms include("../src/interactive.jl")

    println()
    println("Début de mode interactif de seqata.jl ($(ms()))s")
    println()

end


text:
using POMDPToolbox
using POMDPModels

let
    pomdp = BabyPOMDP()
    fwc = FeedWhenCrying()
    rnd = solve(RandomSolver(MersenneTwister(7)), pomdp)

    q = []
    push!(q, Sim(pomdp, fwc, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"feed when crying")))
    push!(q, Sim(pomdp, fwc, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"feed when crying")))
    push!(q, Sim(pomdp, rnd, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"random")))

    println("There should be a warning here:")
    run_parallel(q)

    procs = addprocs(1)
    @everywhere using POMDPToolbox
    @everywhere using POMDPModels
    println("There should not be a warning here:")
    @show run_parallel(q) do sim, hist
        return [:steps=>n_steps(hist), :reward=>discounted_reward(hist)]
    end

    @show data = run_parallel(q)
    @test data[1, :reward] == data[2, :reward]
    @test var(data[:reward][1:2]) == 0.0
    rmprocs(procs)

    mdp = GridWorld()
    q = []
    push!(q, Sim(mdp, RandomPolicy(mdp), max_steps=100))
    run(q)
end

# example from readme
let
    using POMDPToolbox
    using POMDPModels

    pomdp = BabyPOMDP()
    fwc = FeedWhenCrying()
    rnd = solve(RandomSolver(MersenneTwister(7)), pomdp)

    q = [] # vector of the simulations to be run
    push!(q, Sim(pomdp, fwc, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"feed when crying")))
    push!(q, Sim(pomdp, rnd, max_steps=32, rng=MersenneTwister(4), metadata=Dict(:policy=>"random")))

    # this creates two simulations, one with the feed-when-crying policy and one with a random policy

    data = run_parallel(q)

    # by default, the dataframe output contains the reward and the contents of `metadata`
    @show data
    # data = 2×2 DataFrames.DataFrame
    # │ Row │ policy             │ reward   │
    # ├─────┼────────────────────┼──────────┤
    # │ 1   │ "feed when crying" │ -4.5874  │
    # │ 2   │ "random"           │ -27.4139 │

    # to perform additional analysis on each of the simulations one can define a processing function with the `do` syntax:
    data2 = run_parallel(q, progress=false) do sim, hist
        println("finished a simulation - final state was $(last(state_hist(hist)))")
        return [:steps=>n_steps(hist), :reward=>discounted_reward(hist)]
    end

    @show data2
    # 2×3 DataFrames.DataFrame
    # │ Row │ policy             │ reward   │ steps │
    # ├─────┼────────────────────┼──────────┼───────┤
    # │ 1   │ "feed when crying" │ -18.2874 │ 32.0  │
    # │ 2   │ "random"           │ -17.7054 │ 32.0  │
end


text:
struct GroupReserve <: AbstractReservesFormulation end
############################### Reserve Variables` #########################################
"""
This function checks if the variables for reserves were created
"""
function check_activeservice_variables(
    psi_container::PSIContainer,
    contributing_services::Vector{<:PSY.Service},
)
    for service in contributing_services
        get_variable(psi_container, PSY.get_name(service), typeof(service))
    end
    return
end

################################## Reserve Requirement Constraint ##########################
"""
This function creates the requirement constraint that will be attained by the apropriate services
"""
function service_requirement_constraint!(
    psi_container::PSIContainer,
    service::SR,
    ::ServiceModel{SR, GroupReserve},
    contributing_services::Vector{<:PSY.Service},
) where {SR <: PSY.StaticReserveGroup}
    parameters = model_has_parameters(psi_container)
    initial_time = model_initial_time(psi_container)
    @debug initial_time
    time_steps = model_time_steps(psi_container)
    name = PSY.get_name(service)
    constraint = get_constraint(psi_container, make_constraint_name(REQUIREMENT, SR))
    use_slacks = get_services_slack_variables(psi_container.settings)
    reserve_variables = [
        get_variable(psi_container, PSY.get_name(r), typeof(r))
        for r in contributing_services
    ]

    requirement = PSY.get_requirement(service)
    for t in time_steps
        resource_expression = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}()
        for reserve_variable in reserve_variables
            JuMP.add_to_expression!(resource_expression, sum(reserve_variable[:, t]))
        end
        if use_slacks
            resource_expression += slack_vars[t]
        end
        constraint[name, t] =
            JuMP.@constraint(psi_container.JuMPmodel, resource_expression >= requirement)
    end

    return
end


text:
#!/usr/bin/env julia --project

include("simple_box.jl")
ClimateMachine.init(parse_clargs = true)

# Float type
const FT = Float64

# simulation time
const timestart = FT(0)      # s
const timestep = FT(240)     # s
const timeend = FT(86400) # s
timespan = (timestart, timeend)

# DG polynomial order
const N = Int(4)

# Domain resolution
const Nˣ = Int(20)
const Nʸ = Int(20)
const Nᶻ = Int(20)
resolution = (N, Nˣ, Nʸ, Nᶻ)

# Domain size
const Lˣ = 4e6    # m
const Lʸ = 4e6    # m
const H = 1000   # m
dimensions = (Lˣ, Lʸ, H)

BC = (
    OceanBC(Impenetrable(NoSlip()), Insulating()),
    OceanBC(Impenetrable(NoSlip()), Insulating()),
    OceanBC(Penetrable(KinematicStress()), TemperatureFlux()),
)

run_simple_box(
    "ocean_gyre",
    resolution,
    dimensions,
    timespan,
    OceanGyre,
    imex = false,
    Δt = timestep,
    BC = BC,
)


text:
module AutoCorrelationResampling
using Optim
using Base.Iterators

AV  = AbstractVector
AF  = AbstractFloat
VoF = AV{<:AF}

include("filter-coeffs.jl")
include("ac-resample.jl")

export filter_coeffs, ac_resample

end # module


text:
using Discord
include("base.jl")

function cancel(client::Client,m::Message,boss_number::Int)
    global  management_dict

    channel_id_str = string(m.channel_id)

    if !haskey(management_dict,channel_id_str)
        return reply(client,m,"初期設定が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if !haskey(management_dict[channel_id_str],"remain_attack_channel")
        return reply(client,m,"残凸把握板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_channel")
        return reply(client,m,"予約確認板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if boss_number > 5 || boss_number < 1
        return reply(client,m,"不正な入力です．[ボス番号:$(boss_number)は1から5である必要があります]")
    end

    boss_number = string(boss_number)
    reserve_list =  management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"]

    dict_number = 0

    #予約ずみかどうかをチェックする
    for i = 1:length(reserve_list)
        if reserve_list[i]["id"] == m.author.id
            dict_number = i
            break
        end
    end

    #予約していなかっtら怒る
    if dict_number == 0
        return reply(client,m,"凸予約がされていません.")
    end

    deleteat!(management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"],dict_number)
    management_dict[channel_id_str]["reserve"][boss_number]["plan_remain_hp"] = calc_plan_remain_hp(management_dict[channel_id_str]["reserve"][boss_number])

    output_message_content = create_output_message(management_dict[channel_id_str])
    edit_message(client,management_dict[channel_id_str]["output_channel"],management_dict[channel_id_str]["output_message"],content=output_message_content)

    save_management_dict(management_dict)
    create(client, Reaction, m, '👍')

end

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
module Equations
export Equation, EquationChain, Expression, Component, NonAbelian, Operator, Ex, Term, Factor, ╱,Div, Sqrt, Pow,Exp,Log,Cos,Sin,Cosh,Sinh, Der,DerOp, Vec,Cross,Dot,Norm,Mat, Named,Oneable, U,Physical, AbstractTensor,Ten,Up,Alt,Tensor,BraKet,Delta,⊗,TensorProduct, ∧,Wedge,D,Form,Trace,Commutator,Transpose
export equation, solve, expression, ≖, evaluate, simplify, simplify!, componify, componify!, addparse, has, sumnum, sumsym, matches, getarg, findpows, indin, indsin, replace, findsyms, quadratic, complexity, terms, getargs, randeval, expandindices,pushall!,sumconv,duplicates,ps, @equ,@equs
export relations

include("equ.jl")

end # This may seem like the end but is actually a beginning.


text:
immutable System2Intercept{T} <: Model
    β₀::T
    a::T
    ω::T
    δ₁::T
    δ₂::T
end

Base.show(io::IO, ::Type{System2Intercept}) = print(io, "Intercept System-2")

function Base.show(io::IO, m::System2Intercept)
    @printf(
        io,
        "Intercept System-2(β₀ = %.3f, a = %.3f, ω = %.3f, δ₁ = %.3f, δ₂ = %.3f)",
        m.β₀,
        m.a,
        m.ω,
        m.δ₁,
        m.δ₂,
    )
end

function initial_parameters(::Type{System2Intercept})
    return Float64[0.0, log(1.0), logit(0.5), logit(0.25), logit(0.75)]
end

canonical_parameters(m::System2Intercept) = (m.β₀, m.a, m.ω, m.δ₁, m.δ₂)

function System2Intercept(parameters::Vector)
    β₀ = parameters[1]
    a = exp(parameters[2])
    ω = invlogit(parameters[3])
    δ₁ = invlogit(parameters[4])
    δ₂ = invlogit(parameters[5])
    return System2Intercept(β₀, a, ω, δ₁, δ₂)
end

@inline function predict(
    m::System2Intercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, ω, δ₁, δ₂ = m.β₀, m.a, m.ω, m.δ₁, m.δ₂
    u1 = x1 * (ω * δ₁^t1 + (1 - ω) * δ₂^t1)
    u2 = x2 * (ω * δ₁^t2 + (1 - ω) * δ₂^t2)
    z = β₀ + a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::System2Intercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, ω, δ₁, δ₂ = m.β₀, m.a, m.ω, m.δ₁, m.δ₂
    gr = (
        1.0,
        (x2 * (ω * δ₁^t2 + (1 - ω) * δ₂^t2) - x1 * (ω * δ₁^t1 + (1 - ω) * δ₂^t1)) * a,
        a * (x2 * (δ₁^t2 - δ₂^t2) - x1 * (δ₁^t1 - δ₂^t1)) * ω * (1 - ω),
        a * (x2 * ω * t2 * δ₁^(t2 - 1) - x1 * ω * t1 * δ₁^(t1 - 1)) * δ₁ * (1 - δ₁),
        a * (x2 * (1 - ω) * t2 * δ₂^(t2 - 1) - x1 * (1 - ω) * t1 * δ₂^(t1 - 1)) * δ₂ * (1 - δ₂),
    )
    return gr
end


text:
@testset "Trie" begin
    @testset "Core Functionality" begin
        t = Trie{Int}()
        t["amy"] = 56
        t["ann"] = 15
        t["emma"] = 30
        t["rob"] = 27
        t["roger"] = 52
        t["kevin"] = Int8(11)

        @test haskey(t, "roger")
        @test get(t,"rob",nothing) == 27
        @test sort(keys(t)) == ["amy", "ann", "emma", "kevin", "rob", "roger"]
        @test t["rob"] == 27
        @test sort(keys_with_prefix(t,"ro")) == ["rob", "roger"]
    end

    @testset "Constructors" begin
        ks = ["amy", "ann", "emma", "rob", "roger"]
        vs = [56, 15, 30, 27, 52]
        kvs = collect(zip(ks, vs))
        @test isa(Trie(ks, vs), Trie{Int})
        @test isa(Trie(kvs), Trie{Int})
        @test isa(Trie(Dict(kvs)), Trie{Int})
        @test isa(Trie(ks), Trie{Nothing})
    end

    @testset "partial_path iterator" begin
        t = Trie{Int}()
        t["rob"] = 27
        t["roger"] = 52
        t["kevin"] = Int8(11)
        t0 = t
        t1 = t0.children['r']
        t2 = t1.children['o']
        t3 = t2.children['b']
        @test collect(partial_path(t, "b")) == [t0]
        @test collect(partial_path(t, "rob")) == [t0, t1, t2, t3]
        @test collect(partial_path(t, "robb")) == [t0, t1, t2, t3]
        @test collect(partial_path(t, "ro")) == [t0, t1, t2]
        @test collect(partial_path(t, "roa")) == [t0, t1, t2]
    end
end # @testset Trie

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

prob = prob_ode_lotkavoltera
prob = remake(prob, tspan=(0.0, 10.0))
prob = ODEFilters.remake_prob_with_jac(prob)

@testset "Condition numbers of A,Q" begin
    h = 0.1*rand()
    σ = rand()

    d, q = 2, 3

    A!, Q! = ODEFilters.vanilla_ibm(d, q)
    Ah = diagm(0 => ones(d*(q+1)))
    Qh = zeros(d*(q+1), d*(q+1))
    A!(Ah, h)
    Q!(Qh, h, σ^2)

    A_p, Q_p = ODEFilters.ibm(d, q)
    Ah_p = A_p
    Qh_p = Q_p * σ^2

    # First test that they're both equivalent
    P = ODEFilters.preconditioner(d, q)
    P, PI = P(h), inv(P(h))
    @test Qh_p ≈ P * Qh * P'
    @test Ah_p ≈ P * Ah * PI

    # Check that the preconditioning actually helps
    @info "Condition numbers" cond(Qh) cond(Matrix(Qh_p)) cond(Ah) cond(Ah_p)
    @test cond(Qh) > cond(Matrix(Qh_p))
    @test cond(Qh) > cond(Matrix(Qh_p))^2
end


text:
module MedEye3d
import Logging

export  ForDisplayStructs
# export  ForDisplayStructs.TextureSpec
export  SegmentationDisplay

# export  DataStructs.ThreeDimRawDat
# export  DataStructs.DataToScrollDims
# export  DataStructs.FullScrollableDat
# export  ForDisplayStructs.KeyboardStruct
# export  ForDisplayStructs.MouseStruct
# export  ForDisplayStructs.ActorWithOpenGlObjects
export  OpenGLDisplayUtils
# export  DisplayWords.textLinesFromStrings
export  StructsManag
# export  DisplayWords.textLinesFromStrings
# export  StructsManag.getThreeDims

include(joinpath( "display","GLFW","startModules","ModernGlUtil.jl"))

include(joinpath( "structs","BasicStructs.jl"))
include(joinpath( "structs","DataStructs.jl"))
include(joinpath( "structs","ForDisplayStructs.jl"))

include(joinpath( "display","GLFW","DispUtils","StructsManag.jl"))

include(joinpath( "display","GLFW","startModules","PrepareWindowHelpers.jl"))
include(joinpath( "display","GLFW","shadersEtc","CustomFragShad.jl"))

include(joinpath( "display","GLFW","DispUtils","OpenGLDisplayUtils.jl"))
include(joinpath( "display","GLFW","shadersEtc","ShadersAndVerticies.jl"))
include(joinpath( "display","GLFW","shadersEtc","ShadersAndVerticiesForText.jl"))
include(joinpath( "display","GLFW","shadersEtc","Uniforms.jl"))

include(joinpath( "display","GLFW","textRender","DisplayWords.jl"))

include(joinpath( "display","GLFW","DispUtils","TextureManag.jl") )
include(joinpath( "display","GLFW","startModules","PrepareWindow.jl"))

include(joinpath( "display","reactingToMouseKeyboard","ReactToScroll.jl") )
include(joinpath( "display","reactingToMouseKeyboard","ReactOnMouseClickAndDrag.jl") )

include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","KeyboardMouseHelper.jl") )

include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","MaskDiffrence.jl") )
include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","KeyboardVisibility.jl") )
include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","OtherKeyboardActions.jl") )
include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","WindowControll.jl") )
include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","ChangePlane.jl") )
include(joinpath( "display","reactingToMouseKeyboard","reactToKeyboard","reactToKeyboard.jl") )

include(joinpath( "display","reactingToMouseKeyboard","ReactingToInput.jl") )
include(joinpath( "display","GLFW","SegmentationDisplay.jl"))

greet() = print("Hello World!")

end # module


text:
import IntervalRootFinding, LinearAlgebra
using IntervalRootFinding: (..), (×), IntervalBox, interval
export fixedpoints, .., ×, IntervalBox, interval

"""
    fixedpoints(ds::DynamicalSystem, box, p = ds.p; kwargs...) → fp, eigs, stable
Return all fixed points `fp` of the given `ds`
that exist within the state space subset `box` for parameter configuration `p`.
Fixed points are returned as a [`Dataset`](@ref).
For convenience, a vector of the Jacobian eigenvalues of each fixed point, and whether
the fixed points are stable or not, are also returned.

`box` is an appropriate `IntervalBox` from IntervalRootFinding.jl.
E.g. for a 3D system it would be something like
```julia
v, z = -5..5, -2..2   # 1D intervals, can use `interval(-5, 5)` instead
box = v × v × z       # `\\times = ×`, or use `IntervalBox(v, v, z)` instead
```

Internally IntervalRootFinding.jl is used and as a result we are guaranteed to find all
fixed points that exist in `box`, regardless of stability. Since IntervalRootFinding.jl
returns an interval containing a unique fixed point, we return the midpoint of the
interval as the actual fixed point.
Naturally, limitations inherent to IntervalRootFinding.jl apply here.

The output of `fixedpoints` can be used in the [BifurcationKit.jl](https://github.com/rveltz/BifurcationKit.jl)
as a start of a continuation process. See also [`periodicorbits`](@ref).

## Keywords
* `method = IntervalRootFinding.Krawczyk` configures the root finding method,
  see the docs of IntervalRootFinding.jl for all posibilities.
* `tol = 1e-15` is the root-finding tolerance.
"""
function fixedpoints(ds::DynamicalSystem, box, p = ds.p;
        method = IntervalRootFinding.Krawczyk, o = nothing, tol = 1e-15
    )
    if DynamicalSystemsBase.isinplace(ds)
        error("`fixedpoints` works only for out-of-place dynamical systems.")
    end
    # Find roots via IntervalRootFinding.jl
    f = to_root_f(ds, p, o)
    jac = to_root_J(ds, p, o)
    r = IntervalRootFinding.roots(f, jac, box, method, tol)
    D = dimension(ds)
    fp::Dataset{D, Float64} = roots_to_dataset(r, D)
    # Find eigenvalues and stability
    eigs = Vector{Vector{Complex{Float64}}}(undef, length(fp))
    J = Array(jacobian(ds)) # `eigvals` doesn't work with StaticArrays.jl
    for (i, u) in enumerate(fp)
        J .= jacobian(ds, u, p, 0.0) # notice that we use the "pure" jacobian, no -u!
        eigs[i] = LinearAlgebra.eigvals(Array(J))
    end
    stable = Bool[isstable(ds, e) for e in eigs]
    return fp, eigs, stable
end

to_root_f(ds::CDS, p, ::Nothing) = u -> ds.f(u, p, 0.0)
to_root_J(ds::CDS, p, ::Nothing) = u -> ds.jacobian(u, p, 0.0)
to_root_f(ds::DDS, p, ::Nothing) = u -> ds.f(u, p, 0) - u
function to_root_J(ds::DDS{IIP, S}, p, ::Nothing) where {IIP, S}
    c = Diagonal(ones(S))
    return u -> ds.jacobian(u, p, 0) - c
end

# Discrete with periodic order
function to_root_f(ds::DDS, p, o::Int)
    u -> begin
        v = copy(u) # copy is free for StaticArrays
        for _ in 1:o
            v = ds.f(v, p, 0)
        end
        return v - u
    end
end
# TODO: Estimate periodic orbits of period `o` for discrete systems.
# What's left is to create the Jacobian for higher iterates.

function roots_to_dataset(r, D)
    if isempty(r)
        @warn "No fixed points found!"
        return Dataset{D, Float64}()
    end
    if any(root.status != :unique for root in r)
        @warn "Non-unique fixed points found!"
    end
    F = zeros(length(r), D)
    for (j, root) in enumerate(r)
        F[j, :] .= map(i -> (i.hi + i.lo)/2, root.interval)
    end
    return Dataset(F; warn = false)
end

isstable(::CDS, e) = maximum(real(x) for x in e) < 0
isstable(::DDS, e) = maximum(abs(x) for x in e) < 1


text:
function infilt!(I, b, P)
  I  .= P.α.*(b.+P.q*P.f)./(b.+P.q)
end

function infilt(b::Float64, P)
  return P.α*(b+P.q*P.f)/(b+P.q)
end

function infilt(b, P)
  return P.α.*(b.+P.q*P.f)./(b.+P.q)
end

function rhs_h!(ht, h, p, t)
     In, P = p
@.   ht = P.p - In*h + P.dh*$laplacian(h.*h,P.dx)
     nothing
end

function steadyh!(h, In, P)
    probh=SteadyStateProblem(rhs_h!, h, (In, P))
    hh=solve(probh, DynamicSS(Tsit5(), tspan=0.2*P.dx*P.dx/P.dw), save_everystep=false, save_start = false);
    h.=hh.u
    nothing
end

function rhs_stat!(ut, u, p, t)

    global iint

    P, fg, σ, h, bint, wint = p
    b = @view u[:,:,1]
    w = @view u[:,:,2]
    bt = @view ut[:,:,1]
    wt = @view ut[:,:,2]

    I = infilt(b, P)

    # Compute steady h solution
    steadyh!(h, I, P);

    # Compute integrals only every nint steps
    if(mod(iint,P.nint)==0)
        α = initapprox(b, P.η, σ)
        bint .= approxintb(w, fg, α)
        wint .= approxintw(b, fg, α)
        iintb = 0
    end
    iint += 1

@.  bt = P.ν*b*(1-b)*bint - b  + P.db*$laplacian(b, P.dx)
@.  wt = I*h - P.ν*w/(1+P.ρ*b) - P.γ*w*wint + P.dw*$laplacian(w, P.dx)

    nothing
end

max_stars_repo_name: jhardenberg/bwh

max_issues_repo_name: jhardenberg/bwh

max_forks_repo_name: jhardenberg/bwh


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")
include("cylindrical_particle_exact_solution.jl")

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function solve_and_compute_variational_stress_error(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)
    transfstress =
        CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

    analyticalsolution = AnalyticalSolution(
        inradius,
        outradius,
        center,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> -circle_distance_function(x, center, inradius),
        mesh,
    )

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads =
        CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads =
        CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads =
        CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    displacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, width, width),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(
        sysmatrix,
        interfacecondition,
        cutmesh,
    )
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_coherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        displacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onboundary(x, width, width),
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    nodaldisplacement = matrix \ rhs

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_stress_mass_matrix!(sysmatrix, basis, cellquads, cutmesh)
    CutCell.assemble_stress_linear_form!(
        sysrhs,
        basis,
        cellquads,
        stiffness,
        nodaldisplacement,
        cutmesh,
    )
    CutCell.assemble_transformation_stress_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )

    matrix = CutCell.make_sparse_stress_operator(sysmatrix, cutmesh)
    rhs = CutCell.stress_rhs(sysrhs, cutmesh)

    stressvec = matrix \ rhs
    stress = reshape(stressvec, 3, :)

    err = mesh_L2_error(
        stress,
        x -> exact_stress(analyticalsolution, x)[1:3],
        basis,
        cellquads,
        cutmesh,
    )

    return err
end

K1, K2 = 247.0, 192.0
mu1, mu2 = 126.0, 87.0
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)

stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

theta0 = -0.067

width = 1.0

center = [width / 2, width / 2]
inradius = width / 4
outradius = width

polyorder = 2
numqp = required_quadrature_order(polyorder) + 2
penaltyfactor = 1e2

powers = [3,4,5]
nelmts = [2^p + 1 for p in powers]

err = [solve_and_compute_variational_stress_error(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    theta0,
    ne,
    polyorder,
    numqp,
    penaltyfactor,
) for ne in nelmts]

dx = 1.0 ./ nelmts
serr = [[er[i] for er in err] for i = 1:3]
rates = [convergence_rate(v,dx) for v in serr]


text:
module Test_riccati

using Random
using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing algebraic Riccati equation solvers" begin

# only double precision tests are performed

#for (Ty,n,p,m) in ((Float64, 30, 10, 10), (Float32, 5, 3, 3))
for (Ty,n,p,m) in ((Float64, 30, 10, 10),  (Float64, 5, 3, 1))

Random.seed!(21235)

#(Ty,n,p,m) = (Float64, 30, 10, 10)
#(Ty,n,p,m) = (Float64, 2, 1, 1)

ar = randn(Ty,n,n)
er = rand(Ty,n,n)
if m == 1
    br = rand(Ty,n)
    bc = rand(Ty,n)+im*rand(Ty,n)
    sc = rand(Complex{Ty},n)/100
    sr = rand(Ty,n)/100
else
    br = rand(Ty,n,m)
    bc = rand(Ty,n,m)+im*rand(Ty,n,m)
    sc = rand(Complex{Ty},n,m)/100
    sr = rand(Ty,n,m)/100
end
cr = rand(Ty,p,n)
ac = randn(Ty,n,n) + im*randn(Ty,n,n)
ec = er+im*rand(Ty,n,n)
cc = rand(Ty,p,n)+im*rand(Ty,p,n)
bc = rand(Ty,n,m)+im*rand(Ty,n,m)
rr1 = rand(Ty,m,m)
rc1 = rand(Ty,m,m)+im*rand(Ty,m,m)
Ty == Float64 ? rtol = n*n*sqrt(eps(1.)) : rtol = n*sqrt(eps(1.f0))

qc = cc'*cc
Qr = cr'*cr
gc = bc*bc'
gr = br*br'
rr = rr1*rr1'
rc = rc1*rc1'

@testset "Continuous Riccati equation ($Ty, n = $n, m = $m)" begin

@time x, clseig = arec(ar,gr,Qr)
@test norm(ar'*x+x*ar-x*gr*x+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,gr,Qr,as = true)
@test norm(ar'*x+x*ar-x*gr*x+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,gr)
@test norm(ar'*x+x*ar-x*gr*x)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,2I,I)
@test norm(ar'*x+x*ar-x*2*x+I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-2*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-2*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,2,I)
@test norm(ar'*x+x*ar-x*2*x+I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-2*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-2*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,gr,2I)
@test norm(ar'*x+x*ar-x*gr*x+2I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,gr,2)
@test norm(ar'*x+x*ar-x*gr*x+2I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,2I,Qr)
@test norm(ar'*x+x*ar-x*2*x+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-2*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-2*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar',gr,Qr)
@test norm(ar*x+x*ar'-x*gr*x+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar'-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar'-gr*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ac,gc,qc)
@test norm(ac'*x+x*ac-x*gc*x+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-gc*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-gc*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ac',gc,qc)
@test norm(ac*x+x*ac'-x*gc*x+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac'-gc*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac'-gc*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ar,gc,qc)
@test norm(ar'*x+x*ar-x*gc*x+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gc*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gc*x))))/norm(clseig)  < rtol

@time x, clseig = arec(ac,gr,Qr)
@test norm(ac'*x+x*ac-x*gr*x+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-gr*x))))/norm(clseig)  < rtol

try
    arec(zeros(2,2),[1 1;1 1],[1 0;0 0])  # Hamiltonian not dichotomic
    @test false
catch
    @test true
end

try
    arec(rand(2,2),0I,0I)  # no finite solution
    @test false
catch
    @test true
end

end

@testset "Generalized continuous Riccati equation ($Ty, n = $n, m = $m)" begin

@time x, clseig = garec(ar,er,gr,Qr)
@test norm(ar'*x*er+er'*x*ar-er'*x*gr*x*er+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol

@time x, clseig = garec(ar,er,gr,Qr,as = true)
@test norm(ar'*x*er+er'*x*ar-er'*x*gr*x*er+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol

@time x, clseig = garec(ar,er,gr)
@test norm(ar'*x*er+er'*x*ar-er'*x*gr*x*er)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol

@time x, clseig = garec(ar,er,gr,2I)
@test norm(ar'*x*er+er'*x*ar-er'*x*gr*x*er+2I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-gr*x*er,er))))/norm(clseig)  < rtol

@time x, clseig = garec(ar,er,I,2I)
@test norm(ar'*x*er+er'*x*ar-er'*x*x*er+2I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-x*er,er))))/norm(clseig)  < rtol

@time x, clseig = garec(ar,I,1,2)
@test norm(ar'*x+x*ar-x*x+2I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-x))))/norm(clseig)  < rtol

end

@testset "Continuous control Riccati equation ($Ty, n = $n, m = $m)" begin

@time x, clseig, f, z = arec(ar,br,rr,Qr,sr)
@test norm(ar'*x+x*ar-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f, z = arec(ar,br,rr,Qr,sr;orth=true)
@test norm(ar'*x+x*ar-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol
norm(z'*z - I) < rtol

@time x, clseig, f = arec(ar,br,rr,Qr,sr,as=true)
@test norm(ar'*x+x*ar-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,rr,Qr)
@test norm(ar'*x+x*ar-x*br*inv(rr)*br'*x+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,I,Qr)
@test norm(ar'*x+x*ar-x*br*br'*x+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,1,Qr)
@test norm(ar'*x+x*ar-x*br*br'*x+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,rr,0I)
@test norm(ar'*x+x*ar-x*br*inv(rr)*br'*x)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,rr,0)
@test norm(ar'*x+x*ar-x*br*inv(rr)*br'*x)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ac,bc,rc,qc)
@test norm(ac'*x+x*ac-x*bc*inv(rc)*bc'*x+qc)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ac,bc,rc,qc,sc)
@test norm(ac'*x+x*ac-(x*bc+sc)*inv(rc)*(bc'*x+sc')+qc)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ac,bc,rr,Qr,sr)
@test norm(ac'*x+x*ac-(x*bc+sr)*inv(rr)*(bc'*x+sr')+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,gr,rr,Qr,sr)
@test norm(ar'*x+x*ar-x*gr*x-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-gr*x))))/norm(clseig)  < rtol

@time x, clseig, f = arec(ar,br,gr,rr,Qr,sr,as=true)
@test norm(ar'*x+x*ar-x*gr*x-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x))  < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-gr*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-gr*x))))/norm(clseig)  < rtol

g = 2
@time x, clseig, f, z = arec(ar,br,g*I,rr,Qr,sr)
@test norm(ar'*x+x*ar-x*g*x-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-g*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-g*x))))/norm(clseig)  < rtol

g = 2
@time x, clseig, f, z = arec(ar,br,g*I,rr,Qr,sr,as=true)
@test norm(ar'*x+x*ar-x*g*x-(x*br+sr)*inv(rr)*(br'*x+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-g*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-g*x))))/norm(clseig)  < rtol

g = 2; q = 3. *I; r = 1;
@time x, clseig, f, z = arec(ar,br,g*I,r,q,sr,as=true)
@test norm(ar'*x+x*ar-x*g*x-(x*br+sr)*inv(r)*(br'*x+sr')+q)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-g*x))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-g*x))))/norm(clseig)  < rtol

ar1 = rand(Ty,2,2); br1 = rand(Ty,2,2); gr1 = 0*I; r1 = [1.e5 0; 0 1.e-5]; q1 = 0*I; sr1 = zeros(Ty,2,2);
@time x, clseig, f, z = arec(ar1,br1,gr1,r1,q1,sr1)
@test norm(ar1'*x+x*ar1-x*gr1*x-(x*br1+sr1)*inv(r1)*(br1'*x+sr1')+q1)/max(1,norm(x)) < rtol

end

@testset "Generalized continuous control Riccati equation ($Ty, n = $n, m = $m)" begin

@time x, clseig, f = garec(ar,er,br,rr,Qr,sr)
@test norm(ar'*x*er+er'*x*ar-(er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ar,er,br,rr,Qr,sr,as = true)
@test norm(ar'*x*er+er'*x*ar-(er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ar,er,br,gr,rr,Qr,sr)
@test norm(ar'*x*er+er'*x*ar-er'*x*gr*x*er - (er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-gr*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-gr*x*er,er))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ar,er,br,gr,rr,Qr,sr,as = true)
@test norm(ar'*x*er+er'*x*ar-er'*x*gr*x*er - (er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-gr*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-gr*x*er,er))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ar,er,br,rr,2I,sr)
@test norm(ar'*x*er+er'*x*ar-(er'x*br+sr)*inv(rr)*(br'*x*er+sr')+2I)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ac,ec,bc,rc,qc,sc)
@test norm(ac'*x*ec+ec'*x*ac-(ec'x*bc+sc)*inv(rc)*(bc'*x*ec+sc')+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f,ec))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f,ec))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ac,ec,bc,rr,Qr,sr)
@test norm(ac'*x*ec+ec'*x*ac-(ec'x*bc+sr)*inv(rr)*(bc'*x*ec+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f,ec))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f,ec))))/norm(clseig)  < rtol

@time x, clseig, f = garec(ar,er,br,1,rr,Qr,sr)
@test norm(ar'*x*er+er'*x*ar-er'*x*x*er-(er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-x*er,er))))/norm(clseig)  < rtol

g = 2
@time x, clseig, f, z = garec(ar,er,br,g*I,rr,Qr,sr)
@test norm(ar'*x*er+er'*x*ar-er'*x*g*x*er-(er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-g*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-g*x*er,er))))/norm(clseig)  < rtol

g = 2
@time x, clseig, f = garec(ar,er,br,g,rr,Qr,sr,as=true)
@test norm(ar'*x*er+er'*x*ar-er'*x*g*x*er-(er'x*br+sr)*inv(rr)*(br'*x*er+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-g*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-g*x*er,er))))/norm(clseig)  < rtol

g = 2; q = 3. *I; r = 1;
@time x, clseig, f = garec(ar,er,br,g,r,q,sr,as=true)
@test norm(ar'*x*er+er'*x*ar-er'*x*g*x*er-(er'x*br+sr)*inv(r)*(br'*x*er+sr')+q)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f-g*x*er,er))))/norm(clseig)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f-g*x*er,er))))/norm(clseig)  < rtol

end

@testset "Discrete control Riccati equation ($Ty, n = $n, m = $m)" begin
@time x, clseig, f = ared(ar,br,rr,Qr)
@test norm(ar'*x*ar-x-ar'*x*br*inv(rr+br'*x*br*I)*br'*x*ar+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(ar)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(ar)  < rtol

@time x, clseig, f = ared(ar,br,rr,0I)
@test norm(ar'*x*ar-x-ar'*x*br*inv(rr+br'*x*br*I)*br'*x*ar)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(ar)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(ar)  < rtol

@time x, clseig, f = ared(ar,br,rr,0I; as = true)
@test norm(ar'*x*ar-x-ar'*x*br*inv(rr+br'*x*br*I)*br'*x*ar)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(ar)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(ar)  < rtol

@time x, clseig, f = ared(ar,br,rr,Qr,sr)
@test norm(ar'*x*ar-x-(ar'*x*br+sr)*inv(rr+br'*x*br*I)*(br'*x*ar+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(ar)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(ar)  < rtol

@time x, clseig, f = ared(ar,br,rr,Qr,sr,as = true)
@test norm(ar'*x*ar-x-(ar'*x*br+sr)*inv(rr+br'*x*br*I)*(br'*x*ar+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/norm(ar)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/norm(ar)  < rtol

@time x, clseig, f = ared(ac,bc,rc,qc)
@test norm(ac'*x*ac-x-ac'*x*bc*inv(rc+bc'*x*bc*I)*bc'*x*ac+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f))))/norm(ac)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f))))/norm(ac)  < rtol

@time x, clseig, f = ared(ac,bc,rc,qc,sc)
@test norm(ac'*x*ac-x-(ac'*x*bc+sc)*inv(rc+bc'*x*bc*I)*(bc'*x*ac+sc')+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f))))/norm(ac)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f))))/norm(ac)  < rtol

q = 3.f0 *I; r = 1;
@time x, clseig, f = ared(ar,br,r,q,sr,as=true)
@test norm(ar'*x*ar-x-(ar'*x*br+sr)*inv(r*I+br'*x*br)*(br'*x*ar+sr')+q)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f))))/max(opnorm(ar,1),1)  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f))))/max(opnorm(ar,1),1)  < rtol

end

@testset "Generalized discrete control Riccati equation ($Ty, n = $n, m = $m)" begin
@time x, clseig, f = gared(ar,er,br,rr,Qr,sr)
@test norm(ar'*x*ar-er'*x*er-(ar'*x*br+sr)*inv(rr+br'*x*br*I)*(br'*x*ar+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol

@time x, clseig, f = gared(ar,er,br,rr,Qr,sr,as = true)
@test norm(ar'*x*ar-er'*x*er-(ar'*x*br+sr)*inv(rr+br'*x*br*I)*(br'*x*ar+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol

@time x, clseig, f = gared(ar,er,br,rr,0I)
@test norm(ar'*x*ar-er'*x*er-(ar'*x*br)*inv(rr+br'*x*br*I)*(br'*x*ar))/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol

@time x, clseig, f, z = gared(ar,er,br,rr,0I;as = true)
@test norm(ar'*x*ar-er'*x*er-(ar'*x*br)*inv(rr+br'*x*br*I)*(br'*x*ar))/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol

@time x, clseig, f = gared(ac,ec,bc,rc,qc,sc)
@test norm(ac'*x*ac-ec'*x*ec-(ac'*x*bc+sc)*inv(rc+bc'*x*bc*I)*(bc'*x*ac+sc')+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f,ec))))/max(opnorm(ac,1),opnorm(ec,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f,ec))))/max(opnorm(ac,1),opnorm(ec,1))  < rtol

@time x, clseig, f = gared(ac,ec,bc,rc,qc,sc; as = true)
@test norm(ac'*x*ac-ec'*x*ec-(ac'*x*bc+sc)*inv(rc+bc'*x*bc*I)*(bc'*x*ac+sc')+qc)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f,ec))))/max(opnorm(ac,1),opnorm(ec,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f,ec))))/max(opnorm(ac,1),opnorm(ec,1))  < rtol

@time x, clseig, f = gared(ac,ec,bc,rr,Qr,sr)
@test norm(ac'*x*ac-ec'*x*ec-(ac'*x*bc+sr)*inv(rr+bc'*x*bc*I)*(bc'*x*ac+sr')+Qr)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ac-bc*f,ec))))/max(opnorm(ac,1),opnorm(ec,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ac-bc*f,ec))))/max(opnorm(ac,1),opnorm(ec,1))  < rtol

q = 3.f0 *I; r = 1;
@time x, clseig, f = gared(ar,er,br,r,q,sr,as=true)
@test norm(ar'*x*ar-er'*x*er-(ar'*x*br+sr)*inv(r*I+br'*x*br)*(br'*x*ar+sr')+q)/max(1,norm(x)) < rtol &&
norm(sort(real(clseig))-sort(real(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol &&
norm(sort(imag(clseig))-sort(imag(eigvals(ar-br*f,er))))/max(opnorm(ar,1),opnorm(er,1))  < rtol

end

end

end

end


text:
function setup()
	lines = readlines("inputD13.txt")
	dots = Tuple{Int, Int}[]
	folds = Tuple{Char, Int}[]

	for line in lines
		if startswith(line, "fold")
			left, right = split(line, "=")
			left = last(left)
			right = parse(Int, right)
			push!(folds, (left, right))
		elseif line != ""
				x, y = [parse(Int, i) for i in split(line, ",")]
				push!(dots, (x, y))
		end
	end

	return (dots, folds)
end

function folddots(folddir :: Char, foldcoord :: Int, dots :: Vector{Tuple{Int, Int}})
	newdots = Tuple{Int, Int}[]
	if folddir == 'x'
		for dot in dots
			x, y = dot
			x < foldcoord ? push!(newdots, dot) : push!(newdots, (foldcoord * 2 - x, y))
		end

	elseif folddir == 'y'
		for dot in dots
			x, y = dot
			y < foldcoord ? push!(newdots, dot) : push!(newdots, (x, foldcoord * 2 - y))
		end
	end

	return newdots
end

function part1()
	dots, folds = setup()

	folddir, foldcoord = first(folds)
	newdots = folddots(folddir, foldcoord, dots)

	return length(Set(newdots))
end

function part2()
	dots, folds = setup()
	for (folddir, foldcoord) in folds
		dots = folddots(folddir, foldcoord, dots)
	end
	dots = Set(dots)

	highest = maximum(maximum.(dots))
	for i in 0:5
		for j in 0:highest
			print((j, i) in dots ? "#" : " ")
		end
		println()
	end

end

println(part1())
part2()


text:
"""
	Session

This is a web session.

```jldoctest
julia> capabilities = Capabilities("chrome")
Remote WebDriver Capabilities
browserName: chrome
julia> wd = RemoteWebDriver(capabilities, host = ENV["WEBDRIVER_HOST"], port = parse(Int, ENV["WEBDRIVER_PORT"]))
Remote WebDriver
julia> session = Session(wd)
Session
julia> isa(session, Session)
true
julia> delete!(session);

```
"""
struct Session{D<:Object}
    addr::String
    id::String
    attrs::D
    function Session(wd::RemoteWebDriver)
        @unpack addr = wd
        response = HTTP.post(
            "$(wd.addr)/session",
            [("Content-Type" => "application/json")],
            JSON3.write(Dict("desiredCapabilities" => wd.capabilities, wd.kw...)),
        )
        @assert response.status == 200
        json = JSON3.read(response.body)
        new{typeof(json.value)}(addr, json.sessionId, json.value)
    end
end
broadcastable(obj::Session) = Ref(obj)
summary(io::IO, obj::Session) = println(io, "Session")
function show(io::IO, obj::Session)
    print(io, summary(obj))
end


text:
module AnalyticalPlaneSolver

struct PlaneSolver
    L::Any
    W::Any
    R::Any
    lambda1::Any
    mu1::Any
    lambda2::Any
    mu2::Any
    theta0::Any
    C1::Any
    C2::Any
    D::Any
    uI::Any
    function PlaneSolver(L, W, R, lambda1, mu1, lambda2, mu2, theta0)
        C1 = solver_coefficient(lambda1, mu1)
        C2 = solver_coefficient(lambda2, mu2)
        D = 2 * lambda1 * mu1 / (lambda1 + 2mu1)

        uI = interface_displacement(L, R, C1, C2, D, theta0)
        new(L, W, R, lambda1, mu1, lambda2, mu2, theta0, C1, C2, D, uI)
    end
end

function solver_coefficient(lambda, mu)
    C = 4 * mu * (lambda + mu) / (lambda + 2mu)
    return C
end

function interface_displacement(L, R, C1, C2, D, theta0)
    uI = -(C1 + D) / (C2 / R + C1 / (L - R)) * theta0 / 3
    return uI
end

function parent_displacement_field(solver::PlaneSolver,x)
    uI = solver.uI
    R = solver.R
    lambda = solver.lambda2
    mu = solver.mu2

    ux = uI/R*x[1]
    uy = -lambda/(lambda+2mu)*uI/R*x[2]
    return [ux,uy]
end

function product_displacement_field(solver::PlaneSolver,x)
    uI = solver.uI
    L = solver.L
    R = solver.R
    lambda = solver.lambda1
    mu = solver.mu1
    t0 = solver.theta0

    ux = ((L-R)-(x[1]-R))/(L-R)*uI
    uy = (lambda/(lambda+2mu)*uI/(L-R) + (3lambda+2mu)/(lambda+2mu)*t0/3)*x[2]

    return [ux,uy]
end

function displacement_field(solver::PlaneSolver,x)
    R = solver.R
    if x[1] < R
        return parent_displacement_field(solver,x)
    else
        return product_displacement_field(solver,x)
    end
end

end


text:
using Statistics
using ProgressMeter
using MATLAB
include("atn_module.jl")
N = 10_000          # number of nodes
m = 5               # number of active links per active node
mu = 0.1            # probability I -> S
tmax = 3000         # maximum duration of the simulation
Ni = Int(0.01 * N)  # initial number of infects
gamma = -2.22       # exponent of the power law
NTrials = 5         # number of trials for each value of lambda
tolerance = 1e-10   # termination tolerance
TWindow = 500      # time window for evaluating steady state
TMinConvergence = 500

LambdaMin = 0.0
LambdaMax = 1.0
LambdaStep = 0.05
Lambda = LambdaMin:LambdaStep:LambdaMax
LambdaSize = length(Lambda)

EtaMin = 0.0
EtaMax = 15
EtaStep = 0.5
Eta = EtaMin:EtaStep:EtaMax
EtaSize = length(Eta)
IT = zeros(EtaSize, LambdaSize)
@showprogress for (i, lambda) in enumerate(Lambda)
    for (j, eta) in enumerate(Eta)
        ITrend = time_series_eta(Ni, N, m, eta, EtaMax, gamma, lambda, mu, tolerance,
                                 TWindow, TMinConvergence, tmax)
        IT[j, i] = ITrend[end]
    end
end
# mat"plot($ITrend)"


text:
#' @description
#'
#' An AbstractDataArray is an Array whose entries can take on
#' values of type `T` or the value `NA`.
abstract AbstractDataArray{T, N} <: AbstractArray{T, N}

#' @description
#'
#' An AbstractDataVector is an AbstractDataArray of order 1.
typealias AbstractDataVector{T} AbstractDataArray{T, 1}

#' @description
#'
#' An AbstractDataMatrix is an AbstractDataArray of order 2.
typealias AbstractDataMatrix{T} AbstractDataArray{T, 2}

#' @description
#' Determine the type of the elements of an AbstractDataArray.
#'
#' @param ada::AbstractDataArray{T} The AbstractDataArray whose
#'        element type is desired.
#'
#' @returns T::DataType The type of the non-NA elements of `ada`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' T = eltype(dv)
Base.eltype{T, N}(d::AbstractDataArray{T, N}) = T

# Generic iteration over AbstractDataArray's

# TODO: Document
Base.start(x::AbstractDataArray) = 1

# TODO: Document
Base.next(x::AbstractDataArray, state::Integer) = (x[state], state + 1)

# TODO: Document
Base.done(x::AbstractDataArray, state::Integer) = state > length(x)

#' @description
#'
#' Determine if the values of an AbstractArray are `NA`.
#'
#' @param a::AbstractArray{T, N} The AbstractArray whose missingness will
#'        be assessed.
#'
#' @returns na::BitArray{N} Elementwise Boolean whether entry is missing.
#'
#' @examples
#'
#' a = [1, 2, 3]
#' isna(a)
isna{T}(a::AbstractArray{T}) =
    NAtype <: T ? BitArray(map(x->isa(x, NAtype), a)) : falses(size(a)) # -> BitArray

#' @description
#'
#' Safe and type-stable way to determine if element `i` of an
#' AbstractArray is `NA`.
#'
#' @param a::AbstractArray The AbstractArray whose missingness will
#'        be assessed.
#' @param i::Integer The index of the element to be checked for `NA`.
#'
#' @returns na::Bool Is the element `NA` or not?
#'
#' @examples
#'
#' a = [1, 2, 3]
#' isna(a, 1)
isna{T}(a::AbstractArray{T}, i::Real) = NAtype <: T ? isa(a[i], NAtype) : false # -> Bool

#' @description
#'
#' Determine if any of the entries of an AbstractArray are `NA`.
#'
#' @param a::AbstractArray{T, N} The AbstractArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are any of the elements of `a` an `NA` value?
#'
#' @examples
#'
#' a = [1, 2, 3]
#' anyna(a)
anyna(a::AbstractArray) = any(isna(a)) # -> Bool

#' @description
#'
#' Determine if all of the entries of an AbstractArray are `NA`.
#'
#' @param a::AbstractArray{T, N} The AbstractArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are all of the elements of `a` an `NA` value?
#'
#' @examples
#'
#' a = [1, 2, 3]
#' allna(a)
allna(a::AbstractArray) = all(isna(a)) # -> Bool

#' @description
#'
#' NO-OP: Turn a Vector into a Vector. See dropna(dv::DataVector) for
#'        rationale.
#'
#' @param v::Vector{T} Vector that will be converted to a Vector.
#'
#' @returns v::Vector{T} Vector containing all of the values of `v`.
#'
#' @examples
#'
#' v = [1, 2, 3, 4]
#' v = dropna(v)
dropna(v::AbstractVector) = copy(v) # -> AbstractVector

# Iterators
# TODO: Use values()
#       Use DataValueIterator type?

type EachFailNA{T}
    da::AbstractDataArray{T}
end
each_failna{T}(da::AbstractDataArray{T}) = EachFailNA(da)
Base.length(itr::EachFailNA) = length(itr.da)
Base.start(itr::EachFailNA) = 1
Base.done(itr::EachFailNA, ind::Integer) = ind > length(itr)
function Base.next(itr::EachFailNA, ind::Integer)
    if isna(itr.da[ind])
        throw(NAException())
    else
        (itr.da[ind], ind + 1)
    end
end

type EachDropNA{T}
    da::AbstractDataArray{T}
end
each_dropna{T}(da::AbstractDataArray{T}) = EachDropNA(da)
function _next_nonna_ind{T}(da::AbstractDataArray{T}, ind::Int)
    ind += 1
    while ind <= length(da) && isna(da, ind)
        ind += 1
    end
    ind
end
Base.length(itr::EachDropNA) = length(itr.da) - sum(itr.da.na)
Base.start(itr::EachDropNA) = _next_nonna_ind(itr.da, 0)
Base.done(itr::EachDropNA, ind::Int) = ind > length(itr.da)
function Base.next(itr::EachDropNA, ind::Int)
    (itr.da[ind], _next_nonna_ind(itr.da, ind))
end

type EachReplaceNA{S, T}
    da::AbstractDataArray{S}
    replacement::T
end
function each_replacena(da::AbstractDataArray, replacement::Any)
    EachReplaceNA(da, convert(eltype(da), replacement))
end
function each_replacena(replacement::Any)
    x -> each_replacena(x, replacement)
end
Base.length(itr::EachReplaceNA) = length(itr.da)
Base.start(itr::EachReplaceNA) = 1
Base.done(itr::EachReplaceNA, ind::Integer) = ind > length(itr)
function Base.next(itr::EachReplaceNA, ind::Integer)
    item = isna(itr.da, ind) ? itr.replacement : itr.da[ind]
    (item, ind + 1)
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using JLD2
using BenchmarkTools
using Benchmarks
using Plots
pyplot()
# Benchmark parameters

print_system_info()

#set to true to use strong scaling, false to use weak scaling
strong = false
#set to true to use threads, false to use mpi
threaded = false

decomposition = Slab()
# decomposition = Pencil()

#automatically generates thread tuple
#Slab() => min.(2 .^ (0:10), Sys.CPU_THREADS) |> unique
ranks = Dict(
    Slab() => (1, 2, 4),
    Pencil() => (1, 4, 16)
)[decomposition]

if strong
    grid_size(R, decomposition) = (256, 256, 256)
    title = "strong"
else
    grid_size(R, decomposition) = (128, 128, 16R)
    title = "weak"
end

if threaded
    command(julia, R, Nx, Ny, Nz, Rx, Ry, Rz) = `$julia --threads $R --project distributed_nonhydrostatic_model_threaded.jl $Nx $Ny $Nz`
    label = "threads"
    dis_type = "threaded"
    keyop = v -> v[2]
else
    command(julia, R, Nx, Ny, Nz, Rx, Ry, Rz) = `mpiexec -np $R $julia --project distributed_nonhydrostatic_model_mpi.jl $Nx $Ny $Nz $Rx $Ry $Rz`
    label = "ranks"
    dis_type = "mpi"
    keyop = v -> v[2][2]
end

rank_size(R, ::Slab) = (1, R, 1)
rank_size(R, ::Pencil) = Int.((1, √R, √R))

# Run benchmarks
for R in ranks
    Nx, Ny, Nz = grid_size(R, decomposition)
    Rx, Ry, Rz = rank_size(R, decomposition)
    @info string("Benchmarking ", title, " scaling nonhydrostatic model with $(typeof(decomposition)) decomposition [N=($Nx, $Ny, $Nz), ", label, "=($Rx, $Ry, $Rz)]...")
    julia = Base.julia_cmd()
    run(command(julia, R, Nx, Ny, Nz, Rx, Ry, Rz))
end

# Collect and merge benchmarks from all ranks
suite = BenchmarkGroup(["size", label])

for R in ranks
    Nx, Ny, Nz = grid_size(R, decomposition)
    Rx, Ry, Rz = rank_size(R, decomposition)

    if threaded
        case = ((Nx, Ny, Nz), R)
        filename = string("distributed_nonhydrostatic_model_threads$(R).jld2")
        file = jldopen(filename, "r")
        suite[case] = file["trial"]
    else
        case = ((Nx, Ny, Nz), (Rx, Ry, Rz))
        for local_rank in 0:R-1
            filename = string("distributed_nonhydrostatic_model_$(R)ranks_$local_rank.jld2")
            jldopen(filename, "r") do file
                if local_rank == 0
                    suite[case] = file["trial"]
                else
                    merged_trial = suite[case]
                    local_trial = file["trial"]
                    append!(merged_trial.times, local_trial.times)
                    append!(merged_trial.gctimes, local_trial.gctimes)
                end
            end
        end
    end
end

plot_keys = collect(keys(suite))
sort!(plot_keys, by = keyop)
plot_num = length(plot_keys)
rank_num = zeros(Int64, plot_num)
run_times = zeros(Float64, plot_num)
eff_ratio = zeros(Float64, plot_num)
for i in 1:plot_num
    run_times[i] = mean(suite[plot_keys[i]].times) / 1.0e6
    threaded ? rank_num[i] = plot_keys[i][2] : rank_num[i] = plot_keys[i][2][2]
    eff_ratio[i] = median(suite[plot_keys[1]].times) / ((rank_num[i]^strong) * median(suite[plot_keys[i]].times))
end

plt = plot(rank_num, run_times, lw=4, xaxis=:log2, legend=:none,
          xlabel=label, ylabel="Times (ms)", title=string(dis_type, " ", title, " scaling nonhydrostatic times"))
display(plt)
savefig(plt, string(dis_type, "_", title, "_nonhydrostatic_times.png"))

plt2 = plot(rank_num, eff_ratio, lw=4, xaxis=:log2, legend=:none, ylims=(0,1.1),
            xlabel=label, ylabel="Efficiency", title=string(dis_type, " ", title, " scaling nonhydrostatic efficiency"))
display(plt2)
savefig(plt2, string(dis_type, "_", title, "_nonhydrostatic_efficiency.png"))

# Summarize benchmarks

df = benchmarks_dataframe(suite)
sort!(df, Symbol(label))
benchmarks_pretty_table(df, title=string(dis_type, " ", title," scaling nonhydrostatic times"))

base_case = (grid_size(ranks[1], decomposition), threaded ? ranks[1] : rank_size(ranks[1], decomposition))
suite_Δ = speedups_suite(suite, base_case=base_case)
df_Δ = speedups_dataframe(suite_Δ, slowdown=true, efficiency=Symbol(title), base_case=base_case, key2rank=k->prod(k[2]))
sort!(df_Δ, Symbol(label))
benchmarks_pretty_table(df_Δ, title=string(dis_type, " ", title," scaling nonhydrostatic efficiency"))


text:
using Distributed
@everywhere using MCPhylo

## Data
inhalers = Dict{Symbol, Any}(
  :pattern =>
    [1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
     1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4]',
  :Ncum =>
    [ 59 157 173 175 186 253 270 271 271 278 280 281 282 285 285 286
     122 170 173 175 226 268 270 271 278 280 281 281 284 285 286 286]',
  :treat =>
    [ 1 -1
     -1  1],
  :period =>
    [1 -1
     1 -1],
  :carry =>
    [0 -1
     0  1],
  :N => 286,
  :T => 2,
  :G => 2,
  :Npattern => 16,
  :Ncut => 3
)

inhalers[:group] = Array{Int}(undef, inhalers[:N])
inhalers[:response] = Array{Int}(undef, inhalers[:N], inhalers[:T])

i = 1
for k in 1:inhalers[:Npattern], g in 1:inhalers[:G]
  global i
  while i <= inhalers[:Ncum][k, g]
    inhalers[:group][i] = g
    for t in 1:inhalers[:T]
      inhalers[:response][i, t] = inhalers[:pattern][k, t]
    end
    i += 1
  end
end

## Model Specification
model = Model(

  response = Stochastic(2,
    (a1, a2, a3, mu, group, b, N, T) ->
      begin
        a = Float64[a1, a2, a3]
        UnivariateDistribution[(
          eta = mu[group[i], t] + b[i];
          p = ones(4);
          for j in 1:3
            Q = invlogit(-(a[j] + eta))
            p[j] -= Q
            p[j + 1] = Q
          end;
          Categorical(p)) for i in 1:N, t in 1:T
        ]
      end,
    false
  ),

  mu = Logical(2,
    (beta, treat, pi, period, kappa, carry, G, T) ->
      [ beta * treat[g, t] / 2 + pi * period[g, t] / 2 + kappa * carry[g, t]
        for g in 1:G, t in 1:T ],
    false
  ),

  b = Stochastic(1,
    s2 -> Normal(0, sqrt(s2)),
    false
  ),

  a1 = Stochastic(
    a2 -> Truncated(Flat(), -1000, a2)
  ),

  a2 = Stochastic(
    a3 -> Truncated(Flat(), -1000, a3)
  ),

  a3 = Stochastic(
    () -> Truncated(Flat(), -1000, 1000)
  ),

  beta = Stochastic(
    () -> Normal(0, 1000)
  ),

  pi = Stochastic(
    () -> Normal(0, 1000)
  ),

  kappa = Stochastic(
    () -> Normal(0, 1000)
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:response => inhalers[:response], :beta => 0, :pi => 0, :kappa => 0,
       :a1 => 2, :a2 => 3, :a3 => 4, :s2 => 1, :b => zeros(inhalers[:N])),
  Dict(:response => inhalers[:response], :beta => 1, :pi => 1, :kappa => 0,
       :a1 => 3, :a2 => 4, :a3 => 5, :s2 => 10, :b => zeros(inhalers[:N]))
]

## Sampling Scheme
scheme = [AMWG(:b, 0.1),
          Slice([:a1, :a2, :a3], 2.0),
          Slice([:beta, :pi, :kappa, :s2], 1.0, Univariate)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, inhalers, inits, 5000, burnin=1000, thin=2, chains=2)
describe(sim)


text:
module TestsStackBuilding

using Test
using BGEN
using CSV
using TMLEEpistasis
using MLJ
using TOML
using TMLE
using Distributions

include("helper_fns.jl")

GLMClassifier = @load LinearBinaryClassifier pkg=GLM verbosity=0
LogisticClassifier = @load LogisticClassifier pkg=MLJLinearModels verbosity=0
KNNClassifier = @load KNNClassifier pkg=NearestNeighborModels verbosity=0

GLMRegressor = @load LinearRegressor pkg=GLM verbosity=0
LinearRegressor = @load LinearRegressor pkg=MLJLinearModels verbosity=0
KNNRegressor = @load KNNRegressor pkg=NearestNeighborModels verbosity=0

@testset "Categorical target TMLE built from configuration file" begin
    tmle_config = joinpath("config", "tmle_config.toml")
    build_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    tmles =  TMLEEpistasis.estimators_from_toml(TOML.parsefile(tmle_config), queries, TMLEEpistasis.PhenotypeTMLEEpistasis)
    # Test binary target TMLE's Qstack
    tmle = tmles["binary"]
    @test tmle.F isa GLMClassifier
    ## Checking Qstack.metalearner
    @test tmle.Q̅.metalearner isa LogisticClassifier
    @test tmle.Q̅.metalearner.fit_intercept == false
    ## Checking Qstack.resampling
    @test tmle.Q̅.resampling isa StratifiedCV
    @test tmle.Q̅.resampling.nfolds == 2
    ## Checking Qstack XGBoost models
    @test tmle.Q̅.XGBoostClassifier_1.num_round == 10
    ## Checking Qstack  Interaction Logistic models
    @test tmle.Q̅.InteractionLMClassifier_1 isa TMLEEpistasis.InteractionLMClassifier
    @test tmle.Q̅.InteractionLMClassifier_1.interaction_transformer.column_pattern == r"^RS_"
    ## Checking Qstack HAL model
    @test tmle.Q̅.HALClassifier_1.lambda == 10
    @test tmle.Q̅.HALClassifier_1.smoothness_orders == 1
    @test tmle.Q̅.HALClassifier_1.cv_select == false
    @test tmle.Q̅.HALClassifier_1.num_knots == [10, 5]

    # Test binary target TMLE Qstack
    tmle = tmles["continuous"]
    @test tmle.F isa GLMRegressor
    ## Checking Qstack.metalearner
    @test tmle.Q̅.metalearner isa LinearRegressor
    @test tmle.Q̅.metalearner.fit_intercept == false

    ## Checking Qstack.resampling
    @test tmle.Q̅.resampling isa CV
    @test tmle.Q̅.resampling.nfolds == 2
    ## Checking Qstack XGBoost models
    @test tmle.Q̅.XGBoostRegressor_1.num_round == 10
    @test tmle.Q̅.XGBoostRegressor_2.num_round == 20
    ## Checking Qstack Interaction Linear model
    @test tmle.Q̅.InteractionLMRegressor_1.interaction_transformer.column_pattern == r"^RS_"
    ## Checking Qstack HAL model
    @test tmle.Q̅.HALRegressor_1.lambda == 10
    @test tmle.Q̅.HALRegressor_1.smoothness_orders == 1
    @test tmle.Q̅.HALRegressor_1.cv_select == false
    @test tmle.Q̅.HALRegressor_1.num_knots == [10, 5]

    # Both TMLE have the same G Stack
    for (type, tmle) in tmles
        @test tmle.queries == (
            (RSID_10 = ["AG", "GG"], RSID_100 = ["AG", "GG"]),
            (RSID_10 = ["AG", "GG"], RSID_100 = ["AA", "GG"])
            )
        # Checking Gstack
        @test tmle.G isa FullCategoricalJoint
        ## Checking Gstack.metalearner
        @test tmle.G.model.metalearner isa LogisticClassifier
        @test tmle.G.model.metalearner.fit_intercept == false
        ## Checking Gstack.resampling
        @test tmle.G.model.resampling isa StratifiedCV
        @test tmle.G.model.resampling.nfolds == 2
        ## Checking Gstack models
        @test tmle.G.model.LogisticClassifier_1.lambda == 1.0
        @test tmle.G.model.XGBoostClassifier_1.num_round == 10
    end
    rm(queryfile)
end

@testset "Test standard ATE TMLE build" begin
    tmle_config = joinpath("config", "tmle_config.toml")
    build_ate_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    tmles =  TMLEEpistasis.estimators_from_toml(TOML.parsefile(tmle_config), queries, TMLEEpistasis.PhenotypeTMLEEpistasis)
    for (type, tmle) in tmles
        @test tmle.queries == (
            (RSID_10 = ["AG", "GG"],),
            (RSID_10 = ["AA", "GG"],)
            )
        # Checking Gstack
        @test tmle.G isa Stack
    end
    rm(queryfile)
end

@testset "Test cross validation estimators build" begin
    tmle_config = joinpath("config", "tmle_config.toml")
    build_query_file()
    queries = TMLEEpistasis.parse_queries(queryfile)
    libraries =  TMLEEpistasis.estimators_from_toml(TOML.parsefile(tmle_config), queries, TMLEEpistasis.PhenotypeCrossValidation)
    # G settings
    G_settings = libraries["G"]
    @test G_settings[1] isa StratifiedCV
    @test G_settings[1].nfolds == 2
    models =
    @test G_settings[2][:XGBoostClassifier_1].num_round == 10
    @test G_settings[2][:LogisticClassifier_1].fit_intercept == true

    # Qcont settings
    Qcont_settings = libraries["Qcont"]
    @test Qcont_settings[1] isa CV
    @test Qcont_settings[1].nfolds == 2
    @test Qcont_settings[2][:HALRegressor_1].lambda == 10
    @test Qcont_settings[2][:XGBoostRegressor_1].num_round == 10
    @test Qcont_settings[2][:XGBoostRegressor_2].num_round == 20
    @test Qcont_settings[2][:InteractionLMRegressor_1].interaction_transformer.column_pattern == r"^RS_"

    # Qcat settings
    Qcat_settings = libraries["Qcat"]
    @test Qcat_settings[1] isa StratifiedCV
    @test Qcat_settings[1].nfolds == 2
    @test Qcat_settings[2][:HALClassifier_1].lambda == 10
    @test Qcat_settings[2][:XGBoostClassifier_1].num_round == 10
    @test Qcat_settings[2][:InteractionLMClassifier_1].interaction_transformer.column_pattern == r"^RS_"

    rm(queryfile)
end

end;

true


text:
######################################################################
# circle.jl
######################################################################

export circle

"""
`circle(d::Domain, xc::VF, r::Float64)`

Returns a distribution that evenly weights points in a circle centered at `xc` with radius `r`. The distribution is normalized (so it integrates to 1) before being returned to the user.
"""
function circle(d::Domain, xc::Vector{Float64}, radius::Float64)
	phi = zeros(x_cells(d), y_cells(d))
	R2 = radius * radius

	for xi = 1:x_cells(d)
		x = x_min(d) + (xi-0.5)*x_size(d)
		for yi = 1:y_cells(d)
			y = y_min(d) + (yi-0.5)*y_size(d)
			dx = xc[1] - x
			dy = xc[2] - y
			if (dx*dx + dy*dy) < R2
				phi[xi,yi] = 1.0
			end
		end
	end
	normalize!(phi, d.cell_size)
	return phi
end


text:
# Julia wrapper for header: /home/s/juliastuff/tgl/tgl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function TGLC_bn_ctx_new()
    ccall((:TGLC_bn_ctx_new,libtelegram),Ptr{TGLC_bn_ctx},())
end

function TGLC_bn_ctx_free(ctx)
    ccall((:TGLC_bn_ctx_free,libtelegram),Void,(Ptr{TGLC_bn_ctx},),ctx)
end

function TGLC_bn_new()
    ccall((:TGLC_bn_new,libtelegram),Ptr{TGLC_bn},())
end

function TGLC_bn_free(a)
    ccall((:TGLC_bn_free,libtelegram),Void,(Ptr{TGLC_bn},),a)
end

function TGLC_bn_clear_free(a)
    ccall((:TGLC_bn_clear_free,libtelegram),Void,(Ptr{TGLC_bn},),a)
end

function TGLC_bn_cmp(a,b)
    ccall((:TGLC_bn_cmp,libtelegram),Cint,(Ptr{TGLC_bn},Ptr{TGLC_bn}),a,b)
end

function TGLC_bn_is_prime(a,checks::Cint,callback,ctx,cb_arg)
    ccall((:TGLC_bn_is_prime,libtelegram),Cint,(Ptr{TGLC_bn},Cint,Ptr{Void},Ptr{TGLC_bn_ctx},Ptr{Void}),a,checks,callback,ctx,cb_arg)
end

function TGLC_bn_bn2bin(a,to)
    ccall((:TGLC_bn_bn2bin,libtelegram),Cint,(Ptr{TGLC_bn},Ptr{Cuchar}),a,to)
end

function TGLC_bn_bin2bn(s,len::Cint,ret)
    ccall((:TGLC_bn_bin2bn,libtelegram),Ptr{TGLC_bn},(Ptr{Cuchar},Cint,Ptr{TGLC_bn}),s,len,ret)
end

function TGLC_bn_set_word(a,w::Culong)
    ccall((:TGLC_bn_set_word,libtelegram),Cint,(Ptr{TGLC_bn},Culong),a,w)
end

function TGLC_bn_get_word(a)
    ccall((:TGLC_bn_get_word,libtelegram),Culong,(Ptr{TGLC_bn},),a)
end

function TGLC_bn_num_bits(a)
    ccall((:TGLC_bn_num_bits,libtelegram),Cint,(Ptr{TGLC_bn},),a)
end

function TGLC_bn_sub(r,a,b)
    ccall((:TGLC_bn_sub,libtelegram),Void,(Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn}),r,a,b)
end

function TGLC_bn_div(dv,rem,a,d,ctx)
    ccall((:TGLC_bn_div,libtelegram),Cint,(Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn_ctx}),dv,rem,a,d,ctx)
end

function TGLC_bn_mod_exp(r,a,p,m,ctx)
    ccall((:TGLC_bn_mod_exp,libtelegram),Cint,(Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn},Ptr{TGLC_bn_ctx}),r,a,p,m,ctx)
end

function tgl_reopen_binlog_for_writing(TLS)
    ccall((:tgl_reopen_binlog_for_writing,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_replay_log(TLS)
    ccall((:tgl_replay_log,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_print_stat(TLS,s,len::Cint)
    ccall((:tgl_print_stat,libtelegram),Cint,(Ptr{tgl_state},Cstring,Cint),TLS,s,len)
end

function tgl_peer_get(TLS,id::tgl_peer_id_t)
    ccall((:tgl_peer_get,libtelegram),Ptr{tgl_peer_t},(Ptr{tgl_state},tgl_peer_id_t),TLS,id)
end

function tgl_peer_get_by_name(TLS,s)
    ccall((:tgl_peer_get_by_name,libtelegram),Ptr{tgl_peer_t},(Ptr{tgl_state},Cstring),TLS,s)
end

function tgl_message_get(TLS,id)
    ccall((:tgl_message_get,libtelegram),Ptr{tgl_message},(Ptr{tgl_state},Ptr{tgl_message_id_t}),TLS,id)
end

function tgl_peer_iterator_ex(TLS,it,extra)
    ccall((:tgl_peer_iterator_ex,libtelegram),Void,(Ptr{tgl_state},Ptr{Void},Ptr{Void}),TLS,it,extra)
end

function tgl_complete_user_list(TLS,index::Cint,text,len::Cint,R)
    ccall((:tgl_complete_user_list,libtelegram),Cint,(Ptr{tgl_state},Cint,Cstring,Cint,Ptr{Cstring}),TLS,index,text,len,R)
end

function tgl_complete_chat_list(TLS,index::Cint,text,len::Cint,R)
    ccall((:tgl_complete_chat_list,libtelegram),Cint,(Ptr{tgl_state},Cint,Cstring,Cint,Ptr{Cstring}),TLS,index,text,len,R)
end

function tgl_complete_encr_chat_list(TLS,index::Cint,text,len::Cint,R)
    ccall((:tgl_complete_encr_chat_list,libtelegram),Cint,(Ptr{tgl_state},Cint,Cstring,Cint,Ptr{Cstring}),TLS,index,text,len,R)
end

function tgl_complete_peer_list(TLS,index::Cint,text,len::Cint,R)
    ccall((:tgl_complete_peer_list,libtelegram),Cint,(Ptr{tgl_state},Cint,Cstring,Cint,Ptr{Cstring}),TLS,index,text,len,R)
end

function tgl_complete_channel_list(TLS,index::Cint,text,len::Cint,R)
    ccall((:tgl_complete_channel_list,libtelegram),Cint,(Ptr{tgl_state},Cint,Cstring,Cint,Ptr{Cstring}),TLS,index,text,len,R)
end

function tgl_secret_chat_for_user(TLS,user_id::tgl_peer_id_t)
    ccall((:tgl_secret_chat_for_user,libtelegram),Cint,(Ptr{tgl_state},tgl_peer_id_t),TLS,user_id)
end

function tgl_do_send_bot_auth(TLS,code,code_len::Cint,callback,callback_extra)
    ccall((:tgl_do_send_bot_auth,libtelegram),Cint,(Ptr{tgl_state},Cstring,Cint,Ptr{Void},Ptr{Void}),TLS,code,code_len,callback,callback_extra)
end

function tgl_set_binlog_mode(TLS,mode::Cint)
    ccall((:tgl_set_binlog_mode,libtelegram),Void,(Ptr{tgl_state},Cint),TLS,mode)
end

function tgl_set_binlog_path(TLS,path)
    ccall((:tgl_set_binlog_path,libtelegram),Void,(Ptr{tgl_state},Cstring),TLS,path)
end

function tgl_set_auth_file_path(TLS,path)
    ccall((:tgl_set_auth_file_path,libtelegram),Void,(Ptr{tgl_state},Cstring),TLS,path)
end

function tgl_set_download_directory(TLS,path)
    ccall((:tgl_set_download_directory,libtelegram),Void,(Ptr{tgl_state},Cstring),TLS,path)
end

function tgl_set_callback(TLS,cb)
    ccall((:tgl_set_callback,libtelegram),Void,(Ptr{tgl_state},Ptr{tgl_update_callback}),TLS,cb)
end

function tgl_set_rsa_key(TLS,key)
    ccall((:tgl_set_rsa_key,libtelegram),Void,(Ptr{tgl_state},Cstring),TLS,key)
end

function tgl_set_rsa_key_direct(TLS,e::Culong,n_bytes::Cint,n)
    ccall((:tgl_set_rsa_key_direct,libtelegram),Void,(Ptr{tgl_state},Culong,Cint,Ptr{Cuchar}),TLS,e,n_bytes,n)
end

function tgl_set_app_version(TLS,app_version)
    ccall((:tgl_set_app_version,libtelegram),Void,(Ptr{tgl_state},Cstring),TLS,app_version)
end

function tgl_get_peer_type(id::tgl_peer_id_t)
    ccall((:tgl_get_peer_type,libtelegram),Cint,(tgl_peer_id_t,),id)
end

function tgl_get_peer_id(id::tgl_peer_id_t)
    ccall((:tgl_get_peer_id,libtelegram),Cint,(tgl_peer_id_t,),id)
end

function tgl_set_peer_id(_type::Cint,id::Cint)
    ccall((:tgl_set_peer_id,libtelegram),tgl_peer_id_t,(Cint,Cint),_type,id)
end

function tgl_cmp_peer_id(a::tgl_peer_id_t,b::tgl_peer_id_t)
    ccall((:tgl_cmp_peer_id,libtelegram),Cint,(tgl_peer_id_t,tgl_peer_id_t),a,b)
end

function tgl_incr_verbosity(TLS)
    ccall((:tgl_incr_verbosity,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_set_verbosity(TLS,val::Cint)
    ccall((:tgl_set_verbosity,libtelegram),Void,(Ptr{tgl_state},Cint),TLS,val)
end

function tgl_enable_pfs(TLS)
    ccall((:tgl_enable_pfs,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_set_test_mode(TLS)
    ccall((:tgl_set_test_mode,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_set_net_methods(TLS,methods)
    ccall((:tgl_set_net_methods,libtelegram),Void,(Ptr{tgl_state},Ptr{tgl_net_methods}),TLS,methods)
end

function tgl_set_timer_methods(TLS,methods)
    ccall((:tgl_set_timer_methods,libtelegram),Void,(Ptr{tgl_state},Ptr{tgl_timer_methods}),TLS,methods)
end

function tgl_set_ev_base(TLS,ev_base)
    ccall((:tgl_set_ev_base,libtelegram),Void,(Ptr{tgl_state},Ptr{Void}),TLS,ev_base)
end

function tgl_authorized_dc(TLS,DC)
    ccall((:tgl_authorized_dc,libtelegram),Cint,(Ptr{tgl_state},Ptr{tgl_dc}),TLS,DC)
end

function tgl_signed_dc(TLS,DC)
    ccall((:tgl_signed_dc,libtelegram),Cint,(Ptr{tgl_state},Ptr{tgl_dc}),TLS,DC)
end

function tgl_init(TLS)
    ccall((:tgl_init,libtelegram),Cint,(Ptr{tgl_state},),TLS)
end

function tgl_dc_authorize(TLS,DC)
    ccall((:tgl_dc_authorize,libtelegram),Void,(Ptr{tgl_state},Ptr{tgl_dc}),TLS,DC)
end

function tgl_dc_iterator(TLS,iterator)
    ccall((:tgl_dc_iterator,libtelegram),Void,(Ptr{tgl_state},Ptr{Void}),TLS,iterator)
end

function tgl_dc_iterator_ex(TLS,iterator,extra)
    ccall((:tgl_dc_iterator_ex,libtelegram),Void,(Ptr{tgl_state},Ptr{Void},Ptr{Void}),TLS,iterator,extra)
end

function tgl_insert_empty_user(TLS,id::Cint)
    ccall((:tgl_insert_empty_user,libtelegram),Void,(Ptr{tgl_state},Cint),TLS,id)
end

function tgl_insert_empty_chat(TLS,id::Cint)
    ccall((:tgl_insert_empty_chat,libtelegram),Void,(Ptr{tgl_state},Cint),TLS,id)
end

function tgl_free_all(TLS)
    ccall((:tgl_free_all,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_register_app_id(TLS,app_id::Cint,app_hash)
    ccall((:tgl_register_app_id,libtelegram),Void,(Ptr{tgl_state},Cint,Cstring),TLS,app_id,app_hash)
end

function tgl_login(TLS)
    ccall((:tgl_login,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_enable_ipv6(TLS)
    ccall((:tgl_enable_ipv6,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_enable_bot(TLS)
    ccall((:tgl_enable_bot,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_state_alloc()
    ccall((:tgl_state_alloc,libtelegram),Ptr{tgl_state},())
end

function tgl_disable_link_preview(TLS)
    ccall((:tgl_disable_link_preview,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_do_lookup_state(TLS)
    ccall((:tgl_do_lookup_state,libtelegram),Void,(Ptr{tgl_state},),TLS)
end

function tgl_get_allocated_bytes()
    ccall((:tgl_get_allocated_bytes,libtelegram),Clonglong,())
end


text:
# -----------------------------------------------------------------------
# Movement optimization
# -----------------------------------------------------------------------
# OBS.:
# These algorithms were not during the tests!
# It was obseved little impact on the experiments runtime.
# If necessary the user can adapt and use it to determine moving order
# -----------------------------------------------------------------------
# Compute distances
# -----------------------------------------------------------------------
using Distances

"""
    distcylinder(posOld::CylType, posNew::CylType; metric="steps")

## Description
Compute distance between two positions.
## Arguments
- posOld: previous position (from)
- posNew: new position (to)
- metric: metric to compute distance: "euclidean", "manhattan", "steps"
# Examples
```jldoctest
julia> distcylinder(OneCylType(0,0), OneCylType(10,0))
57

julia> distcylinder(OneCylType(0,0), OneCylType(10,10), metric="manhattan")
20.0

julia> distcylinder(ThreeCylType(0,0,5,0,5,0), ThreeCylType(10,5,5,0,5,1))
696
```
"""
function distcylinder(posOld::OneCylType, posNew::OneCylType; metric="steps")
    if metric == "euclidean"
        r = euclidean([posNew.x,posNew.y], [posOld.x,posOld.y])
    elseif metric == "manhattan"
        r = abs(posNew.x-posOld.x) + abs(posNew.y-posOld.y)
    elseif metric == "steps"
        steps = coord2steps(posOld, posNew)[1].step
        r = sum(abs.(steps))
    else
        error("Metric is not valid")
    end
    return r
end

function distcylinder(posOld::TwoCylType, posNew::TwoCylType; metric="steps")
    if metric == "euclidean"
        r1 = euclidean([posNew.x,posNew.y], [posOld.x,posOld.y])
        r2 = euclidean([posNew.x+posNew.dx, posNew.y+posNew.dy],
                        [posOld.x+posOld.dx, posOld.y+posOld.dy])
        r = r1 + r2
    elseif metric == "manhattan"
        r1 = abs(posNew.x-posOld.x) + abs(posNew.y-posOld.y)
        r2 = r1 + abs(posNew.dx-posOld.dx) + abs(posNew.dy-posOld.dy)
        r = r1 + r2
    elseif metric == "steps"
        steps = coord2steps(posOld, posNew)[1].step
        r = sum(abs.(steps))
    else
        error("Metric is not valid")
    end
    return r
end

function distcylinder(posOld::ThreeCylType, posNew::ThreeCylType; metric="steps")
    if metric == "euclidean"
        r1 = euclidean([posNew.x,posNew.y], [posOld.x,posOld.y])
        r2 = euclidean([posNew.x+posNew.dx1,posNew.y+posNew.dy1],
                        [posOld.x+posOld.dx1,posOld.y+posOld.dy1])
        r3 = euclidean([posNew.x+posNew.dx1+posNew.dx2,posNew.y+posNew.dy1+posNew.dy2],
                        [posOld.x+posOld.dx1+posOld.dx2,posOld.y+posOld.dy1+posOld.dy2])
        r = r1 + r2 + r3
    elseif metric == "manhattan"
        r1 = abs(posNew.x-posOld.x) + abs(posNew.y-posOld.y)
        r2 = r1 + abs(posNew.dx1-posOld.dx1) + abs(posNew.dy1-posOld.dy1)
        r3 = r2 + abs(posNew.dx2-posOld.dx2) + abs(posNew.dy2-posOld.dy2)
        r = r1 + r2 + r3
    elseif metric == "steps"
        steps = coord2steps(posOld, posNew)[1].step
        r = sum(abs.(steps))
    end
    return r
end

"""
    greedypath(actualPos::CylType, points::ArrayCylType)

## Description
Greedy algorithm to determine the near optimum path (nearest neighbors).
## Arguments
- actualPos: the initial position of the robot
- points: positions the robot must visit
"""
function greedypath(actualPos::CylType, points::ArrayCylType)

    if typeof(actualPos) != typeof(points[1])
        error("actualPos and points don't have elements with the same type")
    end

    path = [] # initialize path

    n = length(points)
    notVisited = repeat([true],n)
    id = collect(1:n)

    while sum(notVisited) > 0
        minDist = argmin([distcylinder(actualPos, i) for i in points[notVisited]])
        idtemp = id[notVisited][minDist]
        path = push!(path, idtemp)
        notVisited[idtemp] = false
        actualPos = points[idtemp]
    end

    return path
end

"""
    benchoptimizer(n=50; nc=3, optimFun=greedypath)

## Description
Benchmark tool for the path optimizer (compared to random ordered).
## Arguments
- n: number of points
- nc: number of cylinders
- fun: optimization function (only greeypath allowed yet)
"""
function benchoptimizer(n=50; nc=3, fun=greedypath)

    # One cylinder
    initValue = randomInd(nc)
    p = initPop(n,nc=nc)
    s = distcylinder(initValue, p[1])

    id_optim = fun(initValue,p)
    s_optim = distcylinder(initValue, p[id_optim[1]])

    for i in 2:n
        k = i-1
        s += distcylinder(p[k], p[i])
        s_optim += distcylinder(p[id_optim[k]], p[id_optim[i]])
    end

    println("Random order: ", s, " steps")
    println("Optimized: ", s_optim, " steps")
    println("Ratio: ", s_optim/s)

    return s_optim/s
end