text:
using NBodyIPs, JuLIP, Test, IPFitting, DataFrames
using JuLIP.Potentials: evaluate_d
using IPFitting: Dat, LsqDB
using NBodyIPs: BondLengthDesc, BondAngleDesc
Lsq = IPFitting.Lsq
Err = IPFitting.Errors
using LinearAlgebra: norm

# generate random data
function generate_data(species, L, rmax, N, calc)
   data = Dat[]
   for n = 1:N
      at = bulk(species; cubic=true, pbc=true) * L
      rattle!(at, rand() * rmax)
      E = energy(calc, at)
      F = forces(calc, at)
      push!(data, Dat(at, "rand"; E = E, F = F))
   end
   return data
end

r0 = rnn(:Si)
calc = StillingerWeber()
data = generate_data(:Si, 2, 0.33*r0, 100, calc)

TRANSFORM = "exp( - 2 * (r/$r0 - 1.5) )"

rcut2 = cutoff(calc)*2.0
CUTOFF2 = PolyCutSym(2, rcut2)
D2 = BondLengthDesc(TRANSFORM, CUTOFF2)

rcut3 = cutoff(calc)*2.0
# CUTOFF3 = CosCut(rcut3-1, rcut3)
CUTOFF3 = PolyCut(2, rcut3)
# CUTOFF3 = PolyCut2s(2, 0.0, rnn(:Si), rcut3)

rcut3 = cutoff(calc)*2.0
D3 = BondLengthDesc(TRANSFORM, CUTOFF3)

##
err_erms = Float64[]
err_frms = Float64[]
degrees = [4, 6, 8, 10] # [4, 6, 8, 10, 12]
rr = range(0.9*r0, stop=cutoff(calc), length=200)
for deg3 in degrees
   # B = [B1; gen_basis(2, D2, deg2); gen_basis(3, D3, deg3)]
   B = [nbpolys(2, D2, 10); nbpolys(3, D3, deg3)]
   # B = nbpolys(3, D3, deg3)
   @show length(B)
   db = LsqDB("", B, data)
   IP, fitinfo = Lsq.lsqfit( db,
                         E0 = 0.0,
                         configweights = Dict("rand" => 1.0),
                         obsweights   = Dict("E" => 100.0, "F" => 1.0),
                         combineIP = NBodyIP )
   push!(err_erms, fitinfo["errors"]["relrmse"]["set"]["E"])
   push!(err_frms, fitinfo["errors"]["relrmse"]["set"]["F"])
end

##
df = DataFrame( :degrees => degrees,
                :relrms_E => err_erms,
                :relrms_F => err_frms )
display(df)

(@test minimum(err_erms) < 0.001) |> println
(@test minimum(err_frms) < 0.05) |> println

# BOND LENGTH
# │ Row │ degrees │ rms_E       │ rms_F      │
# ├─────┼─────────┼─────────────┼────────────┤
# │ 1   │ 4       │ 0.00134475  │ 0.118625   │
# │ 2   │ 6       │ 0.000244722 │ 0.0362398  │
# │ 3   │ 8       │ 0.000115287 │ 0.0179257  │
# │ 4   │ 10      │ 5.28124e-5  │ 0.00888845 │

# BOND ANGLE
# │ Row │ degrees │ relrms_E    │ relrms_F   │
# ├─────┼─────────┼─────────────┼────────────┤
# │ 1   │ 4       │ 0.000523011 │ 0.070237   │
# │ 2   │ 6       │ 0.000445928 │ 0.0477666  │
# │ 3   │ 8       │ 0.000100544 │ 0.016773   │
# │ 4   │ 10      │ 3.7269e-5   │ 0.00568809 │
# │ 5   │ 12      │ 2.00043e-5  │ 0.00419453 │


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Profiler driver script.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-11

include("hdr.jl")

require("conllx.jl")
require("dep.jl")
require("hybrid.jl")
require("parse.jl")

using CoNLLX
using DepGraph
using Hybrid
using Parse

reportpath = "$repodir/wrk/profile.txt"

open(talpath, "r") do data_f
    sents = collect(Sentence, conllxparse(data_f, useproj=true))
    Profile.init(2^24, 0.001)

    # "Warm up" the code by performing an initial pass.
    next(train(sents, epochs=1))

    @profile begin
        for _ in train(sents, epochs=2)
            nothing
        end
    end

    open(reportpath, "w") do report_f
        Profile.print(report_f)
    end
end

ext: jl

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
using Statistics
using Plots
pyplot()

using SIMD
using BenchmarkTools
using ChipSort
using Unrolled
using JLD

pc(n) = (1:length(n))./length(n)

data = load("chip-medium-bench.jld")
bmk = data["bmk"]

T=Dict("Int32"=>UInt32, "Int64"=>UInt64, "UInt32"=>UInt32, "UInt64"=>UInt64)[data["T"]]
V=data["V"]
J=data["J"]

rj = bmk[2^13][:JuliaStd]
ra = bmk[2^13][:CombSortMix]
rm = bmk[2^13][:ChipSortM]
rx = bmk[2^13][:ChipSortL]

@show ra
@show rx
@show minimum(ra.times)/minimum(rm.times)
@show median(ra.times)/median(rm.times)

plot(size=(600,400))
plot!(ra.times*1e-3, pc(ra.times), l=2, label="ChipSort")
plot!(rj.times*1e-3, pc(rj.times), l=2, label="Julia std")
plot!(rm.times*1e-3, pc(rm.times), l=2, label="ChipSort-M")
plot!(rx.times*1e-3, pc(rx.times), l=2, label="ChipSort-L")
plot!(xlim=(000e0,400e0), ylabel="Fraction of trials", xlabel="Time [μs]", title="Running time eCDF sorting 8k UInt32 values")
rj = bmk[2^18][:JuliaStd]
ra = bmk[2^18][:CombSortMix]
rm = bmk[2^18][:ChipSortM]
rx = bmk[2^18][:ChipSortL]

savefig("chipsort-bench-8k.pdf")

@show ra
@show rx
@show minimum(ra.times)/minimum(rx.times)
@show median(ra.times)/median(rx.times)
plot(size=(600,400), reuse=false)
plot!(ra.times*1e-6, pc(ra.times), l=2, label="ChipSort")
plot!(rj.times*1e-6, pc(rj.times), l=2, label="Julia std")
plot!(rm.times*1e-6, pc(rm.times), l=2, label="ChipSort-M")
plot!(rx.times*1e-6, pc(rx.times), l=2, label="ChipSort-L")
plot!(xlim=(3,17), ylabel="Fraction of trials", xlabel="Time [ms]", title="Running time eCDF sorting 256k UInt32 values")
savefig("chipsort-bench-256k.pdf")


text:
using GLMakie, StatsBase

function plotforest(forest, nx, ny)
    my_colors = ["#ffffff", "#563524", "#e12719", "#30a144"]

    # Create the graph
    heatmap(1:nx, 1:ny, forest, colormap = my_colors, axis = (aspect = 1, xlabel = "x axis", ylabel = "y axis"))
end

max_stars_repo_name: TheSquareRoot/ForestFire

max_issues_repo_name: TheSquareRoot/ForestFire

max_forks_repo_name: TheSquareRoot/ForestFire


text:
module Algames

greet() = print("Hello World!")

using Altro
using BenchmarkTools
using ForwardDiff
using LinearAlgebra
using Parameters
using Printf
using Random
using RecipesBase
using RobotDynamics
using Rotations
using SparseArrays
using StaticArrays
using TrajectoryOptimization

# Utils
export
    add2sub,
    addI2sub,
    sparse_zero!,
    scn

# Dynamics
export
    AbstractGameModel,
    DoubleIntegratorGame,
    UnicycleGame,
    BicycleGame,
    QuadrotorGame,
    dynamics,
    dim

# Struct
export
    Options,
    IBROptions,
    ProblemSize,
    PrimalDualTraj,
    init_traj!,
    update_traj!,
    set_traj!,
    get_traj!,
    Δ_step,
    Statistics,
    record!,
    reset!,
    DynamicsViolation,
    dynamics_violation,
    ControlViolation,
    control_violation,
    StateViolation,
    state_violation,
    OptimalityViolation,
    optimality_violation,
    Regularizer,
    set!,
    mult!

# Constraints
export
    GameConstraintValues,
    ControlBoundConstraint,
    StateBoundConstraint,
    WallConstraint,
    Wall3DConstraint,
    CylinderConstraint,
    add_collision_avoidance!,
    add_spherical_collision_avoidance!,
    add_state_bound!,
    add_velocity_bound!,
    velocity_index,
    add_control_bound!,
    add_circle_constraint!,
    AbstractWall,
    Wall,
    Wall3D,
    CylinderWall,
    add_wall_constraint!,
    set_constraint_params!,
    reset!,
    reset_duals!,
    reset_penalties!,
    penalty_update!,
    dual_update!,
    dual_update,
    evaluate!,
    jacobian!,
    update_active_set!,
    constraint_jacobian_residual!,
    constraint_residual!

# Core
export
    NewtonCore,
    vertical_indices,
    horizontal_indices,
    idx,
    horizontal_idx,
    vertical_idx,
    residual_views,
    jacobian_views,
    dynamics_indices,
    vertical_mask,
    horizontal_mask,
    Stamp,
    VStamp,
    HStamp,
    stampify,
    stampify!,
    valid

# Objective
export
    GameObjective,
    cost_gradient!,
    cost_hessian!,
    add_collision_cost!,
    CollisionCost

# Problem
export
    Penalty,
    GameProblem,
    residual!,
    residual_jacobian!,
    regularize_residual!, # need test
    regularize_residual_jacobian!, # need test
    residual_norm,
    ibr_residual!,
    ibr_residual_jacobian!,
    regularize_ibr_residual!,
    regularize_ibr_residual_jacobian!,
    dynamics_residual,
    ∇dynamics!,
    newton_solve!,
    inner_iteration,
    line_search,
    ibr_newton_solve!,
    ibr_inner_iteration!,
    ibr_line_search

# Active Set
export
    NullSpace,
    add_matrix!,
    ActiveSetCore,
    complete_vertical_indices,
    complete_horizontal_indices,
    complete_residual_views,
    complete_jacobian_views,
    CStamp,
    active,
    active_vertical_mask!,
    active_horizontal_mask!,
    get_collision_conval,
    update_nullspace!

# Plots
export
    recipe_traj,
    recipe_violation

# Utils
include("utils.jl")

# Dynamics
include("dynamics/game_model.jl")
include("dynamics/double_integrator.jl")
include("dynamics/unicycle.jl")
include("dynamics/bicycle.jl")
include("dynamics/quadrotor.jl")

# Struct
include("struct/problem_size.jl")

# Core
include("core/stamp.jl")
include("core/newton_core.jl")

# Struct
include("struct/primal_dual_traj.jl")
include("struct/regularizer.jl")
include("struct/options.jl")

# Constraints
include("constraints/state_bound_constraint.jl")
include("constraints/control_bound_constraint.jl")
include("constraints/wall_constraint.jl")
include("constraints/cylinder_constraint.jl")
include("constraints/game_constraints.jl")
include("constraints/constraints_methods.jl")
include("constraints/velocity_constraint.jl")

# Struct
include("struct/violations.jl")
include("struct/statistics.jl")

# Objective
include("objective/objective.jl")

# Problem
include("problem/problem.jl")
include("problem/local_quantities.jl")
include("problem/global_quantities.jl")
include("problem/solver_methods.jl")

# Constraints
include("constraints/constraint_derivatives.jl")

# Active Set
include("active_set/active_set_stamp.jl")
include("active_set/active_set_core.jl")
include("active_set/active_set_methods.jl")

# Plots
include("plots/solver_plots.jl")

end # module


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

# thermodynamics

@doc """
$(convertext(:temperature,"mass(U,S)*speed(U,S)^2/boltzmann(U,S)"))

Measurement scale for thermodynamic energy or `temperature` (K), unit conversion factor.

```Julia
julia> temperature(Metric,SI2019) # K⋅K⁻¹
$(temperature(Metric,SI2019))

julia> temperature(English,SI2019) # K⋅°R⁻¹
$(temperature(English,SI2019))

julia> temperature(English,Metric) # K⋅°R⁻¹
$(temperature(English,Metric))

julia> temperature(PlanckGauss,Metric) # K⋅TP⁻¹
$(temperature(PlanckGauss,Metric))
```
""" temperature, Θ

@doc """
$(convertext(:entropy,"energy(U,S)/temperature(U,S)"))

Heat capacity or `energy` per `temperature` or `entropy` (J⋅K⁻¹), unit conversion factor.

```Julia
julia> entropy(Metric,SI2019) # K⋅K⁻¹
$(entropy(Metric,SI2019))

julia> entropy(CGS,Metric) # J⋅erg⁻¹
$(entropy(CGS,Metric))

julia> entropy(English,SI2019) # J⋅°R⋅K⁻¹⋅ft⁻¹⋅lb⁻¹
$(entropy(English,SI2019))

julia> entropy(Survey,English) # ftUS²⋅°R⋅°ft⁻²⋅°R⁻¹
$(entropy(Survey,English))
```
""" entropy

@doc """
$(convertext(:specificentropy,"specificenergy(U,S)/temperature(U,S)"))

Specific heat capacity or `specificentropy` (J⋅K⁻¹⋅kg⁻¹), unit conversion factor.

```Julia
julia> specificentropy(Metric,SI2019) # m²⋅K⋅K⁻¹⋅cm⁻²
$(specificentropy(Metric,SI2019))

julia> specificentropy(CGS,Metric) # m²⋅cm⁻²
$(specificentropy(CGS,Metric))

julia> specificentropy(English,Metric) # m²⋅°R⋅K⁻¹⋅ft⁻²
$(specificentropy(English,Metric))

julia> specificentropy(Survey,English) # ft²⋅°R⋅ftUS⁻²⋅°R⁻¹
$(specificentropy(Survey,English))
```
""" specificentropy

@doc """
$(convertext(:volumeheatcapacity,"entropy(U,S)/volume(U,S)"))

The `entropy` per `volume` or `volumeheatcapacity` (J⋅K⁻¹⋅m⁻³), unit conversion factor.

```Julia
julia> volumeheatcapacity(Metric,SI2019) # K⋅K⁻¹
$(volumeheatcapacity(Metric,SI2019))

julia> volumeheatcapacity(CGS,Metric) # J⋅cm³⋅erg⁻¹⋅m⁻³
$(volumeheatcapacity(CGS,Metric))

julia> volumeheatcapacity(English,SI2019) # J⋅ft²⋅°R⋅K⁻¹⋅lb⁻¹⋅m⁻³
$(volumeheatcapacity(English,SI2019))

julia> volumeheatcapacity(Survey,English) # ftUS⁵°R⋅°ft⁻⁵⋅°R⁻¹
$(volumeheatcapacity(Survey,English))
```
""" volumeheatcapacity

@doc """
$(convertext(:thermalconductivity,"force(U,S)/time(U,S)/temperature(U,S)"))

Heat conductivity or `thermalconductivity` (W⋅m⁻¹⋅K⁻¹), unit conversion factor.

```Julia
julia> thermalconductivity(Metric,SI2019) # K⋅K⁻¹
$(thermalconductivity(Metric,SI2019))

julia> thermalconductivity(CGS,Metric) # N⋅dyn⁻¹
$(thermalconductivity(CGS,Metric))

julia> thermalconductivity(English,Metric) # N⋅°R⋅K⁻¹⋅ft⁻¹⋅lb⁻¹
$(thermalconductivity(English,Metric))
```
""" thermalconductivity

@doc """
$(convertext(:thermalconductance,"thermalconductivity(U,S)*length(U,S)"))

Reciprocal of `thermalresistance` (W⋅K⁻¹), unit conversion factor.

```Julia
julia> thermalconductance(Metric,SI2019) # K⋅K⁻¹
$(thermalconductance(Metric,SI2019))

julia> thermalconductance(CGS,Metric) # W⋅s⋅erg⁻¹
$(thermalconductance(CGS,Metric))

julia> thermalconductance(English,Metric) # J⋅°R⋅K⁻¹⋅ft⁻¹⋅lb⁻¹
$(thermalconductance(English,Metric))
```
""" thermalconductance

@doc """
$(convertext(:thermalresistance,"1/thermalconductivity(U,S)/length(U,S)"))

Resistance to heat flow or `thermalresistance` (K⋅W⁻¹), unit conversion factor.

```Julia
julia> thermalresistance(Metric,SI2019) # K⋅K⁻¹
$(thermalresistance(Metric,SI2019))

julia> thermalresistance(CGS,Metric) # erg⋅s⁻¹⋅W⁻¹
$(thermalresistance(CGS,Metric))

julia> thermalresistance(English,Metric) # ft⋅lb⋅K⋅°R⁻¹⋅J⁻¹
$(thermalresistance(English,Metric))
```
""" thermalresistance

@doc """
$(convertext(:thermalresistivity,"1/thermalconductivity(U,S)"))

Resistance to heat flow or `thermalresistance` (K⋅W⁻¹), unit conversion factor.

```Julia
julia> thermalresistance(Metric,SI2019) # K⋅K⁻¹
$(thermalresistance(Metric,SI2019))

julia> thermalresistance(CGS,Metric) # erg⋅s⁻¹⋅W⁻¹
$(thermalresistance(CGS,Metric))

julia> thermalresistance(English,Metric) # ft⋅lb⋅K⋅°R⁻¹⋅J⁻¹
$(thermalresistance(English,Metric))
```
""" thermalresistivity

@doc """
$(convertext(:thermalexpansion,"1/temperature(U,S)"))

Measurement scale for coefficient of `thermalexpansion` (K⁻¹), unit conversion factor.

```Julia
julia> thermalexpansion(Metric,SI2019) # K⋅K⁻¹
$(thermalexpansion(Metric,SI2019))

julia> thermalexpansion(English,SI2019) # °R⋅K⁻¹
$(thermalexpansion(English,SI2019))

julia> thermalexpansion(English,Metric) # °R⋅K⁻¹
$(thermalexpansion(English,Metric))
```
""" thermalexpansion

@doc """
$(convertext(:lapserate,"temperature(U,S)/length(U,S)"))

Temperature gradient over `length` or `lapserate` (K⋅m⁻¹), unit conversion factor.

```Julia
julia> lapserate(Metric,SI2019) # K⋅K⁻¹
$(lapserate(Metric,SI2019))

julia> lapserate(English,SI2019) # K⋅ft⋅°R⁻¹⋅m⁻¹
$(lapserate(English,SI2019))

julia> lapserate(English,Metric) # K⋅ft⋅°R⁻¹⋅m⁻¹
$(lapserate(English,Metric))

julia> lapserate(EnglishUS,English) # °R⋅ftUS⋅°R⁻¹⋅ft⁻¹
$(lapserate(EnglishUS,English))
```
""" lapserate

# molar

@doc """
$(convertext(:molarmass,"molarmass(S)/molarmass(U)"))

Molar mass or `mass` per `mole` (kg⋅mol⁻¹), unit conversion factor.

```Julia
julia> molarmass(CGS,Metric) # kg⋅g⁻¹
$(molarmass(CGS,Metric))

julia> molarmass(Metric,SI2019) # mol⋅mol⁻¹
$(molarmass(Metric,SI2019))
```
""" molarmass(::UnitSystem,::UnitSystem)

@doc """
$(convertext(:molality,"molarmass(U)/molarmass(S)"))

Molality or `mole` per `mass` (mol⋅kg⁻¹), unit conversion factor.

```Julia
julia> molality(CGS,Metric) # kg⋅g⁻¹
$(molality(CGS,Metric))

julia> molality(Metric,SI2019) # mol⋅mol⁻¹
$(molality(Metric,SI2019))
```
""" molality

@doc """
$(convertext(:molaramount,"mass(U,S)*molality(U,S)"))

Amount of molecular substance or `molaramount` (mol), unit conversion factor.

```Julia
julia> molaramount(SI2019,Metric) # mol⋅mol⁻¹
$(molaramount(SI2019,Metric))

julia> molaramount(British,SI2019) # mol⋅slug-mol⁻¹
$(molaramount(English,SI2019))

julia> molaramount(English,SI2019) # mol⋅lb-mol⁻¹
$(molaramount(English,SI2019))
```
""" molaramount, N
#julia> molaramount(English,SI2019)*ft/g₀ # mol⋅lb-mol⁻¹
#$(molaramount(English,SI2019)*ft/g₀)

@doc """
$(convertext(:molarity,"molaramount(U,S)/volume(U,S)"))

Molar concentration or `molaramount` per `volume` (mol⋅m⁻³), unit conversion factor.

```Julia
julia> molarity(CGS,Metric) # cm³⋅m⁻³
$(molarity(CGS,Metric))

julia> molarity(English,SI2019) # ft³⋅m⁻³
$(molarity(English,SI2019))
```
""" molarity

@doc """
$(convertext(:molarvolume,"volume(U,S)/molaramount(U,S)"))

Occupied `volume` per `molaramount` or `molarvolume` (m³⋅mol⁻¹), unit conversion factor.

```Julia
julia> molarvolume(CGS,Metric) # m³⋅cm⁻³
$(molarvolume(CGS,Metric))

julia> molarvolume(English,SI2019) # m³⋅ft⁻³
$(molarvolume(English,SI2019))
```
""" molarvolume

@doc """
$(convertext(:molarentropy,"entropy(U,S)/molaramount(U,S)"))

Molar heat capacity or `entropy` per `molaramount` (J⋅K⁻¹⋅mol⁻¹), unit conversion factor.

```Julia
julia> molarentropy(CGS,Metric) # J⋅erg⁻¹
$(molarentropy(CGS,Metric))

julia> molarentropy(English,SI2019) # J⋅°R⋅lb-mol⋅ft⁻¹⋅lb⁻¹⋅K⁻¹⋅mol⁻¹
$(molarentropy(English,SI2019))
```
""" molarentropy

@doc """
$(convertext(:molarenergy,"energy(U,S)/molaramount(U,S)"))

Gibbs free `energy` per `mole` or `molarenergy` (J⋅mol⁻¹), unit conversion factor.

```Julia
julia> molarenergy(CGS,Metric) # J⋅erg⁻¹
$(molarenergy(CGS,Metric))

julia> molarenergy(English,SI2019) # J⋅slug-mol⋅ft⁻¹⋅lb⁻¹⋅mol⁻¹
$(molarenergy(English,SI2019))
```
""" molarenergy

@doc """
$(convertext(:molarconductivity,"conductivity(U,S)*area(U,S)/molaramount(U,S)"))

Conductivity per `molarvolume` or `molarconductivity` (S⋅m²⋅mol⁻¹), unit conversion factor.

```Julia
julia> molarconductivity(EMU,Metric) # S⋅m²⋅abΩ⋅cm⁻²
$(molarconductivity(EMU,Metric))

julia> molarconductivity(ESU,Metric) # S⋅m²⋅statΩ⋅cm⁻²
$(molarconductivity(ESU,Metric))
```
""" molarconductivity

@doc """
$(convertext(:molarsusceptibility,"specificsusceptibility(U,S)*molarmass(U,S)"))

Magnetic/electric molar mass `susceptibility` (m³⋅mol⁻¹), unit conversion factor.

```Julia
julia> molarsusceptibility(CGS,Metric) # m³⋅cm⁻³
$(molarsusceptibility(CGS,Metric))

julia> molarsusceptibility(Metric,SI2019) # m³⋅mol⋅mol⁻¹⋅cm⁻³
$(molarsusceptibility(Metric,SI2019))
```
""" molarsusceptibility

@doc """
$(convertext(:catalysis,"molaramount(U,S)/time(U,S)"))

Catalytic activity or `molaramount` per `time` (kat, mol⋅s⁻¹), unit conversion factor.

```Julia
julia> catalysis(English,Metric) # kat⋅s⋅lb-mol⁻¹
$(catalysis(English,Metric))
```
""" catalysis

@doc """
$(convertext(:specificity,"volume(U,S)/molaramount(U,S)/time(U,S)"))

Catalytic efficiency or `volume` per `mole` per `time` (m³⋅mol⁻¹⋅s⁻¹), unit conversion factor.

```Julia
julia> specificity(CGS,Metric) # m³⋅cm⁻³
$(specificity(CGS,Metric))

julia> specificity(English,Metric) # m³⋅lb-mol⋅mol⁻¹⋅ft⁻³
$(specificity(English,Metric))
```
""" specificity

# photometrics

@doc """
$(convertext(:luminousflux,"luminousenergy(U,S)*frequency(U,S)"))

Perceived power of light or `luminousflux` (lm, cd⋅rad⋅²), unit conversion factor.
""" luminousflux, J

@doc """
$(convertext(:luminousintensity,"luminousflux(U,S)/solidangle(U,S)"))

Perceived power of light or `luminousintensity` (cd, lm⋅rad⁻²), unit conversion factor.
""" luminousintensity

@doc """
$(convertext(:luminance,"luminousintensity(U,S)/area(U,S)"))

Luminous intensity per `area` or `luminance` (cd⋅m⁻², lm⋅m⁻²⋅rad⁻²), unit conversion factor.

```Julia
julia> luminance(CGS,Metric) # lx⋅ph⁻¹
$(luminance(CGS,Metric))

julia> luminance(IAU,Metric) # lx⋅au²⋅lm⁻¹
$(luminance(IAU,Metric))

julia> luminance(English,Metric) # ft²⋅m⁻²
$(luminance(English,Metric))

julia> 1/10.76 # lx⋅fc⁻¹
$(1/10.76)
```
""" luminance

@doc """
$(convertext(:illuminance,"luminousflux(U,S)/area(U,S)"))

Luminous flux per `area` or `luminance` (lx, lm⋅m⁻², cd⋅m⁻²⋅rad²), unit conversion factor.

```Julia
julia> illuminance(CGS,Metric) # lx⋅ph⁻¹
$(illuminance(CGS,Metric))

julia> illuminance(IAU,Metric) # lx⋅au²⋅lm⁻¹
$(illuminance(IAU,Metric))

julia> illuminance(English,Metric) # ft²⋅m⁻²
$(illuminance(English,Metric))

julia> 1/10.76 # lx⋅fc⁻¹
$(1/10.76)
```
""" illuminance

@doc """
$(convertext(:luminousenergy,"luminousflux(U,S)*time(U,S)"))

Perceived quantity of light or `luminousenergy` (lm⋅s, cd⋅s⋅sr), unit conversion factor.

```Julia
julia> luminousenergy(IAU,Metric) # s⋅day⁻¹
$(luminousenergy(IAU,Metric))
```
""" luminousenergy

@doc """
$(convertext(:luminousexposure,"illuminance(U,S)*time(U,S)"))

Integrated `luminance` along `time` (lx⋅s, lm⋅s⋅m⁻², cd⋅s⋅m⁻²⋅sr), unit conversion factor.

```Julia
julia> luminousexposure(CGS,Metric) # lx⋅ph⁻¹
$(luminousexposure(CGS,Metric))

julia> luminousexposure(IAU,Metric) # s⋅au²⋅day⁻¹⋅m⁻²
$(luminousexposure(IAU,Metric))

julia> luminousexposure(English,Metric) # ft²⋅m⁻²
$(luminousexposure(English,Metric))
```
""" luminousexposure

@doc """
$(convertext(:luminousefficacy,"luminousefficacy(S)/luminousefficacy(U)"))

Ratio of `luminousflux` to `power` or `luminousefficacy` (lm⋅W⁻¹), unit conversion factor.

```Julia
julia> luminousefficacy(CGS,Metric) # erg⋅s⁻¹⋅W⁻¹
$(luminousefficacy(CGS,Metric))

julia> luminousefficacy(English,Metric) # ft⋅lb⋅s⁻¹⋅W⁻¹
$(luminousefficacy(English,Metric))
```
""" luminousefficacy(::UnitSystem,::UnitSystem)


text:
using Base.Order: Ordering
import Base.Order: lt

# Here is just some helper functions to get type stable order functors, cause I don't like
# what Julia base has to offer at the moment.

abstract type Target end

"""
For finding eigenvalues with the Largest Magnitude
"""
struct LM <: Target end

"""
For finding eigenvalues with the Largest Real part
"""
struct LR <: Target end

"""
For finding eigenvalues with the Smallest Real part
"""
struct SR <: Target end

"""
For finding eigenvalues with the Largest Imaginary part
"""
struct LI <: Target end

"""
For finding eigenvalues with the Smallest Imaginary part
"""
struct SI <: Target end

struct OrderPerm{TV<:AbstractVector,O<:Ordering} <: Ordering
    xs::TV
    ord::O
end

struct OrderLt{Tf} <: Ordering
    less::Tf
end

struct OrderBy{Tf,O<:Ordering} <: Ordering
    f::Tf
    ord::O
end

struct OrderReverse{O<:Ordering} <: Ordering
    ord::O
end

const Forward = OrderLt(isless)
const Backward = OrderReverse(Forward)

OrderBy(f) = OrderBy(f, Forward)

@inline lt(o::OrderLt, a, b) = o.less(a, b)
@inline lt(o::OrderBy, a, b) = lt(o.ord, o.f(a), o.f(b))
@inline lt(o::OrderReverse, a, b) = lt(o.ord, b, a)

@inline function lt(o::OrderPerm, i::Integer, j::Integer)
    @inbounds fst = o.xs[i]
    @inbounds snd = o.xs[j]

    # First compare by value, then by key to ensure stable order.
    return lt(o.ord, fst, snd) ? true : (lt(o.ord, snd, fst) ? false : i < j)
end

# Helper functions
get_order(which::LM) = OrderReverse(OrderBy(abs))
get_order(which::LR) = OrderReverse(OrderBy(real))
get_order(which::SR) = OrderBy(real)
get_order(which::LI) = OrderReverse(OrderBy(imag))
get_order(which::SI) = OrderBy(imag)


text:
module CovidPOMDPs

using Reexport
using Random
using ProgressMeter
@reexport using POMDPs
@reexport using POMDPModelTools
@reexport using Distributions
@reexport using Statistics
@reexport using ParticleFilters
using DelimitedFiles
using DataFrames: DataFrame
@reexport using CairoMakie
@reexport using BasicPOMCP

include("init.jl")

include("constants.jl")

include("typedef.jl")
export CovidActionSpace, CovidAction

include("simulate.jl")
export SimHist

include("POMDPsInterface.jl")

include("proportional_control.jl")

include("updater.jl")

include("covid_filter.jl")

include("plots.jl")

include(joinpath("single_strain", "SingleStrain.jl"))
using .SingleStrain
export SingleCovidPOMDP, SingleCovidState

include(joinpath("double_strain", "DoubleStrain.jl"))
using .DoubleStrain
export DoubleCovidPOMDP, DoubleCovidState

export susceptible, infected, recovered

end


text:
using VoronoiCells
using GeometryBasics
import VoronoiDelaunay
using Test

@testset "Points" begin
    @testset "Arithmetic with points" begin
        p = VoronoiCells.IndexablePoint2D(1.0, 1.0)
        q = VoronoiCells.IndexablePoint2D(1.1, 1.1)

        @test p + q == VoronoiDelaunay.Point2D(2.1, 2.1)
        @test q - p ≈ VoronoiDelaunay.Point2D(0.1, 0.1)
        @test 1.5 * p == VoronoiDelaunay.Point2D(1.5, 1.5)
    end

    @testset "Distance between points" begin
        p = VoronoiCells.IndexablePoint2D(1.0, 1.0)
        q = VoronoiCells.IndexablePoint2D(3.0, 2.0)

        @test abs2(p, q) == 5
    end
end


text:
module C
const NAMES = [
    # Heldt et al., PNAS (2018)
    "tRb"             #AU  total level of Rb
    "Skp2"            #-  relative SCF^Skp2 level
    "kSyE2f"          #AU/min  constitutive E2F synthesis
    "kSyE2fE2f"       #AU/min  E2F-dependent E2F synthesis
    "jSyE2f"          #AU  Michealis-Menten F_P for E2F synthesis
    "kAsRbE2f"        #1/(AU*min)  association of Rb and E2F
    "kDsRbE2f"        #1/min  dissociation of Rb:E2F complexes
    "kDeE2f"          #1/min  constiutive E2F degradation
    "kPhRbCd"         #1/(AU*min)  CycD:CDK4/6-mediated Rb phosphorylation
    "kPhRbCe"         #1/(AU*min)  CycE:CDK2-mediated Rb phosphorylation
    "kDpRb"           #1/min  dephosphorylation of Rb
    "kDeP21"          #1/min  constitutive p21 degradation
    "kDeP21Cy"        #1/(AU*min)  cyclin:CDK2-mediated p21 degradation (Skp2-dependent)
    "kSyCe"           #1/min  constitutive CycE synthesis
    "kAsCyP21"        #1/(AU*min)  association of cyclin:CDK2 and p21
    "kDsCyP21"        #1/min  dissociation of cyclin:CDK2:p21 complexes
    "kDeCe"           #1/min  constitutive CycE degradation
    "tP21"
    # Added
    "mitogen"
    "pAkt"
    "DUSP"
    "kfRSK"
    "krRSK"
    "kSyCd"
    "kDeCd"
    "kon"
    "koff"
    "kf"
    "kr"
    "tERK"
    "KmRSK"
    "kfFOS"
    "krFOS"
    "KmFOS"
]

#name2idx
for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module


text:
# FASTA
# =====
#
# Reader and writer of the FASTA file format.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

immutable FASTA <: FileFormat end

"Metadata for FASTA sequence records containing just a `description` field"
type FASTAMetadata
    description::StringField
end

function FASTAMetadata()
    return FASTAMetadata(StringField())
end

@compat function Base.:(==)(a::FASTAMetadata, b::FASTAMetadata)
    return a.description == b.description
end

function Base.copy(metadata::FASTAMetadata)
    return FASTAMetadata(copy(metadata.description))
end

"""
    FASTASeqRecord(name, seq[, description=""])

Create a sequence record for the FASTA file format.
"""
typealias FASTASeqRecord{S} SeqRecord{S,FASTAMetadata}

@compat function (::Type{FASTASeqRecord})(name::AbstractString,
                                          seq::Sequence,
                                          description::AbstractString="")
    return SeqRecord(name, seq, FASTAMetadata(description))
end

function Base.open(filepath::AbstractString, mode::AbstractString, ::Type{FASTA};
                   width::Integer=60)
    io = open(filepath, mode)
    if mode[1] == 'r'
        return open(BufferedInputStream(io), FASTA)
    elseif mode[1] ∈ ('w', 'a')
        return FASTAWriter(io, width)
    end
    error("invalid open mode")
end

function Base.open{S}(input::BufferedInputStream, ::Type{FASTA},
                      ::Type{S}=BioSequence)
    return FASTAParser{S}(input)
end

# Parser
# ------

"A type encapsulating the current state of a FASTA parser"
type FASTAParser{S<:Sequence} <: AbstractParser
    state::Ragel.State
    seqbuf::BufferedOutputStream{BufferedStreams.EmptyStream}

    function FASTAParser(input::BufferedInputStream)
        return new(Ragel.State(fastaparser_start, input),
                   BufferedOutputStream())
    end
end

Base.eltype{S}(::Type{FASTAParser{S}}) = FASTASeqRecord{S}
Base.eof(parser::FASTAParser) = eof(parser.state.stream)

include("fasta-parser.jl")

# Writer
# ------

# Serializer for the FASTA file format.
type FASTAWriter{T<:IO} <: AbstractWriter
    output::T
    # maximum sequence width (no limit when width ≤ 0)
    width::Int
end

function Base.flush(writer::FASTAWriter)
    # TODO: This can be removed on Julia v0.5
    # (because flush will be defined for IOBuffer).
    if applicable(flush, writer.output)
        flush(writer.output)
    end
end

Base.close(writer::FASTAWriter) = close(writer.output)

function Base.write(writer::FASTAWriter, seqrec::FASTASeqRecord)
    output = writer.output
    n = 0

    # header
    n += write(output, '>', seqrec.name)
    if !isempty(seqrec.metadata.description)
        n += write(output, ' ', seqrec.metadata.description)
    end
    n += write(output, '\n')

    # sequence
    w = writer.width
    for x in seqrec.seq
        if writer.width > 0 && w == 0
            n += write(output, '\n')
            w = writer.width
        end
        n += write(output, Char(x))
        w -= 1
    end
    n += write(output, '\n')

    return n
end

function Base.show(io::IO, seqrec::FASTASeqRecord)
    print(io, ">", seqrec.name)
    if !isempty(seqrec.metadata.description)
        print(io, " ", seqrec.metadata.description)
    end
    print(io, '\n')
    print(io, seqrec.seq)
end

# This function is almost deprecated in favor of FASTAWriter.
function Base.write{T}(io::IO, seqrec::SeqRecord{T,FASTAMetadata})
    write(io, ">", seqrec.name)
    if !isempty(seqrec.metadata.description)
        write(io, " ", seqrec.metadata.description)
    end
    write(io, "\n")
    maxchars = 79
    counter = 1
    len = length(seqrec.seq)
    for nt in seqrec.seq
        print(io, nt)
        if counter % maxchars == 0 && counter < len
            write(io, "\n")
        end
        counter += 1
    end
    write(io, "\n")
end


text:
function missing_to_zeros!(X)
    for i in eachindex(X)
        if ismissing(X[i])
            X[i] = 0
        end
    end
    X
end

my_mean(μ, n, x)          = μ + (x - μ) / (n + 1), n + 1
my_mean(μ, n, x::Missing) = μ, n
function agg_mean(x)
    μ, n = zero(eltype(x)), 0
    @inbounds for xi in x
        μ, n = my_mean(μ, n, xi)
    end
    return μ
end

function missing_to_mean(X :: Matrix{Union{T, Missing}}) where T

    # If there are missing values here, then there are too many!
    # this causes quite a lot of allocations, but is still pretty fast
    m = [mean(skipmissing(view(X, :, i))) for i in 1:size(X, 2)]
    # much less allocations, but slow:
    # m = [agg_mean(view(X, :, i)) for i in 1:size(X, 2)]
    # this one is way to slow:
    # m = mapslices(x -> mean(skipmissing(x)), X, dims = 1)
    m2 = mean(m)

    X2 = Array{T}(undef, size(X))
    @inbounds for i in eachindex(X, X2)
        X2[i] = ismissing(X[i]) ? m2 : X[i]
    end

    return X2, m, m2
end

# eachslice requires Julia 1.1!
function missing_to_mean_slice(X :: Matrix{Union{T, Missing}}, dims = 1) where T

    X2 = Matrix{T}(undef, size(X))
    for (s1, s2) in zip(eachslice(X, dims = dims), eachslice(X, dims = dims))
        m = s1 |> skipmissing |> mean
        for i in eachindex(s1, s2)
            @inbounds s1[i] = ismissing(s2[i]) ? m : s2[i]
        end
    end

    X2
end

function remove_random!(
    x::Array{Union{T, Missing}},
    p = 0.05
) where {T <: AbstractFloat}
    for i in eachindex(x)
        if rand(T) < T(p)
            x[i] = missing
        end
    end
    x
end


text:
using DataDeps

register(DataDep(
  "kin8nm",
  "https://www.openml.org/d/189",
  "https://www.openml.org/data/get_csv/3626/dataset_2175_kin8nm.arff",
  "7b9bf0301ac936d88122557a151e1ba8f1ebc278fcf46d9f3c6d462debdbc8ad",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=9,
        feature_indices=1:8,
        header=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module TestCycle3x7

import GraphIdx
import GraphIdx: GraphViz, Const, Graph
import GraphIdx.Grid: GridGraph
import TreeLas: Cycle
import TreeLas: CycleGap

function show(y, lambda = Const(1.0))
    graph = GridGraph(size(y)...)
    cmem = CycleGap.CycMem(y, graph, lambda);
    mem = cmem.gap_mem;
    CycleGap.gap_vec!(mem.gamma, mem.x, mem.alpha, mem.wgraph, -1.0)
    CycleGap.prim_mst_edges(mem.gamma, mem.tree.root, mem.mst)
    cmem.cycles = Cycle.CycleBasis(graph, mem.mst.parent)
    tlam0, tlam = Cycle.extract_rotate(cmem.cycles, lambda)
    cycles = cmem.cycles

    GraphViz.dot(graph, cycles, i -> "[ label=\"$(Int(tlam[i]))\"]")
end

@info "compiled"
@sync for size in [(3, 7), (4, 6), (4, 6), (2, 30), (2, 30), (2, 30)]
    @async show(randn(size))
end

end#module

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
"""
    truncate_sigfigs(data, digits = 8)

Cuts off more than digits significant figures.
Used to reduce tree training degeneracies.
"""
function truncate_sigfigs(data, digits = 7)
    return round.(data, sigdigits = digits)
end

""" Computes the bounds of the union of variables in a BlackBoxLearner and its LinkedLearners. """
function get_joint_bounds(bbl::BlackBoxLearner)
    !isempty(bbl.lls) || throw(OCTHaGOnException("BlackBoxLearner $(bbl.name) must have non-zero LinkedLearners for computing joint bounds."))
    bounds = Dict(var => [] for var in bbl.vars)
    bbl_bounds = flattened_bounds(bbl)
    ll_bounds = [flattened_bounds(ll) for ll in bbl.lls]
    for i = 1:length(bbl.vars)
        var = bbl.vars[i]
        var_bound = [
            minimum([minimum(bbl_bounds[i]),
                    [minimum(ll_bound[i]) for ll_bound in ll_bounds]...])
            maximum([maximum(bbl_bounds[i]),
                    [maximum(ll_bound[i]) for ll_bound in ll_bounds]...])]
        bounds[var] = var_bound
    end
    return bounds
end

"""
    $(TYPEDSIGNATURES)

Uniformly Latin Hypercube samples the variables of GlobalModel, as long as all
lbs and ubs are defined.
"""
function lh_sample(bounds::Dict; lh_iterations::Int64 = 0,
                   n_samples::Int64 = 1000)
    check_bounds(bounds)
    vars = flat(keys(bounds))
    n_dims = length(vars)
    if lh_iterations > 0
        plan, _ = LHCoptim(n_samples, n_dims, lh_iterations);
    else
        plan = rand(n_samples, n_dims)
    end
   X = scaleLHC(plan,[(minimum(bounds[var]), maximum(bounds[var])) for var in vars]);
   return DataFrame(truncate_sigfigs(X), string.(vars))
end

function lh_sample(bbl::BlackBoxLearner; lh_iterations::Int64 = 0,
                   n_samples::Int64 = 1000)
    if isempty(bbl.lls)
        bounds = get_bounds(bbl.vars)
        return lh_sample(bounds; lh_iterations = lh_iterations, n_samples = n_samples)
    else
        bounds = get_joint_bounds(bbl)
        return lh_sample(bounds; lh_iterations = lh_iterations, n_samples = n_samples)
    end
end

function choose(large::Int64, small::Int64)
    return Int64(factorial(big(large)) / (factorial(big(large-small))*factorial(big(small))))
end

"""
    $(TYPEDSIGNATURES)

Samples a BlackBoxLearner on the corners of the variable hypercube.
Samples a subset of corners for learners with large number of variables.
"""
function boundary_sample(bounds::Dict; n_samples::Int64 = 100, fraction::Float64 = 0.5,
                         warn_string::String = "")
    check_bounds(bounds);
    vars = flat(keys(bounds))
    n_vars = length(vars);
    vks = string.(vars);
    lbs = DataFrame(Dict(string(key) => minimum(val) for (key, val) in bounds))
    ubs = DataFrame(Dict(string(key) => maximum(val) for (key, val) in bounds))
    n_comb = 0
    if n_vars <= 12
        n_comb = sum(choose(n_vars, i) for i=0:n_vars);
    else
        n_comb = 1e6
    end
    nX = DataFrame(vks .=> [Float64[] for i in vks])
    sample_indices = [];
    if n_comb >= fraction*n_samples
        n_comb = 2*n_vars+2; # Everything is double because we choose min's and max's
        choosing = 1;
        while n_comb <= fraction*n_samples
            choosing = choosing + 1;
            n_comb += 2*choose(n_vars, choosing);
        end
        choosing = choosing - 1; # Determined maximum 'choose' coefficient
        sample_indices = reduce(vcat,collect(combinations(1:n_vars,i)) for i=0:choosing); # Choose 1 and above
    else
        sample_indices = reduce(vcat,collect(combinations(1:n_vars,i)) for i=0:n_vars); # Choose 1 and above
    end
    for i in sample_indices
        lbscopy = copy(lbs); ubscopy = copy(ubs);
        lbscopy[:, vks[i]] = ubscopy[:, vks[i]];
        append!(nX, lbscopy);
        lbscopy = copy(lbs); ubscopy = copy(ubs);
        ubscopy[:, vks[i]] = lbscopy[:, vks[i]];
        append!(nX, ubscopy);
    end
    return nX # Note: boundary samples are not truncated for significant figures
end

function boundary_sample(bbl::BlackBoxLearner; fraction::Float64 = 0.1)
    if isempty(bbl.lls)
        return boundary_sample(get_bounds(bbl.vars), n_samples = get_param(bbl, :n_samples), fraction = fraction,
                           warn_string = bbl.name)
    else
        bounds = get_joint_bounds(bbl)
        return boundary_sample(bounds, n_samples = get_param(bbl, :n_samples), fraction = fraction,
                            warn_string = bbl.name)
    end
end

"""
    knn_sample(bbl::BlackBoxClassifier; k::Int64 = 10, sample_density = 1e-5, sample_idxs = nothing)

Implements KNN sampling on a BlackBoxLearner.
Note: must have at least one in/feasible point for the algorithm.
"""
function knn_sample(bbl::BlackBoxClassifier; k::Int64 = 10, sample_density = 1e-5, sample_idxs = nothing, sign = 1)
    if bbl.feas_ratio == 0. || bbl.feas_ratio == 1.0
        throw(OCTHaGOnException("Constraint " * string(bbl.name) * " must have at least one feasible or
                            infeasible sample to be KNN-sampled!"))
    end
    vks = string.(bbl.vars)
    df = DataFrame(vks .=> [Float64[] for i in vks])
    build_knn_tree(bbl);
    idxs, dists = find_knn(bbl, k=k);
    feas_class = classify_patches(bbl, idxs);
    negatives = findall(x -> x .< 0, bbl.Y)
    if !isnothing(sample_idxs)
        negatives = intersect(negatives, sample_idxs)
    end
    for i = 1:length(negatives) # This loop is for making sure that every possible root is sampled only once.
        if feas_class[negatives[i]] == "mixed"
            nodes = [idxs[negatives[i]][j] for j=1:length(idxs[negatives[i]])
                            if (bbl.Y[idxs[negatives[i]][j]] >= 0 && dists[negatives[i]][j] >= sample_density)]
            push!(nodes, negatives[i])
            np = secant_method(bbl.X[nodes, :], bbl.Y[nodes, :])
            append!(df, np)
        end
    end
    return truncate_sigfigs(df)
end

"""
    $(TYPEDSIGNATURES)

Samples and evaluates nonlinear constraints using full suite of methods.
"""
function uniform_sample_and_eval!(bbl::BlackBoxLearner;
                          boundary_fraction::Float64 = 0.1,
                          lh_iterations::Int64 = 0, sample_density = 1.e-5)
    @assert size(bbl.X, 1) == 0 # TODO: fix this w.r.t. data-driven constraints.
    vks = string.(bbl.vars)
    n_dims = length(vks);
    check_bounds(get_bounds(bbl))
    df = boundary_sample(bbl, fraction = boundary_fraction)
    eval!(bbl, df)
    df = lh_sample(bbl, lh_iterations = lh_iterations, n_samples = get_param(bbl, :n_samples) - size(df, 1))
    eval!(bbl, df);
    if bbl isa BlackBoxClassifier
        if bbl.feas_ratio == 1.0
            @info(string(bbl.name) * " was not KNN sampled since it has no infeasible samples.")
        elseif bbl.feas_ratio == 0.0
            throw(OCTHaGOnException(string(bbl.name) * " has zero feasible samples. " *
                               "Please find at least one feasible sample, seed the data and KNN sample."))
        else
            df = knn_sample(bbl, k= maximum([10, 2*length(bbl.vars) + 1]), sample_density = sample_density)
            if size(df, 1) > 0
                eval!(bbl, df)
            end
        end
    end
    bbl.max_Y = maximum(filter(!isinf, bbl.Y))
    bbl.min_Y = minimum(filter(!isinf, bbl.Y))
    # TODO: update these mx/min Y more frequently.
    # Setting vague bounds for the dependent variable.
    # This is the only big-M required in the formulation, adding 50% margin to all observed Y values.
    if bbl isa BlackBoxRegressor
        max_Y = bbl.max_Y
        min_Y = bbl.min_Y
        lower_margined_bound = min_Y - (max_Y - min_Y)/2
        upper_margined_bound = max_Y + (max_Y - min_Y)/2
        if !JuMP.has_lower_bound(bbl.dependent_var)
            JuMP.set_lower_bound(bbl.dependent_var, lower_margined_bound)
        end
        if !JuMP.has_upper_bound(bbl.dependent_var)
            JuMP.set_upper_bound(bbl.dependent_var, upper_margined_bound)
        end
        for ll in bbl.lls
            if !JuMP.has_lower_bound(ll.dependent_var)
                JuMP.set_lower_bound(ll.dependent_var, lower_margined_bound)
            end
            if !JuMP.has_upper_bound(ll.dependent_var)
                JuMP.set_upper_bound(ll.dependent_var, upper_margined_bound)
            end
        end
    end
    return
end

function uniform_sample_and_eval!(bbls::Array{BlackBoxLearner}; lh_iterations = 0, sample_density = 1e-5)
    for bbl in bbls
        uniform_sample_and_eval!(bbl, lh_iterations = lh_iterations, sample_density = sample_density)
    end
    return
end

function uniform_sample_and_eval!(gm::GlobalModel)
    @info "Sampling..."
    uniform_sample_and_eval!(gm.bbls; lh_iterations = get_param(gm, :lh_iterations), sample_density = get_param(gm, :sample_density))
    return
end

"""
    $(TYPEDSIGNATURES)

Gets Latin Hypercube samples that fall in the leaf of the last solution.
"""
function last_leaf_sample(bbc::BlackBoxClassifier, n_samples = get_param(bbc, :n_samples))
    isempty(bbc.active_leaves) &&
        throw(OCTHaGOnException("BBC $(bbc.name) needs to be optimized first, to figure out its active leaves."))
    if !bbc.equality
        last_leaf = bbc.active_leaves[1]
        idxs = findall(x -> x .>= 0.5, IAI.apply(bbc.learners[end], bbc.X) .== last_leaf)
        lbs = [minimum(col) for col in eachcol(bbc.X[idxs, :])]
        ubs = [maximum(col) for col in eachcol(bbc.X[idxs, :])]
        plan, _ = LHCoptim(n_samples, length(bbc.vars), 3);
        X = scaleLHC(plan, [(lbs[i], ubs[i]) for i=1:length(lbs)]);
        return DataFrame(truncate_sigfigs(X), string.(bbc.vars))
    else
        for last_leaf in bbc.active_leaves
            idxs = findall(x -> x .>= 0.5, IAI.apply(bbc.learners[end], bbc.X) .== last_leaf)
            lbs = [minimum(col) for col in eachcol(bbc.X[idxs, :])]
            ubs = [maximum(col) for col in eachcol(bbc.X[idxs, :])]
            plan, _ = LHCoptim(Int(ceil(n_samples/2)), length(bbc.vars), 3);
            X = scaleLHC(plan, [(lbs[i], ubs[i]) for i=1:length(lbs)]);
            return DataFrame(truncate_sigfigs(X), string.(bbc.vars))
        end
    end
end

function last_leaf_sample(bbr::BlackBoxRegressor, n_samples = get_param(bbr, :n_samples))
    if isempty(bbr.active_leaves)
        throw(OCTHaGOnException("BBR $(bbr.name) needs to be optimized first, to figure out its active leaves."))
    elseif length(bbr.active_trees) == 2
        upper_leaf, lower_leaf = sort(bbr.active_leaves)
        upper_leafneighbor = [];
        lower_leafneighbor = [];
        for (tree_idx, threshold) in bbr.active_trees
            if threshold.first == "upper"
                append!(upper_leafneighbor, IAI.apply(bbr.learners[tree_idx], bbr.X) .== -upper_leaf)
            elseif threshold.first == "lower"
                append!(lower_leafneighbor, IAI.apply(bbr.learners[tree_idx], bbr.X) .== lower_leaf)
            end
        end
        idxs =  findall(x -> x .>= 0.5, upper_leafneighbor .* lower_leafneighbor)
        if length(idxs) == 0
            @info("No points in $(bbr.name) in the intersection of trees. Widening sampling. ")
            idxs =  findall(x -> x .>= 0.5, upper_leafneighbor + lower_leafneighbor)
        end
        lbs = [minimum(col) for col in eachcol(bbr.X[idxs, :])]
        ubs = [maximum(col) for col in eachcol(bbr.X[idxs, :])]
        plan, _ = LHCoptim(n_samples, length(bbr.vars), 3);
        X = scaleLHC(plan, [(lbs[i], ubs[i]) for i=1:length(lbs)]);
        return DataFrame(truncate_sigfigs(X), string.(bbr.vars))
    elseif length(bbr.active_trees) == 1
        leaf = bbr.active_leaves[1]
        tree = bbr.learners[collect(keys(bbr.active_trees))[1]]
        idxs = findall(x -> x .== leaf, IAI.apply(tree, bbr.X))
        lbs = [minimum(col) for col in eachcol(bbr.X[idxs, :])]
        ubs = [maximum(col) for col in eachcol(bbr.X[idxs, :])]
        plan, _ = LHCoptim(n_samples, length(bbr.vars), 3);
        X = scaleLHC(plan, [(lbs[i], ubs[i]) for i=1:length(lbs)]);
        return DataFrame(truncate_sigfigs(X), string.(bbr.vars))
    else
        throw(OCTHaGOnException("No active trees found in BBR $(bbr.name) while attempting to leaf-sample."))
    end
end

"""
    feasibility_sample(gm::GlobalModel)
    feasibility_sample(bbc::BlackBoxClassifier, n_samples::Int64 = get_param(bbc, :n_samples))

Equally LH samples the leaves of BBC with not enough feasible samples.
"""
function feasibility_sample(bbc::BlackBoxClassifier, n_samples::Int64 = get_param(bbc, :n_samples))
    minbucket = minimum([Int64(round(bbc.feas_ratio .* size(bbc.X, 1))), Int64(floor(0.05*size(bbc.X, 1)))])
    lnr = base_classifier()
    IAI.set_params!(lnr, minbucket = minbucket)
    lnr = learn_from_data!(bbc.X, bbc.Y .>= 0, lnr; fit_classifier_kwargs()...)
    all_leaves = find_leaves(lnr)
    feas_leaves = [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))]
    orig_feasratio = copy(bbc.feas_ratio)
    leaf_idxs = IAI.apply(lnr, bbc.X)
    for leaf in feas_leaves
        idxs = findall(x -> x .>= 0.5, leaf_idxs .== leaf)
        lbs = [minimum(col) for col in eachcol(bbc.X[idxs, :])]
        ubs = [maximum(col) for col in eachcol(bbc.X[idxs, :])]
        plan, _ = LHCoptim(Int64(ceil(n_samples / length(feas_leaves))), length(bbc.vars), 3);
        X = scaleLHC(plan, [(lbs[i], ubs[i]) for i=1:length(lbs)]);
        eval!(bbc, DataFrame(X, string.(bbc.vars)))
    end
    if bbc.feas_ratio >= get_param(bbc, :threshold_feasibility)
        @info("BBC $(bbc.name) has passed threshold feasibility $(round(get_param(bbc, :threshold_feasibility), sigdigits=3)), " *
              "from $(round(orig_feasratio, sigdigits=3)) to $(round(bbc.feas_ratio, sigdigits=3)).")
    elseif bbc.feas_ratio > orig_feasratio
        @info("BBC $(bbc.name) has improved feasibility from $(round(orig_feasratio, sigdigits=3)) to $(round(bbc.feas_ratio, sigdigits=3)), " *
                "but still not over threshold $(round(get_param(bbc, :threshold_feasibility), sigdigits=3)).")
    else
        @info("BBC $(bbc.name) feasibility improvement failed! Please check your functions.")
    end
    return
end

function feasibility_sample(gm::GlobalModel)
    for bbl in gm.bbls
        if bbl isa BlackBoxClassifier && bbl.feas_ratio <= get_param(bbl, :threshold_feasibility)
            feasibility_sample(bbl)
        end
    end
    return
end

"""
    upper_bound_sample(bbr::BlackBoxRegressor)

Samples within the feasible regions of the active upper bounding tree.
"""
function upper_bound_sample(bbr::BlackBoxRegressor)
    @assert !isempty(bbr.active_trees) # We must have optimized at least once over the old tree.
    active_upper_idx = active_upper_tree(bbr)
    lnr = bbr.learners[active_upper_idx]
    @assert lnr isa IAI.OptimalTreeClassifier
    all_leaves = find_leaves(lnr)
    feas_leaves = [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))]
    leaf_idxs = IAI.apply(lnr, bbr.X)
    n_samples = get_param(bbr, :n_samples)
    for leaf in feas_leaves
        idxs = findall(x -> x == leaf, leaf_idxs)
        lbs = [minimum(col) for col in eachcol(bbr.X[idxs, :])]
        ubs = [maximum(col) for col in eachcol(bbr.X[idxs, :])]
        plan, _ = LHCoptim(Int64(ceil(n_samples / length(feas_leaves))), length(bbr.vars), 3);
        X = scaleLHC(plan, [(lbs[i], ubs[i]) for i=1:length(lbs)]);
        eval!(bbr, DataFrame(X, string.(bbr.vars)))
    end
    return
end


text:
"""
    step!(mechanism::Mechanism{T}, z::Vector{T}, u::Vector{T}; opts)

    simulate mechanism for one time step

    mechanism: Mechanism
    z: maximal state
    u: inputs
    opts: SolverOptions
"""
function step!(mechanism::Mechanism{T}, z::Vector{T}, u::Vector{T};
    opts=SolverOptions{T}()) where T

    # set state
    set_maximal_state!(mechanism, z)

    # set control
    set_input!(mechanism, u)

    # solve the 1-step simulation problem
    mehrotra!(mechanism, opts=opts)

    # extract the next state
    z_next = get_next_state(mechanism)

    return z_next
end


text:
using ACE, StaticArrays, BenchmarkTools, Printf

using Zygote: gradient
using ACE: evaluate, evaluate!, evaluate_d, evaluate_d!
using ACEbase: acquire_B!, acquire_dB!

TX = ACE.PositionState{Float64}
B1p = ACE.Utils.RnYlm_1pbasis()
Rn = B1p.bases[1]
cfg = ACEConfig(rand(TX, Rn, 30))

##

# degrees = Dict(2 => [9, 17],
#                3 => [7, 15] )

Pnumprops = [1, 2, 4]

degrees = Dict(2 => [10, 17],
               3 => [9, 15],
               5 => [7, 11] )

wL = 1.5

#zygote tests
site_energy(m,x) = sum(ACE.val.(ACE.evaluate(m, x)))

#loss tests
FS = props -> sum( (1 .+ ACE.val.(props).^2).^0.5 )
sqr(x) = x.rr .^ 2

#full loss
loss(m, x) = (FS(ACE.evaluate(m,x)))^2 + sum(sum(sqr.(gradient(tx->FS(ACE.evaluate(m,tx)), x)[1])))
#only energy loss
lossE(m, x) = (FS(ACE.evaluate(m,x)))^2

##
#Pgroup for benchmarking different number of properties (P)
#fixed maxorder = 2 and maxdeg = 7

fixdeg = 7
fixord = 2

#linear model
Pgroup = BenchmarkGroup()
Pgroup["set_params!"] = BenchmarkGroup()
Pgroup["evaluate"] = BenchmarkGroup()

#only ACE
Pgroup["grad_params"] = BenchmarkGroup()
Pgroup["grad_config"] = BenchmarkGroup()
Pgroup["_rrule_evaluate"] = BenchmarkGroup()
Pgroup["grad_params_config"] = BenchmarkGroup()
Pgroup["adjoint_EVAL_D"] = BenchmarkGroup()

# #zygote calls
# Pgroup["site_energy"] = BenchmarkGroup()
# Pgroup["Zygote_grad_params"] = BenchmarkGroup()
# Pgroup["Zygote_grad_config"] = BenchmarkGroup()
# Pgroup["Zygote_grad_params_config"] = BenchmarkGroup()

# #loss functions
# Pgroup["eval_full_loss"] = BenchmarkGroup()
# Pgroup["eval_energy_loss"] = BenchmarkGroup()
# Pgroup["der_full_loss"] = BenchmarkGroup()
# Pgroup["der_energy_loss"] = BenchmarkGroup()

for numprops in Pnumprops
   local B1p
   Bsel = SparseBasis(; maxorder = fixord, p = 1, default_maxdeg = fixdeg,
                        weight = Dict(:n => 1.0, :l => wL))
   B1p = ACE.Utils.RnYlm_1pbasis(maxdeg = fixdeg,
                                 maxL = ceil(Int, fixdeg / wL),
                                 Bsel = Bsel)

   φ = ACE.Invariant()
   bsis = ACE.SymmetricBasis(φ, B1p, ACE.O3(), Bsel)

   #create a multiple property model
   W = rand(numprops, length(bsis))
   c = [SVector{size(W)[1]}(W[:,i]) for i in 1:size(W)[2]]
   LM = ACE.LinearACEModel(bsis, c, evaluator = :standard)

   Pgroup["set_params!"][numprops] = @benchmarkable ACE.set_params!($LM, $c)
   Pgroup["evaluate"][numprops] = @benchmarkable ACE.evaluate($LM, $cfg)

   #ACE native calls for derivatives
   Pgroup["grad_params"][numprops] = @benchmarkable ACE.grad_params($LM, $cfg)
   Pgroup["grad_config"][numprops] = @benchmarkable ACE.grad_config($LM, $cfg)
   dp = ones(numprops) #the forces pullback input
   Pgroup["_rrule_evaluate"][numprops] = @benchmarkable ACE._rrule_evaluate($dp, $LM, $cfg)
   Pgroup["grad_params_config"][numprops] = @benchmarkable ACE.grad_params_config($LM, $cfg)
   dq = [ACE.DState(rr=rand(SVector{3, Float64})) for _ = 1:length(cfg)] #the adjoint
   Pgroup["adjoint_EVAL_D"][numprops] = @benchmarkable ACE.adjoint_EVAL_D($LM, $LM.evaluator, $cfg, $dq)

   # #Zygote calls for derivatives
   # Pgroup["site_energy"][numprops] = @benchmarkable site_energy($LM, $cfg)
   # Pgroup["Zygote_grad_params"][numprops] = @benchmarkable gradient(x->site_energy(x,$cfg), $LM)
   # Pgroup["Zygote_grad_config"][numprops] = @benchmarkable gradient(x->site_energy($LM,x), $cfg)
   # #sum over the function so we don't compute the jacobian
   # Pgroup["Zygote_grad_params_config"][numprops] = @benchmarkable gradient(m->sum(sum(gradient(x->site_energy(m,x), $cfg)[1]).rr), $LM)

   # #loss function
   # Pgroup["eval_full_loss"][numprops] = @benchmarkable loss($LM, $cfg)
   # Pgroup["eval_energy_loss"][numprops] = @benchmarkable lossE($LM, $cfg)
   # Pgroup["der_full_loss"][numprops] = @benchmarkable gradient(m->loss(m,$cfg), $LM)
   # Pgroup["der_energy_loss"][numprops] = @benchmarkable gradient(m->lossE(m,$cfg), $LM)

end

##
#Benchmarking derivatives for different sizes of basis (B)
#the number of properties is fixed to 2

Nprop = 2

#linear model
Bgroup = BenchmarkGroup()
Bgroup["set_params!"] = BenchmarkGroup()
Bgroup["evaluate"] = BenchmarkGroup()

#only ACE
Bgroup["grad_params"] = BenchmarkGroup()
Bgroup["grad_config"] = BenchmarkGroup()
Bgroup["_rrule_evaluate"] = BenchmarkGroup()
Bgroup["grad_params_config"] = BenchmarkGroup()
Bgroup["adjoint_EVAL_D"] = BenchmarkGroup()

# #zygote calls
# Bgroup["site_energy"] = BenchmarkGroup()
# Bgroup["Zygote_grad_params"] = BenchmarkGroup()
# Bgroup["Zygote_grad_config"] = BenchmarkGroup()
# Bgroup["Zygote_grad_params_config"] = BenchmarkGroup()

# #loss functions
# Bgroup["eval_full_loss"] = BenchmarkGroup()
# Bgroup["eval_energy_loss"] = BenchmarkGroup()
# Bgroup["der_full_loss"] = BenchmarkGroup()
# Bgroup["der_energy_loss"] = BenchmarkGroup()

for ord = keys(degrees), deg in degrees[ord]
   local B1p
   Bsel = SparseBasis(; maxorder = ord, p = 1, default_maxdeg = deg,
                        weight = Dict(:n => 1.0, :l => wL))
   B1p = ACE.Utils.RnYlm_1pbasis(maxdeg = deg,
                                 maxL = ceil(Int, deg / wL),
                                 Bsel = Bsel)

   φ = ACE.Invariant()
   bsis = ACE.SymmetricBasis(φ, B1p, ACE.O3(), Bsel)

   #create a multiple property model
   W = rand(Nprop, length(bsis))
   c = [SVector{size(W)[1]}(W[:,i]) for i in 1:size(W)[2]]
   LM = ACE.LinearACEModel(bsis, c, evaluator = :standard)

   Bgroup["set_params!"][ord, deg] = @benchmarkable ACE.set_params!($LM, $c)
   Bgroup["evaluate"][ord, deg] = @benchmarkable ACE.evaluate($LM, $cfg)

   #ACE native calls for derivatives
   Bgroup["grad_params"][ord, deg] = @benchmarkable ACE.grad_params($LM, $cfg)
   Bgroup["grad_config"][ord, deg] = @benchmarkable ACE.grad_config($LM, $cfg)
   dp = ones(Nprop) #the forces pullback input
   Bgroup["_rrule_evaluate"][ord, deg] = @benchmarkable ACE._rrule_evaluate($dp, $LM, $cfg)
   Bgroup["grad_params_config"][ord, deg] = @benchmarkable ACE.grad_params_config($LM, $cfg)
   dq = [ACE.DState(rr=rand(SVector{3, Float64})) for _ = 1:length(cfg)] #the adjoint
   Bgroup["adjoint_EVAL_D"][ord, deg] = @benchmarkable ACE.adjoint_EVAL_D($LM, $LM.evaluator, $cfg, $dq)

   # #Zygote calls for derivatives
   # Bgroup["site_energy"][ord, deg] = @benchmarkable site_energy($LM, $cfg)
   # Bgroup["Zygote_grad_params"][ord, deg] = @benchmarkable gradient(x->site_energy(x,$cfg), $LM)
   # Bgroup["Zygote_grad_config"][ord, deg] = @benchmarkable gradient(x->site_energy($LM,x), $cfg)
   # #sum over the function so we don't compute the jacobian
   # Bgroup["Zygote_grad_params_config"][ord, deg] = @benchmarkable gradient(m->sum(sum(gradient(x->site_energy(m,x), $cfg)[1]).rr), $LM)

   # #loss function
   # Bgroup["eval_full_loss"][ord, deg] = @benchmarkable loss($LM, $cfg)
   # Bgroup["eval_energy_loss"][ord, deg] = @benchmarkable lossE($LM, $cfg)
   # Bgroup["der_full_loss"][ord, deg] = @benchmarkable gradient(m->loss(m,$cfg), $LM)
   # Bgroup["der_energy_loss"][ord, deg] = @benchmarkable gradient(m->lossE(m,$cfg), $LM)

end

##

linear_suite = BenchmarkGroup()
linear_suite["P"] = Pgroup
linear_suite["B"] = Bgroup

##

# @info("Tune")
# tune!(linear_suite)

# @info("Run")
# results = run(linear_suite, verbose = true)


text:
using DelimitedFiles

function generate_datasets(namedir, application_name, nbpoints)
    """
    This function takes the full dataset and sets the first nbpoints points in a test set,
        then creates an ever growing train dataset by increments of nboints"""

    fulldataset = readdlm("examples/data/consolidated_data_$application_name.csv", ',')
    nb_train = size(fulldataset)[1]

    testdataset = fulldataset[1:nbpoints,:]
    writedlm(namedir*"test_dataset_$application_name.csv", testdataset, ',')

    i = 2
    while nbpoints * i <= nb_train
        writedlm(namedir*"train_dataset_$(application_name)_$(i-1).csv", fulldataset[nbpoints+1:i*nbpoints,:], ',')
        i+=1
    end
end

application_name = "real"
namedir = "examples/data/"
nbpoints = 10000
generate_datasets(namedir, application_name, nbpoints)


text:
using StaticArrays
using SymPy

"""
WENO reconstruction schemes following the lecture notes by Shu (1998) and
the review article by Shu (2009). WENO smoothness indicators are due to
Jiang & Shu (1996).

Shu (1998), "Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
    Schemes for Hyperbolic Conservation Laws", In: Advanced Numerical Approximation
    of Nonlinear Hyperbolic Equations, edited by Cockburn et al., pp. 325–432.
    DOI: https://doi.org/10.1007/BFb0096355

Shu (2009) "High Order Weighted Essentially Nonoscillatory Schemes for Convection
    Dominated Problems", SIAM Review 51(1), pp. 82–126.
    DOI: https://doi.org/10.1137/070679065

Jiang & Shu (1996) "Efficient Implementation of Weighted ENO Schemes", Journal of
    Computational Physics 126, pp. 202–228.
    DOI: https://doi.org/10.1006/jcph.1996.0130
"""

#####
##### ENO reconstruction coefficients and WENO weights on a uniform grid
##### Equation (2.21) from the Shu (1998) lecture notes.
#####

"""
    UΠ(k, r, m, l)

Return the product in the numerator of equation (2.21) of Shu (1998) for the `l`th
term of the `m`th Lagrange basis polynomial of an ENO reconstruction scheme with
order `k` and left shift `r`.
"""
UΠ(k, r, l, m) = prod([r-q+1 for q in 0:k if q ∉ (m, l)])

"""
    U(k, r, m)

Return the numerator in equation (2.21) of Shu (1998) for the `m`th Lagrange basis
polynomial of an ENO reconstruction scheme with order `k` and left shift `r`.
"""
U(k, r, m) = sum([UΠ(k, r, m, l) for l in 0:k if l != m])

"""
    D(k, m)

Return the denominator in equation (2.21) of Shu (1998) for the `m`th Lagrange basis
polynomial of an ENO reconstruction scheme with order `k`.
"""
D(k, m) = prod([m - l for l in 0:k if l != m])

"""
    eno_coefficient(k, r)

Return the `j`th ENO coefficient used to reconstruct a value at the point x(i+½) with
order of accuracy `k` (stencil size) and left shift `r`.
"""
eno_coefficient(k, r, j) = sum([U(k, r, m)//D(k, m) for m in j+1:k])

"""
    eno_coefficients(k, r)

Return an array of ENO coefficients to reconstruct a value at the point x(i+½) with
order of accuracy `k` (stencil size) and left shift `r`.
"""
eno_coefficients(k, r) = [eno_coefficient(k, r, j) for j in 0:k-1]

"""
    eno_coefficients_matrix([FT=Rational], k)

Return a k×k static array containing ENO coefficients to reconstruct a value at the
point x(i+½) with order of accuracy `k` (stencil size) with element type `FT`. Note
that when combined these ENO interpolants produce a WENO scheme of order 2k-1.
"""
eno_coefficients_matrix(FT, k) =
    cat([eno_coefficients(k, r) for r in 0:k-1]..., dims=1) |> SMatrix{k,k,FT}

eno_coefficients_matrix(k) = eno_coefficients_matrix(Rational, k)

"""
    optimal_weno_weights([FT=Rational], k)

Return a static vector containing the optimal weights that can be used to weigh ENO
reconstruction schemes of order `k` to produce a WENO scheme of order 2k-1.
"""
function optimal_weno_weights(FT, k)
    C = zeros(Rational, 2k-1, k)
    b = eno_coefficients(2k-1, k-1)

    for n in 0:k-1
        C[n+1:n+k, n+1] .= eno_coefficients(k, k-1-n)
    end

    γ = C \ b
    γ = rationalize.(γ, tol=√eps(Float64)) |> reverse
    return SVector{k,FT}(γ)
end

optimal_weno_weights(k) = optimal_weno_weights(Rational, k)

#####
##### Jiang & Shu (1996) WENO smoothness indicators β
##### See equation (2.61) of Shu (1998).
#####

x(j) = j  # Assume uniform grid for now.

"""
    ℓ(ξ, k, r, j)

Return a symbolic expression for the `j`th Lagrange basis polynomial ℓ(ξ) for an
ENO reconstruction scheme of order k and left shift `r` with field values `ϕ`.
"""
ℓ(ξ, k, r, j) = prod((ξ - x(m-r)) / (x(j-r) - x(m-r)) for m in 0:k-1 if m != j)

"""
    L(ξ, k, r, ϕ)

Return a symbolic expression for the interpolating Lagrange polynomial p(ξ) for an
ENO reconstruction scheme of order k and left shift `r` with field values `ϕ`.
"""
p(ξ, k, r, ϕ) = sum(ℓ(ξ, k, r, j) * ϕ[j+1] for j in 0:k-1)

"""
    β(k, r, ϕ)

Return a symbolic expression for the smoothness indicator β for a WENO reconstruction
scheme with stencils of size `k` and left shift `r` (WENO scheme of order 2k-1). The
field values are represented by the symbols in `ϕ` which should have length k.
"""
function β(k, r, ϕ)
    @vars ξ
    return sum(integrate(diff(p(ξ, k, r, ϕ), ξ, l)^2, (ξ, Sym(-1//2), Sym(1//2))) for l in 1:k-1)
end

"""
    subscript(n)

Convert the integer `n` to a subscript in the form of a unicode string. Note
that `0x2080` is the unicode encoding for the subscript 0.
"""
subscript(n) = join(Char(0x2080 + parse(Int, d)) for d in string(n))

subscript_sign(n) = n > 0 ? "₊" : n < 0 ? "₋" : ""

subscript_index(n) = n == 0 ? "" : subscript_sign(n) * subscript(abs(n))

"""
    β_coefficients([FT=Rational], k)

Return a k×k×k static array containing the WENO smoothness indicator coefficients
described by Jiang & Shu (1998) for a WENO reconstruction with stencils of size `k`
(WENO scheme of order 2k-1) with element type `FT`.

The `B[m, n, r]` coefficient corresponds to the coefficient of the
`ϕ[r-k+m+1] * ϕ[r-k+n+1]` term where r-k+1 <= m, n <= r and `r` is the left shift of
the ENO interpolants.
"""
function β_coefficients(FT, k)
    B = zeros(Float64, k, k, k)

    for r in 0:k-1
        ϕ = [Sym("ϕᵢ" * subscript_index(n)) for n in r:-1:r-k+1]
        β_symbolic = β(k, r, ϕ) |> expand

        for m in 1:k, n in 1:k
            B[m, n, r+1] = β_symbolic.coeff(ϕ[m] * ϕ[n])
        end
    end

    B = rationalize.(B, tol=√eps(Float64))
    return SArray{Tuple{k,k,k},FT}(B)
end

β_coefficients(k) = β_coefficients(Rational, k)


text:
using Plots, LinearAlgebra, Printf
# helper functions
@views     av(A) = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
@views  av_xa(A) =  0.5*(A[1:end-1,:].+A[2:end,:])
@views  av_ya(A) =  0.5*(A[:,1:end-1].+A[:,2:end])
@views maxloc(A) = max.(A[1:end-2,1:end-2],A[1:end-2,2:end-1],A[1:end-2,3:end],
                        A[2:end-1,1:end-2],A[2:end-1,2:end-1],A[2:end-1,3:end],
                        A[3:end  ,1:end-2],A[3:end  ,2:end-1],A[3:end  ,3:end])
@views   bc2!(A) = begin A[1,:] = A[2,:]; A[end,:] = A[end-1,:]; A[:,1] = A[:,2]; A[:,end] = A[:,end-1] end
# main function
@views function Stokes2D_vep()
    use_vep = true
    # phyiscs
    lx, ly  = 1.0, 1.0
    radi    = 0.1
    η0      = 1.0
    η_reg   = 1.2e-2
    G0      = 1.0
    Gi      = 0.5G0
    τ_y     = 1.6
    sinϕ    = sind(30)
    ebg     = 1.0
    dt      = η0/G0/4.0
    # numerics
    nx, ny  = 63, 63
    nt      = 10
    εnl     = 1e-6
    maxiter = 150max(nx, ny)
    nchk    = 2max(nx, ny)
    Re      = 5π
    r       = 1.0
    CFL     = 0.99/sqrt(2)
    # preprocessing
    dx, dy  = lx/nx, ly/ny
    max_lxy = max(lx, ly)
    vpdτ    = CFL*min(dx, dy)
    xc, yc  = LinRange(-(lx-dx)/2,(lx-dx)/2,nx), LinRange(-(ly-dy)/2,(ly-dy)/2,ny)
    xv, yv  = LinRange(-lx/2,lx/2,nx+1), LinRange(-ly/2,ly/2,ny+1)
    # allocate arrays
    Pr      = zeros(nx  ,ny  )
    τxx     = zeros(nx  ,ny  )
    τyy     = zeros(nx  ,ny  )
    τxy     = zeros(nx+1,ny+1)
    τxyc    = zeros(nx  ,ny  )
    τii     = zeros(nx  ,ny  )
    Eii     = zeros(nx  ,ny  )
    λ       = zeros(nx  ,ny  )
    F       = zeros(nx  ,ny  )
    Fchk    = zeros(nx  ,ny  )
    dQdTxx  = zeros(nx  ,ny  )
    dQdTyy  = zeros(nx  ,ny  )
    dQdTxy  = zeros(nx  ,ny  )
    τxx_o   = zeros(nx  ,ny  )
    τyy_o   = zeros(nx  ,ny  )
    τxyc_o  = zeros(nx  ,ny  )
    τxy_o   = zeros(nx+1,ny+1)
    # Here _r stands for real: physical stress for residual computation
    τxx_r   = zeros(nx  ,ny  )
    τyy_r   = zeros(nx  ,ny  )
    τxyc_r  = zeros(nx  ,ny  )
    τxy_r   = zeros(nx+1,ny+1)
    Vx      = zeros(nx+1,ny  )
    Vy      = zeros(nx  ,ny+1)
    dVx     = zeros(nx-1,ny  )
    dVy     = zeros(nx  ,ny-1)
    Rx      = zeros(nx-1,ny  )
    Ry      = zeros(nx  ,ny-1)
    ∇V      = zeros(nx  ,ny  )
    ρg      = zeros(nx  ,ny  )
    Exx     = zeros(nx  ,ny  )
    Eyy     = zeros(nx  ,ny  )
    Exyc    = zeros(nx  ,ny  )
    Exy     = zeros(nx+1,ny+1)
    Exx_e   = zeros(nx  ,ny  )
    Eyy_e   = zeros(nx  ,ny  )
    Exyc_e  = zeros(nx  ,ny  )
    Exy_e   = zeros(nx+1,ny+1)
    Exx_τ   = zeros(nx  ,ny  )
    Eyy_τ   = zeros(nx  ,ny  )
    Exy_τ   = zeros(nx+1,ny+1)
    Exyc_τ  = zeros(nx  ,ny  )
    η_ve_τ  = zeros(nx  ,ny  )
    η_ve_τv = zeros(nx+1,ny+1)
    η_vem   = zeros(nx  ,ny  )
    η_vev   = zeros(nx+1,ny+1)
    η_vevm  = zeros(nx+1,ny+1)
    dτ_ρ    = zeros(nx  ,ny  )
    dτ_ρv   = zeros(nx+1,ny+1)
    Gdτ     = zeros(nx  ,ny  )
    Gdτv    = zeros(nx+1,ny+1)
    η_vep   = zeros(nx  ,ny  )
    η_vepv  = zeros(nx+1,ny+1)
    η_vec   = zeros(nx  ,ny  )
    # init
    Vx = [ ebg*x for x ∈ xv, _ ∈ yc ]
    Vy = [-ebg*y for _ ∈ xc, y ∈ yv ]
    η  = fill(η0,nx,ny); ηv = fill(η0,nx+1,ny+1)
    G  = fill(G0,nx,ny); Gv = fill(G0,nx+1,ny+1)
    @. G[xc^2 + yc'^2 < radi^2]  = Gi
    @. Gv[xv^2 + yv'^2 < radi^2] = Gi
    η_e   = G.*dt; η_ev = Gv.*dt
    @. η_vec = 1.0/(1.0/η  + 1.0/η_e)
    @. η_vev = 1.0/(1.0/ηv + 1.0/η_ev)
    @. η_vep  = 1.0/(1.0/η  + 1.0/η_e)
    @. η_vepv = 1.0/(1.0/ηv + 1.0/η_ev)
    # action
    t = 0.0; evo_t = Float64[]; evo_τxx = Float64[]; niter = 0
    for it = 1:nt
        τxx_o .= τxx; τyy_o .= τyy; τxy_o .= τxy; τxyc_o .= τxyc
        err  = 2εnl; iter = 0
        while err > εnl && iter < maxiter
            if !use_vep
                η_vem[2:end-1,2:end-1]  .= maxloc(η_ve) ; bc2!(η_vem)
                η_vevm[2:end-1,2:end-1] .= maxloc(η_vev); bc2!(η_vevm)
            else
                η_vem[2:end-1,2:end-1]  .= maxloc(η_vep) ; bc2!(η_vem)
                η_vevm[2:end-1,2:end-1] .= maxloc(η_vepv); bc2!(η_vevm)
            end
            @. dτ_ρ    = vpdτ*max_lxy/Re/η_vem
            @. dτ_ρv   = vpdτ*max_lxy/Re/η_vevm
            @. Gdτ     = vpdτ^2/dτ_ρ/(r+2.0)
            @. Gdτv    = vpdτ^2/dτ_ρv/(r+2.0)
            @. η_ve_τ  = 1.0/(1.0/η + 1.0/η_e + 1.0/Gdτ)
            @. η_ve_τv = 1.0/(1.0/ηv + 1.0/η_ev + 1.0/Gdτv)
            # pressure
            ∇V    .= diff(Vx, dims=1)./dx .+ diff(Vy, dims=2)./dy
            @. Pr -= r*Gdτ*∇V
            # strain rates
            Exx   .= diff(Vx, dims=1)./dx .- 1//3*∇V
            Eyy   .= diff(Vy, dims=2)./dy .- 1//3*∇V
            Exy[2:end-1,2:end-1] .= 0.5.*(diff(Vx[2:end-1,:], dims=2)./dy .+ diff(Vy[:,2:end-1], dims=1)./dx)
            Exyc  .= av(Exy)
            # viscoelastic strain rates
            @. Exx_e  = Exx  + τxx_o /2.0/η_e
            @. Eyy_e  = Eyy  + τyy_o /2.0/η_e
            @. Exy_e  = Exy  + τxy_o /2.0/η_ev
            @. Exyc_e = Exyc + τxyc_o/2.0/η_e
            # viscoelastic pseudo-transient strain rates
            @. Exx_τ  = Exx_e  + τxx /2.0/Gdτ
            @. Eyy_τ  = Eyy_e  + τyy /2.0/Gdτ
            @. Exy_τ  = Exy_e  + τxy /2.0/Gdτv
            @. Exyc_τ = Exyc_e + τxyc/2.0/Gdτ
            # stress update
            @. τxx    = 2.0*η_ve_τ *Exx_τ
            @. τyy    = 2.0*η_ve_τ *Eyy_τ
            @. τxy    = 2.0*η_ve_τv*Exy_τ
            @. τxyc   = 2.0*η_ve_τ *Exyc_τ
            # stress and strain rate invariants
            @. τii    = sqrt(0.5*(τxx^2 + τyy^2) + τxyc*τxyc)
            @. Eii    = sqrt(0.5*(Exx_τ^2 + Eyy_τ^2) + Exyc_τ*Exyc_τ)
            # yield function
            @. F      = τii - τ_y - Pr.*sinϕ
            @. λ      = max(F,0.0)/(η_ve_τ + η_reg)
            @. dQdTxx = 0.5*τxx /τii
            @. dQdTyy = 0.5*τyy /τii
            @. dQdTxy =     τxyc/τii
            # plastic correction
            @. τxx    = 2.0*η_ve_τ *(Exx_τ  - λ*dQdTxx)
            @. τyy    = 2.0*η_ve_τ *(Eyy_τ  - λ*dQdTyy)
            @. τxyc   = 2.0*η_ve_τ *(Exyc_τ - 0.5*λ*dQdTxy)
            τxy[2:end-1,2:end-1] .= 2.0 .* η_ve_τv[2:end-1,2:end-1].*(Exy_τ[2:end-1,2:end-1] .- 0.5 .* av(λ.*dQdTxy))
            @. τii    = sqrt(0.5*(τxx^2 + τyy^2) + τxyc*τxyc)
            @. Fchk   = τii - τ_y - Pr*sinϕ - λ*η_reg
            @. η_vep  = τii / 2.0 / Eii * 19.3 # nx, ny = 63, 63
            # @. η_vep  = τii / 2.0 / Eii * 19.3 * 1.99 # nx, ny = 127, 127
            η_vepv[2:end-1,2:end-1] .= av(η_vep); bc2!(η_vep)
            # velocity update
            dVx .= av_xa(dτ_ρ) .* (.-diff(Pr, dims=1)./dx .+ diff(τxx, dims=1)./dx .+ diff(τxy[2:end-1,:], dims=2)./dy)
            dVy .= av_ya(dτ_ρ) .* (.-diff(Pr, dims=2)./dy .+ diff(τyy, dims=2)./dy .+ diff(τxy[:,2:end-1], dims=1)./dx .+ av_ya(ρg))
            @. Vx[2:end-1,:] = Vx[2:end-1,:] + dVx
            @. Vy[:,2:end-1] = Vy[:,2:end-1] + dVy
            if iter % nchk == 0
                @. τxx_r   = 2.0*η_vec*Exx_e
                @. τyy_r   = 2.0*η_vec*Eyy_e
                @. τxy_r   = 2.0*η_vev*Exy_e
                @. τxyc_r  = 2.0*η_vec*Exyc_e
                @. τii     = sqrt(0.5*(τxx_r^2 + τyy_r^2) + τxyc_r*τxyc_r)
                # yield function
                @. F      = τii - τ_y - Pr.*sinϕ
                @. λ      = max(F,0.0)/(η_vec + η_reg)
                @. dQdTxx = 0.5*τxx_r /τii
                @. dQdTyy = 0.5*τyy_r /τii
                @. dQdTxy =     τxyc_r/τii
                # plastic correction
                @. τxx_r  = 2.0*η_vec*(Exx_e  - λ*dQdTxx)
                @. τyy_r  = 2.0*η_vec*(Eyy_e  - λ*dQdTyy)
                @. τxyc_r = 2.0*η_vec*(Exyc_e - 0.5*λ*dQdTxy)
                τxy_r[2:end-1,2:end-1] .= 2.0 .* η_vev[2:end-1,2:end-1].*(Exy_e[2:end-1,2:end-1] .- 0.5 .* av(λ.*dQdTxy))
                Rx .= .-diff(Pr, dims=1)./dx .+ diff(τxx_r, dims=1)./dx .+ diff(τxy_r[2:end-1,:], dims=2)./dy
                Ry .= .-diff(Pr, dims=2)./dy .+ diff(τyy_r, dims=2)./dy .+ diff(τxy_r[:,2:end-1], dims=1)./dx .+ av_ya(ρg)
                norm_Rx = norm(Rx)/sqrt(length(Rx)); norm_Ry = norm(Ry)/sqrt(length(Ry)); norm_∇V = norm(∇V)/sqrt(length(∇V))
                err = maximum([norm_Rx, norm_Ry, norm_∇V])
                @printf("it = %d, iter = %d, err = %1.2e norm[Rx=%1.2e, Ry=%1.2e, ∇V=%1.2e] (Fchk=%1.2e) \n", it, iter, err, norm_Rx, norm_Ry, norm_∇V, maximum(Fchk))
            end
            iter += 1; niter += 1
        end
        println(norm(Exyc_τ))
        t += dt; push!(evo_t, t); push!(evo_τxx, maximum(τxx))
        p1 = heatmap(xc,yc,τii',aspect_ratio=1,xlims=(-lx/2,lx/2),ylims=(-ly/2,ly/2),title="τii")
        # p3 = heatmap(xc,yc,η_vep',aspect_ratio=1,xlims=(-lx/2,lx/2),ylims=(-ly/2,ly/2),title="τii")
        p2 = plot(evo_t, evo_τxx , legend=false, xlabel="time", ylabel="max(τxx)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
             plot!(evo_t, 2.0.*ebg.*η0.*(1.0.-exp.(.-evo_t.*G0./η0)), linewidth=2.0) # analytical solution for VE loading
             plot!(evo_t, 2.0.*ebg.*η0.*ones(size(evo_t)), linewidth=2.0)            # viscous flow stress
        display(plot(p1,p2))
    end
    println(niter)
    return
end
# action
Stokes2D_vep()


text:
struct PSFDeltaFunction <: AbstractPSF
end
export PSFDeltaFunction


text:
module IEEE802_15_4

export  PacketSink,
        PacketSource,
        Modulation,
        BPSK,
        OQPSK,
        spread,
        CHIP_MAP_BPSK,
        make_packet,
        exec

abstract PackState
type SyncOnZero     <: PackState end
type SFDSearch      <: PackState end
type HeaderSearch   <: PackState end
type PayloadCollect <: PackState end

abstract Modulation
type BPSK  <: Modulation end
type OQPSK <: Modulation end

const MAX_PKT_LEN    = 127
const CHIP_MAP_BPSK  = Uint16[ 0b000100110101111, 0b111011001010000 ]
const CHIP_MASK_BPSK = 0b0011111111111110

type PacketSink{M}
    modulation::Type{M}        # BPSK, OQPSK, etc
    state::Type                # what is our PackState
    chips_per_symbol::Int      # how many chips per demodulated symbol. 15 For BPSK, 32 for OQPSK
    chip_map::Vector{Uint16}   # Chip to bit(s) mapping for modulation T
    sync_sequence::Uint8       # 802.15.4 standard is 4x 0 bytes and 1x0xA7, we will ignore the first byte
    chip_error_threshold::Int  # how many bits may be wrong in sync vector
    chip_shift_reg::Uint16     # chips are shifted in and decoded to look for first 0
    chip_shift_count::Int      # how many chips have we shifted into chip_shift_reg
    last_diff_enc_bit::Int     # previous not-yet-decoded bit, need to keep it to do differential decoding
    packet::Vector{Uint8}      # assembled payload
    packet_byte::Uint8         # byte being assembled
    packet_byte_count::Int     # which bit of d_packet_byte we're working on
    packetlen::Int             # length of packet
    packetlen_cnt::Int         # how many so far
    payload_cnt::Int           # how many bytes in payload
    input_idx::Int             # our location in the input vector
    packet_byte_bit_count::Int # how many bits have we shifted into packet_byte
    differential::Bool         # do differential decoding on the despread chips
    preable_zeros_count::Int   # how many zeros (bits, not byes) received before receiving the SFD byte
    verbosity::Int             # how much debug information to print
    chip_errors::Int           # how many chip errors for the whole received packet
end

function PacketSink( modType, chip_error_threshold; diff_enc = false, verbosity = 0 )

    modulation            = BPSK
    state                 = SyncOnZero
    chips_per_symbol      = 15
    chip_map              = CHIP_MAP_BPSK
    sync_sequence         = 0xA7
    chip_error_threshold  = chip_error_threshold
    chip_shift_reg        = zero( Uint16 )
    chip_shift_count      = 0
    last_diff_enc_bit     = 0
    packet                = zeros( Uint8, MAX_PKT_LEN )
    packet_byte           = zero( Uint8 )
    packet_byte_count     = 0
    packetlen             = 0
    packetlen_cnt         = 0
    payload_cnt           = 0
    input_idx             = 0
    packet_byte_bit_count = 0
    differential          = diff_enc
    preable_zeros_count   = 0
    chip_errors           = 0

    PacketSink( modulation,
                state,
                chips_per_symbol,
                chip_map,
                sync_sequence,
                chip_error_threshold,
                chip_shift_reg,
                chip_shift_count,
                last_diff_enc_bit,
                packet,
                packet_byte,
                packet_byte_count,
                packetlen,
                packetlen_cnt,
                payload_cnt,
                input_idx,
                packet_byte_bit_count,
                differential,
                preable_zeros_count,
                verbosity,
                chip_errors )
end

function count_bit_diffs( a, b, mask )
    count_ones( (a & mask) $ (b & mask) )
end

function chips_to_bit( sink::PacketSink{BPSK}, chips::Integer )
    best_match = 0xFF
    min_errors = 16
    is_valid_seq = false

    for symbol in 0:1
        reference_chips = CHIP_MAP_BPSK[ symbol+1 ]
        error_count     = count_bit_diffs( chips, reference_chips, CHIP_MASK_BPSK )

        if error_count < min_errors
            best_match = symbol
            min_errors = error_count
        end
    end

    is_valid_seq = min_errors <= sink.chip_error_threshold

    return ( is_valid_seq, best_match, min_errors )
end

function set_state( sink::PacketSink{BPSK}, ::Type{SyncOnZero} )
    sink.verbosity > 1 && println( sink.state, " -> SyncOnZero" )
    sink.state               = SyncOnZero
    sink.last_diff_enc_bit   = 0
    sink.chip_shift_reg      = zero( Uint16 )
    sink.chip_shift_count    = 0
    sink.packet_byte         = 0
    sink.preable_zeros_count = 0
end

function set_state( sink::PacketSink{BPSK}, ::Type{SFDSearch} )
    sink.verbosity > 1 && println( sink.state, " -> SFDSearch" )
    sink.state                 = SFDSearch
    sink.chip_shift_reg        = zero( Uint16 )
    sink.chip_shift_count      = 0
    sink.packet_byte           = 0
    sink.packet_byte_bit_count = 0
end

function set_state( sink::PacketSink{BPSK}, ::Type{HeaderSearch} )
    sink.verbosity > 1 && println( sink.state, " -> HeaderSearch" )
    sink.state                 = HeaderSearch
    sink.packet_byte_bit_count = 0
    sink.packet_byte_count     = 0
end

function set_state( sink::PacketSink{BPSK}, ::Type{PayloadCollect} )
    sink.verbosity > 1 && println( sink.state, " -> PayloadCollect" )
    sink.state                 = PayloadCollect
    sink.chip_errors           = 0
    sink.packet_byte           = 0
    sink.packet_byte_bit_count = 0
    sink.packet_byte_count     = 0
end

function synconzero( sink::PacketSink{BPSK}, input::Vector )
    while sink.input_idx <= length( input )
        sink.verbosity > 2 && @printf( "SyncOnZero. input_idx: %d, chip_shift_reg: %s, chip_shift_count: %d, packet_byte_bit_count: %d, packet_byte: %s, %s\n", sink.input_idx, bin(sink.chip_shift_reg, 15),sink.chip_shift_count, sink.packet_byte_bit_count, hex(sink.packet_byte,2), bin(sink.packet_byte,8))
        sink.chip_shift_reg    = uint16( (sink.chip_shift_reg >> 1) | ((input[sink.input_idx] & 1)<<14) )
        sink.input_idx        += 1
        sink.chip_shift_count += 1

        (is_valid_seq, rx_bit, chip_errors) = chips_to_bit( sink, sink.chip_shift_reg )

        if sink.differential
            diff_dec_bit           = rx_bit $ sink.last_diff_enc_bit
            sink.last_diff_enc_bit = rx_bit
            rx_bit                 = diff_dec_bit
        end

        if is_valid_seq && rx_bit == 0
            sink.preable_zeros_count += 1
            set_state( sink, SFDSearch )
            break
        end

        sink.last_diff_enc_bit = 0
    end

end

function sfdsearch( sink::PacketSink{BPSK}, input::Vector )
    while sink.input_idx <= length( input )
        sink.verbosity > 2 && @printf( "SFDSearch. input_idx: %d, chip_shift_reg: %s, chip_shift_count: %d, packet_byte_bit_count: %d, packet_byte: %s, %s\n", sink.input_idx, bin(sink.chip_shift_reg, 15),sink.chip_shift_count, sink.packet_byte_bit_count, hex(sink.packet_byte,2), bin(sink.packet_byte,8))
        sink.chip_shift_reg    = uint16( (sink.chip_shift_reg >> 1) | ((input[sink.input_idx] & 1)<<14) )
        sink.input_idx        += 1
        sink.chip_shift_count += 1

        if sink.chip_shift_count == 15
            (is_valid_seq, rx_bit, chip_errors) = chips_to_bit( sink, sink.chip_shift_reg )
            if !is_valid_seq
                set_state( sink, SyncOnZero )
                break
            end

            if sink.differential
                diff_dec_bit           = rx_bit $ sink.last_diff_enc_bit
                sink.last_diff_enc_bit = rx_bit
                rx_bit                 = diff_dec_bit
            end

            if rx_bit == 0
                sink.preable_zeros_count += 1
            end

            if sink.preable_zeros_count > 4*8+3
                sink.verbosity > 1 && @printf( "Received %d zeros in the preable, going back to SyncOnZero\n", sink.preable_zeros_count )
                set_state( sink, SyncOnZero )
                break
            end

            sink.packet_byte            = uint8( (sink.packet_byte >> 1) | (rx_bit << 7) )
            sink.chip_shift_count       = 0
            sink.chip_shift_reg         = zero( sink.chip_shift_reg )

            if sink.packet_byte == sink.sync_sequence
                sink.verbosity > 0 && println( "Got start of frame delimiter (SFD)")
                set_state( sink, HeaderSearch )
                break
            end
        end
    end

end

function headersearch( sink::PacketSink{BPSK}, input::Vector )
    while sink.input_idx <= length( input )
        sink.verbosity > 2 && @printf( "HeaderSearch. input_idx: %d, chip_shift_reg: %s, chip_shift_count: %d, packet_byte_bit_count: %d, packet_byte: %s, %s\n", sink.input_idx, bin(sink.chip_shift_reg, 15),sink.chip_shift_count, sink.packet_byte_bit_count, hex(sink.packet_byte,2), bin(sink.packet_byte,8))
        sink.chip_shift_reg    = uint16( (sink.chip_shift_reg >> 1) | ((input[sink.input_idx] & 1)<<14) )
        sink.input_idx        += 1
        sink.chip_shift_count += 1

        if sink.chip_shift_count == 15
            (is_valid_seq, rx_bit, chip_errors) = chips_to_bit( sink, sink.chip_shift_reg )
            if !is_valid_seq
                set_state( sink, SyncOnZero )
                break
            end

            if sink.differential
                diff_dec_bit           = rx_bit $ sink.last_diff_enc_bit
                sink.last_diff_enc_bit = rx_bit
                rx_bit                 = diff_dec_bit
            end

            sink.packet_byte            = uint8( (sink.packet_byte >> 1) | (rx_bit << 7) )
            sink.packet_byte_bit_count += 1
            sink.chip_shift_count       = 0
            sink.chip_shift_reg         = zero( sink.chip_shift_reg )

            if sink.packet_byte_bit_count == 8
                sink.packetlen = sink.packet_byte
                sink.verbosity > 0 && println( "Packet lenght = ", sink.packetlen )

                sink.packet = zeros( Uint8, sink.packetlen )

                if  sink.packetlen > MAX_PKT_LEN || sink.packetlen < 1
                    set_state( sink, SyncOnZero )
                    break
                end

                set_state( sink, PayloadCollect )
                break
            end
        end
    end
end

function payloadcollect( sink::PacketSink{BPSK}, input::Vector )
    while sink.input_idx <= length( input )
        sink.verbosity > 2 && @printf( "PayloadCollect. input_idx: %d, chip_shift_reg: %s, chip_shift_count: %d, packet_byte_bit_count: %d, packet_byte: %s, %s\n", sink.input_idx, bin(sink.chip_shift_reg, 15),sink.chip_shift_count, sink.packet_byte_bit_count, hex(sink.packet_byte,2), bin(sink.packet_byte,8))
        sink.chip_shift_reg    = uint16( (sink.chip_shift_reg >> 1) | ((input[sink.input_idx] & 1)<<14) )
        sink.input_idx        += 1
        sink.chip_shift_count += 1

        if sink.chip_shift_count == 15
            (is_valid_seq, rx_bit, chip_errors) = chips_to_bit( sink, sink.chip_shift_reg )
            if !is_valid_seq
                set_state( sink, SyncOnZero )
                break
            end

            if sink.differential
                diff_dec_bit           = rx_bit $ sink.last_diff_enc_bit
                sink.last_diff_enc_bit = rx_bit
                rx_bit                 = diff_dec_bit
            end

            sink.packet_byte            = uint8( (sink.packet_byte >> 1) | (rx_bit << 7) )
            sink.packet_byte_bit_count += 1
            sink.chip_shift_count       = 0
            sink.chip_shift_reg         = zero( sink.chip_shift_reg )
            sink.chip_errors           += chip_errors

            if sink.packet_byte_bit_count == 8
                sink.packet_byte_count += 1
                sink.verbosity > 1 && @printf( "packet[%d] = 0x%s ", sink.packet_byte_count, hex(sink.packet_byte) )
                sink.packet[sink.packet_byte_count] = sink.packet_byte
                sink.packet_byte_bit_count          = 0

                if sink.packet_byte_count + 1 > sink.packetlen
                    returnpacket( sink )
                    set_state( sink, SyncOnZero )
                    break
                end
            end
        end
    end
end

function returnpacket( sink::PacketSink{BPSK} )
    total_chips = sink.packetlen * sink.chips_per_symbol * 8
    lqi = 255-int(sink.chip_errors/total_chips*255)
    println()
    @printf( "LQI: %d\n", int(lqi) )
    print("┏")
    for i in 1:min( sink.packetlen*5+5, 8*5+5 )
        print( "━" )
    end
    println("┓")
    print("┃ 00: ")
    for i in 1:sink.packetlen
        print( "0x", hex(sink.packet[i], 2), " " )
        if mod( i, 8 ) == 0
            println("┃")
            print( "┃ ", dec(i,2), ": " )
        end
    end
    println()
    print("┗")
    for i in 1:min( sink.packetlen*5+5, 8*5+5 )
        print( "━" )
    end
    println("┛")
    println()
    set_state( sink, SyncOnZero )
end

function exec( sink::PacketSink{BPSK}, input::Vector )
    sink.input_idx = 1

    while sink.input_idx <= length(input)

        if sink.state == SyncOnZero
            synconzero( sink, input )
        elseif sink.state == SFDSearch
            sfdsearch( sink, input )
        elseif sink.state == HeaderSearch
            headersearch( sink, input )
        elseif sink.state == PayloadCollect
            payloadcollect( sink, input )
        end

    end
end

type PacketSource
    modulation        # BPSK, OQPSK, etc
    chips_per_symbol  # how many chips per demodulated symbol. 15 For BPSK, 32 for OQPSK
    chip_map          # Chip to bit(s) mapping for modulation T
end

function PacketSource( ::Type{BPSK} )
    PacketSource( BPSK, 15, CHIP_MAP_BPSK )
end

function spread( chip_map::AbstractVector, chips_per_symbol::Integer, packet::AbstractVector{Uint8}; diff_enc = false )

    packet_len        = length( packet )
    frame_len         = packet_len * 8 * chips_per_symbol
    frame             = zeros( Uint8, frame_len )
    last_diff_enc_bit = 0
    frame_idx         = 0

    for packet_idx in 0:packet_len-1
        packet_byte = packet[packet_idx+1]

        for bit_idx in 0:7
            packet_bit = (packet_byte >> bit_idx) & 0x01

            if diff_enc
                packet_bit        = packet_bit $ last_diff_enc_bit
                last_diff_enc_bit = packet_bit
            end

            chips = chip_map[packet_bit+1]

            for chip_idx in 0:chips_per_symbol-1
                frame[frame_idx+1] = (chips >> chip_idx) & 0x01
                frame_idx       += 1
            end

        end
    end

    return frame
end

function make_packet( source::PacketSource, input::Vector{Uint8} )
    payload_len = uint8( length( input ) & 0b0111111 )
    packet      = [ 0x00, 0x00, 0x00, 0x00, 0xA7, payload_len, input ]
end

function exec( source::PacketSource, input::Vector{Uint8}; diff_enc = false )
    packet = make_packet( source, input )
    spread( CHIP_MAP_BPSK, source.chips_per_symbol, packet, diff_enc = diff_enc )
end

end # module


text:
    set_in = SetProg.Sets.PolynomialSublevelSetAtOrigin(d, p[source(s, edge)])
    set_out = SetProg.Sets.PolynomialSublevelSetAtOrigin(d, p[target(s, edge)])
    A = dynamicfort(s, edge, γ)
    @constraint(model, A * set_in ⊆ set_out)


text:
using Lora

DATADIR = "../../data"
SUBDATADIR = "amsmmala"
OUTDIR = "../output"

npars = 4

nchains = 10
nmcmc = 110000
nburnin = 10000

ratio = Array(Float64, nchains)
essizes = Array(Float64, nchains, npars)
results = Dict{Symbol, Any}()

for i in 1:nchains
  ratio[i] = acceptance(readdlm(joinpath(DATADIR, SUBDATADIR, "diagnostics"*lpad(string(i), 2, 0)*".csv"), ',', Bool))
  essizes[i, :] = ess(readdlm(joinpath(DATADIR, SUBDATADIR, "chain"*lpad(string(i), 2, 0)*".csv"), ',', Float64), 2)

  println("Iteration ", i, " of ", nchains, " completed")
end

results[:rate] = mean(ratio)
results[:ess] = mean(essizes, 1)
results[:time] = mean(readdlm(joinpath(DATADIR, SUBDATADIR, "times.csv"), ',', Float64))
results[:efficiency] = minimum(results[:ess])/results[:time]

writedlm(
  joinpath(OUTDIR, "logit_amsmmala_summary.csv"),
  hcat(
    results[:rate],
    results[:ess],
    results[:time],
    results[:efficiency]
  ),
  ','
)

writedlm(
  joinpath(OUTDIR, "logit_amsmmala_summary.txt"),
  Any[
    round(results[:rate], 2)
    [Int64(i) for i in round(results[:ess])]
    round(results[:time], 2)
    round(results[:efficiency], 2)
  ]',
  " & "
)


text:
vars=5
i_train=10_000
i_test=10_000
cc1=10.0
neurons=Int64(floor(i_train/50)) # this automatically
# generates number of neurons== number of clusters (§3.1)
_case_noise_=5

ff(x)=-x[1]+(1/2)x[2].^2-(1/3)x[3].^3+(1/4)x[4].^4-(1/5)x[5].^5
xx_train=0.8*rand(i_train,vars).+0.1
yy_train=zeros(i_train); for i=1:i_train yy_train[i]=ff(xx_train[i,:]) end
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=0.8;yy_train.+=0.1

if _case_noise_==1
    __noise__=rand(i_train)
elseif _case_noise_==2
    __noise__=rand(Normal(0.0, 1.0),i_train)
elseif _case_noise_==3
    __noise__=rand(GeneralizedPareto(0.0,1.0,0.25),i_train)
elseif _case_noise_==4
    __noise__=rand(LogNormal(0.25, 0.8),i_train)
elseif _case_noise_==5
    __noise__=[rand(LogNormal(0.25, 0.8),i_train);10.0.-rand(Exponential(3.0),i_train)]
    __noise__=[__noise__;15.0.+rand(Frechet(10.0,10.0),i_train)]
    __noise__=__noise__[1:3:end]
end
__noise__.-=mean(__noise__)
__noise__./=mean(abs.(__noise__))/(0.05*0.45)
histogram(__noise__,label=sum(__noise__)/sum(abs.(__noise__)),
        title=mean(abs.(__noise__)))
iso=sortperm(yy_train)
plot(yy_train[iso])
yy_train.+=__noise__
display(plot!(yy_train[iso],legend=false))

xx_test=0.8*rand(i_test,vars).+0.1
yy_test=zeros(i_test); for i=1:i_test yy_test[i]=ff(xx_test[i,:]) end
# yy_test.+=noise_perc.*rand(i_test)
yy_test.-=mi;yy_test./=ma;yy_test.*=0.8;yy_test.+=0.1

scatter3d(xx_train[:,1],xx_train[:,2],yy_train)
display(scatter3d!(xx_test[:,1],xx_test[:,2],yy_test))

using DecisionTree, Statistics
num_trees=100;
regr_3 = RandomForestRegressor(n_trees=num_trees);
tmp_xx_train=copy(xx_train);tmp_yy_train=copy(yy_train);
tmp_xx_test=copy(xx_test);tmp_yy_test=copy(yy_test);
DecisionTree.fit!(regr_3, tmp_xx_train, tmp_yy_train[:,1]);
predl1 = DecisionTree.predict(regr_3, tmp_xx_train);
predl1_test = DecisionTree.predict(regr_3, tmp_xx_test);
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot(err[ier],label="Random Forests",color=:black,
    linestyle=:dash,legend=:topleft,linewidth=3))
__err__=string(mean(abs.(err)))

sci_model = fit!(AdaBoostRegressor(n_estimators=100),xx_train,yy_train)
predl1=ScikitLearn.predict(sci_model,xx_train)
predl1_test=ScikitLearn.predict(sci_model,xx_test)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot!(err[ier],label="AdaBoost",color=:black,
    linestyle=:dashdot,legend=:topleft,linewidth=3))
__err__=string(__err__," ",mean(abs.(err)))

num_round = 10
bst = xgboost(copy(xx_train), num_round, label = copy(yy_train), eta = 1, max_depth = 7)
predl1 = XGBoost.predict(bst, copy(xx_train))
predl1_test = XGBoost.predict(bst, copy(xx_test))
mean(abs.(predl1-yy_train))
mae=mean(abs.(predl1_test-yy_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot!(err[ier],label="XGBoost",color=:black,
    linestyle=:dot,legend=:topleft,linewidth=3))
__err__=string(__err__," ",mean(abs.(err)))

##### ANNBN RBF

include(string(path1,"/src/ANNBN.jl"))
inds_all,n_per_part=ANNBN.___clustering(neurons,xx_train,200)

a_all,a_layer1,layer1=ANNBN.train_layer_1_rbf(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train,cc1)
predl1,layer1_train=ANNBN.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,n_per_part,inds_all,xx_train,vars,cc1)
predl1_test,layer1_test=ANNBN.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
maete=mean(abs.(yy_test-predl1_test))
err=yy_test-predl1_test;ier=sortperm(err)
display(plot!(err[ier],label="ANNBN",color=:black,linestyle=:dashdotdot,
    legend=:topleft,linewidth=3))
__err__=string(__err__," ",mean(abs.(err)))

println(__err__)

#


text:
function handle_step(obj::OrekitFixedStepHandler, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "handleStep", void, (SpacecraftState, jboolean), arg0, arg1)
end

function init(obj::OrekitFixedStepHandler, arg0::SpacecraftState, arg1::AbsoluteDate, arg2::jdouble)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate, jdouble), arg0, arg1, arg2)
end


text:
using DiffEqProblemLibrary, Sundials, Base.Test

prob = prob_ode_linear
dt = 1//2^(4)
saveat = float(collect(0:dt:1))
sol = solve(prob,CVODE_BDF())
sol = solve(prob,CVODE_Adams())
@test sol.errors[:l2] < 1e-3
sol = solve(prob,CVODE_Adams(),reltol=1e-5)
@test sol.errors[:l2] < 1e-5
sol = solve(prob,CVODE_Adams(),saveat=saveat)

@test sol.t == saveat

sol = solve(prob,CVODE_Adams(),saveat=dt)

@test sol.t == saveat

sol = solve(prob,CVODE_Adams(),saveat=saveat,save_everystep=true)

@test sol.t != saveat
@test intersect(sol.t,saveat) == saveat

sol = solve(prob,CVODE_Adams(),saveat=saveat,save_everystep=true,save_start=false)

@test sol.t[1] != 0

sol = solve(prob,CVODE_Adams(),tstops=[0.2,0.5,0.7])
@test all(t ∈ sol.t for t in [0.2,0.5,0.7])

prob = prob_ode_2Dlinear
sol = solve(prob,CVODE_BDF())
sol = solve(prob,CVODE_Adams())
sol = solve(prob,CVODE_Adams(),saveat=saveat)

@test sol.t == saveat

sol = solve(prob,CVODE_Adams(),saveat=saveat,save_everystep=false)

@test sol.t == saveat

sol = solve(prob,CVODE_Adams(),tstops=[0.9])

@test 0.9 ∈ sol.t

# Test the other function conversions
k = (t,u,du) -> du[1] = u[1]
prob = ODEProblem(k,[1.0],(0.0,1.0))
sol = solve(prob,CVODE_BDF())
h = (t,u) -> u
u0 = [1.0 2.0
      3.0 2.0]
prob = ODEProblem(h,u0,(0.0,1.0))
sol = solve(prob,CVODE_BDF())

# Test Algorithm Choices
sol1 = solve(prob,CVODE_BDF(method=:Functional))
sol2 = solve(prob,CVODE_BDF(linear_solver=:Band,jac_upper=3,jac_lower=3))
sol3 = solve(prob,CVODE_BDF(linear_solver=:Diagonal))
sol4 = solve(prob,CVODE_BDF(linear_solver=:GMRES))
sol5 = solve(prob,CVODE_BDF(linear_solver=:FGMRES))
sol6 = solve(prob,CVODE_BDF(linear_solver=:PCG))
sol7 = solve(prob,CVODE_BDF(linear_solver=:BCG))
sol8 = solve(prob,CVODE_BDF(linear_solver=:TFQMR))

@test isapprox(sol1[end],sol2[end],rtol=1e-3)
@test isapprox(sol1[end],sol3[end],rtol=1e-3)
@test isapprox(sol1[end],sol4[end],rtol=1e-3)
@test isapprox(sol1[end],sol5[end],rtol=1e-3)
@test isapprox(sol1[end],sol6[end],rtol=1e-3)
@test isapprox(sol1[end],sol7[end],rtol=1e-3)
@test isapprox(sol1[end],sol8[end],rtol=1e-3)
# Backwards
prob = deepcopy(prob_ode_2Dlinear)
prob2 = ODEProblem(prob.f,prob.u0,(1.0,0.0))
sol = solve(prob2,CVODE_BDF())
@test maximum(diff(sol.t)) < 0 # Make sure all go negative