text:
# Model
include_model("hopper")

# modify parameters to match implicit dynamics example
gravity = 9.81 # gravity
μ_world = 0.8 # coefficient of friction
μ_joint = 0.0

mb = 3.0 # body mass
ml = 0.3  # leg mass
Jb = 0.75 # body inertia
Jl = 0.075 # leg inertia

model = Hopper{Discrete, FixedTime}(n, m, d,
			   mb, ml, Jb, Jl,
			   μ_world, gravity,
			   qL, qU,
			   uL, uU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

# Horizon
T = 21
# Time step
h = 0.05
tf = h * (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 10.0 #model_ft.uU
_ul = zeros(model.m)
_ul[model.idx_u] .= -10.0 #model_ft.uL
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0; 0.5; 0.0; 0.5]
# qM = [0.5; 0.5; 0.0; 0.5]
qT = [1.0; 0.5; 0.0; 0.5]
q_ref = [0.5; 0.75; 0.0; 0.25]
x_ref = [q_ref; q_ref]
ql = [-Inf; 0; -Inf; 0.0]
qu = [Inf; Inf; Inf; 1.0]
xl, xu = state_bounds(model, T, [ql; ql], [qu; qu])
xl[1][nq .+ (1:nq)] = q1
xu[1][nq .+ (1:nq)] = q1

# Objective

# # gait 1
Q = [(t == 1 ? 1.0 * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0])
	: t == T ? Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0])
	: 0.1 * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0])) for t = 1:T]
q = [2.0 * Q[t] * x_ref for t = 1:T]
R = [0.1 * Diagonal(ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

# gait 2
Q = [(t == 1 ? 1.0 * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0])
	: t == T ? Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0])
	: 1.0 * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0])) for t = 1:T]
q = [2.0 * Q[t] * x_ref for t = 1:T]
R = [1.0 * Diagonal(ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

# gait 3
Q = [(t == 1 ? 1.0 * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0])
	: t == T ? Diagonal([1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0; 1.0])
	: 0.1 * Diagonal([1.0; 10.0; 1.0; 10.0; 1.0; 10.0; 1.0; 10.0])) for t = 1:T]
q = [2.0 * Q[t] * x_ref for t = 1:T]
R = [1.0e-3 * Diagonal(ones(model.m)) for t = 1:T-1]
r = [zeros(model.m) for t = 1:T-1]

obj_tracking = quadratic_tracking_objective(Q, R, q, r)

obj_contact_penalty = PenaltyObjective(1.0e4, model.m)

obj = MultiObjective([obj_tracking, obj_contact_penalty])

# Constraints
include_constraints(["contact", "stage", "loop"])

function pinned!(c, x, u, t)
    _q1 = view(x, 1:4)
    _q2 = view(x, 4 .+ (1:4))
    c[1:2] = kinematics(model, _q1) - kinematics(model, q1)
    c[3:4] = kinematics(model, _q2) - kinematics(model, q1)
    nothing
end
n_pinned = 4
con_pinned = stage_constraints(pinned!, n_pinned, (1:0), collect(1:1))

function distance_traveled!(c, x, u, t)
    x_travel = 0.5
    _q1 = view(x, 1:4)
    _q2 = view(x, 4 .+ (1:4))
	c[1] = x[1] - x_travel #TODO: add methods for trajectory-wise constraints in order to change to exactly match implicit dynamics example
    c[2] = x[5] - x_travel
    nothing
end

n_distance = 2
con_distance = stage_constraints(distance_traveled!, n_distance, (1:2), collect(T:T))

con_contact = contact_constraints(model, T)

con_loop = loop_constraints(model, collect([(2:model.nq)...,
	(nq .+ (2:model.nq))...]), 1, T)

con = multiple_constraints([con_contact, con_pinned, con_distance, con_loop])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = [[q1; q1] for t = 1:T]
@load joinpath(pwd(), "examples/implicit_dynamics/examples/comparisons/hopper_stand.jld2") u_stand
u0 = [u_stand[1] for t = 1:T-1]

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
# include_snopt()

@time z̄, info = solve(prob, copy(z0),
	nlp = :ipopt,
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

x̄, ū = unpack(z̄, prob)
@show check_slack(z̄, prob)

# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/hopper_gait_1_direct.jld2") x̄ ū
# @load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/hopper_gait_1_direct.jld2") x̄ ū

# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/hopper_gait_2_direct.jld2") x̄ ū
# @load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/hopper_gait_2_direct.jld2") x̄ ū

# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/hopper_gait_3_direct.jld2") x̄ ū
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/hopper_gait_3_direct.jld2") x̄ ū

# compute comparable objective
J = 0.0
x = x̄
u = ū
for t = 1:T
    if t < T
        J += x[t]' * Q[t] * x[t] + q[t]' * x[t] + u[t][1:2]' * R[t][1:2, 1:2] * u[t][1:2] + r[t][1:2]' * u[t][1:2]
    elseif t == T
        J += x[t]' * Q[t] * x[t] + q[t]' * x[t]
    else
        J += 0.0
    end
end
@show J

plot(hcat(ū...)[1:2, :]', linetype = :steppost)

using Plots
t = range(0, stop = h * (T - 1), length = T)
plot(t, hcat(ū..., ū[end])[1:2,:]', linetype=:steppost,
	xlabel="time (s)", ylabel = "control",
	label = ["angle" "length"],
	width = 2.0, legend = :top)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model,
    state_to_configuration(x̄),
	# state_to_configuration([[x̄[1] for i = 1:50]...,x̄..., [x̄[end] for i = 1:50]...]),
	Δt = h,
	scenario = :flip)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
struct HybridQuerySynthesisPQs <: HybridPQs
    qss::QuerySynthesisStrategy
    dist_func::Function
    function HybridQuerySynthesisPQs(model::SVDD.OCClassifier, data::Array{T, 2}; qss_type=nothing, dist_func=:euclidean, other_params...) where T <: Real
        new(initialize_qs(eval(qss_type), model, data, Dict{Symbol, Any}(other_params)), eval(dist_func))
    end
end

"""
get_query_object(qs::HybridQuerySynthesisPQs, data::Array{T, 2}, pools::Vector{Symbol}, global_indices::Vector{Int}, history::Vector{Int})
This query strategy first performs a query synthesis and then selects the unlabeled
observation closest to the synthetic query as final query.
# Arguments
- `query_data`: Subset of the full data set
- `pools`: Labels for `query_data`
- `global_indices`: Indices of the observations in `query_data` relative to the full data set.
- `history`: Indices of previous queries
"""
function get_query_objects(qs::HybridQuerySynthesisPQs,
                          query_data::Array{T, 2},
                          pools::Vector{Symbol},
                          global_indices::Vector{Int},
                          history::Vector{Vector{Int}})::Vector{Int} where T <: Real
    pool_map = MLLabelUtils.labelmap(pools)
    haskey(pool_map, :U) || throw(ArgumentError("No more points that are unlabeled."))
    # Perform query synthesis
    mask = pools .!= :U
    all_history_values = (length(history) > 0) && (collect(Iterators.flatten(history)))
    query_history = [query_data[:, i:i] for i in indexin(all_history_values, global_indices)]
    qss_object = vec(get_query_objects(qs.qss, query_data[:, mask], pools[mask], query_history))
    @assert length(qss_object) == size(query_data, 1)
    # Choose unlabeled observation closest to synthetic query
    candidates = [i for i in pool_map[:U] if global_indices[i] ∉ history]
    distances = [qs.dist_func(qss_object, query_data[:, i]) for i in candidates]
    local_query_index = candidates[argmin(distances)]
    return [global_indices[local_query_index]]
end

ext: jl


text: Dict{String,Any}("albums" => Any[Dict{String,Any}("name" => Dict{String,Any}("value" => "Born to Run","type" => "string"),"songs" => Any[Dict{String,Any}("name" => Dict{String,Any}("value" => "Jungleland","type" => "string")), Dict{String,Any}("name" => Dict{String,Any}("value" => "Meeting Across the River","type" => "string"))]), Dict{String,Any}("name" => Dict{String,Any}("value" => "Born in the USA","type" => "string"),"songs" => Any[Dict{String,Any}("name" => Dict{String,Any}("value" => "Glory Days","type" => "string")), Dict{String,Any}("name" => Dict{String,Any}("value" => "Dancing in the Dark","type" => "string"))])])

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
# load solvers
using HALeqO
using NLPModelsIpopt
using NCL # from https://github.com/JuliaSmoothOptimizers/NCL.jl, using IPOPT
using Percival

# load problems
using CUTEst
probnames = CUTEst.select(
    min_var = 1, # default 1
    min_con = 1, # default 1
    max_var = 100, # default 100
    max_con = 100, # default 100
    only_free_var = true,
    only_equ_con = true,
)

# load tools
using SolverBenchmark
using CSV, DataFrames, Printf
using Plots
include("time_profiles.jl")

# setup benchmarking
filename = "cutest_eq_tol6"
TOL = 1e-6 # default 1e-6
MAXITER = 3000 # default 3000

problems = (CUTEstModel(probname) for probname in probnames)
solvers = Dict{Symbol,Function}(
    :HALeqO => prob -> haleqo(prob; tol = TOL, max_iter = MAXITER),
    :NCL => prob -> NCLSolve(prob, opt_tol = TOL, feas_tol = TOL, max_iter = MAXITER),
    :IPOPT => prob -> ipopt(prob; tol = TOL, print_level = 0, max_iter = MAXITER),
    :Percival => prob -> percival(prob; atol = TOL, rtol = 0.0, ctol = TOL, max_iter = MAXITER),
)

# warm up
nlp = CUTEstModel("BT3")
for solver ∈ keys(solvers)
    out = solvers[solver](nlp)
end
finalize(nlp)

# run solvers!
stats = bmark_solvers(solvers, problems)

# get statistics
for solver ∈ keys(stats)
    @info "$solver statuses" count_unique(stats[solver].status)
end

statuses, avgs = quick_summary(stats; cols = [:iter, :elapsed_time])
for solver ∈ keys(stats)
    @info "statistics for" solver statuses[solver] avgs[solver]
end

num_solved_problems = Dict{Symbol,Int}()
for solver ∈ keys(stats)
    num_solved_problems[solver] = sum(stats[solver].status .== :first_order)
end

# performance and time profiles
cost(df) = (df.status .!= :first_order) * Inf + df.elapsed_time
pprof = performance_profile(stats, cost)
tprof = time_profile(stats, cost)

# store data
for solver ∈ keys(stats)
    CSV.write("data/" * filename * "_" * String(solver) * ".csv", stats[solver], header = true)
end


text:
# run all tests

tests = [ "shapes",
          "diagop",
          "mathfuns",
          "functors",
          "map",
          "vbroadcast",
          "transforms",
          "utils",
          "rkernels",
          "reduce",
          "folddim",
          "reducedim",
          "norms",
          "scan",
          "statistics",
          "wsum"]

srand(6789)

for t in tests
    tf = joinpath("test", "$t.jl")
    println("Running $tf ...")
    include(tf)
end


text:
# This file contains some extra checks that don't fall into the other
# categories, for example because they're very platform-specific.

function check_os_abi(oh::ObjectHandle, p::AbstractPlatform, rest...; verbose::Bool = false, kwargs...)
    if Sys.isfreebsd(p)
        if oh.ei.osabi != ELF.ELFOSABI_FREEBSD
            # The dynamic loader should not have problems in this case, but the
            # linker may not appreciate.  Let the user know about this.
            if verbose
                msg = replace("""
                $(basename(path(oh))) has an ELF header OS/ABI value that is not set to FreeBSD
                ($(ELF.ELFOSABI_FREEBSD)), this may be an issue at link time
                """, '\n' => ' ')
                @warn(strip(msg))
            end
            return false
        end
    elseif call_abi(p) == "eabihf"
        # Make sure the object file has the hard-float ABI.  See Table 4-2 of
        # "ELF for the ARM Architecture" document
        # (https://developer.arm.com/documentation/ihi0044/e/).  Note: `0x000`
        # means "no specific float ABI", `0x400` == EF_ARM_ABI_FLOAT_HARD.
        if header(oh).e_flags & 0xF00 ∉ (0x000, 0x400)
            if verbose
                @error("$(basename(path(oh))) does not match the hard-float ABI")
            end
            return false
        end
    end
    return true
end


text:
"""
    TopParser.samples(content::Union{IO,AbstractString}) -> samplesitr
"""
TopParser.samples
TopParser.samples(content::AbstractString) = SamplesIterator(IOBuffer(content))
TopParser.samples(io::IO) = SamplesIterator(io)

"""
    TopParser.open_samples(f, path)
"""
TopParser.open_samples
function TopParser.open_samples(f, path)
    open(path) do io
        f(TopParser.samples(io))
    end
end

"""
    TopParser.processes(content::Union{IO,AbstractString})
    TopParser.processes(samplesitr)
"""
TopParser.processes
TopParser.processes(content::Union{IO,AbstractString}) =
    TopParser.processes(TopParser.samples(content))
TopParser.processes(samples::SamplesIterator) =
    Iterators.flatten(Iterators.map(get_processes, samples))

get_processes(sample::Sample) = sample.processes

"""
    TopParser.processes_from_file([pred,] path) -> processes::Vector{Process}
"""
TopParser.processes_from_file
TopParser.processes_from_file(path) = TopParser.processes_from_file(always, path)
function TopParser.processes_from_file(pred, path)
    TopParser.open_samples(path) do samples
        collect(Iterators.filter(pred, TopParser.processes(samples)))
    end
end

always(_) = true


text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Run VI for the Sleipner model

## Module loading

using DrWatson
@quickactivate "UncertaintyQuantificationAVP"
using UncertaintyQuantificationAVP
using Flux.Optimise, InvertibleNetworks
using LinearAlgebra, PyPlot, Statistics
using Random; seed = 1; Random.seed!(seed)

## Loading synthetic data and inversion parameters

# Data
data_dict = wload(datadir("Sleipner", "Sleipner_data.bson"))
dat_n = data_dict[:dat_n]
src_pars = data_dict[:src_pars]
abc_geom = data_dict[:abc_geom]
acc = data_dict[:acc_mod]

# Model
model_dict = wload(projectdir()*"/scripts/Sleipner/Sleipner_model.bson")
Mtrue = model_dict[:Mtrue]
Mbg = model_dict[:Mbg]

# Inversion parameters
invpars_dict = wload(datadir("Sleipner", "Sleipner_invpars.bson"))
Mpr = invpars_dict[:Mpr]
M0  = invpars_dict[:M0]
σ2_n = invpars_dict[:σ2_n]
ϵ2 = invpars_dict[:ϵ2]
σ2_pr = invpars_dict[:σ2_pr]
mmin = invpars_dict[:mmin]
mmax = invpars_dict[:mmax]

# MAP result
Mmap = wload(datadir("Sleipner", "Sleipner_MAP.bson"))[:Mmap]

## Negative log density for posterior distribution

# Pre- and post-processing for log-posterior
nz = Mbg.geom.nz
dz = Mbg.geom.dz
function proj_bounds(m::Array{Float32}, mmin, mmax)
    m_ = copy(m)
    m_[m_ .< mmin] .= mmin
    m_[m_ .> mmax] .= mmax
    return m_
end
fact_tune = 1f-8
prec(Δm::Array{Float32, 4}) = Model(Mbg.geom, proj_bounds(Mbg.m.+fact_tune*reshape(Δm, nz, :), mmin, mmax))
# prec(Δm::Array{Float32, 4}) = Model(Mbg.geom, proj_bounds(Mmap.m.+fact_tune*reshape(Δm, nz, :), mmin, mmax))
post(g::Model) = fact_tune*reshape(g.m, (1, 1, nz, :))

# Log-likelihood
F = seismod_fun(dat_n.pars, data_dict[:src_pars]; abc_geom = data_dict[:abc_geom], acc = data_dict[:acc_mod])
negloglike_fun(m::Model; compute_grad::Bool = true) = negLogLikelihood_gauss(dat_n, m, σ2_n, F; compute_grad = compute_grad)

# Log-prior
D = derivative_linop(Mbg.geom)
# neglogpr_fun(m::Model; compute_grad::Bool = true) = negLogPrior_L1smooth(m, σ2_pr, ϵ2; m0 = Mpr, D = D, compute_grad = compute_grad)
neglogpr_fun(m::Model; compute_grad::Bool = true) = negLogPrior_gauss(m, σ2_pr; m0 = Mpr, D = D, compute_grad = compute_grad)

# Log-posterior
neglogpost_fun(m::Model) = negLogPost(m, negloglike_fun, neglogpr_fun; compute_grad = true)

## Setting loss function

function loss(Z::Array{Float32, 4})
    Δm, logdet = T.forward(Z)
    m = prec(Δm)
    nlogp, ∇m = neglogpost_fun(m)
    g = post(∇m)/batchsize
    f = sum(nlogp)/batchsize-logdet
    T.backward(g, Δm)
    return f
end

## Network initialization

# Generator
depth = 5
n_hidden = Mbg.geom.nz
batchsize = 2^3
net_type = NetworkHINT1D
T = net_type(nz, batchsize, n_hidden, depth; permute="full")
# net_type = NetworkGlow1D
# T = net_type(nz, batchsize, n_hidden, depth)
initialize_id!(T)

## Training

# Training data
trainsize = 2^7
Ztrain = randn(Float32, 1, 1, nz, trainsize)

# N of iterations
nepochs = 2^11

# Optimizer
lr = 1f-4
decay = 0.9f0
nbatches = Int64(round(trainsize/batchsize))
decay_step = nbatches*2^6
clip = 0f0
opt = Optimiser(ExpDecay(lr, decay, decay_step, clip), ADAM(lr))

# Train model
fval = Array{Float32, 1}(undef, 0)
time = @elapsed training!(T, fval, loss, Ztrain, nepochs, batchsize, opt; verbose_b=true, gradclip=true)

## Testing

testsize = trainsize
Ztest = randn(Float32, 1, 1, nz, testsize)
Mtest = prec(T.forward(Ztest)[1])

## Saving results

params_exp = Dict(
    :seed => seed, #seed
    :batchsize => batchsize, :trainsize => trainsize, :Ztrain => Ztrain, :nepochs => nepochs, # training pars
    :lr => lr, :decay => decay, :decay_step => decay_step, :clip => clip, # optimizer pars
    :net_type => string(net_type), :depth => depth, :n_hidden => n_hidden, # network arch
    :fact_tune => fact_tune)
results_exp = copy(params_exp)
results_exp[:fval] = fval
results_exp[:net_params] = get_params(T)
results_exp[:Ztest] = Ztest
results_exp[:Ztrain] = Ztrain
results_exp[:Mtest] = Mtest
results_exp[:time] = time
wsave(datadir("Sleipner", savename("Sleipner_VI", params_exp, "bson")), results_exp)


text:
using Printf
using Statistics
using LinearAlgebra
using MDBM
using PyPlot
pygui(true);
using Interpolations
# using Plots

#creating iteration array
function it(A)
    matrdim=size(A,2)-1
    scaleitp=real(A[:,1])
    knots=(scaleitp,)
    ARe=real(A)
    AIm=imag(A)

    imRe=[interpolate(knots, ARe[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imRetemp=[interpolate(knots, ARe[:,j], Gridded(Linear()))]
        imRe=vcat(imRe,imRetemp)
    end
    imIm=[interpolate(knots, AIm[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imImtemp=[interpolate(knots, AIm[:,j], Gridded(Linear()))]
        imIm=vcat(imIm,imImtemp)
    end
    return(hcat(imRe,imIm))
end

#substitution in interpolating function
function sub(it,t)
    subdim=size(it,1)
    out=zeros(ComplexF64,subdim)
    for j = 1:subdim
        out[j]=it[j,1](t)+it[j,2](t)*im
    end
    return(out)
end

#function for finding the greatest eigenvalue
function normmax(a)
    s=maximum(size(a))
    norma=zeros(s)
    for j = 1:s
        norma[j]=norm(a[j])
    end
    return(a[findmax(norma)[2]])
end

#function for finding the smallest eigenvalue
function normmin(a)
    s=maximum(size(a))
    norma=zeros(s)
    for j = 1:s
        norma[j]=norm(a[j])
    end
    return(a[findmin(norma)[2]])
end

#checking the difference between the greatest and the smallest eigenvalue
function checkconv(a,b,err)
    var=false
    amax=norm(normmax(a))
    amin=norm(normmin(a))
    bmax=norm(normmax(b))
    bmin=norm(normmin(b))
    if abs((amax-amin)-(bmax-bmin))/(amax-amin)<err
        var=true
    end
    return(var)
end

#matrix exponential
function expM(M)
    p=10
    n=size(M)[1]
    Mval=zeros(Float64,n,n)
    for j = 0:p
        Mval=Mval+(M^j)/(factorial(j))
    end
    return(Mval)
end
#matrix coefficients
function A(t,vv)
    Amatr=zeros(Float64,dim,dim)
    #vv[1]=kappa #vv[2]=delta #vv[3]=epsilon #vv[4]=omega
    Amatr=[-vv[1] -(vv[2]+vv[3]*cos(vv[4]*t)); 1 0]
    return(Amatr)
end

function B(t,bv)
    Bmatr=zeros(dim*dim)
    Bmatr=[0 bv; 0 0]
    return(Bmatr)
end

###################### ISIM #######################
function ISIM(kappa1,delta1,epsilon1,b1,tau1,k1,n)
        dt=tau1/(n-1)
        v=[kappa1, delta1, epsilon1, 2*pi/(k1*tau1)]
        eval=zeros(ComplexF64,mult,1)
        S=zeros(ComplexF64,n*dim,mult)
        V=zeros(ComplexF64,n*dim,mult)
        Vjnorm=zeros(ComplexF64,n*dim,mult)
        tvect0=collect(-tau1:dt:tau1*k1+dt/2) #discretized time vector
        sol=zeros(ComplexF64,k1*(n-1),mult*dim) #empty solution vector
        sol00=rand(n,mult*dim)*(-2)+2*ones(n,mult*dim) #S0 inital random array
        solreturn=zeros(ComplexF64,1,1+mult*dim)
        for g = 1:i #iteration
        sol0=hcat(tvect0,vcat(sol00,sol)) #constructing solution array
            #time integration
            for m = 2:dim:(dim*mult-(dim-2)) #choosing "one physical system"
            sol0m=hcat(sol0[:,1],sol0[:,m:m+dim-1])
            #choosing numerical method
            if method == "EE"
                for j = 0:k1*(n-1)-1
                    Y=sol0m[n+j,2:2+(dim-1)]+dt*(A(sol0m[n+j,1],v)*sol0m[n+j,2:2+(dim-1)]+B(sol0m[n+j,1],b1)*sol0m[1+j,2:2+(dim-1)])
                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y)
                end
            elseif method == "RK4_int"
                for j = 0:k1*(n-1)-1
                    sol0mittau=it(sol0m[j+1:j+3,:])

                    Ytau0=sub(sol0mittau,-tau1+j*dt)
                    Ytau05=sub(sol0mittau,-tau1+(j+0.5)*dt)
                    Ytau1=sub(sol0mittau,-tau1+(j+1)*dt)

                    Y1=sol0m[j+n,2:dim+1]
                    kk1=(dt)*(A(j*dt,v)*Y1+B(j*dt,b1)*Ytau0)
                    kk2=(dt)*(A((j+0.5)*dt,v)*(Y1+kk1/2)+B((j+0.5)*dt,b1)*Ytau05)
                    kk3=(dt)*(A((j+0.5)*dt,v)*(Y1+kk2/2)+B((j+0.5)*dt,b1)*Ytau05)
                    kk4=(dt)*(A((j+1.0)*dt,v)*(Y1+kk3)+B((j+1.0)*dt,b1)*Ytau1)

                    Y=Y1+(kk1 + 2*kk2 + 2*kk3 + kk4)/6

                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y)
                end
            elseif method == "RK4"
                for j = 0:k1*(n-1)-1
                    Ydtau=0.5*(sol0m[j+1,2:2+(dim-1)]+sol0m[j+2,2:2+(dim-1)])

                    Y1=sol0m[n+j,2:2+(dim-1)]
                    Y2=Y1+(dt/2)*(A(j*dt,v)*Y1+B(j*dt,b1)*sol0m[j+1,2:2+(dim-1)])
                    Y3=Y1+(dt/2)*(A((j+0.5)*dt,v)*Y2+B((j+0.5)*dt,b1)*Ydtau)
                    Y4=Y1+dt*(A((j+1.0)*dt,v)*Y3+B((j+1.0)*dt,b1)*Ydtau)

                    Y=Y1+(dt/6)*((A(j*dt,v)*Y1+B(j*dt,b1)*sol0m[j+1,2:2+(dim-1)])+2*(A((j+0.5)*dt,v)*Y2+B((j+0.5)*dt,b1)*Ydtau)+2*(A((j+0.5)*dt,v)*Y3+B((j+0.5)*dt,b1)*sol0m[j+1,2:2+(dim-1)])+(A((j+1)*dt,v)*Y4+B(j*dt,b1)*sol0m[j+2,2:2+(dim-1)]))

                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y)
                end
            elseif method == "IE"
                for j = 0:k1*(n-1)-1
                    X=sol0m[j+n,2:2+(dim-1)]
                    Yt=sol0m[j+n,2:2+(dim-1)]
                    Ytau=sol0m[j+1,2:2+(dim-1)]
                    for l=1:4
                        X=X-inv(dt*A(j*dt,v)-Matrix{Float64}(I,dim,dim))*(Yt+dt*(A((j+1+n)*dt,v)*X+B((j+1+n)*dt,b1)*Ytau)-X)
                    end
                    Y=X
                   sol0m[n+j+1,2:2+(dim-1)]=transpose(Y) #filling solution array of the physical system
                end
            elseif method == "RK2"
                for j = 0:k1*(n-1)-1

                    Ydtau=0.5*(sol0m[j+1,2:2+(dim-1)]+sol0m[j+2,2:2+(dim-1)])

                    Y1=sol0m[n+j,2:2+(dim-1)]
                    Y2=Y1+(dt/2)*(A((j+0.5)*dt,v)*Y1+B((j+0.5)*dt,b1)*Ydtau)

                    Y=Y1+dt*(A((j+0.5)*dt,v)*Y2+B((j+0.5)*dt,b1)*Ydtau)

                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y)
                end
            end

                sol0[:,m:m+dim-1]=sol0m[:,2:2+(dim-1)] #filling solution array
            end
            #solreturn=vcat(solreturn,sol0)
            #retrieving the solution of last and second-to-last periods
            S1=sol0[1:n,2:mult*dim+1]
            V1=sol0[end-(n-1):end,2:mult*dim+1]
            #restructuring the solution
                for p=1:n
                    for s=1:mult
                        for q=1:dim
                            S[1+(p-1)*dim+(q-1),s]=S1[p,1+(s-1)*dim+(q-1)]
                            V[1+(p-1)*dim+(q-1),s]=V1[p,1+(s-1)*dim+(q-1)]
                        end
                    end
                end
            H=pinv(S)*V #pseudo-inverse calculation
            #H=pinv(S,rtol = sqrt(eps(real(float(one(eltype(S)))))))*V
            (evm,VjT)=eigen(H) #eigenvalue calculationV
            Vj=V*VjT #calculating of new set of eigenvectors
            #normalizing the results
            for h = 1:mult
                Vjnorm[:,h]=normalize(Vj[:,h])
            end
            sol00=zeros(ComplexF64,n,mult*dim) #creating new initial solution array
            for p=1:n
                for s=1:mult
                    for q=1:dim
                             sol00[p,1+(s-1)*dim+(q-1)]=Vjnorm[1+dim*(p-1)+(q-1),s]
                    end
                end
            end
            eval=hcat(eval,evm)     #filling up the eigenvalue array
            #checking convergence
            if g>2 && checkconv(eval[:,end],eval[:,end-1],0.01)
                 break
            end
            #if check(evm) print(g-1); break end
        end
        # print(size(eval)[2]-1)
        # return(norm(normmax(eval[:,end]))-1)
        return eval[:,end]
        #return(eval[:,2:end]) #returning the global eigenvalue array
        #return(solreturn)
end

###################### SDM ########################
function SDM(kappa1, delta1, epsilon1, b1, tau1, k1)
        dt=tau1/m #definition of time-step
        m2=trunc(Int,m/2)
        dim2=trunc(Int,dim/2)
        dimg=dim*(m2+1)
        v=[kappa1, delta1, epsilon1, 2*pi/(k1*tau1)]
        P=zeros(Float64,dim,dim*m*k1) #construction of Pi matrices
        for j=1:m*k1
            P[:,1+dim*(j-1):dim*j]=expM(A((j-1)*dt,v)*dt)
        end
        R=zeros(Float64,dim,dim*m*k1) #construction of Ri matrices
        for j=1:m*k1
            R[:,1+dim*(j-1):dim*j]=0.5*((expM(A((j-1)*dt,v)*dt)-Matrix{Float64}(I, dim, dim))*inv(A((j-1)*dt,v))*B((j-1)*dt,b1))
        end
        Zw=zeros(Float64,dimg,dimg*m*k1) #construction of Zi matrices
        for j=1:m*k1
            P1=P[:,1+(j-1)*dim:j*dim]
            R1=R[:,1+(j-1)*dim:j*dim]
            Zw[:,1+(j-1)*dim*(m2+1):j*dim*(m2+1)]=vcat(hcat(P1,zeros(dim,dim*(m2-1)),R1[:,dim2+1:dim],R1[:,dim2+1:dim]),hcat(zeros(m2*dim,dim2),Matrix{Float64}(I, m2*dim,m2*dim),zeros(m2*dim,dim2)))
        end
        #construction of final Z matrix
        Z=Zw[:,1+dimg*(m-2)+m*dimg*(k1-1):dimg*(m-1)+m*dimg*(k1-1)]*Zw[:,1+dimg*(m-1)+m*dimg*(k1-1):dimg*m+m*dimg*(k1-1)]
        for j=1:(m-2)+(k1-1)*m
            Z=Z*Zw[:,1+dimg*(m-2-j)+m*dimg*(k1-1):dimg*(m-1-j)+m*dimg*(k1-1)]
        end
        #evaluating stability based on largest eigenvalue
        return(norm(normmax(eigvals(Z)))-1)
        #return(P)
end

#parameters of the delayed Mathieu equation
dim=2 #DoF of the system (after Cauchy transcription)
const tau=2*pi
const kappa=0.1
const delta=0.5
const epsilon=1.0
const b=0.2
const k=2
#const T=k*tau
#const omega=2*pi/T
#timestep
#dt=tau/(n-1)

n=5
mus_final=ISIM(kappa, 5.0, 1.0, 2.0 , tau, k,n)
n=5000
mus_final=ISIM(kappa, 5.0, 1.0, 2.0 , tau, k,n)
figure(1)
scatter(real.(mus_final),imag.(mus_final))
plot(sin.(0:0.01:2*pi),cos.(0:0.01:2*pi))
for ni in 2:0.5:14
#
###################### ISIM #######################
mult=16 #Ns=dim*mult
# n=30 #timestep number
#EE: Explicit Euler
#RK4: 4th order Runge-Kutta
#RK4_int: 4th order Runge-Kutta with interpolation function
#IE: Implicit Euler
#RK2: 2nd order Runge-Kutta
method="RK4" #choosing numerical simulation method
i=35 #number of iteration

###################### SDM ########################

mus=ISIM(kappa, 5.0, 1.0, 2.0 , tau, k,Int64(floor(2^ni)))
# figure(1)
# scatter(real.(mus),imag.(mus))
# plot(sin.(0:0.01:2*pi),cos.(0:0.01:2*pi))
figure(2)
# scatter(ones(size(mus)).*m,abs.(mus[1]-mus_final[1]))
println(mus[end])
scatter(ones(size(mus[end])).*ni,log(abs(abs.(mus[end])-abs.(mus_final[end]))))
# scatter(ones(size(mus[end])).*ni,abs.(mus[end]))
pause(0.001)
end

#
# ####### Multi-Dimensional Bisection Method #########
#
# function foo(x,y,z)
#     return(ISIM(kappa, x, y, z , tau, k))
# end
#
#
# # ax1=Axis(-1.0:0.5:5,"delta") # initial grid in x direction
# # ax2=Axis(-1.5:0.5:1.5,"b") # initial grid in y direction
# #
# # ax1=Axis(0.5:0.2:2.0,"tau") # initial grid in x direction
# # ax2=Axis(0.05:0.2:1.5,"b") # initial grid in y direction
# ax1=Axis(0.1:1.0:15.0,"delta") # initial grid in x direction
# ax2=Axis(-5.05:0.5:5.5,"epsilon") # initial grid in y direction
# ax3=Axis(-5.05:0.7:5.5,"b") # initial grid in y direction
#
#
# using PyPlot
# pygui(true)
#
# mymdbm=MDBM_Problem(foo,[ax1,ax2,ax3])
#
# for kinter in 1:1
# MDBM.solve!(mymdbm,1)
#
#
# x_eval,y_eval=getevaluatedpoints(mymdbm)
#
# x_sol,y_sol=getinterpolatedsolution(mymdbm)
#
# fig = figure(11);clf()
# PyPlot.scatter(x_eval,y_eval,s=5)
# PyPlot.scatter(x_sol,y_sol)
#
# pause(0.01)
# end
#
#
# x_eval,y_eval=getevaluatedpoints(mymdbm)
# x_sol,y_sol=getinterpolatedsolution(mymdbm)
# 3
# fig = figure(11);clf()
# PyPlot.scatter(x_eval,y_eval,s=5)
# PyPlot.scatter(x_sol,y_sol,s=5)


text:
export MeterProvider, metrics

using Dates: time

const N_MAX_METRICS = 1_000

Base.@kwdef mutable struct Metric{A<:AbstractAggregation}
    name::String
    description::Union{String,Nothing}
    attribute_keys::Union{Tuple{Vararg{String}},Nothing}
    aggregation::A
    instrument::AbstractInstrument
end

Base.iterate(m::Metric, args...) = iterate(m.aggregation, args...)
Base.getindex(m::Metric, k) = getindex(m.aggregation, k)
Base.length(m::Metric) = length(m.aggregation)
OpenTelemetryAPI.resource(m::Metric) = resource(m.instrument)

function (metric::Metric)(m::Measurement)
    if isnothing(metric.attribute_keys)
        filtered_attributes = StaticAttrs()
    else
        interested_keys = keys(m.attributes) ∩ metric.attribute_keys
        filtered_keys = setdiff(keys(m.attributes), interested_keys)
        m = Measurement(m.value, m.attributes[Tuple(interested_keys)])
        filtered_attributes = m.attributes[Tuple(filtered_keys)]
    end

    span_ctx = span_context()
    if isnothing(span_ctx)
        trace_id = INVALID_TRACE_ID
        span_id = INVALID_SPAN_ID
    else
        trace_id = span_ctx.trace_id
        span_id = span_ctx.span_id
    end

    exemplar = Exemplar(;
        value = m,
        time_unix_nano = UInt(time() * 10^9),
        filtered_attributes = filtered_attributes,
        trace_id = trace_id,
        span_id = span_id,
    )
    metric.aggregation(exemplar)
end

struct MeterProvider <: AbstractMeterProvider
    resource::Resource
    meters::IdDict{Meter,Nothing}
    instrument_linked_metrics::IdDict{AbstractInstrument,IdDict{Metric,Nothing}}
    async_instruments::IdDict{AbstractAsyncInstrument,Nothing}
    views::Vector{View}
    named_view_linked_ins::IdDict{View,AbstractInstrument}
    metrics::IdDict{Metric,Nothing}
    n_max_metrics::UInt
end

OpenTelemetryAPI.resource(p::MeterProvider) = p.resource

"""
    metrics([global_meter_provider()])
"""
metrics() = metrics(global_meter_provider())

metrics(::OpenTelemetryAPI.DummyMeterProvider) = ()
metrics(p::MeterProvider) = keys(p.metrics)

"""
    MeterProvider(;resource = Resource(), views = View[], n_max_metrics = N_MAX_METRICS)

If `views` is empty, a default one ([`View(;instrument_name="*")`](@ref)) will be added to enable all metrics.
"""
function MeterProvider(;
    resource = Resource(),
    views = View[],
    n_max_metrics = N_MAX_METRICS,
)
    if isempty(views)
        push!(views, View(; instrument_name = "*"))
    end

    MeterProvider(
        resource,
        IdDict{Meter,Nothing}(),
        IdDict{AbstractInstrument,IdDict{Metric,Nothing}}(),
        IdDict{AbstractAsyncInstrument,Nothing}(),
        views,
        IdDict{View,AbstractInstrument}(),
        IdDict{Metric,Nothing}(),
        n_max_metrics,
    )
end

function Base.push!(p::MeterProvider, m::Meter)
    p.meters[m] = nothing
    for ins in m.instruments
        push!(p, ins)
    end
end

function Base.push!(p::MeterProvider, ins::AbstractInstrument)
    # 1. Register the meter that the `ins` belongs to
    if !haskey(p.meters, ins.meter)
        p.meters[ins.meter] = nothing
    end

    # 2. Register the `ins`
    drop_views = (v for v in p.views if v.aggregation === DROP)
    valid_views = (v for v in p.views if v.aggregation !== DROP)

    is_drop = false
    for v in drop_views
        if occursin(ins, v.criteria)
            is_drop = true
            break
        end
    end

    # 2.1 Check if the `ins`` is configured to be dropped
    if is_drop
        @info "Metric won't be created since the view for the instrument [$(ins.name)] is configured to DROP."
    else
        for v in valid_views
            # 2.2 try to find the matching views
            if occursin(ins, v.criteria)
                # 2.2.1 for each valid ins => view pair, we're going to
                # create a metric for them So better to check if we've
                # reached the maximum number of metrics configured in the
                # provider
                if length(p.metrics) >= p.n_max_metrics
                    @warn "Maximum number of metrics [$(p.n_max_metrics)] reached. Instrument [$(ins.name)] related metrics are dropped!"
                else
                    metric = Metric(
                        name = something(v.name, ins.name),
                        description = something(v.description, ins.description),
                        attribute_keys = v.attribute_keys,
                        aggregation = something(v.aggregation, default_aggregation(ins)),
                        instrument = ins,
                    )

                    # https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/metrics/sdk.md#view
                    # > In order to avoid conflicts, views which specify
                    # > a name SHOULD have an instrument selector that
                    # > selects at most one instrument. For the
                    # > registration mechanism described above, where
                    # > selection is provided via configuration, the SDK
                    # > MUST NOT allow Views with a specified name to be
                    # > declared with instrument selectors that select
                    # > more than one instrument (e.g. wild card
                    # > instrument name).
                    if !isnothing(v.name)
                        if haskey(p.named_view_linked_ins, v)
                            @info "This named view [$(v.name)] is already registered with an instrument"
                            continue  # !!! important
                        else
                            p.named_view_linked_ins[v] = ins
                        end
                    end

                    get!(p.instrument_linked_metrics, ins, IdDict{Metric,Nothing}())[metric] =
                        nothing
                    p.metrics[metric] = nothing

                    if ins isa AbstractAsyncInstrument
                        p.async_instruments[ins] = nothing
                    end
                end
            end
        end
    end
end

function Base.push!(p::MeterProvider, ins_m::Pair{<:AbstractInstrument,<:Measurement})
    instrument, measurement = ins_m
    if haskey(p.instrument_linked_metrics, instrument)
        for metric in keys(p.instrument_linked_metrics[instrument])
            metric(measurement)
        end
    else
        @debug "Instrument [$(instrument.name)] is not registered in the meter provider."
    end
end

function Base.push!(p::MeterProvider, ins_m::Pair{<:AbstractInstrument})
    ins, m = ins_m
    push!(p, ins => Measurement(m))
end


text:
# using Revise

using Test

using RoME

include(joinpath(@__DIR__,"BeehiveTestUtils.jl"))

@testset "sanity check on Hex example" begin

## start with an empty factor graph object
fg = initfg()

posecount = 0

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)
posecount += 1

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0.0; 0.0; 0.0],
                                           Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=false )

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x0; :l1], p2br, autoinit=false )

## hex 1

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

# writeGraphPdf(fg,engine="neato")

getSolverParams(fg).drawtree = false
getSolverParams(fg).showtree = false
getSolverParams(fg).downsolve = false
getSolverParams(fg).multiproc = false
getSolverParams(fg).async = false
# to disable parent factor sharing
# getSolverParams(fg).devParams[:dontUseParentFactorsInitDown] = ""

# direct solve would be
tree, smt, hist = solveTree!(fg) #, recordcliqs=ls(fg))

@test 55 < sum(-3.0 .< getPoints(getKDE(fg, :x0))[1,:] .< 3.0)
@test 55 < sum(-3.0 .< getPoints(getKDE(fg, :x0))[2,:] .< 3.0)
@test 55 < sum(-0.3 .< getPoints(getKDE(fg, :x0))[3,:] .< 0.3)

@test 55 < sum(7.0 .< getPoints(getKDE(fg, :x1))[1,:] .< 13.0)
@test 55 < sum(-3.0 .< getPoints(getKDE(fg, :x1))[2,:] .< 3.0)
@test 55 < sum(0.7 .< getPoints(getKDE(fg, :x1))[3,:] .< 1.3)

@test 55 < sum(12.0 .< getPoints(getKDE(fg, :x2))[1,:] .< 18.0)
@test 55 < sum(6.0 .< getPoints(getKDE(fg, :x2))[2,:] .< 11.0)
@test 55 < sum(1.8 .< getPoints(getKDE(fg, :x2))[3,:] .< 2.4)

@test 55 < sum(7.0 .< getPoints(getKDE(fg, :x3))[1,:] .< 13.0)
@test 55 < sum(15.0 .< getPoints(getKDE(fg, :x3))[2,:] .< 20.0)
# @test 55 < sum(-0.3 .< getPoints(getKDE(fg, :x3))[3,:] .< 0.3)

@test 55 < sum(-5.0 .< getPoints(getKDE(fg, :x4))[1,:] .< 5.0)
@test 55 < sum(13.0 .< getPoints(getKDE(fg, :x4))[2,:] .< 22.0)
@test 55 < sum(-2.8 .< getPoints(getKDE(fg, :x4))[3,:] .< -1.5)

@test 55 < sum(-8.0 .< getPoints(getKDE(fg, :x5))[1,:] .< -2.0)
@test 55 < sum(6.0 .< getPoints(getKDE(fg, :x5))[2,:] .< 11.0)
@test 55 < sum(-1.3 .< getPoints(getKDE(fg, :x5))[3,:] .< -0.7)

@test 55 < sum(-3.0 .< getPoints(getKDE(fg, :x6))[1,:] .< 3.0)
@test 55 < sum(-3.0 .< getPoints(getKDE(fg, :x6))[2,:] .< 3.0)
@test 55 < sum(-0.3 .< getPoints(getKDE(fg, :x6))[3,:] .< 0.3)

@test 55 < sum(17.0 .< getPoints(getKDE(fg, :l1))[1,:] .< 23.0)
@test 55 < sum(-5.0 .< getPoints(getKDE(fg, :l1))[2,:] .< 5.0)

# #  Do some plotting
# using RoMEPlotting
#
# Gadfly.set_default_plot_size(35cm,20cm)
# pl = drawPosesLandms(fg, meanmax=:max, drawhist=true)
# pl |> PDF("/tmp/caesar/test.pdf");  @async run(`evince /tmp/caesar/test.pdf`)
# pl |> PNG("/tmp/caesar/test.png");
# drawTree(tree, imgs=true)

end # testset

## OR manually do clique by clique
#
# # solve by hand, one cliq at a time
# tree = wipeBuildNewTree!(fg)
# # drawTree(tree,show=true)
#
#
# # solve the first cliq
# smt, hist = solveCliq!(fg, tree, :x0, recordcliq=true)
#
# # solve second clq
# smt, hist = solveCliq!(fg, tree, :x2, cliqHistories=hist, recordcliq=true)
# # resetCliqSolve!(fg, tree, :x2)
#
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x4, cliqHistories=hist, recordcliq=true)
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x6, cliqHistories=hist, recordcliq=true)
# # solve forth clq
# smt, hist = solveCliq!(fg, tree, :x3, cliqHistories=hist, recordcliq=true)


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666],
     y=[0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512],
     color=[0.0770269, 0.685074, 0.988958, 0.2822, 0.824017, 0.69052, 0.077905, 0.953607, 0.0857399, 0.581857],
     Scale.color_asinh)


text:
# This file is a part of TextSearch.jl

import KCenters
export update!, joinmodel
"""
    update!(a::VectorModel, b::VectorModel)

Updates `a` with `b` inplace; returns `a`.
"""
function update!(a::VectorModel, b::VectorModel)
m = a.m
for (token, idfreq) in b.tokens
    x = get(a.tokens, token, nothing)
    if x === nothing
        m += 1
        a.tokens[token] = IdFreq(m, x.freq)
    else
        a.tokens[token] = IdFreq(idfreq.id, idfreq.freq + x.freq)
    end
end

a.maxfreq = max(a.maxfreq, b.maxfreq)
a.n += b.n
a.m = m
a.id2token = Dict(idfreq.id => t for (t, idfreq) in a.tokens)
a
end

"""
    joinmodel(arr::AbstractVector{VectorModel})
    joinmodel(arr::AbstractVector{EntModel})

Joins a list of models into a single one; frequencies and weights are computed as the mean of its non-missing occurrences
"""
function joinmodel(arr::AbstractVector{VectorModel})
    item = first(arr)
    # Vocabulary = Dict{Symbol, IdFreq}
    tokens_ = Dict{Symbol,Int}()
    counts_ = Dict{Symbol,Int32}()

    maxfreq = 0
    for model in arr
        maxfreq = max(maxfreq, model.maxfreq)
        for (k, v) in model.tokens
            tokens_[k] = get(tokens_, k, 0.0) + v.freq # ceil(Int, 0.5 * (
            counts_[k] = get(counts_, k, 0) + 1
        end
    end

    tokens = TextSearch.Vocabulary()
    for (k, v) in tokens_
        tokens[k] = TextSearch.IdFreq(length(tokens) + 1, ceil(Int, v / counts_[k]))
    end

    id2token = Dict(idfreq.id => token for (token, idfreq) in tokens)
    VectorModel(item.config, tokens, id2token, maxfreq, length(tokens_), item.n)
end

function glue(arr::AbstractVector{EntModel})
    item = first(arr)
    tokens_ = Dict{Symbol,Float32}()
    counts_ = Dict{Symbol,Int32}()

    for model in arr
        for (k, v) in model.tokens
            tokens_[k] = get(tokens_, k, 0.0) + v.weight
            counts_[k] = get(counts_, k, 0) + 1
        end
    end

    tokens = TextSearch.WeightedVocabulary()
    for (k, v) in tokens_
        tokens[k] = TextSearch.IdWeight(length(tokens) + 1, v / counts_[k])
    end

    id2token = Dict(w.id => token for (token, w) in tokens)

    EntModel(item.config, tokens, id2token, length(tokens), item.n)
end

max_stars_repo_name: sadit/TextModel

max_issues_repo_name: sadit/TextModel

max_forks_repo_name: sadit/TextModel


text:
using HDF5,JLD,MAT
macro date(_x) :(println("$(now()) "*$(string(_x)));flush(STDOUT);@time $(esc(_x))) end
@date require("KUparser")
@date d = load(ARGS[1])
d = collect(values(d))
@assert length(d)==1 "$ARGS[1] has more than one variable"
d = d[1]
feats = eval(parse("KUparser.Flist.$(ARGS[2])"))
f = fill(feats, length(d))
@date p = pmap(KUparser.oparse, d, f)
x = map(z->z[2], p)
y = map(z->z[3], p)
matwrite(ARGS[3], Dict("x" => x, "y" => y))


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct DisksApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a disk.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskName::String (required)
Param: api_version::String (required)
Param: disk::Disk (required)
Return: Disk
"""
function _swaggerinternal_disksCreateOrUpdate(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, disk; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", Disk, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks/{diskName}", ["azure_auth"], disk)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskName", diskName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksCreateOrUpdate(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, disk; _mediaType=nothing)
    _ctx = _swaggerinternal_disksCreateOrUpdate(_api, subscriptionId, resourceGroupName, diskName, api_version, disk; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksCreateOrUpdate(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, disk; _mediaType=nothing)
    _ctx = _swaggerinternal_disksCreateOrUpdate(_api, subscriptionId, resourceGroupName, diskName, api_version, disk; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes a disk.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_disksDelete(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks/{diskName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskName", diskName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksDelete(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksDelete(_api, subscriptionId, resourceGroupName, diskName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksDelete(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksDelete(_api, subscriptionId, resourceGroupName, diskName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets information about a disk.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskName::String (required)
Param: api_version::String (required)
Return: Disk
"""
function _swaggerinternal_disksGet(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Disk, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks/{diskName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskName", diskName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksGet(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksGet(_api, subscriptionId, resourceGroupName, diskName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksGet(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksGet(_api, subscriptionId, resourceGroupName, diskName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Grants access to a disk.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskName::String (required)
Param: api_version::String (required)
Param: grantAccessData::GrantAccessData (required)
Return: AccessUri
"""
function _swaggerinternal_disksGrantAccess(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, grantAccessData; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", AccessUri, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks/{diskName}/beginGetAccess", ["azure_auth"], grantAccessData)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskName", diskName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksGrantAccess(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, grantAccessData; _mediaType=nothing)
    _ctx = _swaggerinternal_disksGrantAccess(_api, subscriptionId, resourceGroupName, diskName, api_version, grantAccessData; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksGrantAccess(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, grantAccessData; _mediaType=nothing)
    _ctx = _swaggerinternal_disksGrantAccess(_api, subscriptionId, resourceGroupName, diskName, api_version, grantAccessData; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all the disks under a subscription.
Param: subscriptionId::String (required)
Param: api_version::String (required)
Return: DiskList
"""
function _swaggerinternal_disksList(_api::DisksApi, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DiskList, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/disks", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksList(_api::DisksApi, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksList(_api, subscriptionId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksList(_api::DisksApi, response_stream::Channel, subscriptionId::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksList(_api, subscriptionId, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all the disks under a resource group.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Return: DiskList
"""
function _swaggerinternal_disksListByResourceGroup(_api::DisksApi, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", DiskList, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksListByResourceGroup(_api::DisksApi, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksListByResourceGroup(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksListByResourceGroup(_api, subscriptionId, resourceGroupName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Revokes access to a disk.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskName::String (required)
Param: api_version::String (required)
Return: Nothing
"""
function _swaggerinternal_disksRevokeAccess(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks/{diskName}/endGetAccess", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskName", diskName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksRevokeAccess(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksRevokeAccess(_api, subscriptionId, resourceGroupName, diskName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksRevokeAccess(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_disksRevokeAccess(_api, subscriptionId, resourceGroupName, diskName, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates (patches) a disk.
Param: subscriptionId::String (required)
Param: resourceGroupName::String (required)
Param: diskName::String (required)
Param: api_version::String (required)
Param: disk::DiskUpdate (required)
Return: Disk
"""
function _swaggerinternal_disksUpdate(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, disk; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", Disk, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Compute/disks/{diskName}", ["azure_auth"], disk)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "diskName", diskName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function disksUpdate(_api::DisksApi, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, disk; _mediaType=nothing)
    _ctx = _swaggerinternal_disksUpdate(_api, subscriptionId, resourceGroupName, diskName, api_version, disk; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function disksUpdate(_api::DisksApi, response_stream::Channel, subscriptionId::String, resourceGroupName::String, diskName::String, api_version::String, disk; _mediaType=nothing)
    _ctx = _swaggerinternal_disksUpdate(_api, subscriptionId, resourceGroupName, diskName, api_version, disk; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export disksCreateOrUpdate, disksDelete, disksGet, disksGrantAccess, disksList, disksListByResourceGroup, disksRevokeAccess, disksUpdate


max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
"""
	getrange(n)

For a thread, returns respective range within 1:`n` to iterate over.
"""
@inline function getrange(n)
	tid = Threads.threadid()
	nt = Threads.nthreads()
	d , r = divrem(n, nt)
	from = (tid - 1) * d + min(r, tid - 1) + 1
	to = from + d - 1 + (tid ≤ r ? 1 : 0)
	from:to
end

"""
	getindX(supportX,val)

Returns index of `val` within `supportX`. Requires `supportX` to be sorted.
"""
function getindX(supportX,val)
	searchsorted(supportX,val)[1]
end

"""
	estimateΠ(X,supportX)

Returns estimate of transition matrix `Π`, using data `X` and support `supportX`.
"""
function estimateΠ(X,supportX)
	nperiod = size(X,1)
	nfirm = size(X,2)

	Π = zeros(Int,length(supportX),length(supportX))
	Π0 = zeros(Int,length(supportX))

	for t in 2:nperiod, f in 1:nfirm
			_ix = getindX(supportX,X[t-1,f])
			Π[_ix,getindX(supportX,X[t,f])] += 1
			Π0[_ix]+=1
	end

	return Π ./ Π0
end

max_stars_repo_name: rgreminger/DDCModelsExample

max_issues_repo_name: rgreminger/DDCModelsExample

max_forks_repo_name: rgreminger/DDCModelsExample


text:
using Base.Test

using DataReaders

include("test_google_daily.jl")
# include("test_google_quote.jl")  # ToDo
include("test_yahoo_daily.jl")

@test 1 == 1


text:
include("misc.jl")
include("findMin.jl")

function logReg(X,y)

	(n,d) = size(X)

	# Initial guess
	w = zeros(d,1)

	# Function we're going to minimize (and that computes gradient)
	funObj(w) = logisticObj(w,X,y)

	# Solve least squares problem
	w = findMin(funObj,w,derivativeCheck=true)

	# Make linear prediction function
	predict(Xhat) = sign.(Xhat*w)

	# Return model
	return LinearModel(predict,w)
end

function logisticObj(w,X,y)
	yXw = y.*(X*w)
	f = sum(log.(1 + exp.(-yXw)))
	g = -X'*(y./(1+exp.(yXw)))
	return (f,g)
end

# Multi-class one-vs-all version (assumes y_i in {1,2,...,k})
function logRegOnevsAll(X,y)
	(n,d) = size(X)
	k = maximum(y)

	# Each column of 'w' will be a logistic regression classifier
	W = zeros(d,k)

	for c in 1:k
		yc = ones(n,1) # Treat class 'c' as +1
		yc[y .!= c] = -1 # Treat other classes as -1

		# Each binary objective has the same features but different lables
		funObj(w) = logisticObj(w,X,yc)

		W[:,c] = findMin(funObj,W[:,c],verbose=false)
	end

	# Make linear prediction function
	predict(Xhat) = mapslices(indmax,Xhat*W,2)

	return LinearModel(predict,W)
end

# Multi-class softmax version (assumes y_i in {1,2,...,k})
function logRegSoftmax(X,y)
	(n,d) = size(X)
	k = maximum(y)

	# Each column of 'w' will be a logistic regression classifier
	W = zeros(d,k)

	funObj(w) = softmaxObj(w,X,y,k)

	W[:] = findMin(funObj,W[:],derivativeCheck=true,maxIter=500)

	# Make linear prediction function
	predict(Xhat) = mapslices(indmax,Xhat*W,2)

	return LinearModel(predict,W)
end

function softmaxObj(w,X,y,k)
	(n,d) = size(X)

	W = reshape(w,d,k)

	XW = X*W
	Z = sum(exp.(XW),2)

	nll = 0
	G = zeros(d,k)
	for i in 1:n
		nll += -XW[i,y[i]] + log(Z[i])

		pVals = exp.(XW[i,:])./Z[i]
		for c in 1:k
			G[:,c] += X[i,:]*(pVals[c] - (y[i] == c))
		end
	end
	return (nll,reshape(G,d*k,1))
end


text:
using ControlSystems
using Plots

m=1;
b=1;
k=2;

sys = tf([1,0,1], [m,b,m,b,k]);

impulse(sys)

impulseplot(sys)

max_stars_repo_name: akshayhiregoudar/Control-System-Design

max_issues_repo_name: akshayhiregoudar/Control-System-Design

max_forks_repo_name: akshayhiregoudar/Control-System-Design


text:
abstract type AbstractReachSet{SN} end

"""
    ReachSet{SN}

Type that wraps a reach set, which is a set representing (an approximation of)
the reachable states for a given time interval.

### Fields

- `X`       -- set
- `t_start` -- time interval lower bound
- `t_end`   -- time interval upper bound
"""
struct ReachSet{SN} <: AbstractReachSet{SN}
    X::SN
    t_start::Float64
    t_end::Float64
end

function set(rs::ReachSet)
    return rs.X
end

function time_start(rs::ReachSet)
    return rs.t_start
end

function time_end(rs::ReachSet)
    return rs.t_end
end

function substitute(rs::ReachSet; set=set(rs), time_start=time_start(rs),
                    time_end=time_end(rs))
    return ReachSet(set, time_start, time_end)
end

function project(rs::ReachSet, M::AbstractMatrix)
    @assert dim(rs.X) == size(M, 2) "a projection of size $(size(M, 2)) is " *
        "incompatible with a set of dimension $(dim(rs.X))"
    return ReachSet(M * rs.X, rs.t_start, rs.t_end)
end

# =========================
# low-dimensional reach set
# =========================

struct SparseReachSet{SN} <: AbstractReachSet{SN}
    rs::ReachSet{SN}
    dimensions::AbstractVector{Int}
end

# flattened constructor
function SparseReachSet(X, t_start::Float64, t_end::Float64,
                        dimensions::AbstractVector{Int})
    return SparseReachSet(ReachSet(X, t_start, t_end), dimensions)
end

function set(srs::SparseReachSet)
    return set(srs.rs)
end

function time_start(srs::SparseReachSet)
    return time_start(srs.rs)
end

function time_end(srs::SparseReachSet)
    return time_end(srs.rs)
end

function substitute(srs::SparseReachSet; set=set(srs),
                    time_start=time_start(srs), time_end=time_end(srs))
    substituted_rs = substitute(srs.rs, set=set, time_start=time_start, time_end=time_end)
    return SparseReachSet(substituted_rs, srs.dimensions)
end

function project(srs::SparseReachSet, M::AbstractMatrix)
    rs = srs.rs
    M_proj = M[:, srs.dimensions]
    return project(rs, M_proj)
end


text:
using RegressionAndOtherStories
using Test

include(joinpath(@__DIR__, "test_hibbs.jl"))


text:
module ReinforcementLearningDatasets

const RLDatasets = ReinforcementLearningDatasets
export RLDatasets

using DataDeps

include("d4rl/d4rl/register.jl")
include("d4rl/d4rl_pybullet/register.jl")
include("atari/register.jl")
include("rl_unplugged/atari/register.jl")
include("common.jl")
include("init.jl")
include("d4rl/d4rl_dataset.jl")
include("atari/atari_dataset.jl")
include("rl_unplugged/util.jl")
include("rl_unplugged/atari/rl_unplugged_atari.jl")

end


text:
import Base

mutable struct Reshape <: DiffableFunction
    grad_field::GradField
    in_shape :: Tuple
    out_shape :: Tuple
    Reshape(in_shape, out_shape) = new(GradField(), in_shape, out_shape)
end

function forward(f::Reshape, x)
    f.in_shape = size(x)
    return reshape(x, f.out_shape)
end

function backward(f::Reshape, gy)
    return reshape(gy, f.in_shape)
end

Base.reshape(x::Variable, shape) = Reshape(size(x), shape)(x)
is_support(::typeof(Base.reshape)) = true


text:
# This script provides implementations of
# multi-step methods of different orders
# for solving a system of differential equations
# and compares them to the exact solutions and
# other methods while solving a 1-body problem

println('\n', " "^4, "> Loading the packages...")

using LinearAlgebra
using NLsolve
using Printf

# Define the value of ϰ
const ϰ = -1

# Integrate equations of motion using the
# Euler's method, return the values of position
# and velocity on the last step
function euler(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Prepare the output vectors
    r = copy(r₀)
    v = copy(v₀)
    # Compute the solution
    for _ in 1:n
        a = ϰ * r / norm(r)^3
        r += h * v
        v += h * a
    end
    return r, v
end

# Integrate equations of motion using the
# two-step Adams–Bashforth's method, return
# the values of position and velocity on the
# last step
function _ab2(
    r₀::Vector{F},
    v₀::Vector{F},
    r::Vector{F},
    v::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Define a couple of independent coefficients
    k₁ = 3 / 2 * h
    k₂ = -1 / 2 * h
    # Compute the rest in two steps
    for _ in 2:n
        a₁ = k₁ * ϰ * r / norm(r)^3
        a₂ = k₂ * ϰ * r₀ / norm(r₀)^3
        r₀ = r
        r += k₁ * v + k₂ * v₀
        v₀ = v
        v += a₁ + a₂
    end
    return r, v
end

# Integrate equations of motion using the
# two-step Adams–Bashforth's method (with
# the Euler method as a starter), return
# the values of position and velocity on
# the last step
function ab2(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the one-step Euler's method
    r, v = euler(r₀, v₀, h, UInt(1))
    # Integrate and return results
    return _ab2(r₀, v₀, r, v, h, n)
end

# Integrate equations of motion using the
# two-step Adams–Bashforth's method (with
# the 4th-order Runge-Kutta method as a starter),
# return the values of position and velocity on
# the last step
function ab2_rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the 4th-order Runge-Kutta method
    r, v = rk4(r₀, v₀, h, UInt(1))
    # Integrate and return results
    return _ab2(r₀, v₀, r, v, h, n)
end

# Integrate equations of motion using the
# three-step Adams–Bashforth's method, return
# the values of position and velocity on the
# last step
function _ab3(
    r₀::Vector{F},
    v₀::Vector{F},
    r₁::Vector{F},
    v₁::Vector{F},
    r::Vector{F},
    v::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Define independent coefficients
    k₁ = 23 / 12 * h
    k₂ = -4 / 3 * h
    k₃ = 5 / 12 * h
    # Compute the rest in three steps
    for _ in 3:n
        a₁ = k₁ * ϰ * r / norm(r)^3
        a₂ = k₂ * ϰ * r₁ / norm(r₁)^3
        a₃ = k₃ * ϰ * r₀ / norm(r₀)^3
        r₀, r₁ = r₁, r
        r += k₁ * v + k₂ * v₁ + k₃ * v₀
        v₀, v₁ = v₁, v
        v += a₁ + a₂ + a₃
    end
    return r, v
end

# Integrate equations of motion using the
# three-step Adams–Bashforth's method (with
# the Euler method and the two-step
# Adams–Bashforth's method as starters),
# return the values of position and velocity
# on the last step
function ab3(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the one-step Euler's method
    r₁, v₁ = euler(r₀, v₀, h, UInt(1))
    # Compute the third value of the solution
    # using the two-step Adams–Bashforth's method
    r, v = ab2(r₀, v₀, h, UInt(2))
    return _ab3(r₀, v₀, r₁, v₁, r, v, h, n)
end

# Integrate equations of motion using the
# three-step Adams–Bashforth's method (with
# the 4th-order Runge-Kutta method as a starter),
# return the values of position and velocity
# on the last step
function ab3_rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second and third values of the
    # solution using the 4th-order Runge-Kutta method
    r₁, v₁ = rk4(r₀, v₀, h, UInt(1))
    r, v = rk4(r₀, v₀, h, UInt(2))
    return _ab3(r₀, v₀, r₁, v₁, r, v, h, n)
end

# Integrate equations of motion using the
# four-step Adams–Bashforth's method, return
# the values of position and velocity on the
# last step
function _ab4(
    r₀::Vector{F},
    v₀::Vector{F},
    r₁::Vector{F},
    v₁::Vector{F},
    r₂::Vector{F},
    v₂::Vector{F},
    r::Vector{F},
    v::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Define independent coefficients
    k₁ = 55 / 24 * h
    k₂ = -59 / 24 * h
    k₃ = 37 / 24 * h
    k₄ = -9 / 24 * h
    # Compute the rest in three steps
    for _ in 3:n
        a₁ = k₁ * ϰ * r / norm(r)^3
        a₂ = k₂ * ϰ * r₂ / norm(r₂)^3
        a₃ = k₃ * ϰ * r₁ / norm(r₁)^3
        a₄ = k₄ * ϰ * r₀ / norm(r₀)^3
        r₀, r₁, r₂ = r₁, r₂, r
        r += k₁ * v + k₂ * v₂ + k₃ * v₁ + k₄ * v₀
        v₀, v₁, v₂ = v₁, v₂, v
        v += a₁ + a₂ + a₃ + a₄
    end
    return r, v
end

# Integrate equations of motion using the
# four-step Adams–Bashforth's method (with
# the Euler method, two-step and three-step
# Adams–Bashforth's methods as starters),
# return the values of position and velocity
# on the last step
function ab4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the one-step Euler's method
    r₁, v₁ = euler(r₀, v₀, h, UInt(1))
    # Compute the third value of the solution
    # using the two-step Adams–Bashforth's method
    r₂, v₂ = ab2(r₀, v₀, h, UInt(2))
    # Compute the third value of the solution
    # using the three-step Adams–Bashforth's method
    r, v = ab3(r₀, v₀, h, UInt(3))
    return _ab4(r₀, v₀, r₁, v₁, r₂, v₂, r, v, h, n)
end

# Integrate equations of motion using the
# four-step Adams–Bashforth's method (with
# the 4th-order Runge-Kutta method as a starter),
# return the values of position and velocity
# on the last step
function ab4_rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second, third, and fourth values of
    # the solution using the 4th-order Runge-Kutta method
    r₁, v₁ = rk4(r₀, v₀, h, UInt(1))
    r₂, v₂ = rk4(r₀, v₀, h, UInt(2))
    r, v = rk4(r₀, v₀, h, UInt(3))
    return _ab4(r₀, v₀, r₁, v₁, r₂, v₂, r, v, h, n)
end

# Integrate equations of motion using the
# 4th-order Runge-Kutta's method, return
# the values of position and velocity on
# the last step
function rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Prepare the output vectors
    r = copy(r₀)
    v = copy(v₀)
    # Compute the solution
    for _ in 1:n
        k₁ = v
        a₁ = ϰ * r / norm(r)^3
        k₂ = v + h * a₁ / 2
        a₂ = ϰ * (r + h * k₁ / 2) / norm(r + h * k₁ / 2)^3
        k₃ = v + h * a₂ / 2
        a₃ = ϰ * (r + h * k₂ / 2) / norm(r + h * k₂ / 2)^3
        k₄ = v + h * a₃
        a₄ = ϰ * (r + h * k₃) / norm(r + h * k₃)^3
        r += h / 6 * (k₁ + 2 * k₂ + 2 * k₃ + k₄)
        v += h / 6 * (a₁ + 2 * a₂ + 2 * a₃ + a₄)
    end
    return r, v
end

# Integrate equations of motion using the
# backward Euler's method, return the values
# of position and velocity on the last step
function b_euler(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Prepare the output vectors
    r = copy(r₀)
    v = copy(v₀)
    # Compute the solution
    for _ in 1:n
        function f!(F, x)
            F[1:N] = -x[1:N] + r + h * x[N+1:end]
            F[N+1:end] = -x[N+1:end] + v + h * ϰ * x[1:N] / norm(x[1:N])^3
        end
        res = nlsolve(f!, [r; v], autodiff=:forward, method=:newton, ftol=1e-16).zero
        r, v = res[1:N], res[N+1:end]
    end
    return r, v
end

# Integrate equations of motion using the
# trapezoidal rule, return the values
# of position and velocity on the last step
function trapezoidal(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Prepare the output vectors
    r = copy(r₀)
    v = copy(v₀)
    # Compute the solution
    for _ in 1:n
        ρₖ = norm(r)^3
        function f!(F, x)
            F[1:N] = -x[1:N] + r + h / 2 * (x[N+1:end] + v)
            F[N+1:end] = -x[N+1:end] + v + h / 2 * ϰ * (x[1:N] / norm([x[1:N]])^3 + r / ρₖ)
        end
        res = nlsolve(f!, [r; v], autodiff=:forward, method=:newton, ftol=1e-16).zero
        r, v = res[1:N], res[N+1:end]
    end
    return r, v
end

# Integrate equations of motion using the
# two-step Adams–Moulton's method, return
# the values of position and velocity on the
# last step
function _am2(
    r₀::Vector{F},
    v₀::Vector{F},
    r::Vector{F},
    v::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Define independent coefficients
    k₁ = 5 / 12 * h
    k₂ = 2 / 3 * h
    k₃ = -1 / 12 * h
    # Compute the rest in two steps
    for _ in 2:n
        k₄ = k₂ * ϰ / norm(r)^3
        k₅ = k₃ * ϰ / norm(r₀)^3
        function g!(F, x)
            F[1:N] = -x[1:N] + r + k₁ * x[N+1:end] + k₂ * v + k₃ * v₀
            F[N+1:end] = -x[N+1:end] + v +
                k₁ * ϰ * x[1:N] / norm(x[1:N])^3 +
                k₄ * r + k₅ * r₀
        end
        res = nlsolve(g!, [r; v]; autodiff=:forward, method=:newton, ftol=1e-16).zero
        r₀, v₀ = r, v
        r, v = res[1:N], res[N+1:end]
    end
    return r, v
end

# Integrate equations of motion using the
# two-step Adams–Moulton's method (with
# the trapezoidal rule as a starter), return
# the values of position and velocity on the
# last step
function am2(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the trapezoidal rule
    r, v = trapezoidal(r₀, v₀, h, UInt(1))
    return _am2(r₀, v₀, r, v, h, n)
end

# Integrate equations of motion using the
# two-step Adams–Moulton's method (with the
# 4th-order Runge-Kutta method as a starter),
# return the values of position and velocity
# on the last step
function am2_rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the 4th-order Runge-Kutta method
    r, v = rk4(r₀, v₀, h, UInt(1))
    return _am2(r₀, v₀, r, v, h, n)
end

# Integrate equations of motion using the
# three-step Adams–Moulton's method, return
# the values of position and velocity on the
# last step
function _am3(
    r₀::Vector{F},
    v₀::Vector{F},
    r₁::Vector{F},
    v₁::Vector{F},
    r::Vector{F},
    v::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Determine the length of the input vectors
    N = length(r₀)
    # Define independent coefficients
    k₁ = 9 / 24 * h
    k₂ = 19 / 24 * h
    k₃ = -5 / 24 * h
    k₄ = 1 / 24 * h
    # Compute the rest in two steps
    for _ in 3:n
        # Define dependent coefficients
        k₅ = k₂ * ϰ / norm(r)^3
        k₆ = k₃ * ϰ / norm(r₁)^3
        k₇ = k₄ * ϰ / norm(r₀)^3
        function g!(F, x)
            F[1:N] = -x[1:N] + r + k₁ * x[N+1:end] + k₂ * v + k₃ * v₁ + k₄ * v₀
            F[N+1:end] = -x[N+1:end] + v +
                   k₁ * ϰ * x[1:N] / norm(x[1:N])^3 +
                   k₅ * r + k₆ * r₁ + k₇ * r₀
        end
        res = nlsolve(g!, [r; v]; autodiff=:forward, method=:newton, ftol=1e-16).zero
        r₀, r₁, v₀, v₁ = r₁, r, v₁, v
        r, v = res[1:N], res[N+1:end]
    end
    return r, v
end

# Integrate equations of motion using the
# three-step Adams–Moulton's method (with
# the trapezoidal rule and the two-step
# Adams–Moulton's method as starters), return
# the values of position and velocity on the
# last step
function am3(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second value of the solution
    # using the trapezoidal rule
    r₁, v₁ = trapezoidal(r₀, v₀, h, UInt(1))
    # Compute the third value of the solution
    # using the two-step Adams–Moulton's method
    r, v = am2(r₀, v₀, h, UInt(2))
    return _am3(r₀, v₀, r₁, v₁, r, v, h, n)
end

# Integrate equations of motion using the
# three-step Adams–Moulton's method (with the
# 4th-order Runge-Kutta method as a starter),
# return the values of position and velocity
# on the last step
function am3_rk4(
    r₀::Vector{F},
    v₀::Vector{F},
    h::F,
    n::I,
)::Tuple{Vector{F},Vector{F}} where
{F <: AbstractFloat,I <: Unsigned}
    # Compute the second and third values of the
    # solution using the 4th-order Runge-Kutta method
    r₁, v₁ = rk4(r₀, v₀, h, UInt(1))
    r, v = rk4(r₀, v₀, h, UInt(2))
    return _am3(r₀, v₀, r₁, v₁, r, v, h, n)
end

println(" "^4, "> Integrating the 1-body problem...")

# Define the initial values of the position and velocity
r₀ = [1.0, 0.0]
v₀ = [0.0, 0.5]

# Calculate integrals and some of the orbit parameters
E₀ = 1 / 2 * norm(v₀)^2 - ϰ^2 / norm(r₀)

# Print the initial values
println(
    '\n',
    " "^6, "r₀: ", r₀[1], " ", r₀[2], '\n',
    " "^6, "v₀: ", v₀[1], " ", v₀[2], '\n',
    " "^6, "E₀: ", E₀,
)

# Define the output directories
const data_dir = "$(@__DIR__)/../data/1-body"
const files = joinpath.(
    data_dir,
    [
        "euler.dat",
        "ab2.dat",
        "ab2_rk4.dat",
        "ab3.dat",
        "ab3_rk4.dat",
        "ab4.dat",
        "ab4_rk4.dat",
        "rk4.dat",
        "b_euler.dat",
        "trapezoidal.dat",
        "am2.dat",
        "am2_rk4.dat",
        "am3.dat",
        "am3_rk4.dat",
    ]
)

# Prepare a list of methods
const methods = (
    euler,
    ab2,
    ab2_rk4,
    ab3,
    ab3_rk4,
    ab4,
    ab4_rk4,
    rk4,
    b_euler,
    trapezoidal,
    am2,
    am2_rk4,
    am3,
    am3_rk4,
)

# Integrate equations of motion using all available
# methods and print the values of position and velocity
# on the last step
function integrate(h::Vector{F}, n::Vector{I}) where {F <: AbstractFloat,I <: Unsigned}
    # Open the data files
    io = open.(files, "a")

    # Prepare results vectors
    r = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    v = [ Vector{Vector{F}}() for _ in 1:length(io) ]
    E = [ Vector{F}() for _ in 1:length(io) ]
    ΔE = [ Vector{F}() for _ in 1:length(io) ]

    # For each pair of parameters
    for i in eachindex(h)
        # Integrate and get the results of the last steps
        for j in eachindex(io)
            rᵢ, vᵢ = methods[j](r₀, v₀, h[i], n[i])

            push!(r[j], rᵢ)
            push!(v[j], vᵢ)

            Eᵢ = 1 / 2 * norm(vᵢ)^2 - ϰ^2 / norm(rᵢ)
            ΔEᵢ = abs(E₀ - Eᵢ)

            push!(E[j], Eᵢ)
            push!(ΔE[j], ΔEᵢ)
        end
    end

    # Get the format string based on the index of the iteration
    function get_format(i)
        Printf.Format(
            "\$ 10^{-$(length("$(n[i])") - 1)} \$ " *
            "& \$ $(n[i] % 10 == 0 ? "10^$(length("$(n[i])") - 1)" : n[i]) \$ " *
            "& \$ %.14f \$ " *
            "& \$ %.14f \$ " *
            "\\\\\n"
        )
    end

    # Print the data in the files in the LaTeX format

    for j in eachindex(io)
        println(io[j], "# Positions")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), r[j][i][1], r[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Velocities")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i),  v[j][i][1], v[j][i][2])
        end
    end

    for j in eachindex(io)
        println(io[j], "\n# Total energy")
        for i in eachindex(h)
            Printf.format(io[j], get_format(i), E[j][i], ΔE[j][i])
        end
        println(io[j])
        close(io[j])
    end
end

# Truncate the previous results
open.(files; truncate=true)

# Integrate n = 2^m iterations
integrate(
    [1e-2, 1e-3, 1e-4, 1e-5], #, 1e-6, 1e-7],
    UInt.([ 10^i for i in 2:5 ]),
)

# Integrate the full cycle
integrate(
    [1e-2, 1e-3, 1e-4, 1e-5], #, 1e-6, 1e-7],
    UInt.([ 271, 2714, 27141, 271408]), #, 2714081, 27140809 ]),
)

println()


text:
"""
    writeOutput(model::MoominModel, fileName; nSol=1, format="standard", type="output", asString=true)

  Function for exporting outputs of the algorithm into a text file in various formats.

# INPUTS

- `model`:     A MoominModel structure that is an output of `runMoomin`
- `fileName`:  String to determine were to write the results

# OPTIONAL INPUTS

- `nSol`:      Int, number of the solution to be printed (default: 1)
- `format`:    String, one of the three options (default: "standard"):

               "standard": output is a delimited file with two columns
                listing the ID and the selected output for each
                reaction (see later)
               "full": a table that lists various outputs as well as other information
                for each reaction
               "json": a json file with the selected output for each reaction

- `type`:       String, one of the four options (default: "output"):

                "output": output colour from one solution
                "input": a priori colours ie colours based solely on data
                "combined": attempted consensus between all enumerated solutions
                 If colour differs between solutions, 6/yellow (indeterminate)
                 is entered
                "frequency": how often a reaction is coloured in all enumerated
                 solutions

- `asString`:    Boolean, if true, output is given in actual colours (eg "red")
                (default: true)

"""
function writeOutput(model::MoominModel, fileName; nSol=1, format="standard", type="output", asString=true)
  if isempty(model.reactions.outputColours)
    @warn "The model contains no solutions."
    return
  elseif nSol > size(model.reactions.outputColours, 2)
    @warn "`nSol` exceeds number of solutions in the model."
    return
  end

  if format=="full"
    toWrite = writeFullOutput(model, nSol, asString)
  else
    if type == "output"
      solution = model.reactions.outputColours[:, nSol]
    elseif type == "input"
      solution = model.reactions.inputColours
    elseif type == "combined"
      solution = model.reactions.combinedOutput
    elseif type == "frequency"
      solution = model.reactions.outputFrequency
    else
      @warn "Cannot recognise type given."
      return
    end
    if asString & (type != "frequency")
      solution = numberToColour.(solution)
    end
    solution = string.(solution)
    if format=="standard"
      toWrite = ["ID"; model.reactions.IDs]
      toWrite = [toWrite [type; solution]]
    elseif format=="json"
      toWrite = "{\n"
      for rInd in 1:nReactions(model)-1
        toWrite = toWrite * "\"" * model.reactions.IDs[rInd] * "\": " *
                    solution[rInd] * ",\n"
      end
      toWrite = toWrite * "\"" * model.reactions.IDs[end] * "\": " *
                  solution[end] * "\n}"
      open(fileName, "w") do f
        write(f, toWrite)
      end
      return
    else
      @warn "Cannot recognise format given."
    end
  end
  writedlm(fileName, toWrite)
end

function numberToColour(n)
  if n == 2
    return "r.red"
  elseif n==1
    return "red"
  elseif n==0
    return "grey"
  elseif n==-1
    return "blue"
  elseif n==-2
    return "r.blue"
  elseif n==6
    return "yellow"
  else
    return "potato"
  end
end

function writeFullOutput(model::MoominModel, nSol, asString)
  toWrite = ["ID"; model.reactions.IDs]
  if asString
    toWrite = [toWrite ["input"; numberToColour.(model.reactions.inputColours)]]
    toWrite = [toWrite ["sol#$nSol"; numberToColour.(model.reactions.outputColours[:, nSol])]]
    toWrite = [toWrite ["combined"; numberToColour.(model.reactions.combinedOutput)]]
  else
    toWrite = [toWrite ["input"; model.reactions.inputColours]]
    toWrite = [toWrite ["sol#$nSol"; model.reactions.outputColours[:, nSol]]]
    toWrite = [toWrite ["combined"; model.reactions.combinedOutput]]
  end
  toWrite = [toWrite ["frequency"; model.reactions.outputFrequency]]
  toWrite = [toWrite ["subsystem"; model.reactions.subsystems]]
  genes = ["NA"; model.genes.IDs]
  PPDEs = ["NA"; model.genes.PPDE]
  FCs = ["NA"; model.genes.logFC]
  adjustedInds = model.reactions.leadingGenes.+1
  toWrite = [toWrite ["leadingGene"; genes[adjustedInds]]]
  toWrite = [toWrite ["PPDE"; PPDEs[adjustedInds]]]
  toWrite = [toWrite ["logFC"; FCs[adjustedInds]]]
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
taskid(t=current_task()) = hex(hash(t) & 0xffff, 4)

macro debug(n::Int, s)
    DEBUG_LEVEL >= n ? :(println("DEBUG: ", taskid(), " ", $(esc(s)))) :
                       :()
end

macro debugshow(n::Int, s)
    DEBUG_LEVEL >= n ? :(println("DEBUG: ", taskid(), " ",
                                 $(sprint(Base.show_unquoted, s)), " = ",
                                 sprint(io->show(io, "text/plain",
                                                 begin value=$(esc(s)) end)))) :
                       :()

end

macro debugshort(n::Int, s)
    DEBUG_LEVEL >= n ? :(println("DEBUG: ", taskid(), " ",
                                 sprint(showcompact, $(esc(s))))) :
                       :()
end

printlncompact(x) = println(sprint(showcompact, x))

@noinline function precondition_error(msg, frame)
    msg = string(#=sprint(StackTraces.show_spec_linfo,
                        StackTrace0.lookup(frame)[2])=# "function",
                 " requires ", msg)
    return ArgumentError(msg)
end

"""
    @require precondition [message]
Throw `ArgumentError` if `precondition` is false.
"""
macro require(condition, msg = string(condition))
    esc(:(if ! $condition throw(precondition_error($msg, backtrace()[1])) end))
end

@noinline function postcondition_error(msg, frame, ls="", l="", rs="", r="")
    msg = string(#= sprint(StackTraces.show_spec_linfo,
                        StackTraces.lookup(frame)[2]) =# "function",
                 " failed to ensure ", msg)
    if ls != ""
        msg = string(msg, "\n", ls, " = ", sprint(show, l),
                          "\n", rs, " = ", sprint(show, r))
    end
    return AssertionError(msg)
end

# Copied from stdlib/Test/src/Test.jl:get_test_result()
iscondition(ex) = isa(ex, Expr) &&
                  ex.head == :call &&
                  length(ex.args) == 3 &&
                  first(string(ex.args[1])) != '.' &&
                  (!isa(ex.args[2], Expr) || ex.args[2].head != :...) &&
                  (!isa(ex.args[3], Expr) || ex.args[3].head != :...) &&
                  (ex.args[1] === :(==) ||
                       Base.operator_precedence(ex.args[1]) ==
                           Base.operator_precedence(:(==)))

"""
    @ensure postcondition [message]
Throw `ArgumentError` if `postcondition` is false.
"""
macro ensure(condition, msg = string(condition))

    if DEBUG_LEVEL < 0
        return :()
    end

    if iscondition(condition)
        l,r = condition.args[2], condition.args[3]
        ls, rs = string(l), string(r)
        return esc(quote
            if ! $condition
                # FIXME double-execution of condition l and r!
                throw(postcondition_error($msg, backtrace()[1],
                                          $ls, $l, $rs, $r))
            end
        end)
    end

    esc(:(if ! $condition throw(postcondition_error($msg, backtrace()[1])) end))
end

# FIXME
# Should this have a branch-prediction hint? (same for @assert?)
# http://llvm.org/docs/BranchWeightMetadata.html#built-in-expect-instructions

#=
macro src()
    @static if VERSION >= v"0.7-" && length(:(@test).args) == 2
        esc(quote
            (__module__,
             __source__.file == nothing ? "?" : String(__source__.file),
             __source__.line)
        end)
    else
        esc(quote
            (current_module(),
             (p = Base.source_path(); p == nothing ? "REPL" : p),
             Int(unsafe_load(cglobal(:jl_lineno, Cint))))
        end)
    end
end
=#


text:
module FakePackageTests
# this is used to test `load` keyword of retest, but also to manually check that Revise
# works as expected: e.g. do `using Revise, FakePackageTests` and `retest(FakePackageTests)`
# multiple times, while changing the testcode in this file or the included one "sub.jl"

using ReTest
__revise_mode__ = :eval

@testset "toplevel" begin
    @test true
end

module Sub
using ReTest
__revise_mode__ = :eval

@testset "Sub" begin
    @test true
end

include("./"*"sub.jl")

end # Sub
end # FakePackageTests

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
#push!(LOAD_PATH, "../MocosSim")

module MocosSimLauncher

using ArgParse
using Base.Threads
using CodecZlib
using DataFrames
using FileIO
using JLD2
using JSON
using ProgressMeter
using Random
using Setfield
using TOML

import MocosSim
import MocosSim: ContactKind, NoContact, contactkind, time

const OptTimePoint = Union{Missing, MocosSim.TimePoint}
optreal2float32(optreal::Union{Missing,T} where T<:Real) = ismissing(optreal) ? NaN32 : Float32(optreal)

include("cmd_parsing.jl")
include("callback.jl")
include("load_params.jl")
include("outputs.jl")

export launch

function string2enum(enum_group::Type{<:Enum{T}}, str::AbstractString) where {T<:Integer}
    sym = Symbol(str)
    for val in instances(enum_group)
        if sym == Symbol(val)
            return val
        end
    end
    error("not found $str in $enum_group")
end

dict2kwargs(dict::Dict{S, Any} where S<:AbstractString) = NamedTuple{Tuple(Symbol.(keys(dict)))}(values(dict))

function launch(args::AbstractVector{T} where T<:AbstractString)
  @info "Stated" nthreads()
  if nthreads() == 1
    @warn "using single thread, set more threads by passing --threads agrument to julia or setting JULIA_NUM_THREADS environment variable"
  end

  cmd_args = parse_commandline(args)
  @info "Parsed args" cmd_args
  config_path = cmd_args["CONFIG"]
  config = endswith(config_path, ".json") ? JSON.parsefile(config_path) : TOML.parsefile(config_path)

  max_num_infected = config["stop_simulation_threshold"] |> Int
  time_limit = get(config, "stop_simulation_time", typemax(MocosSim.TimePoint)) |> MocosSim.TimePoint
  num_trajectories = config["num_trajectories"] |> Int
  params_seed = get(config, "params_seed", 0)

  @info "loading population and setting up parameters" params_seed
  rng = MersenneTwister(params_seed)
  GC.gc()
  params = read_params(config, rng)
  GC.gc()
  num_individuals =  MocosSim.numindividuals(params)

  immunization = nothing
  immune = nothing
  if haskey(config, "initial_conditions")
    initial_conditions = config["initial_conditions"]
    if haskey(initial_conditions, "immunization")
      immunization_cfg = initial_conditions["immunization"]

      if haskey(immunization_cfg, "order_file")
        immunization = load(immunization_cfg["order_file"], "immunization")::MocosSim.Immunization
        enqueue_immunizations = get(immunization_cfg, "enqueue", true) |> Bool
      end

      # keeping legacy immunization for a while
      if haskey(immunization_cfg, "level")
        @warn "legacy immunization in use"
        immunization_ordering::AbstractVector{T} where T<: Integer = load(immunization_cfg["order_data"], "ordering")
        immunization_level::Real = immunization_cfg["level"]

        num_immune = round(UInt, immunization_level * num_individuals)
        immune_ids = @view immunization_ordering[begin : num_immune]
        immune = falses(num_individuals)
        immune[immune_ids] .= true
      end
    end
  end
  @info "allocating simulation states"
  states = [MocosSim.SimState(num_individuals) for _ in 1:nthreads()]
  callbacks = [DetectionCallback(num_individuals, max_num_infected, time_limit) for _ in 1:nthreads()]
  outputs = make_outputs(cmd_args, num_trajectories)

  for o in outputs
    beforetrajectories(o, params)
  end

  @info "starting simulation" num_trajectories
  writelock = ReentrantLock()
  progress = ProgressMeter.Progress(num_trajectories)
  GC.gc()

  outside_case_imports = MocosSim.AbstractOutsideCases[]
  if haskey(config, "imported_cases")
    for outside_import in config["imported_cases"]
      name = outside_import["function"]
      params_dict = get(outside_import, "params", Dict{String,Any}())
      if haskey(params_dict, "strain")
        strain_str = params_dict["strain"] * "Strain"
        params_dict["strain"] = string2enum(MocosSim.StrainKind, strain_str)
      end
      outside_cases = MocosSim.make_imported_cases(name; dict2kwargs(params_dict)...)
      push!(outside_case_imports, outside_cases)
    end
  else
    error("the import function was not used!")
  end

  @threads for trajectory_id in 1:num_trajectories
    state = states[threadid()]
    MocosSim.reset!(state, trajectory_id)

    for outside_fun in outside_case_imports
      outside_fun(state, params)
    end
    if params.screening_params !== nothing
      MocosSim.add_screening!(state, params)
    end

    if immune !== nothing
      immune::AbstractVector{Bool}
      for i in 1:num_individuals
        if !immune[i]
          continue
        end
        individual = state.individuals[i]
        state.individuals[i] = @set individual.health = MocosSim.Recovered
      end
    end

    if immunization !== nothing
      immunization::MocosSim::Immunization
      MocosSim.immunize!(state, immunization, enqueue_immunizations)
    end

    callback = callbacks[threadid()]
    reset!(callback)
    try
      MocosSim.simulate!(state, params, callback)
      for o in outputs
        pushtrajectory!(o, trajectory_id, writelock, state, params, callback)
      end
    catch err
      @warn "Failed on thread " threadid() trajectory_id err
      foreach(x -> println(stderr, x), stacktrace(catch_backtrace()))
    end

    ProgressMeter.next!(progress) # is thread-safe
  end

  for o in outputs
    aftertrajectories(o, params)
  end
end

function julia_main()::Cint
  try
    launch(ARGS)
  catch err
    Base.invokelatest(Base.display_error, Base.catch_stack())
    return 1
  end
  return 0
end

precompile(MocosSim.simulate!, (MocosSim.SimState, MocosSim.SimParams, DetectionCallback))
precompile(launch, (Vector{String},))
precompile(julia_main, ())

end


text:
_check_jld2_nobs(nobs) = all(==(first(nobs)), nobs[2:end])

"""
    JLD2Dataset(file::AbstractString, paths)
    JLD2Dataset(fid::JLD2.JLDFile, paths::Union{String, Vector{String}})

Wrap several JLD2 datasets (`paths`) as a single dataset container.
Each dataset `p` in `paths` should be accessible as `fid[p]`.
Calling `getobs` on a `JLD2Dataset` is equivalent to mapping `getobs` on
each dataset in `paths`.
See [`close(::JLD2Dataset)`](@ref) for closing the underlying JLD2 file pointer.
"""
struct JLD2Dataset{T<:JLD2.JLDFile, S<:Tuple} <: AbstractDataContainer
    fid::T
    paths::S

    function JLD2Dataset(fid::JLD2.JLDFile, paths)
        _paths = Tuple(map(p -> fid[p], paths))
        nobs = map(numobs, _paths)
        _check_jld2_nobs(nobs) ||
            throw(ArgumentError("Cannot create JLD2Dataset for datasets with mismatched number of observations (got $nobs)."))

        new{typeof(fid), typeof(_paths)}(fid, _paths)
    end
end

JLD2Dataset(file::JLD2.JLDFile, path::String) = JLD2Dataset(file, (path,))
JLD2Dataset(file::AbstractString, paths) = JLD2Dataset(jldopen(file, "r"), paths)

Base.getindex(dataset::JLD2Dataset{<:JLD2.JLDFile, <:NTuple{1}}, i) = getobs(only(dataset.paths), i)
Base.getindex(dataset::JLD2Dataset, i) = map(Base.Fix2(getobs, i), dataset.paths)
Base.length(dataset::JLD2Dataset) = numobs(dataset.paths[1])

"""
    close(dataset::JLD2Dataset)

Close the underlying JLD2 file pointer for `dataset`.
"""
Base.close(dataset::JLD2Dataset) = close(dataset.fid)