text:
using DrWatson, Test

d = (a = 0.153456453, b = 5.0, mode = "double")
@test savename(d; digits = 4) == "a=0.1535_b=5_mode=double"
@test savename("n", d) == "n_a=0.153_b=5_mode=double"
@test savename("n/", d) == "n/a=0.153_b=5_mode=double"
@test savename("n/", d; connector = "-") == "n/a=0.153-b=5-mode=double"
@test savename("n\\", d) == "n\\a=0.153_b=5_mode=double"
@test savename(d, "n") == "a=0.153_b=5_mode=double.n"
@test savename("n", d, "n") == "n_a=0.153_b=5_mode=double.n"
@test savename("n", d, "n"; connector = "-") == "n-a=0.153-b=5-mode=double.n"
@test savename(d, allowedtypes = (String,)) == "mode=double"

rick = (never = "gonna", give = "you", up = "!");
@test savename(rick) == "give=you_never=gonna_up=!"

x = 3; y = 5.0;
d = Dict(:x => x, :y => y)
n = (x = x, y = y)

@test d == @dict x y
@test Dict("x" => x, "y" => y) == @strdict x y
@test n == @ntuple x y

z = "lala"
d2 = Dict(:x => x, :y => y, :z => z)
n2 = (x = x, y = y, z= z)

@test d2 == @dict x y z
@test n2 == @ntuple x y z

@test savename(@dict x y) == "x=3_y=5"
@test savename(@ntuple x y) == "x=3_y=5"
w = rand(50)
@test savename(@dict x y w) == savename(@dict x y)
@test savename(@ntuple x y w) == savename(@dict x y)

@test ntuple2dict(@ntuple x y) == @dict x y
@test sort(collect(keys(dict2ntuple(@dict x y)))) == sort(collect(keys(@ntuple x y)))

a = 3; b = 4
c = @ntuple a b
d = 5; e = @dict c d

@test DrWatson.access(e, :c, :a) == a
ff = dict2ntuple(e)
@test DrWatson.access(ff, :c, :a) == a
@test ff.c.a == a

# Expand tests:
a = 3; b = 4
c = @dict a b
d = 5; e = @dict c d

s = savename(e; allowedtypes = (Any,), expand = ["c"])
@test '(' ∈ s
@test ')' ∈ s
@test occursin("a=3", s)
@test occursin("b=4", s)


text:
using Printf

# Functions can have default parameters
function can_i_vote(age=16)
    if age > 18
        println("Can Vote")
    else
        println("Can't Vote")
    end
end

# Parameters are always call-by-value
function try_to_change_value(some_val)
    println("At start of try_to_change_value() some_val is ", some_val)
    some_val = some_val * 2
    println("At end of try_to_change_value() some_val is ", some_val)
end

# Functions can have a variable number of parameters
function get_sum(args...)
    sum = 0
    for a in args
        sum += a
    end
    println(sum)
end

# Functions can return values
function squared(some_int)
    return some_int * some_int
end

# Functions can return tuples
function squared_and_cubed(some_int)
    return (some_int * some_int,
            some_int * some_int * some_int)
end

# Functions can return functions!
function get_multiplier(num)
    return function (x)
        return x * num
    end
end

# Functions can be overloaded
function get_sum(num1::Number, num2::Number)
    return num1 + num2
end

# Functions can be overloaded
function get_sum(num1::String, num2::String)
    return parse(Int32, num1) + parse(Int32, num2)
end

# Functions can be nested
function outer_function()
    function inner_function()
        println("inside inner_function()")
    end

    println("start outer_function()")
    inner_function()
    println("finished outer_function()")
end

# We can create functions..
function doubler(x)
    x * 2
end
# ..and more functions..
function tripler(x)
    x * 3
end
# ..and then pass them as parameters
function increase(n, f)
    f(n)
end

function function_examples()
    println("--------------------------------------")
    println("Function examples")

    println()

    println("We can create single expression functions")
    getSum(x,y) = x + y
    x, y = 1, 2
    @printf("%d + %d = %d\n", x, y, getSum(x, y))

    println("Functions can have default parameters")
    can_i_vote()
    can_i_vote(21)

    println("Arguments are always passed by value")
    some_int = 10
    println("Before calling try_to_change_value() some_int is ", some_int)
    try_to_change_value(some_int)
    println("After calling try_to_change_value() some_int is ", some_int)

    println("Functions can have a variable number of parameters");
    get_sum(1, 2, 3, 4)
    get_sum(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)

    println("Functions can return values")
    println("5 squared is ", squared(5))

    println("Functions can return multiple values")
    println("5 squared and cubed is ", squared_and_cubed(5))

    println("Functions can return multiple values and we can get them individually")
    five_squared, five_cubed = squared_and_cubed(5)
    println("5 squared is ", five_squared)
    println("5 cubed is ", five_cubed)

    println("Functions can return functions")
    times_three = get_multiplier(3)
    println(times_three(10))

    println("Functions can be overloaded")
    println(get_sum(10, 50))
    println(get_sum("10", "50"))

    println("Functions can be nested")
    outer_function()

    println("We can pass functions as parameters");
    println(increase(5, doubler))
    println(increase(5, tripler))
end

ext: jl

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text: function f() end

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:
const IntChar = UInt8
const IntSmallString = UInt64
const maxTagLength = 8
const SmallStringStorage = SVector{maxTagLength,IntChar}
const MSmallStringStorage = MVector{maxTagLength,IntChar} # Mutable SmallString storage

struct SmallString
  data::SmallStringStorage
  length::Int
  SmallString(sv::SmallStringStorage,l::Int) = new(sv,l)
  function SmallString()
    sv = SmallStringStorage(ntuple(_ -> IntChar(0),Val(length(SmallStringStorage))))
    return new(sv, 0)
  end
end

function SmallString(sv::SmallStringStorage)
  sv[1] == IntChar(0) && return SmallString(sv,0)
  len = 1
  while len < length(SmallStringStorage) && @inbounds sv[len+1] ≠ IntChar(0)
    len += 1
  end
  SmallString(sv,len)
end

Base.length(s::SmallString) = s.length
Base.getindex(s::SmallString,n::Integer) = getindex(s.data,n)
function Base.setindex(s::SmallString,val,n::Integer)
  len = length(s)
  if n > len
    len = n
  end
  return SmallString(setindex(s.data,val,n),len)
end

isNull(s::SmallString) = @inbounds s[1] == IntChar(0)

function StaticArrays.push(s::SmallString,val)
  len = length(s)
  return SmallString(setindex(s.data,val,len+1),len+1)
end

# Cast to SmallString:
function SmallString(i::IntSmallString)
  SmallString(unsafe_load(convert(Ptr{SmallStringStorage},pointer_from_objref(MVector{1,IntSmallString}(ntoh(i))))))
end

# Cast to IntSmallString:
function cast_to_uint64(a)
  return ntoh(unsafe_load(convert(Ptr{IntSmallString},pointer_from_objref(MSmallStringStorage(a)))))
end

#isint(i::IntSmallString) = isint(SmallString(i))

function isint(s::SmallString)
  ndigits = 1
  cur_char = Char(s[ndigits])
  !isdigit(cur_char) && return false
  while ndigits < length(s)
    ndigits += 1
    cur_char = Char(s[ndigits])
    !isdigit(cur_char) && return false
  end
  return true
end

#Base.parse(::Type{Int}, i::IntSmallString) = parse(Int,SmallString(i))

function Base.parse(::Type{Int}, s::SmallString)
  n = length(s)
  int = 0
  for j = 1:n
    int = int*10 + parse(Int,Char(s[j]))
  end
  return int
end

# Here we use the StaticArrays comparison
Base.:(==)(s1::SmallString,s2::SmallString) = (s1.data == s2.data)
Base.isless(s1::SmallString,s2::SmallString) = isless(s1.data,s2.data)

########################################################
# Here are alternative SmallString comparison implementations
#

#Base.isless(a::SmallString,b::SmallString) = cast_to_uint64(a) < cast_to_uint64(b)
#Base.:(==)(a::SmallString,b::SmallString) = cast_to_uint64(a) == cast_to_uint64(b)

# Here we use the c-function memcmp (used in Julia string comparison):
#function Base.cmp(a::SmallString, b::SmallString)
#  return ccall(:memcmp, Int32, (Ptr{IntChar}, Ptr{IntChar}, IntChar),
#               a.data, b.data, IntChar(8))
#end
#Base.isless(a::SmallString,b::SmallString) = cmp(a, b) < 0
#Base.:(==)(a::SmallString,b::SmallString) = cmp(a, b) == 0

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

function Base.String(s::SmallString)
  res = ""
  for n=1:length(s)
    res *= Char(s[n])
  end
  return res
end

function Base.show(io::IO, s::SmallString)
  for n=1:length(s)
    print(io,Char(s[n]))
  end
end


text:
# For the conversion to go through ForwardDiff, I need
#Base.convert(FDT::Type{ForwardDiff.Dual{T,V,N}}, p::AIBECS.Parameters) where {T,V,N} = AIBECS.Parameters(convert(Vector{FDT}, vec(p))...)

@testset "Derivatives" begin
    nt = length(T_all)
    n = nt * nb
    @unpack xgeo = p
    x = xgeo * ones(n)
    testp = AIBECS.reconstruct(typeof(p), 2vec(p)) # other value for testing p
    @testset "∇ₓF" begin
        @test ForwardDiff.jacobian(x -> F(x, p), x) ≈ ∇ₓF(x, p) rtol = 1e-14
        @test ForwardDiff.jacobian(x -> F(x, p), 2x) ≈ ∇ₓF(2x, p) rtol = 1e-14
        @test ForwardDiff.jacobian(x -> F(x, p), -x) ≈ ∇ₓF(-x, p) rtol = 1e-14
        @test ForwardDiff.jacobian(x -> F(x, testp), x) ≈ ∇ₓF(x, testp) rtol = 1e-14
    end
    @testset "∇ₚf" begin
        fAD = generate_objective(ωs, μx, σ²x, v, ωp, typeof(p))
        @test fAD(x,p) ≈ f(x,p)
        @test ForwardDiff.jacobian(p -> [fAD(x, p)], p) ≈ ∇ₚf(x, p) rtol = 1e-14
        @test ForwardDiff.jacobian(p -> [fAD(x, p)], testp) ≈ ∇ₚf(x, testp) rtol = 1e-14
        @test ForwardDiff.jacobian(p -> [fAD(2x, p)], p) ≈ ∇ₚf(2x, p) rtol = 1e-14
    end
    @testset "∇ₓf" begin
        @test ForwardDiff.jacobian(x -> [f(x, p)], x) ≈ ∇ₓf(x, p) rtol = 1e-14
        @test ForwardDiff.jacobian(x -> [f(x, testp)], x) ≈ ∇ₓf(x, testp) rtol = 1e-14
        @test ForwardDiff.jacobian(x -> [f(x, p)], 2x) ≈ ∇ₓf(2x, p) rtol = 1e-14
    end
end


text:
module LoadSave # can no longer be called `IO` because it conflicts with `Base.IO`

using ..Utils
using ...LogicCircuits

include("vtree_line_compiler.jl")
include("vtree_loaders.jl")
include("vtree_savers.jl")

include("circuit_line_types.jl")
include("circuit_line_compiler.jl")
include("circuit_loaders.jl")
include("circuit_savers.jl")

include("data_loaders.jl")

include("plot.jl")
end


text:
"""
Algorithm for solving Optimal Stopping Problems.

```julia
NeuralPDE.NNStopping(chain, opt, sdealg, ensemblealg )
```
Arguments:
- `chain`: A Chain neural network with an N-dimensional output according to N stopping times and the last layer (softmax function).
- `opt`: The optimizer to train the neural network. Defaults to `ADAM(0.1)`.
- `sdealg`: The algorithm used to solve the discretized SDE according to the process that X follows. Defaults to `EM()`.
- `ensemblealg`: The algorithm used to solve the Ensemble Problem that performs Ensemble simulations for the SDE. Defaults to `EnsembleThreads()`. See
  the [Ensemble Algorithms](https://diffeq.sciml.ai/stable/features/ensemble/#EnsembleAlgorithms-1)
  documentation for more details.

[1]Becker, Sebastian, et al. "Solving high-dimensional optimal stopping problems using deep learning." arXiv preprint arXiv:1908.01602 (2019).
"""
struct NNStopping{C,O,S,E} <: NeuralPDEAlgorithm
    chain::C
    opt::O
    sdealg::S
    ensemblealg::E
end
NNStopping(chain  ; opt=Flux.ADAM(0.1) , sdealg = EM() , ensemblealg = EnsembleThreads()) = NNStopping(chain , opt , sdealg , ensemblealg)

function DiffEqBase.solve(
    prob::SDEProblem,
    alg::NNStopping;
    abstol = 1f-6,
    verbose = false,
    maxiters = 300,
    trajectories = 1000,
    save_everystep = false,
    dt,
    kwargs...
    )

    tspan = prob.tspan
    sigma = prob.g
    μ = prob.f
    g = prob.kwargs.data.g
    u0 = prob.u0
    ts = tspan[1]:dt:tspan[2]
    N = size(ts)[1]
    T = tspan[2]

    m  = alg.chain
    opt    = alg.opt
    sdealg = alg.sdealg
    ensemblealg = alg.ensemblealg

    prob = SDEProblem(μ,sigma,u0,tspan)
    ensembleprob = EnsembleProblem(prob)
    sim = solve(ensembleprob, sdealg, ensemblealg, dt=dt,trajectories=trajectories,adaptive=false)
    payoff = []
    times = []
    iter = 0
    # for u in sim.u
    un = []
    function Un(n , X )
        if size(un)[1] >= n
            return un[n]
        else
            if(n == 1)
                  ans =  first(m(X[1])[1])
                  un = [ans]
                  return ans
              else
                  ans =  max(first(m(X[n])[n]) , n + 1 - size(ts)[1])*(1 - sum(Un(i , X ) for i in 1:n-1))
                  un = vcat( un , ans)
                  return ans
              end
          end
    end

    function loss()
        reward = 0.00
        for u in sim.u
            X = u.u
            reward = reward + sum(Un(i , X )*g(ts[i] , X[i]) for i in 1 : size(ts)[1])
            un = []
        end
        return 10000 - reward
    end
    dataset = Iterators.repeated(() , maxiters)

    callback = function ()
        l = loss()
        un = []
        println("Current loss is: $l")
    end
    Flux.train!(loss, Flux.params(m), dataset, opt; cb = callback)

    Usum = 0
    ti = 0
    Xt = sim.u[1].u
    for i in 1:N
          un = []
          Usum = Usum + Un(i , Xt)
          if Usum >= 1 - Un(i , Xt)
            ti = i
            break
          end
    end
    for u in sim.u
        X = u.u
        price = g(ts[ti] , X[ti])
        payoff = vcat(payoff , price)
        times = vcat(times, ti)
        iter = iter + 1
        # println("SUM : $sump")
        # println("TIME : $ti")
    end
    sum(payoff)/size(payoff)[1]
 end #solve


text:
using StochasticDiffEq, DiffEqDevTools, ParameterizedFunctions, DiffEqProblemLibrary
using Plots; gr()
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_additive,
            prob_sde_linear, prob_sde_wave
const N = 1000

prob = prob_sde_additive

reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = reltols#[0.0 for i in eachindex(reltols)]
setups = [
          Dict(:alg=>EM(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1))
          Dict(:alg=>RKMil(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRIW1(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRA1(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRA1())
          Dict(:alg=>SRIW1())
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;numruns_error=N,
                      save_everystep = false,
                      parallel_type = :none,
                      error_estimate=:weak_final)#
plot(wp)

sample_size = Int[10;1e2;1e3;1e4]
se = get_sample_errors(prob,setups[6],numruns=sample_size,
                                      sample_error_runs = 100_000,solution_runs=100)

times = [wp[i].times for i in 1:length(wp)]
times = [minimum(minimum(t) for t in times),maximum(maximum(t) for t in times)]
plot!([se[end];se[end]],times,color=:red,linestyle=:dash,label="Sample Error: 10000",lw=3)

prob = prob_sde_additive

reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = reltols#[0.0 for i in eachindex(reltols)]
setups = [
          Dict(:alg=>SRA1())
          Dict(:alg=>SRA2())
          Dict(:alg=>SRA3())
          Dict(:alg=>SOSRA())
          Dict(:alg=>SOSRA2())
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;numruns_error=N,
                      save_everystep = false,
                      maxiters = 1e7,
                      parallel_type = :none,
                      error_estimate=:weak_final)
plot(wp)

sample_size = Int[10;1e2;1e3;1e4]
se = get_sample_errors(prob,setups[4],numruns=sample_size,
                                      sample_error_runs = 100_000,solution_runs=100)

times = [wp[i].times for i in 1:length(wp)]
times = [minimum(minimum(t) for t in times),maximum(maximum(t) for t in times)]
plot!([se[end];se[end]],times,color=:red,linestyle=:dash,label="Sample Error: 10000",lw=3)

prob = prob_sde_linear

reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = reltols#[0.0 for i in eachindex(reltols)]

setups = [Dict(:alg=>SRIW1())
          Dict(:alg=>EM(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1))
          Dict(:alg=>RKMil(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRIW1(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;numruns_error=N,
                      save_everystep = false,
                      maxiters = 1e7,
                      parallel_type = :none,
                      error_estimate=:weak_final)
plot(wp)

sample_size = Int[10;1e2;1e3;1e4]
se = get_sample_errors(prob,setups[1],numruns=sample_size,
                                      sample_error_runs = 100_000,solution_runs=100)

times = [wp[i].times for i in 1:length(wp)]
times = [minimum(minimum(t) for t in times),maximum(maximum(t) for t in times)]
plot!([se[end];se[end]],times,color=:red,linestyle=:dash,label="Sample Error: 10000",lw=3)

prob = prob_sde_linear

reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = reltols#[0.0 for i in eachindex(reltols)]

setups = [Dict(:alg=>EM(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 2))
          Dict(:alg=>RKMil(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 2),:adaptive=>false)
          Dict(:alg=>SRI())
          Dict(:alg=>SRIW1())
          Dict(:alg=>SRIW2())
          Dict(:alg=>SOSRI())
          Dict(:alg=>SOSRI2())
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;numruns_error=N,
                      save_everystep = false,
                      maxiters = 1e7,
                      parallel_type = :none,
                      error_estimate=:weak_final)
plot(wp)

sample_size = Int[10;1e2;1e3;1e4]
se = get_sample_errors(prob,setups[6],numruns=sample_size,
                                      sample_error_runs = 100_000,solution_runs=100)

times = [wp[i].times for i in 1:length(wp)]
times = [minimum(minimum(t) for t in times),maximum(maximum(t) for t in times)]
plot!([se[end];se[end]],times,color=:red,linestyle=:dash,label="Sample Error: 10000",lw=3)

prob = prob_sde_wave

reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = reltols#[0.0 for i in eachindex(reltols)]

setups = [
          Dict(:alg=>EM(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1))
          Dict(:alg=>RKMil(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRIW1(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 1),:adaptive=>false)
          Dict(:alg=>SRIW1())
          ]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;numruns_error=N,
                      save_everystep = false,
                      maxiters = 1e7,
                      parallel_type = :none,
                      error_estimate=:weak_final)
plot(wp)

sample_size = Int[10;1e2;1e3;1e4]
se = get_sample_errors(prob,setups[4],numruns=sample_size,
                                      sample_error_runs = 100_000,solution_runs=100)

times = [wp[i].times for i in 1:length(wp)]
times = [minimum(minimum(t) for t in times),maximum(maximum(t) for t in times)]
plot!([se[end];se[end]],times,color=:red,linestyle=:dash,label="Sample Error: 10000",lw=3)

prob = prob_sde_wave

reltols = 1.0 ./ 10.0 .^ (1:5)
abstols = reltols#[0.0 for i in eachindex(reltols)]

setups = [Dict(:alg=>EM(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 2))
          Dict(:alg=>RKMil(),:dts=>1.0./5.0.^((1:length(reltols)) .+ 2),:adaptive=>false)
          Dict(:alg=>SRI())
          Dict(:alg=>SRIW1())
          Dict(:alg=>SRIW2())
          Dict(:alg=>SOSRI())
          Dict(:alg=>SOSRI2())
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;numruns_error=N,
                      save_everystep = false,
                      maxiters = 1e7,
                      parallel_type = :none,
                      error_estimate=:weak_final)
plot(wp)

sample_size = Int[10;1e2;1e3;1e4]
se = get_sample_errors(prob,setups[6],numruns=sample_size,
                                      sample_error_runs = 100_000,solution_runs=100)

times = [wp[i].times for i in 1:length(wp)]
times = [minimum(minimum(t) for t in times),maximum(maximum(t) for t in times)]
plot!([se[end];se[end]],times,color=:red,linestyle=:dash,label="Sample Error: 10000",lw=3)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
## (c) 2016 David A. van Leeuwen
## tests for ../src/index.jl

import Base.indices

print("getindex, ")
## getindex
## Test Integer indices up to 5 dimensions, as well as CartesianIndexes
for i in 1:5
	dims = fill(3, i)
	n1 = NamedArray(rand(dims...))
	for i in CartesianRange(tuple(dims...))
		@test n1[i.I...] == n1.array[i.I...]
		@test n1[i] == n1.array[i]
		@test n1[[string(ii) for ii in i.I]...] == n1.array[i]
		n1[i.I...] +=  1 ## test setindex
	end
end

include("init-namedarrays.jl")

@test n[:] == n.array[:]
n1 = NamedArray(rand(10))
@test n1[:] == n1
@test [x for x in n] == [x for x in n.array]

## more indexing
first = n.array[1,:]
@test convert(Array, n["one", :]) == first
@test n[Not("two"), :].array == n.array[1:1,:]
@test n[:, ["b", "d"]] == n[:, [2, 4]]

if VERSION < v"0.5.0-dev"
    @test names(n["one", :],1) == ["one"]
end
@test names(n[Not("one"), :],1) == ["two"]
@test names(n[1:1, Not("a")], 2) == ["b", "c", "d"]
## indexing by pair
@test n[:B=>"a", :A=>"two"] == n.array[2, 1]
@test n[:A=>"one", :B=>"d"] == n.array[1, 4]
## https://github.com/nalimilan/FreqTables.jl/issues/10
bi = [false, true, false, true] ## Array{Bool}
for i in 1:2
    @test n[:, bi].array == n.array[:, bi]
    @test names(n[:, bi], 2) == ["b", "d"]
    bi = BitArray(bi)
end


text:
using Parameters
using Statistics

"""
Object for storing useful debugging information about potentially problematic species
tol is the tolerance analyzed at, G is the Gibbs free energy, dy is
the flux to that species and ratio is the ratio between dy and the flux scale
index is the index of the species
"""
@with_kw struct DebugSpecies
    name::String
    G::Float64 = 0.0
    ratio::Float64 = 0.0
    dy::Float64 = 0.0
    tol::Float64 = 0.0
    index::Int64 = 0
end

"""
Object for storing useful debugging information about potentially problematic reactions
tol is the tolerance analyzed at, kf and krev are the forward and reverse rate coefficients
rt is the rate of the reaction, ratio is the ratio between rt and the rate scale
T and P are the temperature and pressure and index is the index of the reaction
"""
@with_kw struct DebugReaction
    reactants::Array{DebugSpecies,1}
    products::Array{DebugSpecies,1}
    rxnstring::String
    tol::Float64 = 0.0
    kf::Float64 = 0.0
    krev::Float64 = 0.0
    rt::Float64 = 0.0
    ratio::Float64 = 0.0
    T::Float64 = 0.0
    P::Float64 = 0.0
    index::Int64 = 0
end

struct DebugMech
    spcs::Array{DebugSpecies,1}
    rxns::Array{DebugReaction,1}
end

function getdebugreaction(rxn::ElementaryReaction;tol=0.0,kf=0.0,krev=0.0,rt=0.0,ratio=0.0,T=0.0,P=0.0,index=0)
    return DebugReaction([getdebugspecies(spc,T) for spc in rxn.reactants],
            [getdebugspecies(spc,T) for spc in rxn.products],
            getrxnstr(rxn),tol,kf,krev,rt,ratio,T,P,index)
end

export getdebugreaction

function getdebugspecies(spc::Species,T::Float64;dy=0.0,ratio=0.0,tol=0.0,index=0)
    G = getGibbs(spc.thermo,T)
    return DebugSpecies(spc.name,G,ratio,dy,tol,index)
end

export getdebugspecies

"""
This function prints a report analyzing rates and dydt to determine reactions
and thermochemistry worth looking at
sim should be a simulation object created from a solution object corresponding to a crash
tol is the ratio relative to the median absolute rate or dn_i/dt value above which reactions
and species will be reported
"""
function analyzecrash(sim::Simulation;tol=1e6)
    rxns = Array{DebugReaction,1}()
    spcs = Array{DebugSpecies,1}()
    t = sim.sol.t[end]
    T = getT(sim,t)
    rts = rates(sim,t)
    rmedian = median(abs.([rt for rt in rts if !isnan(rt) && rt != 0.0]))
    dydt = zeros(length(sim.sol.u[end]))
    dydtreactor!(dydt,sim.sol.u[end],0.0,sim.domain,sim.interfaces;p=sim.p)
    dymedian = median(abs.([dy for dy in dydt if !isnan(dy) && dy != 0.0]))
    y = sim.sol(t)
    p = sim.domain.p
    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(sim.domain,y,t,p)
    for (i,rt) in enumerate(rts)
        if isnan(rt)
            push!(rxns,getdebugreaction(sim.reactions[i];tol=tol,kf=kfs[i],krev=krevs[i],rt=rt,T=T,ratio=NaN,index=i))
        elseif abs(rt/rmedian) > tol
            ratio = abs(rt/rmedian)
            push!(rxns,getdebugreaction(sim.reactions[i];tol=tol,kf=kfs[i],krev=krevs[i],rt=rt,T=T,ratio=ratio,index=i))
        end
    end
    for (i,dy) in enumerate(dydt)
        if i > length(sim.species)
            continue
        end
        if isnan(dy)
            push!(spcs,getdebugspecies(sim.species[i],T;dy=dy,ratio=NaN,tol=tol,index=i))
        elseif abs(dy/dymedian) > tol
            ratio = abs(dy/dymedian)
            push!(spcs,getdebugspecies(sim.species[i],T;dy=dy,ratio=ratio,tol=tol,index=i))
        end
    end
    return DebugMech(spcs,rxns)
end

"""
This function prints a report analyzing rates and dydt to determine reactions
and thermochemistry worth looking at
sim should be a simulation object created from a solution object corresponding to a crash
tol is the ratio relative to the median absolute rate or dn_i/dt value above which reactions
and species will be reported
"""
function analyzecrash(sim::SystemSimulation;tol=1e6)
    rxns = Array{DebugReaction,1}()
    spcs = Array{DebugSpecies,1}()
    t = sim.sol.t[end]
    rts = rates(sim,t)
    rmedian = median(abs.([rt for rt in rts if !isnan(rt) && rt != 0.0]))
    dydt = zeros(length(sim.sol.u[end]))
    dydtreactor!(dydt,sim.sol.u[end],0.0,tuple([s.domain for s in sim.sims]),sim.interfaces;p=sim.p)
    dymedian = median(abs.([dy for dy in dydt if !isnan(dy) && dy != 0.0]))
    y = sim.sol(t)
    p = sim.domain.p
    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR = calcthermo(sim.domain,y,t,p)
    for (i,rt) in enumerate(rts)
        if isnan(rt)
            push!(rxns,getdebugreaction(sim.reactions[i];tol=tol,kf=kfs[i],krev=krevs[i],rt=rt,ratio=NaN,index=i))
        elseif abs(rt/rmedian) > tol
            ratio = abs(rt/rmedian)
            push!(rxns,getdebugreaction(sim.reactions[i];tol=tol,kf=kfs[i],krev=krevs[i],rt=rt,ratio=ratio,index=i))
        end
    end
    for (i,dy) in enumerate(dydt)
        if i > length(sim.species)
            continue
        end
        if isnan(dy)
            push!(spcs,getdebugspecies(sim.species[i];dy=dy,ratio=NaN,tol=tol,index=i))
        elseif abs(dy/dymedian) > tol
            ratio = abs(dy/dymedian)
            push!(spcs,getdebugspecies(sim.species[i];dy=dy,ratio=ratio,tol=tol,index=i))
        end
    end
    return DebugMech(spcs,rxns)
end
export analyzecrash

function getdebugmechstring(mech::DebugMech)
    s = "Crash Analysis Report\n"
    for rxn in mech.rxns
        s *= getdebugrxnstring(rxn)
    end
    for spc in mech.spcs
        s *= getdebugspeciesstring(spc)
    end
    return s
end

export getdebugmechstring

function getdebugrxnstring(rxn)
    s = ""
    rt = rxn.rt
    tol = rxn.tol
    ratio = rxn.ratio
    kf = rxn.kf
    krev = rxn.krev
    index = rxn.index
    if isnan(rt)
        s *= "NaN for Reaction:\n"
    elseif ratio > tol
        s *= "rt/rmedian > $tol, rt=$rt, rt/rmedian=$ratio \nReaction:\n"
    end
    s *= "Index: $index\n"
    s *= rxn.rxnstring * "\n"
    if length(rxn.reactants) == 1
        kunitsf = "s^-1"
    elseif length(rxn.reactants) == 2
        kunitsf = "m^3/(mol*s)"
    elseif length(rxn.reactants) == 3
        kunitsf = "m^6/(mol^2*s)"
    else
        error("cannot accomodate reactants>3")
    end
    if length(rxn.products) == 1
        kunitsr = "s^-1"
    elseif length(rxn.products) == 2
        kunitsr = "m^3/(mol*s)"
    elseif length(rxn.products) == 3
        kunitsr = "m^6/(mol^2*s)"
    else
        error("cannot accomodate products>3")
    end
    s *= "kf: $kf $kunitsf\n"
    s *= "krev: $krev $kunitsr\n"
    s *= "Reactants: \n"
    for reactant in rxn.reactants
        s *= getdebugspeciesstring(reactant)
    end
    s *= "Products: \n"
    for product in rxn.products
        s *= getdebugspeciesstring(product)
    end
    return s
end

export getdebugreactionstring

function getdebugspeciesstring(spc)
    s = ""
    dy = spc.dy
    ratio = spc.ratio
    tol = spc.tol
    index = spc.index
    if isnan(dy)
        s *= "NaN for Species net flux:\n"
    elseif ratio> tol
        s *= "dydt/dydtmedian > $tol, dydt=$dy, dydt/dydtmedian=$ratio \nSpecies:\n"
    end
    if index != 0
        s *= "Index: $index\n"
    end
    G = spc.G/4184.0
    name = spc.name
    s *= "\t$name G(T): $G kcal/mol\n"
    return s
end

export getdebugspeciesstring

function printcrashanalysis(mech::DebugMech)
    println(getdebugmechstring(mech))
end

function printcrashanalysis(sim::Simulation;tol=1e6)
    printcrashanalysis(analyzecrash(sim;tol=tol))
end

export printcrashanalysis

"""
This calculates the collision limit of H + H -> H2 as an upper bound
"""
function calccollisionlimit(T)
    mu = 0.00050397
    sigma = 2.05e-10
    epsilon = 1205.6
    Tr = T*kB*Na/epsilon
    collintegral = 1.16145*Tr^(-0.14874)+0.52487*exp(-0.7732 * Tr)+2.16178*exp(-2.437887*Tr)
    kcoll = sqrt(8.0*pi*kB*T*Na/mu)*sigma^2*collintegral*Na
    return kcoll
end

export calccollisionlimit

"""
Compares bimolecular and trimolecular reactions with the collision limit for H+H->H2
to check whether they are physical, a report with only a title and kcoll means no violators
"""
function analyzecolllimit(phase,Tmin,Tmax,Pmin,Pmax)
    println("Collision Limit Report (Comparisons are done with collision limit for H+H->H2)")
    println("kcoll=$kcoll")
    for T in [Tmin,Tmax]
        for P in [Pmin,Pmax]
            kcoll = calccollisionlimit(T)
            cpdivR,hdivRT,sdivR = calcHSCpdless(phase.vecthermo,T)
            Gs = (hdivRT.-sdivR)*(R*T)
            if phase.diffusionlimited
                diffs = getfield.(phase.species,:diffusion)(T=T,mu=0.0,P=P)
            else
                diffs = Array{Float64,1}()
            end
            kfs,krevs = getkfkrevs(phase,T,P,P/(R*T),1.0,Gs,diffs,V=1.0/C)
            for (i,rxn) in enumerate(phase.reactions)
                boo = false
                if len(rxn.reactants) > 2
                    if kfs[i] > kcoll
                        if !boo
                            println(getrxnstr(rxn))
                            boo = true
                        end
                        kf = kfs[i]
                        ratio = kfs[i]/kcoll
                        println("Violated in forward direction kf=$kf ratio=$ratio T=$T P=$P")
                    end
                end
                if len(rxn.products) > 2
                    if krevs[i] > kcoll
                        if !boo
                            println(getrxnstr(rxn))
                            boo = true
                        end
                        krev = krevs[i]
                        ratio = krevs[i]/kcoll
                        println("Violated in forward direction krev=$krev ratio=$ratio T=$T P=$P")
                    end
                end
            end
        end
    end
end

export analyzecolllimit


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using Test

# Random.seed!(233)
import ADCME: SparseAssembler

sparse_accumulator = load_op_and_grad("./build/libSparseAccumulator","sparse_accumulator", multiple=false)
sparse_accumulator_add = load_op_and_grad("./build/libSparseAccumulator","sparse_accumulator_add", multiple=false)
sparse_accumulator_copy = load_op_and_grad("./build/libSparseAccumulator","sparse_accumulator_copy", multiple=true)

"""
    SparseAssembler(handle::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, tol::Union{PyObject, <:Real}=0.0)

Creates a SparseAssembler for accumulating `row`, `col`, `val` for sparse matrices.
- `handle`: an integer handle for creating a sparse matrix. If the handle already exists, `SparseAssembler` return the existing sparse matrix handle. If you are creating different sparse matrices, the handles should be different.
- `n`: Number of rows of the sparse matrix.
- `tol` (optional): Tolerance. `SparseAssembler` will treats any values less than `tol` as zero.

# Example
```julia
handle = SparseAssembler(100, 5, 1e-8)
op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
J = assemble(5, 5, [op1;op2])
```
`J` will be a [`SparseTensor`](@ref) object.
"""
function SparseAssembler(handle::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, tol::Union{PyObject, <:Real}=0.0)
    n = convert_to_tensor(n, dtype=Int32)
    tol = convert_to_tensor(tol, dtype=Float64)
    handle = convert_to_tensor(handle, dtype=Int32)
    sparse_accumulator(tol, n, handle)
end

"""
    accumulate(handle::PyObject, row::Union{PyObject, <:Integer}, cols::Union{PyObject, Array{<:Integer}}, vals::Union{PyObject, Array{<:Real}})

Accumulates `row`-th row. It adds the value to the sparse matrix
```julia
for k = 1:length(cols)
    A[row, cols[k]] += vals[k]
end
```
`handle` is the handle created by [`SparseAssembler`](@ref).

See [`SparseAssembler`](@ref) for an example.

!!! Note
    `accumulate` returns a `op::PyObject`. Only when `op` is executed, the nonzero values are populated into the sparse matrix.
"""
function accumulate(handle::PyObject, row::Union{PyObject, <:Integer}, cols::Union{PyObject, Array{<:Integer}},
    vals::Union{PyObject, Array{<:Real}})
    row = convert_to_tensor(row, dtype=Int32)
    cols = convert_to_tensor(cols, dtype=Int32)
    vals = convert_to_tensor(vals, dtype=Float64)
    return sparse_accumulator_add(acc, row, cols, vals)
end

"""
    assemble(m::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, ops::PyObject)

Assembles the sparse matrix from the `ops` created by [`accumulate`](@ref). `ops` is either a single output from `accumulate`, or concated from several `ops`
```julia
op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
op = [op1;op2] # equivalent to `vcat([op1, op2]...)`
```
`m` and `n` are rows and columns of the sparse matrix.

See [`SparseAssembler`](@ref) for an example.
"""
function assemble(m::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, ops::PyObject)
    if length(size(ops))==0
        ops = reshape(ops, 1)
    end
    ii, jj, vv = sparse_accumulator_copy(ops)
    return SparseTensor(ii, jj, vv, m, n)
end

handle = SparseAssembler(100, 5, 1e-8)
op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
J = assemble(5, 5, [op1;op2])
# J = assemble(acc, 5, 5)
sess = Session(); init(sess)
run(sess, J)

m = 20
n = 100
handle = SparseAssembler(100, m, 0.0)
op = PyObject[]
A = zeros(m, n)
for i = 1:1
    ncol = rand(1:n, 10)
    row = rand(1:m)
    v = rand(10)
    for (k,val) in enumerate(v)
        @show k
        A[row, ncol[k]] += val
    end
    @show v
    push!(op, accumulate(handle, row, ncol, v))
end
op = vcat(op...)
J = assemble(m, n, op)
B = run(sess, J)
@test norm(A-B)<1e-8

handle = SparseAssembler(100, 5, 1.0)
op1 = accumulate(handle, 1, [1;2;3], [2.0;0.5;0.5])
J = assemble(5, 5, op1)
B = run(sess, J)
@test norm(B-[2.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0])<1e-8

handle = SparseAssembler(100, 5, 0.0)
op1 = accumulate(handle, 1, [1;1], [1.0;1.0])
op2 = accumulate(handle, 1, [1;2], [1.0;1.0])

J = assemble(5, 5, [op1;op2])
B = run(sess, J)
@test norm(B-[3.0  1.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0])<1e-8


text:
# This file is a part of EasyPkg, licensed under the MIT License (MIT).

called_from_repl() = isa(Base.source_path(), Void)

include_base_path() = called_from_repl() ? pwd() : dirname(Base.source_path())

current_src_file() = called_from_repl() ? "" : Base.source_path()

export include_sources

"""
`include_sources(name...)`

Includes all listed Julia sources. Source names are interpreted to be relative
to the path of the calling source file (or the current directory, when run
directly from the REPL).

Names may be ".jl" files or directories. For directories (e.g. "foo"),
`include_sources` will include the Julia source file with the same name inside
that directory (e.g. "foo/foo.jl").
"""

function include_sources(names::AbstractString...)
    base_path = include_base_path()

    for fname in names
        absfname = joinpath(base_path, fname)
        if isdir(absfname)
            subfname = joinpath(absfname, "$(basename(absfname)).jl")
            if isfile(subfname)
                include(subfname)
            else
                error("Error while including directory \"$absfname\", no source file \"$subfname\" found.")
            end
        elseif isfile(absfname) && endswith(absfname, ".jl")
            include(absfname)
        elseif endswith(absfname, ".jl")
            error("No such file: \"$absfname\"")
        else
            error("Don't know how to handle source name \"$absfname\"")
        end
    end
end


text:
using Random, BenchmarkTools, UnsafeArrays, Statistics
using BenchmarkTools: prettytime, prettypercent

push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Shapes

include("benchmarks.jl")

function bgroup!(suite_or_group, name, bench, std)
    g = suite_or_group[name] = BenchmarkGroup()
    g["bench"] = bench
    g["std"] = std
    g
end

# Define a parent BenchmarkGroup to contain our suite
suite = BenchmarkGroup()

suite["sv"] = BenchmarkGroup()
bgroup!(suite["sv"], "vs_assign", sv_vs_assign()...)
bgroup!(suite["sv"], "ms_assign", sv_ms_assign()...)
bgroup!(suite["sv"], "ms_nested_assign", sv_ms_nested_assign()...)

summarize_results(suite_or_group) = summarize_results(suite_or_group, "")
function summarize_results(suite_or_group, name)
    if haskey(suite_or_group, "bench")
        println()
        println("Benchmark $name")
        bench, std = median.(values(suite_or_group))
        if bench.allocs != 0
            @warn "Non-zero allocs in bench: $(bench.allocs)"
        end
        if std.allocs != 0
            @error "Non-zero allocs in std: $(std.allocs)"
        end
        bench_t = time(bench)
        std_t = time(std)
        ratio = (bench_t - std_t) / std_t

        str = "Bench: $(prettytime(bench_t))"
        str *= " Std: $(prettytime(std_t))"
        str *= " Percent slower: $(prettypercent(ratio))"
        println(str)
    else
        for (k, v) in suite_or_group
            summarize_results(v, k)
        end
    end
end

function run_benchmarks()
    tune!(suite)
    results = run(suite)
    summarize_results(results)
    results
end

run_benchmarks()


text:
@testset "Data" begin
  @testset "CurveData" begin
    c = CurveData(Dict(:z => 1:10), [j for i in 1:3, j in 1:10])
    @test coordnames(c) == (:x1, :x2, :x3)
    @test variables(c) == Dict(:z => Int)
    @test npoints(c) == 10

    if visualtests
      gr(size=(800,800))
      c1 = CurveData(Dict(:z => 1:10), [j for i in 1:1, j in 1:10])
      c2 = CurveData(Dict(:z => 1:10), [j for i in 1:2, j in 1:10])
      c3 = CurveData(Dict(:z => 1:10), [j for i in 1:3, j in 1:10])
      @plottest plot(c1,ms=4) joinpath(datadir,"CurveData1D.png") !istravis
      @plottest plot(c2,ms=4) joinpath(datadir,"CurveData2D.png") !istravis
      @plottest plot(c3,ms=4) joinpath(datadir,"CurveData3D.png") !istravis
    end
  end

  @testset "GeoDataFrame" begin
    # basic checks
    data3D   = readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t')
    X, z = valid(data3D, :value)
    @test coordnames(data3D) == (:x, :y, :z)
    @test variables(data3D) == Dict(:value => Float64)
    @test npoints(data3D) == 100
    @test size(X) == (3, 100)
    @test length(z) == 100

    # missing data and NaN
    missdata = readgeotable(joinpath(datadir,"missing.tsv"), delim='\t', coordnames=[:x,:y])
    X, z = valid(missdata, :value)
    @test size(X) == (2,1)
    @test length(z) == 1

    # show methods
    rawdata = DataFrame(x=[1,2,3],y=[4,5,6])
    sdata = GeoDataFrame(rawdata, [:x,:y])
    @test sprint(show, sdata) == "3×2 GeoDataFrame (x and y)"
    @test sprint(show, MIME"text/plain"(), sdata) == "3×2 GeoDataFrame (x and y)\n\n│ Row │ x     │ y     │\n│     │ Int64 │ Int64 │\n├─────┼───────┼───────┤\n│ 1   │ 1     │ 4     │\n│ 2   │ 2     │ 5     │\n│ 3   │ 3     │ 6     │"
    @test sprint(show, MIME"text/html"(), sdata) == "3×2 GeoDataFrame (x and y)\n

	x	y

	Int64	Int64

1	1	4

2	2	5

3	3	6

"

    if visualtests
      gr(size=(800,800))
      df = DataFrame(x=[25.,50.,75.],y=[25.,75.,50.],z=[1.,0.,1.])
      sdata = GeoDataFrame(df, [:x,:y])
      @plottest plot(sdata) joinpath(datadir,"GeoDataFrame.png") !istravis
    end
  end

  @testset "PointSetData" begin
    # basic checks
    data3D = readgeotable(joinpath(datadir,"data3D.tsv"), delim='\t')
    X, z = valid(data3D, :value)
    ps = PointSetData(Dict(:value => z), X)
    X, z = valid(ps, :value)
    @test coordnames(ps) == (:x1, :x2, :x3)
    @test variables(ps) == Dict(:value => Float64)
    @test npoints(ps) == 100
    @test size(X,2) == 100
    @test length(z) == 100

    # show methods
    ps = PointSetData(Dict(:value => [1,2,3]), [1. 0. 1.; 0. 1. 1.])
    @test sprint(show, ps) == "3 PointSetData{Float64,2}"
    @test sprint(show, MIME"text/plain"(), ps) == "3 PointSetData{Float64,2}\n  variables\n    └─value (Int64)"

    if visualtests
      gr(size=(800,800))
      sdata = PointSetData(Dict(:z => [1.,0.,1.]), [25. 50. 75.; 25. 75. 50.])
      @plottest plot(sdata) joinpath(datadir,"PointSetData.png") !istravis
    end
  end

  @testset "RegularGridData" begin
    # basic checks
    g = RegularGridData{Float64}(Dict(:value => rand(100,100)))
    X, z = valid(g, :value)
    @test size(g) == (100, 100)
    @test origin(g) == (0., 0.)
    @test spacing(g) == (1., 1.)
    @test coordnames(g) == (:x1, :x2)
    @test variables(g) == Dict(:value => Float64)
    @test npoints(g) == 10000
    @test size(X) == (2, 10000)
    @test length(z) == 10000

    # show methods
    g = RegularGridData(Dict(:z => [1 2; 3 4]), (0.,0.), (1.,1.))
    @test sprint(show, g) == "2×2 RegularGridData{Float64,2}"
    @test sprint(show, MIME"text/plain"(), g) == "2×2 RegularGridData{Float64,2}\n  origin:  (0.0, 0.0)\n  spacing: (1.0, 1.0)\n  variables\n    └─z (Int64)"

    if visualtests
      gr(size=(800,800))
      sdata = RegularGridData(Dict(:z => [1 2; 3 4]), (0.,0.), (0.1,0.1))
      @plottest plot(sdata) joinpath(datadir,"RegularGridData.png") !istravis
    end
  end

  @testset "StructuredGridData" begin
    X = readdlm(joinpath(datadir,"HurricaneX.dat"))
    Y = readdlm(joinpath(datadir,"HurricaneY.dat"))
    P = readdlm(joinpath(datadir,"HurricaneP.dat"))
    g = StructuredGridData(Dict(:precipitation => P), X, Y)

    # basic checks
    @test size(g) == (221, 366)
    @test coordnames(g) == (:x1, :x2)
    @test variables(g) == Dict(:precipitation => Float64)
    @test npoints(g) == 221*366

    # missing values
    X, z = valid(g, :precipitation)
    @test size(X) == (2, 64416)
    @test length(z) == 64416

    # show methods
    @test sprint(show, g) == "221×366 StructuredGridData{Float64,2}"
    @test sprint(show, MIME"text/plain"(), g) == "221×366 StructuredGridData{Float64,2}\n  variables\n    └─precipitation (Float64)"

    if visualtests
      gr(size=(800,800))
      @plottest plot(g,ms=0.5) joinpath(datadir,"StructuredGridData.png") !istravis
    end
  end
end


text:
module Dance

import REPL

include("Configuration.jl")
include("Logger.jl")
include("utils/URIUtils.jl")
include("Router.jl")
include("utils/PayloadUtils.jl")
include("renderers/HTMLRenderer.jl")
include("renderers/JSONRenderer.jl")
include("renderers/StaticRenderer.jl")
include("renderers/CoreRenderer.jl")
include("engines/CoreEngine.jl")

export start_project

"""
    launch(start_server::Bool)

Launch either (depending on `start_server` param):

- Web server via selected engine
- REPL
"""
function launch(start_server::Bool) :: Union{REPL.REPLBackend, Nothing}
    if start_server
        server_host::String = Configuration.Settings[:server_host]
        server_port::Int64 = Configuration.Settings[:server_port]

        @info "Web server started at $server_host:$server_port"
        CoreEngine.start_server(server_host, server_port)
    else
        # https://github.com/JuliaLang/julia/blob/master/stdlib/REPL/src/REPL.jl
        terminal = REPL.Terminals.TTYTerminal("", stdin, stdout, stderr)
        repl = REPL.LineEditREPL(terminal, true)
        REPL.run_repl(repl)
    end
end

"""
    populate_load_path(file_path::String; ignore_dirs::Array{String, 1})

- Read Settings dict for routes filepath
- Get static dir by parsing routes.jl

- Optional: array of other paths to ignore
"""
function populate_load_path(file_path::String; ignore_dirs::Array{String, 1}) :: Array{String, 1}
    dirs_ignore_array::Array{String, 1} = [".git", ".hg", "node_modules"]
    load_path_array::Array{String, 1} = [abspath(file_path)]

    function _get_static_dir() :: String
        static_dir::String = ""

        for line in readlines(joinpath(file_path, Configuration.Settings[:routes_filename]*".jl"))
            if occursin("static_dir(", line)
                static_dir = replace(
                    strip(
                        split(
                            split(
                                split(line, ",")[2],
                                ")"
                            )[1],
                            "/"
                        )[1]
                    ),
                    "\"" => ""
                )
            end
        end

        return static_dir
    end

    function _populate_load_path() :: Nothing
        for item in readdir(file_path)
            if isdir(item) && !(item in filter(x -> !occursin("/", x), dirs_ignore_array))
                for (root, dirs, files) in walkdir(item)
                    push!(load_path_array, abspath(root))
                    for dir in dirs
                        if size(split(root, "/"))[1]==1
                            for ignore_dir in dirs_ignore_array
                                if !occursin(abspath(ignore_dir), abspath(joinpath(root)))
                                    if !(abspath(root) in load_path_array)
                                        push!(load_path_array, abspath(root))
                                    end
                                end
                            end
                        end

                        ignore_dirs_contains_dir_counter::Int8 = 0
                        for ignore_dir in dirs_ignore_array
                            if occursin(abspath(ignore_dir), abspath(joinpath(root, dir)))
                                ignore_dirs_contains_dir_counter = 1
                                break
                            end
                        end
                        if ignore_dirs_contains_dir_counter==0
                            if !(abspath(joinpath(root, dir)) in load_path_array)
                                push!(load_path_array, abspath(joinpath(root, dir)))
                            end
                        end
                    end
                end
            end
        end
        nothing
    end

    if _get_static_dir()!=""
        push!(dirs_ignore_array, _get_static_dir())
    end
    dirs_ignore_array = vcat(dirs_ignore_array, ignore_dirs)
    _populate_load_path()

    return load_path_array
end

"""
    populate_router(file_path::String)

Populate Router.ROUTES OrderedDict (when including/compiling Julia files)
"""
function populate_router(file_path::String) :: Bool
    return Router.populate(file_path)
end

"""
    populate_settings(file_path::String)

Populate Configuration.Settings dict
"""
function populate_settings(file_path::String) :: Bool
    return Configuration.populate(file_path)
end

"""
    start_project(project_name::String, path::String=".")

- Specify new project name
- Copy `files` dir to root of new project
"""
function start_project(project_name::String, path::String=".") :: Nothing
    project_directory::String = joinpath(abspath(path), project_name)

    mkdir(project_directory)
    cp(joinpath(@__DIR__, "../files"), project_directory; force=true)
    if Sys.iswindows()
        run(`icacls.exe $project_name /reset /T /Q`)
        run(`icacls "$project_name" /reset /T /Q`)
        username::String = read(run(`whoami`), String)
        run(`icacls.exe $project_name /grant $username:F /T /Q`)
        run(`icacls "$project_name" /grant "$username:(OI)(CI)(IO)F" /T /Q`)
        run(`icacls "$project_name" /grant "$username:F" /Q`)

        #= TODO: fix Windows PowerShell escape issue
        See: https://discourse.julialang.org/t/how-to-quote-special-characters-in-run-function/21359/4 =#
        @info "Please run REPL command `run(`attrib –r \"$project_name\\*.*\" /s`)` to finish (due to Julia <-> Windows PowerShell issue)"
    else
        run(`chmod -R 755 $project_directory`)
    end
    @info "Project files created for $project_name"
end

end


text:
function shareParams!(net::Array{Layer})
  for i=1:size(sharedParams,1)
    layerIndices=sharedParams[i,1]
    vars=[sharedParams[i,2];]
    for var in vars
       meanValue=1/length(layerIndices)*reduce(.+,map(layerIndex->net[layerIndex].(var),layerIndices))
       map(layerIndex->net[layerIndex].(var)=meanValue,layerIndices)
    end
  end
end


text:
# classical central limit theorem

using Distributions
using Random
using Gnuplot

function clt(maxn)
    μ = 10
    D = Normal(μ, 1)
    savg = zeros(Float64, maxn)
    for n = 1:maxn
        xbar = mean(rand(D, n))
        savg[n] = sqrt(n)*(xbar - μ)
    end
    savg
end

function plothist(x)
    h = hist(x, bs=0.05)
    bs = fill(h.binsize, length(h.bins));
    println(h.counts)
    @gp h.bins h.counts./2 bs./2 h.counts./2 "w boxxy notit fs solid 0.4"
end

max_stars_repo_name: affans/notebooks

max_issues_repo_name: affans/notebooks

max_forks_repo_name: affans/notebooks


text:
export linear,∂linear

function linear(rᵢ)
  rₒ=rᵢ
end
function ∂linear(rᵢ)
  rₒ=ones(rᵢ)
end


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

mutable struct GatewayRoute <: SwaggerModel
    localAddress::Any # spec type: Union{ Nothing, String } # spec name: localAddress
    network::Any # spec type: Union{ Nothing, String } # spec name: network
    nextHop::Any # spec type: Union{ Nothing, String } # spec name: nextHop
    sourcePeer::Any # spec type: Union{ Nothing, String } # spec name: sourcePeer
    origin::Any # spec type: Union{ Nothing, String } # spec name: origin
    asPath::Any # spec type: Union{ Nothing, String } # spec name: asPath
    weight::Any # spec type: Union{ Nothing, Int32 } # spec name: weight

    function GatewayRoute(;localAddress=nothing, network=nothing, nextHop=nothing, sourcePeer=nothing, origin=nothing, asPath=nothing, weight=nothing)
        o = new()
        validate_property(GatewayRoute, Symbol("localAddress"), localAddress)
        setfield!(o, Symbol("localAddress"), localAddress)
        validate_property(GatewayRoute, Symbol("network"), network)
        setfield!(o, Symbol("network"), network)
        validate_property(GatewayRoute, Symbol("nextHop"), nextHop)
        setfield!(o, Symbol("nextHop"), nextHop)
        validate_property(GatewayRoute, Symbol("sourcePeer"), sourcePeer)
        setfield!(o, Symbol("sourcePeer"), sourcePeer)
        validate_property(GatewayRoute, Symbol("origin"), origin)
        setfield!(o, Symbol("origin"), origin)
        validate_property(GatewayRoute, Symbol("asPath"), asPath)
        setfield!(o, Symbol("asPath"), asPath)
        validate_property(GatewayRoute, Symbol("weight"), weight)
        setfield!(o, Symbol("weight"), weight)
        o
    end
end # type GatewayRoute

const _property_map_GatewayRoute = Dict{Symbol,Symbol}(Symbol("localAddress")=>Symbol("localAddress"), Symbol("network")=>Symbol("network"), Symbol("nextHop")=>Symbol("nextHop"), Symbol("sourcePeer")=>Symbol("sourcePeer"), Symbol("origin")=>Symbol("origin"), Symbol("asPath")=>Symbol("asPath"), Symbol("weight")=>Symbol("weight"))
const _property_types_GatewayRoute = Dict{Symbol,String}(Symbol("localAddress")=>"String", Symbol("network")=>"String", Symbol("nextHop")=>"String", Symbol("sourcePeer")=>"String", Symbol("origin")=>"String", Symbol("asPath")=>"String", Symbol("weight")=>"Int32")
Base.propertynames(::Type{ GatewayRoute }) = collect(keys(_property_map_GatewayRoute))
Swagger.property_type(::Type{ GatewayRoute }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_GatewayRoute[name]))}
Swagger.field_name(::Type{ GatewayRoute }, property_name::Symbol) =  _property_map_GatewayRoute[property_name]

function check_required(o::GatewayRoute)
    true
end

function validate_property(::Type{ GatewayRoute }, name::Symbol, val)
end


text:
@testset "AbstractArray" begin
    unsorted_dv = @data [2, 1, missing]

    # TODO: Make this work
    # tiedrank(dv)

    @test first(unsorted_dv) == 2
    @test ismissing(last(unsorted_dv))

    # ismissing with AbstractArray
    a = [1, 2, 3]
    @test ismissing.(a) == falses(3)
    a = Any[1, 2, missing, 3]
    @test ismissing.(a) == [false, false, true, false]
    for i = 1:length(a)
        @test ismissing(a, i) == ismissing.(a)[i]
    end
end


text:
# a basic create robot type node example

using Caesar, IncrementalInference
using CloudGraphs

using KernelDensityEstimate

# switch IncrementalInference to use CloudGraphs (Neo4j) data layer
# dbaddress = "localhost"
# dbusr = ""
# dbpwd = ""
# mongoaddress = "localhost"
session = "SESSTEST"
robot = "TESTROBOT"
user = "TESTUSER"
include(joinpath(dirname(@__FILE__) ,"blandauth.jl"))
Nparticles = 200
println("Attempting to solve session $(session)...")

configuration = CloudGraphs.CloudGraphConfiguration(dbaddress, 7474, dbusr, dbpwd, mongoaddress, 27017, false, "", "");
cloudGraph = connect(configuration, IncrementalInference.encodePackedType, Caesar.getpackedtype, IncrementalInference.decodePackedType);
# register types of interest in CloudGraphs
# registerGeneralVariableTypes!(cloudGraph)
Caesar.usecloudgraphsdatalayer!()

# Uncomment out for command line operation
# cloudGraph, addrdict = standardcloudgraphsetup(drawdepth=true)
# session = addrdict["session"]
# Nparticles = addrdict["num particles"]

# this is being replaced by cloudGraph, added here for development period
fg = initfg(sessionname=session, robotname=robot, username=user, cloudgraph=cloudGraph)

# Robot navigation and inference type stuff
N=Nparticles
doors = reshape([-100.0;0.0;100.0;300.0],1,4)
cov = 3.0*ones(1,1)

#TODO: WIP
# robot style, add first pose vertex
v1 = addVariable!(fg,:x1,ContinuousScalar,N=N, labels=["POSE"])

# add a prior for the initial position of the robot
fct = Obsv2(doors, cov, [1.0])
# TODO: Fix here - what's up with it failing to encode?
#WIP
f1  = addFactor!(fg,[v1], fct)

# add second pose vertex
tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addVariable!(fg, :x2, ContinuousScalar, N=N,labels=["POSE"])

# now add the odometry factor between them
f1 = addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

v3=addVariable!(fg,:x3,4.0*randn(1,N)+getVal(v2)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors, cov', [1.0]))

v4=addVariable!(fg,:x4,2.0*randn(1,N)+getVal(v3)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

# if true
l1=addVariable!(fg, :l1, 0.5*randn(1,N)+getVal(v3)+64.0, N=N,labels=["LANDMARK"])
addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))
# end

v5=addVariable!(fg,:x5,2.0*randn(1,N)+getVal(v4)+50.0, N=N,labels=["POSE"])
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

v6=addVariable!(fg,:x6,1.25*randn(1,N)+getVal(v5)+40.0, N=N,labels=["POSE"])
addFactor!(fg,[v5;v6],Odo([40.0]',[1.25]',[1.0]))

v7=addVariable!(fg,:x7,2.0*randn(1,N)+getVal(v6) +60.0, N=N,labels=["POSE"])
addFactor!(fg,[v6;v7],Odo([60.0]',[2.0]',[1.0]))

f3 = addFactor!(fg,[v7], Obsv2(doors, cov', [1.0]))

# Now operate with the data in the DB

tree = wipeBuildNewTree!(fg,drawpdf=true)
# run(`evince bt.pdf`)

# gc()
# recursive solving (single process, easy stack trace for debugging)
println("Starting inference over tree")
inferOverTreeR!(fg, tree)

# Get neighbors
# neighs = CloudGraphs.get_neighbors(fg.cg, cv1r)

gt = Dict{Symbol, Array{Float64,2}}()
gt[:x1]=([0.0;1.97304 ]')' # -0.0342366
gt[:x2]=([50.0; 2.83153 ]')' # 49.8797
gt[:x3]=([100.0; 1.65557 ]')' # 99.8351
gt[:x4]=([150.0; 1.64945 ]')' # 148.637
gt[:x5]=([200.0; 1.77992 ]')' # 198.62
gt[:x6]=([240.0; 2.20466 ]')' # 238.492
gt[:x7]=([300.0; 2.14353 ]')' # 298.467
gt[:l1]=([165.0; 1.17284 ]')' # 164.102

# get vertex back from DB
x1neoID = fg.cgIDs[fg.IDs[:x1]]
cv1r = CloudGraphs.get_vertex(fg.cg, x1neoID, false)
v1r = cloudVertex2ExVertex(cv1r)
plotKDE(getKDE(v1r))

plotKDE( getVertKDE(fg, :x2) )

# draw all beliefs
using Gadfly

DOYTICKS = false
xx,ll = ls(fg)
msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
evalstr = ""
for i in 1:length(msgPlots)
    evalstr = string(evalstr, ",msgPlots[$(i)]")
end
pl = eval(parse(string("vstack(",evalstr[2:end],")")));
Gadfly.draw(PDF("cg4doors.pdf",15cm,20cm),pl) # can also do PNG


text:
module Utilities

export ComponentValue, ComponentType, Component, ComponentChild, Enumerated,
        FlagList, RangeValue, Size, Count, Normalized, Color, FontSize, Font,
        Side, AxisRef, XAnchor, YAnchor, Margin, Title, ColorScale,
        ColorScaleSet, Orientation, Easing, Ordering, Transition, NameLength,
        HoverLabel, Calendar, TickMode, Ticks, TickFormatStops, ShowTickPrefix,
        ShowTickSuffix, ExponentFormat, ShowExponent, Dash, default

abstract type ComponentValue end
abstract type ComponentType end
abstract type Component <: ComponentType end
abstract type ComponentChild <: ComponentType end

const default = Dict{String,Component}()

diff(ref::ComponentValue,obj::ComponentValue) = diff(ref.value,obj.value)

function diff(ref,obj)
    ref != obj && return obj
    return
end

struct Enumerated{E,T} <: ComponentValue
    value::Union{String,Bool}

    function Enumerated{E,T}(v) where {E,T}
        Symbol(v) in T || v in T ||
        error("$(E): \"$(v)\" not found in $(tuple([isbits(t) ? t : string(t) for t in T]...))")
        return new{E,T}(v)
    end
end

struct FlagList{S,T} <: ComponentValue
    value::Union{String,Bool}

    function FlagList{S,T}(s::Union{String,Bool}) where {S,T}
        flags = isa(s,String) ? Symbol.(split(s,'+')) : Bool[s]
        length(flags) == 1 && (Symbol(flags[1]) in T || flags[1] in T) && return new(s)
        for flag in flags
            Symbol(flag) in S || flag in S || error("$(flag) not found in $(S)")
        end
        return new(s)
    end
end

struct RangeValue{N,R,T} <: ComponentValue
    value::N

    function RangeValue{N,R,T}(v::N) where {N,R,T}
        errmsg = "$(v) is outside the range $(T[1] ? "[" : "(")$(R[1]),$(R[2])$(T[2] ? "]" : ")")"
        (v < R[1]) && error(errmsg)
        (v > R[2]) && error(errmsg)
        (!T[1] && v == R[1]) && error(errmsg)
        (!T[2] && v == R[2]) && error(errmsg)
        return new(v)
    end
end
RangeValue(range,value) = RangeValue{typeof(value),range,(true,true)}(value)

const Size{S} = RangeValue{S,(0,Inf),(true,false)}

const Count = RangeValue{Int,(1,Inf),(true,false)}

const Normalized = RangeValue{Float64,(0,1),(true,true)}

struct Color <: ComponentValue
    value::String
end

const FontSize = RangeValue{Int,(1,Inf),(true,false)}

struct Font <: ComponentType
    family::Union{String,Vector{String}}
    size::Union{FontSize,Vector{FontSize}}
    color::Union{Color,Vector{Color}}

    function Font(;
        family::Union{String,Vector{String}}="",
        size::Union{Int,Vector{Int}}=6,
        color::Union{String,Vector{String}}="#000")
        return new(family,FontSize.(size),Color.(color))
    end
end

const Side = Enumerated{:Side,(:left,:right,:top,:bottom)}

const AxisRef = Enumerated{:AxisRef,(:container,:paper)}

const XAnchor = Enumerated{:XAnchor,(:auto,:left,:center,:right)}

const YAnchor = Enumerated{:YAnchor,(:auto,:top,:middle,:bottom)}

struct Margin <: ComponentType
    t::Size{Int}
    r::Size{Int}
    b::Size{Int}
    l::Size{Int}
    pad::Size{Int}
    autoexpand::Bool

    function Margin(;
        t::Int = 0,
        r::Int = 0,
        b::Int = 0,
        l::Int = 0,
        pad::Int = 0,
        autoexpand::Bool = true)
        return new(
            Size{Int}(t),
            Size{Int}(r),
            Size{Int}(b),
            Size{Int}(l),
            Size{Int}(pad),
            autoexpand)
    end
end

struct Title <: ComponentType
    text::String
    font::Font
    side::Side
    xref::AxisRef
    yref::AxisRef
    x::Normalized
    y::Normalized
    xanchor::XAnchor
    yanchor::YAnchor
    pad::Margin

    function Title(;
        text::String = "",
        font::Font = Font(),
        side::String = "top",
        xref::String = "container",
        yref::String = "container",
        x = .5,
        y = .5,
        xanchor::String = "auto",
        yanchor::String = "auto",
        pad = Margin())
        return new(
            text,
            font,
            Side(side),
            AxisRef(xref),
            AxisRef(yref),
            Normalized(Float64(x)),
            Normalized(Float64(y)),
            XAnchor(xanchor),
            YAnchor(yanchor),
            pad)
    end
end
Title(t::Title) = t

struct ColorScale <: ComponentType
    colors::Vector{Tuple{Normalized,Color}}
end

colorscalesets = (:Greys,:YlGnBu,:Greens,:YlOrRd,:Bluered,:RdBu,:Reds,:Blues,:Picnic,:Rainbow,:Portland,:Jet,:Hot,:Blackbody,:Earth,:Electric,:Viridis,:Cividis)
const ColorScaleSet = Enumerated{:colorscaleset,colorscalesets}

const Orientation = Enumerated{:Orientation,(:v,:h)}

easings = Symbol.(("linear","quad","cubic","sin","exp","circle","elastic","back","bounce","linear-in","quad-in","cubic-in","sin-in","exp-in","circle-in","elastic-in","back-in","bounce-in","linear-out","quad-out","cubic-out","sin-out","exp-out","circle-out","elastic-out","back-out","bounce-out","linear-in-out","quad-in-out","cubic-in-out","sin-in-out","exp-in-out","circle-in-out","elastic-in-out","back-in-out","bounce-in-out"))
const Easing = Enumerated{:Easing,easings}

const Ordering = Enumerated{:Ordering,Symbol.(("layout first","traces first"))}

struct Transition
    duration::Size{Int}
    easing::Easing
    ordering::Ordering

    function Transition(;
        duration::Int = 500,
        easing::String = "cubic-in-out",
        ordering::String = "layout first")
        return new(
            Size{Int}(duration),
            Easing(easing),
            Ordering(ordering)
        )
    end
end

const NameLength = RangeValue{Int,(-1,Inf),(true,false)}

struct HoverLabel <: ComponentType
    bgcolor::Union{Color,Vector{Color}}
    bordercolor::Union{Color,Vector{Color}}
    font::Font
    namelength::Union{NameLength,Vector{NameLength}}

    function HoverLabel(;
        bgcolor::Union{String,Vector{String}} = "#fff",
        bordercolor::Union{String,Vector{String}} = "#fff",
        font = Font(),
        namelength::Union{Int,Vector{Int}} = -1)
        return new(
            Color.(bgcolor),
            Color.(bordercolor),
            font,
            NameLength.(namelength))
    end
end

const Calendar = Enumerated{:Calendar,(:gregorian,:chinese,:coptic,:discworld,:ethiopian,:hebrew,:islamic,:mayan,:nanakshahi,:nepali,:persian,:jalali,:taiwan,:thai,:ummalqura)}

const TickMode = Enumerated{:TickMode,(:auto,:linear,:array)}

const TickLocation = Enumerated{:TickLocation,(:outside,:inside,Symbol(""))}

const TicksOn = Enumerated{:TicksOn,(:labels,:boundaries)}

const Mirror = Enumerated{:Mirror,(true,:ticks,false,:all,:allticks)}

struct TickFormatStops <: ComponentType
    enable::Bool
    dtickrange::Vector{NamedTuple{(:min,:max),NTuple{2,Float64}}}
    value::String

    function TickFormatStops(;
        enable::Bool = true,
        dtickrange = Vector{NamedTuple{(:min,:max),NTuple{2,Float64}}}(),
        value::String = "")
        return new(enable,dtickrange,value)
    end
end

const ShowTickPrefix = Enumerated{:ShowTickPrefix,(:all,:first,:last,:none)}

const ShowTickSuffix = Enumerated{:ShowTickSuffix,(:all,:first,:last,:none)}

const ExponentFormat = Enumerated{:ExponentFormat,(:none,:e,:E,:power,:SI,:B)}

const ShowExponent = Enumerated{:ShowExponent,(:all,:first,:last,:none)}

struct Ticks <: ComponentChild
    tickmode::TickMode
    nticks::Size{Int}
    tick0::Any
    dtick::Any
    tickvals::Vector{Float64}
    ticktext::Vector{String}
    ticks::TickLocation
    tickson::TicksOn
    mirror::Mirror
    ticklen::Size{Int}
    tickwidth::Size{Int}
    tickcolor::Color
    showticklabels::Bool
    automargin::Bool
    tickfont::Font
    tickangle::Float64
    tickprefix::String
    showtickprefix::ShowTickPrefix
    ticksuffix::String
    showticksuffix::ShowTickSuffix
    showexponent::ShowExponent
    exponentformat::ExponentFormat
    separatethousands::Bool
    tickformat::String
    tickformatstops::TickFormatStops

    function Ticks(;
        tickmode::String = "auto",
        nticks::Int = 0,
        tick0::Any = "",
        dtick::Any = "",
        tickvals::Vector{Float64} = Float64[],
        ticktext::Vector{String} = String[],
        ticks::String = "inside",
        tickson::String = "labels",
        mirror::Union{String,Bool} = false,
        ticklen::Int = 5,
        tickwidth::Int = 1,
        tickcolor::String = "#444",
        showticklabels::Bool = true,
        automargin::Bool = false,
        tickfont::Font = Font(),
        tickangle::Float64 = 0.,
        tickprefix::String = "",
        showtickprefix::String = "all",
        ticksuffix::String = "",
        showticksuffix::String = "all",
        showexponent::String = "all",
        exponentformat::String = "B",
        separatethousands::Bool = false,
        tickformat::String = "",
        tickformatstops::TickFormatStops = TickFormatStops())
        return new(
            TickMode(tickmode),
            Size{Int}(nticks),
            tick0,
            dtick,
            tickvals,
            ticktext,
            TickLocation(ticks),
            TicksOn(tickson),
            Mirror(mirror),
            Size{Int}(ticklen),
            Size{Int}(tickwidth),
            Color(tickcolor),
            showticklabels,
            automargin,
            tickfont,
            tickangle,
            tickprefix,
            ShowTickPrefix(showtickprefix),
            ticksuffix,
            ShowTickSuffix(showticksuffix),
            ShowExponent(showexponent),
            ExponentFormat(exponentformat),
            separatethousands,
            tickformat,
            tickformatstops)
    end
end

const Dash = Enumerated{:Dash,(:solid,:dot,:dash,:longdash,:dashdot,:longdashdot)}

end


text:
function rna(dna)
end

max_stars_repo_name: benallan/lovelace-website

max_issues_repo_name: benallan/lovelace-website

max_forks_repo_name: benallan/lovelace-website


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using DynamicalSystems, InteractiveDynamics
import GLMakie
using AbstractPlotting
using LinearAlgebra
using Statistics

# Code for 3D animation
ds = Systems.lorenz()
u0 = trajectory(ds, 10; Ttr = 100)[end]
dt = 0.005
N = 100
framerate = 30
systemtitle = "Lorenz63"
azimuth = 2.24
elevation = 0.38

ds = Systems.towel()
u0 = trajectory(ds, 10; Ttr = 100)[end]
dt = 1
N = 15 # total steps to take
framerate = 2
systemtitle = "towelmap"

elevation = 0.5
azimuth = 6.02

# %%
tinteg = tangent_integrator(ds, 3; u0)
Y = get_deviations(tinteg)

# Now I define functions that make an ellipsoid given the point vectors
function ellipsoid_transform(Y)
    svddecomp = svd(Y)
    d = svddecomp.S
    s = d ./ maximum(d)
    Q = svddecomp.U * svddecomp.Vt # rotation matrix
    θx = atan(Q[3,2],Q[3,3])
    θy = atan(-Q[3,1], sqrt(Q[3,2]^2+Q[3,3]^2))
    θz = atan(Q[2,1], Q[1,1])
    M = AbstractPlotting.rotationmatrix_x(θx) *
        AbstractPlotting.rotationmatrix_y(θy) *
        AbstractPlotting.rotationmatrix_z(θz) *
        AbstractPlotting.scalematrix(Vec3f0(s...))
    return M, s
end
function transform_mesh(msh, mat4x4)
    pos_trans = Point3f0.(Ref(mat4x4) .* to_ndim.(Point4f0, spheremesh.position, 0))
    AbstractPlotting.GeometryBasics.Mesh(pos_trans, AbstractPlotting.GeometryBasics.faces(msh))
end
function ellipsoid_from_pointcloud(v, msh = spheremesh)
    M, s = ellipsoid_transform(v)
    return transform_mesh(msh, M), s
end
spheremesh = AbstractPlotting.GeometryBasics.mesh(Sphere(Point3f0(0), 1))

ellipsoid, s = ellipsoid_from_pointcloud(Y, spheremesh)
colors = [RGBAf0((abs.(u)./norm(u))..., 1) for u in spheremesh.position]

ellipsobs = Observable(ellipsoid)

# Setup figure and plot everything

# fig = Figure(resolution = (1200, 600))
fig = Figure(resolution = (1200, 600)); display(fig)
# Plot trajectory
t = Observable(0.0)
tstr = lift(x -> systemtitle *", t = $(round(x; digits = 5))", t)
axtr = Axis3(fig[1,1]; title = tstr, titlealign = :left)
tr = trajectory(ds, 100N*dt, u0; dt)
dotobs = Observable([u0])
trobs = Observable([u0])
if ds isa ContinuousDynamicalSystem
    lines!(axtr, tr.data; linewidth = 0.2, color = RGBAf0(0,0,0,0.5), transparency = true )
    lines!(axtr, trobs; linewidth = 4.0, color = to_color(JULIADYNAMICS_COLORS[1]), overdraw = true)
else
    scatter!(axtr, tr.data; markersize = 1000, color = RGBAf0(0,0,0,0.5), strokewidth=0, transparency = true )
    scatter!(axtr, trobs; markersize = 2000, color = to_color(JULIADYNAMICS_COLORS[1]), overdraw = true)
end
scatter!(axtr, dotobs; markersize = 4000, marker=:diamond,
color = to_color(JULIADYNAMICS_COLORS[3]), overdraw = true)
axtr.azimuth = azimuth
axtr.elevation = elevation

# Plot ellipsoid
ax3D = Axis3(fig[1, 2]; aspect = :data)
mesh!(ellipsobs, color = colors)

ax3D.xticklabelsvisible = false
ax3D.yticklabelsvisible = false
ax3D.zticklabelsvisible = false
ax3D.xlabelvisible = false
ax3D.ylabelvisible = false
ax3D.zlabelvisible = false
ax3D.azimuth = azimuth
ax3D.elevation = elevation

ax2 = Axis(fig[1, 3]; width = 300) # this axis has values of principal lengths
ax2.title = "ellipsoid axes"
tightlimits!(ax2, Bottom())
sobs = Observable(s)
barplot!(ax2, 1:3, sobs; color = to_color.(COLORSCHEME[1:3]))

display(fig)

# %%
record(
        fig, joinpath(@__DIR__, "volume_growth_$(systemtitle).mp4"), 1:N;
        framerate
        ) do i
    DynamicalSystems.step!(tinteg, dt)
    Y = get_deviations(tinteg)
    ellipsobs[], sobs[] = ellipsoid_from_pointcloud(Y)
    # s[] = principals(newus)
    # autolimits!(ax2)
    # must set new limits centered on u
    # maxε = maximum(norm(u) for u in newus)
    # maxε = maximum(s)
    # ax3D.limits = map(j -> (-maxε, +maxε), (1,2,3))
    # Update trajectory and title and elevation
    u = get_state(tinteg)
    dotobs[] = [u]
    trobs[] = push!(trobs[], u)
    t[] = tinteg.t
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Test, ScientificTypes, ScientificTypesBase, Random
using Tables, CategoricalArrays, DataFrames
using ColorTypes
using Dates
# using CSV # dropped until julia release new LTS as issue for 1.0
import Distributions

const Arr  = AbstractArray
const CArr = CategoricalArray
const Cat  = CategoricalValue
const Vec  = AbstractVector
const Dist = Distributions

include("type_tests.jl")
include("basic_tests.jl")
include("extra_tests.jl")
include("autotype.jl")
include("extra_coercion.jl")

include("coverage.jl")


text:
abstract type AbstractSplitCritereon end

struct GiniSplitCritereon <: AbstractSplitCritereon end
struct NumericAttributeObserver{T<:AbstractSplitCritereon}
    distribution_per_class::Dict{Int, GaussianEstimator}
    num_split_suggestions::Int64
    split_critereon::T
end

abstract type AbstractConditionalTest end

struct NumericBinaryTest <: AbstractConditionalTest
    attr_idx::Int64
    split_value::Float64
end

struct AttributeSplitSuggestion
    split_test::AbstractConditionalTest
    resulting_class_distribution::Tuple{Dict{Int64, Float64}, Dict{Int64, Float64}}
    merit::Float64
    num_splits::Int64
end

function NumericAttributeObserver()
    return NumericAttributeObserver(Dict{Int, GaussianEstimator}(), 10, GiniSplitCritereon())
end

function update(o::NumericAttributeObserver, x::Float64, y::Int64)
    if !haskey(o.distribution_per_class, y)
        o.distribution_per_class[y] = GaussianEstimator()
    end
    o.distribution_per_class[y] = add_weight(o.distribution_per_class[y], x, 1.0)
    return
end

function range(o::NumericAttributeObserver)
    min_seen = typemax(Float64)
    max_seen = typemin(Float64)
    for estimator in values(o.distribution_per_class)
        min_seen = min(min_seen, estimator.min_seen)
        max_seen = max(max_seen, estimator.max_seen)
    end
    return min_seen, max_seen
end

"""
Calculate the gini impurity of a split as the weighted sum
of the gini impurity values of each child node.
"""
function get_split_merit(critereon::GiniSplitCritereon, pre_split_distribution, post_split_distribution)
    child_weights = [sum(collect(values(child_dist))) for (child_idx, child_dist) in enumerate(post_split_distribution)]
    total_weight = sum(child_weights)

    gini = 0.0
    for (child_idx, child_dist) in enumerate(post_split_distribution)
        child_relative_weight = child_weights[child_idx] / total_weight
        child_gini = compute_gini(child_dist, child_weights[child_idx])
        gini += child_relative_weight * child_gini
    end

    # Return Inverse to a low gini gets a high merit.
    return 1.0 - gini
end

"""
Compute the gini impurity of a given distribution of class weights.
The squared proportion of each class is subtracted from 1,
so a gini impurity of 0.0 is acheived when only one class is present.
"""
function compute_gini(class_weight_distribution, total_weight)
    gini = 1.0
    if total_weight > 0.0
        for class_weight in values(class_weight_distribution)
            class_proportion = class_weight / total_weight
            gini -= class_proportion^2
        end
    end
    return gini
end

function get_range_of_merit(critereon::GiniSplitCritereon)
    return 1.0
end

"""
Propose values an attribute could be split at, between the min and max seen.
"""
function get_split_point_suggestions(o::NumericAttributeObserver)
    min_seen, max_seen = range(o)
    split_difference = (max_seen - min_seen) / (o.num_split_suggestions + 1)
    split_values = Vector{Float64}()
    for i in 1:o.num_split_suggestions
        split_value = min_seen + (split_difference * i)
        if split_value > min_seen
            push!(split_values, split_value)
        end
    end
    return split_values
end

function get_binary_split_class_distribution(o::NumericAttributeObserver, split_value::Float64)
    lhs_dist = Dict{Int64, Float64}()
    rhs_dist = Dict{Int64, Float64}()
    for (class, estimator) in o.distribution_per_class
        if split_value < estimator.min_seen
            rhs_dist[class] = estimator.weight_sum
        elseif split_value > estimator.max_seen
            lhs_dist[class] = estimator.weight_sum
        else
            l_dist, e_dist, r_dist = estimated_weight_split(estimator, split_value)
            lhs_dist[class] = l_dist+e_dist
            rhs_dist[class] = r_dist
        end
    end
    return (lhs_dist, rhs_dist)
end

function get_best_split_suggestion(o::NumericAttributeObserver, pre_split_distribution::Dict{Int64, Float64}, att_idx::Int64)
    best_suggestion = (typemin(Float64), 0.0, (Dict{Int64, Float64}(), Dict{Int64, Float64}()))
    suggested_split_values = get_split_point_suggestions(o)
    for split_value in suggested_split_values
        post_split_dist = get_binary_split_class_distribution(o, split_value)
        merit = get_split_merit(o.split_critereon, pre_split_distribution, post_split_dist)
        if merit > best_suggestion[1]
            best_suggestion = (merit, split_value, post_split_dist)
        end
    end
    split_test = NumericBinaryTest(att_idx, best_suggestion[2])
    return AttributeSplitSuggestion(split_test, best_suggestion[3], best_suggestion[1], 2)
end


text:
using LDUFacts: piv_to_perm, perm_to_piv, LDUPerm

res(A, la) = la.P' * Matrix(A) * la.P - la.L * la.D * la.U

@testset "throwing" begin
   A = [0.0 1; 1 0]
   @test_throws SingularException ldu(A, NoPivot())
   @test_throws SingularException ldu(A, DiagonalPivot())
   @test rank(ldu(A, FullPivot())) == 2
end

@testset "ldu(3,2) NoPivot" begin
    A = Matrix(Symmetric([1. 0 2; 0 1 1; 0 0 1]))
    la = ldu(A, NoPivot(); check = false)
    @test issuccess(la)
    @test rank(la) == 3
    @test la.p == [1, 2, 3]
    @test norm(res(A, la)) <= 1e-13
    @test Matrix(la.P) == I(3)
    @test propertynames(la) ⊇ (:U, :L, :D, :d, :p, :P)
end
@testset "ldu(3,2) DiagonalPivot" begin
    A = Symmetric([1. 0 2; 0 1 1; 0 0 1])
    la = ldu(A, DiagonalPivot(); check = false)
    @test issuccess(la)
    @test rank(la) == 3
    @test la.p == [1, 3, 2]
    @test norm(res(A, la)) <= 1e-13
    @test Matrix(la.P) == [1. 0 0; 0 0 1; 0 1 0]
end
@testset "ldu(3,2) FullPivot" begin
    A = Hermitian([1. 0 2; 0 1 1; 0 0 1])
    la = ldu(A, FullPivot(); check = false)
    @test issuccess(la)
    @test rank(la) == 3
    @test la.p == [1, 2, 3]
    @test norm(res(A, la)) <= 1e-13
    @test Matrix(la.P) == [1. 0 0; 0 1 0; 1 0 1]
end

@testset "ldu!(7,10), $ps" for ps in (NoPivot(), DiagonalPivot(), FullPivot())
    R = [
         0  -1   0   0   0  -1  -1
        -1   1   0   0  -1  -1   0
         1   1   1   1  -1   1   0
         1   0   1   0   0   0   1
         0   0   0  -1   0   0   0
         0  -1  -1   0  -1  -1   0
         0   0   0   1   1   0  -1
        -1   1  -1  -1   1   0  -1
         0   1  -1  -1  -1   0   0
        -1  -1   0   1   0   1  -1
    ]
    R = float.(R)
    A = R * Diagonal([1,1,1,1,-1,-1,-1]) * R'
    la = ldu(copy(A), ps, tol=-1, check=false)
    @test norm(res(A, la)) < 1e-14
    ps isa NoPivot || @test issuccess(la)
    @test rank(la) == 7
end

@testset "ldu!(rational)" begin
    A = [20 4 11 -23 20; 4 -38 -25 -46 -1; 11 -25 12 -7 27; -23 -46 -7 10 14; 20 -1 27 14 -3] .// 4
    la = ldu!(copy(A), FullPivot(), tol = 0)
    @test rank(la) == 5
    @test issuccess(la)
    @test iszero(res(A, la))
end

@testset "simple complex" begin
    A = [
        -3//1+0//1*im   3//2-2//1*im  -1//1+2//1*im
         3//2+2//1*im   2//1+0//1*im  -3//2-1//2*im
        -1//1-2//1*im  -3//2+1//2*im  -3//2+0//1*im
    ]
    la = ldu(A, DiagonalPivot(), tol = 0)
    @test issuccess(la)
    @test rank(la) == 3
    @test iszero(res(A, la))
    @test iszero(inv(A) - inv(la))
    @test det(A) == det(la)
end

@testset "advanced complex and rational" begin
    A = [
        2//1+0//1*im   2//1-3//2*im   0//1-2//1*im   3//1+3//1*im  -1//1+1//1*im
        2//1+3//2*im   1//1+0//1*im  -1//2+1//2*im  -3//2+2//1*im   1//1+0//1*im
        0//1+2//1*im  -1//2-1//2*im  -1//1+0//1*im   1//1-3//2*im  -3//2-3//2*im
        3//1-3//1*im  -3//2-2//1*im   1//1+3//2*im   3//1+0//1*im   1//1+1//2*im
       -1//1-1//1*im   1//1+0//1*im  -3//2+3//2*im   1//1-1//2*im  -1//1+0//1*im
    ]
    A = big.(A)
    A = Hermitian(A, :L)
    la = ldu(A, FullPivot(), tol = 0)
    @test issuccess(la)
    @test rank(la) == 5
    @test iszero(res(A, la))
    @test inv(la) == inv(A)
    @test logdet(la) ≈ logdet(A)
    b = big.([1+0im, 2, 3, 4, 5])
    @test la \ b == A \ b
    @test b' / la == b' / A
end

@testset "PivotLike" begin
    A = [20 4 11 -23 20; 4 -38 -25 -46 -1; 11 -25 12 -7 27; -23 -46 -7 10 14; 20 -1 27 14 -3] .// 4
    la = ldu(A, FullPivot())
    ps = PivotLike(la)
    lb = ldu(A, ps)
    @test la.piv == lb.piv
    @test la.pam == lb.pam
    @test la.ame == lb.ame
end

@testset "perm_to_piv" begin
    perm = [2, 4, 3, 9, 5, 1, 10, 7, 6, 8]
    piv = perm_to_piv(perm)
    @test piv_to_perm(piv) == perm
end

@testset "mut_by_revperm" begin
    P = LDUPerm([1, 2], [2, 0], [1+2im, 0])
    b = [10+0im, 20]
    MP = Matrix(P)
    @test MP == [1 0; 1+2im 1]
    @test P * b == MP * b
    @test P * [b b] == MP * [b b]
    @test [b'; b'] * P' == [b'; b'] * MP'
end


text:
# global config

const DATA_PATH = "../resources/datasets"
const MODEL_PATH = "../resources/models"
const INIT_LR = 0.003
const CHANNELS = 3

const NUM_CLASSES = 10


text:
function _default_tween_creator()
    # println("Pool: used default tween creator")
    Tween()
end

mutable struct TweenPool
    tweens::Array{AbstractTween,1}
    creator::Function # f() -> AbstractTween

    function TweenPool(creator::Function = _default_tween_creator)
        o = new()
        o.tweens = Array{AbstractTween,1}()
        o.creator = creator
        o
    end
end

function count(pool::TweenPool)
    length(pool.tweens)
end

function get_from!(pool::TweenPool)
    tween = if count(pool) == 0
        pool.creator()
    else
        pop!(pool.tweens)
    end

    reset!(tween)

    tween
end

function put_to!(pool::TweenPool, tween::AbstractTween)
    push!(pool.tweens, tween)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
module MDMKP
include("../src/JOH.jl")
using JuMP #used to create the MOI model

export MDMKP_Prob, MDMKP_Sol, load_folder, create_MIPS_model

struct Problem_ID <: JOH.ProblemID
	id::Int
	dataset::Int
	instance::Int
	case::Int
	tightness::Float64
	n_vars::Int
	n_demands::Int
	n_dimensions::Int
	mixed_obj::Bool
	loosened::Bool
	category::Char
end

struct MDMKP_Prob <: JOH.Problem
    objective::Vector{Int}
    upper_bounds::Vector{Tuple{Vector{Int},Int}}
    lower_bounds::Vector{Tuple{Vector{Int},Int}}
	id::Problem_ID
end

struct MDMKP_Sol <: JOH.Solution
	problem::MDMKP_Prob
	value::BitArray
	score::Number
	objective::Number
	infeasibility::Number
end

include("load_folder.jl")

include("decision_tree.jl")

"""Select slices of MDMKP problems"""
function slice_select(problems::Vector{MDMKP_Prob};
		datasets=1:9, cases=1:6, instances=1:15)::Vector{MDMKP_Prob}
	extracted = []
	for prob in problems
		if prob.id.dataset in datasets &&
				prob.id.case in cases &&
				prob.id.instance in instances
			push!(extracted, prob)
		end
	end
	extracted
end

"""Solution Constructor """
function MDMKP_Sol(bitlist::AbstractArray{Bool}, problem::MDMKP_Prob)
    upper_bounds_totals = [sum(coeffs .* bitlist) for (coeffs, bound) in
		problem.upper_bounds]
    lower_bounds_totals = [sum(coeffs .* bitlist) for (coeffs, bound) in
		problem.lower_bounds]

	#add dimension constraint infeasibility
    infeasibility = sum([
		upper_bounds_totals[i] > upper_bound ?
			upper_bounds_totals[i] - upper_bound : 0
    	for (i, (constraint_coeffs, upper_bound)) in enumerate(
				problem.upper_bounds)
	])

	#add demand constraint infeasibility
    infeasibility += sum([
		lower_bounds_totals[i] < lower_bound ?
			lower_bound - lower_bounds_totals[i] : 0
        for (i, (contraint_coeffs, lower_bound)) in enumerate(
				problem.lower_bounds)
	])

	#calculate the objective function
    objective_value = sum(problem.objective .* bitlist)

	#the score is the objective function if feasible, else infeasibility
    score = infeasibility > 0 ? -infeasibility : objective_value

    MDMKP_Sol(
		problem,
        bitlist,
        score,
        objective_value,
        infeasibility
    )
end

"""Accept an MDMKP problem, and return a formulation that includes heavily
penalized artificial variables to make the discovery of a feasibile solution
trivial. """
function create_MIPS_model(problem::MDMKP_Prob, optimizer;
		time_limit=20,
		weight=1000,
		num_threads=6)::Model
	model = Model(optimizer)
	_create_MIPS_model(problem, model, time_limit, weight, num_threads)
end
function create_MIPS_model(problem::MDMKP_Prob;
		time_limit=20,
		weight=1000,
		num_threads=6)::Model
	model = Model()
	_create_MIPS_model(problem, model, time_limit, weight, num_threads)
end
function direct_model(problem::MDMKP_Prob, model;
		time_limit=20,
		weight=1000000,
		num_threads=6)::Model
	_create_MIPS_model(problem, model, time_limit, weight, num_threads)
end

function _create_MIPS_model(problem, model, time_limit, weight, num_threads)::Model
	JOH.Matheur.set_threads!(model, num_threads)
	JOH.Matheur.set_time!(model, time_limit)

    #make the problem variables with a Binary constraint
    @variable(model, x[1:length(problem.objective)], Bin)

	#make the artificial variables to fix infeasibility
	@variable(model, s[1:length(problem.upper_bounds)] <= 0)
	@variable(model, ss[1:length(problem.lower_bounds)] >= 0)

	# objective is the normal objective value minus the artificial values
	# needed to make feasible
    @objective(model, Max,
		sum(problem.objective .* x) - weight * (-sum(s)+sum(ss)))

	# add dimension constraints
    for (i, ub) in enumerate(problem.upper_bounds)
        @constraint(model, sum(ub[1] .* x) + s[i] <= ub[2])
    end

	# add demand constraints
    for (i, lb) in enumerate(problem.lower_bounds)
        @constraint(model, sum(lb[1] .* x) + ss[i] >= lb[2])
    end

    model
end

end

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
const DEFAULT_COMMANDS  = Dict("help"     => Help.help,
                               "type"     => Keyboard.type,
                               "email"    => Email.email,
                               "internet" => Internet.internet)

"""
    start()
    start(<keyword arguments>)

Start offline, low latency, highly accurate and secure speech to command translation.

# Keyword arguments
- `commands::Dict{String, <:Any}=DEFAULT_COMMANDS`: the commands to be recognized with their mapping to a function or to a keyboard key or shortcut.
- `subset::AbstractArray{String}=nothing`: a subset of the `commands` to be recognised and executed (instead of the complete `commands` list).
- `max_speed_subset::AbstractArray{String}=nothing`: a subset of the `commands` for which the command names (first word of a command) are to be recognised with maxium speed rather than with maximum accuracy. Forcing maximum speed is usually desired for single word commands that map to functions or keyboard shortcuts that should trigger immediate actions as, e.g., mouse clicks or page up/down (in general, actions that do not modify content and can therefore safely be triggered at maximum speed). Note that forcing maximum speed means not to wait for a certain amount of silence after the end of a command as normally done for the full confirmation of a recognition. As a result, it enables a minimal latency between the saying of a command name and its execution. Note that it is usually possible to define very distinctive command names, which allow for a safe command name to shortcut mapping at maximum speed (to be tested case by case).
- `modeldirs::Dict{String, String}=DEFAULT_MODELDIRS`: the directories where the unziped speech recognition models to be used are located. Models are downloadable from here: https://alphacephei.com/vosk/models
- `noises::Dict{String, <:AbstractArray{String}}=DEFAULT_NOISES`: for each model, an array of strings with noises (tokens that are to be ignored in the speech as, e.g., "huh").
- `audio_input_cmd::Cmd=nothing`: a command that returns an audio stream to replace the default audio recorder. The audio stream must fullfill the following properties: `samplerate=$SAMPLERATE`, `channels=$AUDIO_IN_CHANNELS` and `format=Int16` (signed 16-bit integer).

# Submodules for command name to function mapping
- [`Help`](@ref)
- [`Keyboard`](@ref)
- [`Mouse`](@ref)
- [`Email`](@ref)
- [`Internet`](@ref)

To see a description of a submodule, type `?<modulename>`.

# Default `commands`
```
$(pretty_dict_string(DEFAULT_COMMANDS))
```

# Default `modeldirs`
```
$(pretty_dict_string(DEFAULT_MODELDIRS))
```

# Default `noises`
```
$(pretty_dict_string(DEFAULT_NOISES))
```

# Examples

#### Define `subset`
```
# Listen only to the commands "help" and "type".
using JustSayIt
start(subset=["help", "type"])
```

#### Define custom `commands` - functions and keyboard shortcuts - and a `max_speed_subset`
```
using JustSayIt
commands = Dict("help"      => Help.help,
                "type"      => Keyboard.type,
                "ma"        => Mouse.click_left,
                "select"    => Mouse.press_left,
                "okay"      => Mouse.release_left,
                "middle"    => Mouse.click_middle,
                "right"     => Mouse.click_right,
                "double"    => Mouse.click_double,
                "triple"    => Mouse.click_triple,
                "copy"      => (Key.ctrl, 'c'),
                "cut"       => (Key.ctrl, 'x'),
                "paste"     => (Key.ctrl, 'v'),
                "undo"      => (Key.ctrl, 'z'),
                "redo"      => (Key.ctrl, Key.shift, 'z'),
                "upwards"   => Key.page_up,
                "downwards" => Key.page_down,
                );
start(commands=commands, max_speed_subset=["ma", "select", "okay", "middle", "right", "double", "triple", "copy", "upwards", "downwards"])
```

#### Define custom `modeldirs`
```
using JustSayIt
modeldirs = Dict(DEFAULT_MODEL_NAME => "$(homedir())/mymodels/vosk-model-small-en-us-0.15",
                 TYPE_MODEL_NAME    => "$(homedir())/mymodels/vosk-model-en-us-daanzu-20200905")
start(modeldirs=modeldirs)
```

#### Define `audio_input_cmd`
```
# Use a custom command to create the audio input stream - instead of the default recorder (the rate, channels and format must not be chosen different!)
using JustSayIt
audio_input_cmd = `arecord --rate=$SAMPLERATE --channels=$AUDIO_IN_CHANNELS --format=S16_LE`
start(audio_input_cmd=audio_input_cmd)
```
"""
function start(; commands::Dict{String, <:Any}=DEFAULT_COMMANDS, subset::Union{Nothing, AbstractArray{String}}=nothing, max_speed_subset::Union{Nothing, AbstractArray{String}}=nothing, modeldirs::Dict{String,String}=DEFAULT_MODELDIRS, noises::Dict{String,<:AbstractArray{String}}=DEFAULT_NOISES, audio_input_cmd::Union{Cmd,Nothing}=nothing) where N
    if (!isnothing(subset) && !issubset(subset, keys(commands))) @IncoherentArgumentError("'subset' incoherent: the obtained command name subset ($(subset)) is not a subset of the command names ($(keys(commands))).") end
    if (!isnothing(max_speed_subset) && !issubset(max_speed_subset, keys(commands))) @IncoherentArgumentError("'max_speed_subset' incoherent: the obtained max_speed_subset ($(max_speed_subset)) is not a subset of the command names ($(keys(commands))).") end
    if !isnothing(subset) commands = filter(x -> x[1] in subset, commands) end
    if isnothing(max_speed_subset) max_speed_subset = String[] end
    max_speed_token_subset = map(max_speed_subset) do cmd_name # Only the first tokens are used to decide if max speed is used.
       string(first(split(cmd_name)))
    end
    max_speed_multiword_cmds = [x for x in keys(commands) if any(startswith.(x, [x for x in max_speed_token_subset if x ∉ max_speed_subset]))]
    incoherent_subset = [x for x in max_speed_multiword_cmds if x ∉ max_speed_subset]
    if !isempty(incoherent_subset) @IncoherentArgumentError("'max_speed_subset' incoherent: the following commands are not part of 'max_speed_subset', but start with the same word as a command that is part of it: \"$(join(incoherent_subset,"\", \"", " and "))\". Adjust the 'max_speed_subset' to prevent this.") end

    # Initializations
    @info "JustSayIt: I am initializing (say \"sleep JustSayIt\" to put me to sleep; press CTRL+c to terminate)..."
    init_jsi(commands, modeldirs, noises)
    start_recording(; audio_input_cmd=audio_input_cmd)

    # Interprete commands
    @info "Listening for commands..."
    valid_cmd_names = command_names()
    cmd_name = ""
    is_sleeping = false
    use_max_speed = true
    try
        while true
            if is_sleeping
                if cmd_name == COMMAND_NAME_AWAKE
                    if is_confirmed() is_sleeping = false end
                    if (is_sleeping) @info("I heard \"awake\". If you want to awake me, say \"awake JustSayIt\".")
                    else             @info("... awake again. Listening for commands...")
                    end
                end
            else
                if cmd_name in valid_cmd_names
                    cmd = command(string(cmd_name))
                    if (t0_latency() > 0.0) @debug "Latency of command `$cmd_name`: $(round(toc(t0_latency()),sigdigits=2))." end
                    try
                        latency_msg = use_max_speed ? " (latency: $(round(Int,toc(t0_latency())*1000)) ms)" : ""
                        @info "Starting command: $(pretty_cmd_string(cmd))"*latency_msg
                        execute(cmd)
                    catch e
                        if isa(e, InsecureRecognitionException)
                            @info("Command `$cmd_name` aborted: insecure command argument recognition.")
                        elseif isa(e, InterruptException)
                            @info "Command `$cmd_name` aborted (with CTRL+C)."
                        else
                            rethrow(e)
                        end
                    end
                elseif cmd_name == COMMAND_NAME_SLEEP
                    if is_confirmed() is_sleeping = true end
                    if (!is_sleeping) @info("I heard \"sleep\". If you want me to sleep, say \"sleep JustSayIt\".")
                    else              @info("Going to sleep... (To awake me, just say \"awake JustSayIt\".)")
                    end
                elseif cmd_name != ""
                    @debug "Invalid command: $cmd_name."
                end
            end
            try
                force_reset_previous(recognizer(COMMAND_RECOGNIZER_ID))
                use_max_speed = _is_next(max_speed_token_subset, recognizer(COMMAND_RECOGNIZER_ID), _noises(DEFAULT_MODEL_NAME); use_partial_recognitions=true, ignore_unknown=false)
                cmd_name = next_token(recognizer(COMMAND_RECOGNIZER_ID), _noises(DEFAULT_MODEL_NAME); use_partial_recognitions = use_max_speed, ignore_unknown=false)
                if cmd_name == UNKNOWN_TOKEN # For increased recognition security, ignore the current word group if the unknown token was obtained as command name (achieved by doing a full reset). This will prevent for example "text right" or "text type text" to trigger an action, while "right" or "type text" does so.
                    reset_all()
                    cmd_name = ""
                end
                while (cmd_name != "") && (cmd_name ∉ valid_cmd_names) && any(startswith.(valid_cmd_names, cmd_name))
                    token = next_token(recognizer(COMMAND_RECOGNIZER_ID), _noises(DEFAULT_MODEL_NAME); use_partial_recognitions = use_max_speed, ignore_unknown=false)
                    if token == UNKNOWN_TOKEN # For increased recognition security, ignore the current word group if the unknown token was obtained as command name (achieved by doing a full reset). This will prevent for example "text right" or "text type text" to trigger an action, while "right" or "type text" does so.
                        reset_all()
                        cmd_name = ""
                    else
                        cmd_name = cmd_name * " " * token
                    end
                end
            catch e
                if isa(e, InsecureRecognitionException)
                    if !is_sleeping @debug(e.msg) end
                    cmd_name = ""
                else
                    rethrow(e)
                end
            end
        end
    catch e
        if isa(e, InterruptException)
            @info "Terminating JustSayIt..."
        else
            rethrow(e)
        end
    finally
        stop_recording()
        @info "Stopped listening for commands."
    end

    finalize_jsi()
end

function is_confirmed()
    try
        _is_confirmed()
    catch e
        if isa(e, InsecureRecognitionException)
            return false
        else
            rethrow(e)
        end
    end
end

@voiceargs words=>(valid_input=["just say it"], vararg_timeout=2.0, vararg_max=3) function _is_confirmed(words::String...)
    return (join(words, " ") == "just say it")
end

execute(cmd::Function)                  = cmd()
execute(cmd::PyKey)                     = Keyboard.press_keys(cmd)
execute(cmd::NTuple{N,PyKey} where {N}) = Keyboard.press_keys(cmd...)


text:
export Solution, OrderedSolution, PolicySolution
export get_action, best_action, rand_action

"Abstract solution type."
abstract type Solution end

"Return action for step `t` at `state`."
get_action(sol::Solution, t::Int, state::State) = error("Not implemented.")

"Null solution that indicates no plan was found."
struct NullSolution <: Solution end

## Ordered solutions ##

"Abstract type for ordered planner solutions."
abstract type OrderedSolution <: Solution end

"Return action for step `t`."
get_action(sol::OrderedSolution, t::Int) = error("Not implemented.")
get_action(sol::OrderedSolution, t::Int, ::State) = get_action(sol, t)

# Ordered solutions should support indexing and iteration over action terms
Base.iterate(::OrderedSolution) = error("Not implemented.")
Base.iterate(::OrderedSolution, iterstate) = error("Not implemented.")
Base.getindex(::OrderedSolution, ::Int) = error("Not implemented.")
Base.eltype(::Type{OrderedSolution}) = Term

include("solutions/ordered_plan.jl")

## Policy-based solutions ##

"Abstract type for policy-based solutions."
abstract type PolicySolution <: Solution end

"Return action for the given state."
get_action(sol::PolicySolution, state::State) = best_action(sol, state)
get_action(sol::PolicySolution, ::Int, state::State) = get_action(sol, state)

"Returns the best action for the given state."
best_action(::PolicySolution, ::State) = error("Not implemented.")

"Samples an action according to the policy for the given state."
rand_action(::PolicySolution, ::State) = error("Not implemented.")

"Return value (i.e. expected future reward) of the given state (and action)."
get_value(::PolicySolution, ::State) = error("Not implemented.")
get_value(::PolicySolution, ::State, ::Term) = error("Not implemented.")

"Return Q-values of actions for the given state as an iterator over pairs."
get_action_values(::PolicySolution, ::State) = error("Not implemented.")

include("solutions/policy_value.jl")
include("solutions/random_policy.jl")
include("solutions/boltzmann_policy.jl")
include("solutions/epsilon_greedy.jl")