text:
using ModelingToolkit, Unitful, OrdinaryDiffEq, DiffEqJump
using Test
MT = ModelingToolkit
@parameters τ [unit = u"ms"]
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
D = Differential(t)

@test MT.get_unit(t) == u"ms"
@test MT.get_unit(E) == u"kJ"
@test MT.get_unit(τ) == u"ms"

@parameters β [unit = u"°"] α [unit = u"°C"] γ [unit = 1u"s"]
@test_throws MT.ValidationError MT.get_unit(β)
@test_throws MT.ValidationError MT.get_unit(α)
@test_throws MT.ValidationError MT.get_unit(γ)

unitless = Unitful.unit(1)
@test MT.get_unit(0.5) == unitless
@test MT.get_unit(t) == u"ms"
@test MT.get_unit(P) == u"MW"
@test MT.get_unit(τ) == u"ms"

@test MT.get_unit(τ^-1) == u"ms^-1"
@test MT.equivalent(MT.get_unit(D(E)),u"MW")
@test MT.equivalent(MT.get_unit(E/τ), u"MW")
@test MT.get_unit(2*P) == u"MW"
@test MT.get_unit(t/τ) == unitless
@test MT.equivalent(MT.get_unit(P - E/τ),u"MW")
@test MT.equivalent(MT.get_unit(D(D(E))),u"MW/ms")

@test MT.get_unit(1.0^(t/τ)) == unitless
@test MT.get_unit(exp(t/τ)) == unitless
@test MT.get_unit(sin(t/τ)) == unitless
@test MT.get_unit(sin(1u"rad")) == unitless
@test MT.get_unit(t^2) == u"ms^2"

@test !MT.validate(E^1.5 ~ E^(t/τ))
@test MT.validate(E^(t/τ) ~ E^(t/τ))

eqs = [D(E) ~ P - E/τ
        0 ~ P]
@test MT.validate(eqs)
@named sys = ODESystem(eqs)
@named sys = ODESystem(eqs, t, [P, E], [τ])

@test !MT.validate(D(D(E)) ~ P)
@test !MT.validate(0 ~ P + E*τ)

#Unit-free
@variables x y z u
@parameters σ ρ β
eqs = [0 ~ σ*(y - x)]
@test MT.validate(eqs)

#Array variables
@variables t x[1:3,1:3](t)
D = Differential(t)
eqs = D.(x) .~ x
ODESystem(eqs,name=:sys)

# Array ops
using Symbolics: unwrap, wrap
using LinearAlgebra
@variables t
sts = @variables x[1:3](t) y(t)
ps = @parameters p[1:3] = [1, 2, 3]
D = Differential(t)
eqs = [
       collect(D.(x) ~ x)
       D(y) ~ norm(x)*y
      ]
ODESystem(eqs, t, [sts...;], [ps...;],name=:sys)

#= Not supported yet b/c iterate doesn't work on unitful array
# Array ops with units
@variables t [unit =u"s"]
sts = @variables x[1:3](t) [unit = u"kg"] y(t) [unit = u"kg"]
ps = @parameters b [unit = u"s"^-1]
D = Differential(t)
eqs = [
       collect(D.(x) ~ b*x)
       D(y) ~ b*norm(x)
      ]
ODESystem(eqs, t, [sts...;], [ps...;])

#Array variables with units
@variables t [unit = u"s"] x[1:3,1:3](t) [unit = u"kg"]
@parameters a [unit = u"s"^-1]
D = Differential(t)
eqs = D.(x) .~ a*x
ODESystem(eqs)
=#

#Difference equation with units
@parameters t [unit = u"s"] a [unit = u"s"^-1]
@variables x(t) [unit = u"kg"]
δ = Differential(t)
D = Difference(t; dt = 0.1u"s")
eqs = [
    δ(x) ~ a*x
]
de = ODESystem(eqs, t, [x, y], [a],name=:sys)

@parameters t
@variables y[1:3](t)
@parameters k[1:3]
D = Differential(t)

eqs = [D(y[1]) ~ -k[1]*y[1] + k[3]*y[2]*y[3],
       D(y[2]) ~  k[1]*y[1] - k[3]*y[2]*y[3] - k[2]*y[2]^2,
       0 ~  y[1] + y[2] + y[3] - 1]

@named sys = ODESystem(eqs,t,y,k)

# Nonlinear system
@parameters a [unit = u"kg"^-1]
@variables x [unit = u"kg"]
eqs = [
    0 ~ a*x
]
@named nls = NonlinearSystem(eqs, [x], [a])

# SDE test w/ noise vector
@parameters τ [unit = u"ms"] Q [unit = u"MW"]
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
D = Differential(t)
eqs = [D(E) ~ P - E/τ
      P ~ Q]

noiseeqs = [0.1u"MW",
            0.1u"MW"]
@named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
# With noise matrix
noiseeqs = [0.1u"MW" 0.1u"MW"
            0.1u"MW" 0.1u"MW"]
@named sys = SDESystem(eqs,noiseeqs, t, [P, E], [τ, Q])

noiseeqs = [0.1u"MW" 0.1u"MW"
            0.1u"MW" 0.1u"s"]
@test !MT.validate(eqs,noiseeqs)

#Test non-trivial simplifications
@variables t [unit = u"s"] V(t) [unit = u"m"^3] L(t) [unit = u"m"]
@parameters v [unit = u"m/s"] r [unit =u"m"^3/u"s"]
D = Differential(t)
eqs = [D(L) ~ v,
       V ~ L^3]
@named sys = ODESystem(eqs)
sys_simple = structural_simplify(sys)

eqs = [D(V) ~ r,
       V ~ L^3]
@named sys = ODESystem(eqs)
sys_simple = structural_simplify(sys)

@variables V [unit = u"m"^3] L [unit = u"m"]
@parameters v [unit = u"m/s"] r [unit = u"m"^3/u"s"] t [unit = u"s"]
eqs = [V ~ r*t,
       V ~ L^3]
@named sys = NonlinearSystem(eqs, [V, L], [t, r])
sys_simple = structural_simplify(sys)

eqs = [L ~ v*t,
       V ~ L^3]
@named sys = NonlinearSystem(eqs, [V,L], [t,r])
sys_simple = structural_simplify(sys)

#Jump System
@parameters β [unit = u"(mol^2*s)^-1"] γ [unit = u"(mol*s)^-1"] t [unit = u"s"] jumpmol [unit = u"mol"]
@variables S(t) [unit = u"mol"] I(t) [unit = u"mol"] R(t) [unit = u"mol"]
rate₁   = β*S*I
affect₁ = [S ~ S - 1*jumpmol, I ~ I + 1*jumpmol]
rate₂   = γ*I
affect₂ = [I ~ I - 1*jumpmol, R ~ R + 1*jumpmol]
j₁      = ConstantRateJump(rate₁, affect₁)
j₂      = VariableRateJump(rate₂, affect₂)
js      = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ],name=:sys)

affect_wrong = [S ~ S - jumpmol, I ~ I + 1]
j_wrong      = ConstantRateJump(rate₁, affect_wrong)
@test_throws MT.ValidationError JumpSystem([j_wrong, j₂], t, [S, I, R], [β, γ],name=:sys)

rate_wrong   = γ^2*I
j_wrong     = ConstantRateJump(rate_wrong, affect₂)
@test_throws MT.ValidationError JumpSystem([j₁, j_wrong], t, [S, I, R], [β, γ],name=:sys)

# mass action jump tests for SIR model
maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])
maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
@named js3  = JumpSystem([maj1, maj2], t, [S,I,R], [β,γ])

#Test unusual jump system
@parameters β γ t
@variables S(t) I(t) R(t)

maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
maj2 = MassActionJump(γ, [S => 1], [S => -1])
@named js4  = JumpSystem([maj1, maj2], t, [S], [β, γ])


text:
@testset "Intensity" begin

    @testset "Constant Intensity" begin
        bg = 0.5
        kappa = 0.5
        kern(x) = pdf.(Distributions.Exponential(1/0.5), x)
        testEvents = [1, 2, 3]
        inferTests && @inferred HawkesProcesses.intensity(0, testEvents, bg, kappa, kern)
        @test HawkesProcesses.intensity(0, testEvents, bg, kappa, kern) == bg
        @test HawkesProcesses.intensity(1, testEvents, bg, kappa, kern) == bg
        @test HawkesProcesses.intensity(1000, testEvents, bg, kappa, kern) == bg
        @test HawkesProcesses.intensity([0, 1000], testEvents, bg, kappa, kern) == [bg, bg]
        @test HawkesProcesses.intensity(testEvents, testEvents, bg, kappa, kern) ≈ [bg, bg + kappa * kern(1), bg + kappa * kern(1) + kappa * kern(2)] atol=0.01
    end

    @testset "Kernel Distribution Intensity" begin
        bg = 0.5
        kappa = 0.5
        kern = Distributions.Exponential(1/0.5)
        testEvents = [1, 2, 3]
        inferTests && @inferred HawkesProcesses.intensity([0], testEvents, bg, kappa, kern)

        int = HawkesProcesses.intensity([1,2,3], testEvents, bg, kappa, kern)
        @test length(int) == 3
    end

    @testset "Functional Background Intensity" begin
        bg(x) = 0.5*x
        kappa = 0.5
        kern = Distributions.Exponential(1/0.5)
        testEvents = [1, 2, 3]
        inferTests && @inferred HawkesProcesses.intensity([0], testEvents, bg, kappa, kern)

        int = HawkesProcesses.intensity([1,2,3], testEvents, bg, kappa, kern)
        @test length(int) == 3
    end

end

ext: jl


text:
# Base Quality
# ============
#
# A representation of positions-specific integer quality scores, as in FASTQ.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

immutable QualityEncoding
    # valid ASCII character range
    ascii::NTuple{2,UInt8}

    # valid base quality range
    qual::NTuple{2,Int8}

    function QualityEncoding(ascii::Range{Char}, qual::Range{Int})
        if length(ascii) != length(qual)
            throw(ArgumentError("the ranges of ASCII and base quality don't match"))
        end
        return new((first(ascii), last(ascii)), (first(qual), last(qual)))
    end
end

# Return the offset of converting a base quality to an ASCII code (quality + offset).
function ascii_offset(encoding::QualityEncoding)
    return encoding.ascii[1] - encoding.qual[1]
end

"Sanger (Phred+33) quality score encoding"
const SANGER_QUAL_ENCODING     = QualityEncoding('!':'~', 0:93)

"Solexa (Solexa+64) quality score encoding"
const SOLEXA_QUAL_ENCODING     = QualityEncoding(';':'~', -5:62)

"Illumina 1.3 (Phred+64) quality score encoding"
const ILLUMINA13_QUAL_ENCODING = QualityEncoding('@':'~', 0:62)

"Illumina 1.5 (Phred+64) quality score encoding"
const ILLUMINA15_QUAL_ENCODING = QualityEncoding('B':'~', 2:62)

"Illumina 1.8 (Phred+33) quality score encoding"
const ILLUMINA18_QUAL_ENCODING = QualityEncoding('!':'~', 0:93)

# Check quality encoding of `input[start:stop]`.
function check_quality_string(encoding, input, start, stop)
    ascii_lo, ascii_hi = encoding.ascii
    for i in start:stop
        @inbounds a = input[i]
        if !(ascii_lo ≤ a ≤ ascii_hi)
            error("base quality '$(Char(a))' is out of range")
        end
    end
end

# Decode an ASCII string `input[start:stop]` into base quality scores `output`.
function decode_quality_string!(encoding, input, output, start, stop)
    resize!(output, stop - start + 1)
    offset = ascii_offset(encoding)
    for i in 1:(stop - start + 1)
        @inbounds output[i] = input[start + i - 1] - offset
    end
    return output
end

# Encode base quality scores `input` into a ASCII string `output`.
function encode_quality_string!(encoding, input, output, start, stop)
    resize!(output, stop - start + 1)
    offset = ascii_offset(encoding)
    for i in 1:(stop - start + 1)
        @inbounds output[i] = input[start + i - 1] + offset
    end
    return output
end


text:
#Demo 1: Simple sine wave plot
#-------------------------------------------------------------------------------
using InspectDR
using Colors
import NumericIO: UEXPONENT_SI

#==Input
===============================================================================#
#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
tcycle = 1e-9
npts = 2000 #Relatively large dataset
ncycles = 4
tmax = tcycle*ncycles
Δ = tmax/(npts-1)
t = collect(0:Δ:tmax)
y = sin(t*(ncycles*2pi/tmax))

t[end] = t[1] #x-values are no longer ordered

#Lower resolution t-vector:
t_lres = collect(0:(tmax/10):tmax)

#==Generate plot
===============================================================================#
plot = InspectDR.Plot2D(title="Sample Plot (λ)")
plot.layout = InspectDR.Layout(fontname="monospace", fontscale=1.5)
plot.layout.grid = grid(vmajor=true, vminor=true, hmajor=true)
plot.layout.legend.enabled = true
plot.layout.legend.width = 150
plot.layout.showtimestamp = true
plot.layout.xlabelformat.expdisplay = UEXPONENT_SI #Use SI notation on x-axis

style = :dashdot #solid/dashdot/...
wfrm = add(plot, t, y+1, id="sin(2πt)+1")
	wfrm.line = line(color=blue, width=1, style=style)
wfrm = add(plot, t, y-1, id="sin(2πt)-1")
	wfrm.line = line(color=red, width=5, style=style)
wfrm = add(plot, t_lres, -2+4*(t_lres./tmax), id="-2+4t/tmax")
	wfrm.line = line(color=blue, width=3, style=:dash)
	wfrm.glyph = glyph(shape=:*, size=10)
#= Supported shapes:
	:square, :diamond,
	:uarrow, :darrow, :larrow, :rarrow, #usually triangles
	:cross, :+, :diagcross, :x,
	:circle, :o, :star, :*,
=#

a = plot.annotation
a.xlabel = "Time (s)"
a.ylabel = "Signal Voltage (V)"

gplot = display(InspectDR.GtkDisplay(), plot)

InspectDR.write_png("export_plot.png", plot)
InspectDR.write_svg("export_plot.svg", plot)
InspectDR.write_eps("export_plot.eps", plot)
InspectDR.write_pdf("export_plot.pdf", plot)

:DONE


text:
module CrossSections
using Unitful
const UF = Unitful
using Unitful: MeV, NoUnits, cm, m, J, s, kg
using Random
using Distributions

include("constants.jl")
include("compton.jl")

end # module


text:
// Our struct to use in the array
typedef struct Foo_t{
  int bar;
} * Foo;

int main(){
  // Init struct array
  Foo[] foos = new Foo[10];

  // Init all values (default value is null)
  int i = 0;
  while(i < foos.length){
    foos[i] = new Foo_t;
    foos[i]->bar = i;
    i++;
  }

  // Print all values
  printInt(foos.length);
  for(Foo foo : foos){
     printInt(foo->bar);
  }
  return 0;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
mutable struct Particles{T<: AbstractFloat, J<: Integer}
    x::Vector{Vector{T}} # particle positions
    v::Vector{Vector{T}} # particle velocities
    p::Vector{Vector{T}} # best particle positions
    g::Vector{Vector{T}} # best neighborhood swarm position

    lb::Vector{T} #lower bound
    ub::Vector{T} #uper bound

    swarmsize::J #swarm size
    particlelength::J # particle length

    fx::Vector{T} # current particle function values
    fs::Vector{Bool} # feasibility of each particle
    fp::Vector{T} # best particle function values
    fg::Vector{T} # best neighborhood swarm position value

    order::Vector{J} # neighborhood order
    best_position::Vector{T} # best position
    best_value::T # best position value
    history::Vector{T} # best position value history

    function Particles(swarmsize::V, lb::J, ub::J) where {V<:Integer, J<:Vector}
        @assert length(ub) == length(lb) "ub and lb must have same dimension"
        @assert all(ub .> lb) "Each value of 'ub' must be greater than 'lb'"

        S = length(lb)
        x = [ rand(S) for i = 1:swarmsize]
        v = [ rand(S) for i = 1:swarmsize]

        vhigh = abs.(ub .- lb)
        vlow = -vhigh
        x = map(x-> lb .+ x .* (ub .- lb), x) # particle positions
        v = map(x-> vlow .+ x .* (vhigh .- vlow), v) # particle velocities
        p = copy(x)  # best particle positions
        g = copy(x)  # best neighborhood swarm position

        new{typeof(rand()), typeof(swarmsize) }(x, v, p, g, lb, ub, swarmsize, S, fill(Inf,swarmsize),
         fill(false,swarmsize), fill(Inf,swarmsize), fill(Inf,swarmsize), collect(1:swarmsize),
          [Inf for i=1:S], Inf, [] );
    end
end


text:
using Test, AutoARIMA

@testset "Forecast" begin

    @test forecast(MAModel{1}(0.0, [-0.7], 1.0), [3.0,8,2,5,6]) ≈ 1.01941

    @test forecast(ARModel{2}(1.0, [0.5, 0.5], 1.0), [3.0,2]) ≈ 3.5

    @test forecast(ARModel{2}(10.0, [0.5, 0.5], 1.0), [3.0,2]) ≈ 12.5

end


text:
export KernelInterface
"""
Objects can be selected and deselected for certain operations (in the future). Supertype of [`AbstractComposite`](@ref).
"""
abstract type Selectable end

"""
Interface which `KERNEL` types implement. See [`AbstractComposite`](@ref) and
[`DataFrameSystem`](@ref) which both implement `KernelInterface`.
"""
abstract type KernelInterface <: Selectable end

function select(x::Selectable)
    x.selected_ = true
end

function deselect(x::Selectable)
    x.selected_ = false
end

function isSelected(x::Selectable)
    return x.selected_
end


text:
function blockMatching(X::Array{Float32,3}, Par::PAR, Neighbor_arr::Array{Int32,2}, Num_arr::Array{Int32,1}, SelfIndex_arr::Array{Int32,1})::Array{Int32,2}

    L = length(Num_arr)

    Init_Index = zeros(Int32, Par.patnum, L)

    for  i in 1 : L
        Patch = X[:, SelfIndex_arr[i], :]

        Neighbors = X[:, Neighbor_arr[1:Num_arr[i], i], :]

        Dist = zeros(Float32,size(Neighbors)[2])
        for j in 1:size(Neighbors)[2]
            Dist[j] = sum(float32.(Neighbors[:,j,:] - Patch).^2)
        end

        index = Int32.(sortperm(Dist, alg = QuickSort))

        Init_Index[:, i] = Neighbor_arr[index[1:Par.patnum], i]
    end

    Init_Index
end


text:
"""
    Metropolis(y, x, β₀, β, M)

A Metropolis algorithm for probit regression.

# Arguments
- `y::Vector{Bool}`: boolean vector of outcomes
- `x::Matrix{<:Real}`: real valued vector or matrix of inputs
- `β₀::MvNormal`: a univariate or multivariate normal distribution for the prior
- `β::MvNormal`: a univariate or multivariate normal distribution for the candidate
- `M::Integer`: the number of draws

# Output
- `chain::Matrix{<:Float}`: matrix of draws from the target distribution
- `accept::Vector{Bool}`: boolean vector of acceptance indicators

The function can be called without 'x' in which case a constant-only model is estimated.
"""

function estimate(Time::Vector{T} where T<:Real,
                  Indicator::BitVector,
                  X::Matrix{T} where T<:Real,
                  β::MvNormal,
                  α::Exponential,
                  draws = 10_000::Integer)

    N, J = size(X)
    @assert length(Time) == N "Time and X must have the same number of rows."
    @assert length(Indicator) == N "Indicator and X must have the same number of rows."
    @assert J == length(β.μ) "The number of columns in x must match the dimension of β."

    # Initial values
    b = zeros(draws, J)
    a = zeros(draws)

    b[1, :] = β.μ
    a[1] = α.θ

    Ts = copy(Time)
    Ss = findall(Indicator .== 1)
    for i in 2:draws
        # Propose new values
        b[i, :], a[i] = proposal(b[i-1, :], a[i-1], β, α, Ts, X)

        # Time - update missing values (i.e. censored values)
        for n in Ss
            λ = exp(-dot(X[n, :], b[i, :]) * a[i])
            Ts[n] = rand(Truncated(Weibull(a[i], λ), Time[n], Inf))
        end
    end
    return (b, a)
end

function lambda(b::Vector, a::AbstractFloat, X::Matrix)
    m = X * b
    λ = exp.(-m * a)
    return λ
end

function loglik(b::Vector, a::AbstractFloat, Time::Vector, X::Matrix)
    λ = lambda(b, a, X)
    ll = sum(logpdf.(Weibull.(a, λ), Time))
    return ll
end

function logprior(prior, x)
    return sum(logpdf(prior, x))
end

function proposal(b, a, β, α, Time, X)
    # Existing
    pr = logprior(α, a) + logprior(β, b) # Log-prior
    ll = loglik(b, a, Time, X)           # Log-likelihood
    po = pr + ll                         # Log-posterior

    # Proposal
    a1 = exp(rand(Normal(log(a), 0.05))) # Proposal for alpha
    b1 = rand.(Normal.(b, 0.05))         # Proposal for beta

    pr1 = logprior(α, a1) + logprior(β, b1) # Log-prior
    ll1 = loglik(b1, a1, Time, X)           # Log-likelihood
    po1 = pr1 + ll1                         # Log-posterior

    if log(rand()) < (po1 - po)
        return (b1, a1)
    else
        return (b, a)
    end
end


text:
using DeepForest, Flux
using StatsBase: mode, mean
using Test

@testset "Leaf precompute" begin
    l = DeepForest.Leaf()
    x = Float32.([1.0 2.0; 3.0 4.0])
    y = [1, 1]
    DeepForest.precompute!(l, x, y)
    @test l.pred == 1
    y = [2, 2]
    DeepForest.precompute!(l, x, y)
    @test l.pred == 2
end

@testset "Leaf predict" begin
    l = DeepForest.Leaf()
    x = Float32.([1.0 2.0; 3.0 4.0])
    y = [1, 1]
    DeepForest.precompute!(l, x, y)
    @test DeepForest.predict(l, x) == [1, 1]
end

@testset "Node init" begin
    feat_sub = Dict(
        1 => [1, 2]
    )

    n = Node(feat_sub, 5, 1, 1)
    @test n.subset == [1, 2]
    @test typeof(n.left) == DeepForest.Leaf

    feat_sub = Dict(
        1 => [1, 2],
        2 => [2, 3],
        3 => [1, 3]
    )

    n = Node(feat_sub, 5, 2, 1)
    @test n.subset == [1, 2]
    @test n.left.subset == [2, 3]
    @test n.right.subset == [1, 3]
end

@testset "Node parameters" begin
    feat_sub = Dict(
        1 => [1, 2],
        2 => [2, 3],
        3 => [1, 3]
    )

    n = Node(feat_sub, 5, 2, 1)
    @test params(n) != Flux.Params([])
end

@testset "Node precompute" begin
    feat_sub = Dict(
        1 => [1, 2],
        2 => [2, 3],
        3 => [1, 3]
    )

    n = Node(feat_sub, 5, 2, 1)
    x = rand(Float32, 3, 4)
    y = [1, 2, 2, 2]
    DeepForest.precompute!(n, x, y)
    @test n.best.first != 1 || n.best.second != 1
end

@testset "Node predict" begin
    feat_sub = Dict(
        1 => [1, 2],
        2 => [2, 3],
        3 => [1, 3]
    )

    n = Node(feat_sub, 5, 2, 1)
    x = rand(Float32, 3, 4)
    y = [1, 2, 2, 2]
    DeepForest.precompute!(n, x, y)
    @test length(predict(n, x)) == 4
    @test length(unique(predict(n, x))) <= 2
end

@testset "Node loss" begin
    feat_sub = Dict(
        1 => [1, 2],
        2 => [2, 3],
        3 => [1, 3]
    )

    n = Node(feat_sub, 5, 2, 1)
    x = rand(Float32, 3, 4)
    y = [1, 2, 2, 2]
    DeepForest.precompute!(n, x, y)

    @test loss!(n, x, y) != 0

    gs = gradient(() -> loss!(n, x, y), params(n))
    @test gs != IdDict()
end

@testset "Node training" begin
    x = rand(Float32, 2, 1000);
    x[1, :] .*= 2*π;
    x[1, :] .-= π;
    x[2, :] .*= 3;
    x[2, :] .-= 1.5;
    y = Int64.(x[2, :] .< sin.(x[1, :])) .+ 1;

    features = Dict(
        7 => [1, 2],
        6 => [1, 2],
        5 => [1, 2],
        4 => [1, 2],
        3 => [1, 2],
        2 => [1, 2],
        1 => [1, 2]
    )

    model = Node(features, 10, 2, 1)

    tree_train!(1000, model, x, y)
    @test mean(predict(model, x) .== y) > .9
end

@testset "Node importances" begin
    x = rand(Float32, 2, 1000);
    x[1, :] .*= 2*π;
    x[1, :] .-= π;
    x[2, :] .*= 3;
    x[2, :] .-= 1.5;
    y = Int64.(x[2, :] .< sin.(x[1, :])) .+ 1;

    features = Dict(
        7 => [1, 2],
        6 => [1, 2],
        5 => [1, 2],
        4 => [1, 2],
        3 => [1, 2],
        2 => [1, 2],
        1 => [1, 2]
    )

    model = Node(features, 10, 2, 1)

    tree_train!(1000, model, x, y)
    imp = importance(model, x)
    @test imp[2] >= imp[1]
end


text:
using NLSolvers, StaticArrays, Test
@testset "mixed optimization problems" begin
    function theta(x)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end
    f(x) = 100.0 * ((x[3] - 10.0 * theta(x))^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2

    function f∇f!(x, ∇f)
        if !(∇f == nothing)
            if (x[1]^2 + x[2]^2 == 0)
                dtdx1 = 0
                dtdx2 = 0
            else
                dtdx1 = -x[2] / (2 * pi * (x[1]^2 + x[2]^2))
                dtdx2 = x[1] / (2 * pi * (x[1]^2 + x[2]^2))
            end
            ∇f[1] =
                -2000.0 * (x[3] - 10.0 * theta(x)) * dtdx1 +
                200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[1] / sqrt(x[1]^2 + x[2]^2)
            ∇f[2] =
                -2000.0 * (x[3] - 10.0 * theta(x)) * dtdx2 +
                200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[2] / sqrt(x[1]^2 + x[2]^2)
            ∇f[3] = 200.0 * (x[3] - 10.0 * theta(x)) + 2.0 * x[3]
        end

        fx = f(x)
        objective_return(fx, ∇f)
    end

    function f∇f(x, ∇f)
        if !(∇f == nothing)
            ∇f = similar(x)
        end
        fx, ∇f = f∇f!(gx, x)
        objective_return(fx, ∇f)
    end
    function f∇fs(x, ∇f)
        if !(∇f == nothing)
            if (x[1]^2 + x[2]^2 == 0)
                dtdx1 = 0
                dtdx2 = 0
            else
                dtdx1 = -x[2] / (2 * pi * (x[1]^2 + x[2]^2))
                dtdx2 = x[1] / (2 * pi * (x[1]^2 + x[2]^2))
            end

            s1 =
                -2000.0 * (x[3] - 10.0 * theta(x)) * dtdx1 +
                200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[1] / sqrt(x[1]^2 + x[2]^2)
            s2 =
                -2000.0 * (x[3] - 10.0 * theta(x)) * dtdx2 +
                200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[2] / sqrt(x[1]^2 + x[2]^2)
            s3 = 200.0 * (x[3] - 10.0 * theta(x)) + 2.0 * x[3]
            ∇f = @SVector [s1, s2, s3]
        end
        objective_return(f(x), ∇f)
    end
    obj_inplace = OnceDiffed(f∇f!)
    obj_outofplace = OnceDiffed(f∇f)
    obj_static = OnceDiffed(f∇fs)

    x0 = [-1.0, 0.0, 0.0]
    xopt = [1.0, 0.0, 0.0]
    x0s = @SVector [-1.0, 0.0, 0.0]

    println("Starting  from: ", x0)
    println("Targeting     : ", xopt)

    res = minimize(obj_inplace, copy(x0), NelderMead(), MinOptions())
    print("NN  $(summary(NelderMead()))         ")
    @printf(
        "%2.2e  %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.nm_obj, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), SimulatedAnnealing(), MinOptions())
    print("NN  $(summary(SimulatedAnnealing()))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = HZ()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=HZ())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = CD()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=CD())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = HS()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=HS())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = FR()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=FR())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = PRP()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=PRP(;plus=false))))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(
        obj_inplace,
        copy(x0),
        ConjugateGradient(update = PRP(plus = true)),
        MinOptions(),
    )
    print("NN  $(summary(ConjugateGradient(update=PRP(;plus=true))))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = VPRP()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=VPRP())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = LS()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=LS())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(update = DY()), MinOptions())
    print("NN  $(summary(ConjugateGradient(update=DY())))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )

    for _method in (GradientDescent, LBFGS, BFGS, DBFGS, DFP, SR1)
        # for _method in (GradientDescent, BFGS, DBFGS, DFP, SR1)
        methodtxt = summary(_method())
        for m in (Inverse(), Direct())
            mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
            if _method == LBFGS && m isa Inverse
                res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
            elseif _method !== LBFGS
                res = minimize(obj_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN  $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize(obj_static, x0s, LineSearch(_method(m)), MinOptions())
                print("NN  $_method(S) $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                if m isa Direct && !(_method == GradientDescent)
                    println("Trust region: NWI")
                    res = minimize!(
                        obj_inplace,
                        copy(x0),
                        TrustRegion(_method(Direct()), NWI()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e  %d\n",
                        norm(res.info.solution - xopt, Inf),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                    println("Trust region: NTR")
                    res = minimize!(
                        obj_inplace,
                        copy(x0),
                        TrustRegion(_method(Direct()), NTR()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e  %d\n",
                        norm(res.info.solution - xopt, Inf),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                end
            end
        end
    end
    println()

    res = minimize(obj_inplace, x0, LineSearch(BFGS(Inverse())), MinOptions())
    @test res.info.iter == 30
    @printf(
        "NN  BFGS    (inverse): %2.2e  %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.∇fz, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize!(
        obj_inplace,
        copy(x0),
        LineSearch(BFGS(Inverse()), Backtracking()),
        MinOptions(),
    )
    @test res.info.iter == 30
    @printf(
        "NN! BFGS    (inverse): %2.2e  %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.∇fz, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize!(
        obj_inplace,
        copy(x0),
        LineSearch(BFGS(Inverse()), Backtracking(interp = FFQuadInterp())),
        MinOptions(),
    )
    @test res.info.iter == 30
    @printf(
        "NN! BFGS    (inverse, quad): %2.2e %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.∇fz, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_static, x0s, LineSearch(BFGS(Inverse())), MinOptions())
    @test res.info.iter == 30
    @printf(
        "NN  BFGS(S) (inverse): %2.2e %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.∇fz, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(
        obj_static,
        x0s,
        LineSearch(BFGS(Inverse()), Backtracking(interp = FFQuadInterp())),
        MinOptions(),
    )
    @test res.info.iter == 30
    @printf(
        "NN  BFGS(S) (inverse, quad): %2.2e %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.∇fz, Inf),
        res.info.minimum,
        res.info.iter
    )

    x0 = rand(3)
    x0s = SVector{3}(x0)
    println("\nFrom a random point: ", x0)
    res = minimize(obj_inplace, copy(x0), NelderMead(), MinOptions())
    print("NN  $(summary(NelderMead()))         ")
    @printf(
        "%2.2e  %2.2e %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        norm(res.info.nm_obj, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), SimulatedAnnealing(), MinOptions())
    print("NN  $(summary(SimulatedAnnealing()))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )
    res = minimize(obj_inplace, copy(x0), ConjugateGradient(), MinOptions())
    print("NN  $(summary(ConjugateGradient()))         ")
    @printf(
        "%2.2e  %2.2e %d\n",
        norm(res.info.solution - xopt, Inf),
        res.info.minimum,
        res.info.iter
    )

    for _method in (GradientDescent, LBFGS, BFGS, DBFGS, DFP, SR1)
        methodtxt = summary(_method())
        for m in (Inverse(), Direct())
            mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
            if _method == LBFGS && m isa Inverse
                res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
            elseif _method !== LBFGS
                res = minimize(obj_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN  $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize!(obj_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize(obj_static, x0s, LineSearch(_method(m)), MinOptions())
                print("NN  $_method(S) $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e %d\n",
                    norm(res.info.solution - xopt, Inf),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                if m isa Direct && !(_method == GradientDescent)
                    println("Trust region: NWI")
                    res = minimize!(
                        obj_inplace,
                        copy(x0),
                        TrustRegion(_method(Direct()), NWI()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e %d\n",
                        norm(res.info.solution - xopt, Inf),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                    println("Trust region: NTR")
                    res = minimize!(
                        obj_inplace,
                        copy(x0),
                        TrustRegion(_method(Direct()), NTR()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e %d\n",
                        norm(res.info.solution - xopt, Inf),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                end
            end
        end
    end
    println()

    function himmelblau!(x, ∇f)
        if !(∇f == nothing)
            ∇f[1] =
                4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 -
                14.0
            ∇f[2] =
                2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 -
                28.0 * x[2]
        end

        fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
        objective_return(fx, ∇f)
    end

    function himmelblaus(x, ∇f)
        fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
        if !(∇f == nothing)
            ∇f1 =
                4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 -
                14.0
            ∇f2 =
                2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 -
                28.0 * x[2]
            ∇f = @SVector([∇f1, ∇f2])
        end
        objective_return(fx, ∇f)
    end

    function himmelblau(x, ∇f)
        g = ∇f == nothing ? ∇f : similar(x)

        return himmelblau!(x, g)
    end

    him_inplace = OnceDiffed(himmelblau!)
    him_static = OnceDiffed(himmelblaus)
    him_outofplace = OnceDiffed(himmelblau)

    println("\nHimmelblau function")
    x0 = [3.0, 1.0]
    x0s = SVector{2}(x0)
    minimizers =
        [[3.0, 2.0], [-2.805118, 3.131312], [-3.779310, -3.283186], [3.584428, -1.848126]]
    for _method in (GradientDescent, LBFGS, BFGS, DBFGS, DFP, SR1)
        # for _method in (GradientDescent, BFGS, DBFGS, DFP, SR1)
        methodtxt = summary(_method())
        for m in (Inverse(), Direct())
            mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
            if _method == LBFGS && m isa Inverse
                res = minimize!(him_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
            elseif _method !== LBFGS
                res =
                    minimize(him_outofplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN  $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize!(him_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize(him_static, x0s, LineSearch(_method(m)), MinOptions())
                print("NN  $_method(S) $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                if m isa Direct && !(_method == GradientDescent)
                    println("Trust region: NWI")
                    res = minimize!(
                        him_inplace,
                        copy(x0),
                        TrustRegion(_method(Direct()), NWI()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e  %d\n",
                        minimum([
                            norm(res.info.solution - xopt, Inf) for xopt in minimizers
                        ]),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                    println("Trust region: NTR")
                    res = minimize!(
                        him_inplace,
                        copy(x0),
                        TrustRegion(_method(Direct()), NTR()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e  %d\n",
                        minimum([
                            norm(res.info.solution - xopt, Inf) for xopt in minimizers
                        ]),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                end
            end
        end
    end
    println()
    println()
    xrand = rand(2)
    xrands = SVector{2}(xrand)
    println("\nFrom a random point: ", xrand)

    for _method in (GradientDescent, LBFGS, BFGS, DBFGS, SR1)
        # for _method in (GradientDescent, BFGS, DBFGS, DFP, SR1)
        methodtxt = summary(_method())
        for m in (Inverse(), Direct())
            mtxt = m isa Inverse ? "(inverse): " : "(direct):  "
            if _method == LBFGS && m isa Inverse
                res = minimize!(him_inplace, copy(x0), LineSearch(_method(m)), MinOptions())
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
            elseif _method !== LBFGS
                res = minimize(
                    him_outofplace,
                    copy(xrand),
                    LineSearch(_method(m)),
                    MinOptions(),
                )
                print("NN  $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize!(
                    him_inplace,
                    copy(xrand),
                    LineSearch(_method(m)),
                    MinOptions(),
                )
                print("NN! $_method    $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                res = minimize(him_static, xrands, LineSearch(_method(m)), MinOptions())
                print("NN  $_method(S) $mtxt")
                @printf(
                    "%2.2e  %2.2e %2.2e  %d\n",
                    minimum([norm(res.info.solution - xopt, Inf) for xopt in minimizers]),
                    norm(res.info.∇fz, Inf),
                    res.info.minimum,
                    res.info.iter
                )
                if m isa Direct && !(_method == GradientDescent)
                    println("Trust region: NWI")
                    res = minimize!(
                        him_inplace,
                        copy(xrand),
                        TrustRegion(_method(Direct()), NWI()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e  %d\n",
                        minimum([
                            norm(res.info.solution - xopt, Inf) for xopt in minimizers
                        ]),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                    println("Trust region: NTR")
                    res = minimize!(
                        him_inplace,
                        copy(xrand),
                        TrustRegion(_method(Direct()), NTR()),
                        MinOptions(),
                    )
                    print("NN! $_method    $mtxt")
                    @printf(
                        "%2.2e  %2.2e %2.2e  %d\n",
                        minimum([
                            norm(res.info.solution - xopt, Inf) for xopt in minimizers
                        ]),
                        norm(res.info.∇fz, Inf),
                        res.info.minimum,
                        res.info.iter
                    )
                end
            end
        end
    end
    println()

end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
# downsample.jl

using MIRT: downsample_dim1
using MIRT: downsample1
using MIRT: downsample2
using MIRT: downsample3
using MIRT: downsample
import MIRT: downsample3_loop, downsample3_perm

using Test: @test, @inferred

"""
    downsample_dim1_test()
"""
function downsample_dim1_test()
	@inferred downsample_dim1(1:4, 2)
	x = reshape(2:2:48, 4, 6)
#	y = @inferred downsample_dim1(x, 2) # todo: fails
	y = downsample_dim1(x, 2)
	@test y == reshape(3:4:47, 2, 6)
	true
end

"""
downsample1_test()
"""
function downsample1_test()
	x = 2:2:48
	y = @inferred downsample1(x, 2)
	@test y == 3:4:47
	@inferred downsample1(1:5, 2, warn=false) # warn
	true
end

"""
downsample2_test()
"""
function downsample2_test()
	x = reshape(0:2:46, 4, 6)
	y = @inferred downsample2(x, 2)
	@test y == [5 21 37; 9 25 41]
	true
end

"""
downsample3_test()
"""
function downsample3_test()
	x = [6, 5, 2]
	x = reshape(2*(1:prod(x)), x...)
	for down = 1:2
		y = downsample3(x, down)
	#	y = @inferred downsample3(x, down) # todo fails
		if down == 1
			@test y == x
		elseif down == 2
			@test y == reshape([39 63; 43 67; 47 71], 3,2,1) # squeeze
		end
	end

	x = ones(4, 6, 8)
	@test downsample3_loop(x, (2, 2, 2)) == downsample3_perm(x, (2, 2, 2))

	true
end

#=
if has_aspire # todo
	filex = [test_dir filesep 'testx.fld']
	filey = [test_dir filesep 'testy.fld']
	fld_write(filex, x)
	delete(filey)
	com = ['op sample3 mean ' filey ' ' filex ' %d %d %d']
	com = sprintf(com, down, down, down)
	os_run(com)
	z = fld_read(filey)
	!isequal(y, z) && throw("aspire/matlab mismatch")
end
=#

@test downsample_dim1_test()
@test downsample1_test()
@test downsample2_test()
@test downsample3_test()


text:
# -*- mode: julia; -*-
module WAVPlay
import WAV.wavplay
using Compat
immutable PulseSample
    format::Int32
    rate::UInt32
    channels::UInt8
end

function wavplay(data, fs)
    ss = PulseSample(5,fs,size(data,2))
    data = convert(Array{Float32}, data)

    err = Cint[0]
    s = ccall((:pa_simple_new, "libpulse-simple"),
        Ptr{Void}, (Ptr{Void}, Ptr{UInt8}, Cint, Ptr{Void}, Ptr{UInt8},
        Ptr{PulseSample}, Ptr{Void}, Ptr{Void}, Ptr{Cint}),
        0,"hey",1,0,"playback",&ss,0,0,err)
    assert(s != C_NULL)

    assert(0 == ccall((:pa_simple_write, "libpulse-simple"), Int32,
        (Ptr{Void}, Ptr{Void}, Cssize_t, Ptr{Cint}),
        s, data, sizeof(data), err))

    assert(0 == ccall((:pa_simple_drain, "libpulse-simple"), Int32,
        (Ptr{Void}, Ptr{Cint}), s, err))
end
end # module


text:
module PyVirtualenv

using Base: PkgId, UUID
using Libdl

const _pycall_deps_jl = Ref{Module}()

""" Get path to `PyCall/deps/deps.jl` """
function pycall_deps_jl_path()
    modpath = Base.locate_package(Base.identify_package("PyCall"))
    if modpath === nothing
        error("PyCall not found")
    end
    @debug "PyCall.jl found at: $modpath"
    return joinpath(dirname(modpath), "..", "deps", "deps.jl")
end
# e.g., ~/.julia/dev/PyCall/deps/deps.jl

""" Load `PyCall/deps/deps.jl` namespace. """
function pycall_deps_jl()
    if isdefined(_pycall_deps_jl, 1)
        return _pycall_deps_jl[]
    end
    mod = Module(:__PyCall_deps__)
    Base.include(mod, pycall_deps_jl_path())
    _pycall_deps_jl[] = mod
    @debug "PyCall/deps/deps.jl:\n$(sprint(show_deps_jl, mod))"
    return mod
end

function show_deps_jl(io, mod)
    for name in sort(names(mod; all=true))
        val = try
            getproperty(mod, name)
        catch
            continue
        end
        if val isa Union{String, Bool}
            println(name, " =\t", val)
        end
    end
end

_clength(x::Cstring) = ccall(:strlen, Csize_t, (Cstring,), x) + 1
_clength(x) = length(x)

function __leak(::Type{T}, x) where T
    n = _clength(x)
    ptr = ccall(:malloc, Ptr{T}, (Csize_t,), n * sizeof(T))
    unsafe_copyto!(ptr, pointer(x), n)
    return ptr
end

"""
    _leak(T::Type, x::AbstractString) :: Ptr
    _leak(x::Array) :: Ptr

Leak `x` from Julia's GCer.
"""
_leak(x::Union{Cstring, Array}) = __leak(eltype(x), x)
_leak(T::Type, x::AbstractString) =
    _leak(Base.unsafe_convert(T, Base.cconvert(T, x)))
_leak(::Type{Cwstring}, x::AbstractString) =
    _leak(Base.cconvert(Cwstring, x))

function pythonhome_of(pyprogramname::AbstractString)
    if Sys.iswindows()
        script = """
        import sys
        sys.stdout.write(sys.exec_prefix)
        """
        # [[~/.julia/dev/PyCall/deps/build.jl::PYTHONHOME]]
    else
        script = """
        import sys
        sys.stdout.write(sys.prefix)
        sys.stdout.write(":")
        sys.stdout.write(sys.exec_prefix)
        """
        # https://docs.python.org/3/using/cmdline.html#envvar-PYTHONHOME
    end
    cmd = `$pyprogramname -c $script`
    @debug "Trying to find out PYTHONHOME." cmd

    # For Windows:
    env = copy(ENV)
    env["PYTHONIOENCODING"] = "UTF-8"
    cmd = setenv(cmd, env)

    return read(cmd, String)
end

function libpython_handle()
    libpython = pycall_deps_jl().libpython
    return Libdl.dlopen(libpython,
                        Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)
end

function Py_IsInitialized(handle = libpython_handle())
    return ccall(Libdl.dlsym(handle, :Py_IsInitialized), Cint, ()) != 0
end

function is_same_env(pyprogramname::AbstractString,
                     PYTHONHOME::AbstractString = pythonhome_of(pyprogramname))
    if !Py_IsInitialized()
        @debug "Py_IsInitialized() = false"
        return false
    end
    @debug "Loading PyCall..."
    @eval import PyCall
    return Base.invokelatest(_is_same_env, PyCall, pyprogramname, PYTHONHOME)
end

function _is_same_env(PyCall, pyprogramname, PYTHONHOME)
    sys = PyCall.pyimport("sys")
    loaded_executable = sys[:executable]
    prefix = sys[:prefix]
    exec_prefix = sys[:exec_prefix]
    if Sys.iswindows()
        loaded_pythonhome = "$exec_prefix"
    else
        loaded_pythonhome = "$prefix:$exec_prefix"
    end
    @debug """
    Given:
        pyprogramname = $pyprogramname
        PYTHONHOME    = $PYTHONHOME
    From Python interpreter:
        sys.executable  = $loaded_executable
        sys.prefix      = $prefix
        sys.exec_prefix = $exec_prefix
    """
    if Sys.iswindows()
        # Py_SetProgramName does not work on Windows?
        return loaded_pythonhome == PYTHONHOME
    end
    return loaded_executable == pyprogramname &&
        loaded_pythonhome == PYTHONHOME
end

"""
    activate(pyprogramname::String, [PYTHONHOME::String])

Activate virtual environment of Python program at `pyprogramname`; it MUST
be compatible with the Python program configured for PyCall.jl.
"""
function activate(pyprogramname::AbstractString,
                  PYTHONHOME::AbstractString = pythonhome_of(pyprogramname))
    handle = libpython_handle()

    fp(f::Symbol) = Libdl.dlsym(handle, f)

    already_inited = Py_IsInitialized(handle)
    if already_inited
        if is_same_env(pyprogramname, PYTHONHOME)
            @debug "Re-activating the same environment (which is OK)."
            return
        end
        error("Re-activation not supported.")
    end

    is_py2 = pycall_deps_jl().pyversion_build.major < 3

    @debug "Py_SetPythonHome($PYTHONHOME)"
    if is_py2
        ccall(fp(:Py_SetPythonHome), Cvoid, (Cstring,),
              _leak(Cstring, PYTHONHOME))
    else
        ccall(fp(:Py_SetPythonHome), Cvoid, (Ptr{Cwchar_t},),
              _leak(Cwstring, PYTHONHOME))
    end

    @debug "Py_SetProgramName($pyprogramname)"
    if is_py2
        ccall(fp(:Py_SetProgramName), Cvoid, (Cstring,),
              _leak(Cstring, pyprogramname))
    else
        ccall(fp(:Py_SetProgramName), Cvoid, (Ptr{Cwchar_t},),
              _leak(Cwstring, pyprogramname))
    end

    @debug "Py_InitializeEx()"
    ccall(fp(:Py_InitializeEx), Cvoid, (Cint,), 0)
end

function pipenv_python(path::AbstractString)
    dir = isdir(path) ? path : dirname(path)
    return rstrip(read(setenv(`pipenv --py`; dir=dir), String))
end

"""
    activate_pipenv([path::String = pwd()])
    activate_pipenv(module::Module)
    activate_pipenv(package::Symbol)

**WARNING**  Same restriction as `activate` applies here.

# Examples
```
julia> PyVirtualenv.activate_pipenv("PATH/TO/PROJECT/Pipfile")
julia> PyVirtualenv.activate_pipenv("PATH/TO/PROJECT/")  # equivalent
```

If you have `Pipfile` in the same directory as in `Project.toml`, you can
pass the module object to this function:

```
julia> using MyModule
julia> PyVirtualenv.activate_pipenv(MyModule)
```

However, if you use `PyCall` inside `MyModule`, above call to `activate_pipenv`
is too late.  One way of working around this problem is:

```
module MyModule

using PyVirtualenv

function __init__()
    PyVirtualenv.activate_pipenv(@__MODULE__)
    @eval using PyCall
end

end
```

Another way is to manually run `activate_pipenv(:MyModule)` (note that the
argument is a `Symbol`) before importing `MyModule`.

```
julia> PyVirtualenv.activate_pipenv(:MyModule)
julia> using MyModule
```
"""
function activate_pipenv(path::AbstractString = pwd())
    activate(pipenv_python(path))
end

activate_pipenv(m::Union{Module, Symbol}) = activate_pipenv(_pathof(m))

function _pathof(m)
    path = __pathof(m)
    if path === nothing
        error("Cannot find the path of module $m.")
    end
    return path
end

__pathof(m::Module) = pathof(m)
__pathof(m::Symbol) = Base.find_package(string(m))

const pycall_pkgid =
    PkgId(UUID("438e738f-606a-5dbb-bf0a-cddfbfd45ab0"), "PyCall")

is_pycall_loaded() = haskey(Base.loaded_modules, pycall_pkgid)

function loaded_pycall()
    if !is_pycall_loaded()
        error("PyCall not loaded.")
    end
    return Base.loaded_modules[pycall_pkgid]
end

"""
    pyexecfile(path, globals, [locals])

Works like Python 2 `execfile` except that `globals` is always required.
"""
pyexecfile(path::AbstractString, args...) =
    _pyexecfile(loaded_pycall(), path, args...)

function _pyexecfile(PyCall, path::AbstractString,
                     globals::AbstractDict,
                     locals::AbstractDict = globals)
    globals = convert(PyCall.PyDict, globals)
    locals = convert(PyCall.PyDict, locals)
    if PyCall.pyversion.major < 3
        PyCall.pybuiltin("execfile")(path, globals, locals)
    else
        src = read(path, String)
        code = PyCall.pybuiltin("compile")(src, path, "exec")
        PyCall.pybuiltin("exec")(code, globals, locals)
    end
end

"""
    activate_this(activate_this_py::String)

See:
https://virtualenv.pypa.io/en/stable/userguide/#using-virtualenv-without-bin-python

# Examples
```
julia> PyVirtualenv.activate_this("VIRTUAL_ENV/bin/activate_this.py")
```
"""
function activate_this(activate_this_py::AbstractString)
    pyexecfile(activate_this_py, Dict("__file__" => activate_this_py))
    return nothing
end

end  # module


text:
using JuMP, Gurobi, Cbc
import XLSX
using BenchmarkTools

# make sure u set the working directory to current files folder
data = XLSX.readxlsx("../data/part1/data_numP100_numC25_10_05_23_17.xlsx")

#General
General = data["General"];
l = General["B1"]
num_providers = General["B2"]
num_clients = General["B3"]

#Providers
Providers = data["Providers"];
provCostString  = string("A2:A", num_providers+1)
provStartString = string("B2:B", num_providers+1)
f = Providers[provCostString]
w = Providers[provStartString]

#Providers
Clients = data["Clients"];
clientStartString  = string("A2:A", num_clients+1)
clientDurationString  = string("B2:B", num_clients+1)
s_times = Clients[clientStartString]
d = Clients[clientDurationString]

#Distances
Distances = data["Distances"];
t = Distances[Distances.dimension]

# need to convert to floats to use 'sortperm' below
w = convert(Array{Float64,1}, w[:,1])
s_times = convert(Array{Float64,1}, s_times[:,1])
d = convert(Array{Float64,1}, d[:,1])
f = convert(Array{Float64,1}, f[:,1])

client_sorted = sortperm(s_times)  # client start times sorted (in index)
arr = [[f[i], w[i]] for i in 1:length(w)]
provider_sorted = sortperm(f) # providers cost and start time sorted

# heuristic function that computes feasible routes (but likely not optimal, i.e. minimum cost)
function heuristics_result(provider_sorted, client_sorted, w, s_times, t, d)
    past_client = zeros(0)
    num_served = zeros(num_providers, 2)
    num_served[:,1] = provider_sorted

    for m in provider_sorted  # for each index in sorted provider, and then for each client
        current_location = 1  # 1 is at depot(14HS) currently
        for n in client_sorted
            # only proceed if this client has not been assigned to a provider yet.
            if n ∉ past_client
                currentTime_provider = w[m]
                currentTime_client = s_times[n]
                travelTime = t[current_location, n+1] # n+1 is nth client in the t matrix

                #time to travel from provider's current location to the client plus the time when the provider is free.
                time_to_arrive = travelTime + currentTime_provider

                #Time to serve client "n".
                duration = d[n]
                #Time finished serving client "n".
                finish_time = currentTime_client + duration

                #Ensure the provider can arrive in time (at or earlier than service start times)
                # Must return back the origin (14HS) by hour l at latest.
                if currentTime_client >= time_to_arrive && finish_time <= l - t[n+1, 1]
                        println("Provider: ", m, " is assigned to client: ", n)
                        num_served[findfirst(isequal(m), num_served[:,1]), 2] += 1
                        w[m] = finish_time # update the time of provider m
                        current_location = n+1  # n+1 is nth client in the t-matrix
                        append!(past_client, n) # record service provided for client n
                end
            end
        end
        # end the algorithm if we have served all clients. so O(N*<M)?
        if length(past_client) == length(client_sorted)
            break
        end
    end
    println("\n")
    return convert(Array{Int,2}, num_served), convert(Array{Int,1}, past_client)
end

function compute_cost(num_served)
    cost = 0
    for i in 1:length(num_served[:,2])
        if num_served[:,2][i] > 0
            cost += f[num_served[:,1][i]]
        end
    end
    return cost
end

function is_feasible(past_client)
    return length(past_client) == length(client_sorted)
end

function get_heuristics_path(num_served, past_client, k)
    path = zeros(1)
    k_ind = findfirst(isequal(k), num_served[:,1])
    if k_ind == 1
        for i in 1:num_served[1,2]
            append!(path, past_client[i])
        end
    end
    if k_ind != 1
        start_ind = sum(num_served[1:k_ind-1,2]) + 1
        for i in start_ind:start_ind+num_served[k_ind,2]-1
            append!(path, past_client[i])
        end
    end
    return convert(Array{Int,1}, append!(path, zeros(1))) # convert to integer array
end

num_served, past_client = heuristics_result(provider_sorted, client_sorted, copy(w), s_times, t, d)
num_served
println(client_sorted)

#@btime heuristics_result(provider_sorted, client_sorted, copy(w), s_times, t, d)

length(past_client)
is_feasible(past_client)
compute_cost(num_served)

function get_path(num_served, past_client, k)
    path = get_heuristics_path(num_served, past_client, k)
    path += convert(Array{Int,1}, ones(length(get_heuristics_path(num_served, past_client, k))))
    return (path)
end

function get_workers(num_served)
    workers = zeros(0)
    for i in 1:size(num_served)[1]
        if num_served[i,:][2] != 0
            append!(workers, num_served[i,:][1])
        end
    end
    return workers
end

path_dict = Dict(trunc(Int, get_workers(num_served)[i]) => get_path(num_served, past_client, get_workers(num_served)[i]) for i in 1:length(get_workers(num_served)))

################### MIP ######################################
################### MIP ######################################
################### MIP ######################################
################### MIP ######################################
################### MIP ######################################
using JuMP, Gurobi, Cbc
import XLSX
using BenchmarkTools

#### read data ####
#xf = XLSX.readxlsx("../data/part1/data_numP10_numC10_09_05_13_29.xlsx")
xf = XLSX.readxlsx("../data/part1/data_numP200_numC150_11_05_34_00.xlsx")

l = xf["General!B1"]     # work day duration
N = xf["General!B3"] + 1  # number of clients + 1
M = xf["General!B2"]  # number of providers
# + 1 is for including the headquarter (HS14)

# Clients
s = xf["Clients!A2:A$N"]   # client service start times
d = vcat([0], xf["Clients!B2:B$N"])   # durations for clients, set d[1] = 0 for HS14

# Providers
w = xf["Providers!B2:B$(M+1)"]  # start time of providers
f = xf["Providers!A2:A$(M+1)"]  # hiring cost of providers

# make separate s vector for each provider, inserting into the first element the provider's work start time
s_test = ones(M, N)
for k in 1:M
    global s_test[k,:] = vcat([w[k]], s)
end
s = s_test

# Distances
t = xf["Distances"][:]  # 2-d array of travel times
##########

model = Model();
set_optimizer(model, Gurobi.Optimizer);
set_time_limit_sec(model, 300.0);

#(13)
@variable(model, x[1:N, 1:N, 1:M], Bin)

#################### setting start values from heursitic solution ###########

for i in 1:N
    for j in 1:N
        for k in 1:M
            set_start_value(x[i,j,k], 0)
        end
    end
end

for key in keys(path_dict)
    path = path_dict[key]
    for i in 1:length(path)-1
        set_start_value(x[path[i], path[i+1], key], 1)
    end
end

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

# subtour constraint variables using u
@variable(model, 1 <= u[1:N, 1:M] <= N-1)

######################## routing constraints ######################

# (2) each client served once by one provider
@constraint(model, served[j in 2:N], sum(x[i,j,k] for i in 1:N for k in 1:M) == 1)

# (3) route continuity
@constraint(model, continuity[p in 2:N, k in 1:M], sum(x[i,p,k] for i in 1:N) == sum(x[p,j,k] for j in 1:N) )

# (4) not all providers need to leave 14HS (and start working)
@constraint(model, leave_depot[k in 1:M], sum(x[1,j,k] for j in 2:N) <= 1)

# (7) at least 1 provider must leave 14HS and start working to serve clients
@constraint(model, sum(x[1,j,k] for j in 2:N for k in 1:M) >= 1)

#(9) subtour elimination for each provider
for k in 1:M
    for i in 2:N
        for j in 2:N
            @constraint(model, u[i,k] - u[j,k] + (N-1) * x[i,j,k] <= (N-1) - 1 )
        end
    end
end

######################## end routing constraints ######################

######################## time constraints ######################

# (10) initial client service times
@constraint(model, start[j in 2:N, k in 1:M], x[1,j,k]*(s[k,1] + t[1,j] + d[1]) <= s[k,j])

# (11) in between start and end of work, service times
@constraint(model, bet[i in 2:N, j in 2:N, k in 1:M], x[i,j,k]*(s[k,i] + t[i,j] + d[i]) <= s[k,j])

# (12) come back to 14HS by hour l.
@constraint(model, endofwork[i in 2:N, k in 1:M], x[i,1,k]*(s[k,i] + t[i,1] + d[i]) <= l)

######################## end time constraints ######################

# newer objective to minimise total cost based on hours of work.
### Note. the max() function inside the objective can only be solved by Gurobi solver (Cbc cannot solve)
#@objective(model, Min, sum(x[1,j,k]*f[k]*(sum(x[p,q,k]*(t[p,q] + d[p] + max(0, s[k,q] - s[k,p] - t[p,q] - d[p]) ) for p in 1:N for q in 1:N)) for j in 2:N for k in 1:M))

# original non-hourly objective for part 1
@objective(model, Min, sum(f[k]*x[1,j,k] for j in 1:N for k in 1:M))

optimize!(model)

@btime optimize!(model)

x[1,:,1]
# see which providers are chosen to provide service.
for k in 1:M
    if sum(value.(x[1,:,k])) == 1
        println("Provider $k provides service.")
        println("Provider $k costs $(f[k]) per hour.")
        # come back hr
        for i in 2:N
            if value.(x[i,1,k]) == 1
                println("Provider $k comes back to 14HS at $(t[i,1] + s[k,i] + d[i])hr. \n")
                break # added
            end
        end
    end
end

# Given N: number of clients+1, and provider number k, find the path
### So if you found that providers 2, 3, 4 and 5 are chosen to work, run get_path(N, 2), get_path(N, 3) etc.
function get_path(N, k)
    path = []
    for j in 2:N
        if value.(x[1,j,k]) == 1
            println("\nstart of path")
            println(1, " => ")
            path = vcat(path, [1])
            println(j, " => ")
            path = vcat(path, [j])
            i = j
            while value.(x[i,1,k]) == 0
                for m in 2:N
                    if value.(x[i,m,k]) == 1
                        println(m, " => ")
                        path = vcat(path, [m])
                        i = m
                        break
                    end
                end
            end
            println(1)
            path = vcat(path, [1])
            println("end of path.... or in client index below (not location)")
        end
    end

    for j in 2:N
        if value.(x[1,j,k]) == 1
            println("\nstart of path")
            println(0, " => ")
            println(j-1, " => ")
            i = j
            while value.(x[i,1,k]) == 0
                for m in 2:N
                    if value.(x[i,m,k]) == 1
                        println(m-1, " => ")
                        i = m
                        break
                    end
                end
            end
            println(0)
            println("end of path (in client index)")
        end
    end
    return path
end

k = 8
path = get_path(N, k)

function overview(path, k)
    println("$(s[k,1])hr: Becomes available, leaves to location $(path[2]).")
    for i in 2:length(path)-1
        println("$(s[k,path[i-1]]+t[path[i-1], path[i]] + d[path[i-1]])hr: Arrived at location $(path[i]).")
        if s[k,path[i-1]]+t[path[i-1], path[i]] + d[path[i-1]] != s[k,path[i]]
            println("   Arrived early... Wait until service start time.")
        end
        println("$(s[k,path[i]])hr: Serve for $(d[path[i]]) hrs.")
        println("$(s[k,path[i]]+d[path[i]])hr: Finished service. Leave to location $(path[i+1]).")
    end
    println("$(s[k,path[length(path)-1]] + d[path[length(path)-1]] + t[path[length(path)-1], 1])hr: Arrived at location 1 (14HS).\n")
end

# print the overview of the path for provider k, where 'path' is found above using get_path(N,k)
path = [1 7 1]
k = 2
overview(path, k)

# find the number of worked hours = provider come back time (to 14HS) -  provider start time (w_k)
worked_hr = s[k,path[length(path)-1]] + d[path[length(path)-1]] + t[path[length(path)-1], 1] - s[k,1]
println("Provider $k worked for $worked_hr hours.")


text:
export continuum_model

"Compute the mean continuum flux as a function of wavelength (x) using polynomial model (coeff)."
function continuum_model(x::T1, coeff::A2; center = zero(typeof(x)) ) where { T1<:Number, T2<:Number, A2<:AbstractVector{T2} }
    #@assert 2 <= length(coeff) <= 10
    y = coeff[end]
    for i in (length(coeff)-1):-1:1      # a range starting a degree and running backwards to 1
        y *= (x.-center)
        y += coeff[i]
    end
    return y
end

function continuum_model(x::A1, coeff::A2; center = zero(typeof(x)) ) where { T1<:Number, T2<:Number, A1<:AbstractArray{T1}, A2<:AbstractVector{T2} }
    f(y) = continuum_model(y,coeff,center=center)  # make helper function so that broadcasting over x is unambiguous
    @inbounds f.(x)
end


text:
using Plots
include_ddp()

# Model
include_model("rocket3D")

function fd(model::Rocket3D{Midpoint, FixedTime}, x, u, w, h, t)
	if t == 1
		_h = h + u[4]
		# _hh = u[4]
	else
		_h = x[13]
	end
	# h = 1.0 / (1.0 + exp(-h))
	return [view(x, 1:model.n) + _h * f(model, view(x, 1:model.n) + 0.5 * _h * f(model, view(x, 1:model.n), view(u, 1:model.m), w), view(u, 1:model.m), w);
	        _h]
end

n = [t == 1 ? model.n : model.n + 1 for t = 1:T]
m = [t == 1 ? model.m + 1 : model.m for t = 1:T]

fd(model, rand(model.n), rand(model.m + 1), zeros(model.d), 1.0, 1)
fd(model, rand(model.n + 1), rand(model.m), zeros(model.d), 1.0, 2)

# Time
T = 201
h = 0.01

# Initial conditions, controls, disturbances
x1 = zeros(model.n)
x1[1] = 1.0
x1[2] = 1.0
x1[3] = 10.0
mrp = MRP(RotY(-0.5 * π) * RotX(0.25 * π))
x1[4:6] = [mrp.x; mrp.y; mrp.z]

xT = zeros(model.n)
# xT[1] = 2.5
# xT[2] = 0.0
xT[3] = model.length

u_ref = [0.0; 0.0; 0.0]
ū = [t == 1 ? [ū_fixed_time[t]; 0.0] : ū_fixed_time[t] for t = 1:T-1]
# ū = [t == 1 ? [[1.0e-2; 1.0e-2; 1.0e-2] .* randn(model.m); 0.0] : [1.0e-2; 1.0e-2; 1.0e-2] .* randn(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# # Objective
# Q = [(t == 1 ? Diagonal([1.0e-3 * ones(3); 0.0 * ones(3); 1.0e-3 * ones(3); 10.0 * ones(3)])
#         : (t == T ? Diagonal([0.0 * ones(model.n); 0.0])
# 		: Diagonal([1.0e-3 * ones(3); 0.0 * ones(3); 1.0e-3 * ones(3); 10.0 * ones(3); 0.0]))) for t = 1:T]
# q = [(t == 1 ? -2.0 * Q[t] * xT
#  	 : -2.0 * Q[t] * [xT; 0.0]) for t = 1:T]
#
# R = [(t == 1 ? Diagonal([100.0; 100.0; 1.0; 1.0e-6])
# 	 : Diagonal([100.0; 100.0; 1.0])) for t = 1:T-1]
# r = [(t == 1 ? [-2.0 * R[t][1:model.m, 1:model.m] * u_ref; -2.0 * 1.0 * 0.0]
# 	 : -2.0 * R[t] * u_ref)  for t = 1:T-1]
#
# obj = StageCosts([QuadraticCost(Q[t], q[t],
# 	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)
#
# function g(obj::StageCosts, x, u, t)
# 	T = obj.T
#     if t < T
# 		Q = obj.cost[t].Q
# 		q = obj.cost[t].q
# 	    R = obj.cost[t].R
# 		r = obj.cost[t].r
#         return x' * Q * x + q' * x + u' * R * u + r' * u
#     elseif t == T
# 		Q = obj.cost[T].Q
# 		q = obj.cost[T].q
#         return x' * Q * x + q' * x
#     else
#         return 0.0
#     end
# end

# Objective
Q = [(t < T ? 1.0 * Diagonal([1.0e-1 * ones(3); 0.0 * ones(3); 1.0e-1 * ones(3); 1000.0 * ones(3)])
        : 0.0 * Diagonal(0.0 * ones(model.n))) for t = 1:T]
q = [-2.0 * Q[t] * xT for t = 1:T]

R = [Diagonal([10000.0; 10000.0; 100.0]) for t = 1:T-1]
r = [-2.0 * R[t] * u_ref  for t = 1:T-1]

obj = StageCosts([NonlinearCost() for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		t == 1 ? (_h = h + u[end]) : (_h = x[end])
		view(u, 1:model.m)
        return _h * (view(x, 1:model.n)' * Q[t] * view(x, 1:model.n) + q[t]' * view(x, 1:model.n) + view(u, 1:model.m)' * R[t] * view(u, 1:model.m) + r[t]' * view(u, 1:model.m))
		_h = x[end]
        return _h * (view(x, 1:model.n)' * Q[t] * view(x, 1:model.n) + q[t]' * view(x, 1:model.n))
    else
        return 0.0
    end
end

# g(obj, x̄[1], ū[1], 1)
# g(obj, x̄[2], ū[2], 2)
# g(obj, x̄[T], nothing, T)

# Constraints
p = [t < T ? (t == 1 ? 2 * model.m + 2 : 2 * model.m) : model.n for t = 1:T]
info_1 = Dict(:ul => [-5.0; -5.0; 0.0; 0.0], :uu => [5.0; 5.0; 100.0; 0.0], :inequality => (1:(2 * model.m + 2)))
info_t = Dict(:ul => [-5.0; -5.0; 0.0], :uu => [5.0; 5.0; 100.0], :inequality => (1:2 * model.m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? (t == 1 ? info_1 : info_t) : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		c .= [ul - u; u - uu]
	elseif t == T
		xT = cons.con[T].info[:xT]
		c .= view(x, 1:model.n) - xT
	else
		nothing
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T,
	n = n, m = m)

prob.m_data.obj.ρ[1][4] *= 1000.0
prob.m_data.obj.ρ[1][4]

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	con_tol = 1.0e-3,
	ρ_init = 1.0, ρ_scale = 10.0)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# Time step
@show u[1][end]
x[2][end]

# Trajectories
plot(hcat([ut[1:model.m] for ut in u]...)', linetype = :steppost)
plot(hcat([xt[1:model.n] for xt in x]...)[1:3, :]', linetype = :steppost)

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
# open(vis)
visualize!(vis, model, x̄, Δt = ū[1][4])

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using UserNSSandbox_jll

# Use User Namespaces to provide isolation on Linux hosts whose kernels support it
export UserNamespacesExecutor, UnprivilegedUserNamespacesExecutor, PrivilegedUserNamespacesExecutor

abstract type UserNamespacesExecutor <: SandboxExecutor; end

# Because we can run in "privileged" or "unprivileged" mode, let's treat
# these as two separate, but very similar, executors.
struct UnprivilegedUserNamespacesExecutor <: UserNamespacesExecutor
end
struct PrivilegedUserNamespacesExecutor <: UserNamespacesExecutor
end

Base.show(io::IO, exe::UnprivilegedUserNamespacesExecutor) = write(io, "Unprivileged User Namespaces Executor")
Base.show(io::IO, exe::PrivilegedUserNamespacesExecutor) = write(io, "Privileged User Namespaces Executor")

function executor_available(::Type{T}; verbose::Bool=false) where {T <: UserNamespacesExecutor}
    return check_kernel_version(;verbose) && probe_executor(T(); test_read_only_map=true, test_read_write_map=true, verbose)
end

function check_kernel_version(;verbose::Bool = false)
    # Don't bother to do anything on non-Linux
    if !Sys.islinux()
        return true
    end
    kernel_version = get_kernel_version()

    # If we were unable to parse any part of the version number, then warn and exit.
    if kernel_version === nothing
        @warn("Unable to check version number; assuming kernel version >= 3.18")
        return true
    end

    # Otherwise, we have a kernel version and if it's too old, we should freak out.
    if kernel_version < v"3.18"
        if verbose
            @warn("Kernel version too old: detected $(kernel_version), need at least 3.18!")
        end
        return false
    end

    if verbose
        @info("Parsed kernel version \"$(kernel_version)\"")
    end
    return true
end

function build_executor_command(exe::UserNamespacesExecutor, config::SandboxConfig, user_cmd::Cmd)
    # Check to make sure that
    # While we would usually prefer to use the `executable_product()` function to get a
    # `Cmd` object that has all of the `PATH` and `LD_LIBRARY_PATH` environment variables
    # set properly so that the executable product can be run, we are careful to ensure
    # that `sandbox` has no dependencies (as much as that is possible).
    cmd_string = String[UserNSSandbox_jll.sandbox_path]

    # Enable verbose mode on the sandbox wrapper itself
    if config.verbose
        push!(cmd_string, "--verbose")
    end

    # Extract the rootfs, as it's treated specially
    append!(cmd_string, ["--rootfs", config.read_only_maps["/"]])

    # Add our `--cd` command
    append!(cmd_string, ["--cd", config.pwd])

    # Add in read-only mappings (skipping the rootfs)
    for (dst, src) in config.read_only_maps
        if dst == "/"
            continue
        end
        append!(cmd_string, ["--map", "$(src):$(dst)"])
    end

    # Add in read-write mappings
    for (dst, src) in config.read_write_maps
        append!(cmd_string, ["--workspace", "$(src):$(dst)"])
    end

    # If we're running in privileged mode, we need to add `sudo` (or `su`, if `sudo` doesn't exist)
    if isa(exe, PrivilegedUserNamespacesExecutor)
        # Next, prefer `sudo`, but allow fallback to `su`. Also, force-set our
        # environmental mappings with sudo, because many of these are often  lost
        # and forgotten due to `sudo` restrictions on setting `LD_LIBRARY_PATH`, etc...
        if sudo_cmd()[1] == "sudo"
            sudo_envs = vcat([["-E", "$k=$(config.env[k])"] for k in keys(config.env)]...)
            if user_cmd.env !== nothing
                append!(sudo_envs, vcat([["-E", pair] for pair in user_cmd.env]...))
            end
            prepend!(cmd_string, String[sudo_cmd()..., sudo_envs...])
        else
            prepend!(sudo_cmd, sudo_cmd())
        end
    end

    # Finally, append the user-requested command string
    push!(cmd_string, "--")
    append!(cmd_string, user_cmd.exec)

    # Construct a `Cmd` object off of those, with the SandboxConfig's env (if this is an unprivileged runner):
    sandbox_cmd = Cmd(cmd_string)
    if isa(exe, UnprivilegedUserNamespacesExecutor)
        sandbox_cmd = setenv(sandbox_cmd, config.env)

        # If the user has provided an environment with their command, merge that in as well
        if user_cmd.env !== nothing
            sandbox_cmd = addenv(sandbox_cmd, user_cmd.env)
        end
    end

    # If the user has asked that this command be allowed to fail silently, pass that on
    if user_cmd.ignorestatus
        sandbox_cmd = ignorestatus(sandbox_cmd)
    end

    return sandbox_cmd
end


text:
function initlearning(N,T,J,Choice,y,x)

    bigN    = zeros(Int64,J)
    bstart  = zeros(size(x,2),J)
    resid   = Array{Float64}[]
    abilsub = Array{Float64}[]
    Resid   = zeros(size(y))
    Csum    = zeros(N*S,J)
    tresid  = zeros(N*S,J)
    for j=1:J
        bigN[j]                  = sum(Choice.==j) # later sum needs to be weighted by PTypesl
        flag                     = y[:,j].!=999
        bstart[:,j]              = x[flag,:,j]\y[flag,j]
        push!(resid  ,y[flag,j].-x[flag,:,j]*bstart[:,j])
        push!(abilsub,y[flag,j].-x[flag,:,j]*bstart[:,j])
        Resid[vec(Choice.==j),j] = resid[j]
        Csum[:,j]                = (sum(reshape(Choice.==j,(T,N*S));dims=1))'
        tresid[:,j]              = (sum(reshape(Resid[:,j],(T,N*S));dims=1))'
    end

    abil = tresid./(Csum.+eps())
    abil = kron(abil,ones(T,1))

    Psi1 = deepcopy(Csum)

    for j=1:J
        resid[j]               .+= abil[vec(Choice.==j),j]
        Resid[vec(Choice.==j),j] = resid[j]
    end

    sig2=zeros(500,J)
    cov2=zeros(500,convert(Int64,(J+1)*J/2))

    sigtemp = zeros(J,1)
    covtemp = zeros(J)

    return bigN,bstart,resid,abilsub,Resid,Csum,tresid,Psi1,sig2,cov2,sigtemp,covtemp
end


text:
function SimpleExponentialAtmosphere(arg0::BodyShape, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return SimpleExponentialAtmosphere((BodyShape, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function get_density(obj::SimpleExponentialAtmosphere, arg0::AbsoluteDate, arg1::Vector3D, arg2::Frame)
    return jcall(obj, "getDensity", jdouble, (AbsoluteDate, Vector3D, Frame), arg0, arg1, arg2)
end

function get_density(obj::SimpleExponentialAtmosphere, arg0::FieldAbsoluteDate, arg1::FieldVector3D, arg2::Frame)
    return jcall(obj, "getDensity", RealFieldElement, (FieldAbsoluteDate, FieldVector3D, Frame), arg0, arg1, arg2)
end

function get_frame(obj::SimpleExponentialAtmosphere)
    return jcall(obj, "getFrame", Frame, ())
end


text:
##==============================================================================
## Accessors
##==============================================================================

##==============================================================================
## GenericFunctionNodeData
##==============================================================================

## COMMON
# getSolveInProgress
# isSolveInProgress

#TODO  getFactorFunction = getFactorType
"""
    $SIGNATURES

Return reference to the user factor in `<:AbstractDFG` identified by `::Symbol`.
"""
getFactorFunction(fcd::GenericFunctionNodeData) = fcd.fnc.usrfnc!
getFactorFunction(fc::DFGFactor) = getFactorFunction(getSolverData(fc))
function getFactorFunction(dfg::G, fsym::Symbol) where G <: AbstractDFG
  getFactorFunction(getFactor(dfg, fsym))
end

"""
    $SIGNATURES

Return user factor type from factor graph identified by label `::Symbol`.

Notes
- Replaces older `getfnctype`.
"""
getFactorType(data::GenericFunctionNodeData) = data.fnc.usrfnc!
getFactorType(fct::DFGFactor) = getFactorType(getSolverData(fct))
function getFactorType(dfg::G, lbl::Symbol) where G <: AbstractDFG
  getFactorType(getFactor(dfg, lbl))
end

##==============================================================================
## Factors
##==============================================================================
# |                   | label | tags | timestamp | solvable | solverData |
# |-------------------|:-----:|:----:|:---------:|:--------:|:----------:|
# | SkeletonDFGFactor |   X   |   x  |           |          |            |
# | DFGFactorSummary  |   X   |   X  |     X     |          |            |
# | DFGFactor         |   X   |   X  |     X     |     X    |      X     |

##------------------------------------------------------------------------------
## label
##------------------------------------------------------------------------------

## COMMON
# getLabel

##------------------------------------------------------------------------------
## tags
##------------------------------------------------------------------------------

## COMMON
# getTags
# setTags!

##------------------------------------------------------------------------------
## timestamp
##------------------------------------------------------------------------------

## COMMON
# getTimestamp

function setTimestamp(f::DFGFactor, ts::DateTime)
    return DFGFactor(f.label, ts, f.nstime, f.tags, f.solverData, f.solvable, getfield(f,:_variableOrderSymbols))
end

function setTimestamp(f::DFGFactorSummary, ts::DateTime)
    return DFGFactorSummary(f.label, ts, f.tags, f._variableOrderSymbols)
end

##------------------------------------------------------------------------------
## solvable
##------------------------------------------------------------------------------

## COMMON
# getSolvable
# setSolvable!
# isSolvable

##------------------------------------------------------------------------------
## solvable
##------------------------------------------------------------------------------

## COMMON

##------------------------------------------------------------------------------
## _variableOrderSymbols
##------------------------------------------------------------------------------

#TODO perhaps making _variableOrderSymbols imutable (NTuple) will be a save option
"""
$SIGNATURES

Get the variable ordering for this factor.
Should be equivalent to getNeighbors unless something was deleted in the graph.
"""
getVariableOrder(fct::DFGFactor)::Vector{Symbol} = fct._variableOrderSymbols
getVariableOrder(dfg::AbstractDFG, fct::Symbol) = getVariableOrder(getFactor(dfg, fct))

##------------------------------------------------------------------------------
## solverData
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Retrieve solver data structure stored in a factor.
"""
function getSolverData(f::F) where F <: DFGFactor
  return f.solverData
end

setSolverData!(f::DFGFactor, data::GenericFunctionNodeData) = f.solverData = data

##------------------------------------------------------------------------------
## utility
##------------------------------------------------------------------------------

"""
    $SIGNATURES

Return `::Bool` on whether given factor `fc::Symbol` is a prior in factor graph `dfg`.
"""
function isPrior(dfg::G, fc::Symbol)::Bool where G <: AbstractDFG
  fco = getFactor(dfg, fc)
  getFactorType(fco) isa FunctorSingleton
end

##==============================================================================
## Layer 2 CRUD (none) and Sets
##==============================================================================

##==============================================================================
## TAGS - See CommonAccessors
##==============================================================================


text:
using YAML, DataFrames

const _dbFolder = joinpath(Pkg.dir("GlassCatalog"), "assets","refractiveindex.info-database","database")
const _libFile = joinpath(_dbFolder, "library.yml")

# Prepare for Julia 0.5 unified strings by overriding deprecated Base.String
typealias String UTF8String

strip_html(s::AbstractString) = replace(s, r"<[^>]*>", "")

_library = DataFrame(Name=String[], Alias=String[], Category=String[], CatAlias=String[], Shelf=String[], Path=String[])
_shelf_info = Dict{String,String}()

#= Form of the YAML library.yml file
 - SHELF: main
  name: "MAIN - simple inorganic materials"
  content:
    - DIVIDER: "Ag - Silver"
    - BOOK: Ag
      name: "Ag (Silver)"
      content:
        - DIVIDER: "optional sub-grouping of pages"
        - PAGE: Rakic
          name: "Rakić 1998: n,k 0.2066-12.40 µm"
          path: "main/Ag/Rakic.yml"
        - PAGE: McPeak
          name: "McPeak et al. 2015: Thin film; n,k 0.3-1.7 µm"
          path: "main/Ag/McPeak.yml"
        ...
=#

# returns an array of strings representing a row in the dataframe
function _parse_page(shelf, shelfdiv, book, bookdiv, bookname, page, pagename, path)
  shelf = lowercase(shelf)
  local spath = split(path, '/')
  local lc_shelfdiv = lowercase(shelfdiv)
  if shelf == "main"
    # shelfdiv form: "Al - Aluminium and aluminates"
    # parse as: "$div_formula - $div_name"
    m = match(r"^([^\-]+)\s\-\s(.+)$", shelfdiv)
    @assert m !== nothing && length(m.captures) >= 2 "Bad shelfdiv in $shelf: $shelfdiv"
    (div_formula, div_name) = m.captures

    # book form: "MgAl2O4"
    # bookname form: "MgAl2O4 (Magnesium aluminate, spinel)"
    # parse as: "$bk_formula ($bk_name)"
    m = match(r"^(.+)\s\((.+)\)$", bookname)
    @assert m !== nothing && length(m.captures) >= 2 "Bad bookname in $shelf: $bookname"
    (bk_formula, bk_name) = m.captures

    # page form: "Tropf"
    # pagename form: "Tropf and Thomas 1991: n 0.35-5.5 µm"
    # leave unparsed

    # path form: "glass/schott/FK5HTi.yml"
    # parse as: ".../$(file_name).yml"
    @assert length(spath) > 0
    file_name = split(spath[end], '.')[1]

    entry_name = isempty(bookdiv) ? "$book-$file_name" : "$bookdiv-$file_name"

    # return (name, alias, category, catalias, shelf, path)
    return (entry_name, "$bk_name - $pagename", div_formula, div_name, shelf, path)

  elseif shelf == "organic"
    # shelfdiv form: "Benzene and its derivatives"
    # leave unparsed

    # book form: "benzene"
    # bookname form: "C6H6 (Benzene)"
    # parse as: "$bk_formula ($bk_name)"
    m = match(r"^(.+)\s\((.+)\)$", bookname)
    @assert m !== nothing && length(m.captures) >= 2 "Bad bookname in $shelf: $bookname"
    (bk_formula, bk_name) = m.captures

    # page form: "Moutzouris"
    # pagename form: "Moutzouris et al. 2013: n 0.450-1.551 µm"
    # leave unparsed

    # return (name, alias, category, catalias, shelf, path)
    return ("$book-$page", "$bk_name - $pagename", shelfdiv, shelfdiv, shelf, path)

  elseif shelf == "glass" && startswith(lc_shelfdiv, "optical")
    # shelfdiv form: "Optical glass - SCHOTT"
    # leave unparsed

    # book form: "SCHOTT-FK"
    # bookname form: "SCHOTT - FK (Fluorite crown)"
    #   OR "SCHOTT - Inquiry glasses"
    # parse as: "$bk_mfg - $bk_class ($bk_classname)"
    #   OR "$bk_mfg - $bk_classname"
    if ismatch(r"^(.+)\s-\s(.+)\s\((.+)\)$", bookname)
      m = match(r"^(.+)\s-\s(.+)\s\((.+)\)$", bookname)
      @assert m !== nothing && length(m.captures) >= 3 "Bad bookname in $shelf: $bookname"
      (bk_mfg, bk_class, bk_classname) = m.captures
    else
      m = match(r"^(.+)\s-\s(.+)$", bookname)
      @assert m !== nothing && length(m.captures) >= 2 "Bad bookname in $shelf: $bookname"
      (bk_mfg, bk_classname) = m.captures
      m = nothing #match(r"^(.+) glasses$", bk_classname)
      if m === nothing
        bk_class = bk_classname
      else
        (bk_class) = m.captures
      end
    end

    # page form: "FK5HTi"
    # pagename form: "FK5HTi"
    # leave unparsed

    # return (name, alias, category, catalias, shelf, path)
    return ("$page-$bk_mfg", "$pagename - $bk_mfg", "$bk_mfg - $bk_class", "$shelfdiv - $bk_classname", shelf, path)

  elseif shelf == "glass" && startswith(lc_shelfdiv, "popular")
    # popular glasses are listed twice. Modify the entry name to start with popular
    # return (name, alias, category, catalias, shelf, path)
    return ("popular-$book-$page", "$bookname - $pagename (popular)", shelfdiv, shelfdiv, shelf, path)

  elseif shelf == "glass"
          # shelfdiv form: "Multi purpose glass"
          # leave unparsed

          # book form: "soda-lime"
          # bookname form: "Soda lime glass"
          #   OR "SCHOTT - multiple purpose"
          # parse as: "$bk_mfg - $bk_class ($bk_classname)"
          #   OR "$bk_mfg - $bk_classname"
          if ismatch(r"^(.+)\s-\s(.+)$", bookname)
            m = match(r"^(.+)\s-\s(.+)$", bookname)
            @assert m !== nothing && length(m.captures) >= 2 "Bad bookname in $shelf: $bookname"
            (bk_mfg, bk_class) = m.captures
            bk_classname = bk_class
          else
            bk_mfg = ""
            bk_class = bookname
            bk_classname = bk_class
          end

          # page form: "Rubin-clear"
          #   OR "ZERODUR"
          # leave unparsed

          if length(bk_mfg) == 0
            # return (name, alias, category, catalias, shelf, path)
            return ("$book-$page", "$bookname - $pagename", shelfdiv, shelfdiv, shelf, path)
          else
            # return (name, alias, category, catalias, shelf, path)
            return ("$page-$bk_mfg", "$pagename - $bk_mfg", "$bk_mfg - $bk_class", "$shelfdiv - $bk_classname", shelf, path)
          end

  else
    if ismatch(r"^(.+)\s\((.+)\)$", bookname)
      m = match(r"^(.+)\s\((.+)\)$", bookname)
      @assert m !== nothing && length(m.captures) >= 2 "Bad bookname in $shelf: $bookname"
      (bk_formula, bk_name) = m.captures

      # return (name, alias, category, catalias, shelf, path)
      return ("$bk_formula-$page", "$bk_name - $pagename", shelfdiv, shelfdiv, shelf, path)
    else
      if isempty(shelfdiv)
        shelfdiv = _shelf_info[shelf]
      end
      # return (name, alias, category, catalias, shelf, path)
      return ("$book-$page", "$bookname - $pagename", shelfdiv, shelfdiv, shelf, path)
    end
  end
end

function _parse_library()
  # load the .yml file as a Dict hierarchy
  lib_d = YAML.load_file(_libFile)

  for shelf_d in lib_d
    (shelf, shelfname) = (shelf_d["SHELF"], shelf_d["name"])
    _shelf_info[shelf] = shelfname
    #if lowercase(shelf) != "main"
    #  continue
    #end
    shelfdiv::String = ""
    for book_d in shelf_d["content"]
      if haskey(book_d,"DIVIDER")
        shelfdiv = book_d["DIVIDER"]
      else
        (book, bookname) = (book_d["BOOK"], strip_html(book_d["name"]))
        bookdiv::String = ""
        for page_d in book_d["content"]
          if haskey(page_d,"DIVIDER")
            bookdiv = page_d["DIVIDER"]
          end
          if haskey(page_d, "PAGE")
            (page, pagename, path) = (page_d["PAGE"], strip_html(page_d["name"]), page_d["path"])
            entry = _parse_page(shelf, shelfdiv, book, bookdiv, bookname, page, pagename, path)
            push!(_library, entry)
          end
        end
      end
    end
  end
end

using DataFrames

function as_dataframe()
  _parse_library()
  return _library
end

writetable(joinpath(homedir(),"Devel/sandbox/julia/df.csv"), as_dataframe())


text:
using .Turing, .MCMCChains
export pointwise_log_likelihoods

const TURING_LOOP_WARN = """
**Important Note:** The posterior log-likelihood must be computed with a `for` loop inside a
Turing model; broadcasting will result in all observations being treated as if they are a
single point.
"""
const CHAINS_ARG = """
`chains::Chains`: A chain object from MCMCChains.
"""
const TURING_MODEL_ARG = """
`model`: A Turing model with data in the form of `model(data)`.
"""

"""
    pointwise_log_likelihoods(model::DynamicPPL.Model, chain::Chains)

Compute the pointwise log-likelihoods from a Turing model.

$TURING_LOOP_WARN

# Arguments
  - $TURING_MODEL_ARG
  - $CHAINS_ARG

# Returns
  - `Array`: A three dimensional array of pointwise log likelihoods. This array should be
    indexed using `array[data, sample, chain]`.
"""
function pointwise_log_likelihoods(model::DynamicPPL.Model, chain::Chains)

    @info TURING_LOOP_WARN

    # subset of chain for mcmc samples
    chain_params = MCMCChains.get_sections(chain, :parameters)
    # compute the pointwise log likelihoods
    log_like_dict = DynamicPPL.pointwise_loglikelihoods(model, chain_params)
    # Size of array (n_steps, n_chains) using first parameter
    dims = size(last(first(log_like_dict)))
    # parse "var[i]" -> i
    ind_from_string(x) = parse(Int, split(split(x, "[")[2], "]")[1])
    # collect variable names
    sorted_keys = sort(collect(keys(log_like_dict)); by=ind_from_string)
    # Convert from dictionary to 3d array
    array = [reshape(log_like_dict[i], 1, dims...) for i in sorted_keys]
    return reduce(vcat, array)
end

"""
    function psis_loo(
        ll_fun::Function,
        chains::Chains
        [, args...; kwargs...]
    ) -> PsisLoo

Use Pareto-Smoothed Importance Sampling to calculate the leave-one-out cross validation
score from an MCMCChain object and a Turing model.

$TURING_LOOP_WARN

# Arguments

  - $CHAINS_ARG
  - $TURING_MODEL_ARG
  - $ARGS [`psis`](@ref).
  - $KWARGS [`psis`](@ref).

See also: [`psis`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).
"""
function psis_loo(model::DynamicPPL.Model, chain::Chains, args...; kwargs...)
    @info TURING_LOOP_WARN
    pointwise_log_likes = pointwise_log_likelihoods(model, chain)
    return psis_loo(pointwise_log_likes, args...; kwargs...)
end

"""
    psis(
        model::DynamicPPL.Model,
        chains::Chains,
        args...;
        kwargs...
    ) -> Psis

Generate samples using Pareto smoothed importance sampling (PSIS).

$TURING_LOOP_WARN

# Arguments
  - $TURING_MODEL_ARG
  - $CHAINS_ARG
  - $ARGS [`psis`](@ref).
  - $KWARGS [`psis`](@ref).

See also: [`psis`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).
"""
function psis(model::DynamicPPL.Model, chain::Chains, args...; kwargs...)
    @info TURING_LOOP_WARN
    log_ratios = pointwise_log_likelihoods(model, chain)
    return psis(-log_ratios, args...; kwargs...)
end