text:
"""
   `interp_chunk_to_shifted_grid_sinc!( flux_out, var_out, chunk_of_spectrum, wavelengths, boost_factor )`
Return spectra interpolated onto a grid of points using sinc interpolation.
# Arguments:
- flux_out: (results stored into this array)
- var_out: (results stored into this array)
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
- boost_factor: divide wavelengths by boost_factor
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_shifted_grid_sinc!( flux_out::AbstractArray{T1,1}, var_out::AbstractArray{T2,1}, chunk::AC, grid::AR, boost_factor::Real; Filter::AbstractVector = zeros(0) ) where {
    T1<:Real, T2<:Real, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    if issorted(chunk.λ)
        use_λ = chunk.λ
        use_flux = chunk.flux
        use_var = chunk.var
    else  # Deal with ocassion issue near boundary of pixels with LFC and non-LFC calibration, where pixel wavelengths may not be monotonic.
        mask = trues(length(chunk.λ))
        λ_last = first(chunk.λ)
        for i in 2:length(mask)
            if chunk.λ[i]<=λ_last
                mask[i] = false
            else
                λ_last = chunk.λ[i]
            end
        end
        use_λ = chunk.λ[mask]
        use_flux = chunk.flux[mask]
        use_var = chunk.var[mask]
        @assert issorted(use_λ)
    end

    if length(Filter) >= 1
        flux_out .= SincInterpolation.spectra_interpolate(grid./ boost_factor, use_λ,use_flux, Filter=Filter)
        var_out .= SincInterpolation.spectra_interpolate(grid./ boost_factor, use_λ ,use_var, Filter=Filter)
    else
        flux_out .= SincInterpolation.spectra_interpolate(grid./ boost_factor, use_λ ,use_flux)
        var_out .= SincInterpolation.spectra_interpolate(grid ./ boost_factor, use_λ,use_var)
    end
    return flux_out
end

"""
   `interp_chunk_to_shifted_grid_sinc( chunk_of_spectrum, wavelengths, boost_factor )`
Return spectra interpolated onto a grid of points using sinc interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
- boost_factor: divide wavelengths by boost_factor
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_shifted_grid_sinc( chunk::AC, grid::AR, boost_factor::Real; Filter::AbstractVector = zeros(0) ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    if length(Filter) >= 1
        interp_chunk_to_shifted_grid_sinc!(flux_out, var_out, chunk, grid, boost_factor, Filter=Filter)
    else
        interp_chunk_to_shifted_grid_sinc!(flux_out, var_out, chunk, grid, boost_factor)
    end
    return (flux=flux_out, var=var_out)
end

"""
   `interp_chunk_to_grid_sinc!( flux_out, var_out, chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using sinc interpolation.
# Arguments:
- flux_out: (results stored into this array)
- var_out: (results stored into this array)
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
Optional Arguments:
- Filter: Vector with pre-allocated workspace (if length>=1)
# Returns
- flux_out
"""
function interp_chunk_to_grid_sinc!( flux_out::AbstractArray{T1,1}, var_out::AbstractArray{T2,1}, chunk::AC, grid::AR; Filter::AbstractVector = zeros(0) ) where {
    T1<:Real, T2<:Real, AC<:AbstractChunkOfSpectrum, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    @assert size(flux_out) == size(var_out)
    @assert size(flux_out) == size(grid)
    if length(Filter) >= 1
        flux_out .= SincInterpolation.spectra_interpolate(grid,chunk.λ,chunk.flux,Filter=Filter)
        var_out .= SincInterpolation.spectra_interpolate(grid,chunk.λ,chunk.var,Filter=Filter)
    else
        flux_out .= SincInterpolation.spectra_interpolate(grid,chunk.λ,chunk.flux)
        var_out .= SincInterpolation.spectra_interpolate(grid,chunk.λ,chunk.var)
    end
    return flux_out
end

"""
   `interp_chunk_to_grid_sinc( chunk_of_spectrum, wavelengths )`
Return spectra interpolated onto a grid of points using sinc interpolation.
# Arguments:
- chunk_of_spectrum
- wavelengths: AbstractRange or AbstractArray of locations where chunk is to be interpolated to
# Returns
- flux_out
"""
function interp_chunk_to_grid_sinc( chunk::AC, grid::AR ) where {  AC<:AbstractChunkOfSpectrum, T2<:Real, AR<:Union{AbstractRange,AbstractArray{T2,1}} }
    flux_out = Array{Float64,1}(undef,length(grid))
    var_out = Array{Float64,1}(undef,length(grid))
    interp_chunk_to_grid_sinc!(flux_out, var_out, chunk, grid)
    return (flux=flux_out, var=var_out)
end


text:
include("MainModule.jl")
import Gtk
function main()
  win = MainModule.MainWindow()
  Gtk.showall(win)
  if !isinteractive()
    c = Condition()
    Gtk.signal_connect(win, :destroy) do widget
      notify(c)
    end
    wait(c)
  end
end

ext: jl

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui


text:
module Mechanic

using Modia

Flange = Model(
    phi = potential,
    tau = flow
)

Fixed = Model(
    flange = Flange,
    equations = :[
        flange.phi = 0.0]
)

PU2Abs = Model(
    flange_pu = Flange,
    flange_SI = Flange,
    Speedrad = 1.0,
    Tbase = 1.0,
    NoPolepairs = 1,
    speedpu = Var(),
    mechpower = Var(),
    mechPower = Var(),
    equations = :[
        flange_SI.phi = Speedrad * flange_pu.phi
        0.0 = Tbase * flange_pu.tau + flange_SI.tau
        speedpu = der(flange_pu.phi)
        mechpower = flange_pu.tau * speedpu
        mechPower = mechpower * Tbase * Speedrad
    ]
)

Inertia = Model(
    flange_a = Flange,
    flange_b = Flange,
    J = 1.0,
    phi = Var(start=0.0),
    w = Var(),
    equations = :[
        phi = flange_a.phi
        phi = flange_b.phi
        w = der(phi)
        a = der(w)
        J * a = flange_a.tau + flange_b.tau ]
)

end

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

module PerfectPlasticModule

using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters

import ..AbstractMaterial, ..AbstractMaterialState
import ..Utilities: Symm2, Symm4, isotropic_elasticity_tensor, IS, ID, lame
import ..integrate_material!  # for method extension

# parametrically polymorphic for any type representing ℝ
export GenericPerfectPlastic, GenericPerfectPlasticDriverState, GenericPerfectPlasticParameterState, GenericPerfectPlasticVariableState

# specialization for Float64
export PerfectPlastic, PerfectPlasticDriverState, PerfectPlasticParameterState, PerfectPlasticVariableState

@with_kw mutable struct GenericPerfectPlasticDriverState{T <: Real} <: AbstractMaterialState
    time::T = zero(T)
    strain::Symm2{T} = zero(Symm2{T})
end

"""Parameter state for perfect plastic material.

- youngs_modulus
- poissons_ratio
- yield_stress
"""
@with_kw struct GenericPerfectPlasticParameterState{T <: Real} <: AbstractMaterialState
    youngs_modulus::T = zero(T)
    poissons_ratio::T = zero(T)
    yield_stress::T = zero(T)
end

"""Problem state for perfect plastic material.

- `stress`: stress tensor
- `plastic_strain`: plastic part of strain tensor
- `cumeq`: cumulative equivalent plastic strain (scalar, ≥ 0)
- `jacobian`: ∂σij/∂εkl
"""
@with_kw struct GenericPerfectPlasticVariableState{T <: Real} <: AbstractMaterialState
    stress::Symm2{T} = zero(Symm2{T})
    plastic_strain::Symm2{T} = zero(Symm2{T})
    cumeq::T = zero(T)
    jacobian::Symm4{T} = zero(Symm4{T})
end

@with_kw mutable struct GenericPerfectPlastic{T <: Real} <: AbstractMaterial
    drivers::GenericPerfectPlasticDriverState{T} = GenericPerfectPlasticDriverState{T}()
    ddrivers::GenericPerfectPlasticDriverState{T} = GenericPerfectPlasticDriverState{T}()
    variables::GenericPerfectPlasticVariableState{T} = GenericPerfectPlasticVariableState{T}()
    variables_new::GenericPerfectPlasticVariableState{T} = GenericPerfectPlasticVariableState{T}()
    parameters::GenericPerfectPlasticParameterState{T} = GenericPerfectPlasticParameterState{T}()
    dparameters::GenericPerfectPlasticParameterState{T} = GenericPerfectPlasticParameterState{T}()
end

PerfectPlastic = GenericPerfectPlastic{Float64}
PerfectPlasticDriverState = GenericPerfectPlasticDriverState{Float64}
PerfectPlasticParameterState = GenericPerfectPlasticParameterState{Float64}
PerfectPlasticVariableState = GenericPerfectPlasticVariableState{Float64}

"""
    integrate_material!(material::GenericPerfectPlastic)

Perfect plastic material: no hardening. The elastic region remains centered on the
origin, and retains its original size.

This is a standard basic plasticity model; see a textbook, such as:

    J. C. Simo, T. J. R. Hughes. Computational Inelasticity. Interdisciplinary
    Applied Mathematics volume 7. Springer. 1998. http://dx.doi.org/10.1007/b98904

    The notation in the book somewhat differs from ours; see:
    https://github.com/JuliaFEM/Materials.jl/pull/66#issuecomment-674786955
"""
function integrate_material!(material::GenericPerfectPlastic{T}) where T <: Real
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers

    E = p.youngs_modulus
    nu = p.poissons_ratio
    lambda, mu = lame(E, nu)
    R0 = p.yield_stress

    # @unpack strain, time = d  # not needed for this material
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, plastic_strain, cumeq, jacobian = v

    jacobian = isotropic_elasticity_tensor(lambda, mu)  # dσ/dε, i.e. ∂σij/∂εkl
    stress += dcontract(jacobian, dstrain)  # add the elastic stress increment, get the elastic trial stress
    seff_dev = dev(stress)
    f = sqrt(1.5)*norm(seff_dev) - R0  # von Mises yield function; f := J(seff_dev) - Y

    if f > 0.0
        # see Simo & Hughes ch. 3.1: radial return mapping, eq. (3.37)
        dp = 1.0/(3.0*mu) * f
        n = sqrt(1.5)*seff_dev/norm(seff_dev)  # a (tensorial) unit direction, s.t. 2/3 * (n : n) = 1

        plastic_strain += dp*n
        cumeq += dp  # cumulative equivalent plastic strain (note dp ≥ 0)

        # Perfect plastic material: the stress state cannot be outside the yield surface.
        # Project it back to the yield surface.
        stress -= dcontract(jacobian, dp*n)

        # Compute ∂σij/∂εkl, accounting for the plastic contribution.
        # EE = IS + dp/R0 * (∂σ/∂ε)_e : ((3/2) ID - n ⊗ n)
        EE = IS(T) + dp/R0 * dcontract(jacobian, 1.5*ID(T) - otimes(n,n))  # using the elastic jacobian
        ED = dcontract(inv(EE), jacobian)
        # J = ED - (ED : n) ⊗ (n : ED) / (n : ED : n)
        jacobian = ED - otimes(dcontract(ED, n), dcontract(n, ED)) / dcontract(dcontract(n, ED), n)
    end
    variables_new = GenericPerfectPlasticVariableState(stress=stress,
                                                       plastic_strain=plastic_strain,
                                                       cumeq=cumeq,
                                                       jacobian=jacobian)
    material.variables_new = variables_new
    return nothing
end

end


text:
export BLUE, cov_to_cor

"""
    cov_to_cor(cov::Array{<:Real, 2})

Convert a covariance matrix `cov` to a correlation matrix and a vector of uncertainty values.

Returns a matrix and a vector.
Throws a warning when the covariance matrix is not positive definite.

Example:

```julia
cor, unc = cov_to_cor(cov)
```
"""
function cov_to_cor(cov::Array{<:Real, 2})
    if !isposdef(cov)
        @warn "The covariance matrix $cov is not positive definite!"
    end

    cor = StatsBase.cov2cor(cov, sqrt.(diag(cov)))
    unc = sqrt.(diag(cov))

    return cor, unc
end

"""
    BLUE(model::EFTfitterModel)

Calculates the best linear unbiased estimator (BLUE) for multiple
measurements of the same observable, according to https://www.sciencedirect.com/science/article/pii/0168900288900186.

Note: Works only for an `EFTfitterModel` where all measurements have the same observable.
If this is not the case, an error is thrown.

Returns a `NamedTuple` with the fields:
* `:value`: BLUE value
* `:unc`: BLUE uncertainty
* `:weights`: Array with the weights for each measurement

Example:

```julia
blue = BLUE(model)
println(blue.value, blue.unc, blue.weights)
```
"""
function BLUE(m::EFTfitterModel)
    if !all_observables_equal(m)
        @error "The measurements in the EFTfitterModel have different observables.
        Cannot calculate the best linear unbiased estimator (BLUE)."
        return
    end

    meas = [meas.value for meas in values(m.measurements)]
    nmeas = length(meas)
    covariance = get_total_covariance(m)

    u = ones(nmeas)
    α = inv(covariance)*u/(u'*inv(covariance)*u)
    τ_BLUE = dot(α, meas)
    σ_BLUE = dot(α, covariance*α)

    return (value = τ_BLUE, unc = σ_BLUE, weights = α)
end

function all_observables_equal(model::EFTfitterModel)
    observable_functions = [m.observable.func for m in values(model.measurements)]
    all(y->y==observable_functions[1], observable_functions)
end


text:
using OrdinaryDiffEq, Sundials, DiffEqDevTools, Plots, ODEInterfaceDiffEq, ODE, LSODA
using Random
Random.seed!(123)
gr()
# 2D Linear ODE
function f(du,u,p,t)
  @inbounds for i in eachindex(u)
    du[i] = 1.01*u[i]
  end
end
function f_analytic(u₀,p,t)
  u₀*exp(1.01*t)
end
tspan = (0.0,10.0)
prob = ODEProblem(ODEFunction(f,analytic=f_analytic),rand(100,100),tspan)

abstols = 1.0 ./ 10.0 .^ (3:13)
reltols = 1.0 ./ 10.0 .^ (0:10);

setups = [Dict(:alg=>DP5())
          Dict(:alg=>ode45())
          Dict(:alg=>dopri5())
          Dict(:alg=>ARKODE(Sundials.Explicit(),etable=Sundials.DORMAND_PRINCE_7_4_5))
          Dict(:alg=>Tsit5())]
solnames = ["OrdinaryDiffEq";"ODE";"ODEInterface";"Sundials ARKODE";"OrdinaryDiffEq Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,save_everystep=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5(),:dense=>false)
          Dict(:alg=>ode45(),:dense=>false)
          Dict(:alg=>dopri5()) # dense=false by default: no nonlinear interpolation
          Dict(:alg=>ARKODE(Sundials.Explicit(),etable=Sundials.DORMAND_PRINCE_7_4_5),:dense=>false)
          Dict(:alg=>Tsit5(),:dense=>false)]
solnames = ["OrdinaryDiffEq";"ODE";"ODEInterface";"Sundials ARKODE";"OrdinaryDiffEq Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>ode45())
          Dict(:alg=>dopri5())
          Dict(:alg=>ARKODE(Sundials.Explicit(),etable=Sundials.DORMAND_PRINCE_7_4_5))
          Dict(:alg=>Tsit5())]
solnames = ["OrdinaryDiffEq";"ODE";"ODEInterface";"Sundials ARKODE";"OrdinaryDiffEq Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>BS3())
          Dict(:alg=>BS5())
          Dict(:alg=>Tsit5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Vern6())
          Dict(:alg=>TanYam7())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>DP8())
          Dict(:alg=>Vern9())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Vern7())
          Dict(:alg=>dop853())
          Dict(:alg=>ode78())
          Dict(:alg=>odex())
          Dict(:alg=>lsoda())
          Dict(:alg=>ddeabm())
          Dict(:alg=>ARKODE(Sundials.Explicit(),order=8))
          Dict(:alg=>CVODE_Adams())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;save_everystep=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          #Dict(:alg=>ode45())
          Dict(:alg=>Tsit5())]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:L2,dense_errors=true,numruns=100)
plot(wp)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Vern7())
          #Dict(:alg=>ode78())
          ]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;error_estimate=:L2,dense_errors=true,numruns=100)
plot(wp)

abstols = 1.0 ./ 10.0 .^ (3:13)
reltols = 1.0 ./ 10.0 .^ (0:10);
dts = [1,1/2,1/4,1/10,1/20,1/40,1/60,1/80,1/100,1/140,1/240]
setups = [Dict(:alg=>DP5())
          Dict(:alg=>ode45())
          Dict(:alg=>dopri5())
          Dict(:alg=>RK4(),:dts=>dts)
          Dict(:alg=>Tsit5())]
solnames = ["DifferentialEquations";"ODE";"ODEInterface";"DifferentialEquations RK4";"DifferentialEquations Tsit5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>Tsit5())
          Dict(:alg=>Vern9())
          Dict(:alg=>VCABM())
          Dict(:alg=>AitkenNeville(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))]
solnames = ["Tsit5","Vern9","VCABM","AitkenNeville","Midpoint Deuflhard","Midpoint Hairer Wanner"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointDeuflhard(min_order=1, max_order=9, init_order=9, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=false))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :romberg, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, sequence = :bulirsch, threading=true))]
solnames = ["Deuflhard","No threads","standard","Romberg","Bulirsch"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

setups = [Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=11, init_order=4, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=11, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=2, max_order=15, init_order=10, threading=true))
          Dict(:alg=>ExtrapolationMidpointHairerWanner(min_order=5, max_order=7, init_order=6, threading=true))]
solnames = ["1","2","3","4","5"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;names=solnames,
                      save_everystep=false,verbose=false,numruns=100)
plot(wp)

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


text:
# Addition
#
# commutativity
# S1 + S2 == S2 + S1
# S + C == C + S
# C1 + C2 == C2 + C1
# S + U - U == S (for sorted S)
#
# associativity
# (S1 + S2) + S3 == S1 + (S2 + S3)
# (S1 + S2) + C == S1 + (S2 + C)
function +(S::T, U::T) where {T<:GphysData}
  Ω = deepcopy(S)
  append!(Ω, U)
  merge!(Ω, purge_only=true)
  return Ω
end
+(S::SeisData, C::SeisChannel) = +(S, SeisData(C))
+(C::SeisChannel, S::SeisData) = +(S, SeisData(C))
+(C::SeisChannel, D::SeisChannel) = +(SeisData(C), SeisData(D))

# Subtraction
-(S::GphysData, i::Int)          = (U = deepcopy(S); deleteat!(U,i); return U)  # By channel #
-(S::GphysData, J::Array{Int,1}) = (U = deepcopy(S); deleteat!(U,J); return U)  # By array of channel #s

# Multiplication
# distributivity: (S1+S2)*S3) == (S1*S3 + S2*S3)
*(S::SeisData, U::SeisData) = merge(Array{SeisData,1}([S,U]))
*(S::SeisData, C::SeisChannel) = merge(S, SeisData(C))
function *(C::SeisChannel, D::SeisChannel)
  s1 = deepcopy(C)
  s2 = deepcopy(D)
  S = merge(SeisData(s1),SeisData(s2))
  return S
end

# Division will happen eventually; for S/U to be logical, we need to extract
# time ranges of data from S that are not in U. This will work like unix `diff`


text:
## TYPES AND CONSTRUCTOR

mutable struct DeterministicTunedModel{T,M<:Deterministic} <: MLJBase.Deterministic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,Vector{<:Real}}
    operation
    range
    selection_heuristic
    train_best::Bool
    repeats::Int
    n::Union{Int,Nothing}
    acceleration::AbstractResource
    acceleration_resampling::AbstractResource
    check_measure::Bool
    cache::Bool
end

mutable struct ProbabilisticTunedModel{T,M<:Probabilistic} <: MLJBase.Probabilistic
    model::M
    tuning::T  # tuning strategy
    resampling # resampling strategy
    measure
    weights::Union{Nothing,AbstractVector{<:Real}}
    operation
    range
    selection_heuristic
    train_best::Bool
    repeats::Int
    n::Union{Int,Nothing}
    acceleration::AbstractResource
    acceleration_resampling::AbstractResource
    check_measure::Bool
    cache::Bool
end

const EitherTunedModel{T,M} =
    Union{DeterministicTunedModel{T,M},ProbabilisticTunedModel{T,M}}

#todo update:
"""
    tuned_model = TunedModel(; model=nothing,
                             tuning=Grid(),
                             resampling=Holdout(),
                             measure=nothing,
                             weights=nothing,
                             repeats=1,
                             operation=predict,
                             range=nothing,
                             selection_heuristic=NaiveSelection(),
                             n=default_n(tuning, range),
                             train_best=true,
                             acceleration=default_resource(),
                             acceleration_resampling=CPU1(),
                             check_measure=true,
                             cache=true)

Construct a model wrapper for hyperparameter optimization of a
supervised learner.

Calling `fit!(mach)` on a machine `mach=machine(tuned_model, X, y)` or
`mach=machine(tuned_model, X, y, w)` will:

- Instigate a search, over clones of `model`, with the hyperparameter
  mutations specified by `range`, for a model optimizing the specified
  `measure`, using performance evaluations carried out using the
  specified `tuning` strategy and `resampling` strategy.

- Fit an internal machine, based on the optimal model
  `fitted_params(mach).best_model`, wrapping the optimal `model`
  object in *all* the provided data `X`, `y`(, `w`). Calling
  `predict(mach, Xnew)` then returns predictions on `Xnew` of this
  internal machine. The final train can be supressed by setting
  `train_best=false`.

The `range` objects supported depend on the `tuning` strategy
specified. Query the `strategy` docstring for details. To optimize
over an explicit list `v` of models of the same type, use
`strategy=Explicit()` and specify `model=v[1]` and `range=v`.

The number of models searched is specified by `n`. If unspecified,
then `MLJTuning.default_n(tuning, range)` is used. When `n` is
increased and `fit!(mach)` called again, the old search history is
re-instated and the search continues where it left off.

If `measure` supports weights (`supports_weights(measure) == true`)
then any `weights` specified will be passed to the measure. If more
than one `measure` is specified, then only the first is optimized
(unless `strategy` is multi-objective) but the performance against
every measure specified will be computed and reported in
`report(mach).best_performance` and other relevant attributes of the
generated report.

Specify `repeats > 1` for repeated resampling per model
evaluation. See [`evaluate!`](@ref) options for details.

*Important.* If a custom `measure` is used, and the measure is
a score, rather than a loss, be sure to check that
`MLJ.orientation(measure) == :score` to ensure maximization of the
measure, rather than minimization. Override an incorrect value with
`MLJ.orientation(::typeof(measure)) = :score`.

In the case of two-parameter tuning, a Plots.jl plot of performance
estimates is returned by `plot(mach)` or `heatmap(mach)`.

Once a tuning machine `mach` has bee trained as above, then
`fitted_params(mach)` has these keys/values:

key                 | value
--------------------|--------------------------------------------------
`best_model`        | optimal model instance
`best_fitted_params`| learned parameters of the optimal model

The named tuple `report(mach)` includes these keys/values:

key                 | value
--------------------|--------------------------------------------------
`best_model`        | optimal model instance
`best_history_entry`| corresponding entry in the history, including performance estimate
`best_report`       | report generated by fitting the optimal model to all data
`history`           | tuning strategy-specific history of all evaluations

plus other key/value pairs specific to the `tuning` strategy.

### Summary of key-word arguments

- `model`: `Supervised` model prototype that is cloned and mutated to
  generate models for evaluation

- `tuning=Grid()`: tuning strategy to be applied (eg, `RandomSearch()`)

- `resampling=Holdout()`: resampling strategy (eg, `Holdout()`, `CV()`),
  `StratifiedCV()`) to be applied in performance evaluations

- `measure`: measure or measures to be applied in performance
  evaluations; only the first used in optimization (unless the
  strategy is multi-objective) but all reported to the history

- `weights`: sample weights to be passed the measure(s) in performance
  evaluations, if supported.

- `repeats=1`: for generating train/test sets multiple times in
  resampling; see [`evaluate!`](@ref) for details

- `operation=predict`: operation to be applied to each fitted model;
  usually `predict` but `predict_mean`, `predict_median` or
  `predict_mode` can be used for `Probabilistic` models, if
  the specified measures are `Deterministic`

- `range`: range object; tuning strategy documentation describes
  supported types

- `selection_heuristic`: the rule determining how the best model is
  decided. According to the default heuristic,
  `NaiveSelection()`, `measure` (or the first
  element of `measure`) is evaluated for each resample and these
  per-fold measurements are aggregrated. The model with the lowest
  (resp. highest) aggregate is chosen if the measure is a `:loss`
  (resp. a `:score`).

- `n`: number of iterations (ie, models to be evaluated); set by
  tuning strategy if left unspecified

- `train_best=true`: whether to train the optimal model

- `acceleration=default_resource()`: mode of parallelization for
  tuning strategies that support this

- `acceleration_resampling=CPU1()`: mode of parallelization for
  resampling

- `check_measure=true`: whether to check `measure` is compatible with the
  specified `model` and `operation`)

- `cache=true`: whether to cache model-specific representations of
  user-suplied data; set to `false` to conserve memory. Speed gains
  likely limited to the case `resampling isa Holdout`.

"""
function TunedModel(; model=nothing,
                    tuning=Grid(),
                    resampling=MLJBase.Holdout(),
                    measures=nothing,
                    measure=measures,
                    weights=nothing,
                    operation=predict,
                    ranges=nothing,
                    range=ranges,
                    selection_heuristic=NaiveSelection(),
                    train_best=true,
                    repeats=1,
                    n=nothing,
                    acceleration=default_resource(),
                    acceleration_resampling=CPU1(),
                    check_measure=true,
                    cache=true)

    range === nothing && error("You need to specify `range=...`.")
    model == nothing && error("You need to specify model=... .\n"*
                              "If `tuning=Explicit()`, any model in the "*
                              "range will do. ")

    if model isa Deterministic
        tuned_model = DeterministicTunedModel(model, tuning, resampling,
                                              measure, weights, operation,
                                              range, selection_heuristic,
                                              train_best, repeats, n,
                                              acceleration,
                                              acceleration_resampling,
                                              check_measure,
                                              cache)
    elseif model isa Probabilistic
        tuned_model = ProbabilisticTunedModel(model, tuning, resampling,
                                              measure, weights, operation,
                                              range, selection_heuristic,
                                              train_best, repeats, n,
                                              acceleration,
                                              acceleration_resampling,
                                              check_measure,
                                              cache)
    else
        error("Only `Deterministic` and `Probabilistic` "*
              "model types supported.")
    end

    message = clean!(tuned_model)
    isempty(message) || @info message

    return tuned_model

end

function MLJBase.clean!(tuned_model::EitherTunedModel)
    message = ""
    if tuned_model.measure === nothing
        tuned_model.measure = default_measure(tuned_model.model)
        if tuned_model.measure === nothing
            error("Unable to deduce a default measure for specified model. "*
                  "You must specify `measure=...`. ")
        else
            message *= "No measure specified. "*
            "Setting measure=$(tuned_model.measure). "
        end
    end

    if !supports_heuristic(tuned_model.tuning, tuned_model.selection_heuristic)
        message *= "`selection_heuristic=$(tuned_model.selection_heuristic)` "*
        "is not supported by $(tuned_model.tuning). Resetting to "*
        "`NaiveSelectionment()`."
        tuned_model.selection_heuristic = NaiveSelection()
    end

    if (tuned_model.acceleration isa CPUProcesses &&
        tuned_model.acceleration_resampling isa CPUProcesses)
        message *=
        "The combination acceleration=$(tuned_model.acceleration) and"*
        " acceleration_resampling=$(tuned_model.acceleration_resampling) is"*
        "  not generally optimal. You may want to consider setting"*
        " `acceleration = CPUProcesses()` and"*
        " `acceleration_resampling = CPUThreads()`."
    end

    if (tuned_model.acceleration isa CPUThreads &&
        tuned_model.acceleration_resampling isa CPUProcesses)
        message *=
        "The combination acceleration=$(tuned_model.acceleration) and"*
        " acceleration_resampling=$(tuned_model.acceleration_resampling) isn't"*
        " supported. \n Resetting to"*
        " `acceleration = CPUProcesses()` and"*
        " `acceleration_resampling = CPUThreads()`."

        tuned_model.acceleration = CPUProcesses()
        tuned_model.acceleration_resampling = CPUThreads()
    end

    tuned_model.acceleration =
        _process_accel_settings(tuned_model.acceleration)

    return message
end

## FIT AND UPDATE METHODS

# A *metamodel* is either a `Model` instance, `model`, or a tuple
# `(model, s)`, where `s` is extra data associated with `model` that
# the tuning strategy implementation wants available to the `result`
# method for recording in the history.

_first(m::MLJBase.Model) = m
_last(m::MLJBase.Model) = nothing
_first(m::Tuple{Model,Any}) = first(m)
_last(m::Tuple{Model,Any}) = last(m)

# returns a (model, result) pair for the history:
function event(metamodel,
               resampling_machine,
               verbosity,
               tuning,
               history,
               state)
    model = _first(metamodel)
    metadata = _last(metamodel)
    resampling_machine.model.model = model
    verb = (verbosity >= 2 ? verbosity - 3 : verbosity - 1)
    fit!(resampling_machine, verbosity=verb)
    E = evaluate(resampling_machine)
    entry0 = (model       = model,
              measure     = E.measure,
              measurement = E.measurement,
              per_fold    = E.per_fold,
              metadata    = metadata)
    entry = merge(entry0, extras(tuning, history, state, E))
    if verbosity > 2
        println("hyperparameters: $(params(model))")
    end

    if verbosity > 1
        println("measurement: $(E.measurement[1])")
    end
    return entry
end

function assemble_events(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         acceleration::CPU1)

     n_metamodels = length(metamodels)

     p = Progress(n_metamodels,
         dt = 0,
         desc = "Evaluating over $(n_metamodels) metamodels: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 25,
         color = :yellow)

    verbosity !=1 || update!(p,0)

    entries = map(metamodels) do m
        r = event(m, resampling_machine, verbosity, tuning, history, state)
        verbosity < 1 || begin
                  p.counter += 1
                  ProgressMeter.updateProgress!(p)
                end
        r
      end

    return entries
end

function assemble_events(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         acceleration::CPUProcesses)

    n_metamodels = length(metamodels)

    entries = @sync begin
        channel = RemoteChannel(()->Channel{Bool}(min(1000, n_metamodels)), 1)
        p = Progress(n_metamodels,
                     dt = 0,
                     desc = "Evaluating over $n_metamodels metamodels: ",
                     barglyphs = BarGlyphs("[=> ]"),
                     barlen = 25,
                     color = :yellow)

        # printing the progress bar
        verbosity < 1 || begin
            update!(p,0)
            @async while take!(channel)
                p.counter +=1
                ProgressMeter.updateProgress!(p)
            end
        end

        ret = @distributed vcat for m in metamodels
            r = event(m, resampling_machine, verbosity, tuning, history, state)
            verbosity < 1 || begin
                put!(channel, true)
            end
            r
        end
        verbosity < 1 || put!(channel, false)
        ret
    end

    return entries
end

@static if VERSION >= v"1.3.0-DEV.573"
# one machine for each thread; cycle through available threads:
function assemble_events(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         acceleration::CPUThreads)

    if Threads.nthreads() == 1
        return assemble_events(metamodels,
                         resampling_machine,
                         verbosity,
                         tuning,
                         history,
                         state,
                         CPU1())
   end

    n_metamodels = length(metamodels)
    ntasks = acceleration.settings
    partitions = chunks(1:n_metamodels, ntasks)
    #tasks = Vector{Task}(undef, length(partitions))
    entries = Vector(undef, length(partitions))
    p = Progress(n_metamodels,
         dt = 0,
         desc = "Evaluating over $(n_metamodels) metamodels: ",
         barglyphs = BarGlyphs("[=> ]"),
         barlen = 25,
         color = :yellow)
    ch = Channel{Bool}(min(1000, length(partitions)) )

    @sync begin
        # printing the progress bar
        verbosity < 1 || begin
            update!(p,0)
            @async while take!(ch)
                p.counter +=1
                ProgressMeter.updateProgress!(p)
            end
        end
        # One resampling_machine per task
         machs = [resampling_machine,
                 [machine(Resampler(
                     model= resampling_machine.model.model,
                     resampling    = resampling_machine.model.resampling,
                     measure       = resampling_machine.model.measure,
                     weights       = resampling_machine.model.weights,
                     operation     = resampling_machine.model.operation,
                     check_measure = resampling_machine.model.check_measure,
                     repeats       = resampling_machine.model.repeats,
                     acceleration  = resampling_machine.model.acceleration,
                     cache         = resampling_machine.model.cache),
                          resampling_machine.args...; cache=false) for
                  _ in 2:length(partitions)]...]

        @sync for (i, parts) in enumerate(partitions)
            Threads.@spawn begin
                entries[i] =  map(metamodels[parts]) do m
                    r = event(m, machs[i],
                              verbosity, tuning, history, state)
                    verbosity < 1 || put!(ch, true)
                    r
                end
            end
        end
        verbosity < 1 || put!(ch, false)
    end
    reduce(vcat, entries)
end

end # of if VERSION ...

# history is intialized to `nothing` because it's type is not known.
_vcat(history, Δhistory) = vcat(history, Δhistory)
_vcat(history::Nothing, Δhistory) = Δhistory
_length(history) = length(history)
_length(::Nothing) = 0

# builds on an existing `history` until the length is `n` or the model
# supply is exhausted (method shared by `fit` and `update`). Returns
# the bigger history:
function build(history,
               n,
               tuning,
               model,
               state,
               verbosity,
               acceleration,
               resampling_machine)
    j = _length(history)
    models_exhausted = false
    while j < n && !models_exhausted
        metamodels, state  = models(tuning,
                                    model,
                                    history,
                                    state,
                                    n - j,
                                    verbosity)
        Δj = _length(metamodels)
        Δj == 0 && (models_exhausted = true)
        shortfall = n - Δj
        if models_exhausted && shortfall > 0 && verbosity > -1
            @info "Only $j (of $n) models evaluated.\n"*
            "Model supply exhausted. "
        end
        Δj == 0 && break
        shortfall < 0 && (metamodels = metamodels[1:n - j])
        j += Δj

        Δhistory = assemble_events(metamodels,
                                   resampling_machine,
                                   verbosity,
                                   tuning,
                                   history,
                                   state,
                                   acceleration)
        history = _vcat(history, Δhistory)
    end
    return history, state
end

# given complete history, pick out best model, fit it on all data and
# generate report and cache (meta_state):
function finalize(tuned_model, history, state, verbosity, rm, data...)
    model = tuned_model.model
    tuning = tuned_model.tuning

    user_history = map(history) do entry
        delete(entry, :metadata)
    end

    entry =  best(tuned_model.selection_heuristic, history)
    best_model = entry.model
    best_history_entry = delete(entry, :metadata)
    fitresult = machine(best_model, data...)

    report0 = (best_model         = best_model,
               best_history_entry = best_history_entry,
               history            = user_history)

    if tuned_model.train_best
        fit!(fitresult, verbosity=verbosity - 1)
        report1 = merge(report0, (best_report=MLJBase.report(fitresult),))
    else
        report1 = merge(report0, (best_report=missing,))
    end

    report = merge(report1, tuning_report(tuning, history, state))
    meta_state = (history, deepcopy(tuned_model), state, rm)

    return fitresult, meta_state, report
end

function MLJBase.fit(tuned_model::EitherTunedModel{T,M},
                     verbosity::Integer, data...) where {T,M}
    tuning = tuned_model.tuning
    model = tuned_model.model
    range = tuned_model.range
    n = tuned_model.n === nothing ?
        default_n(tuning, range) : tuned_model.n

    verbosity < 1 || @info "Attempting to evaluate $n models."

    acceleration = tuned_model.acceleration

    state = setup(tuning, model, range, tuned_model.n, verbosity)

    # instantiate resampler (`model` to be replaced with mutated
    # clones during iteration below):
    resampler = Resampler(model=model,
                          resampling    = tuned_model.resampling,
                          measure       = tuned_model.measure,
                          weights       = tuned_model.weights,
                          operation     = tuned_model.operation,
                          check_measure = tuned_model.check_measure,
                          repeats       = tuned_model.repeats,
                          acceleration  = tuned_model.acceleration_resampling,
                          cache         = tuned_model.cache)
    resampling_machine = machine(resampler, data...; cache=false)
    history, state = build(nothing, n, tuning, model, state,
                           verbosity, acceleration, resampling_machine)

    rm = resampling_machine
    return finalize(tuned_model, history, state, verbosity, rm, data...)

end

function MLJBase.update(tuned_model::EitherTunedModel,
                        verbosity::Integer,
                        old_fitresult, old_meta_state, data...)

    history, old_tuned_model, state, resampling_machine = old_meta_state
    acceleration = tuned_model.acceleration

    tuning = tuned_model.tuning
    range = tuned_model.range
    model = tuned_model.model

    # exclamation points are for values actually used rather than
    # stored:
    n! = tuned_model.n === nothing ?
        default_n(tuning, range) : tuned_model.n

    old_n! = old_tuned_model.n === nothing ?
        default_n(tuning, range) : old_tuned_model.n

    if MLJBase.is_same_except(tuned_model, old_tuned_model, :n) &&
        n! >= old_n!

        verbosity < 1 || @info "Attempting to add $(n! - old_n!) models "*
        "to search, bringing total to $n!. "

        history, state = build(history, n!, tuning, model, state,
                               verbosity, acceleration, resampling_machine)

        rm = resampling_machine
        return finalize(tuned_model, history, state, verbosity, rm, data...)
    else
        return  fit(tuned_model, verbosity, data...)
    end
end

MLJBase.predict(tuned_model::EitherTunedModel, fitresult, Xnew) =
    predict(fitresult, Xnew)

function MLJBase.fitted_params(tuned_model::EitherTunedModel, fitresult)
    if tuned_model.train_best
        return (best_model=fitresult.model,
                best_fitted_params=fitted_params(fitresult))
    else
        return (best_model=fitresult.model,
                best_fitted_params=missing)
    end
end

## SUPPORT FOR MLJ ITERATION API

MLJBase.iteration_parameter(::Type{<:EitherTunedModel}) = :n
MLJBase.supports_training_losses(::Type{<:EitherTunedModel}) = true

function MLJBase.training_losses(tuned_model::EitherTunedModel, _report)
    _losses = MLJTuning.losses(tuned_model.selection_heuristic, _report.history)
    MLJTuning._length(_losses) == 0 && return nothing

    ret = similar(_losses)
    lowest = first(_losses)
    for i in eachindex(_losses)
        current = _losses[i]
        lowest = min(current, lowest)
        ret[i] = lowest
    end
    return ret
end

## METADATA

MLJBase.is_wrapper(::Type{<:EitherTunedModel}) = true
MLJBase.supports_weights(::Type{<:EitherTunedModel{<:Any,M}}) where M =
    MLJBase.supports_weights(M)
MLJBase.load_path(::Type{<:ProbabilisticTunedModel}) =
    "MLJTuning.ProbabilisticTunedModel"
MLJBase.load_path(::Type{<:DeterministicTunedModel}) =
    "MLJTuning.DeterministicTunedModel"
MLJBase.package_name(::Type{<:EitherTunedModel}) = "MLJTuning"
MLJBase.package_uuid(::Type{<:EitherTunedModel}) =
    "03970b2e-30c4-11ea-3135-d1576263f10f"
MLJBase.package_url(::Type{<:EitherTunedModel}) =
    "https://github.com/alan-turing-institute/MLJTuning.jl"
MLJBase.package_license(::Type{<:EitherTunedModel}) = "MIT"
MLJBase.is_pure_julia(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.is_pure_julia(M)
MLJBase.input_scitype(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.input_scitype(M)
MLJBase.target_scitype(::Type{<:EitherTunedModel{T,M}}) where {T,M} =
    MLJBase.target_scitype(M)


text:
"""
    $(SIGNATURES)

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.

DevNotes
- TODO, standardize around ::MIME="image/svg", see JuliaRobotics/DistributedFactorGraphs.jl#640
"""
function plotLocalProduct(fgl::AbstractDFG,
                          lbl::Symbol;
                          solveKey::Symbol=:default,
                          N::Int=100,
                          dims::Vector{Int}=Int[],
                          levels::Int=1,
                          show::Bool=false,
                          dirpath="/tmp/",
                          mimetype::AbstractString="svg",
                          sidelength=20cm,
                          title::String="Local product ($solveKey), ",
                          xmin=nothing, xmax=nothing, ymin=nothing, ymax=nothing  )
  #
  @warn "not showing partial constraints, but included in the product"
  arr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getBelief(getVariable(fgl, lbl), solveKey))
  push!(lbls, "curr")
  pl = nothing
  pp, parr, lb, sinfd = IIF.localProduct(fgl, lbl, N=N, solveKey=solveKey)

  partials = []

  # another sanity check
  xmin !== nothing && xmax !== nothing && xmin == xmax ? error("xmin must be less than xmax") : nothing
  ymin !== nothing && ymax !== nothing && ymin == ymax ? error("ymin must be less than ymax") : nothing

  # helper functions
  function plotDirectProducts()
      if pp != parr[1]
        push!(arr,pp)
        push!(lbls, "prod")
        for a in parr
          push!(arr, a)
        end
        @show lb, lbls
        lbls = union(lbls, string.(lb))
      end
      dims = length(dims) > 0 ? dims : collect(1:Ndim(pp))
      colors = getColorsByLength(length(arr))
      plotKDE(arr, dims=dims, levels=levels, c=colors, title=string(title,lbl), legend=string.(lbls)) #
  end
  function plotPartialProducts()
      # stack 1d plots to accomodate all the partials
      PL = []
      lbls = String["prod";"curr";string.(lb)]
      pdims = sort(collect(keys(partials)))
      for dimn in pdims
        vals = partials[dimn]
        proddim = marginal(pp, [dimn])
        colors = getColorsByLength(length(vals)+2)
        newbel_ = getBelief(getVariable(fgl, lbl),solveKey)
        newbel = newbel_ isa ManifoldKernelDensity ? newbel_.belief : newbel_
        pl = plotKDE( [proddim;newbel;vals], dims=[1;], levels=levels, c=colors, title=string("Local product, dim=$(dimn), ",lbl))
        push!(PL, pl)
      end
      Gadfly.vstack(PL...)
  end

  @show length(parr), length(partials)
  if length(parr) > 0 && length(partials) == 0
    pl = plotDirectProducts()
  elseif length(parr) == 0 && length(partials) > 0
    pl = plotPartialProducts()
  else
    return error("plotLocalProduct not built for lengths parr, partials = $(length(parr)), $(length(partials)) yet.")
  end

  # # set coordinates accordingly
  # co = Coord.Cartesian(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax)
  # pl.coord = co

  # now let's export:
  backend = getfield(Gadfly, Symbol(uppercase(mimetype)))
  Gadfly.draw(backend(dirpath*"test_$(lbl).$(mimetype)",sidelength,0.75*sidelength), pl)
  driver = mimetype in ["pdf"] ? "evince" : "eog"
  show ? (@async run(`$driver $(dirpath)test_$(lbl).$(mimetype)`)) : nothing

  return pl
end

"""
    $(SIGNATURES)

Plot---for the cylindrical manifold only---the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference.  The linear and circular belief products are returned as a tuple.
"""
function plotLocalProductCylinder(fgl::G,
                                  lbl::Symbol;
                                  N::Int=100,
                                  levels::Int=1,
                                  show=false,
                                  dirpath="/tmp/",
                                  mimetype::AbstractString="svg",
                                  sidelength=30cm,
                                  scale::Float64=0.2 ) where G <: AbstractDFG
  #
  @warn "not showing partial constraints, but included in the product"
  arr = Array{BallTreeDensity,1}()
  carr = Array{BallTreeDensity,1}()
  lbls = String[]
  push!(arr, getBelief(getVariable(fgl, lbl)))
  push!(lbls, "curr")
  pl = nothing
  plc = nothing
  pp, parr, partials, lb = IncrementalInference.localProduct(fgl, lbl, N=N)
  if length(parr) > 0 && length(partials) == 0
    if pp != parr[1]
      push!(arr, marginal(pp,[1]))
      push!(lbls, "prod")
      push!(carr, marginal(pp, [2]))
      for a in parr
        push!(arr, marginal(a,[1]))
        push!(carr, marginal(a,[2]))
      end
      lbls = union(lbls, lb)
    end
    colors = getColorsByLength(length(arr))
    pll = plotKDE(arr, dims=[1;], levels=levels, c=colors, legend=lbls, title=string("Local product, ",lbl))
    plc = AMP.plotKDECircular(carr, c=colors, scale=scale)
    # (pll, plc)
  elseif length(parr) == 0 && length(partials) > 0
    # stack 1d plots to accomodate all the partials
    PL = []
    PLC = []
    lbls = ["prod";"curr";lb]
    pdims = sort(collect(keys(partials)))
    for dimn in pdims
      vals = partials[dimn]
      proddim = marginal(pp, [dimn])
      colors = getColorsByLength(length(vals)+2)
      if dimn == 1
        pll = plotKDE([proddim;getBelief(getVariable(fgl, lbl));vals], dims=[1;], levels=levels, c=colors, title=string("Local product, dim=$(dimn), ",lbl))
        push!(PL, pl)
      else
        plc = AMP.plotKDECircular([proddim; marginal(getBelief(getVariable(fgl, lbl)),[2]);vals], c=colors, scale=scale)
        push!(PLC, plc)
      end
    end
    pl = Gadfly.vstack(PL..., PLC...)
  else
    return error("plotLocalProduct not built for lengths parr, partials = $(length(parr)), $(length(partials)) yet.")
  end

  # now let's export:
  # backend = getfield(Gadfly, Symbol(uppercase(mimetype)))
  # Gadfly.draw(backend(dirpath*"test_$(lbl).$(mimetype)",sidelength,sidelength), pl)
  # driver = mimetype in ["pdf"] ? "evince" : "eog"
  # show ? (@async run(`$driver $(dirpath)test_$(lbl).$(mimetype)`)) : nothing

  return pll, plc
end

"""
    $SIGNATURES

Plot the proposal belief from neighboring factors to `lbl` in the factor graph (ignoring Bayes tree representation),
and show with new product approximation for reference. String version is obsolete and will be deprecated.
"""
plotLocalProduct(fgl::AbstractDFG, lbl::AbstractString; N::Int=100, dims::Vector{Int}=Int[] ) = plotLocalProduct(fgl, Symbol(lbl), N=N, dims=dims)

"""
    $SIGNATURES

Development function to plot the same variable from both factor graphs together.
"""
function plotPairVariables(dfg1::G1,
                          dfg2::G2,
                          sym::Symbol;
                          dims=nothing,
                          levels::Int=3,
                          title::String="") where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  X1 = getBelief(dfg1, sym)
  X2 = getBelief(dfg2, sym)

  plotKDE([X1;X2], c=["black";"red"], legend=["dfg1";"dfg2"], dims=dims, levels=levels, title=title*" $sym")
end

function plotPairPose2(dfg1::G1,
                          dfg2::G2,
                          sym::Symbol;
                          dims=nothing,
                          levels::Int=3,
                          title::String="") where {G1 <: AbstractDFG, G2 <: AbstractDFG}
  #
  X1 = getBelief(dfg1, sym)
  X2 = getBelief(dfg2, sym)

  plotPose(Pose2(), [X1;X2], title*" $sym", c=["black";"red"], levels=levels)
  # , legend=["dfg1";"dfg2"]
end

"""
    $SIGNATURES

Plot new proposal (convolution) for factor x of variable y given all other -- i.e. y|X

Notes
- plot new result on `towards` as first color.
- Plot other variables in factor on 2nd to m colors.
"""
function plotVariableGivenFactor(dfg::G,
                                fct::Symbol,
                                towards::Symbol;
                                levels::Int=2,
                                dims=nothing  ) where G <: AbstractDFG
  #
  pts = approxConv(dfg,fct,towards)
  mani = getManifold(dfg, towards)
  res = manikde!(mani, pts)

  lie = ls(dfg, fct)
  setdiff!(lie, [towards])

  otr = map(x->getBelief(dfg,x),lie)
  lbls = string.([towards;lie])

  pl = plotKDE([res;otr],dims=dims,levels=levels,legend=lbls)

  return pl
end

#


text:
using TreeShakeTest

TreeShakeTest.foo()


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ d7e3ef54-4307-11eb-221c-e354f095486e
md" ## _Equilibrium Extension: Acids and Bases_
"

# ╔═╡ 45a52da0-4308-11eb-1765-9fc9cc901eb8
md" ## Acid and Base Theories
Originally, the arrhenius definition was first used to understand acid and bases. It said that acids donate H+ ions and bases donate OH- ions. However, people figured out that this did not apply for all bases, like NH3+. So, the definition was revised to with the Bronsted-Lowry definition which says:

Acids donate protons (h+ ions)

Bases accept protons (h+ ions)

**Conjugate acid-base pair**: this a pair that differ with the prescence of an h+ ion. Some examples of these pairs are: $\ce{NH3 and NH4+, H2PO4- and H3PO4}$. The conjugate pair is the resulting species after the h+ is either donated or accepted. The formation of different species dictates whether the solution is acidic or basic. If H3O+ is formed, then the solution is acidic. If OH- is formed, then the solution is neutral. If equal amounts of H+ and OH- are formed, then the solution is neutral.

*Acids only donate one proton at a time*
"

# ╔═╡ 4cea74b2-4308-11eb-2c18-a9ffb737a3fd
md" ## Strength of Acids and Bases
"

# ╔═╡ 5eac15ac-4308-11eb-0970-0db13817d006
md" ## Auto-Ionization
"

# ╔═╡ 76169ee2-4308-11eb-06a3-8d1cf54f2856
md" ## pH and pOH
"

# ╔═╡ 7dce9ee6-4308-11eb-0f9b-69e65e4c2c1b
md"## Acid-Base Properties of Salt
"

# ╔═╡ 952f1e30-4308-11eb-0a9d-4911391ca300

# ╔═╡ Cell order:
# ╠═d7e3ef54-4307-11eb-221c-e354f095486e
# ╠═45a52da0-4308-11eb-1765-9fc9cc901eb8
# ╠═4cea74b2-4308-11eb-2c18-a9ffb737a3fd
# ╠═5eac15ac-4308-11eb-0970-0db13817d006
# ╠═76169ee2-4308-11eb-06a3-8d1cf54f2856
# ╠═7dce9ee6-4308-11eb-0f9b-69e65e4c2c1b
# ╠═952f1e30-4308-11eb-0a9d-4911391ca300


text:
@testset "Quoted Symbol Syntax" begin
    @test (Node(:a) → Node(:b)) == (StenoGraphs.@variable_as_node a → b)
    b = Node(:c)
    @test (Node(:a) → Node(:c)) == (StenoGraphs.@variable_as_node a → _(b))
    let err = nothing
        try
            eval(:(StenoGraphs.@variable_as_node a → _(b, c)))
        catch err
        end

        @test err isa Exception
        @test  occursin("Unqote only a single argument. Right: `_(x)`, Wrong: `_(x, y)`.", sprint(showerror, err))
    end
end

@testset "Addition as hcat" begin
    @test StenoGraphs.addition_to_vector!(:(a * b + c)) == :([a * b c])
    @test StenoGraphs.addition_to_vector!(:(a * (b + c))) == :(a * [b c])
end

@testset "Broadcasting Edges" begin
    @test Edge(Node(:a), [Node(:b) Node(:c)]) == vec([Edge(Node(:a), Node(:b)) Edge(Node(:a), Node(:c))])
    @test Edge([Node(:a) Node(:b)], Node(:c)) == vec([Edge(Node(:a), Node(:c)) Edge(Node(:b), Node(:c))])
    @test Edge([Node(:a) Node(:b)], [Node(:c) Node(:d)]) == vec([Edge(Node(:a), Node(:c)) Edge(Node(:a), Node(:d)) Edge(Node(:b), Node(:c)) Edge(Node(:b), Node(:d))])
end

@testset "Singleline Arrows" begin
    @test [Edge(Node(:a), Node(:c)), Edge(Node(:b), Node(:c))] == @StenoGraph [a b] → c
end

@testset "Multiline Arrows" begin
    @test [Edge(Node(:a), Node(:c)), Edge(Node(:b), Node(:c)), Edge(Node.(:f), Node.(:e))] == @StenoGraph begin
        [a b] → c
        e ← f
    end
end

@testset "Multiline Mix Arrow Edge" begin
    @test [Edge(Node(:a), Node(:c)), Edge(Node(:b), Node(:c)), Edge(Node(:f), Node(:e))] ==
    @StenoGraph begin
        [a b] → c
        Edge(f, e)
    end
end


text:
# # Feature Extraction

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/feature_extraction.ipynb)

# In this example we demonstrate how to extract feature from trajectories using AutomotiveSimulator

# **Load a Dataset**
# First let's create a synthetic dataset.

using AutomotiveSimulator
using AutomotiveVisualization
AutomotiveVisualization.colortheme["background"] = colorant"white"; # hide
using Random

roadway = roadway = gen_straight_roadway(3, 1000.0)
veh_state = VehicleState(Frenet(roadway[LaneTag(1,2)], 0.0), roadway, 10.)
veh1 = Entity(veh_state, VehicleDef(), "bob")
veh_state = VehicleState(Frenet(roadway[LaneTag(1,2)], 20.0), roadway, 2.)
veh2 = Entity(veh_state, VehicleDef(), "alice")

dt = 0.5
n_steps = 10
models = Dict{String, DriverModel}()
models["bob"] = Tim2DDriver(mlane=MOBIL())
set_desired_speed!(models["bob"], 10.0)
models["alice"] = Tim2DDriver(mlane=MOBIL())
set_desired_speed!(models["alice"], 2.0)

scene = Scene([veh1, veh2])
scenes = simulate(scene, roadway, models, n_steps, dt)
camera = SceneFollowCamera()
update_camera!(camera, scene)
snapshot = render([roadway, scene], camera=camera)
#md write("feature_initial.svg", snapshot) # hide
#md # ![initial state of feature extraction scenario](feature_initial.svg)

#md nothing # hide

# One can also load the data from the Stadium tutorial

#md load data from stadium tutorial
#md scenes = open("2Dstadium_listrec.txt", "r") do io
#md     read(io, Vector{EntityScene{VehicleState, VehicleDef, String}})
#md end

# **Extract features from a recorded trajectory**
# Recorded trajectories are expected to be vectors of `Scene`s where each element correspond to one time step.
# To extract features, one can use the `extract_features` function which takes as input
# a list of feature we want to extract and the list of vehicle ids for which we want those features.
# For this example, let's first query two features, the longitudinal and lateral position of Bob, and whether or not Bob is colliding:

dfs = extract_features((posfs, posft, iscolliding), roadway, scenes, ["bob"])
dfs["bob"]

# To query features for all traffic participants we can just add their ID to the list:

dfs = extract_features((posfs, posft, iscolliding), roadway, scenes, ["bob", "alice"])
dfs["alice"]

# The output is a dictionary mapping ID to dataframes. To learn more about DataFrames visit [DataFrames.jl](https://github.com/JuliaData/DataFrames.jl).

# For the list of all possible features available see the documentation.
# Features are generally just function. AutomotiveSimulator provides some convenience
# to automatically generate feature function like `distance_to_$x`
# The `distance_to` function takes as input a vehicle ID and returns a function to extract
# the distance between the queried vehicle and the vehicle given to `distance_to`

distance_to("alice")

# we can use this newly generated funciton in the feature extraction pipeline

dfs = extract_features((distance_to_alice, posfs), roadway, scenes, ["bob"])
dfs["bob"].distance_to_alice[1] # distance between Bob and Alice in the first scene.


text:
using Documenter

using Pkg
Pkg.activate("../")
using JMcDM

push!(LOAD_PATH,"../src/")
makedocs(
	sitename="JMcDM",
	pages=[
		"mcdms.md", "game.md", "dataenvelop.md", "scdm.md", "utility.md"
	]
)

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
@testset "Infeasible Problem" begin
prob, opts = Problems.DubinsCar(:escape)
model_inf = Altro.InfeasibleModel.(prob.model)
@test RD.dims(model_inf) == (fill(3,prob.N), fill(5,prob.N))
Z0 = prob.Z
Z = Altro.infeasible_trajectory(model_inf, Z0)
@test all(z->length(z) == 8, Z)
conSet = TO.change_dimension(get_constraints(prob), 3, 5, 1:3, 1:2)
inf = Altro.InfeasibleConstraint(model_inf[1])
TO.add_constraint!(conSet, inf, 1:prob.N-1)
obj = Altro.infeasible_objective(prob.obj, 1.0)
@test RD.dims(obj) == (fill(3, prob.N), fill(5, prob.N))
prob_inf = Problem(model_inf, obj, conSet, prob.x0, prob.xf, Z, prob.N,
                   TO.get_initial_time(prob), TO.get_final_time(prob))
ilqr = Altro.iLQRSolver(prob_inf, opts, use_static=Val(true))
@test Altro.usestatic(ilqr)
ilqr = Altro.iLQRSolver(prob_inf, opts, use_static=Val(false))
@test Altro.usestatic(ilqr) == false
end


text:
using .CUDA, Adapt

struct CuSnpArray{T} <: AbstractMatrix{UInt8}
    data::CuMatrix{UInt8}
    m::Int
    model::Union{Val{1}, Val{2}, Val{3}}
    center::Bool
    scale::Bool
    impute::Bool
    μ::CuVector{T}
    σinv::CuVector{T}
    storagev1::CuVector{T}
    storagev2::CuVector{T}
end

"""
    CuSnpArray{T}(s; model=ADDITIVE_MODEL, center=false, scale=false, impute=true)

Copy a `SnpArray` to a CUDA GPU to perform linear algebera operations.

# Arguments
- s: a `SnpArray`.
- model: one of `ADDITIVE_MODEL`(default), `DOMINANT_MODEL`, `RECESSIVE_MODEL`.
- center: whether to center (default: false).
- scale: whether to scale to standard deviation 1 (default: false).
- impute: whether to impute missing value with column mean (default: true).
"""
function CuSnpArray{T}(s::SnpArray;
    model = ADDITIVE_MODEL,
    center::Bool = false,
    scale::Bool = false, impute::Bool = false) where T <: AbstractFloat

    data = adapt(CuMatrix{UInt8}, s.data)
    if model == ADDITIVE_MODEL
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=ADDITIVE_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j] / 2))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    elseif model == DOMINANT_MODEL
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=DOMINANT_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j]))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    elseif model == RECESSIVE_MODEL
        if center || scale
            μ = Vector{T}(undef, size(s, 2))
            μ[:] = mean(s; dims=1, model=RECESSIVE_MODEL)
        else
            μ = T[]
        end
        if scale
            σinv = Vector{T}(undef, size(s, 2))
            @inbounds @simd for j in 1:size(s, 2)
                σinv[j] = sqrt(μ[j] * (1 - μ[j]))
                σinv[j] = σinv[j] > 0 ? inv(σinv[j]) : one(T)
            end
        else
            σinv = T[]
        end
    else
        throw(ArgumentError("unrecognized model $model"))
    end
    μ = adapt(CuVector{T}, μ)
    σinv = adapt(CuVector{T}, σinv)
    storagev1 = CuVector{T}(undef, size(s, 1))
    storagev2 = CuVector{T}(undef, size(s, 2))
    CuSnpArray{T}(data, s.m, model, center, scale, impute, μ, σinv, storagev1, storagev2)
end

Base.size(s::CuSnpArray) = s.m, size(s.data, 2)
eltype(s::CuSnpArray) = eltype(s.μ)

for (_ftn!, expr) in [
        (:_snparray_cuda_ax_additive!, :(((Aij >= 2) + (Aij >= 3)) * v[j])),
        (:_snparray_cuda_ax_dominant!, :((Aij >= 2) * v[j])),
        (:_snparray_cuda_ax_recessive!, :((Aij >= 3) * v[j])),
        (:_snparray_cuda_ax_additive_meanimpute!, :(((Aij >= 2) + (Aij >= 3) + (Aij == 1) * μ[j]) * v[j])),
        (:_snparray_cuda_ax_dominant_meanimpute!, :(((Aij >= 2) + (Aij == 1) * μ[j]) * v[j])),
        (:_snparray_cuda_ax_recessive_meanimpute!, :(((Aij >= 3) + (Aij == 1) * μ[j]) * v[j]))
    ]
    @eval begin
        function $_ftn!(out, s, v)
            packedstride = size(s, 1)
            index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
            stride_x = blockDim().x * gridDim().x
            for i = index_x:stride_x:length(out)
                outi = zero(eltype(out))
                for j = 1:length(v)
                    k = 2 * ((i-1) & 3)
                    block = s[(j-1) * packedstride + ((i-1) >> 2) + 1]
                    Aij = (block >> k) & 3
                    @inbounds outi += $expr
                end
                out[i] = outi
            end
        end
    end
end

for (_ftn!, expr) in [
        (:_snparray_cuda_atx_additive!, :((((Aij >= 2) + (Aij >= 3))) * v[j])),
        (:_snparray_cuda_atx_dominant!, :((Aij >= 2) * v[j])),
        (:_snparray_cuda_atx_recessive!, :((Aij >= 3) * v[j])),
        (:_snparray_cuda_atx_additive_meanimpute!, :(((Aij >= 2) + (Aij >= 3) + (Aij == 1) * μ[i]) * v[j])),
        (:_snparray_cuda_atx_dominant_meanimpute!, :(((Aij >= 2) + (Aij == 1) * μ[i]) * v[j])),
        (:_snparray_cuda_atx_recessive_meanimpute!, :(((Aij >= 3) + (Aij == 1) * μ[i]) * v[j]))
    ]
    @eval begin
        function $_ftn!(out, s, v)
            packedstride = size(s, 1)
            index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
            stride_x = blockDim().x * gridDim().x
            for i = index_x:stride_x:length(out)
                outi = zero(eltype(out))
                for j = 1:length(v)
                    k = 2 * ((j-1) & 3)
                    block = s[(i-1) * packedstride + ((j-1) >> 2) + 1]
                    Aij = (block >> k) & 3
                    @inbounds outi += $expr
                end
                out[i] = outi
            end
        end
    end
end

"""
    LinearAlgebra.mul!(out::CuVector{T}, s::CuSnpArray{T}, v::CuVector{T})

In-place matrix-vector multiplication on a GPU.
"""
function mul!(
    out::CuVector{T},
    s::CuSnpArray{T},
    v::CuVector{T}) where T <: AbstractFloat
    @assert length(out) == size(s, 1) && length(v) == size(s, 2)
    fill!(out, zero(T))
    if s.scale
        s.storagev2 .= s.σinv .* v
        w = s.storagev2
    else
        w = v
    end
    numblocks = ceil(Int, length(out)/256)
    CUDA.@sync begin
        if s.model == ADDITIVE_MODEL
            @cuda threads=256 blocks=numblocks _snparray_cuda_ax_additive!(out, s.data, w)
        elseif s.model == DOMINANT_MODEL
            @cuda threads=256 blocks=numblocks _snparray_cuda_ax_dominant!(out, s.data, w)
        else
            @cuda threads=256 blocks=numblocks _snparray_cuda_ax_recessive!(out, s.data, w)
        end
    end

    if s.center
        return out .-= dot(s.μ, w)
    else
        return out
    end
end

"""
    LinearAlgebra.mul!(out::CuVector{T}, s::Union{Transpose{T, CuSnpArray{T}}, Adjoint{T, CuSnpArray{T}}}, v::CuVector{T})

In-place matrix-vector multiplication on a GPU, with transposed CuSnpArray.
"""
function mul!(
    out::CuVector{T},
    st::Union{Transpose{T, CuSnpArray{T}}, Adjoint{T, CuSnpArray{T}}},
    v::CuVector{T}) where T <: AbstractFloat
    @assert length(out) == size(st, 1) && length(v) == size(st, 2)
    fill!(out, zero(T))
    s = st.parent
    numblocks = ceil(Int, length(out)/256)
    CUDA.@sync begin
        if !s.impute
            if s.model == ADDITIVE_MODEL
                @cuda threads=256 blocks=numblocks _snparray_cuda_atx_additive!(out, s.data, v)
            elseif s.model == DOMINANT_MODEL
                @cuda threads=256 blocks=numblocks _snparray_cuda_atx_dominant!(out, s.data, v)
            else
                @cuda threads=256 blocks=numblocks _snparray_cuda_atx_recessive!(out, s.data, v)
            end
        else
            if s.model == ADDITIVE_MODEL
                @cuda threads=256 blocks=numblocks _snparray_cuda_atx_additive_meanimpute!(out, s.data, v, s.μ)
            elseif s.model == DOMINANT_MODEL
                @cuda threads=256 blocks=numblocks _snparray_cuda_atx_dominant_meanimpute!(out, s.data, v, s.μ)
            else
                @cuda threads=256 blocks=numblocks _snparray_cuda_atx_recessive_meanimpute!(out, s.data, v, s.μ)
            end
        end
    end

    if s.center
        out .-= sum(v) .* s.μ
    end
    if s.scale
        return out .*= s.σinv
    else
        return out
    end
end


text:
# Read span and put spikes into poi_samples collection.

function read_spike_samples(samples::SampleData{T}, model::ModelData, file::String, span::Int64) where {T <: Number}
    synapses = Model.synapses(model)
    span_time = Model.span_time(model)

    # println("Loading new samples: ", file)

    # Load samples
    idx = 1
    open(file, "r") do f
        syn_samples = readlines(f)
        for syn_sample in syn_samples
            # Format:
            # id 1010001011010...::

            # Parse out "id" field
            idx_range = findfirst(" ", syn_sample)
            id = parse(Int64, SubString(syn_sample, 1, idx_range[1]))
            # print(id, " ")

            # EOL marker
            bits_end = findfirst("::", syn_sample)

            # Extract just spike data
            bits = SubString(syn_sample, idx_range[1] + 1, bits_end[1] - 1)

            # load row of spikes into samples
            t = samples.t
            for bit in bits
                samples.samples[t] = parse(T, bit)
                t += 1
            end
            idx += 1
        end
    end

    # Prepare for next span my moving "t" to the start of the next span
    # position within the full duration of samples.
    samples.t += span_time
end

function read_poi_samples(samples::Samples, model::ModelData, span::Int64)
    # Where to access fresh samples
    path = Model.data_output_path(model)

    # source file.
    file = path * Model.poisson_files(model) * string(span) * ".data"

    synapses = Model.synapses(model)
    span_time = Model.span_time(model)

    # println("Loading new samples: ", file)

    # Load samples
    idx = 1
    open(file, "r") do f
        syn_samples = readlines(f)
        for syn_sample in syn_samples
            # Format:
            # id 1010001011010...::

            # Parse out "id" field
            idx_range = findfirst(" ", syn_sample)
            id = parse(Int64, SubString(syn_sample, 1, idx_range[1]))
            # print(id, " ")

            # EOL marker
            bits_end = findfirst("::", syn_sample)

            # Extract just spike data
            bits = SubString(syn_sample, idx_range[1] + 1, bits_end[1] - 1)

            # load row of spikes into samples
            t = samples.poi_t
            for bit in bits
                samples.poi_samples[idx, t] = parse(UInt8, bit)
                t += 1
            end
            idx += 1
        end
    end

    # Prepare for next span my moving "t" to the start of the next span
    # position within the full duration of samples.
    samples.poi_t += span_time
end

function read_stimulus_samples(samples::Samples, model::ModelData, span::Int64)
    # Where to access fresh samples
    path = Model.data_output_path(model)

    # source file.
    file = path * Model.output_stimulus_files(model) * string(span) * ".data"

    span_time = Model.span_time(model)

    # println("Loading new samples: ", file)

    # Load samples
    idx = 1
    open(file, "r") do f
        syn_samples = readlines(f)
        for syn_sample in syn_samples
            # Format:
            # id 1010001011010...::

            # Parse out "id" field
            idx_range = findfirst(" ", syn_sample)
            id = parse(Int64, SubString(syn_sample, 1, idx_range[1]))
            # print(id, " ")

            # EOL marker
            bits_end = findfirst("::", syn_sample)

            # Extract just spike data
            bits = SubString(syn_sample, idx_range[1] + 1, bits_end[1] - 1)
            # println(bits)
            # load row of spikes into samples
            t = samples.stim_t
            for bit in bits
                samples.stimulus_samples[idx, t] = parse(UInt8, bit)
                t += 1
            end
            idx += 1
        end
    end

    # Prepare for next span my moving "t" to the start of the next span
    # position within the full duration of samples.
    samples.stim_t += span_time
end

function read_float_samples(data_sams::SampleData, model::ModelData, file::String, span::Int64)
    span_time = Model.span_time(model)

    # println("Loading apFast samples: ", file)

    # Load samples
    open(file, "r") do f
        data = readlines(f)
        t = data_sams.t
        for sample in data
            # Format:
            # float
            # float
            # ...
            v = parse(Float64, sample)
            data_sams.min = min(data_sams.min, v)
            data_sams.max = max(data_sams.max, v)

            data_sams.samples[t] = v
            t += 1
        end
    end

    # println("min: ", data_sams.min, ", max: ", data_sams.max)
    # Prepare for next span my moving "t" to the start of the next span
    # position within the full duration of samples.
    data_sams.t += span_time
end

function read_synaptic_float_samples(syn_samples::SynapticSamples, model::ModelData, file::String, span::Int64)
    span_time = Model.span_time(model)

    # Format:
    # id float,float,float...
    # id float,float,float...
    # Nth float,float,float...
    open(file, "r") do f
        synapses = Model.synapses(model)

        for syn_id in 1:synapses
            t = syn_samples.t # reset time index for next synapse

            line_data = readline(f)

            # Fetch samples bucket
            syn_data = syn_samples.data[syn_id]

            # Parse out "id" field
            idx_range = findfirst(" ", line_data)
            id = parse(Int64, SubString(line_data, 1, idx_range[1]))

            # Extract just weight data
            sub_line = SubString(line_data, idx_range[1] + 1, length(line_data))
            weights = split(sub_line, ' ')

            for weight in weights
                w = parse(Float64, weight)
                syn_data.min = min(syn_data.min, w)
                syn_data.max = max(syn_data.max, w)

                syn_data.samples[t] = w
                t += 1
            end
        end
    end

    syn_samples.t += span_time
end


text:
module InformationDistances

using Distances: PreMetric
import CodecXz
import TranscodingStreams

using LibDeflate

export NormalizedCompressionDistance, LibDeflateCompressor, CodecCompressor,
        compressed_length, compressed_lengths

include("compressors.jl")

"""
    NormalizedCompressionDistance{<: AbstractCompressor} <: Distances.PreMetric

A `normalized compression distance` metric between two strings.

The metric is defined by ``d(x, y) := \\frac{Z(xy) - \\min(Z(x), Z(y))} {\\max(Z(x), Z(y))}``

where `Z(x)` is the length when compressing the string `x` with a certain compression codec.

----

    NormalizedCompressionDistance(, [compressor::AbstractCompressor])

Create a `NormalizedCompressionDistance`.

# Arguments
- `compressor` The compressor to use. If not specified,
    `CodecCompressor{CodecXz.XzCompressor}(;level=9; check=CodecXz.LZMA_CHECK_NONE)` is used.

# Examples
```jldoctest
julia> d1 = NormalizedCompressionDistance()
NormalizedCompressionDistance{CodecCompressor{CodecXz.XzCompressor}}(CodecCompressor{CodecXz.XzCompressor}(Base.Iterators.Pairs{Symbol,Signed,Tuple{Symbol,Symbol},NamedTuple{(:level, :check),Tuple{Int64,Int32}}}(:level => 9,:check => 0)))

julia> d1("hello", "world")
0.07142857142857142

julia> d2 = NormalizedCompressionDistance(LibDeflateCompressor())
NormalizedCompressionDistance{LibDeflateCompressor}(LibDeflateCompressor(12))

julia> d2("hello", "world")
0.5
```
"""
struct NormalizedCompressionDistance{C <: AbstractCompressor} <: PreMetric

    compressor::C

    function NormalizedCompressionDistance(compressor=
            CodecCompressor{CodecXz.XzCompressor}(;level=9, check=CodecXz.LZMA_CHECK_NONE))

        return new{typeof(compressor)}(compressor)
    end
end

function (dist::NormalizedCompressionDistance)(x::AbstractString, y::AbstractString)

    x == y && return 0.0

    compressor = dist.compressor

    len_xy, len_x, len_y = compressed_lengths(compressor, (x * y, x, y))

    return clamp(Float64((len_xy - min(len_x, len_y)) / max(len_x, len_y)), 0.0, 1.0)
end

end # module


text:
struct PCloudError{T} <: Exception
    msg::T
end

function PCloudClient(; auth_token = "", user = "", password = "", set_endpoints = true)
    client = PCloudClient("https://api.pcloud.com/",
                          "https://api.pcloud.com/",
                          "https://api.pcloud.com/",
                          Dict{String, String}())
    if !isempty(auth_token)
        authorize!(client, auth_token)
    elseif !isempty(user) & !isempty(password)
        authorize!(client, user, password)
    end

    if set_endpoints
        update_endpoints!(client)
    end

    return client
end

function authorize!(client::PCloudClient, auth_token)
    client.cookies["auth"] = auth_token
end

function authorize!(client::PCloudClient, user, password)
    digest_json = getdigest(client)
    digest = digest_json.digest
    pwddigest = bytes2hex(sha1(password * bytes2hex(sha1(user)) * digest))
    res = userinfo(client; getauth = 1, logout = 1, username = user,
                   digest = digest, passworddigest = pwddigest)
    client.cookies["auth"] = res.auth
end

function authtoken(client::PCloudClient)
    return client.cookies["auth"]
end

function update_endpoints!(client::PCloudClient)
    res = API.getapiserver(client)
    client.apiep = "https://" * first(res.api) * "/"
    client.binep = "https://" * first(res.binapi) * "/"

    return client
end

function reset_endpoints!(client::PCloudClient)
    client.apiep = client.rootep
    client.binep = client.rootep

    return client
end

function useglobally!(client::PCloudClient)
    OPTS.client = client
end


text:
using Mass2
using ReferenceMicrocalFiles
using ZMQ
using Base.Test

function setup_channel(ljh_filename, noise_filename)
	ljh = LJHGroup(ljh_filename)

	mc=MassChannel()

	mc[:pretrig_mean] = RunningVector(Float32)
	mc[:pretrig_rms] = RunningVector(Float32)
	mc[:pulse_average] = RunningVector(Float32)
	mc[:pulse_rms] = RunningVector(Float32)
	mc[:rise_time] = RunningVector(Float32)
	mc[:postpeak_deriv] = RunningVector(Float32)
	mc[:peak_index] = RunningVector(UInt16)
	mc[:peak_value] = RunningVector(UInt16)
	mc[:min_value] = RunningVector(UInt16)
	mc[:filt_value] = RunningVector(Float32)
	mc[:filt_value_dc] = RunningVector(Float32)
	mc[:filt_phase] = RunningVector(Float32)
	mc[:energy] = RunningVector(Float32)
	mc[:pulse] = RunningVector(Vector{UInt16})
	mc[:rowcount] = RunningVector(Int)
	mc[:timestamp_posix_usec] = RunningVector(Int)
	mc[:selection_good] = RunningSumBitVector()
	mc[:filt_value_hist] = Histogram(0:1:20000)
	mc[:filt_value_dc_hist] = Histogram(0:1:20000)
	mc[:energy_hist] = Histogram(0:1:20000)

	mc[:ljh] = LJHGroup(ljh_filename)
	mc[:pretrig_nsamples] = LJH.pretrig_nsamples(mc[:ljh])
	mc[:samples_per_record] = LJH.record_nsamples(mc[:ljh])
	mc[:frametime] = LJH.frametime(mc[:ljh])
	mc[:channel_number] = LJH.channel(mc[:ljh])

	mc[:f_3db] = 10000
	mc[:known_energies] = [5898.801,6490.59] # FeKAlpha and FeKBeta

	#metadata
	mc[:calibration_nextra] = 1 # when finding peaks, how many peaks other than the largest n to include when assigning peaks to energies

	mc[:noise_filename]=ascii(noise_filename)
	mc[:ljh_filename]=ascii(ljh_filename) # h5py can't read UFT8String written by julia
	mc[:name] = "summarize and filter test"
	mc[:hdf5_filename] = "$(splitext(ljh_filename)[1])_jl.hdf5"
	mc[:oncleanfinish] = markhdf5oncleanfinish
	isfile(mc[:hdf5_filename]) && rm(mc[:hdf5_filename])

	steps = AbstractStep[
	GetPulsesStep(mc[:ljh], [:pulse, :rowcount, :timestamp_posix_usec], 0,100,0,false) # the false is for analyzing files on disk that aren't being written
	@perpulse pretrig_mean, pretrig_rms, pulse_average, pulse_rms, rise_time, postpeak_deriv, peak_index, peak_value, min_value = compute_summary(pulse, pretrig_nsamples, frametime)
	@threshold pretrig_rms_criteria, postpeak_deriv_criteria = estimate_pretrig_rms_and_postpeak_deriv_criteria(noise_filename, pretrig_nsamples) when length(pulse) > 100
	@threshold peak_index_criteria = estimate_peak_index_criteria(peak_index) when length(peak_index)>100
	@perpulse selection_good = selectfromcriteria(pretrig_rms, pretrig_rms_criteria, peak_index, peak_index_criteria, postpeak_deriv, postpeak_deriv_criteria)
	@threshold average_pulse = compute_average_pulse(pulse, selection_good) when sum(selection_good) > 100
	@threshold rise_tau_s, fall_tau_s = fit_pulse_two_exponential(average_pulse, pretrig_nsamples, frametime) when length(filt_value)>10
	@threshold noise_autocorr = compute_noise_autocorr(noise_filename, samples_per_record) when sum(selection_good) > 100
	@threshold filter, vdv = compute_filter(average_pulse, noise_autocorr, f_3db, frametime) when sum(selection_good) > 100
	@perpulse filt_value, filt_phase = filter5lag(filter, pulse)
	@threshold pretrigger_mean_drift_correction_params = calc_dc(pretrig_mean, filt_value, selection_good) when sum(selection_good)>2500
	@perpulse filt_value_dc = applycorrection(pretrigger_mean_drift_correction_params, pretrig_mean, filt_value)
	@threshold calibration = calibrate_nofit(filt_value_dc_hist, known_energies, calibration_nextra) when counted(filt_value_dc_hist)>2500
	@perpulse energy = apply_calibration(calibration, filt_value_dc)
	@histogram update_histogram!(filt_value_hist, selection_good, filt_value)
	@histogram update_histogram!(energy_hist, selection_good, energy)
	@histogram update_histogram!(filt_value_dc_hist, selection_good, filt_value_dc)
	ToJLDStep([:filt_value, :filt_value_dc, :energy, :filt_phase, :pretrig_rms, :postpeak_deriv, :rise_time, :peak_index,
	:pretrig_mean, :pulse_average, :pulse_rms, :peak_value, :min_value, :rowcount, :timestamp_posix_usec],
	Pair[:filter=>"filter/filter", :f_3db=>"filter/f_3db", :frametime=>"filter/frametime", :noise_autocorr=>"filter/noise_autocorr", :average_pulse=>"filter/average_pulse",
	:average_pulse=>"average_pulse",
	:samples_per_record=>"samples_per_record", :frametime=>"frametime", :pretrig_nsamples=>"pretrig_nsamples",
	:ljh_filename=>"ljh_filename", :noise_filename=>"noise_filename",
	:peak_index_criteria=>"selection_criteria/peak_index", :pretrig_rms_criteria=>"selection_criteria/pretrig_rms", :postpeak_deriv_criteria=>"selection_criteria/postpeak_deriv",
	:pretrigger_mean_drift_correction_params=>"pretrigger_mean_drift_correction_params"],
	mc[:hdf5_filename])
	MemoryLimitStep(Int(4e6)) # throw error if mc uses more than 4 MB
	FreeMemoryStep()
	] # end steps
	# write a verification function that makes sure all inputs either exist, or are the output of another step
	setsteps!(mc, steps)
	preptasks!(mc)
	mc
end

masschannels, watcher_exitchannel, watcher_task = schedule_MATTER_watcher(setup_channel,1)

rmf = ReferenceMicrocalFiles.dict["good_mnka_mystery"]

t = @schedule begin
getopenfilelimit() = parse(Int,split(split(readall(`ulimit -a`),"\n")[6])[end])
#if getopenfilelimit()>=1000
	sleep(2)
	LJHUtil.write_sentinel_file(rmf.noise_filename,false)
	sleep(2)
	LJHUtil.write_sentinel_file(rmf.filename,false)
	sleep(2)
	LJHUtil.write_sentinel_file(rmf.noise_filename,false)
	sleep(2)
	put!(watcher_exitchannel,1)
	@schedule begin
		tasks = [mc.task.value for (ch, mc) in masschannels]
		for task in tasks # wait for all tasks to finish
			try wait(task) end # wait rethrows errors from failed tasks, we can always look at them later
		end
		ndone = sum([task.state==:done for task in tasks])
		nfailed = sum([task.state==:failed for task in tasks])
		info("all tasks done!! Yay. $ndone :done and $nfailed :failed")
	end;
#else
#	println("open file limit is too low, try ulimit -n 1000")
#end
end # write to MATTER sentinel file to simulate matter writing various files

wait(t)
mc=masschannels[13];

wait(mc)
@test get(mc.task).state == :done
@test seen(mc[:energy_hist])==length(mc[:selection_good])
@test seen(mc[:filt_value_hist])==length(mc[:selection_good])
@test seen(mc[:filt_value_dc_hist])==length(mc[:selection_good])
@test 0.8*sum(mc[:selection_good])<counted(mc[:energy_hist])<sum(mc[:selection_good])
@test 0.8*sum(mc[:selection_good])<counted(mc[:filt_value_hist])<sum(mc[:selection_good])
@test 0.8*sum(mc[:selection_good])<counted(mc[:filt_value_dc_hist])<sum(mc[:selection_good])

h5open(mc[:hdf5_filename]) do h5
	@test read(h5["selection_criteria/peak_index"]) == mc[:peak_index_criteria]
	@test read(h5["selection_criteria/pretrig_rms"]) == mc[:pretrig_rms_criteria]
	@test read(h5["selection_criteria/postpeak_deriv"]) == mc[:postpeak_deriv_criteria]
	@test length(h5["filt_value"]) == length(mc[:filt_value])
end

# failing because I'm watching the LJH file for changes, and it doesn't change


text:
#=
Constants
=#

global OPEN_EPHYS = false

global EMG = false

function set_open_ephys()
    global OPEN_EPHYS = true
    nothing
end

function set_emg()
    global EMG = true
    nothing
end

#Constant parameters
@static if is_linux()
    base_path=string(dirname(Base.source_path()),"/../lib/")
    intan_lib = string(base_path,"libokFrontPanel_0118.so")
    #const lib = Libdl.dlopen(intan_lib,Libdl.RTLD_NOW)
end

@static if is_apple()
    base_path=string(dirname(Base.source_path()),"/../lib/")
    intan_lib = string(base_path,"libokFrontPanel.dylib")
end

@static if is_windows()
    base_path=string(dirname(Base.source_path()),"\\..\\lib\\")
    intan_lib = string(base_path,"okFrontPanel.dll")
end

const bit = string(base_path,"main.bit")
#const usb3bit = string(base_path,"XEM6310_512ch.bit")
const usb3bit_open_ephys = string(base_path,"rhd2000_usb3_0518.bit")
#const usb3bit_open_ephys = string(base_path,"rhd2000_usb3_oe.bit")
const usb3bit = string(base_path,"rhd2000_usb3.bit")

const USB_BUFFER_SIZE = 2400000
const RHYTHM_BOARD_ID = 500
const RHYTHM_BOARD_ID_OPENEPHYS = 600
const MAX_NUM_DATA_STREAMS = 16
const FIFO_CAPACITY_WORDS = 67108864

#const SAMPLES_PER_DATA_BLOCK = 600
#const SAMPLES_PER_DATA_BLOCK = 256
const SAMPLES_PER_DATA_BLOCK = 512*2
const USB3_BLOCK_SIZE = 1024

const DDR_BURST_LENGTH = 32

const RHD2000_HEADER_MAGIC_NUMBER_OPEN_EPHYS = 0xc691199927021942

const RHD2000_HEADER_MAGIC_NUMBER = 0xd7a22aaa38132a53

const WireInResetRun = 0x00

#Open Ephys
const WireInMaxTimeStepLsb = 0x01
const WireInMaxTimeStepMsb = 0x02

#Intan USB3
const WireInMaxTimeStep = 0x01
const WireInSerialDigitalInCntl = 0x02

const WireInDataFreqPll = 0x03
const WireInMisoDelay = 0x04
const WireInCmdRamAddr = 0x05
const WireInCmdRamBank = 0x06
const WireInCmdRamData = 0x07
const WireInAuxCmdBank1 = 0x08
const WireInAuxCmdBank2 = 0x09
const WireInAuxCmdBank3 = 0x0a

#OPEN EPHYS
const WireInAuxCmdLength1 = 0x0b
const WireInAuxCmdLength2 = 0x0c
const WireInAuxCmdLength3 = 0x0d
const WireInAuxCmdLoop1 = 0x0e
const WireInAuxCmdLoop2 = 0x0f #These conflict with TTL generation on Intan
const WireInAuxCmdLoop3 = 0x10 #Conflicts with Stimulation
const WireInLedDisplay_openephys = 0x11 #Conflicts with Stimulation
const WireInDataStreamSel1234 = 0x12 #conflicts with stimulation
const WireInDataStreamSel5678 = 0x13 #conflicts with stimulation

#Intan USB 3
const WireInAuxCmdLength = 0x0b
const WireInAuxCmdLoop = 0x0c
const WireInDacReref = 0x0e
const WireInLedDisplay = 0x0d

#Intan USB 3 Stimulation
const WireInStimCmdMode = 0x0f
const WireInStimRegAddr = 0x10
const WireInStimRegWord = 0x11

#Manual TTL triggers (0-15) and manual 16-bit dack value (31-16)
const WireInManualTriggers = 0x12

const WireInTtlOut = 0x13
const TrigInRamAddrReset = 0x43

const WireInDataStreamEn = 0x14
#const WireInTtlOut = 0x15
const WireInDacSource1 = 0x16
const WireInDacSource2 = 0x17
const WireInDacSource3 = 0x18
const WireInDacSource4 = 0x19
const WireInDacSource5 = 0x1a
const WireInDacSource6 = 0x1b
const WireInDacSource7 = 0x1c
const WireInDacSource8 = 0x1d
const WireInDacManual = 0x1e
const WireInMultiUse = 0x1f

const TrigInSpiStart = 0x41

#Intan USB3
const TrigInConfig = 0x40
const TrigInDacConfig = 0x42

#OpenEphys
const TrigInDcmProg = 0x40
const TrigInRamWrite = 0x42
const TrigInDacThresh = 0x43
const TrigInDacHpf = 0x44
const TrigInExtFastSettle = 0x45
const TrigInExtDigOut = 0x46
const TrigInOpenEphys = 0x5a

#Open Ephys
const WireOutNumWordsLsb = 0x20
const WireOutNumWordsMsb = 0x21

#Intan USB3
const WireOutNumWords = 0x20
const WireOutSerialDigitalIn = 0x21

const WireOutSpiRunning = 0x22
const WireOutTtlIn = 0x23
const WireOutDataClkLocked = 0x24
const WireOutBoardMode = 0x25
const WireOutBoardId = 0x3e
const WireOutBoardVersion = 0x3f

const PipeOutData = 0xa0

#For 32 channel amps
const PortA1 = 0
const PortA2 = 1
const PortB1 = 2
const PortB2 = 3
const PortC1 = 4
const PortC2 = 5
const PortD1 = 6
const PortD2 = 7

#For 64 channel amps
const PortA1Ddr = 8
const PortA2Ddr = 9
const PortB1Ddr = 10
const PortB2Ddr = 11
const PortC1Ddr = 12
const PortC2Ddr = 13
const PortD1Ddr = 14
const PortD2Ddr = 15

const speedOfLight = 299792458.0
const xilinxLvdsOutputDelay=1.9e-9
const xilinxLvdsInputDelay=1.4e-9
const rhd2000Delay=9.0e-9
const misoSettleTime=6.7e-9
const cableVelocity=.555*speedOfLight

# Adapted from https://github.com/JuliaAttic/Color.jl/issues/75#issuecomment-68073631
const jet_r,jet_g,jet_b=zeros(UInt32,256),zeros(UInt32,256),zeros(UInt32,256)
for i = 0:255
  n=4*i/256
  jet_r[i+1]=round(UInt32,255*min(max(min(n-1.5,-n+4.5),0),1));
  jet_g[i+1]=round(UInt32,255*min(max(min(n-0.5,-n+3.5),0),1));
  jet_b[i+1]=round(UInt32,255*min(max(min(n+0.5,-n+2.5),0),1));
end


text:
using JSON, Printf
import Base: merge, length

abstract type AbstractExtractor end
abstract type JSONEntry end
StringOrNumber = Union{AbstractString,Number}
max_keys = 10_000

"""
	updatemaxkeys!(n::Int)

limits the maximum number of keys in statistics of leaves in JSON. Default value is `10_000`.
"""
function updatemaxkeys!(n::Int)
	global max_keys = n
end

max_len = 10_000

"""
	updatemaxlen!(n::Int)

limits the maximum length of string values in statistics of nodes in JSON. Default value is `10_000`.
Longer strings will be trimmed and their length and hash will be appended to retain the uniqueness.
This is due to some strings being very long and causing the schema to be even order of magnitute larger than needed.
"""
function updatemaxlen!(n::Int)
	global max_len = n
end

function safe_update!(s::JSONEntry, d; path = "")
	success = update!(s, d)
	isnothing(success) && return nothing
	success && return s
	@info "In path $path: Instability in the schema detected. Using multiple representation."
	s = MultiEntry([s], s.updated)
	update!(s, d, path=path)
	return s
end
safe_update!(::Nothing, d; path = "") = newentry!(d)
update!(s::JSONEntry, d; path = "") = false

include("scalar.jl")
include("dict.jl")
include("array.jl")
include("multientry.jl")
include("makeschema.jl")

# todo: otestovat si schema pro inty a floaty v různém pořadí
# todo : pro multirepresentation mít inty a stringy
# todo: pro stringy to natvrdo přetypovat na abstractstringy, aby to pobralo šimonovu issue
# todo: otestovat když mám na stejném fieldu stringy a dictionary, co se stane když bude první dict nebo string, otestovat a popsat

updated(s::T) where {T<:JSONEntry} = s.updated
merge(combine::typeof(merge), es::JSONEntry...) = merge(es...)
merge(::Nothing, e::JSONEntry) = e

function Mill.reflectinmodel(sch::JSONEntry, ex::AbstractExtractor, fm=d->Flux.Dense(d, 10), fa=d->SegmentedMean(d); fsm = Dict(), fsa = Dict(),
			   single_key_identity=true, single_scalar_identity=true)
	# we do catobs of 3 samples here, because we want full representative sample to build whole "supersample" from which some samples are just subset
	# we also want sample with leaves missing in places where leaves can be actually missing in order to have imputation in all correct places
	full_sample = ex(sample_synthetic(sch, empty_dict_vals=false))
	leaves_empty_sample = ex(sample_synthetic(sch, empty_dict_vals=true))
	specimen = catobs(full_sample, leaves_empty_sample)
	reflectinmodel(specimen, fm, fa, fsm=fsm, fsa=fsa, single_key_identity=single_key_identity, single_scalar_identity=single_scalar_identity)
end

# this can probably be
make_selector(s) = e->select_by_string(s, e)
select_by_string(s::AbstractString, e::DictEntry) = e.childs[Symbol(s)]
select_by_string(s::AbstractString, e::ArrayEntry) = s == "[]" ? e.items : @error "wrong selector"
select_by_string(s::AbstractString, e::MultiEntry) = e.childs[parse(Int, s)]

"""
Deletes `field` at the specified `path` from the schema `sch`.
For instance, the following:
	`delete!(schema, ".field.subfield.[]", "x")`
deletes the field `x` from `schema` at:
	`schema.childs[:field].childs[:subfield].items.childs`
"""
function Base.delete!(sch::JSONEntry, path::AbstractString, field::AbstractString)
	@assert field != "[]"
	selectors = split(path, ".")[2:end]
	item = reduce((s, f) -> f(s), make_selector.(selectors), init=sch)
	delete!(item.childs, Symbol(field))
end

"""
	prune_json(json, schema)

Removes keys from `json` which are not part of the `schema`.

# Example
```jldoctest
julia> using JSON

julia> j1 = JSON.parse("{\\"a\\": 4, \\"b\\": {\\"a\\":1, \\"b\\": 1}}");

julia> j2 = JSON.parse("{\\"a\\": 4, \\"b\\": {\\"a\\":1}}");

julia> sch = JsonGrinder.schema([j1,j2])
[Dict] (updated = 2)
  ├── a: [Scalar - Int64], 1 unique values, updated = 2
  └── b: [Dict] (updated = 2)
           ├── a: [Scalar - Int64], 1 unique values, updated = 2
           └── b: [Scalar - Int64], 1 unique values, updated = 1

julia> j3 = Dict("a" => 4, "b" => Dict("a"=>1), "c" => 1, "d" => 2)
Dict{String,Any} with 4 entries:
  "c" => 1
  "b" => Dict("a"=>1)
  "a" => 4
  "d" => 2

julia> JsonGrinder.prune_json(j3, sch)
Dict{Any,Any} with 2 entries:
  "b" => Dict{Any,Any}("a"=>1)
  "a" => 4
```
so the `JsonGrinder.prune_json` removes keys `c` and `d`.
"""
prune_json(json, sch::Entry) = json

prune_json(json, sch::ArrayEntry) = map(json) do el
	prune_json(el, sch.items)
end

function prune_json(json, sch::DictEntry)
    out = Dict()
    for (k,v) in children(sch)
        String(k) ∈ keys(json) && (out[String(k)] = prune_json(json[String(k)], v))
    end
    out
end


text:
using Flux
using Random
using Distributions

mutable struct ActorCriticNetwork{M, Π, V, D}
    base::M
    π::Π
    v::V
    dist::D
end

mutable struct ActorCriticAgent{M, TN, O, LU, AP<:AbstractValuePolicy, Φ, ERP} <: AbstractAgent
    network::M
    target_network::TN
    opt::O
    lu::LU
    ap::AP
    er::ExperienceReplay{ERP}
    γ::Float32
    batch_size::Int64
    tn_counter_init::Int64
    target_network_counter::Int64
    action::Int64
    prev_s::Φ
end

function ActorCriticAgent(model, target_network, opt, lu, ap, size_buffer, γ, batch_size, tn_counter_init, s)

    model = ActorCriticNetwork(
        Chain(Dense(2, 128, Flux.relu; initW=Flux.glorot_uniform)
              Dense(128, 128, Flux.relu; initW=Flux.glorot_uniform),
              Dense(128, 32, Flux.relu; initW=Flux.glorot_uniform)),
        [Dense(32, 1), Dense(32, 1, Flux.relu)],
        Dense(32, 1),
        Normal)
    # model = Chain(Dense(length(s), 128, Flux.relu; initW=Flux.glorot_uniform),
    #               Dense(128, 128, Flux.relu; initW=Flux.glorot_uniform),
    #               Dense(128, 32, Flux.relu; initW=Flux.glorot_uniform),
    #               Dense(32, num_actions; initW=Flux.glorot_uniform))

    ActorCriticAgent(model,
                     target_network,
                     opt,
                     lu,
                     ap,
                     ExperienceReplay(100000,
                                      (Array{Float32, 1}, Int64, Array{Float32, 1}, Float32, Bool),
                                      (:s, :a, :sp, :r, :t)),
                     γ,
                     batch_size,
                     tn_counter_init,
                     tn_counter_init,
                     0,
                     s)
end

function JuliaRL.start!(agent::ActorCriticAgent, env_s_tp1, rng::AbstractRNG=Random.GLOBAL_RNG; kwargs...)
    # Start an Episode
    agent.action = sample(agent.ap,
                          agent.model(env_s_tp1),
                          rng)
    return agent.action
end

function JuliaRL.step!(agent::ActorCriticAgent, env_s_tp1, r, terminal, rng=Random.GLOBAL_RNG; kwargs...)

    add!(agent.er, (Float32.(agent.prev_s), agent.action, copy(Float32.(env_s_tp1)), r, terminal))

    if size(agent.er)[1] > 1000
        e = sample(agent.er, agent.batch_size; rng=rng)
        update_params!(agent, e)
    end

    agent.prev_s .= env_s_tp1
    agent.action = sample(agent.ap,
                          agent.model(agent.prev_s),
                          rng)

    return agent.action
end

function update_params!(agent::ActorCriticAgent, e)

    # update!(agent.model, agent.lu, agent.opt, e.s, e.a, e.sp, e.r, e.t, agent.target_network)

    s_t = hcat(e.s...)
    s_tp1 = hcat(e.sp...)

    b_tp1 = Flux.data(base(s_tp1))
    v_tp1 = Flux.data(agent.model.v(b_tp1))
    target = e.r .+ agent.γ.*(.!e.t).*v_tp1

    # update value function

    gs = Flux.Tracker.Gradient(params(agent.model.base, agent.model.π, agent.model.v)) do
        b_t = base(s_t)
        δ = target .- agent.model.v(b_t)
        μ, σ = agent.model.π[1](b_t)
    end

    # update policy

    if agent.target_network_counter == 1
        agent.target_network_counter = agent.tn_counter_init
        agent.target_network = mapleaves(
            Flux.Tracker.data,
            deepcopy(agent.model))
    else
        agent.target_network_counter -= 1
    end

    return nothing

end


text:
# Autogenerated wrapper script for ngspice_jll for x86_64-apple-darwin
export libngspice, ngspice

JLLWrappers.@generate_wrapper_header("ngspice")
JLLWrappers.@declare_library_product(libngspice, "@rpath/libngspice.0.dylib")
JLLWrappers.@declare_executable_product(ngspice)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libngspice,
        "lib/libngspice.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        ngspice,
        "bin/ngspice",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()