text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Human1")

# ------------------------------------------------------------------
# ARGS
import ArgParse: ArgParseSettings, @add_arg_table!, parse_args

set = ArgParseSettings()
@add_arg_table! set begin
    "--no-restart"
        help = "Do not restart env"
        action = :store_true
    "--comparativeFVA"
        help = "run comparativeFVA_humanModels function"
        action = :store_true
    "--predictGRate"
        help = "run predict_cellLines_gRates script"
        action = :store_true
end

if isinteractive()
    no_restart_flag = false
    comparativeFVA_flag = false
    predictGRate_flag = false
else
    parsed_args = parse_args(set)
    no_restart_flag = parsed_args["no-restart"]
    comparativeFVA_flag = parsed_args["comparativeFVA"]
    predictGRate_flag = parsed_args["predictGRate"]
end

# --------------------------------------------------------------------

import Chemostat_Human1
const H1 = Chemostat_Human1
import Chemostat_Human1: Chemostat
const ChU = Chemostat.Utils

import MATLAB
import MATLAB: @mat_str, eval_string, @mget, mxcall, @show
import ProgressMeter: ProgressUnknown, next!, finish!

## --------------------------------------------------------------------
# Systems
const PATHS_TO_DEL = Set{String}()
atexit() do
    for path in PATHS_TO_DEL
        try; rm(path, force = true, recursive = true) catch end
    end
    try; string_eval("path('$(MATPATH0)')") catch end # restore initial path
end

# flushing buffers
@async while true
    flush.([stderr, stdout])
    sleep(5 + 10 * rand())
end

## --------------------------------------------------------------------
# testing MATLAB connection
println("Testing MATLAB connection")
try; mat"disp(['Hello from MATLAB ' version])"
    @eval string("Connection OK")
catch err
    error("ERROR CALLING MATLAB:\n"*
          "I use MATLAB.jl, see https://github.com/JuliaInterop/MATLAB.jl for help:\n"*
          ChU.err_str(err)
        )
end
# Initial path
const MATPATH0 = let
    mat"str_path = path"
    str_path = @mget str_path
    mat"clear str_path"
    str_path
end;

## --------------------------------------------------------------------
# MAT TOOLS
function get_matpath()
    mat"str_path = path"
    str_path = @mget str_path
    mat"clear str_path"
    split(str_path, ":")
end

matwhich(ider) = mxcall(:which, 1, string(ider));
matcd(dir) = (mxcall(:cd, 1, string(dir)); mat"disp(['pwd: ' pwd])")
addto_matpath(path; addsubs = true) = addsubs ?
    eval_string("path([genpath('$path') pathsep path])") :
    eval_string("path('$path', path)")
mxcall_ans(x...) = (mxcall(x...); @mget ans)
meval_ans(exp::String) = (eval_string(exp); @mget ans)

function plant_matlocator(path)
    path = endswith(path, ".m") ? path : path * ".m"
    fun_name = replace(basename(path), ".m" => "")
    src = """
    function [ dir ] = $fun_name()
        % return the directory where this function is placed
            script_file = mfilename('fullpath');
            parts = strsplit(script_file, filesep);
            dir  = strjoin(parts(1:end-1), filesep);
        end
    """
    write(path, src)
    addto_matpath(path |> dirname; addsubs = false)
    push!(PATHS_TO_DEL, path)
    @assert abspath(path |> dirname) == mxcall_ans(Symbol(fun_name), 0)
    path
end

## --------------------------------------------------------------------
# SETTING HUMAN1 DEPS
println("Setting up Human1 deps")
HERE = @__DIR__
cd(HERE); matcd(HERE)
mat"clear"

## --------------------------------------------------------------------
println("Setting up COBRA Toolbox")
cobra_dir = joinpath(H1.DEPS_DIR, "cobratoolbox")
if !isdir(cobra_dir)
    run(`git clone --depth=1 https://github.com/opencobra/cobratoolbox.git $cobra_dir`)
    @info string("cobratoolbox cloned at: ", relpath(cobra_dir))
else
    @info string("cobratoolbox found at: ", relpath(cobra_dir))
    cd(cobra_dir); run(`git checkout -- .`); cd(HERE)
end
# run initCobraToolbox
@info string("running initCobraToolbox")
addto_matpath(cobra_dir; addsubs = false)
matcd(cobra_dir)
mat"initCobraToolbox(false)"
matcd(HERE)

## --------------------------------------------------------------------
println("Setting up GECKO Toolbox")
gecko_dir = joinpath(H1.DEPS_DIR, "GECKO")
if !isdir(gecko_dir)
    # run(`git clone https://github.com/SysBioChalmers/GECKO $gecko_dir`)
    run(`git clone /Users/Pereiro/Documents/MATLAB/GECKO $gecko_dir`) # dev
    @info string("GECKO cloned at: ", relpath(gecko_dir))
else
    @info string("GECKO found at: ", relpath(gecko_dir))
    cd(gecko_dir); run(`git checkout -- .`); cd(HERE)
end
addto_matpath(gecko_dir; addsubs = true)

## --------------------------------------------------------------------
println("Setting up RAVEN Toolbox")
raven_dir = joinpath(H1.DEPS_DIR, "RAVEN")
if !isdir(raven_dir)
    run(`git clone https://github.com/SysBioChalmers/RAVEN $raven_dir`)
    # run(`git clone /Users/Pereiro/Documents/MATLAB/RAVEN $raven_dir`) # Dev
    @info string("RAVEN cloned at: ", relpath(raven_dir))
else
    @info string("RAVEN found at: ", relpath(raven_dir))
    cd(raven_dir); run(`git checkout -- .`); cd(HERE)
end
addto_matpath(raven_dir; addsubs = true)

## --------------------------------------------------------------------
# Prepare playground
temp_dir = joinpath(H1.PROJ_ROOT, "temp")
if !no_restart_flag
    push!(PATHS_TO_DEL, temp_dir)
    rm(temp_dir, force = true, recursive = true)
    cp(H1.HUMAN1_PROBLICATION_RAW_DATA_DIR, temp_dir, force = true);
    addto_matpath(temp_dir; addsubs = true)
    cd(HERE)
end
@assert isdir(temp_dir)

## --------------------------------------------------------------------
# Plant markers
plant_matlocator(joinpath(H1.HUMAN1_MODELS_DIR, "get_models_folder"))
plant_matlocator(joinpath(H1.CACHE_DIR, "get_cache_folder"))
plant_matlocator(joinpath(H1.HUMAN1_RESULTS_DIR, "get_results_folder"))
plant_matlocator(joinpath(temp_dir, "get_proj_folder"))

## --------------------------------------------------------------------
matscripts_dir = joinpath(H1.PROJ_ROOT, "scripts", "matlab_scripts")

## --------------------------------------------------------------------
# Run coparative_FVA_humanModels
if comparativeFVA_flag
    # prepare comparativeFVA_humanModels
    cd(matscripts_dir)
    try
        cp("my_comparativeFVA_humanModels", matwhich("comparativeFVA_humanModels"); force = true);
        @info "comparativeFVA_humanModels updated"
    catch; rethrow() end

    # prepare GECKO
    cd(gecko_dir); wait(run(`git checkout 45804e1`))
    cd(matscripts_dir)
    try
        cp("my_comparativeFVA", matwhich("comparativeFVA"); force = true);
        @info "comparativeFVA updated"
    catch; rethrow() end

    # prepare RAVEN
    cd(raven_dir); wait(run(`git checkout -- .`))
    cd(matscripts_dir)
    try
        cp("my_setExchangeBounds", matwhich("setExchangeBounds"); force = true);
        @info "setExchangeBounds updated"
    catch; rethrow() end
    cd(HERE)

    @assert isdir(matscripts_dir)
    # run simulation
    mat"clear"
    mat"cellLine = 'HOP62'"
    mat"comparativeFVA_humanModels(cellLine)"
end

## --------------------------------------------------------------------
# run predict cellLines
if predictGRate_flag
    # prepare GECKO
    cd(gecko_dir); wait(run(`git checkout 45804e1`))
    # prepare RAVEN
    cd(raven_dir); wait(run(`git checkout -- .`))
    cd(matscripts_dir)
    try
        cp("my_setExchangeBounds", matwhich("setExchangeBounds"); force = true);
        @info "setExchangeBounds updated"
    catch; rethrow() end
    # prepare comparativeFVA_humanModels
    try
        cp("my_predict_cellLines_gRates", matwhich("predict_cellLines_gRates"); force = true);
        @info "predict_cellLines_gRates updated"
    catch; rethrow() end
    try
        cp("my_ExchFluxesComparison_NCI60", matwhich("ExchFluxesComparison_NCI60"); force = true);
        @info "ExchFluxesComparison_NCI60 updated"
    catch; rethrow() end

    # run simulation
    sim_folder = joinpath(temp_dir, "ec_GEMs/ComplementaryScripts")
    matcd(sim_folder); cd(sim_folder)
    mat"clear"
    mat"predict_cellLines_gRates"
    cd(HERE)
end


text:
function f(n)
    while n < 10
        println(n)
        n = n+1
    end
end

f(5)

ext: jl

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
module GlobalApproximationValueIteration

using Base.Iterators: repeated

# Stdlib imports
using LinearAlgebra
using Random
using Printf

# For function approximation
using MultivariateStats
using Flux

# POMDPs imports
using POMDPs
using POMDPPolicies
using POMDPModelTools
using POMDPLinter: @POMDP_require, @warn_requirements

export
    GlobalFunctionApproximator,
    fit!,
    compute_value,
    LinearGlobalFunctionApproximator,
    NonlinearGlobalFunctionApproximator

export
    GlobalApproximationValueIterationSolver,
    GlobalApproximationValueIterationPolicy,
    convert_featurevector,
    sample_state

function sample_state end

include("global_function_approximation.jl")
include("linear_gfa.jl")
include("nonlinear_gfa.jl")
include("global_approximation_vi.jl")

end # module


text:
libmx = C_NULL
libeng = C_NULL
libmat = C_NULL

# Determine MATLAB library path and provide facilities to load libraries with
# this path

function get_paths()
    global matlab_homepath = get(ENV, "MATLAB_HOME", "")
    global default_startcmd = C_NULL
    global matlab_library_path = nothing

    if matlab_homepath == ""
        if OS_NAME == :Linux
            matlab_homepath = dirname(dirname(realpath(chomp(readall(`which matlab`)))))
        elseif OS_NAME == :Darwin
            apps = readdir("/Applications")
            filter!(app -> ismatch(r"^MATLAB_R[0-9]+[ab]\.app$", app), apps)
            if ~isempty(apps)
                matlab_homepath = joinpath("/Applications", minimum(apps))
            end
        elseif OS_NAME == :Windows
            default_dir = Int == Int32 ? "C:\\Program Files (x86)\\MATLAB" : "C:\\Program Files\\MATLAB"
            if isdir(default_dir)
                dirs = readdir(default_dir)
                filter!(dir -> ismatch(r"^R[0-9]+[ab]$", dir), dirs)
                if ~isempty(dirs)
                    matlab_homepath = joinpath(default_dir, minimum(dirs))
                end
            end
        end
    end

    if matlab_homepath == ""
        error("The MATLAB path could not be found. Set the MATLAB_HOME environmental variable to specify the MATLAB path.")
    end

    if OS_NAME != :Windows
        default_startcmd = joinpath(matlab_homepath, "bin", "matlab")
        if !isexecutable(default_startcmd)
            error("The MATLAB path is invalid. Set the MATLAB_HOME evironmental variable to the MATLAB root.")
        end
        default_startcmd = "exec $(Base.shell_escape(default_startcmd)) -nosplash"
    elseif OS_NAME == :Windows
        default_startcmd = joinpath(matlab_homepath, "bin", (Int == Int32 ? "win32" : "win64"), "MATLAB.exe")
        if !isfile(default_startcmd)
            error("The MATLAB path is invalid. Set the MATLAB_HOME evironmental variable to the MATLAB root.")
        end
        default_startcmd *= " -nosplash"
    end

    # Get path to MATLAB libraries
    matlab_library_path = nothing
    if OS_NAME == :Linux
        matlab_library_path = joinpath(matlab_homepath, "bin", (Int == Int32 ? "glnx86" : "glnxa64"))
    elseif OS_NAME == :Darwin
        matlab_library_path = joinpath(matlab_homepath, "bin", (Int == Int32 ? "maci" : "maci64"))
    elseif OS_NAME == :Windows
        matlab_library_path = joinpath(matlab_homepath, "bin", (Int == Int32 ? "win32" : "win64"))
    end

    if matlab_library_path != nothing && !isdir(matlab_library_path)
        matlab_library_path = nothing
    end
end
get_paths()

matlab_library(lib::String) =
    matlab_library_path == nothing ? lib : joinpath(matlab_library_path, lib)

# libmx (loaded when the module is imported)

function load_libmx()
    global libmx
    if libmx == C_NULL
        libmx = dlopen(matlab_library("libmx"), RTLD_GLOBAL | RTLD_LAZY)

        if libmx == C_NULL
            error("Failed to load libmx.")
        end
    end
end

load_libmx()

# libmat (loaded when the module is imported)

function load_libmat()
    global libmat
    if libmat == C_NULL
        libmat = dlopen(matlab_library("libmat"), RTLD_GLOBAL | RTLD_LAZY)

        if libmat == C_NULL
            error("Failed to load libmat.")
        end
    end
end

load_libmat()

# libeng (loaded when needed)

function load_libeng()
    global libeng
    if libeng == C_NULL
        libeng = dlopen(matlab_library("libeng"), RTLD_GLOBAL | RTLD_LAZY)
        if libeng == C_NULL
            error("Failed to load libeng.")
        end
    end
end

engfunc(fun::Symbol) = dlsym(libeng::Ptr{Void}, fun)
mxfunc(fun::Symbol) = dlsym(libmx::Ptr{Void}, fun)
matfunc(fun::Symbol) = dlsym(libmat::Ptr{Void}, fun)


text:
abstract type Objective end

"""
    empty objective
"""
struct EmptyObjective <: Objective end
objective(Z, obj::EmptyObjective, idx, T) = 0.0
objective_gradient!(∇J, Z, obj::EmptyObjective, idx, T) = nothing

"""
    multiple objectives
"""
struct MultiObjective <: Objective
    obj::Vector{Objective}
end

function objective(Z, obj::MultiObjective, idx, T)
    return sum([objective(Z, o, idx, T) for o in obj.obj])
end

function objective_gradient!(∇J, Z, obj::MultiObjective, idx, T)
    for o in obj.obj
        objective_gradient!(∇J, Z, o, idx, T)
    end
    return nothing
end

function include_objective(str::String)
    include(joinpath(pwd(), "src/objectives", str * ".jl"))
end

function include_objective(strs::Vector{String})
    for s in strs
        include_objective(s)
    end
end


text:
import Distributions
import DataFrames

println("# of cores: ", Sys.CPU_CORES, "; # of workers: ", nworkers())

# load simulation functions
@everywhere include("simulation.jl")

## PARAMETER VALUES #####################################################
tic()

# parse input arguments
p = Int(eval(parse(ARGS[1])))  # population number
r = Int(eval(parse(ARGS[2])))  # replicate number

# varying parameters ----------------------------------------------------
H = [0.1, 0.5, 0.9]                        # heritabilities                     TODO: run beetle values
ρ = [-0.9, -0.5, -0.1, 0.0, 0.1, 0.5, 0.9] # covariances                        TODO: run beetle values

# fixed parameters ------------------------------------------------------
n = 20               # number of individuals
M = [1.63, 2.74]     # phenotype means
V = [0.40, 0.35]     # total phenotypic variances                               TODO: add additional fixed variances for appendix
b = 4.05/10.0        # Beverton-Holt parameter (re-scaled for simulation)       TODO: add additional values for appendix
ngens = 20           # number of generations to simulate
bat_time = time()    # batch processing time
c = 0                # loop counter

# Vectors to store parameter combinations
vl = length(H)^2 * length(ρ)^2

P0 = Array{Any}(vl)  # popmatrix
P1 = Array{Any}(vl)  # ngens
P2 = Array{Any}(vl)  # phenotype means
P3 = Array{Any}(vl)  # variances
P4 = Array{Any}(vl)  # correlations
P5 = Array{Any}(vl)  # heritabilities
P6 = Array{Any}(vl)  # Beverton-Holt parameter
P7 = Array{Any}(vl)  # population
P8 = Array{Any}(vl)  # rep
PT = Array{Any}(vl)  # bat_time

for h1 in H
	for h2 in H
		for ρ1 in ρ
			for ρ2 in ρ

				c += 1

				# set seed so that founding population can be replicated later
				srand(plant_seed([ρ1, ρ2], [h1, h2], r, p))

				P0[c] = init_inds(n, M, V, [ρ1, ρ2], [h1, h2]) # founding pop'n
				P1[c] = ngens        # ngens
				P2[c] = M            # phenotype means
				P3[c] = V            # variances
				P4[c] = [ρ1, ρ2]     # correlations
				P5[c] = [h1, h2]     # heritabilities
				P6[c] = b            # Beverton-Holt parameter
				P7[c] = p            # population
				P8[c] = r            # rep
				PT[c] = bat_time     # bat_time
			end
		end
	end
end

## SIMULATION ##################################################################

# run all parameter combinations in parallel
R = pmap((p0,p1,p2,p3,p4,p5,p6,p7,p8,PT) ->
   runsim(p0,p1,p2,p3,p4,p5,p6,p7,p8,PT),
	        P0,P1,P2,P3,P4,P5,P6,P7,P8,PT);

## RECORD OUTPUT ###############################################################

# convert R to a single array (could be a better way to do this)
out = Array{Any}(size(R, 1) * ngens * 12, size(R[1], 2))

# fill 'out' with data from R
for i in 1:size(R, 1)
	out[((i -1) * ngens * 12 + 1) : i * ngens * 12, :] = R[i]
end

# write output
writedlm(outname(p, r), out, '\,')
toc()

max_stars_repo_name: bochocki/correlatedtraits

max_issues_repo_name: bochocki/correlatedtraits

max_forks_repo_name: bochocki/correlatedtraits


text:
using AbstractNumbers, SpecialFunctions
using Test, Aqua

Aqua.test_ambiguities([AbstractNumbers, Base, Core])
Aqua.test_unbound_args(AbstractNumbers)
Aqua.test_undefined_exports(AbstractNumbers)
Aqua.test_project_extras(AbstractNumbers)
Aqua.test_stale_deps(AbstractNumbers)
Aqua.test_deps_compat(AbstractNumbers)
Aqua.test_project_toml_formatting(AbstractNumbers)

struct MyNumber{T} <: AbstractNumbers.AbstractNumber{T}
    number::T
end
const SF = SpecialFunctions

AbstractNumbers.number(x::MyNumber) = x.number
AbstractNumbers.basetype(::Type{<: MyNumber}) = MyNumber

x = MyNumber(1)
@test convert(Number, x) === x

all_funcs = (
    :~, :conj, :abs, :sin, :cos, :tan, :sinh, :cosh, :tanh, :asin, :acos, :atan,
    :asinh, :acosh, :atanh, :sec, :csc, :cot, :asec, :acsc, :acot, :sech, :csch,
    :coth, :asech, :acsch, :acoth, :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
    :sind, :tand, :acosd, :acotd, :acscd, :asecd, :asind, :atand, :rad2deg,
    :deg2rad, :log, :log2, :log10, :log1p, :exponent, :exp, :exp2, :expm1,
    :cbrt, :sqrt, :ceil, :floor, :trunc, :round, :significand,
    :frexp, :ldexp, :modf, :real, :imag, :!, :identity,
    :zero, :one, :<<, :>>, :abs2, :sign, :sinpi, :cospi, :exp10,
    :iseven, :ispow2, :isfinite, :isinf, :isodd, :isinteger, :isreal,
    :isnan, :isempty, :iszero, :transpose, :copysign, :flipsign, :signbit,
    :+, :-, :*, :/, :\, :^, :(==), :(!=), :<, :(<=), :>, :(>=), :min, :max,
    :div, :fld, :rem, :mod, :mod1, :cmp, :&, :|, :xor,
    :clamp
)

special_funcs = (
    :airyai,
    :airyaiprime,
    :airybi,
    :airybiprime,
    :airyaix,
    :airyaiprimex,
    :airybix,
    :airybiprimex,
    :besselh,
    :besselhx,
    :besseli,
    :besselix,
    :besselj,
    :besselj0,
    :besselj1,
    :besseljx,
    :besselk,
    :besselkx,
    :bessely,
    :bessely0,
    :bessely1,
    :besselyx,
    :dawson,
    :erf,
    :erfc,
    :erfcinv,
    :erfcx,
    :erfi,
    :erfinv,
    :eta,
    :digamma,
    :invdigamma,
    :polygamma,
    :trigamma,
    :hankelh1,
    :hankelh1x,
    :hankelh2,
    :hankelh2x,
    :zeta
)

nancompare(a, b) = a == b
function nancompare(a::Number, b::Number)
    isnan(a) && isnan(b) && return true
    a == b
end

function myrand(::Type{T}) where T <: Integer
    x = rand(T(-10):T(10))
    x == 0 ? T(1) : x
end
myrand(::Type{T}) where T <: Unsigned = rand(T(1):T(20))
myrand(::Type{T}) where T = rand(T)
@testset "AbstractNumbers" begin
    for (mod, funcs) in ((Base, all_funcs),), f in funcs
        println(f)
        func = getfield(mod, f)
        @testset "testing $f" begin
            for i = 1:4
                for T in (Float32, Float64, ComplexF64, Int, UInt)
                    args = ntuple(x-> myrand(T), i)
                    if i == 3 && func === cmp
                        continue
                    end
                    if T <: Complex
                        (func in (
                            cotd, cosd, cscd, secd, sind, tand, acosd, acotd, acscd,
                            asecd, asind, atand, rem, modf, (<), (>), (<=), (>=),
                            min, max, cmp, &, |, xor, clamp, div, fld
                        )) && continue
                    end
                    if T <: AbstractFloat
                        (func in (
                            &, |, xor
                        )) && continue
                    end
                    # i have no idea what these functions are - seem to throw exceptions when not careful with input values
                    # also, they're not usable since they're deprecated, even if you get them from specialfunctions
                    if applicable(func, args...)
                        res = try
                            a = func(args...)
                            b = func(MyNumber.(args)...)
                            @testset "$T $args" begin
                                @test nancompare(a, b)
                            end
                        catch e
                            isa(e, DomainError) || rethrow(e)
                        end
                    end
                end
            end
        end
    end

    # isapprox
    for T in (Float32, Float64, ComplexF64, Int, UInt)
        T = Float64
        args = ntuple(x-> myrand(T), 2)
        a = ≈(args...)
        b = ≈(MyNumber.(args)...)
        a == b
    end

    @testset "bessel functions" begin
        bessel_funcs = [(SF.bessely0, SF.bessely1, SF.bessely), (SF.besselj0, SF.besselj1, SF.besselj)]
        @testset "$z, $o" for (z, o, f) in bessel_funcs
            @test z(Float32(2.0)) ≈ z(Float64(2.0))
            @test o(Float32(2.0)) ≈ o(Float64(2.0))
            @test z(MyNumber(2)) ≈ z(MyNumber(2.0))
            @test o(MyNumber(2)) ≈ o(MyNumber(2.0))
            @test z(MyNumber(2.0 + im)) ≈ f(MyNumber(0), MyNumber(2.0 + im))
            @test o(MyNumber(2.0 + im)) ≈ f(MyNumber(1), MyNumber(2.0 + im))
        end
        @testset "besselj error throwing" begin
            @test_throws MethodError SF.besselj(MyNumber(1.2), MyNumber(big(1.0)))
            @test_throws MethodError SF.besselj(MyNumber(1), MyNumber(complex(big(1.0))))
            @test_throws MethodError SF.besseljx(MyNumber(1), MyNumber(big(1.0)))
            @test_throws MethodError SF.besseljx(MyNumber(1), MyNumber(complex(big(1.0))))
        end
    end

end

nothing


text:
struct FlattenTree

    feature::Vector{Int}
    threshold::Vector{Float32}
    left::Vector{Int}
    right::Vector{Int}

end

function tree_assign_ids!(node, id=1)

    node.id = id

    if !node.is_leaf
        if node.split.n_left >= node.split.n_right
            id = tree_assign_ids!(node.left, id+1)
            id = tree_assign_ids!(node.right, id+1)
        else
            id = tree_assign_ids!(node.right, id+1)
            id = tree_assign_ids!(node.left, id+1)
        end
    end

    return id

end

function tree_converter(root)

    n = tree_assign_ids!(root)

    feature = Vector{Int}(undef, n)
    threshold = Vector{Float32}(undef, n)
    left = Vector{Int}(undef, n)
    right = Vector{Int}(undef, n)

    stack = Node[root]
    while length(stack) > 0
        node = pop!(stack)
        if node.is_leaf
            feature[node.id] = -1
            threshold[node.id] = node.probability
        else
            feature[node.id] = node.split.feature
            threshold[node.id] = node.split.threshold
            left[node.id] = node.left.id
            right[node.id] = node.right.id
            push!(stack, node.left, node.right)
        end
    end

    return FlattenTree(feature, threshold, left, right)

end

struct FlattenModel

    trees::Vector{FlattenTree}
    options::Options
    metadata::Metadata
    oob_samples::Vector{BitArray{1}}

    function FlattenModel(model::Model)

        return new([tree_converter(tree) for tree in model.trees], model.options, model.metadata, model.oob_samples)

    end

end

FlattenModel(X, Y;
    n_trees::Int=1,
    n_subfeat::Int=0,
    n_thresholds::Int=1,
    max_depth::Int=-1,
    min_samples_leaf::Int=1,
    min_samples_split::Int=2,
    rand_thresholds::Bool=true,
    soft_leaf_score::Bool=true,
    bagging::Float64=0.0,
    extra_pos_bagging::Float64=0.0,
    exclude_samples::Set{Int}=Set{Int}(),
    description::String="none"
    ) =
    FlattenModel(
        Model(X, Y,
            Options(
                n_trees,
                (1 <= n_subfeat <= size(X, 2)) ? n_subfeat : round(Int, sqrt(size(X, 2))),
                n_thresholds,
                max_depth,
                min_samples_leaf,
                min_samples_split,
                rand_thresholds,
                soft_leaf_score,
                bagging,
                extra_pos_bagging,
                exclude_samples
                ),
            description
        )
    )

Base.show(io::IO, ::MIME"text/plain", model::FlattenModel) = print(io,
    "Model: ExtraTrees [flatten array-based representation]\n",
    "Trees: ", length(model.trees), "\n",
    "Data features: ", model.metadata.n_features, "\n",
    "Average depth: ", round(model.metadata.avg_tree_depth, digits=1), "\n",
    "Training time: ", round(model.metadata.trainingtime, digits=2), " sec\n",
    "Training pos samples: ", model.metadata.n_pos_samples, "\n",
    "Training neg samples: ", model.metadata.n_neg_samples, "\n",
    "Custom description: ", model.metadata.description, "\n"
    )

function (tree::FlattenTree)(sample::Vector{Float32})

    id = 1
    while tree.feature[id] > 0
        id = sample[tree.feature[id]] < tree.threshold[id] ? tree.left[id] : tree.right[id]
    end

    return tree.threshold[id]

end

function (tree::FlattenTree)(sample::Vector{Float32}, level::Bool)

    id = 1
    level = 1
    while tree.feature[id] > 0
        id = sample[tree.feature[id]] < tree.threshold[id] ? tree.left[id] : tree.right[id]
        level += 1
    end

    return (tree.threshold[id], level)

end

function (model::FlattenModel)(sample::Vector{Float32})

    @assert model.metadata.n_features == length(sample)

    probability = 0.0f0

    for tree in model.trees
        probability += tree(sample)
    end

    return probability / length(model.trees)

end

function (model::FlattenModel)(X::Matrix{Float32})

    @assert model.metadata.n_features == size(X, 2)

    return [model(X[ss, :]) for ss in 1:size(X, 1)]

end

function (model::FlattenModel)(XS::SharedArray{Float32,2}, range::UnitRange{Int64})

    X = [XS[ss, :] for ss in range]

    probability = zeros(Float32, length(X))

    for tree in model.trees
        for ss in 1:length(X)
            probability[ss] += tree(X[ss])
        end
    end

    return probability ./ length(model.trees)

end

function oob_model(model::FlattenModel, XS::SharedArray{Float32,2}, range::UnitRange{Int64})

    X = [XS[ss, :] for ss in range]

    scores = zeros(Float32, length(X))
    tree_count = zeros(Int, length(X))

    for (tree, oob_samples) in zip(model.trees, model.oob_samples)
        for xx in 1:length(X)
            if oob_samples[range[xx]]
                scores[xx] += tree(X[xx])
                tree_count[xx] += 1
            end
        end
    end

    return scores ./ tree_count

end

function (model::FlattenModel)(XS::SharedArray{Float32,2}, range::UnitRange{Int64}, oob_samples::Bool)

    return oob_samples ? oob_model(model, XS, range) : model(XS, range)

end

function (model::FlattenModel)(XS::SharedArray{Float32,2}, oob_samples::Bool=false)

    @assert model.metadata.n_features == size(XS, 2)

    if oob_samples
        @assert length(model.oob_samples[1]) == size(XS, 1)
    end

    n_samples = size(XS, 1)
    n_workers = length(workers())
    step_size = ceil(Int, n_samples / n_workers)

    ranges = [ start:((stop=start+step_size-1) < n_samples ? stop : n_samples) for start in 1:step_size:n_samples ]

    results = pmap(range -> (first(range), model(XS, range, oob_samples)), ranges)

    return vcat(last.(sort(results, by=x->x[1]))...)

end


text:
include("monitor.jl")

containerID, containerStr, netIO = stats()

taintedFlag = map(containerID, containerStr) do x,y
	z = String.(split(y, '-'))
	println(z[1])
	hash = open(f->read(f, String), z[1])
	hashlog(x) != hash
end


text:
export saveRecoParams, loadRecoParams, defaultRecoParams, defaultOnlineRecoParams

function defaultRecoParams()
  params = Dict{Symbol,Any}()
  params[:lambd] = 1e-2
  params[:iterations] = 4
  params[:SNRThresh] = 2.0
  params[:minFreq] = 80e3
  params[:maxFreq] = 1.25e6
  params[:sortBySNR] = false
  params[:nAverages] = 1
  params[:repetitionTime] = 0.0
  params[:denoiseWeight] = 0
  params[:loadas32bit] = true
  params[:loadasreal] = false
  params[:sparseTrafo] = nothing
  params[:redFactor] = 0.0
  params[:solver] = "kaczmarz"
  params[:emptyMeasPath] = nothing
  params[:frames] = 1
  params[:spectralCleaning] = true
  params[:recChannels] = [1,2,3]
  params[:reconstructor] = get(ENV,"USER","default")
  params[:firstFrameBG] = 1
  params[:lastFrameBG] = 1

  return params
end

function defaultRecoParamsOld()
  params = Dict{Symbol,Any}()
  params[:lambd] = "1e-2"
  params[:iterations] = "4"
  params[:SNRThresh] = "2"
  params[:minFreq] = "80e3"
  params[:maxFreq] = "1.25e6"
  params[:sortBySNR] = "false"
  params[:nAverages] = "1"
  params[:repetitionTime] = "0.0"
  params[:denoiseWeight] = "0"
  params[:loadas32bit] = "true"
  params[:loadasreal] = "false"
  params[:maxload] = "100"
  params[:sparseTrafo] = "nothing"
  params[:redFactor] = "0.0"
  params[:solver] = "kaczmarz"
  params[:bEmpty] = "nothing"

  return params
end

function defaultOnlineRecoParams()
  params = defaultRecoParams()
  params[:iterations] = 1
  params[:SNRThresh] = 5

  return params
end

function saveRecoParams(filename::AbstractString, params)
  ini = Inifile()
  for (key,value) in params
    set(ini, string(key), string(value) )
  end
  open(filename,"w") do fd
    write(fd, ini)
  end
end

to_bool(s::AbstractString) = (lowercase(s) == "true") ? true : false
to_bool(b::Bool) = b

function loadRecoParams(filename::AbstractString)
  ini = Inifile()

  if isfile(filename)
    read(ini, filename)
  end

  params = defaultRecoParamsOld()

  for key in [:lambd, :SNRThresh, :minFreq, :maxFreq, :repetitionTime, :denoiseWeight, :redFactor]
    params[key] = parse(Float64,get(ini,"","$key", params[key]))
  end

  for key in [:iterations, :nAverages, :maxload]
    params[key] = parse(Int,get(ini,"","$key", params[key]))
  end

  for key in [:sortBySNR, :loadas32bit, :loadasreal]
    params[key] = to_bool(get(ini,"","$key", params[key]))
  end

  sparseTrafo = get(ini,"","sparseTrafo", params[:sparseTrafo])
  params[:sparseTrafo] = (sparseTrafo == "nothing") ? nothing : sparseTrafo
  params[:solver] = get(ini,"","solver", params[:solver])

  params[:bEmpty] = get(ini,"bEmpty")
  if params[:bEmpty] == :notfound
    params[:bEmpty] = nothing
  end

  params[:SFPath] = get(ini,"SFPath")
  if params[:SFPath] == :notfound
    params[:SFPath] = nothing
  else
    params[:SFPath] = String[strip(path) for path in split(params[:SFPath],",")]
  end

  params[:SFPathFreq] = get(ini,"SFPathFreq")
  if params[:SFPathFreq] == :notfound
    params[:SFPathFreq] = nothing
  else
    params[:SFPathFreq] = String[strip(path) for path in split(params[:SFPathFreq],",")]
  end

  recChanStr = get(ini,"recChannels")
  if recChanStr != :notfound
    params[:recChannels] = [parse(Int,chan) for chan in split(recChanStr,",")]
  end

  return params
end


text:
"""
    ZeroOneBridge{T}

The `ZeroOneBridge` splits a `MOI.SingleVariable`-in-`MOI.ZeroOne` constraint
into a `MOI.SingleVariable`-in-`MOI.Integer` constraint
and a `MOI.SingleVariable`-in-`MOI.Interval(0, 1)` constraint.
"""
struct ZeroOneBridge{T} <: AbstractBridge
    interval_index::MOI.ConstraintIndex{MOI.SingleVariable,MOI.Interval{T}}
    integer_index::MOI.ConstraintIndex{MOI.SingleVariable,MOI.Integer}
end

function bridge_constraint(
    ::Type{ZeroOneBridge{T}},
    model::MOI.ModelLike,
    f::MOI.SingleVariable,
    ::MOI.ZeroOne,
) where {T<:Real}
    interval_index =
        MOI.add_constraint(model, f, MOI.Interval{T}(zero(T), one(T)))
    integer_index = MOI.add_constraint(model, f, MOI.Integer())
    return ZeroOneBridge{T}(interval_index, integer_index)
end

function MOIB.added_constraint_types(::Type{<:ZeroOneBridge{T}}) where {T}
    return [
        (MOI.SingleVariable, MOI.Interval{T}),
        (MOI.SingleVariable, MOI.Integer),
    ]
end

function MOIB.added_constrained_variable_types(::Type{<:ZeroOneBridge})
    return Tuple{DataType}[]
end

function concrete_bridge_type(
    ::Type{<:ZeroOneBridge{T}},
    ::Type{MOI.SingleVariable},
    ::Type{MOI.ZeroOne},
) where {T}
    return ZeroOneBridge{T}
end

function MOI.supports_constraint(
    ::Type{<:ZeroOneBridge},
    ::Type{MOI.SingleVariable},
    ::Type{MOI.ZeroOne},
)
    return true
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::ZeroOneBridge,
)
    return MOI.ZeroOne()
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::ZeroOneBridge,
)
    return MOI.get(model, attr, bridge.interval_index)
end

function MOI.delete(model::MOI.ModelLike, bridge::ZeroOneBridge)
    MOI.delete(model, bridge.interval_index)
    return MOI.delete(model, bridge.integer_index)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::ZeroOneBridge,
)
    return MOI.get(model, attr, bridge.interval_index)
end

function MOI.supports(
    ::MOI.ModelLike,
    ::MOI.ConstraintPrimalStart,
    ::Type{<:ZeroOneBridge},
)
    return true
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::ZeroOneBridge{T},
    value,
) where {T}
    MOI.set(model, attr, bridge.integer_index, value)
    return MOI.set(model, attr, bridge.interval_index, value)
end

# Attributes, Bridge acting as a model
function MOI.get(
    bridge::ZeroOneBridge{T},
    ::MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Interval{T}},
) where {T}
    return 1
end

function MOI.get(
    bridge::ZeroOneBridge,
    ::MOI.NumberOfConstraints{MOI.SingleVariable,MOI.Integer},
)
    return 1
end

function MOI.get(
    bridge::ZeroOneBridge,
    ::MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Interval{T}},
) where {T}
    return [bridge.interval_index]
end

function MOI.get(
    bridge::ZeroOneBridge,
    ::MOI.ListOfConstraintIndices{MOI.SingleVariable,MOI.Integer},
)
    return [bridge.integer_index]
end


text:
module pretrained
import ..pretrained_dir
include(joinpath(pretrained_dir, "v2_tiny", "voc2007", "v2_tiny_voc.jl"))
include(joinpath(pretrained_dir, "v2", "voc2007", "v2_voc.jl"))
end #module


text:
using BinDeps
using Compat

@BinDeps.setup

const libsilo = library_dependency("libsilo", aliases=["libsilo"])
const silover = "4.9.1"
const silofilebase = "silo-$(silover)-bsd"

provides(Sources,
         URI("https://wci.llnl.gov/content/assets/docs/simulation/computer-codes/silo/silo-$(silover)/$(silofilebase).tar.gz"),
         libsilo)

println(BinDeps.depsdir(libsilo))

prefix = joinpath(BinDeps.depsdir(libsilo), "usr")
srcdir = joinpath(BinDeps.depsdir(libsilo), "src", silofilebase)
println("Installing Silo source to ", srcdir)

provides(SimpleBuild, (@build_steps begin
                         GetSources(libsilo)
                         CreateDirectory(srcdir)
                         @build_steps begin
                           ChangeDirectory(srcdir)
                           `./configure --disable-silex --disable-browser --prefix=$prefix --enable-shared`
                           `make`
                           `make install`
                         end
                       end), libsilo)

@BinDeps.install @compat Dict(:libsilo => :libsilo)


text:
@doc """
### linearMap

Take in one coordinate in (x,y,z) space and convert it to parametric (s,t,u)
space. The symbols used in this function have been taken from the paper
"Free Form Deformation of Solid Geometric Models, Sederberg & Parry, 1986". It
may change at a later date

**Arguments**

*  `map` : Object of mapping type
*  `box` : BoundingBox object
*  `X`   : Point coordinate in (x,y,z) space
*  `pX`   : corresponding coordinates in (s,t,u) space

"""->

function linearMap(map::AbstractMappingType, box::AbstractBoundingBox,
                   X, pX)

  # The assumption of the mapping for the bounding box presently is that
  # coordinate transformation from physical x,y,z coordinate transformation to
  # parametric s,t,u coordinate involved only translation and scaling. There is
  # no rotation.

  # Get the x,y,z coordinates for the origin of the s,t,u system
  origin = box.origin
  S = box.unitVector[:,1]*box.lengths[1]
  T = box.unitVector[:,2]*box.lengths[2]
  U = box.unitVector[:,3]*box.lengths[3]

  XmX0 = X - origin

  # calculate s
  # Division by lengths[i] to ensure s,t,u lie between [0,1]
  TcrossU = cross(T,U)
  s = dot(TcrossU,XmX0)/dot(TcrossU,S)

  # Calculate t
  ScrossU = cross(S,U)
  t = dot(ScrossU,XmX0)/dot(ScrossU,T)

  # calculate u
  ScrossT = cross(S,T)
  u = dot(ScrossT,XmX0)/dot(ScrossT,U)

  pX[:] = [s,t,u]

  return nothing
end

@doc """
### nonlinearMap

Computes the (s,t,u) parametric coordinates of the a point in the geometry
embedded in the FFD box of any aribitrary shape. This is done by doing a Newton
solve.

**Arguments**

*  `map` : Object of Mapping type
*  `box` : bounding box object. In this case the bounding box can be any shape
*  `X`   : (x,y,z) coordinates of the point in geometry
*  `pX`  : (s,t,u) coordinates of the point. An initial guess of pX must be
           supplied. This is needed by the Newton's solve

"""->

function nonlinearMap(map::AbstractMappingType{Tffd},
                      box::AbstractBoundingBox{Tffd}, X, pX) where Tffd

  origin = box.origin

  # Compute the residual
  res = zeros(Tffd, box.ndim)
  pointVal = zeros(Tffd, box.ndim)
  xi = zeros(Tffd, box.ndim)
  xi_new = zeros(xi)
  xi[:] = pX

  # Do the newton solve to get the (s,t,u coordinates)
  for itr = 1:50
    # Compute residual
    fill!(pointVal, 0.0)
    evalVolumePoint(map, xi, pointVal)
    res = X - pointVal
    # Construct jacobian
    J = zeros(box.ndim, box.ndim)
    jderiv = zeros(Int, box.ndim)

    for i = 1:box.ndim
      fill!(jderiv, 0)
      jderiv[i] = 1
      Jrow =sview(J,i,:)
      calcdXdxi(map, xi, jderiv, Jrow)
    end
    xi_new = xi + J\res
    if norm(xi_new - xi, 2) < 1e-15
      pX[:] = xi_new[:]
      break
    else
      xi[:] = xi_new[:]
    end

  end

  return nothing
end

@doc """
### calcParametricMappingLinear

Creates a linear mapping for an array of nodes in the (x,y,z) space to the
(s,t,u) space.

**Arguments**

*  `map` : Object of Mapping type
*  `box` : BoundingBox object
*  `nodes_xyz` : (x,y,z) coordinates of the nodes of the mesh
"""->

function calcParametricMappingLinear(map::Mapping, box,
                                     nodes_xyz::AbstractArray{AbstractFloat,4})

  X = zeros(map.ndim)
  for k = 1:map.numnodes[3]
    for j = 1:map.numnodes[2]
      for i = 1:map.numnodes[1]
        X[:] = nodes_xyz[i,j,k,:]
        pX =sview(map.xi,i,j,k,:)
        linearMap(map, box, X, pX)
      end
    end
  end

  return nothing
end  # End function calcParametricLinear

function calcParametricMappingLinear(map::PumiMapping{Tffd},
                               box::PumiBoundingBox, mesh::AbstractCGMesh) where Tffd

  @assert false "CG meshes not supported"

  if mesh.dim == 2
    X = zeros(Tffd,3)
    for i = 1:mesh.numEl
      for j = 1:mesh.numNodesPerElement
        X[1:2] = mesh.coords[:,j,i]
        pX =sview(map.xi,:,j,i)
        linearMap(map, box, X, pX)
      end
    end
  else
    for i = 1:mesh.numEl
      for j = 1:mesh.numNodesPerElement
        X =sview(mesh.coords,:,j,i)
        pX =sview(map.xi,:,j,i)
        linearMap(map, box, X, pX)
      end
    end
  end

  return nothing
end

function calcParametricMappingLinear(map::PumiMapping{Tffd},
                               box::PumiBoundingBox, mesh::AbstractDGMesh) where Tffd

  @assert false "this function is broken"
  if mesh.dim == 2
    X = zeros(Tffd,3)
    for i = 1:mesh.numEl
      for j = 1:size(mesh.vert_coords,2) # 1:mesh.numNodesPerElement
        X[1:2] = mesh.vert_coords[:,j,i]
        pX =sview(map.xi,:,j,i)
        linearMap(map, box, X, pX)
      end
    end
  else
    for i = 1:mesh.numEl
      for j = 1:size(mesh.vert_coords,2) # 1:mesh.numNodesPerElement
        X =sview(mesh.vert_coords,:,j,i)
        pX =sview(map.xi,:,j,i)
        linearMap(map, box, X, pX)
      end
    end
  end

  return nothing
end

function calcParametricMappingLinear(map::PumiMapping{Tffd},
                               box::PumiBoundingBox, mesh::AbstractCGMesh,
                               bc_nums::AbstractArray{Int,1}) where Tffd

  @assert false "CG meshes not supported"
  if mesh.dim == 2
    x = zeros(Tffd,3)
    for (idx, itr) in enumerate(bc_nums)
      start_index = mesh.bndry_offsets[itr]
      end_index = mesh.bndry_offsets[itr+1]
      idx_range = start_index:end_index
      bndry_facenums = sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          fill!(x, 0.0)
          x[1:2] = mesh.coords[:,vtx_arr[j],bndry_i.element]
          pX = sview(map.xi[idx], :, j, i)
          linearMap(map, box, x, pX)
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)
  else
    for (idx, itr) in enumerate(bc_nums)
      start_index = mesh.bndry_offsets[itr]
      end_index = mesh.bndry_offsets[itr+1]
      idx_range = start_index:end_index
      bndry_facenums = sview(mesh.bndryfaces, start_index:(end_index - 1))
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          fill!(x, 0.0)
          X =sview(mesh.coords,:,vtx_arr,bndry_i.elements)
          pX =sview(map.xi[idx], :, j, i)
          linearMap(map, box, X, pX)
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)

  end  # End if mesh.dim == 2

  return nothing
end

function calcParametricMappingLinear(map::PumiMapping{Tffd},
                               box::PumiBoundingBox, mesh::AbstractDGMesh,
                               bc_nums::AbstractArray{Int,1}) where Tffd

  @assert false "This function is broken"
  # Check if the knot vectors are for Bezier Curves with Bernstein polynomial
  # basis functions
  for i = 1:length(map.edge_knot)
    ctr = 0
    for j = 1:length(map.edge_knot[i])
      if map.edge_knot[i][j] != 0.0 && map.edge_knot[i][j] != 1.0
        ctr += 1
      end
    end
    @assert ctr == 0 "Linear mapping works only for Bezier Curves with Bernstein polynomial basis functions"
  end

  x = zeros(Tffd,3)
  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      # get the local index of the vertices
      vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
      for j = 1:length(vtx_arr)
        fill!(x, 0.0)
        for k = 1:mesh.dim
          x[k] = mesh.vert_coords[k,vtx_arr[j],bndry_i.element]
        end
        pX =sview(map.xi[idx], :, j, i)
        linearMap(map, box, x, pX)
      end  # End for j = 1:length(vtx_arr)
    end    # End for i = 1:nfaces
  end      # End for itr = 1:length(geomfaces)

  return nothing
end

@doc """
### calcParametricMappingNonlinear

Creates a non linear mapping for an array of nodes in the (x,y,z) space to the
(s,t,u) space.

**Arguments**

*  `map` : Object of Mapping type
*  `box` : BoundingBox object
*  `nodes_xyz` : (x,y,z) coordinates of the points in the embedded geometry

"""->

function calcParametricMappingNonlinear(map::Mapping, box,
                                        nodes_xyz::AbstractArray{AbstractFloat,4})

  X = zeros(map.ndim)
  pX = zeros(map.ndim)
  for k = 1:map.numnodes[3]
    for j = 1:map.numnodes[2]
      for i = 1:map.numnodes[1]
        X[:] = nodes_xyz[i,j,k,:]
        pX[:] = [0.,0.,0.]
        nonlinearMap(map, box, X, pX)
        map.xi[i,j,k,:] = pX[:]
      end
    end
  end

  return nothing
end

# unused now that full_geom is no more?
function calcParametricMappingNonlinear(map::PumiMapping{Tffd},
                                  box::PumiBoundingBox, mesh::AbstractCGMesh) where Tffd

  @assert false "CG meshes not supported"
  if mesh.dim == 2
    X = zeros(Tffd,3)
    for i = 1:mesh.numEl
      for j = 1:mesh.numNodesPerElement
        X[1:2] = mesh.coords[:,j,i]
        pX =sview(map.xi,:,j,i)
        nonlinearMap(map, box, X, pX)
      end
    end
  else
    for i = 1:mesh.numEl
      for j = 1:mesh.numNodesPerElement
        X =sview(mesh.coords,:,j,i)
        pX =sview(map.xi,:,j,i)
        nonlinearMap(map, box, X, pX)
      end
    end
  end

  return nothing
end

# unused now that full_geom is no more?
function calcParametricMappingNonlinear(map::PumiMapping{Tffd},
                                  box::PumiBoundingBox, mesh::AbstractDGMesh) where Tffd

  @assert false "This function is broken"
  X = zeros(Tffd,3)
  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      for k = 1:mesh.dim
        X[k] = mesh.vert_coords[k,j,i]
      end
      pX =sview(map.xi,:,j,i)
      nonlinearMap(map, box, X, pX)
    end
  end

  return nothing
end

function calcParametricMappingNonlinear(map::PumiMapping{Tffd},
                                     box::PumiBoundingBox, mesh::AbstractCGMesh,
                                     bc_nums::AbstractArray{Int,1}) where Tffd
  @assert false "CG Meshes not supported"
  if mesh.dim == 2
    x = zeros(Tffd,3)
    for (idx, itr) in enumerate(bc_nums)
      start_index = mesh.bndry_offsets[itr]
      end_index = mesh.bndry_offsets[itr+1]
      idx_range = start_index:(end_index-1)
      bndry_facenums =sview(mesh.bndryfaces, idx_range)
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          fill!(x, 0.0)
          x[1:2] = mesh.coords[:,vtx_arr[j],bndry_i.element]
          pX =sview(map.xi[idx], :, j, i)
          nonlinearMap(map, box, x, pX)
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)
  else
    for (idx, itr) in enumerate(bc_nums)
      start_index = mesh.bndry_offsets[itr]
      end_index = mesh.bndry_offsets[itr+1]
      idx_range = start_index:(end_index-1)
      bndry_facenums =sview(mesh.bndryfaces, idx_range)
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        for j = 1:length(vtx_arr)
          X =sview(mesh.coords,:,vtx_arr[j],bndry_i.element)
          pX =sview(map.xi, :, j, i)
          nonlinearMap(map, box, X, pX)
        end  # End for j = 1:length(vtx_arr)
      end    # End for i = 1:nfaces
    end      # End for itr = 1:length(geomfaces)

  end  # End if mesh.dim == 2

  return nothing
end

function calcParametricMappingNonlinear(map::PumiMapping{Tffd},
                                     box::PumiBoundingBox, mesh::AbstractDGMesh,
                                     bc_nums::AbstractArray{Int,1}) where Tffd
# this works for quadratic coordinate fields

  x_real = zeros(Float64, 3)
  for i=1:map.numFacePts
    v_i = map.face_verts[i]

    #apf.getPoint(mesh.m_ptr, v_i, 0, x_real)
    getCoords(mesh, v_i, 0, x_real)
    pX = sview(map.xi, :, i)
    nonlinearMap(map, box, sview(x_real, :), pX)
  end
#=
  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
    nfaces = length(bndry_facenums)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      # get the local index of the vertices
      vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
      for j = 1:length(vtx_arr)
        fill!(x, 0.0)
        for k = 1:mesh.dim
          x[k] = mesh.vert_coords[k,vtx_arr[j],bndry_i.element]
        end
        pX =sview(map.xi[idx], :, j, i)
        nonlinearMap(map, box, x, pX)
      end  # End for j = 1:length(vtx_arr)
    end    # End for i = 1:nfaces
  end      # End for itr = 1:length(geomfaces)

  return nothing
end
=#
  return nothing
end
@doc """
### calcdXdxi

It calculates the partial derivative of the mapping (including the 0th order)

(x,y,z) = [x(xi,eta,zeta), y(xi,eta,zeta), z(xi,eta,zeta)]

**Arguments**

*  `map`    : Object of Mapping type
*  `xi`     : Mapping parametric coordinates (s,t,u) to be evaluated
*  `jderiv` : Derivative indices. jderi[mdi] = n means take the nth derivative
              in the xi[mdi] direction.
*  `dX`     : Derivative of X w.r.t the 3 parametric variables. length = 3

REFERENCE: Carl de Boor, 'Pratical Guide to Splines', pgs 138-149, function
           BVALUE

Notes: (taken from Mapping_Mod.f90)
1) de Boor points out (pg 149) that for many points (as we may
   have) it is faster to compute the piecewise polys and differentiate
   them.  Something to consider for the future.
2) Since direction indices are used for both the physical and
   mathematical coordinates, di will be reserved for the physical
   and mdi for the mathematical coordinates in this function.

"""->

function calcdXdxi(map, xi, jderiv, dX)

  # initialize intermediate arrays
  km1 = zeros(Int, 3)
  jcmax = zeros(Int, 3)
  jcmin = zeros(Int, 3)
  imk = zeros(Int, 3)
  left = zeros(Int, 3)

  # the derivative is zero if any of jderiv(mdi) >= order(mdi)
  for mdi = 1:3
    if jderiv[mdi] >= map.order[mdi]
      return
    end
  end

  # find the spatially varying knot vector
  # calcKnot(map)

  # find the left(3) array such that
  #   knot(mdi,left(mdi)) <= xi(mdi) <= knot(mdi,left(mdi)+1)
  for mdi = 1:3
    left[mdi] = findSpan(xi[mdi], map.edge_knot[mdi], map.order[mdi], map.nctl[mdi])
  end
  # the calculations involving the knot vectors are independent,
  # so loop through them sequentially

    for mdi = 1:3

      # we will store the (order) b-spline coefficients relevant
      # to the knot interval [knot(mdi,left),knot(mdi,left+1)]
      # in aj[mdi,1],...,aj[mdi,order]
      # and compute dl[mdi,j] = xi[mdi] - knot[mdi][left-j]
      #             dr[mdi,j] = knot[mdi][left+j] - xi[mdi]
      # for all j = 1,...,order[mdi]-1

      km1[mdi] = map.order[mdi] - 1
      jcmin[mdi] = 1
      imk[mdi] = left[mdi] - map.order[mdi]
      if imk[mdi] < 0
        # we are close to the left end of the knot interval, so some
        # of the aj will be set to zero later
        jcmin[mdi] = 1 - imk[mdi]
        for j = 1:left[mdi]
          map.dl[j,mdi] = xi[mdi] - map.edge_knot[mdi][left[mdi]+1-j]
        end
        for j = left[mdi]:km1[mdi]
          map.dl[j,mdi] = map.dl[left[mdi],mdi]
        end
      else
        for j = 1:km1[mdi]
          map.dl[j,mdi] = xi[mdi] - map.edge_knot[mdi][left[mdi]+1-j]
        end
      end

      jcmax[mdi] = map.order[mdi]
      n = map.nctl[mdi]
      nmi = n - left[mdi]
      if nmi < 0
        # we are close to the right end of the knot interval, so some
        # of the aj will be set to zero later
        jcmax[mdi] = map.order[mdi] + nmi
        for j = 1:jcmax[mdi]
          map.dr[j,mdi] = map.edge_knot[mdi][left[mdi]+j] - xi[mdi]
        end
        for j = jcmax[mdi]:km[mdi]
          map.dr[j,mdi] = map.dr[jcmax[mdi],mdi]
        end
      else
        for j = 1:km1[mdi]
          map.dr[j,mdi] = map.edge_knot[mdi][left[mdi]+j] - xi[mdi]
        end
      end
    end # mdi loop
    # set all elements of aj[:,:,1] to zero, in case we are close to
    # the ends of the knot vector
    map.aj[:,:,1] = 0.0

    for jc1 = jcmin[1]:jcmax[1]
      p = imk[1] + jc1

      # set all the elements of aj[:,:,2] to zero, in case we are
      # close to the ends of the knot vector
      map.aj[:,:,2] = 0.0

      for jc2 = jcmin[2]:jcmax[2]
        q = imk[2] + jc2

        # set all the elements of aj[:,:,3] to zero, in case we are
        # close to the ends of the knot vector
        map.aj[:,:,3] = 0.0

        for jc3 = jcmin[3]:jcmax[3]
          map.aj[:,jc3,3] = map.cp_xyz[:,p,q,imk[3]+jc3]
        end

        if jderiv[3] != 0
          # derivative: apply the recursive formula X.12b from de Boor
          for j = 1:jderiv[3]
            kmj = map.order[3] - j
            ilo = kmj
            for jj = 1:kmj
              map.aj[:,jj,3] = (map.aj[:,jj+1,3] - map.aj[:,jj,3]) *
                               kmj / (map.dl[ilo,3] + map.dr[jj,3])
              ilo = ilo - 1
            end
          end
        end
        #println("map.aj[:,:,3] = \n", map.aj[:,:,3])
        if jderiv[3] != km1[3]
          # if jderiv != order - 1, we need to apply the recursive
          # formula from de Boor
          for j = jderiv[3]+1:km1[3]
            kmj = map.order[3] - j
            ilo = kmj
            for jj = 1:kmj
              map.aj[:,jj,3] = (map.aj[:,jj+1,3]*map.dl[ilo,3] +
                   map.aj[:,jj,3]*map.dr[jj,3]) / (map.dl[ilo,3] + map.dr[jj,3])
              ilo = ilo - 1
            end
          end
        end
        # println("map.aj[:,:,3] = \n", map.aj[:,:,3])

        map.aj[:,jc2,2] = map.aj[:,1,3]
      end # get next element of aj(:,2)

      if jderiv[2] != 0
        # derivative: apply the recursive formula X.12b from de Boor
        for j = 1:jderiv[2]
          kmj = map.order[2] - j
          ilo = kmj
          for jj = 1:kmj
            map.aj[:,jj,2] = (map.aj[:,jj+1,2] - map.aj[:,jj,2]) * kmj /
                             (map.dl[ilo,2] + map.dr[jj,2])
            ilo = ilo - 1
          end
        end
      end

      if jderiv[2] != km1[2]
        # if jderiv /= order - 1, we need to apply the recursive
        # formula from de Boor
        for j = jderiv[2]+1:km1[2]
          kmj = map.order[2] - j
          ilo = kmj
          for jj = 1:kmj
            map.aj[:,jj,2] = (map.aj[:,jj+1,2]*map.dl[ilo,2] + map.aj[:,jj,2]*
                             map.dr[jj,2]) / (map.dl[ilo,2] + map.dr[jj,2])
            ilo = ilo - 1
          end
        end
      end

      map.aj[:,jc1,1] = map.aj[:,1,2]
      # println("map.aj[:,:,2] = \n", map.aj[:,:,2])

    end # get next element of map.aj[:,:,1]

    if jderiv[1] != 0
      # derivative: apply the recursive formula X.12b from de Boor
      for j = 1:jderiv[1]
        kmj = map.order[1] - j
        ilo = kmj
        for jj = 1:kmj
          map.aj[:,jj,1] = (map.aj[:,jj+1,1] - map.aj[:,jj,1]) * kmj /
                           (map.dl[ilo,1] + map.dr[jj,1])
          ilo = ilo - 1
        end
      end
    end

    if jderiv[1] != km1[1]
      # if jderiv /= order - 1, we need to apply the recursive
      # formula from de Boor
      for j = jderiv[1]+1:km1[1]
        kmj = map.order[1] - j
        ilo = kmj
        for jj = 1:kmj
          map.aj[:,jj,1] = (map.aj[:,jj+1,1]*map.dl[ilo,1] + map.aj[:,jj,1]*
                           map.dr[jj,1]) / (map.dl[ilo,1] + map.dr[jj,1])
          ilo = ilo - 1
        end
      end
    end
    # println("map.aj[:,:,1] = \n", map.aj[:,:,2])

    dX[:] = map.aj[:,1,1]

  return nothing
end  # End function calcdXdxi(map, xi, jderiv)


text:
function _bibtorem(bibkey::AbstractString; add_pdfs = false)

    # get ref
    ref = get(bibrefs(), bibkey) do
        error("'", bibkey, "' not found")
    end

    # create rem
    tags = join([remtag(ref.keywords)..., remtag("Document"), remtag("Stub"), remtag("Read Report"), remtag("Read Report")])
    remname = string(bibkey, tags)
        # Meta
            title = _if_nonempty(ref.title, string("\t\t", ref.title, remtag("Title")))
            authors = _if_nonempty(ref.authors, string("\t\t", ref.authors, remtag("Authors")))
            year = _if_nonempty(ref.year, string("\t\t", remref(ref.year), remtag("Issued")))
            doi = _if_nonempty(ref.doi, string("\t\t", ref.doi, remtag("DOI")))
        # Abstract
            abstract = _if_nonempty(ref.abstract, string("\t\t", ref.abstract, remtag("Abstract")))
        # Files
            files = filter(_getfiles(ref)) do path
                endswith(path, ".pdf")
            end
            pdfs = _if_nonempty(files,
                string(
                    "\t\t_pdfs:_",
                    remtag("Edit Later"),
                    "\n\t\t\t",
                    join(files, "\n\t\t\t")
                )
            )
            pdfs = add_pdfs ? pdfs : ""

    return join(
        [
            remname,
                "\tMeta",
                    title,
                    authors,
                    year,
                    doi,
                    pdfs,
                "\tAbstract",
                    abstract
        ]
        , "\n"
    )
end


text:
JSONRPC.parse_params(::Type{Val{Symbol("textDocument/hover")}}, params) = TextDocumentPositionParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/hover")},TextDocumentPositionParams}, server)
    doc = getdocument(server, URI2(r.params.textDocument.uri))
    x = get_expr1(getcst(doc), get_offset(doc, r.params.position))
    documentation = get_hover(x, "", server)
    documentation = get_fcall_position(x, documentation)
    documentation = sanitize_docstring(documentation)

    return Hover(MarkupContent(documentation), missing)
end

function get_hover(x, documentation, server) documentation end

function get_hover(x::EXPR, documentation, server)
    if parentof(x) isa EXPR  && (kindof(x) === CSTParser.Tokens.END || kindof(x) === CSTParser.Tokens.RPAREN || kindof(x) === CSTParser.Tokens.RBRACE || kindof(x) === CSTParser.Tokens.RSQUARE)
        documentation = "Closes `$(typof(parentof(x)))` expression."
    elseif CSTParser.isidentifier(x) && StaticLint.hasref(x)
        if refof(x) isa StaticLint.Binding
            documentation = get_hover(refof(x), documentation, server)
        elseif refof(x) isa SymbolServer.SymStore
            documentation = string(documentation, refof(x).doc)
        end
    end
    return documentation
end

function get_hover(b::StaticLint.Binding, documentation, server)
    if b.val isa EXPR
        if CSTParser.defines_function(b.val)
            while true
                if b.val isa EXPR
                    if parentof(b.val) isa EXPR && typof(parentof(b.val)) === CSTParser.MacroCall && length(parentof(b.val).args) == 3 && typof(parentof(b.val).args[1]) === CSTParser.GlobalRefDoc && CSTParser.isstring(parentof(b.val).args[2])
                        # Binding has preceding docs so use them..
                        documentation = string(documentation, Expr(parentof(b.val).args[2]))
                    elseif CSTParser.defines_function(b.val)
                        documentation = string(documentation, "```julia\n", Expr(CSTParser.get_sig(b.val)), "\n```\n")
                    elseif CSTParser.defines_datatype(b.val)
                        documentation = string(documentation, "```julia\n", Expr(b.val), "\n```\n")
                    end
                elseif b.val isa SymbolServer.SymStore
                    documentation = string(documentation, b.val.doc)
                else
                    break
                end
                if b.prev isa StaticLint.Binding && b.prev != b && (b.prev.type == StaticLint.CoreTypes.Function || b.prev.type == StaticLint.CoreTypes.DataType || b.prev.val isa Union{SymbolServer.FunctionStore,SymbolServer.DataTypeStore})
                    b = b.prev
                else
                    break
                end
            end
        else
            documentation = string(documentation, "```julia\n", Expr(b.val), "\n```\n")
        end
    elseif b.val isa SymbolServer.SymStore
        documentation = string(documentation, b.val.doc)
    elseif b.val isa StaticLint.Binding
        documentation = get_hover(b.val, documentation, server)
    end
    return documentation
end

get_fcall_position(x, documentation) = documentation
function get_fcall_position(x::EXPR, documentation)
    while parentof(x) isa EXPR
        if typof(parentof(x)) === CSTParser.Call
            call_counts = StaticLint.call_nargs(parentof(x))
            call_counts[1] < 5 && return documentation
            arg_i = 0
            for i = 1:length(parentof(x).args)
                arg = parentof(x).args[i]
                if arg == x
                    arg_i = div(i-1, 2)
                end
            end
            arg_i == 0 && return documentation
            fname = CSTParser.get_name(parentof(x))
            if StaticLint.hasref(fname) &&
                (refof(fname) isa StaticLint.Binding && refof(fname).val isa EXPR && CSTParser.defines_struct(refof(fname).val) && StaticLint.struct_nargs(refof(fname).val)[1] == call_counts[1])
                dt_ex = refof(fname).val

                args = CSTParser.defines_mutable(dt_ex) ? dt_ex.args[4] : dt_ex.args[3]
                _fieldname = CSTParser.str_value(CSTParser.get_arg_name(args.args[arg_i]))
                documentation = string("Datatype field `$_fieldname` of $(CSTParser.str_value(CSTParser.get_name(dt_ex)))", "\n", documentation)
            elseif StaticLint.hasref(fname) && (refof(fname) isa SymbolServer.DataTypeStore || refof(fname) isa StaticLint.Binding && refof(fname).val isa SymbolServer.DataTypeStore)
                dts = refof(fname) isa StaticLint.Binding ? refof(fname).val : refof(fname)
                if length(dts.fields) == call_counts[1]
                    documentation = string("Datatype field `$(dts.fields[arg_i])`", "\n", documentation)
                end
            else
                documentation = string("Argument $arg_i of $(call_counts[1]) in call to `", CSTParser.str_value(fname), "`\n", documentation)
            end
            return documentation
        end
        x = parentof(x)
    end
    return documentation
end


text:
using ChaosTools, Test, Statistics, DynamicalSystemsBase

println("\nTesting transit time statistics...")
@testset "Transit time statistics" begin

@testset "Standard map (exact)" begin
# INPUT
ds = Systems.standardmap()
T = 10000 # maximum time

# period 3 of standard map
u0 = SVector(0.8121, 1.6243)
εs = sort!([4.0, 2.0, 0.01]; rev=true)
# radius of 4.0 covers the first 2 points of the period 3 while 2.0 covers only
# the first point. Therefore:
# all first return times must be 1, and all second return times must be 2,
# and 3rd must be same as second (because it doesn't matter how close are)
exits, entries = exit_entry_times(ds, u0, εs, T)
transits, returns = transit_return(exits, entries)

@test all(x -> length(x) > 5, exits)
@test all(x -> length(x) > 5, entries)
@test all(issorted, exits)
@test all(issorted, entries)
@test all(isequal(2), transits[1])
@test all(isequal(1), transits[2])
@test transits[2] == transits[3]
@test all(isequal(1), returns[1])
@test all(isequal(2), returns[2])
@test returns[2] == returns[3]

τ, c = mean_return_times(ds, u0, εs, T)
@test τ == mean.(returns)
@test length(unique(c)) == 1
@test c[1] == T÷3 # should return exactly 1 every 3 steps.

### quasiperiodic around period 3:
u0 = SVector(0.877, 1.565)
εs = sort!([4.0, 0.5, 0.1]; rev=true)
exits, entries = exit_entry_times(ds, u0, εs, T)
transits, returns = transit_return(exits, entries)

@test all(issorted, exits)
@test all(issorted, entries)
@test all(x -> length(x) > 5, exits)
@test all(x -> length(x) > 5, entries)

# For ε=4.0, nothing changes with the before
@test all(isequal(1), returns[1])
@test all(isequal(2), transits[1])

# Similarly, 0.5 should be the same as before
@test all(isequal(1), transits[2])
@test all(isequal(2), returns[2])

# But now, the third entry is different, because it has the size of the quasiperiodic
# stability island torous
@test returns[3] ≠ returns[2]
@test transits[3] ≠ transits[2]
@test any(>(3), returns[3])
@test all(isequal(1), transits[3]) # still need only one step to exit

τ, c = mean_return_times(ds, u0, εs, T)
@test τ == mean.(returns)
@test length(unique(c)) == 2
@test c[1] == T÷3
@test c[2] == T÷3
@test c[3] < T÷3

end

@testset "Towel map (boxes)" begin

to = Systems.towel()
tr = trajectory(to, 5000; Ttr = 10)
u0 = tr[3000]

# With these boxes, in the first 5 steps, the trajectory enters the y and z range
# but not the x range. Therefore it should NOT enter the box. See figure!
εs = [
    SVector(0.05, 0.05, 0.125),
    SVector(0.005, 0.005, 0.025),
]

# Visual guidance
# using PyPlot
# tr0 = trajectory(to, 5, u0)
# fig, axs = subplots(1,3)
# comb = ((1, 2), (1, 3), (2, 3))
# for i in 1:3
#     j, k = comb[i]
#     ax = axs[i]
#     ax.scatter(tr[:, j], tr[:, k], s = 2, color = "C$(i-1)")
#     ax.scatter([u0[j]], [u0[k]], s = 20, color = "k")
#     ax.plot(tr0[:, j], tr0[:, k], color = "k")
#     for l in 1:length(εs)
#         rect = matplotlib.patches.Rectangle(
#         u0[[j, k]] .- εs[l][[j, k]], 2εs[l][j], 2εs[l][k],
#         alpha = 0.25, color = "k"
#         )
#         ax.add_artist(rect)
#     end
# end

exits, entries = exit_entry_times(to, u0, εs, 10000)
transits, returns = transit_return(exits, entries)

@test all(issorted, exits)
@test all(issorted, entries)
@test length(exits[1]) > length(exits[2])
@test returns[1][1] > 5
@test mean(returns[1]) < mean(returns[2])

τ, c = mean_return_times(to, u0, εs, 10000)
@test τ == mean.(returns)
@test c[1] > c[2] > 0

end
#
@testset "Continuous Roessler" begin
# %%
# using OrdinaryDiffEq: Tsit5
alg = DynamicalSystemsBase.DEFAULT_SOLVER

ro = Systems.roessler(ones(3))
u0 = SVector(
    4.705494942754781,
    -10.221120945130545,
    0.06186563933318555
)
εs = sort!([1.0, 0.1, 0.01]; rev=true)
avg_period = 6.0

# # Visual guidance
# using PyPlot
# tr = trajectory(ro, 5000, u0; alg)
# tr0 = trajectory(ro, avg_period, u0; alg)
# tr2 = trajectory(ro, 3avg_period, u0; alg)
# fig, axs = subplots(1,3)
# comb = ((1, 2), (1, 3), (2, 3))
# for i in 1:3
#     j, k = comb[i]
#     ax = axs[i]
#     # ax.plot(tr[:, j], tr[:, k], lw = 2.0, color = "C$(i-1)", alpha = 0.5, marker = "o", ms = 2)
#     ax.plot(tr[:, j], tr[:, k], lw = 2.0, color = "C$(i-1)")
#     ax.scatter([u0[j]], [u0[k]], s = 20, color = "k")
#     ax.plot(tr2[:, j], tr2[:, k], color = "C4", lw = 1.0, ls = "--")
#     ax.plot(tr0[:, j], tr0[:, k], color = "k", lw = 1.0, ls = "--")
#     if eltype(εs[1]) <: Vector
#         for l in 1:length(εs)
#             rect = matplotlib.patches.Rectangle(
#             u0[[j, k]] .- εs[l][[j, k]], 2εs[l][j], 2εs[l][k],
#             alpha = 0.1, color = "k"
#             )
#             ax.add_artist(rect)
#         end
#     else
#         for l in 1:length(εs)
#             circ = matplotlib.patches.Circle(
#                 u0[[j, k]], εs[l]; alpha = 0.1, color = "k"
#             )
#             ax.add_artist(circ)
#         end
#     end
# end

# %%
# We know the average period of the Roessler system. Therefore the return times
# cannot be possibly smaller than it (because u0 is in the xy plane)
# We also know (from the plot) that the first returns after around 3 periods.
# We aaaalso know (by zooming in the plot) that the innermost ball is recurred exactly twice,
# however one of the two crossings is grazing and thus likely to not be spotted by the
# algorithm
τ, c = mean_return_times(ro, u0, εs, 3avg_period; diffeq = (alg = alg,), i=20)
@test c[1] == 1
@test c[2] == c[3] == 0
@test 2avg_period < τ[1] < 3avg_period
@test issorted(c; rev=true)

τ, c = mean_return_times(ro, u0, εs, 5000.0; diffeq = (alg = alg,), i=20)
@test all(τ .> avg_period/2)
@test 0 < c[3] ≤ 2
@test issorted(c; rev=true)

x = sort!(MathConstants.e .^ (-4:0.5:-1); rev = true)
Ts = 10.0 .^ range(3, 6, length = 7)
is = range(10; step = 4, length = 7)
τd, cd_ = mean_return_times(ro, u0, x, Ts; i=is, dmin =10.0)
@test issorted(τd)
@test all(z -> z > 0, cd_)

# figure()
# plot(log.(x), log.(τd); marker ="o")
# d = -ChaosTools.slope(log.(x), log.(τd))
# The slope of the above plot should approximate the fractal dimension
# but I don't find this to be true unfortunately...

# Test with hyper rectangles. Both same outcome because we are anyway in the flat part
εs = [
    SVector(0.1, 0.1, 0.5),
    SVector(0.1, 0.1, 0.05),
]
τ2, c2 = mean_return_times(ro, u0, εs, 5000.0; diffeq = (alg = alg,), i=20)

@test τ2[1] < τ[2]

end
end


text:
using MonteCarloTreeSearch
using Test

@testset "MonteCarloTreeSearch.jl" begin
    # Write your tests here.
end


text:
module LocalBinaryPatterns

using TiledIteration: EdgeIterator
using Interpolations
using ImageCore
using ImageCore.OffsetArrays
using StaticArrays
using IntegralArrays

export
    local_binary_pattern,
    LBP,
    center_mode,
    average_mode

include("local_binary_pattern.jl")

include("block_modes.jl")
include("bitrotate_encoding.jl")
include("uniform_encoding.jl")
include("utils.jl")
include("compat.jl")

end


text:
"""
    Const{T, TSubCon<:Construct{T}, VT} <: Validator{T}

Field enforcing a constant.
"""
struct Const{T, TSubCon<:Construct{T}, VT} <: Validator{T}
    subcon::TSubCon
    value::VT
end

Const(::Type{T}, value::VT) where {T, VT} = Const(Construct(T), value)
Const(value::T) where {T} = Const(Construct(T), value)
Const(value::AbstractArray{V, N}) where {V, N} = Const(SizedArray(V, size(value)...), value)

subcon(wrapper::Const) = wrapper.subcon
function validate(cons::Const{T, TSubCon, VT}, obj::T; contextkw...) where {T, TSubCon, VT}
    (cons.value == obj) === true ? ValidationOK : ValidationError("$obj mismatch the const value $(cons.value).")
end


text:
#using Run a number of regression tests for RunTests.jl
#
# NOTE: the baseline files can be regenerated by passing the --rebaseline
# argument at the command line
#

module StandardizeRunTestsOutput
  import RunTests

  RunTests.set_show_progress(false)
  RunTests.set_tty_cols(80)
  RunTests.set_show_backtraces(false)
end

module RegressionTest
  import RunTests
  import Base.Test.@test

  export regression_test, @test

  const REBASELINE_DFLT = ("--rebaseline" in ARGS)

  function regression_test(fn::Function, baseline_file::AbstractString; rebaseline::Bool=REBASELINE_DFLT)
    result, data, err = RunTests.capture_output() do
      fn()
    end

    if err!=nothing
       rethrow(err)
    end

    if rebaseline
      open(baseline_file, "w") do f
        serialize(f, result)
        serialize(f, data)
      end
      println("Regenerated baseline file: $baseline_file")
    else
      open(baseline_file, "r") do f
        @test result==deserialize(f)
        @test data==deserialize(f)
      end
    end
  end
end

module EachTestFolderOneAtATime
  using RegressionTest
  using RunTests
  import Base.Test.@test

  const PREFIX = VERSION < v"0.5.0-" ? VERSION > v"0.3.0-" ? "0.4_" : "0.2_" : ""

  println("test1")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test1.out")) do
    run_tests("test1")
  end

  println("test2")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test2.out")) do
    run_tests("test2")
  end

  println("test3")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test3.out")) do
    run_tests("test3")
  end

  println("test4")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test4.out")) do
    run_tests("test4")
  end

  println("test5")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test5.out")) do
    run_tests("test5")
  end

  println("test6")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test6.out")) do
    run_tests("test6")
  end

  println("test7")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)test7.out")) do
    run_tests("test7")
  end
end

module AllTheTestsAtOnce
  using RegressionTest
  using RunTests
  import Base.Test.@test

  const PREFIX = VERSION < v"0.5.0-" ? VERSION > v"0.3.0-" ? "0.4_" : "0.2_" : ""

  println("testAll")
  regression_test(joinpath(dirname(@__FILE__), "$(PREFIX)testAll.out")) do
    run_tests()
  end
end


text:
negbinmodel_stan = "data {
  int<lower=1> N;
  int<lower=0> y[N];
}
parameters {
  real<lower=0> alpha;
  real<lower=0> beta;
}
model {
  alpha ~ cauchy(0,10);
  beta ~ cauchy(0,10);
  for (i in 1:N)
    y[i] ~ neg_binomial(alpha, beta);
}
"

@model negbinmodel(y) = begin
  α ~ Cauchy(0,10)
  β ~ Cauchy(0,10)
  for i = 1:length(y)
    y[i] ~ NegativeBinomial(α, β)  # α > 0, 0 < β < 1
  end
end


text:
#
#
using Query
using DataFrames

function make_people_frame( n :: Integer ) ::DataFrame

    pid :: Vector(BigInteger

)


text:
"""
LX_TOKENS

List of names of latex tokens. (See markdown tokens)
"""
const LX_TOKENS = (:LX_BRACE_OPEN, :LX_BRACE_CLOSE, :LX_COMMAND)

"""
$(TYPEDEF)

Structure to keep track of the definition of a latex command declared via a
`\newcommand{\name}[narg]{def}`.
"""
struct LxDef
    name::String
    narg::Int
    def ::AS
    # location of the definition
    from::Int
    to  ::Int
end
# if offset unspecified, start from basically -∞ (configs etc)
function LxDef(name::String, narg::Int, def::AS)
    o = FD_ENV[:OFFSET_LXDEFS] += 5  # precise offset doesn't matter, jus
    LxDef(name, narg, def, o, o + 3) # just forward a bit
end

from(lxd::LxDef) = lxd.from
to(lxd::LxDef)   = lxd.to

"""
pastdef(λ)

Convenience function to return a copy of a definition indicated as having
already been earlier in the context i.e.: earlier than any other definition
appearing in the current chunk.
"""
pastdef(λ::LxDef) = LxDef(λ.name, λ.narg, λ.def)

"""
$(TYPEDEF)

A `LxCom` has a similar content as a `Block`, with the addition of the definition and a vector of brace blocks.
"""
struct LxCom <: AbstractBlock
    ss    ::SubString                 # \\command
    lxdef ::Union{Nothing,Ref{LxDef}} # definition of the command
    braces::Vector{OCBlock}           # relevant {...} with the command
end
LxCom(ss, def)   = LxCom(ss, def, Vector{OCBlock}())
from(lxc::LxCom) = from(lxc.ss)
to(lxc::LxCom)   = to(lxc.ss)

"""
For a given `LxCom`, retrieve the definition attached to the corresponding
`LxDef` via the reference.
"""
function getdef(lxc::LxCom)::Union{Nothing,AS}
    if isnothing(lxc.lxdef)
        return nothing
    end
    return getindex(lxc.lxdef).def
end

"""
For a given `LxCom`, retrieve the name of the command via the reference.
Example: `\\cite` --> `cite`.
"""
function getname(lxc::LxCom)::String
    if isnothing(lxc.lxdef)
        s = String(lxc.ss)
        j = findfirst('{', s)
        return lxc.ss[2:prevind(s, j)]
    end
    return String(getindex(lxc.lxdef).name)[2:end]
end
# XXX if nothing lxc.lxdef then extract from lxc.ss but extract before {}


text:
"""
Supertype for linear minimization oracles.

All LMOs must implement `compute_extreme_point(lmo::LMO, direction)`
and return a vector `v` of the appropriate type.
"""
abstract type LinearMinimizationOracle end

"""
    compute_extreme_point(lmo::LinearMinimizationOracle, direction; kwargs...)

Computes the point `argmin_{v ∈ C} v ⋅ direction`
with `C` the set represented by the LMO.
Most LMOs feature `v` as a keyword argument that allows for an in-place computation whenever `v` is dense.
All LMOs should accept keyword arguments that they can ignore.
"""
function compute_extreme_point end

"""
    CachedLinearMinimizationOracle{LMO}

Oracle wrapping another one of type lmo.
Subtypes of `CachedLinearMinimizationOracle` contain a cache of
previous solutions.

By convention, the inner oracle is named `inner`.
Cached optimizers are expected to implement `Base.empty!` and `Base.length`.
"""
abstract type CachedLinearMinimizationOracle{LMO<:LinearMinimizationOracle} <:
              LinearMinimizationOracle end

# by default do nothing and return the LMO itself
Base.empty!(lmo::CachedLinearMinimizationOracle) = lmo
Base.length(::CachedLinearMinimizationOracle) = 0

"""
    SingleLastCachedLMO{LMO, VT}

Caches only the last result from an LMO and stores it in `last_vertex`.
Vertices of `LMO` have to be of type `VT` if provided.
"""
mutable struct SingleLastCachedLMO{LMO,A} <: CachedLinearMinimizationOracle{LMO}
    last_vertex::Union{Nothing,A}
    inner::LMO
end

# initializes with no cache by default
SingleLastCachedLMO(lmo::LMO) where {LMO<:LinearMinimizationOracle} =
    SingleLastCachedLMO{LMO,AbstractVector}(nothing, lmo)

function compute_extreme_point(
    lmo::SingleLastCachedLMO,
    direction;
    v = nothing,
    threshold=-Inf,
    store_cache=true,
    kwargs...,
)
    if lmo.last_vertex !== nothing && isfinite(threshold)
        if fast_dot(lmo.last_vertex, direction) ≤ threshold # cache is a sufficiently-decreasing direction
            return lmo.last_vertex
        end
    end
    v = compute_extreme_point(lmo.inner, direction, kwargs...)
    if store_cache
        lmo.last_vertex = v
    end
    return v
end

function Base.empty!(lmo::SingleLastCachedLMO)
    lmo.last_vertex = nothing
    return lmo
end

Base.length(lmo::SingleLastCachedLMO) = Int(lmo.last_vertex !== nothing)

"""
    MultiCacheLMO{N, LMO, A}

Cache for a LMO storing up to `N` vertices in the cache, removed in FIFO style.
`oldest_idx` keeps track of the oldest index in the tuple, i.e. to replace next.
`VT`, if provided, must be the type of vertices returned by `LMO`
"""
mutable struct MultiCacheLMO{N,LMO<:LinearMinimizationOracle,A} <:
               CachedLinearMinimizationOracle{LMO}
    vertices::NTuple{N,Union{A,Nothing}}
    inner::LMO
    oldest_idx::Int
end

function MultiCacheLMO{N,LMO,A}(lmo::LMO) where {N,LMO<:LinearMinimizationOracle,A}
    return MultiCacheLMO{N,LMO,A}(ntuple(_ -> nothing, Val{N}()), lmo, 1)
end

function MultiCacheLMO{N}(lmo::LMO) where {N,LMO<:LinearMinimizationOracle}
    return MultiCacheLMO{N,LMO,AbstractVector}(ntuple(_ -> nothing, Val{N}()), lmo, 1)
end

# arbitrary default to 10 points
function MultiCacheLMO(lmo::LMO) where {LMO<:LinearMinimizationOracle}
    return MultiCacheLMO{10}(lmo)
end

# type-unstable
function MultiCacheLMO(n::Integer, lmo::LMO) where {LMO<:LinearMinimizationOracle}
    return MultiCacheLMO{n}(lmo)
end

function Base.empty!(lmo::MultiCacheLMO{N}) where {N}
    lmo.vertices = ntuple(_ -> nothing, Val{N}())
    lmo.oldest_idx = 1
    return lmo
end

Base.length(lmo::MultiCacheLMO) = count(!isnothing, lmo.vertices)

"""
Compute the extreme point with a multi-vertex cache.
`store_cache` indicates whether the newly-computed point should be stored in cache.
`greedy` determines if we should return the first point with dot-product
below `threshold` or look for the best one.
"""
function compute_extreme_point(
    lmo::MultiCacheLMO{N},
    direction;
    v = nothing,
    threshold=-Inf,
    store_cache=true,
    greedy=false,
    kwargs...,
) where {N}
    if isfinite(threshold)
        best_idx = -1
        best_val = Inf
        best_v = nothing
        # create an iteration order to visit most recent vertices first
        iter_order = if lmo.oldest_idx > 1
            Iterators.flatten((lmo.oldest_idx-1:-1:1, N:-1:lmo.oldest_idx))
        else
            N:-1:1
        end
        for idx in iter_order
            if lmo.vertices[idx] !== nothing
                v = lmo.vertices[idx]
                new_val = fast_dot(v, direction)
                if new_val ≤ threshold # cache is a sufficiently-decreasing direction
                    # if greedy, stop and return point
                    if greedy
                        # println("greedy cache sol")
                        return v
                    end
                    # otherwise, keep the index only if better than incumbent
                    if new_val < best_val
                        best_idx = idx
                        best_val = new_val
                        best_v = v
                    end
                end
            end
        end
        if best_idx > 0 # && fast_dot(best_v, direction) ≤ threshold
            # println("cache sol")
            return best_v
        end
    end
    # no interesting point found, computing new
    # println("LP sol")
    v = compute_extreme_point(lmo.inner, direction, kwargs...)
    if store_cache
        tup = Base.setindex(lmo.vertices, v, lmo.oldest_idx)
        lmo.vertices = tup
        # if oldest_idx was last, we get back to 1, otherwise we increment oldest index
        lmo.oldest_idx = lmo.oldest_idx < N ? lmo.oldest_idx + 1 : 1
    end
    return v
end

"""
    VectorCacheLMO{N, LMO, VT}

Cache for a LMO storing an unbounded number of vertices of type `VT` in the cache.
`VT`, if provided, must be the type of vertices returned by `LMO`
"""
mutable struct VectorCacheLMO{LMO<:LinearMinimizationOracle,VT} <:
               CachedLinearMinimizationOracle{LMO}
    vertices::Vector{VT}
    inner::LMO
end

function VectorCacheLMO{LMO,VT}(lmo::LMO) where {VT,LMO<:LinearMinimizationOracle}
    return VectorCacheLMO{LMO,VT}(VT[], lmo)
end

function VectorCacheLMO(lmo::LMO) where {LMO<:LinearMinimizationOracle}
    return VectorCacheLMO{LMO,Vector{Float64}}(AbstractVector[], lmo)
end

function Base.empty!(lmo::VectorCacheLMO)
    empty!(lmo.vertices)
    return lmo
end

Base.length(lmo::VectorCacheLMO) = length(lmo.vertices)

function compute_extreme_point(
    lmo::VectorCacheLMO,
    direction;
    v = nothing,
    threshold=-Inf,
    store_cache=true,
    greedy=false,
    kwargs...,
)
    if isempty(lmo.vertices)
        v = compute_extreme_point(lmo.inner, direction)
        if store_cache
            push!(lmo.vertices, v)
        end
        return v
    end
    best_idx = -1
    best_val = Inf
    best_v = nothing
    for idx in reverse(eachindex(lmo.vertices))
        @inbounds v = lmo.vertices[idx]
        new_val = fast_dot(v, direction)
        if new_val ≤ threshold
            # stop, store and return
            if greedy
                return v
            end
            # otherwise, compare to incumbent
            if new_val < best_val
                best_v = v
                best_val = new_val
                best_idx = idx
            end
        end
    end
    v = best_v
    if best_idx < 0
        v = compute_extreme_point(lmo.inner, direction)
        if store_cache
            # note: we do not check for duplicates. hence you might end up with more vertices,
            # in fact up to number of dual steps many, that might be already in the cache
            # in order to reach this point, if v was already in the cache is must not meet the threshold (otherwise we would have returned it)
            # and it is the best possible, hence we will perform a dual step on the outside.
            #
            # note: another possibility could be to test against that in the if statement but then you might end you recalculating the same vertex a few times.
            # as such this might be a better tradeoff, i.e., to not check the set for duplicates and potentially accept #dualSteps many duplicates.
            push!(lmo.vertices, v)
        end
    end
    return v
end

"""
    ProductLMO(lmos...)

Linear minimization oracle over the Cartesian product of multiple LMOs.
"""
struct ProductLMO{N,TL<:NTuple{N,LinearMinimizationOracle}} <: LinearMinimizationOracle
    lmos::TL
end

function ProductLMO{N}(lmos::TL) where {N,TL<:NTuple{N,LinearMinimizationOracle}}
    return ProductLMO{N,TL}(lmos)
end

function ProductLMO(lmos::Vararg{LinearMinimizationOracle,N}) where {N}
    return ProductLMO{N}(lmos)
end

"""
    compute_extreme_point(lmo::ProductLMO, direction::Tuple; kwargs...)

Extreme point computation on Cartesian product, with a direction `(d1, d2, ...)` given as a tuple of directions.
All keyword arguments are passed to all LMOs.
"""
function compute_extreme_point(lmo::ProductLMO, direction::Tuple; kwargs...)
    return compute_extreme_point.(lmo.lmos, direction; kwargs...)
end

"""
    compute_extreme_point(lmo::ProductLMO, direction::AbstractArray; direction_indices, storage=similar(direction))

Extreme point computation, with a direction array and `direction_indices` provided such that:
`direction[direction_indices[i]]` is passed to the i-th LMO.
The result is stored in the optional `storage` container.

All keyword arguments are passed to all LMOs.
"""
function compute_extreme_point(
    lmo::ProductLMO{N},
    direction::AbstractArray;
    storage=similar(direction),
    direction_indices,
    kwargs...,
) where {N}
    for idx in 1:N
        storage[direction_indices[idx]] .=
            compute_extreme_point(lmo.lmos[idx], direction[direction_indices[idx]]; kwargs...)
    end
    return storage
end


text:
using PyPlot
PyPlot.ioff()
push!(LOAD_PATH, "/home/berceanu/Development/topo-photon/anc/scripts")
import BP
import HH
using PyCall
@pyimport matplotlib.gridspec as gspec

# matplotlib parameters
matplotlib["rcParams"][:update](Dict("axes.labelsize" => 14,
                                     "axes.titlesize" => 20,
                                     "font.size" => 18,
                                     "legend.fontsize" => 14,
                                     "axes.linewidth" => 1.5,
                                     "font.family" => "serif",
                                     "font.serif" => "Computer Modern Roman",
                                     "xtick.labelsize" => 12,
                                     "xtick.major.size" => 5.5,
                                     "xtick.major.width" => 1.5,
                                     "ytick.labelsize" => 12,
                                     "ytick.major.size" => 5.5,
                                     "ytick.major.width" => 1.5,
                                     "text.usetex" => true,
                                     "figure.autolayout" => true))

# system parameters
sN=11 # true system size
# where we pump
n=5
m=0
# inverse lifetime
sγ=0.05
# strength of the trap
sκ = 0.2
q = 7
r = 11 # points in MBZ
# N should be an odd multiple of q
N = r*q # zero-padded system size
l = div(N-1,2)
x = -l:l
δk = 2π/N #resolution in mom space
k = x * δk
#generate all p values inside MBZ
xmbz = -div(r-1,2):div(r-1,2)
kxmbz = xmbz * δk

# exact spectrum, first 15 eigenvalues
exexp = BP.ExactStates(15, :landau, sN, 1/q, sκ)

# calculate HH energies E₀ and E₁
ladder = Array(Float64, (25,25,q))
HH.hhladder!(ladder, 1) # p = 1

#we filter state η
η = 8
β =  η-1
#at energy
sω0=exexp.νs[η]
#with a δ-like pump
δpmp(n₀::Int,m₀::Int) = BP.δpmp(sN; n0=n₀, m0=m₀)
P = δpmp(n,m)

# compute spectrum
sω1=-3.1
sω2=-1.0
dδ=0.001
#spectral range
ν =  sω1:dδ:sω2
sp = BP.Spectrum(collect(ν), P, :landau, 1/q, sγ, sκ)

# no pumping or dissipation
state = BP.getstate(exexp, η)
ψrex = abs2(state)
ψkex = abs2(BP.myfft2(state, k, k)) #77x77 Array{Float64,2} (0.0, 21.3)
ψkmbzex = BP.mbz(ψkex, r, q, kxmbz, k) #77x11 Array{Float64,2} (0., 1.24)
# with pumping
X = BP.getstate(sp, sω0)
ψr = abs2(X)
ψk = abs2(BP.myfft2(X, k, k))
ψkmbz = BP.mbz(ψk, r, q, kxmbz, k)

#analytical w.f. in MBZ
function getχ(Np,q,β)
    α = 1/q
    ydata = linspace(-π,π,Np)
    v = Array(Complex{Float64}, Np)
    for (i,y) in enumerate(ydata)
        v[i] = BP.χ(0.,y,α,β)
    end
    radical = sqrt(sqrt(2/q) / (2π*2^β * factorial(β) * 2π*α))
    radical .* v
end

χ = getχ(N,q,β)

# plotting
ics = 4
t1 = 0.97
b3 = 0.08
uai = 1.2
alfa = 1.5
el = (t1-b3)/(2 + 2alfa + 1/ics + 2/uai)
elp = alfa*el
b1 = t1 - el
t2 = b1 - el/uai
b2 = t2 - 2elp - el/ics
t3 = b2 - el/uai

fig = plt[:figure]()
#plot spectrum
gs1 = gspec.GridSpec(1, 1)
gs1[:update](top=t1, bottom=b1, left=0.225, right=0.8)
ax1 = plt[:subplot](get(gs1, (0, 0)))
ax1[:plot](sp.νs, sp.intensity, "k", linewidth = 2.0)
k₀ = 9 # the first k₀ states are in the n=0 HH band
for (ka, ν) in enumerate(exexp.νs)
    if (ka > k₀) && (mod(ka-k₀, 2) == 1)
        ax1[:axvline](x = ν, color="green", linestyle="-.", lw=2.0)
    else
        ka != 8 && ax1[:axvline](x = ν, color="orange",
           linestyle="dashed", lw=2.0)
        # we excluded the value where we filter
    end
end
ax1[:axvline](x = sω0, color="k", ls=":", lw=2.0)
ax1[:set_xlim](sp.νs[1], sp.νs[end])
ax1[:yaxis][:set_ticks]([0.,1.7,3.5])
ax1[:set_xlabel](L"$\omega_0/J$")
ax1[:set_ylabel](L"$\sum_{m,n} |a_{m,n}|^2$ (arb. units)", fontsize=10)
#plot w.f. in real and mom space
gs2 = gspec.GridSpec(2, 3)
gs2[:update](top=t2, bottom=b2, left=0.225, right=0.8)
ax2 = plt[:subplot](get(gs2, (0, 0)))
ax3 = plt[:subplot](get(gs2, (0, 1)))
ax4 = plt[:subplot](get(gs2, (0, 2)))
ax5 = plt[:subplot](get(gs2, (1, 0)))
ax6 = plt[:subplot](get(gs2, (1, 1)))
ax7 = plt[:subplot](get(gs2, (1, 2)))
axes = [ax2 ax3 ax4;
        ax5 ax6 ax7]
#real space
for (i,ψ) in enumerate((ψrex, ψr))
    ax = axes[i,1]
    img = ax[:imshow](ψ, origin="upper", ColorMap("gist_heat_r"),
                     interpolation="none",
                     extent=[-5.5, 5.5, -5.5, 5.5], aspect=1,
                     vmin=0, vmax=0.1)
    ax[:set_ylabel](L"$n$", labelpad=-5)
    ax[:set_xticks]([-4,0,4])
    ax[:set_yticks]([-4,0,4])
    if i == 2
        ax[:set_xlabel](L"$m$")
    else
        ax[:set_xticklabels]([])
    end
end
#momentum space
for (i,ψ) in enumerate((ψkex, ψk))
    ax = axes[i,2]
    img = ax[:imshow](ψ, origin="upper", ColorMap("gist_heat_r"),
                     interpolation="none",
                     extent=[-π, π, -π, π], aspect=1,
                     vmin=0, vmax=14)
    ax[:set_ylabel](L"$k_y$", labelpad=-8)
    ax[:set_xticks]([-π,0,π])
    if i == 2
        ax[:set_xlabel](L"$k_x$")
        ax[:set_xticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    else
        ax[:set_xticklabels]([])
    end
    ax[:set_yticks]([-π,0,π])
    ax[:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
end
#MBZ
for (i,ψ) in enumerate((ψkmbzex, ψkmbz))
    ax = axes[i,3]
    img = ax[:imshow](ψ, origin="upper", ColorMap("gist_heat_r"),
                      interpolation="none",
                      extent=[-π/q, π/q, -π, π], aspect=1/q,
                      vmin=0, vmax=1)
    # vertical line to indicate slice
    ax[:axvline](x = 0.0, color = "k", linestyle = "-.", linewidth = 2.5)
    ax[:set_ylabel](L"$k_y$", labelpad=-8)
    ax[:set_xticks]([-π/q,0,π/q])
    if i == 2
        ax[:set_xlabel](L"$k_x^0$")
        ax[:set_xticklabels]([L"$-\pi/7$",L"$0$",L"$\pi/7$"])
    else
        ax[:set_xticklabels]([])
    end
    ax[:set_yticks]([-π,0,π])
    ax[:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
end
gs3 = gspec.GridSpec(1, 1)
gs3[:update](top=t3, bottom=b3, left=0.225, right=0.8)
ax8 = plt[:subplot](get(gs3, (0, 0)))
ax8[:plot](k, abs2(χ), "blue", ls="dotted", linewidth = 2.0)
ax8[:plot](k, reverse(ψkmbz[:,6]), "k", linewidth = 2.0)
ax8[:plot](k, reverse(ψkmbzex[:,6]), color="orange", ls="--", linewidth = 2.0)
ax8[:set_xlim](-π, π)
ax8[:set_xticks]([-π,-π/2,0,π/2,π])
ax8[:set_xticklabels]([L"$-\pi$",L"$-\pi/2$",L"$0$",L"$\pi/2$",L"$\pi$"])
ax8[:set_ylabel](L"$|\chi_7(0,k_y)|^2$")
ax8[:set_xlabel](L"$k_y$")
ax8[:yaxis][:set_ticks]([0.,0.2,0.4])
fig[:savefig]("../../figures/exp_fig.pdf", pad_inches=0.0, transparent=true, bbox_inches="tight")
plt[:close](fig)

max_stars_repo_name: berceanu/topo-photon

max_issues_repo_name: berceanu/topo-photon

max_forks_repo_name: berceanu/topo-photon


text:
import Test: @test, @testset
include("constant.jl")

@testset "sample functions" begin
    include("sample.jl")
    Size=64
    idx, prop=newIdx(Size)
    @test prop≈1.0/delIdx(idx, Size)

    β=10.0
    tau, prop=newTau(β)
    @test prop≈1.0/delTau(tau, β)

    Kf=1.0
    K=zero(Mom)
    prop=newK!(K, Kf)
    OldK=copy(K)
    @test prop≈1.0/delK(K, Kf)
    @test K≈OldK #make sure oldK will not be changed

    NewK=zero(Mom)
    #mulitple attempts to test all internal sampling possibilities
    for i in 1:100
        shiftK!(OldK, NewK, Kf)
    end
    @test K≈OldK #make sure shiftK! will not change OldK
end

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
using BinDeps

@BinDeps.setup

libQuasiDefinite = library_dependency("libQuasiDefinite")

provides(SimpleBuild,
         (@build_steps begin
             `make`
             `make install`
         end),
         libQuasiDefinite)

@BinDeps.install [:libQuasiDefinite => :libQuasiDefinite]


text:
using PyPlot2Igor
using Documenter

DocMeta.setdocmeta!(PyPlot2Igor, :DocTestSetup, :(using PyPlot2Igor); recursive=true)

makedocs(;
    modules=[PyPlot2Igor],
    authors="MLackner <lacknersmichael@gmail.com> and contributors",
    repo="https://github.com/MLackner/PyPlot2Igor.jl/blob/{commit}{path}#{line}",
    sitename="PyPlot2Igor.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MLackner.github.io/PyPlot2Igor.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MLackner/PyPlot2Igor.jl",
    devbranch="main",
)


text:
using PyPlot
include("RiemannSolver.jl")

N = 1000
xmin = -0.5
xmax = 0.5
x = collect(LinRange(xmin, xmax, N))
save_figs = true #if this is false the figures will use show(), otherwise they will be saved locally in 'figs' dir

#create directory for figures
if save_figs && !isdir("figs")
	mkdir("figs")
end

struct TestCase
	name::String
	t::Float64
	left::HydroStatus
	right::HydroStatus
end

function calculate_and_plot(test::TestCase)
	name = test.name

	println("Time for the test $name:")

	#dummy run to compile all code
	sample_riemann(x, test.t, test.left, test.right)

	@time profiles = sample_riemann(x, test.t, test.left, test.right)
	density = [status.rho for status in profiles]
	velocity = [status.u for status in profiles]
	pressure = [status.p for status in profiles]
	energy = pressure ./ density ./ (test.left.gamma - 1.)

	figure()
	suptitle(name)
	subplot(221)
	plot(x, density)
	title("density")
	xlabel("x")
	grid(true)
	subplot(222)
	plot(x, velocity)
	title("velocity")
	xlabel("x")
	grid(true)
	subplot(223)
	plot(x, pressure)
	title("pressure")
	xlabel("x")
	grid(true)
	subplot(224)
	plot(x, energy)
	title("energy")
	xlabel("x")
	grid(true)
	tight_layout()

	if !save_figs
		show()
	else

		savefig("figs/$name.pdf",format="pdf",dpi=200)
		savefig("figs/$name.png",format="png",dpi=200)
	end
end

gamma = 7. / 5.

#The first 5 case are taken from Toro's book: (Table 4.1)
#Toro, E.F., Riemann Solvers and Numerical Methods for Fluid Dynamics, A Practical Introduction, 3nd ed., Springer, Berlin, 2009.
#cases 6-13 are takeb from:
#Kamm et. al. 'Enhanced Verification Test Suite for Physics Simulation Codes' SAND2008-7813 2009

#case 1:
sod = TestCase("SodShockTube",  0.25, HydroStatus(1., 0., 1., gamma),
									  HydroStatus(0.125, 0., 0.1, gamma))
calculate_and_plot(sod)

#case 2:
einfeldt = TestCase("Einfeldt123",  0.15, HydroStatus(1., -2., 0.4, gamma),
									      HydroStatus(1., 2., 0.4, gamma))
calculate_and_plot(einfeldt)

#case 3:
toro3 =  TestCase("ToroTest3",  0.012, HydroStatus(1., 0., 1000., gamma),
									   HydroStatus(1., 0., 0.1, gamma))
calculate_and_plot(toro3)

#case 4:
toro4 =  TestCase("ToroTest4",  0.035, HydroStatus(1., 0., 0.01, gamma),
									   HydroStatus(1., 0., 100., gamma))
calculate_and_plot(toro4)

#case 5:
toro5 =  TestCase("ToroTest5",  0.035, HydroStatus(5.99924, 19.5975, 460.894, gamma),
									   HydroStatus(5.99242, -6.19633, 46.0950, gamma))
calculate_and_plot(toro5)

#case 6:
vaccum_left =  TestCase("VaccumExpansionLeft",  0.75, HydroStatus(0., 0., 0., gamma),
													  HydroStatus(1., 0., 1., gamma))
calculate_and_plot(vaccum_left)

#case 7:
vaccum_right =  TestCase("VaccumExpansionRight",  0.75, HydroStatus(1., 0., 1., gamma),
														HydroStatus(0., 0., 0., gamma))
calculate_and_plot(vaccum_right)

#case 8:
rcvcr =  TestCase("RCVCR",  0.75, HydroStatus(1., -4., 0.4, gamma),
								  HydroStatus(1., 4., 0.4, gamma))
calculate_and_plot(rcvcr)

#case 9:
modified_sod = TestCase("ModifiedSod",  0.2, HydroStatus(1., 0.75, 1., gamma),
											 HydroStatus(0.125, 0., 0.1, gamma))
calculate_and_plot(modified_sod)

#case 9:
stream = TestCase("StreamCollision",  0.8, HydroStatus(1., 2., 0.1, gamma),
										   HydroStatus(1., -2., 0.1, gamma))
calculate_and_plot(stream)

#case 10:
leblanc = TestCase("LeBlanc",  0.5, HydroStatus(1., 0., (2. / 3.)*1.e-1, gamma),
									HydroStatus(1.e-3, 0., (2. / 3.)*1.e-10, gamma))
calculate_and_plot(leblanc)

#case 11:
peak = TestCase("PeakProblem",  3.9e-3, HydroStatus(0.1261192, 8.9047029, 782.92899, gamma),
									    HydroStatus(6.591493, 2.2654207, 3.1544874, gamma))
calculate_and_plot(peak)

#case 12:
slow = TestCase("SlowShock",  2., HydroStatus(3.857143, -0.810631, 10.33333, gamma),
								  HydroStatus(1.0, -3.44, 1.0, gamma))
calculate_and_plot(slow)

#case 13:
stationary = TestCase("StationaryContact",  0.012, HydroStatus(1.0, -19.59745, 1.e3, gamma),
												   HydroStatus(1.0, -19.59745, 1.e-2, gamma))
calculate_and_plot(stationary)

##if you want to profile the tests - remove the following comments
##you will also need a larger N than used before
#
#using Profile
#N = 1000000
#x = collect(LinRange(xmin, xmax, N))
#chosen_test = rcvcr
#sample_riemann(x, test.t, test.left, test.right)
#Profile.clear()
#Profile.clear_malloc_data()
#@time sample_riemann(x, test.t, test.left, test.right)
#@profile sample_riemann(x, test.t, test.left, test.right)
#Profile.print()

max_stars_repo_name: ryarazi/ExactRiemannProblemSolver

max_issues_repo_name: ryarazi/ExactRiemannProblemSolver

max_forks_repo_name: ryarazi/ExactRiemannProblemSolver