text:
"""
Abstract supertype for organ parameters.
Extend to add additional components to organ parameters.
"""
abstract type AbstractParams end

"""
Model parameters that vary between organs
"""
@default_kw @flattenable @selectable struct Params{As,Sh,Al,Ma,Tr,Re,Ge,Pr} <: AbstractParams
    # Field               | Default                    | _     | Selectable Types
    assimilation_pars::As | ConstantCAssim()           | _     | Union{Nothing,AbstractAssim}
    scaling_pars::Sh      | nothing                    | _     | Union{Nothing,AbstractScaling}
    allometry_pars::Al    | nothing                    | _     | Union{Nothing,AbstractAllometry}
    maturity_pars::Ma     | nothing                    | _     | Union{Nothing,AbstractMaturity}
    activetrans_pars::Tr  | nothing                    | _     | Union{Nothing,ActiveTranslocation}
    passivetrans_pars::Re | LosslessPassiveTranslocation() | _ | PassiveTranslocation
    germination_pars::Ge  | nothing                    | _     | Union{Nothing,AbstractGermination}
    production_pars::Pr   | nothing                    | _     | Union{Nothing,AbstractProduction}
end

for fn in fieldnames(Params)
    @eval $fn(p::Params) = p.$fn
end

"""
Astract supertype for shared parameters.
Extend to change the components that are shared.
"""
abstract type AbstractSharedParams end

"""
    SharedParams(su_pars, core_pars, resorption_pars, tempcorr_pars, catabolism_pars)

Model parameters shared between organs.

# FieldMetadata macros
- `@default_kw` provides default values and constructors,
- `@selectable` provides the set of types that can be used for the parameter,
  so that they can be selected in a live interface.
"""
@default_kw @selectable struct SharedParams{SU,Co,Fe,Te,Ca} <: AbstractSharedParams
    # Field               | Default                   | Selectable Types
    su_pars::SU           | ParallelComplementarySU() | AbstractSynthesizingUnit
    core_pars::Co         | DEBCore()                 | _
    resorption_pars::Fe   | nothing                   | Union{Nothing,AbstractResorption}
    tempcorr_pars::Te     | nothing                   | Union{Nothing,AbstractTemperatureCorrection}
    catabolism_pars::Ca   | CatabolismCN()            | AbstractCatabolism
end

for fn in fieldnames(SharedParams)
    @eval $fn(p::SharedParams) = p.$fn
end

###########################################################################################
# Variables

"""
Model variables.
Allow storing and accessing variables for use by multiple components.
"""
abstract type AbstractVars end

"""
    PlottableVars()

Plottable model variables. These are vectors witih values for each time-step,
to allow plotting and model introspection.
"""
@udefault_kw @units @plottable struct PlottableVars{F,R,E,T,WP,H,TS} <: AbstractVars
    scaling::F        | [0.0] | _                 | _
    rate::R           | [0.0] | mol*mol^-1*d^-1   | _
    E_ctb::E          | [0.0] | mol/hr            | _
    θE::F             | [0.0] | _                 | _
    temp::T           | [0.0] | K                 | _
    tempcorrection::F | [0.0] | _                 | _
    swp::WP           | [0.0] | kPa               | _
    soilcorrection::F | [0.0] | _                 | _
    height::H         | [0.0] | m                 | _
    tstep::TS         | [1]   | _                 | false
end

"""
    Vars()

Mutable struct to allow storing variables
for use by multiple components.
"""
@udefault_kw @units mutable struct Vars{F,R,E,T,WP,H} <: AbstractVars
    scaling::F        | 0.0   | _
    rate::R           | 0.0   | mol*mol^-1*d^-1
    θE::F             | 0.0   | _
    E_ctb::E          | 0.0   | mol/hr
    temp::T           | 0.0   | K
    tempcorrection::F | 0.0   | _
    swp::WP           | 0.0   | kPa
    soilcorrection::F | 0.0   | _
    height::H         | 0.0   | m
end

tstep(v::PlottableVars) = v.tstep[1]
tstep(v::Vars) = nothing

set_tstep!(v::PlottableVars, val) = v.tstep[1] = val
set_tstep!(v::Vars, val) = nothing

depth(v) = height(v)

build_vars(vars::Vars, tspan) = vars
build_vars(vars::PlottableVars, tspan) = begin
    len = length(tspan)
    len <= length(vars.rate) && return vars

    for fname in fieldnames(typeof(vars))
        varvec = getfield(vars, fname)
        append!(varvec, fill(getfield(vars, fname)[1], len - length(varvec)))
    end
    vars
end

###########################################################################################
# Organs and Organisms

"""
Abstract supertype for organs. Inherit from it if you need to difine
behaviour different to that or [`Organ`](@ref).
"""
abstract type AbstractOrgan{P,S} end

params(o::AbstractOrgan) = o.params
shared(o::AbstractOrgan) = o.shared
vars(o::AbstractOrgan) = o.vars
flux(o::AbstractOrgan) = o.J

#= κ functions allow modular addition of destinations for catabolised flux,
by default maturity and active translocation. If those components are
not used the flux fraction will be allocated to κsoma. =#

κtra(o::AbstractOrgan) = κtra(activetrans_pars(o))
κtra(o::Nothing) = 0.0

κmat(o::AbstractOrgan) = κmat(maturity_pars(o))
κmat(::Nothing) = 0.0

κsoma(o::AbstractOrgan) = oneunit(κtra(o)) - κtra(o) - κmat(o)

# Define `scaling` and `setscaling` etc. methods
for fn in fieldnames(Vars)
    setfn = Symbol.(:set_, fn, :!)
    @eval $fn(o::AbstractOrgan) = $fn(vars(o))
    @eval $setfn(o::AbstractOrgan, x) = $setfn(vars(o), x)
    @eval @inline ($fn)(vars::PlottableVars) = vars.$fn[tstep(vars)]
    @eval @inline ($setfn)(vars::PlottableVars, val) = vars.$fn[tstep(vars)] = val
    @eval @inline ($fn)(vars::Vars) = vars.$fn
    @eval @inline ($setfn)(vars::Vars, val) = vars.$fn = val
end

#= Forward variable and parameter getter/setter methods
so they can be acessesed/set directly from the organ,
without knowning where they are actually stored. =#
tstep(o::AbstractOrgan) = tstep(vars(o))
set_tstep!(o::AbstractOrgan, t) = set_tstep!(vars(o), t)

for fn in fieldnames(DEBCore)
    @eval $fn(p::SharedParams) = $fn(core_pars(p))
    @eval $fn(o::AbstractOrgan{<:Any,<:SharedParams}) = $fn(shared(o))
end

for fn in fieldnames(SharedParams)
    @eval $fn(o::AbstractOrgan{<:Any,<:SharedParams}) = $fn(shared(o))
end

for fn in fieldnames(Params)
    @eval $fn(o::AbstractOrgan{<:Params}) = $fn(params(o))
end

"""
    Organ(params, shared, vars, J)

Basic model components. For a plants, organs might be roots, stem and leaves
"""
struct Organ{P,S,V,F} <: AbstractOrgan{P,S}
    params::P
    shared::S
    vars::V
    J::F
end
"""
    Organ(params, shared, records)

Construct an organ from parameters, shared parameters and
views into records arrays for vaiable and flux matrices.
"""
Organ(params::AbstractParams, shared::AbstractSharedParams, records) = begin
    vars = records.vars
    set_tstep!(vars, 1)
    J = view(records.J, Ti(1))
    Organ(params, shared, vars, J)
end

# Records hold vars and flux, allowing storage of each
# timestep for plotting if PlottableVars are used.
abstract type AbstractRecords end

@plottable struct PlottableRecords{V,F} <: AbstractRecords
    vars::V | true
    J::F    | false
end
PlottableRecords(vars::PlottableVars, J::AbstractArray) =
    PlottableRecords{map(typeof,(vars,J))...}(vars, J)
PlottableRecords(vars::PlottableVars, fluxval::Number, fluxaxes::Tuple, tspan::AbstractRange) = begin
    vars = build_vars(vars, tspan)
    J = build_flux(fluxval, fluxaxes..., tspan)
    PlottableRecords(vars, J)
end

@plottable struct Records{V,F} <: AbstractRecords
    vars::V | false
    J::F    | false
end
Records(vars::Vars, J::AbstractArray) =
    Records{map(typeof,(vars,J))...}(vars, J)
Records(vars::Vars, fluxval::Number, fluxaxes::Tuple) = begin
    J = build_flux(fluxval, fluxaxes...)
    Records{map(typeof, (vars,J))...}(vars, J)
end

build_flux(fluxval, x::Tuple, y::Tuple) = begin
    dims = X(Val(x)), Y(Val(y))
    A = zeros(typeof(fluxval), map(length, dims)...)
    DimensionalArray(A, dims)
end
build_flux(fluxval, x::Tuple, y::Tuple, time::AbstractRange) = begin
    dims = X(Val(x)), Y(Val(y)), Ti(time)
    A = zeros(typeof(fluxval), map(length, dims)...)
    DimensionalArray(A, dims)
end

"""
A an Organism is a model object for modelling the growth
of a single organism.

It can be run as a functor:

```julia
(o::AbstactOrganism)(du, u, p, t) =
...
```

So that it can be passed into DifferentialEquations.jl solvers
if required.

Where `du` is the flux to be updated, `u` is the current state,
`p` are new model paremeters or `nothing`, and `t` is the current
time step. See model.jl for implemntations.

When julia 1.5 is released this will be implemented for any
`AbstactOrganism`, but for now it is only implemented for `Plant`
due to current limitations in Julia.
"""
abstract type AbstractOrganism end

params(o::AbstractOrganism) = o.params
shared(o::AbstractOrganism) = o.shared
records(o::AbstractOrganism) = o.records
organs(o::AbstractOrganism) = o.organs
environment(o::AbstractOrganism) = o.environment
environment_start(o::AbstractOrganism) = o.environment_start
dead(o::AbstractOrganism) = o.dead[]
set_dead!(o::AbstractOrganism, val) = o.dead[] = val

"""
    define_organs(o::AbstractOrganism, t)

Organs are constructed with views of Records and J Arrays at time t
"""
define_organs(o::AbstractOrganism, t) =
    define_organs(params(o), shared(o), records(o), t)
define_organs(params::Tuple, shared, records::Tuple, t) =
    map((p, r) -> Organ(p, shared, r), params, records)

update_organs(organs, t) = update_organs(vars(first(organs)), organs, t)
update_organs(::Vars, organs, t) = organs
update_organs(::PlottableVars, organs, t) = begin
    i = floor(Int, ustrip(t))
    map(organs) do organ
        organ.vars.tstep[1] = i
        J = organ.J
        data = J.data
        @set! data.indices[3] = i
        @set! data.offset1 = (i - 1) * length(data.indices[1]) * length(data.indices[2])
        # Setfield.jl isn't type-stable, DD `rebuild` is better
        @set! organ.J = DimensionalData.rebuild(J, data)
        organ
    end
end

"""
    Plant(params, shared, records, environment, environment_start, dead)
    Plant(states=(:V, :C, :N),
          transformations=(:asi, :gro, :mai, :rej, :res),
          params=(ShootParamsCN(), RootParamsCN()),
          vars=(Vars(), Vars()),
          shared=SharedParams(),
          records=nothing,
          environment=nothing,
          time=0.0hr:1.0hr:8760.0hr,
          environment_start=Ref(1.0hr),
          dead=Ref(false))

Plant model.

`params` are a tuple of Parameters objects, one for each organ.

`shared` is a struct of parameters shared between organs.

`vars`: can be `Vars` or `PlottableVars` or custom struct with additional variables.
  `PlottableVars` will be stored for each timestep for plotting.

`environment` can be `nothing`, or a `MicroclimPoint` or `MicroclimControl`
from Microclimates.jl

`time` determines the timespan over which `PlottableVars` will be constructed.
it isn't used with regular `Vars`.

If 4 state model is used, pass in
```julia
states=(:V, :E, :C, :N),
```

Similarly, if additional components like Maturity or ActiveTranslocation
are used, their state label needs to be passed in:

```julia
transformations=(:asi, :gro, :mai, :mat, :rej, :tra, :res),
```
"""
@flattenable @description mutable struct Plant{P,S,R,O,E,ES,D} <: AbstractOrganism
    params::P             | true  | "Model parameters"
    shared::S             | true  | "Parameters shared between organs"
    records::R            | false | "Plotable variables stored in arrays and sliced for each timestep on demand"
    organs::O             | false | ""
    environment::E        | false | "Environment object, provides environmental variables for each timestep"
    environment_start::ES | false | "Start index of environmental data"
    dead::D               | false | "`Bool` flag: has the plant died"
    Plant(params::P, shared::S, records::R, organs::O, environment::E, environment_start::ES, dead::D) where {P,S,R,O,E,ES,D} = begin
        organs = define_organs(params, shared, records, 0)
        new{P,S,R,typeof(organs),E,ES,D}(params, shared, records, organs, environment, environment_start, dead)
    end
end
Plant(; states=(:V, :C, :N),
        transformations=(:asi, :gro, :mai, :rej, :tra, :res),
        params=(
            Params(assimilation_pars=KooijmanSLAPhotosynthesis()),
            Params(assimilation_pars=NAssim()),
        ),
        vars=(Vars(), Vars()),
        shared=SharedParams(),
        environment=nothing,
        time=0.0hr:1.0hr:8760.0hr,
        records=nothing,
        environment_start=Ref(0.0hr),
        dead=Ref(false)
      ) = begin
    fluxaxes = states, transformations
    fluxval = 1.0mol/hr
    records = build_records(records, vars, fluxval, fluxaxes, time)
    Plant(params, shared, records, nothing, environment, environment_start, dead)
end

build_records(records::AbstractRecords, args...) = records
build_records(records::Nothing, vars::Tuple, fluxval, fluxaxes, tspan) =
    map(vars) do v
        build_records(v, fluxval, fluxaxes, tspan)
    end
build_records(vars::Vars, fluxval::Number, fluxaxes::Tuple, tspan::AbstractRange) =
    Records(vars, fluxval, fluxaxes)
build_records(vars::PlottableVars, fluxval::Number, fluxaxes::Tuple, tspan::AbstractRange) =
    PlottableRecords(vars, fluxval, fluxaxes, tspan)

# Define a SubArray constructor so we can increment the time index
ConstructionBase.constructorof(::Type{A}) where A<:SubArray{T,N,P,I,L} where {T,N,P,I,L} =
    (parent, indices, offset1, stride1) -> begin
        SubArray{eltype(parent),ndims(parent),typeof(parent),typeof(indices),L}(parent, indices, offset1, stride1)
    end


text:
# Script for running DA and ADA using the same traninig data

# go to Ricker model folder
try
  cd("Ricker model")
catch
 warn("Already in the Ricker model folder")
end

# load algorithms
include("rickermodel.jl")

# packages for storing data
using JLD
using HDF5

################################################################################
###      set up problem                                                      ###
################################################################################

problem = set_up_gp_problem(ploton = false)

# set adaptive updating
#problem.adaptive_update = AMUpdate(eye(3), 2.4/sqrt(3), 1., 0.7, 50)

#problem.adaptive_update = noAdaptation(2.4/sqrt(3)*eye(3))

# or, use AM gen alg for adaptive updating
#problem.adaptive_update = AMUpdate_gen(eye(3), 2.4/sqrt(3), 0.2, 1., 0.8, 25)
problem.adaptive_update = AMUpdate_gen(eye(3), 2.4/sqrt(3), 0.3, 1., 0.8, 25)

# set algorithm parameters
problem.alg_param.N = 1000 # nbr particels
problem.alg_param.R = 50000 # nbr iterations
problem.alg_param.burn_in = 0 # burn in
problem.alg_param.length_training_data = 2000
problem.alg_param.alg = "MCWM" # we should only! use the MCWM algorithm
problem.alg_param.compare_GP_and_PF = false
problem.alg_param.noisy_est = false
problem.alg_param.pred_method = "sample"
problem.alg_param.print_interval = 10000 # problem.alg_param.R#
problem.alg_param.beta_MH = 0.15 # "local_loglik_approx" # "max_loglik"
problem.alg_param.lasso = false

#problem.data.y = Array(readtable("y.csv"))[:,1]
#problem.data.y = Array(readtable("y_data_set_1.csv"))[:,1]
problem.data.y = Array(readtable("y_data_set_2.csv"))[:,1]

################################################################################
###      generate traning data                                               ###
################################################################################

# set up training problem

#accelerated_da = true

problem_training = set_up_problem(ploton = false)

length_training_data = 2000
length_test_data = 2000
burn_in = 2000

problem_training.alg_param.N = 1000 # nbr particels
problem_training.alg_param.R = length_training_data + length_test_data + burn_in # nbr iterations
problem_training.alg_param.burn_in = burn_in # burn_in
problem_training.data.y = Array(readtable("y_data_set_2.csv"))[:,1] #Array(readtable("y.csv"))[:,1]
problem_training.alg_param.print_interval = 1000

# test starting at true parameters
#problem.model_param.theta_0 = problem.model_param.theta_true

# PMCMC
problem_training.alg_param.alg = "MCWM"

# use AM alg for adaptive updating
#problem.adaptive_update = AMUpdate(eye(3), 2.4/sqrt(3), 1., 0.7, 25)

#problem.adaptive_update = noAdaptation(2.4/sqrt(3)*eye(3))

# or, use AM gen alg for adaptive updating
#problem_training.adaptive_update = AMUpdate_gen(eye(3), 2.4/sqrt(3), 0.2, 1., 0.8, 25)
problem_training.adaptive_update = AMUpdate_gen(eye(3), 2.4/sqrt(3), 0.4, 1., 0.8, 25)

load_training_data = true

if !load_training_data

  # generate training data
  tic()
  # collect data
  res_training, Theta_star_training, loglik_star_training,Theta_old_training,loglik_old_training, cov_matrix = mcmc(problem_training, true, true)

  time_pre_er = toc()

  # write outputs
  res = res_training[1]

  Theta = res.Theta_est
  loglik = res.loglik_est
  accept_vec = res.accept_vec
  prior_vec = res.prior_vec

  loglik_avec_priorvec = zeros(3, length(loglik))
  loglik_avec_priorvec[1,:] = loglik
  loglik_avec_priorvec[2,:] = accept_vec
  loglik_avec_priorvec[3,:] = prior_vec

  algorithm_parameters = zeros(10, 2)

  algorithm_parameters[1,1] = problem_training.alg_param.burn_in
  algorithm_parameters[2:4,1] = problem_training.model_param.theta_true
  algorithm_parameters[5:7,1] = problem_training.model_param.theta_0
  algorithm_parameters[8:end,:] = problem_training.prior_dist.prior_parameters

  writetable("Results/Theta_training.csv", convert(DataFrame, Theta))
  writetable("Results/loglik_avec_priorvec_training.csv", convert(DataFrame, loglik_avec_priorvec))
  writetable("Results/algorithm_parameters_training.csv", convert(DataFrame, algorithm_parameters))

  # split tranining and test data

  Theta_test_star = Theta_star_training[:,(end-length_test_data+1):end]
  loglik_test_star = loglik_star_training[(end-length_test_data+1):end]

  Theta_test_old = Theta_old_training[:,(end-length_test_data+1):end]
  loglik_test_old = loglik_old_training[(end-length_test_data+1):end]

  data_test_star = [Theta_test_star; loglik_test_star']
  data_test_old = [Theta_test_old; loglik_test_old']

  Theta_training_star = Theta_star_training[:,1:length_training_data]
  loglik_training_star = loglik_star_training[1:length_training_data]

  Theta_training_old = Theta_old_training[:,1:length_training_data]
  loglik_training_old = loglik_old_training[1:length_training_data]

  data_training_star = [Theta_training_star; loglik_training_star']
  data_training_old = [Theta_training_old; loglik_training_old']

  save("gp_training_and_test_data.jld",
        "res_training", res_training,
        "data_training_star", data_training_star,
        "data_training_old", data_training_old,
        "data_test_star", data_test_star,
        "data_test_old", data_test_old,
        "cov_matrix",cov_matrix)

else

  @load "gp_training_and_test_data.jld"

end

################################################################################
###     fit gp model                                                         ###
################################################################################

# fit gp model
tic()

# create gp object
gp = GPModel("est_method",zeros(6), zeros(4),
eye(problem.alg_param.length_training_data-20), zeros(problem.alg_param.length_training_data-20),zeros(2,problem.alg_param.length_training_data-20),
collect(1:10))

data_training = data_training_star
data_test = data_test_star

# fit GP model
if problem.alg_param.est_method == "ml"
  # fit GP model using ml

  perc_outlier = 0.1 # used when using PMCMC for trainig data 0.05
  tail_rm = "left"
  problem.alg_param.lasso = false

  ml_est(gp, data_training,"SE", problem.alg_param.lasso,perc_outlier,tail_rm)
else
  error("The two stage estimation method is not in use")
  #two_stage_est(gp, data_training)
end

time_fit_gp = toc()

# save fitted gp model
# @save "gp_fitted_model.jld" gp
################################################################################
###     DA-GP-MCMC                                                           ###
################################################################################

accelerated_da = false

problem.model_param.theta_0 = mean(res_training[1].Theta_est[:,problem_training.alg_param.burn_in+1:problem_training.alg_param.burn_in+length_training_data],2)

res = @time dagpmcmc(problem_training, problem, gp, cov_matrix)

# calc res
mcmc_results = Result(res[1].Theta_est, res[1].loglik_est, res[1].accept_vec, res[1].prior_vec)

# write output
Theta = mcmc_results.Theta_est
loglik = mcmc_results.loglik_est
accept_vec = mcmc_results.accept_vec
prior_vec = mcmc_results.prior_vec

loglik_avec_priorvec = zeros(3, length(loglik))
loglik_avec_priorvec[1,:] = loglik
loglik_avec_priorvec[2,:] = accept_vec
loglik_avec_priorvec[3,:] = prior_vec

algorithm_parameters = zeros(10, 2)

algorithm_parameters[1,1] = problem.alg_param.burn_in + 1
algorithm_parameters[2:4,1] = problem.model_param.theta_true
algorithm_parameters[5:7,1] = problem.model_param.theta_0
algorithm_parameters[8:end,:] = problem.prior_dist.prior_parameters

if !accelerated_da
  writetable("Results/Theta_dagpmcmc_lunarc.csv", convert(DataFrame, Theta))
  writetable("Results/loglik_avec_priorvec_dagpmcmc_lunarc.csv", convert(DataFrame, loglik_avec_priorvec))
  writetable("Results/algorithm_parameters_dagpmcmc_lunarc.csv", convert(DataFrame, algorithm_parameters))
else
  writetable("Results/Theta_adagpmcmc_lunarc.csv", convert(DataFrame, Theta))
  writetable("Results/loglik_avec_priorvec_adagpmcmc_lunarc.csv", convert(DataFrame, loglik_avec_priorvec))
  writetable("Results/algorithm_parameters_adagpmcmc_lunarc.csv", convert(DataFrame, algorithm_parameters))
end

################################################################################
###     A-DA-GP-MCMC                                                         ###
################################################################################

accelerated_da = true

problem.model_param.theta_0 = mean(res_training[1].Theta_est[:,problem_training.alg_param.burn_in+1:problem_training.alg_param.burn_in+length_training_data],2)

################################################################################
##  Create features for classification models                                                                            ##
################################################################################

n = size(data_training,2)
n_burn_in = problem_training.alg_param.burn_in

dim = length(problem.model_param.theta_true)
data_gp_loglik_star_old = zeros(4,n)
data_gp_loglik_star_old[3,:] = data_training_star[dim+1,:]
data_gp_loglik_star_old[4,:] = data_training_old[dim+1,:]

# data_gp_loglik_star_old = [gp_star, gp_old, ll_star, ll_old]

std_pred_gp_star = zeros(n)

noisy_pred = problem.alg_param.noisy_est

for i = 1:n
  (loglik_est_star, var_pred_ml_star, prediction_sample_ml_star) = predict(data_training_star[1:dim,i],gp,noisy_pred)
  (loglik_est_old, var_pred_ml, prediction_sample_ml_old) = predict(data_training_old[1:dim,i],gp,noisy_pred)
  data_gp_loglik_star_old[1,i] = prediction_sample_ml_star[1]
  data_gp_loglik_star_old[2,i] = prediction_sample_ml_old[1]
  std_pred_gp_star[i] = sqrt(var_pred_ml_star[1])
end

nbr_GP_star_geq_GP_old = zero(Int64)
nbr_case_1 = zero(Int64)
nbr_case_2 = zero(Int64)

targets_case_1_and_3 = []
data_case_1_and_3 = []

targets_case_2_and_4 = []
data_case_2_and_4 = []

for i = 1:n
  if data_gp_loglik_star_old[1,i] > data_gp_loglik_star_old[2,i]
    nbr_GP_star_geq_GP_old += 1
    data_case_1_and_3 = vcat(data_case_1_and_3, [data_training[1:dim,i]; data_gp_loglik_star_old[1,i]/data_gp_loglik_star_old[2,i]; std_pred_gp_star[i]])
    if data_gp_loglik_star_old[3,i] > data_gp_loglik_star_old[4,i]
      append!(targets_case_1_and_3, 1)
      nbr_case_1 += 1
    else
      append!(targets_case_1_and_3, 0)
    end
  else
    data_case_2_and_4 = vcat(data_case_2_and_4, [data_training[1:dim,i]; data_gp_loglik_star_old[1,i]/data_gp_loglik_star_old[2,i]; std_pred_gp_star[i]])
    if data_gp_loglik_star_old[3,i] < data_gp_loglik_star_old[4,i]
      append!(targets_case_2_and_4, 1)
      nbr_case_2 += 1
    else
      append!(targets_case_2_and_4, 0)
    end
  end
end

# tansform features and set input data

# convert matricies to floats
data_case_1_and_3 = convert(Array{Float64,2},reshape(data_case_1_and_3, (dim+2, length(targets_case_1_and_3))))
data_case_2_and_4 = convert(Array{Float64,2},reshape(data_case_2_and_4, (dim+2, length(targets_case_2_and_4))))

targets_case_1_and_3 = convert(Array{Float64,1}, targets_case_1_and_3)
targets_case_2_and_4 = convert(Array{Float64,1}, targets_case_2_and_4)

################################################################################
##   set case model                                                          ###
################################################################################

select_case_model = "dt" # logisticregression or dt

nbr_GP_star_led_GP_old = n-nbr_GP_star_geq_GP_old

prob_case_1 = nbr_case_1/nbr_GP_star_geq_GP_old
prob_case_2 = (nbr_case_2)/nbr_GP_star_led_GP_old
prob_case_3 = 1-prob_case_1
prob_case_4 = (nbr_GP_star_led_GP_old-nbr_case_2)/nbr_GP_star_led_GP_old
prob_cases = [prob_case_1;prob_case_2;prob_case_3;prob_case_4]

println("Est prob:")
println(prob_cases)
#prob_cases = [0.2;0.2;0.8;0.8]

if select_case_model == "biasedcoin"

  casemodel = BiasedCoin(prob_cases)

elseif select_case_model == "logisticregression"

  mean_posterior = mean(theta_training,2)[:]

  input_data_case_1_and_3 = zeros(length(targets_case_1_and_3), dim+3)
  input_data_case_1_and_3[:,1] = sqrt((mean_posterior[1] - data_case_1_and_3[1,:]).^2)
  input_data_case_1_and_3[:,2] = sqrt((mean_posterior[2] - data_case_1_and_3[2,:]).^2)
  input_data_case_1_and_3[:,3] = sqrt((mean_posterior[3] - data_case_1_and_3[3,:]).^2)
  input_data_case_1_and_3[:,4] = sqrt(sum((repmat(mean_posterior', size(data_case_1_and_3,2))'-data_case_1_and_3[1:3,:]).^2,1))
  input_data_case_1_and_3[:,5] = data_case_1_and_3[4,:]
  input_data_case_1_and_3[:,end] = targets_case_1_and_3

  input_data_case_1_and_3 = DataFrame(input_data_case_1_and_3)

  input_data_case_2_and_4 = zeros(length(targets_case_2_and_4), dim+3)
  input_data_case_2_and_4[:,1] = sqrt((mean_posterior[1] - data_case_2_and_4[1,:]).^2)
  input_data_case_2_and_4[:,2] = sqrt((mean_posterior[2] - data_case_2_and_4[2,:]).^2)
  input_data_case_2_and_4[:,3] = sqrt((mean_posterior[3] - data_case_2_and_4[3,:]).^2)
  input_data_case_2_and_4[:,4] = sqrt(sum((repmat(mean_posterior', size(data_case_2_and_4,2))'-data_case_2_and_4[1:3,:]).^2,1))
  input_data_case_2_and_4[:,5] = data_case_2_and_4[4,:]
  input_data_case_2_and_4[:,end] = targets_case_2_and_4

  input_data_case_2_and_4 = DataFrame(input_data_case_2_and_4)

  log_reg_model_case_1_and_3 = glm(@formula(x6 ~ x1 + x2 + x3 + x4 + x5), input_data_case_1_and_3, Binomial(), LogitLink())

  log_reg_model_case_2_and_4 = glm(@formula(x6 ~ x1 + x2 + x3 + x4 + x5), input_data_case_2_and_4, Binomial(), LogitLink())

  β_for_model1or3 = coef(log_reg_model_case_1_and_3)
  β_for_model2or4 = coef(log_reg_model_case_2_and_4)

  casemodel = LogisticRegression(β_for_model1or3, β_for_model2or4, mean_posterior)

elseif select_case_model == "dt"

  #standardization!(data_case_1_and_3)

  input_data_case_1_and_3 = zeros(length(targets_case_1_and_3), dim+2)

  input_data_case_1_and_3[:,1] =  data_case_1_and_3[1,:]
  input_data_case_1_and_3[:,2] =  data_case_1_and_3[2,:]
  input_data_case_1_and_3[:,3] =  data_case_1_and_3[3,:]
  input_data_case_1_and_3[:,4] =  data_case_1_and_3[4,:]
  input_data_case_1_and_3[:,end] = targets_case_1_and_3

  #standardization!(data_case_2_and_4)

  input_data_case_2_and_4 = zeros(length(targets_case_2_and_4), dim+2)

  input_data_case_2_and_4[:,1] =  data_case_2_and_4[1,:]
  input_data_case_2_and_4[:,2] =  data_case_2_and_4[2,:]
  input_data_case_2_and_4[:,3] =  data_case_2_and_4[3,:]
  input_data_case_2_and_4[:,4] =  data_case_2_and_4[4,:]
  input_data_case_2_and_4[:,end] = targets_case_2_and_4

  # 4 features model
  features_1_and_3 = convert(Array, input_data_case_1_and_3[:, 1:dim+1])

  labels_1_and_3 = convert(Array, input_data_case_1_and_3[:, end])

  labels_1_and_3 = Array{String}(size(features_1_and_3,1))

  for i = 1:length(labels_1_and_3)
    if input_data_case_1_and_3[i,end] == 0
      labels_1_and_3[i] = "case 3"
    else
      labels_1_and_3[i] = "case 1"
    end
  end

  decisiontree1or3 = build_tree(labels_1_and_3, features_1_and_3)

  decisiontree1or3 = prune_tree(decisiontree1or3, 0.9)

  # 4 features model
  features_case_2_and_4 = convert(Array, input_data_case_2_and_4[:, 1:dim+1])

  labels_case_2_and_4 = Array{String}(size(features_case_2_and_4,1))

  for i = 1:length(labels_case_2_and_4)
    if input_data_case_2_and_4[i,end] == 0
      labels_case_2_and_4[i] = "case 4"
    else
      labels_case_2_and_4[i] = "case 2"
    end
  end

  # train full-tree classifier
  decisiontree2or4 = build_tree(labels_case_2_and_4, features_case_2_and_4)

  decisiontree2or4 = prune_tree(decisiontree2or4, 0.9)

  casemodel = DT(decisiontree1or3, decisiontree2or4)

end

# run ADA

res = @time adagpmcmc(problem_training, problem, gp, casemodel, cov_matrix)

# write results
mcmc_results = Result(res[1].Theta_est, res[1].loglik_est, res[1].accept_vec, res[1].prior_vec)

# write output
Theta = mcmc_results.Theta_est
loglik = mcmc_results.loglik_est

accept_vec = mcmc_results.accept_vec
prior_vec = mcmc_results.prior_vec

loglik_avec_priorvec = zeros(3, length(loglik))
loglik_avec_priorvec[1,:] = loglik
loglik_avec_priorvec[2,:] = accept_vec
loglik_avec_priorvec[3,:] = prior_vec

algorithm_parameters = zeros(10, 2)

algorithm_parameters[1,1] = problem.alg_param.burn_in + 1
algorithm_parameters[2:4,1] = problem.model_param.theta_true
algorithm_parameters[5:7,1] = problem.model_param.theta_0
algorithm_parameters[8:end,:] = problem.prior_dist.prior_parameters

if !accelerated_da
  writetable("Results/Theta_dagpmcmc_lunarc.csv", convert(DataFrame, Theta))
  writetable("Results/loglik_avec_priorvec_dagpmcmc_lunarc.csv", convert(DataFrame, loglik_avec_priorvec))
  writetable("Results/algorithm_parameters_dagpmcmc_lunarc.csv", convert(DataFrame, algorithm_parameters))
else
  writetable("Results/Theta_adagpmcmc_lunarc.csv", convert(DataFrame, Theta))
  writetable("Results/loglik_avec_priorvec_adagpmcmc_lunarc.csv", convert(DataFrame, loglik_avec_priorvec))
  writetable("Results/algorithm_parameters_adagpmcmc_lunarc.csv", convert(DataFrame, algorithm_parameters))
end

ext: jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
using Base, OrthogonalPolynomialsQuasi, ContinuumArrays, QuasiArrays, FillArrays,
        LazyArrays, BandedMatrices, LinearAlgebra, FastTransforms, IntervalSets,
        InfiniteLinearAlgebra, Test
using ForwardDiff, SemiseparableMatrices, SpecialFunctions, LazyBandedMatrices
import ContinuumArrays: BasisLayout, MappedBasisLayout
import OrthogonalPolynomialsQuasi: jacobimatrix, ∞, ChebyshevInterval, LegendreWeight,
            Clenshaw, forwardrecurrence!, singularities
import LazyArrays: ApplyStyle, colsupport, MemoryLayout, arguments
import SemiseparableMatrices: VcatAlmostBandedLayout
import QuasiArrays: MulQuasiMatrix
import Base: OneTo
import InfiniteLinearAlgebra: KronTrav, Block
import FastTransforms: clenshaw!

@testset "singularities" begin
    x = Inclusion(ChebyshevInterval())
    @test singularities(x) == singularities(exp.(x)) == singularities(x.^2) ==
        singularities(x .+ 1) == singularities(1 .+ x) == singularities(x .+ x) ==
        LegendreWeight()
    @test singularities(exp.(x) .* JacobiWeight(0.1,0.2)) ==
        singularities(JacobiWeight(0.1,0.2) .* exp.(x)) ==
        JacobiWeight(0.1,0.2)

    x = Inclusion(0..1)
    @test singularities(x) == singularities(exp.(x)) == singularities(x.^2) ==
        singularities(x .+ 1) == singularities(1 .+ x) == singularities(x .+ x) ==
        legendreweight(0..1)
end

@testset "weight with coeffs" begin
    x = Inclusion(ChebyshevInterval())
    w = π*x.^2 .+ 2
    Pw = LanczosPolynomial(w)
    J = Pw\(x.*Pw)
    x = 1
    n = 10
    x*Pw[x,1] ≈ J[1,1]*Pw[x,1] + J[1,2]*Pw[x,2]
    x*Pw[x,n] ≈ J[n,n-1]*Pw[x,n-1] + J[n,n]*Pw[x,n] + J[n,n+1]*Pw[x,n+1]
end

@testset "Chebyshev weight times polynomial" begin
    x = Inclusion(ChebyshevInterval())
    w = ChebyshevTWeight()
    ϕ = x.^2
    P = LanczosPolynomial(w.*ϕ)
    J = P\(x.*P)
    x = 1
    n = 10
    x*P[x,1] ≈ J[1,1]*P[x,1] + J[1,2]*P[x,2]
    x*P[x,n] ≈ J[n,n-1]*P[x,n-1] + J[n,n]*P[x,n] + J[n,n+1]*P[x,n+1]
end

include("test_chebyshev.jl")
include("test_legendre.jl")
include("test_ultraspherical.jl")
include("test_jacobi.jl")
include("test_fourier.jl")
include("test_odes.jl")
include("test_normalized.jl")
include("test_lanczos.jl")

@testset "Auto-diff" begin
    U = Ultraspherical(1)
    C = Ultraspherical(2)

    f = x -> ChebyshevT{eltype(x)}()[x,5]
    @test ForwardDiff.derivative(f,0.1) ≈ 4*U[0.1,4]
    f = x -> ChebyshevT{eltype(x)}()[x,5][1]
    @test ForwardDiff.gradient(f,[0.1]) ≈ [4*U[0.1,4]]
    @test ForwardDiff.hessian(f,[0.1]) ≈ [8*C[0.1,3]]

    f = x -> ChebyshevT{eltype(x)}()[x,1:5]
    @test ForwardDiff.derivative(f,0.1) ≈ [0;(1:4).*U[0.1,1:4]]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ b2f02100-576f-11eb-3453-ff071cb563df
begin
	using Pkg; Pkg.add(["Images", "ImageIO", "ImageMagick"])
	using Images
end

# ╔═╡ 19b34ee0-576f-11eb-350f-2bbb4050b1db
element = 1 // 1

# ╔═╡ 33005be0-576f-11eb-06d7-0d2df9299f71
fill(element, 3, 4)

# ╔═╡ 41dd7170-576f-11eb-1ae2-cb1e4b6a63bc
typeof(element)

# ╔═╡ 491d1080-576f-11eb-2ce6-fd4f20efecda
keeptrack = (typeof(1), typeof(1.0), typeof("one"), typeof(1 // 1))

# ╔═╡ 5a3b1560-576f-11eb-2584-ad4e0fa643b8
typeof(keeptrack)

# ╔═╡ a84b36e0-576f-11eb-208a-e7b829b847a9
A1 = rand(1:9, 3, 4)

# ╔═╡ 8f3617e0-5771-11eb-23ff-b1f1f5f85951

# ╔═╡ Cell order:
# ╠═b2f02100-576f-11eb-3453-ff071cb563df
# ╠═19b34ee0-576f-11eb-350f-2bbb4050b1db
# ╠═33005be0-576f-11eb-06d7-0d2df9299f71
# ╠═41dd7170-576f-11eb-1ae2-cb1e4b6a63bc
# ╠═491d1080-576f-11eb-2ce6-fd4f20efecda
# ╠═5a3b1560-576f-11eb-2584-ad4e0fa643b8
# ╠═a84b36e0-576f-11eb-208a-e7b829b847a9
# ╠═8f3617e0-5771-11eb-23ff-b1f1f5f85951

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
module TerminalSpinners

export Spinner, autospin, spin!, @spin

ESC = "\e"
CSI = "\e["
DEC_RST  = "l"
DEC_SET  = "h"
DEC_TCEM = "?25"

const TTY_HIDE_CURSOR = CSI * DEC_TCEM * DEC_RST
const TTY_SHOW_CURSOR = CSI * DEC_TCEM * DEC_SET
const TTY_CLEAR_LINE = CSI * "2K"

include("spinners.jl")

mutable struct SpinnerState
    idx::Int
end

struct Spinner{IO_t <: IO}
    frames::Vector{String}
    interval::Float64 # [s]
    stream::IO_t
    msg::Union{String, Nothing}
    channel::Channel{Symbol}
    state::SpinnerState
end

function Spinner(msg::Union{String, Nothing}=nothing; format::Symbol=:classic, output=stderr)
    format = get(FORMATS, format, nothing)
    format == nothing && error("unknown format :$format")
    frames, frequency = format[:frames], format[:frequency]
    return Spinner(frames, 1/frequency, output, msg, Channel{Symbol}(1), SpinnerState(1))
end

const SPINNER_TOKEN = ":spinner"

function spin!(s::Spinner)
    s.state.idx += 1
    s.state.idx > length(s.frames) && (s.state.idx = 1)
    print(s.stream, "\r")
    v = s.frames[s.state.idx]
    if s.msg === nothing
        print(s.stream, v)
    else
        print(s.stream, replace(s.msg, SPINNER_TOKEN => v))
    end
    return
end

function _autospin(s::Spinner)
    while true
        t = @elapsed begin
            if isready(s.channel)
                should_stop = handle_command(s.channel)
                should_stop && break
            end
            spin!(s)
        end
        sleep_time = max(0, s.interval - t)
        sleep(sleep_time)
    end
end

autospin(s::Spinner) = Threads.@spawn _autospin(s)

function handle_command(c::Channel)
    while true
        command = take!(c)
        if command !== :pause && command !== :resume && command !== :stop
            error("unknown command :$command")
        end
        command == :pause && continue
        command == :resume && return false
        command == :stop && return true
    end
end

pause!(s::Spinner)  = put!(s.channel, :pause)
resume!(s::Spinner) = put!(s.channel, :resume)
stop!(s::Spinner)   = put!(s.channel, :stop)

macro spin(spinner, task)
    return quote
        t = autospin($spinner)
        try
            return $task
        finally
            stop!($spinner)
            print($spinner.stream, TTY_SHOW_CURSOR)
            wait(t)
        end
    end
end

end # module


text:
# Tests automatic differentiation on sn.jl:
include("../src/sn_jacobian.jl")

@testset "sn_jacobian" begin

r0 = [0.01,100.0]
nb = 20
#l_max = 20
l_max = 10
n_max = l_max^2+2*l_max

# Now, carry out finite-difference derivative computation:
function sn_jac_num(l_max::Int64,r::T,b::T) where {T <: Real}
  dq = big(1e-24)
  n_max = l_max^2+2*l_max
# Allocate an array for s_n:
  sn_big = zeros(BigFloat,n_max+1)
# Make BigFloat versions of r & b:
  r_big = big(r); b_big = big(b)
# Compute s_n to BigFloat precision:
  s_n!(l_max,r_big,b_big,sn_big)
# Now, compute finite differences:
  sn_jac_big= zeros(BigFloat,n_max+1,2)
  sn_plus = zeros(BigFloat,n_max+1)
  s_n!(l_max,r_big+dq,b_big,sn_plus)
  sn_minus = zeros(BigFloat,n_max+1)
  s_n!(l_max,r_big-dq,b_big,sn_minus)
  sn_jac_big[:,1] = (sn_plus-sn_minus)*.5/dq
  s_n!(l_max,r_big,b_big+dq,sn_plus)
  s_n!(l_max,r_big,b_big-dq,sn_minus)
  sn_jac_big[:,2] = (sn_plus-sn_minus)*.5/dq
return convert(Array{Float64,2},sn_jac_big)
end

#sn_jacobian
#convert(Array{Float64,2},sn_jac_big)-sn_jacobian

using PyPlot
fig,axes = subplots(1,2)
get_cmap("plasma")
epsilon = 1e-12; delta = 1e-3
i=1
for i=1:2
  r=r0[i]
  if r < 1.0
    b = [linspace(1e-8,epsilon,nb); linspace(epsilon,delta,nb); linspace(delta,r-delta,nb);
#     linspace(r-delta,r-epsilon,nb); linspace(r-epsilon,r,nb); linspace(r,r+epsilon,nb); linspace(r+epsilon,r+delta,nb);
     linspace(r-delta,r-epsilon,nb); linspace(r-epsilon,r+epsilon,nb); linspace(r+epsilon,r+delta,nb);
     linspace(r+delta,1-r-delta,nb); linspace(1-r-delta,1-r-epsilon,nb); linspace(1-r-epsilon,1-r+epsilon,nb);
     linspace(1-r+epsilon,1-r+delta,nb); linspace(1-r+delta,1+r-delta,nb); linspace(1+r-delta,1+r-epsilon,nb);linspace(1+r-epsilon,1+r-1e-8,nb)]
  else
    b = [linspace(r-1+1e-8,r-1+epsilon,nb); linspace(r-1+epsilon,r-1+delta,nb); linspace(r-1+delta,r-delta,nb);
#     linspace(r-delta,r-epsilon,nb); linspace(r-epsilon,r,nb); linspace(r,r+epsilon,nb); linspace(r+epsilon,r+delta,nb);
     linspace(r-delta,r-epsilon,nb); linspace(r-epsilon,r+epsilon,nb); linspace(r+epsilon,r+delta,nb);
     linspace(r+delta,r+1-delta,nb); linspace(r+1-delta,r+1-epsilon,nb); linspace(r+1-epsilon,r+1-1e-8,nb)]
  end
  k2 = (1.-(r-b).^2)./(4.*b.*r)
  k2H = ((b+1).^2-r.^2)./(4b)
  igrid=linspace(1,length(b),length(b))-1
  sn_jac_grid = zeros(length(b),n_max+1,2)
  sn_grid = zeros(length(b),n_max+1)
  sn_jac_grid_num = zeros(length(b),n_max+1,2)
  for j=1:length(b)
    println("r: ",r," b: ",b[j])
    sn,sn_jacobian= sn_jac(l_max,r,b[j])
    sn_grid[j,:]=sn
    sn_jac_grid[j,:,:]=sn_jacobian
    sn_jac_grid_num[j,:,:]= sn_jac_num(l_max,r,b[j])
  end
# Now, make plots:
  ax = axes[i]
  m=0;l=0
  for n=0:n_max
    if m==0
#      ax[:plot](abs.(sn_jac_grid[:,n+1,1]-sn_jac_grid_num[:,n+1,1]),color=cmap(l/(l_max+2)),lw=1,label=string("l=",l))
#      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-sn_jac_grid_num[:,n+1,1]),lw=1,label=string("l=",l))
      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1,label=string("l=",l))
    else
#      ax[:plot](abs.(sn_jac_grid[:,n+1,1]-sn_jac_grid_num[:,n+1,1]),color=cmap(l/(l_max+2)),lw=1)
#      ax[:semilogy](abs.(sn_jac_grid[:,n+1,1]-sn_jac_grid_num[:,n+1,1]),lw=1)
      ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,1])-asinh.(sn_jac_grid_num[:,n+1,1])),lw=1)
    end
#    ax[:plot](abs.(sn_jac_grid[:,n+1,2]-sn_jac_grid_num[:,n+1,2]),color=cmap(l/(l_max+2)),lw=1)
#    ax[:semilogy](abs.(sn_jac_grid[:,n+1,2]-sn_jac_grid_num[:,n+1,2]),lw=1)
    ax[:semilogy](abs.(asinh.(sn_jac_grid[:,n+1,2])-asinh.(sn_jac_grid_num[:,n+1,2])),lw=1)
    m +=1
    if m > l
      l += 1
      m = -l
    end
  end
  ax[:legend](loc="upper right",fontsize=6)
  ax[:set_xlabel]("b values")
  ax[:set_ylabel]("Derivative Error")
  ax[:axis]([0,length(b),1e-16,1])
end
#read(STDIN,Char)

# Loop over n and see where the differences between the finite-difference
# and AutoDiff are greater than the derivative value:
l=0; m=0
for n=0:n_max
  diff = sn_jac_grid[:,n+1,:]-sn_jac_grid_num[:,n+1,:]
  for ib=1:length(b)
    @test isapprox(sn_jac_grid[ib,n+1,1],sn_jac_grid_num[ib,n+1,1])
    @test isapprox(sn_jac_grid[ib,n+1,2],sn_jac_grid_num[ib,n+1,2])
  end
  end
  mask = (abs.(diff) .> 1e-3*abs.(sn_jac_grid[:,n+1,:])) .& (abs.(sn_jac_grid[:,n+1,:]) .> 1e-5)
  if sum(mask) > 0 || mod(l-m,4) == 0
    println("n: ",n," max dsn/dr: ",maximum(abs.(diff)))
    println("b: ",b[mask[:,1]],b[mask[:,2]]," k^2_H: ",k2H[mask[:,1]],k2H[mask[:,2]])
# end
#end
# Now, plot derivatives for each value of n:
#for n=0:n_max
    clf()
#    plot(sn_grid[:,n+1])
#    plot(sn_jac_grid[:,n+1,1])
#    plot(igrid[mask[:,1]],sn_jac_grid[mask[:,1],n+1,1],"o")
#    plot(sn_jac_grid[:,n+1,2])
#    plot(igrid[mask[:,2]],sn_jac_grid[mask[:,2],n+1,2],"o")
#    plot(igrid,sn_jac_grid_num[:,n+1,1],linestyle="--",linewidth=2)
#    plot(igrid,sn_jac_grid_num[:,n+1,2],linestyle="--",linewidth=2)
    plot(b,sn_grid[:,n+1])
    plot(b,sn_jac_grid[:,n+1,1])
    plot(b[mask[:,1]],sn_jac_grid[mask[:,1],n+1,1],"o")
    plot(b,sn_jac_grid[:,n+1,2])
    plot(b[mask[:,2]],sn_jac_grid[mask[:,2],n+1,2],"o")
    plot(b,sn_jac_grid_num[:,n+1,1],linestyle="--",linewidth=2)
    plot(b,sn_jac_grid_num[:,n+1,2],linestyle="--",linewidth=2)
    println("n: ",n," l: ",l," m: ",m," mu: ",l-m," nu: ",l+m)
#    read(STDIN,Char)
  end
  m +=1
  if m > l
    l += 1
    m = -l
  end
end

#end

#loglog(abs.(reshape(sn_jac_grid,length(b)*(n_max+1)*2)),abs.(reshape(sn_jac_grid-sn_jac_grid_num,length(b)*(n_max+1)*2)),".")
#loglog(abs.(reshape(sn_jac_grid,length(b)*(n_max+1)*2)),abs.(reshape(sn_jac_grid_num,length(b)*(n_max+1)*2)),".")
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
@testset "Controller: Linearized Approximations" begin
	names = ["particle", "quadruped"]

	ss = [get_simulation("particle", "quadratic_bowl_3D_lc", "quadratic"),
	      get_simulation("quadruped", "flat_2D_lc", "flat")]

	for s in ss
		nz = ContactImplicitMPC.num_var(s.model, s.env)
		nθ = ContactImplicitMPC.num_data(s.model)
		z = rand(nz)
		θ = rand(nθ)
		κ = 1.0e-5

		r0 = rand(nz)
		rz0 = similar(s.rz, Float64)

		lin = ContactImplicitMPC.LinearizedStep(s, z, θ, κ)

		@test norm(lin.z - z) < 1.0e-8
		@test norm(lin.θ - θ) < 1.0e-8
		@test norm(lin.κ - κ) < 1.0e-8

		# Test r!: the linear linearization is exact at the linearization point
		s.res.r!(r0, z, θ, κ)
		@test norm(r0 - lin.r, Inf) < 1.0e-8

		s.res.rz!(rz0, z, θ)
		@test norm(rz0 - lin.rz, Inf) < 1.0e-8

		#TODO: add more tests
	end
end


text:
export CDS_MVp_MT_subfunc
"""
wrapper around Julia threads to compute a multi-threaded matrix vector product in CDS formulated
This is a subfunction of CDS_MVp_MT.jl
"""
function CDS_MVp_MT_subfunc(
        R  ::Array{TF,2},
        x  ::Vector{TF},
        y  ::Vector{TF},
        r0 ::Int,
        c0 ::Int,
        r1 ::Int,
        i  ::Int) where {TF<:Real}

  @Threads.threads for r = r0 : r1
      c = r - r0 + c0 #original
      @inbounds y[r] = y[r] + R[r,i] * x[c]#original
    end
  return y
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using RandomMatrices
using Base.Test

srand(4)

# Excursions in formal power series (fps)
MaxSeriesSize=8
MaxRange = 20
MatrixSize=10
TT=Int64

(n1, n2, n3) = rand(1:MaxSeriesSize, 3)

X = FormalPowerSeries{TT}(rand(1:MaxRange, n1))
Y = FormalPowerSeries{TT}(rand(1:MaxRange, n2))
Z = FormalPowerSeries{TT}(rand(1:MaxRange, n3))

c = rand(1:MatrixSize) #Size of matrix representation to generate

nzeros = rand(1:MaxSeriesSize)
@test X == trim(X)
XX = deepcopy(X)
for i=1:nzeros
    idx = rand(1:MaxRange)
    if !haskey(XX.c, idx)
        XX.c[idx] = convert(TT, 0)
    end
end
@test trim(XX) == X

#Test addition, p.15, (1.3-4)
@test X+X == 2X
@test X+Y == Y+X
@test MatrixForm(X+Y,c) == MatrixForm(X,c)+MatrixForm(Y,c)

#Test subtraction, p.15, (1.3-4)
@test X-X == 0X
@test X-Y == -(Y-X)
@test MatrixForm(X-Y,c) == MatrixForm(X,c)-MatrixForm(Y,c)

#Test multiplication, p.15, (1.3-5)
@test X*Y == Y*X
@test MatrixForm(X*X,c) == MatrixForm(X,c)*MatrixForm(X,c)
@test MatrixForm(X*Y,c) == MatrixForm(X,c)*MatrixForm(Y,c)
@test MatrixForm(Y*X,c) == MatrixForm(Y,c)*MatrixForm(X,c)
@test MatrixForm(Y*Y,c) == MatrixForm(Y,c)*MatrixForm(Y,c)

@test X.*Y == Y.*X

#The reciprocal series has associated matrix that is the matrix inverse
#of the original series
#Force reciprocal to exist
X.c[0] = 1
discrepancy = if Base.VERSION < v"0.5-"
    (norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c), 0:c-1)))
else
    (norm(inv(float(MatrixForm(X,c)))[1, :][:, 1] - tovector(reciprocal(X, c), 0:c-1)))
end
tol = c*√eps()
if discrepancy > tol
    error(string("Error ", discrepancy, " exceeds tolerance ", tol))
end
#@test norm(inv(float(MatrixForm(X,c)))[1, :]'[:, 1] - tovector(reciprocal(X, c),c)) < c*sqrt(eps())

#Test differentiation
XX = derivative(X)
for (k, v) in XX.c
    k==0 && continue
    @test X.c[k+1] == v/(k+1)
end

#Test product rule [H, Sec.1.4, p.19]
@test derivative(X*Y) == derivative(X)*Y + X*derivative(Y)

#Test right distributive law of composition [H, Sec.1.6, p.38]
@test compose(X,Z)*compose(Y,Z) == compose(X*Y, Z)

#Test chain rule [H, Sec.1.6, p.40]
@test derivative(compose(X,Y)) == compose(derivative(X),Y)*derivative(Y)


text:
using Random
using LinearAlgebra
using GenericSVD

using ArbNumerics  # ArbFloat
using Decimals     # Decimal
using DecFP        # Decimal128
using DoubleFloats # Double64
using Quadmath     # Float128
using MultiFloats  # Float64x2

setextrabits(0)
setprecision(BigFloat, 106)
setprecision(ArbFloat, 106)

function test_qr(::Type{T}, n::Int) where {T}
    Random.seed!(0)
    A = T.(rand(Float64, n, n))
    result, elapsed_time, _, _, _ = @timed qr(A)
    num_correct_digits = -log10(Float64(sum(abs.(result.Q * result.R - A))))
    println(rpad(T, 30), " |  qr  | ", round(num_correct_digits, digits=1), " | ", elapsed_time)
end

function test_pinv(::Type{T}, m::Int, n::Int) where {T}
    Random.seed!(0)
    A = T.(rand(Float64, m, n))
    Apinv, elapsed_time, _, _, _ = @timed pinv(A)
    num_correct_digits = -log10(Float64(sum(abs.(A * Apinv * A - A))))
    println(rpad(T, 30), " | pinv | ", round(num_correct_digits, digits=1), " | ", elapsed_time)
end

const N = 40
const M = 25

function main()
    while true
        test_qr(Float64x2, N)
        test_qr(BigFloat, N)
        test_qr(ArbFloat, N)
        test_qr(Dec128, N)
        test_qr(Double64, N)
        test_qr(Float128, N)
        test_pinv(Float64x2, N, M)
        test_pinv(BigFloat, N, M)
        test_pinv(ArbFloat, N, M)
        test_pinv(Double64, N, M)
        test_pinv(Float128, N, M)
    end
end

main()


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin

    using ProgressMeter
    using Base.Threads

    import Chemostat
    import Chemostat.MetNets
    import Chemostat.MetLP
    import Chemostat.MetEP

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
    const H1 = ChR.Human1
    const HG = H1.HumanGEM
    const AG = H1.HartGEMs
end

## ---------------------------------------------------------------------
# DESCRIPTION
# It use the contextualized Hart GEMs and produce a set of other GEMs
# contextualized using Rath data

## ---------------------------------------------------------------------
# base models
let
    modelid = "base"
    for tissue in AG.TISSUES
        (tissue != "GBM") && continue # Test

        AG.check_modelfile(;modelid, tissue) && continue

        raw_model = AG.load_model(;modelid = :raw, tissue, uncompress = true)
        raw_model = MetNets.force_dims(raw_model)

        hart_model = ChR.prepare_metnet(AG, raw_model; inf_medium = true);
        hart_model = MetNets.force_dims(hart_model)

        println()
        @info("Saving", modelid, tissue)
        AG.save_model(hart_model; modelid, tissue)

        println()
        @info("Done!", tissue)

        println()
    end # for tissue
end

## ---------------------------------------------------------------------
# scaled models
let
    modelid = "scaled"
    lift_bound = 1e8
    scale_factor = 50.0
    sglob(AG; lift_bound, scale_factor)

    for tissue in AG.TISSUES
        (tissue != "GBM") && continue # Test

        println("-"^50)

        @info("Scaling",
            modelid, tissue,
            scale_factor,
            lift_bound
        )
        AG.check_modelfile(;modelid, tissue) && continue

        base_model = AG.load_model(;modelid = "base", tissue, uncompress = true)
        @info("Processing", tissue,
            base_model = size(base_model),
            nzrange = MetNets.nzabs_range(base_model.S)
        )
        scaled_model = MetNets.well_scaled_model(base_model, scale_factor; lift_bound)
        scaled_model = MetNets.force_dims(scaled_model)
        @info("Done", tissue,
            scaled_model = size(scaled_model),
            nzrange = MetNets.nzabs_range(scaled_model.S)
        )

        println()
        @info("Comparing with experiments")
        ChR.compare_with_experimets(scaled_model)

        println()
        @info("Saving", tissue)
        AG.save_model(scaled_model; modelid, tissue)
        println()

    end # for tissue
end

## ---------------------------------------------------------------------
# Fva models
let
    for baseid in ["base", "scaled"], tissue in AG.TISSUES
        (tissue != "GBM") && continue # Test

        base_model = AG.load_model(;modelid = baseid, tissue, uncompress = true)
        modelid = "fva_$(baseid)"

        println("-"^50)
        @info("FVA processing", modelid, tissue,
            base_model = size(base_model),
        )
        AG.check_modelfile(;modelid, tissue) && continue

        # This run in parallel
        fva_model = MetLP.fva_preprocess(
            base_model;
            check_obj = HG.HUMAN_BIOMASS_IDER
        ) |> MetNets.force_dims

        println()
        @info("Done!!", tissue, modelid,
            fva_model = size(fva_model)
        )

        println()
        @info("Comparing with experiments")
        ChR.compare_with_experimets(fva_model)

        println()
        @info("Saving", modelid, tissue)
        AG.save_model(fva_model; modelid, tissue)

        println()
        @info("Done!", modelid, tissue)

        println()

    end
end

## ---------------------------------------------------------------------
# stst models
let
    for stst in Rd.STSTS, baseid in ["scaled"], tissue in AG.TISSUES
        (tissue != "GBM") && continue # Test
        (stst != "A") && continue # Test

        modelid = "stst_$(baseid)"
        fva_model = AG.load_model(;modelid = "fva_$(baseid)", tissue, uncompress = true)

        println("-"^50)
        @info("STST processing",
            modelid, tissue, stst,
            fva_model = size(fva_model),
        )
        AG.check_modelfile(;modelid, tissue, stst) && continue

        # stst contextualization
        intake_info = HG.stst_base_intake_info(stst)
        xi = Rd.val(:ξ, stst)
        Chemostat.apply_bound!(fva_model, xi, intake_info;
            emptyfirst = true, ignore_miss = true
        )

        # This run in parallel
        stst_model = MetLP.fva_preprocess(fva_model;
            check_obj = HG.HUMAN_BIOMASS_IDER
        ) |> MetNets.force_dims

        println()
        @info("Done!!", tissue, modelid,
            stst_model = size(stst_model)
        )

        println()
        @info("Comparing with experiments")
        ChR.compare_with_experimets(stst_model)

        println()
        @info("Saving", modelid, tissue)
        AG.save_model(stst_model; modelid, tissue, stst)

        println()
        @info("Done!", modelid, tissue, stst)

        println()
        return

    end
end


text:
degtorad(deg) = deg*π/180

# adapted from Distances.jl
@inline function haversine(λ₁, φ₁, λ₂, φ₂)
    Δλ = λ₂ - λ₁  # longitudes
    Δφ = φ₂ - φ₁  # latitudes

    # haversine formula
    a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2

    # # distance on the sphere
    2 * 6_371_000 * asin( min(√a, one(a)) ) # take care of floating point errors
end

function dist_pairwise(long, lat)
    dim = size(long, 1)
    dist = Array{Float64}(undef, dim, dim)
    @fastmath @inbounds for i ∈ axes(dist, 1), j ∈ 1:i
        dist[i, j] = haversine(long[i], lat[i], long[j], lat[j])
    end
    return Symmetric(dist, :L)
end

function dist_pointwise!(out, long_point, lat_point, long_vec, lat_vec)
    @fastmath @inbounds for i ∈ eachindex(out)
        out[i] = haversine(long_point, lat_point, long_vec[i], lat_vec[i])
    end
    return out
end

function fill_diag!(matrix, value)
    foreach(i -> matrix[i, i] = value, 1:size(matrix, 1))
end


text:
"""
    fastindices(i...)

Wrap ```ComponentArray``` symbolic indices in ```Val```s for type-stable indexing.

# Examples
```julia-repl
julia> using ComponentArrays, BenchmarkTools

julia> ca = ComponentArray(a=1, b=[2, 1, 4], c=(a=2, b=[1, 2]))
ComponentArray{Float64}(a = 1.0, b = [2.0, 1.0, 4.0], c = (a = 2.0, b = [1.0, 2.0]))

julia> ca2 = ca .* ca'
7×7 ComponentArray{Tuple{Axis{(a = 1, b = 2:4, c = (5:7, (a = 1, b = 2:3)))},Axis{(a = 1, b = 2:4, c = (5:7, (a = 1, b =
2:3)))}},Float64,2,Array{Float64,2}}:
 1.0  2.0  1.0   4.0  2.0  1.0  2.0
 2.0  4.0  2.0   8.0  4.0  2.0  4.0
 1.0  2.0  1.0   4.0  2.0  1.0  2.0
 4.0  8.0  4.0  16.0  8.0  4.0  8.0
 2.0  4.0  2.0   8.0  4.0  2.0  4.0
 1.0  2.0  1.0   4.0  2.0  1.0  2.0
 2.0  4.0  2.0   8.0  4.0  2.0  4.0

julia> _a, _b, _c = fastindices(:a, :b, :c)
(Val{:a}(), Val{:b}(), Val{:c}())

julia> @btime \$ca2[:c,:c];
  12.199 μs (2 allocations: 80 bytes)

julia> @btime \$ca2[\$_c, \$_c];
  14.728 ns (2 allocations: 80 bytes)
```
"""
fastindices(i...) = toval.(i)
fastindices(i::Tuple) = toval.(i)

# Make a Val if input isn't already one
toval(x::Val) = x
toval(x) = Val(x)
toval(x::String) = Val(Symbol(x))

# Get value from Val type
getval(::Val{x}) where x = x
getval(::Type{Val{x}}) where x = x

# Split an array up into partitions where the Ns are partition sizes on each dimension
partition(A) = A
partition(a, N...) = partition(a, N)
# Faster method for vectors and matrices
partition(v, N) = (view(v, i:i+N-1) for i in firstindex(v):N:lastindex(v))
function partition(m, N1, N2)
    ax = axes(m)
    firsts = firstindex.(ax)
    lasts = lastindex.(ax)
    return (view(m, i:i+N1-1, j:j+N2-1) for i in firsts[1]:N1:lasts[1], j in firsts[2]:N2:lasts[2])
end
# Slower fallback for higher dimensions
function partition(a::A, N::Tuple) where A<:AbstractArray
    ax = axes(a)
    offs = firstindex.(ax)
    return (view(a, (:).((I.I .- 1) .* N .+ offs, ((I.I .- 1) .* N .+ N .- 1 .+ offs))...) for I in CartesianIndices(div.(size(a), N)))
end
# partition(a::A, N::Tuple) where A<:AbstractVector = reshape(view(a, :), N)

# Faster filtering of tuples by type
filter_by_type(::Type{T}, args...) where T = filter_by_type(T, (), args...)
filter_by_type(::Type{T}, part::Tuple) where T = part
filter_by_type(::Type{T}, part::Tuple, ax, args...) where T = filter_by_type(T, part, args...)
filter_by_type(::Type{T}, part::Tuple, ax::T, args...) where T = filter_by_type(T, (part..., ax), args...)

# Flat length of an arbitrarily nested named tuple
recursive_length(x) = length(x)
recursive_length(a::AbstractArray{T,N}) where {T<:Number,N} = length(a)
recursive_length(a::AbstractArray) = recursive_length.(a) |> sum
recursive_length(nt::NamedTuple) = values(nt) .|> recursive_length |> sum
recursive_length(::Union{Nothing, Missing}) = 1


text:
@testset "Typical Section" begin

    submodel = Submodel(Section())

    # test provided jacobians
    submodel_jacobian_tests(submodel)

    # test consistency of input/output functions
    submodel_io_tests(submodel)
end


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

import Test
Test.@testset "Package H3DDetectorSystems" begin

# test code goes here

end # testset


text:
## Array.jl ##

for (T, suff) in ((Float64, ""), (Float32, "f"))

    # 1 arg functions
    for f in (:ceil,:floor,:sqrt,:rsqrt,:rec,
              :exp,:exp2,:expm1,:log,:log1p,:log2,:log10,
              :sin,:sinpi,:cos,:cospi,:tan,:tanpi,:asin,:acos,:atan,
              :sinh,:cosh,:tanh,:asinh,:acosh,:atanh)
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T})
                out = Array($T,size(X))
                ($f!)(out, X)
            end
            function ($f!)(out::Array{$T}, X::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{Cint}),out,X,&length(X))
                out
            end
        end
    end

    # renamed 1 arg functions
    for (f,fa) in ((:trunc,:int),(:round,:nint),(:exponent,:logb),
                   (:abs,:fabs))
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T})
                out = Array($T,size(X))
                ($f!)(out, X)
            end
            function ($f!)(out::Array{$T}, X::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{Cint}),out,X,&length(X))
                out
            end
        end
    end

    # 2 arg functions
    for f in (:copysign,:atan2)
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = Array($T,size(X))
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out,X,Y,&length(X))
                out
            end
        end
    end

    # for some bizarre reason, vvpow/vvpowf reverse the order of arguments.
    for f in (:pow,)
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = Array($T,size(X))
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out,Y,X,&length(X))
                out
            end
        end
    end

    # renamed 2 arg functions
    for (f,fa) in ((:rem,:fmod),(:fdiv,:div))
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T}, Y::Array{$T})
                size(X) == size(Y) || throw(DimensionMismatch("Arguments must have same shape"))
                out = Array($T,size(X))
                ($f!)(out, X, Y)
            end
            function ($f!)(out::Array{$T}, X::Array{$T}, Y::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out,X,Y,&length(X))
                out
            end
        end
    end

    # two-arg return
    for f in (:sincos,)
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T})
                out1 = Array($T,size(X))
                out2 = Array($T,size(X))
                ($f!)(out1, out2, X)
            end
            function ($f!)(out1::Array{$T}, out2::Array{$T}, X::Array{$T})
                ccall(($(string("vv",f,suff)),libacc),Void,
                      (Ptr{$T},Ptr{$T},Ptr{$T},Ptr{Cint}),out1,out2,X,&length(X))
                out1, out2
            end
        end
    end

    # complex return
    for (f,fa) in ((:cis,:cosisin),)
        f! = symbol("$(f)!")
        @eval begin
            function ($f)(X::Array{$T})
                out = Array(Complex{$T},size(X))
                ($f!)(out, X)
            end
            function ($f!)(out::Array{Complex{$T}}, X::Array{$T})
                ccall(($(string("vv",fa,suff)),libacc),Void,
                      (Ptr{Complex{$T}},Ptr{$T},Ptr{Cint}),out,X,&length(X))
                out
            end
        end
    end
end

# Functions over single vectors that return scalars/tuples
for (T, suff) in ((Float32, ""), (Float64, "D"))

    for (f, fa) in ((:maximum, :maxv), (:minimum, :minv), (:mean, :meanv),
                    (:meansqr, :measqv), (:meanmag,  :meamgv), (:sum, :sve))
        @eval begin
            function ($f)(X::Vector{$T})
                val = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa, suff), libacc)),  Void,
                      (Ptr{$T}, Int64,  Ref{$T}, UInt64),
                      X, 1, val, length(X))
                return val[]
            end
        end
    end

    for (f, fa) in ((:findmax, :maxvi), (:findmin, :minvi))
        @eval begin
            function ($f)(X::Vector{$T})
                index = Ref{Int}(0)
                val = Ref{$T}(0.0)
                ccall(($(string("vDSP_", fa, suff), libacc)),  Void,
                      (Ptr{$T}, Int64,  Ref{$T}, Ref{Int}, UInt64),
                      X, 1, val, index, length(X))
                return (val[], index[]+1)
            end
        end
    end
end

# Element-wise operations over two vectors
for (T, suff) in ((Float32, ""), (Float64, "D"))

    for (f, name) in ((:vadd, "addition"),  (:vsub, "subtraction"),
                      (:vdiv, "division"), (:vmul, "multiplication"))
        f! = symbol("$(f)!")

        @eval begin
            @doc """
            `$($f!)(result::Vector{$($T)}, X::Vector{$($T)}, Y::Vector{$($T)})`

            Implements element-wise **$($name)** over two **Vector{$($T)}** and overwrites
            the result vector with computed value. *Returns:* **Vector{$($T)}** `result`
            """ ->
            function ($f!)(result::Vector{$T}, X::Vector{$T}, Y::Vector{$T})
                ccall(($(string("vDSP_", f, suff), libacc)),  Void,
                      (Ptr{$T}, Int64, Ptr{$T},  Int64, Ptr{$T}, Int64,  UInt64),
                      Y, 1, X, 1, result, 1, length(result))
                return result
            end
        end

        @eval begin
            @doc """
            `$($f)(X::Vector{$($T)}, Y::Vector{$($T)})`

            Implements element-wise **$($name)** over two **Vector{$($T)}**. Allocates
            memory to store result. *Returns:* **Vector{$($T)}**
            """ ->
            function ($f)(X::Vector{$T}, Y::Vector{$T})
                result = similar(X)
                ($f!)(result, X, Y)
                return result
            end
        end
    end
end


text:
using NeuralQuantum, Test
using NeuralQuantum: set_index!
num_types   = [Float32, Float64]
atol_types  = [1e-5, 1e-8]

re_machines = Dict()
im_machines = Dict()

ma = (T, hi) -> RBMSplit(T, hi, 2)
im_machines["RBMSplit"] = ma

ma = (T, hi) -> RBM(T, hi, 2, NeuralQuantum.logℒ)
im_machines["RBM_softplus"] = ma

ma = (T, hi) -> RBM(T, N, 2, NeuralQuantum.logℒ2)
im_machines["RBM_cosh"] = ma

ma = (T, hi) -> NDM(T, hi, 1, 2, NeuralQuantum.logℒ)
re_machines["NDM_softplus"] = ma

ma = (T, hi) -> NDM(T, hi, 1, 2, NeuralQuantum.logℒ2)
re_machines["NDM_cosh"] = ma

ma = (T, hi) -> NDMSymm(T, hi, 1, 2, translational_symm_table(HyperCube([N], periodic=true)), NeuralQuantum.logℒ2)
re_machines["NDMSymm_cosh"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    PureStateAnsatz(Chain(Dense(T, N, N*2), Dense(T, N*2, N*3), WSum(T, N*3)), hi)
end
re_machines["ChainKet"] = ma

ma = (T, hi) -> begin
    N =nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_softplus), sum_autobatch)
    return PureStateAnsatz(ch, hi)
end
im_machines["chain_pure_softplus"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_logcosh), Dense(T, 2*N, N, af_logcosh), sum_autobatch)
    return PureStateAnsatz(ch, hi)
end
im_machines["chain_pure_cosh_2"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(DenseSplit(T, N, 2*N, af_softplus), sum_autobatch)
    return MixedStateAnsatz(ch, hi)
end
im_machines["chain_mixed_softplus"] = ma

all_machines = merge(re_machines, im_machines)

N = 4

@testset "Test Properties $name" for name=keys(all_machines)
    for T=num_types
        if name ∈ keys(im_machines)
            T = Complex{T}
        end
        name == "ChainKet" && T != Float32 && continue
        hilb = HomogeneousFock(N, 2)

        net = all_machines[name](T,hilb)

        @test NeuralQuantum.out_type(net) == Complex{real(T)}
        @test NeuralQuantum.is_analytic(net)
    end
end

@testset "Test cached evaluation $name" for name=keys(all_machines)
    for T=num_types
        name == "ChainKet" && T != Float32 && continue
        hilb = HomogeneousFock(N, 2)

        net = all_machines[name](T, hilb)
        cnet = cached(net)

        if net isa NeuralQuantum.MatrixNet
            hilb = SuperOpSpace(hilb)
        end

        v = state(T, hilb, net)
        # compute exact
        vals = []; cvals = [];
        for i=1:spacedimension(hilb)
            set!(v, hilb, i)
            push!(vals, net(v))
            push!(cvals, cnet(v))
        end
        # Test evaluation
        @test vals ≈ cvals
        # Test cached vals have the same type
        all(typeof.(cvals) .== typeof(first(cvals)))
        # and same type as exact computation
        if all(imag(vals) .== 0)
            @test real(typeof(first(vals))) == real(typeof(first(cvals)))
        else
            @test typeof(first(vals)) == typeof(first(cvals))
        end
    end
end

@testset "Test cached gradient $name" for name=keys(all_machines)
    for (T,  atol)=zip(num_types, atol_types)
        name == "ChainKet" && T != Float32 && continue
        name == "NDMSymm_cosh" && continue
        hilb = HomogeneousFock(N, 2)

        net = all_machines[name](T, hilb)
        cnet = cached(net)

        if net isa NeuralQuantum.MatrixNet
            hilb = SuperOpSpace(hilb)
        end

        cder = grad_cache(net)

        # compute exact
        grads = [];
        for v=states(T,hilb)
            der_ad = ∇logψ(net,  v)
            ∇logψ!(cder, cnet, v)
            #=for f=propertynames(der_ad)
                ∇  = getproperty(der_ad,  f)
                c∇ = getproperty(cder, f)
                push!(grads, ∇ ≈ c∇)
            end=#
            isgood = isapprox(der_ad, cder, atol=atol)
            push!(grads, isgood)
        end
        @test all(grads)
    end
end


text:
abstract type AbstractNote{X,I} end

struct Noted{X,I}
    value :: X
    info :: I
end

function Base.show(io::IO, n::Noted)
    print(io, "Noted(")
    show(io, n.value)
    print(io, ", ")
    show(io, n.info)
    print(io, ")")
end

export value
value(n::Noted) = n.value
value(nt::NamedTuple) = nt.value
value(other) = other

export info
info(n::Noted) = n.info
info(other) = other

# s = StructArray{Noted}((randn(10),Fill(5,(10,))))

# r = StructArray{NamedTuple{(:a, :b, :c),Tuple{Float64,Float64,Float64}}}((randn(10),randn(10),randn(10)))

# StructArray{NamedTuple{(:a, :b),Tuple{Int64,NamedTuple{(:b1, :b2),Tuple{Int64,Int64}}}}}((rand(1:10,4),rand(1:10,4),rand(1:10,4)))

# StructArray{NamedTuple{(:a, :b),Tuple{Int64,NamedTuple{(:b1, :b2),Tuple{Int64,Int64}}}}}((rand(

# StructArray{NamedTuple{(:b1, :b2),Tuple{Int64,Int64}}}((rand(1:10,4),rand(1:10,4)))

# StructArray((1,2,3); names=(:a,:B,:c))