text:
# ## model
include("trajopt_model_v2.jl")
include(joinpath(@__DIR__, "..", "..", "..", "src/dynamics/centroidal_quadruped/visuals.jl"))

vis = Visualizer()
open(vis)

# ## horizon
T = 11
h = 0.1

# ## centroidal_quadruped
s = get_simulation("centroidal_quadruped", "flat_3D_lc", "flat")
model = s.model
env = s.env
nx = 2 * model.nq
nc = 4 # model.nc
nu = model.nu + nc + 4 * nc + nc + 4 * nc + 1
nθ = 53

# ## model
d1 = DTO.Dynamics((y, x, u, w) -> centroidal_quadruped_dyn1(
	model, env, [h], y, x, u, w), nx + nθ + nx, nx, nu)
dt = DTO.Dynamics((y, x, u, w) -> centroidal_quadruped_dynt(
	model, env, [h], y, x, u, w), nx + nθ + nx, nx + nθ + nx, nu)

dyn = [d1, [dt for t = 2:T-1]...]

# ## initial conditions
body_height = 0.3
foot_x = 0.17
foot_y = 0.17

function nominal_configuration(model::CentroidalQuadruped)
    [
        0.0; 0.0; body_height;
        0.0; 0.0; 0.0;
        foot_x; foot_y; 0.0;
        foot_x;-foot_y; 0.0;
       -foot_x; foot_y; 0.0;
       -foot_x;-foot_y; 0.0;
    ]
end

q1 = nominal_configuration(model)
qM = nominal_configuration(model)
qT = nominal_configuration(model)
q_ref = nominal_configuration(model)

x1 = [q1; q1]
xM = [qM; qM]
xT = [qT; qT]
x_ref = [q_ref; q_ref]

# ## objective
function obj1(x, u, w)
	J = 0.0
	J += 0.5 * transpose(x[1:nx] - x_ref) * Diagonal(ones(nx)) * (x[1:nx] - x_ref)
	J += 0.5 * transpose(u) * Diagonal([ones(model.nu); zeros(nu - model.nu)]) * u
    J += 1000.0 * u[end] # slack
    J += 0.5 * transpose(u[model.nu + 4 .+ (1:20)]) * Diagonal(1.0 * ones(20)) * u[model.nu + 4 .+ (1:20)]

	return J
end

function objt(x, u, w)
	J = 0.0

    u_prev = x[nx .+ (1:53)]
    w = (u - u_prev) ./ h
    J += 0.5 * 10.0 * dot(w[1:end-1], w[1:end-1])

	J += 0.5 * transpose(x[1:nx] - x_ref) * Diagonal(ones(nx)) * (x[1:nx] - x_ref)
	J += 0.5 * transpose(u) * Diagonal([ones(model.nu); zeros(nu - model.nu)]) * u
    J += 1000.0 * u[end] # slack
    J += 0.5 * transpose(u[model.nu + 4 .+ (1:20)]) * Diagonal(1.0 * ones(20)) * u[model.nu + 4 .+ (1:20)]

	return J
end

function objT(x, u, w)
	J = 0.0
	J += 0.5 * transpose(x[1:nx] - x_ref) * Diagonal(ones(nx)) * (x[1:nx] - x_ref)
    return J
end

c1 = DTO.Cost(obj1, nx, nu)
ct = DTO.Cost(objt, nx + nθ + nx, nu)
cT = DTO.Cost(objT, nx + nθ + nx, 0)
obj = [c1, [ct for t = 2:T-1]..., cT];

# ## constraints
ql = q1
qu = q1

# initial condition
xl1 = [q1; q1]
xu1 = [q1; q1]
xlt = [-Inf * ones(nx); -Inf * ones(nθ); -Inf * ones(nx)]
xut = [Inf * ones(nx); Inf * ones(nθ); Inf * ones(nx)]

# final condition
xlT = [q1; q1; -Inf * ones(nθ); -Inf * ones(nx)]
xuT = [q1; q1; Inf * ones(nθ); Inf * ones(nx)]

ul = [-Inf * ones(model.nu); zeros(nu - model.nu)]
uu = [Inf * ones(model.nu); Inf * ones(nu - model.nu)]

bnd1 = DTO.Bound(nx, nu, state_lower=xl1, state_upper=xu1, action_lower=ul, action_upper=uu)
bndt = DTO.Bound(nx + nθ + nx, nu, state_lower=xlt, state_upper=xut, action_lower=ul, action_upper=uu)
bndT = DTO.Bound(nx + nθ + nx, 0, state_lower=xlT, state_upper=xuT)
bnds = [bnd1, [bndt for t = 2:T-1]..., bndT];

function constraints_1(x, u, w)
    [
     # equality (16)
     contact_constraints_equality(model, env, h, x, u, w);
     # inequality (28)
     contact_constraints_inequality_1(model, env, h, x, u, w);
     x[6 .+ (1:12)] - q1[6 .+ (1:12)];
     x[18 + 6 .+ (1:12)] - q1[6 .+ (1:12)];
    ]
end

function constraints_t(x, u, w)
    [
     # equality (16)
     contact_constraints_equality(model, env, h, x, u, w);
     # inequality (32)
     contact_constraints_inequality_t(model, env, h, x, u, w);
     x[6 .+ (1:12)] - q1[6 .+ (1:12)];
     x[18 + 6 .+ (1:12)] - q1[6 .+ (1:12)];
    ]
end

function constraints_T(x, u, w)
    [
     # inequality (8)
     contact_constraints_inequality_T(model, env, h, x, u, w);
     x[6 .+ (1:12)] - q1[6 .+ (1:12)];
     x[18 + 6 .+ (1:12)] - q1[6 .+ (1:12)];
    ]
end

con1 = DTO.Constraint(constraints_1, nx, nu, idx_ineq=collect(16 .+ (1:28)))
cont = DTO.Constraint(constraints_t, nx + nθ + nx, nu, idx_ineq=collect(16 .+ (1:32)))
conT = DTO.Constraint(constraints_T, nx + nθ + nx, nu, idx_ineq=collect(0 .+ (1:8)))
cons = [con1, [cont for t = 2:T-1]..., conT];

# ## problem
direct_solver = DTO.solver(dyn, obj, cons, bnds,
    options=DTO.Options(
        tol=1.0e-3,
        constr_viol_tol=1.0e-3,
        ))

# ## initialize
x_interpolation = [x1, [[x1; zeros(nθ); zeros(nx)] for t = 2:T]...]
u_guess = [1.0e-4 * rand(nu) for t = 1:T-1] # may need to run more than once to get good trajectory
DTO.initialize_states!(direct_solver, x_interpolation)
DTO.initialize_controls!(direct_solver, u_guess)

# ## solve
@time DTO.solve!(direct_solver)

# ## solution
x_sol, u_sol = DTO.get_trajectory(direct_solver)
@show x_sol[1]
@show x_sol[T]
maximum([u[nu] for u in u_sol[1:end-1]])

# ## visualize
# q_sol = state_to_configuration([x[1:nx] for x in x_sol])
visualize!(vis, model, x_sol, Δt=h)

q_opt = [x_sol[1][1:model.nq], [x[model.nq .+ (1:model.nq)] for x in x_sol]...]
v_opt = [(x[model.nq .+ (1:model.nq)] - x[0 .+ (1:model.nq)]) ./ h for x in x_sol]
u_opt = [u[1:model.nu] for u in u_sol]
λ_opt = [u[model.nu .+ (1:4)] for u in u_sol]
b_opt = [u[model.nu + 4 .+ (1:16)] for u in u_sol]

q_opt = [x_sol[1][1:model.nq], [x[model.nq .+ (1:model.nq)] for x in x_sol]...]
v_opt = [(x[model.nq .+ (1:model.nq)] - x[0 .+ (1:model.nq)]) ./ h for x in x_sol]
u_opt = [u[1:model.nu] for u in u_sol]
γ_opt = [u[model.nu .+ (1:4)] for u in u_sol]
b_opt = [u[model.nu + 4 .+ (1:16)] for u in u_sol]
ψ_opt = [u[model.nu + 4 + 16 .+ (1:4)] for u in u_sol]
η_opt = [u[model.nu + 4 + 16 + 4 .+ (1:16)] for u in u_sol]

qm = q_opt
vm = v_opt
um = u_opt
γm = γ_opt
bm = b_opt
ψm = ψ_opt
ηm = η_opt

μm = model.μ_world
hm = h
timesteps = range(0.0, stop=(h * (length(qm) - 2)), length=(length(qm) - 2))
plot(timesteps, hcat(qm[2:end-1]...)', labels="")
plot(timesteps, hcat(um...)', labels="")
plot(timesteps, hcat(γm...)', labels="")
plot(timesteps, hcat(bm...)', labels="")
plot(timesteps, hcat(ψm...)', labels="")
plot(timesteps, hcat(ηm...)', labels="")

using JLD2
@save joinpath(@__DIR__, "stand_10Hz.jld2") qm um γm bm ψm ηm μm hm
@load joinpath(@__DIR__, "stand_10Hz.jld2") qm um γm bm ψm ηm μm hm


text:
###### State and States ######
mutable struct CarSt
    s::Float64
    v::Float64
    r::UInt16
end
CarSt() = CarSt(0.0, 0.0, UInt16(1))
Base.:(==)(a::CarSt, b::CarSt) = a.s == b.s && a.v == b.v && a.r == b.r
Base.hash(St::CarSt, h::UInt) = hash(St.s, hash(St.v, hash(St.r, h)))

struct Sts
    Ego::CarSt
    Other::CarSt
end
Sts() = Sts(CarSt(), CarSt())
Base.:(==)(A::Sts, B::Sts) = A.Ego == B.Ego && A.Other == B.Other
Base.hash(sts::Sts, h::UInt) = hash(sts.Ego, hash(sts.Other, h))

###### Observation ######
mutable struct CarOb
    v::Float64
    x::Float64
    y::Float64
end
CarOb() = CarOb(0.0, 0.0, 0.0)
Base.:(==)(a::CarOb, b::CarOb) = a.v == b.v && a.x == b.x && a.y == b.y
Base.hash(ob::CarOb, h::UInt) = hash(ob.v, hash(ob.x, hash(ob.y, h)))

mutable struct Route
    Geos::Vector{Geometry}
    intersect_Infos::Dict{UInt16, Vector{Tuple{String, Point2D{Float64}, Float64}}} # UInt16: the index of the car whoes route intersects with that of ego car;
    Vref::Vector{Float64}
    Aref::Vector{Float64}
    Length::Float64

    function Route(Geos::Vector{Geometry})  # constructor
        intersect_Infos = Dict{UInt16, Vector{Tuple{String, Point2D{Float64}, Float64}}}()
        Length = round(sum(map(x->x.Length, Geos)), 1.0e-8)
        return new(Geos, intersect_Infos, [], [], Length)
    end
end

# POMDP{State, Action, Observation}; action is the acceleration of ego car.
struct DrivePOMDP <: POMDP{Sts, Symbol, CarOb}
        ##### constants #####
        discount_factor::Float64
        Δt::Float64 # time step
        Δs::Float64
        Δv::Float64
        ##### ranges #####
        sV::Vector{Float64}
        vV::Vector{Float64}
        rV::Vector{UInt16}
        Aset::NamedTuple{(:min, :max, :comfort),Tuple{Float64, Float64, Float64}}

        Routes::Vector{Route}
        SSpace::Vector{Sts}
        OSpace::Vector{CarOb}
        SsInit::Sts

        # used only in IDM model
        Stopline::Float64
        Smin::Float64
end

function round(x::Float64, pre::Union{Float64, Int64}) # pre: precision
    n = floor(Int, x/pre)
    m = x - n*pre
    if m >= pre/2
        return trunc((n+1)*pre; digits= 13)
    else
        return trunc(n*pre; digits= 13)
    end
end

ceil(x::Float64, pre::Union{Float64, Int64}) = sign(x)*pre*ceil(Int, abs(x)/pre)

function Rmat(α::Union{Float64, Irrational})
    return [cos(α) -sin(α); sin(α) cos(α)]
end

ext: jl


text:
include("../required.jl")

outcome = "bart"

## Prepare data

# Load presence-absence data for all species
@load joinpath("data", "jld2", "$(outcome)-distributions.jld2") distributions
# Create matrix Y
Y = Ymatrix(distributions; transform=false)
# Create matrix of observed sites only
Yobs = _Yobs(Y)
# Get richness
richness = calculate_richness(Y, distributions[1])
# Get LCBD values
lcbd_rel = calculate_lcbd(Y, distributions[1]; transform=true, relative=true)
# Get non-relative values
lcbd_abs = calculate_lcbd(Y, distributions[1]; transform=true, relative=false)
round.(Float64.(extrema(lcbd_abs)); sigdigits=4)
lcbd = lcbd_abs

# Create a simpler MWE
mwe = rand(10, 10) / 2e5

## Option 1: Default plots
gr()
heatmap(mwe; c=:viridis, colorbar_title="LCBD value")
heatmap(lcbd; c=:viridis, colorbar_title="LCBD value")
plot_layer(lcbd; c=:viridis, colorbar_title="LCBD value")

#= Problems:
    - Colorbar title overlaps with colorbar ticks
   Pros:
    - That's what I have already
=#

## Option 2: Add scale in colorbar title
# Suggestion from https://discourse.julialang.org/t/plots-jl-clean-scientific-formatting/38112/2
gr()
heatmap(mwe * 1e5; c=:viridis, colorbar_title="LCBD value (x 10⁻⁵)")
layer1e5 = replace(similar(lcbd), 0.0 => 1e5)
heatmap(lcbd * layer1e5; c=:viridis, colorbar_title="LCBD value (x 10⁻⁵)")
plot_layer(lcbd * layer1e5; c=:viridis, colorbar_title="LCBD value (x 10⁻⁵)")

# Alternative to rescale
lcbd_scaled = lcbd_abs .* 1e5
extrema(lcbd_scaled)
lcbd_expo = rescale(lcbd_abs, (extrema(lcbd_scaled)))
collect(lcbd_expo)
isapprox(collect(lcbd_expo), lcbd_scaled)

#= Problems:
    - Still a bit of overlap, but could be ok
    - Not sure if it's considered a valid presentation form
   Pros:
    - Easier to fit with what I already have
=#

## Option 2: PyPlot
pyplot()
heatmap(mwe; c=:viridis, colorbar_title="LCBD value")
heatmap(lcbd; c=:viridis, colorbar_title="LCBD value")
heatmap(replace(lcbd.grid, nothing => NaN); c=:viridis, colorbar_title="LCBD value")
plot_layer(lcbd; c=:viridis, colorbar_title="LCBD value")

# Customize colorbar ticks
# New features from https://github.com/JuliaPlots/Plots.jl/pull/3346
# But only for PyPlot for now...
# GR extension discussed in https://github.com/JuliaPlots/Plots.jl/issues/3174
heatmap(rand(10, 10); colorbar_ticks=([0.3, 0.6], ["0.300", "0.60"]))
# How to check these possible attributes
plotattr(:Subplot)
plotattr("colorbar")

#= Problems:
    - super slow for layers (long to display)
    - looks terrible for layers
    - some weird lines on the heatmap, like a bad printer...
    - doesn't work with plot_layer
    - python
   Pros:
    - has scientific notation out-of-the-box
    - can customize colorbar ticks
=#

## Option 3: PlotlyJS (or Plotly)
plotlyjs()
heatmap(mwe; c=:viridis, colorbar_title="LCBD value")
heatmap(lcbd; c=:viridis, colorbar_title="LCBD value")
heatmap(lcbd; c=:viridis, colorbar_title="LCBD value", leftmargin=100px)
plot_layer(lcbd; c=:viridis, colorbar_title="LCBD value")
savefig("test-plotlyjs.png")
savefig("test-plotlyjs.pdf")

#= Problems:
    - micro sign, not sure how to remove
    - space between plot & colorbar (because of aspect_ratio)
        - can improve with leftmargin
    - terrible png quality, but pdf is OK and not too big
   Pros:
    - overall fast & still nice
    - works with plot_layer
=#

## Option 4: PlotlyJS.jl (not the same as the plotlyjs backend)
# From https://discourse.julialang.org/t/scientific-notation-in-plots-jl-colorbar/57793/9
using PlotlyJS: PlotlyJS
trace = PlotlyJS.heatmap(;
    z=replace(lcbd.grid, nothing => NaN)',
    colorscale="Viridis",
    colorbar_title="LCBD value",
    colorbar_exponentformat="power",
)
PlotlyJS.plot(trace)

#= Problems:
    - doesn't work with plot_layer
    - works differently from Plots, so I don't know how to tweak it
    - not sure if it can be used to combine plots either, which I need for subplots
   Pros:
    - scientific notation
=#

## Option 5: CairoMakie
using CairoMakie: CairoMakie
CairoMakie.activate!()
CairoMakie.heatmap(
    replace(lcbd.grid, nothing => NaN)'; c=:viridis, colorbar_title="LCBD value"
)

#= Problems:
    - don't really know how it works
    - no scientific notation
   Pros:
    - scientific notation
=#

## Option 6: Play with layouts
# https://discourse.julialang.org/t/single-colorbar-for-heatmap-grid-layout/50609/3
# https://discourse.julialang.org/t/plots-jl-shared-colorbar-with-subplots/47269/4
gr()
begin
    p1 = plot(lcbd; c=:viridis)
    p2 = plot(; frame=:none)
    annotate!(p2, 0.5, 0.5, text("LCBD value", 11, :center, 90.0))
    l = @layout [a b{0.01w}]
    plot(p1, p2; layout=l, leftmargin=[0.0mm -7.0mm])
end
begin
    p1 = plot_layer(lcbd; c=:viridis)
    p2 = plot(; frame=:none)
    annotate!(p2, 0.5, 0.5, text("LCBD value", 11, :center, 90.0))
    l = @layout [a b{0.01w}]
    plot(p1, p2; layout=l, leftmargin=[0.0mm -7.0mm])
end
# Comparison with relative values
begin
    p1 = plot_layer(lcbd_rel; c=:viridis)
    p2 = plot(; frame=:none)
    annotate!(p2, 0.5, 0.5, text("LCBD value", 11, :center, 90.0))
    l = @layout [a b{0.01w}]
    plot(p1, p2; layout=l, leftmargin=[0.0mm -20.0mm])
end
plot_layer(lcbd_rel; c=:viridis, colorbar_title="LCBD value")
# Just need to play with margins

#= Problems:
    - cannot format colorbar to similar units
   Pros:
    - Should work with everything I have if I simply define a new plot_layer function
=#

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
"Loads csv with timeslots formatted as in 'test/Sample data/timeslottest_output.csv' "
function loadtimeslots(path)
    colparser = Dict(:startTime => dateformat"H:M:S",:endTime => dateformat"H:M:S")
    res = loadtable(path,colparsers=colparser)|> @map((resourceID = _.resourceID, type= _.type,dayID = _.dayID,timeslotID = _.timeslotID,startTime = Time(_.startTime),endTime= Time(_.endTime),booked = _.booked == "true"  )) |> collect
    return table(res,pkey = [:resourceID,:dayID,:timeslotID])
end

"Saves any object to a CSV file"
function save(object,path::String)
    CSVFiles.save(path, object)
end

function loadresources(path)
    loadtable(path,indexcols= [:intID])
end
"Loads csv with patientdata formatted as in 'test/Sample data/patienttest_output.csv' "
function loadpatients(path)
        loadtable(path,indexcols= [:intID])
end
"Loads csv with \"time delta\" which is the number of days that need to be in between two types of visits. Data should be formatted as in 'test/Sample data/GUCHamb_timeDelta.csv' "
function loadTimeDelta(path_TimeDelta::String)
    loadndsparse(path_TimeDelta,indexcols=1:2)
end

"Loads csv with visitdata formatted as in 'test/Sample data/visittest_output.csv' "
function loadvisits(path)
    loadtable(path,indexcols= [:intID])
end

"Loads csv with a solution formatted as in 'test/Sample data/solutiontest_input.csv' "
function loadsolution(path)
    loadtable(path;indexcols=[:visitID])
end


text:
function bad()
    println(div(1, 0))
    false
end
function good()
    println("good")
    false
end

println(false && bad())
println(true && good())
println(good() && bad())

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# ---
# title: Generate your cover on the fly
# cover: assets/logo.svg
# description: this demo shows you how to generate cover on the fly
# date: 2020-09-13
# author: Jun Tian
# ---

# There're many reasons that you don't want to mannually manage the cover image.
# DemoCards.jl allows you to generate the card cover on the fly for demos written
# in julia.
#
# Let's do this with a simple example

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:BasicDQN},
    ::Val{:CartPole},
    ::Nothing;
    seed = 123,
    save_dir = nothing,
)
    rng = StableRNG(seed)
    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))

    policy = Agent(
        policy = QBasedPolicy(
            learner = BasicDQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = Chain(
                        Dense(ns, 128, relu; init = glorot_uniform(rng)),
                        Dense(128, 128, relu; init = glorot_uniform(rng)),
                        Dense(128, na; init = glorot_uniform(rng)),
                    ) |> cpu,
                    optimizer = ADAM(),
                ),
                batch_size = 32,
                min_replay_history = 100,
                loss_func = huber_loss,
                rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = 500,
                rng = rng,
            ),
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
        ),
    )
    stop_condition = StopAfterStep(10_000)
    hook = ComposedHook(TotalRewardPerEpisode(), TimePerStep())
    Experiment(policy, env, stop_condition, hook, "")
end

# Now let's get start!

#+ tangle=false
run(RL.Experiment(Val(:JuliaRL), Val(:BasicDQN), Val(:CartPole), nothing))

1 + 2


text:
@testset "different distance transform implementations" begin

    box = [20   -20;           # RA
          -10     10] * degree  # DEC
    shape, wcs = geometry(Pixell.WCS.WCSTransform, box, (1/2) * degree)

    for kk in 1:300
        m = Enmap(ones(shape), wcs)
        # m[2,2] = 0.0
        for i in 1:30
            m[rand(2:(size(m,1)-1)), rand(2:(size(m,2)-1))] = 0.0
        end

        distmap = distance_transform(ExactSeqSDT(), m)
        distmap_bf = distance_transform(BruteForceSDT(), m)
        distmap_approx = distance_transform(ApproxSeqSDT(), m)

        @test sum(abs.((distmap_bf .- distmap))) < 1e-13
        @test (sum(distmap_bf .!= distmap_approx) / prod(size(distmap))) < 0.2
    end

end

##
@testset "distance transform metric" begin
    box = [20   -20;           # RA
           0     10] * degree  # DEC
    shape, wcs = geometry(Pixell.WCS.WCSTransform, box, (1/2) * degree)

    m = Enmap(ones(shape), wcs)
    m[1,1] = 0.0
    distmap = distance_transform(ExactSeqSDT(), m)

    αs = pix2sky(m, collect(1:size(m,1)), ones(size(m,1)))[1]
    δs = pix2sky(m, ones(size(m,2)), collect(1:size(m,2)))[2]

    for i in axes(m,1)
        @test (αs[1] - αs[i]) ≈ distmap[i,1]
    end
    for j in axes(m,2)
        @test (δs[j] - δs[1]) ≈ distmap[1,j]
    end
end


text:
using FlightGNC
using UnPack, ComponentArrays, Transducers
using LinearAlgebra, DifferentialEquations
using Plots

function main(x₀::Float64, h₀::Float64, γ₀::Float64, V₀::Float64, Δt_plan_update::Float64, s_guidance::IASCG)
    @unpack p, env_model = s_guidance
    @unpack m, S, C_D_0, K, ρ₀, H, g = env_model

    # Environments
    env     = VerticalPlaneDynamicMissile(m, S, C_D_0, K, ρ₀, H, g)

    # Initial condition
    s₀      = State(env)(x₀, h₀, γ₀, V₀)
    R₀      = 1.0
    c₀      = zeros(p)
    params₀ = Params(env)(R₀, c₀)

    # Simulation parameters
    Δt_save        = 0.01
    t_sim_f        = 100.0

    # callbacks
    function guidance_update(s_guidance::IASCG)
        return function (integrator)
            # println("t: $(integrator.t)")
            # println(integrator.p)
            if abs(integrator.u[1] - s_guidance.x_f) > 1E3
                integrator.p = IASCG_Planner(s_guidance, integrator.u, integrator.p)
            end
        end
    end
    cb_plan_update = PeriodicCallback(guidance_update(s_guidance), Δt_plan_update; initial_affect = true)

    function condition_stop(s_guidance::IASCG)
        @unpack x_f = s_guidance
        return function (u, t, integrator)
            u.x - x_f
        end
    end
    affect!(integrator) = terminate!(integrator)
    cb_stop = ContinuousCallback(condition_stop(s_guidance), affect!)

    cb = CallbackSet(cb_plan_update, cb_stop)

    # Execute Simulation
    # prob: DE problem, df: DataFrame
    simulator = Simulator(s₀,  # initial condition
                        apply_inputs(Dynamics!(env);
                                    u = IASCG_cmd(s_guidance)),
                        params₀;
                        Problem = :ODE,
                        solver = Tsit5(),
                        tf = t_sim_f)

    # Non-interactive simulation: solve approach (automatically reinitialised)
    @time df = solve(simulator;
                    savestep = Δt_save,
                    callback = cb,
                    reltol = 1e-8,
                    abstol = 1e-8)

    # # Execute Simulation
    # # prob: DE problem, df: DataFrame
    # @time prob, df = FSimBase.sim(
    #         s₀,  # initial condition
    #         apply_inputs(Dynamics!(env);
    #                      u = IASCG_cmd(s_guidance)),
    #         params₀;  # dynamics!; apply_inputs is exported from FS and is so useful for systems with inputs
    #         tf = t_sim_f,
    #         savestep = Δt_save,  # savestep is NOT simulation step
    #         solver = Tsit5(),
    #         callback = cb,
    #         reltol = 1e-8,
    #         abstol = 1e-8
    #        )

    t_sim   = df.time
    x_sim   = df.sol |> Map(datum -> datum.x)  |> collect
    h_sim   = df.sol |> Map(datum -> datum.h)  |> collect
    γ_sim   = rad2deg.( df.sol |> Map(datum -> datum.γ)  |> collect )
    V_sim   = df.sol |> Map(datum -> datum.V)  |> collect
    A_L_sim = df.sol |> Map(datum -> datum.u)  |> collect
    R_sim   = df.sol |> Map(datum -> datum.R)  |> collect
    c_sim   = df.sol |> Map(datum -> datum.c)  |> collect

    c_sim  = hcat(c_sim...)'

    return df, t_sim, x_sim, h_sim, γ_sim, V_sim, A_L_sim, R_sim, c_sim
end

function plan_path(x₀::Float64, h₀::Float64, γ₀::Float64, V₀::Float64, s_guidance::IASCG)
    @unpack p, Δx_pred, env_model = s_guidance
    @unpack g = env_model

    s₀     = State(env_model)(x₀, h₀, γ₀, V₀)
    R₀     = 1.0
    c₀     = zeros(p)
    param₀ = Params(env_model)(R₀, c₀)
    param_plan = IASCG_Planner(s_guidance, s₀, param₀)

    x_plan = x₀:Δx_pred:x_f
    order = p-1:-1:0
    path_s = x_plan |> Map(x -> begin
                            ϕ  = (x / x_ref).^order
                            ϕ′ = order .* [ϕ[2:end]; 0] / x_ref
                            ϕ″ = order .* (order .- 1) .* [ϕ[3:end]; zeros(2)] / x_ref^2
                            # path_deriv = h_ref * [ϕ'; ϕ′'; ϕ″'] * param_plan[:c]
                            path = [h_ref * dot(ϕ, param_plan[:c]) atand(h_ref * dot(ϕ′, param_plan[:c])) h_ref * dot(ϕ″, param_plan[:c])]
                        end) |> tcollect
    path_s  = vcat(path_s...)
    h_plan  = path_s[:,1]
    γ_plan  = path_s[:,2]
    h″_plan = path_s[:,3]

    prob    = ODEProblem(Speed_Dynamics(s_guidance, param_plan[:c]), V₀, (x₀, x_f))
    sol     = DifferentialEquations.solve(prob, Tsit5(), reltol = 1e-6, abstol = 1e-6, saveat = x_plan)
    V_plan  = sol.u

    A_L_plan = ( (h″_plan .* ((cos.(deg2rad.(γ_plan))).^2) .* (V_plan.^2)) .+ g ) .* cos.(deg2rad.(γ_plan))

    return x_plan, h_plan, γ_plan, V_plan, A_L_plan, param_plan[:R], param_plan[:c]
end

# -----------------------------------------------------------------------------------------
# Simulation Input Default Values
# -----------------------------------------------------------------------------------------
# Algorithm parameters
p       = 10
x_ref   = 40E3
h_ref   = 3E3
h_des   = 1E3
Approx_mode = 1
n       = 1.5
Δx_pred = 1
x_f     = x_ref
h_f_d   = 0E3
γ_f_d   = deg2rad(-30)
V_f_d   = 400.0

# environmental parameters
m       = 544
S       = 0.258
C_D_0   = 0.126
K       = 0.370
ρ₀      = 1.225
H       = 8435
g       = 9.805

env_model  = VerticalPlaneDynamicMissile(m, S, C_D_0, K, ρ₀, H, g)

# initial conditions
R   = 1.0
x₀  = 0.0
h₀  = 3E3
γ₀  = deg2rad(5)
V₀  = 1000.0

s_guidance = IASCG(p, x_ref, h_ref, h_des, Approx_mode, n, x₀, V₀, Δx_pred, x_f, h_f_d, γ_f_d, V_f_d, env_model)


text:
export  ProPolyHermite, Cpro, degree, ProPolyH, prohermite_coeffmatrix,
        FamilyProPolyHermite, FamilyScaledProPolyHermite,
        derivative,
        evaluate!, evaluate,
        vander!, vander

"""
ProPolyHermite <: Hermite

An immutable structure for probabilistic Hermite polynomials

## Fields
-  `m` : order of the polynomial
-  `P` : probabilistic Hermite polynomial of order m
- `scaled` : with rescaling to have unitary norm
"""
struct ProPolyHermite <: Hermite
    m::Int64
    P::ImmutablePolynomial{Float64}
    scaled::Bool
end

function Base.show(io::IO, P::ProPolyHermite)
    println(io, string(P.m)*"-th order probabilistic Hermite polynomial "*string(P.P)*", scaled = "*string(P.scaled))
end

# Hen(x)  = (-1)ⁿ*exp(x²/2)dⁿ/dxⁿ exp(-x²/2)
# Hen′(x) = n*Hen-1(x)
# Hen″(x) = n*(n-1)*Hen-1(x)

"""
$(TYPEDSIGNATURES)

Normalization constant of the m-th order probabilistic Hermite polynomial
"""
# Cpro(m::Int64) =sqrt(sqrt(2*π) * gamma(m+1))
Cpro(m::Int64) =sqrt(gamma(m+1)) #sqrt(2π) is already accounted for by sampling

"""
$(TYPEDSIGNATURES)

Normalization constant of the m-th order probabilistic Hermite polynomial
"""
Cpro(P::ProPolyHermite) = Cpro(P.m)

"""
$(TYPEDSIGNATURES)

Return the degree of the probabilistic Hermite polynomial `P`.
"""
degree(P::ProPolyHermite) = P.m

# Adapted https://people.sc.fsu.edu/~jburkardt/m_src/hermite_polynomial/h_polynomial_coefficients.m

#    Input, integer N, the highest order polynomial to compute.
#    Note that polynomials 0 through N will be computed.
#
#    Output, real C(1:N+1,1:N+1), the coefficients of the Hermite
#    polynomials.

# Recurrence relation:
# Hen+1 = x*Hen - n*Hen-1

"""
$(TYPEDSIGNATURES)

Computes the coefficients of the probabilistic Hermite polynomials up to order m.
"""
function prohermite_coeffmatrix(m::Int64)

    if m < 0
        return Float64[]
    end

    coeff = zeros(m+1, m+1)
    coeff[1, 1] = 1.0

    if m == 0
        return coeff
    end

    coeff[2, 2] = 1.0

    if m==1
        return coeff
    end

    for i  = 2:m
        coeff[i+1, 1]      =  -      (i - 1) * coeff[i-1, 1]
        coeff[i+1, 2:i-1]  =                   coeff[i  , 1:i-2] -
                                     (i - 1) * coeff[i-1, 2:i-1]
        coeff[i+1, i]      =                   coeff[i  , i-1]
        coeff[i+1, i+1]    =                   coeff[i  , i]
    end
    return coeff
end

const ProPolyH = prohermite_coeffmatrix(CstMaxDegree)

function ProPolyHermite(m::Int64;scaled::Bool= false)
    @assert m>=0 "The order of the polynomial should be >=0"
    if scaled ==false
        return ProPolyHermite(m, ImmutablePolynomial(view(ProPolyH,m+1,1:m+1)), scaled)
    else
        C = 1/Cpro(m)
        return ProPolyHermite(m, ImmutablePolynomial(C*view(ProPolyH,m+1,1:m+1)), scaled)
    end
end

"""
$(TYPEDSIGNATURES)

Evaluates the probabilistic Hermite polynomial `P`.
"""
(P::ProPolyHermite)(x) = P.P(x)

const FamilyProPolyHermite = ntuple(i->ProPolyHermite(i-1), CstMaxDegree+1)
const FamilyScaledProPolyHermite = ntuple(i->ProPolyHermite(i-1; scaled = true), CstMaxDegree+1)

"""
$(TYPEDSIGNATURES)

Computes the k-th derivative of the probabilistic Hermite polynomial `P`.
"""
function derivative(P::ProPolyHermite, k::Int64)
    m = P.m
    @assert k>=0 "This function doesn't compute anti-derivatives of Hermite polynomials"
    if m>=k
        factor = exp(loggamma(m+1) - loggamma(m+1-k))
        if P.scaled == false
            return ImmutablePolynomial(factor*ProPolyH[m-k+1,1:m-k+1])
        else
            C = 1/Cpro(m)
            return ImmutablePolynomial(factor*C*ProPolyH[m-k+1,1:m-k+1])
        end
    else
        return ImmutablePolynomial((0.0))
    end
end

"""
$(TYPEDSIGNATURES)

Evaluates in-place the probabilistic Hermite polynomial `P`.
"""
function evaluate!(dV, P::ProPolyHermite, x)
    m = P.m
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    # H_0(x) = 1, H_1(x) = x
    col = view(dV,:,1)
    @avx @. col = 1.0

    if P.scaled
        rmul!(col, 1/Cpro(0))
    end
    if m == 0
        return dV
    end

    col = view(dV,:,2)
    @avx @. col = x

    if P.scaled
        rmul!(col, 1/Cpro(1))
    end

    if m == 1
        return dV
    end

    if P.scaled
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = Cpro(i-1)* x * col - Cpro(i-2)*(i-1)*colm1
            # dV[:,i+1] = 2.0*Cphy(i-1)*x .* dV[:,i] - 2.0*Cphy(i-2)*i*dV[:,i-1]
            rmul!(colp1, 1/Cpro(i))
        end
    else
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = x * col - (i-1)*colm1
        end
    end
    return dV
end

"""
$(TYPEDSIGNATURES)

Evaluates the probabilistic Hermite polynomial `P`.
"""
evaluate(P::ProPolyHermite, x::Array{Float64,1}) = evaluate!(zeros(size(x,1), P.m+1), P, x)

# He_{n}^(k)(x) =  n!/(n-k)! H_{n-k}(x)

"""
$(TYPEDSIGNATURES)

Computes in-place the Vandermonde matrix of the k-th derivative of the probabilistic Hermite polynomial `P`.
"""
function vander!(dV, P::ProPolyHermite, k::Int64, x)
    m = P.m

    if k==0
        evaluate!(dV, P, x)
    # Derivative
    elseif k==1
        # Use recurrence relation H′_{n+1} = (1 + 1/n)*(x*H′_{n} - n*H′_{n-1})
        N = size(x,1)
        @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

        col = view(dV,:,1)
        @avx @. col = 0.0
        if m==0
            return dV
        end

        col = view(dV,:,2)
        @avx @. col = 1.0

        if m==1
            return dV
        end

        col = view(dV,:,3)
        @avx @. col = 2.0*x

        if m==2
            return dV
        end

        @inbounds for i=3:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)

            @avx @. colp1 = (1.0 + 1.0/(i-1.0))*(x * col - (i-1.0) * colm1)
        end

        if P.scaled
            @inbounds for i=1:m
                colp1 = view(dV,:,i+1)
                factor = exp(loggamma(i+1) - loggamma(i+1-k))
                rmul!(colp1, 1/(Cpro(i)*sqrt(factor)))
            end

        end
        return dV
    # Second Derivative
    elseif k==2
        # Use recurrence relation He″_{n+1} = (n+1)/(n-1)*(x*He″_{n} - 2*n*He″_{n-1})
        N = size(x,1)
        @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

        col = view(dV,:,1)
        @avx @. col = 0.0
        if m==0
            return dV
        end

        col = view(dV,:,2)
        @avx @. col = 0.0

        if m==1
            return dV
        end

        col = view(dV,:,3)
        @avx @. col = 2.0

        if m==2
            return dV
        end

        col = view(dV,:,4)
        @avx @. col = 6.0 * x

        if m==3
            return dV
        end

        @inbounds for i=4:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)

            @avx @. colp1 = ((i)/(i-2.0))*(x * col - (i-1.0) * colm1)
        end

        if P.scaled
            @inbounds for i=1:m
                colp1 = view(dV,:,i+1)
                factor = exp(loggamma(i+1.0) - loggamma(i+1.0-k))
                @avx rmul!(colp1, 1/(Cpro(i)*sqrt(factor)))
            end

        end
        return dV
    end
end

"""
$(TYPEDSIGNATURES)

Computes the Vandermonde matrix of the k-th derivative of the probabilistic Hermite polynomial `P`.
"""
vander(P::ProPolyHermite, k::Int64, x::Array{Float64,1}) = vander!(zeros(size(x,1), P.m+1), P, k, x)


text:
function getPlotArr(arr)
		intervalArray=arr   #as we refactor code from above code,don't change internal code structure, use temp variable
     arrlength=length(intervalArray)
     xsArr=[0:i:2π for i in intervalArray]

    plotArr=[]
      for i in  1:arrlength
       ys=[sin(x) for x in xsArr[i]]
       p=scatter(xsArr[i],ys,color=:green,markersize=2,title="interval=$(intervalArray[i])",label=false,framestyle=:origin)
      push!(plotArr,p)
      end
    return  plotArr
end

max_stars_repo_path: Mathematical Chart /introduction.jl

max_stars_repo_name: phpsmarter/plutonote

max_issues_repo_path: Mathematical Chart /introduction.jl

max_issues_repo_name: phpsmarter/plutonote

max_forks_repo_path: Mathematical Chart /introduction.jl

max_forks_repo_name: phpsmarter/plutonote


text:
module FrankieMae

import Pandas
import Dates

# ReadFiles
fmopen(parentpath,startrange,endrange) = begin
## for daterange, create an array where 1:4 is repeated endrange-(startrange+1) and with startrange quarter:4 appended to front and 1:endrange quarter is appended to the end
 daterange =

[Pandas.read_csv(
 join(parentpath, daterange[n]),header=-1,delimiter="|") for n=1:length(daterange)]
end # module

max_stars_repo_name: alejandrorbraun/FrankieMae

max_issues_repo_name: alejandrorbraun/FrankieMae

max_forks_repo_name: alejandrorbraun/FrankieMae


text:
function attach_dirichlet(cellmatvec,cellvals,cellmask=Fill(true,length(cellvals)))
  k = AttachDirichletMap()
  lazy_map(k,cellmatvec,cellvals,cellmask)
end

struct AttachDirichletMap <: Map
  muladd::MulAddMap{Int}
  AttachDirichletMap() = new(MulAddMap(-1,1))
end

function Arrays.return_cache(k::AttachDirichletMap,matvec::Tuple,vals,mask)
  mat, vec = matvec
  return_cache(k.muladd,mat,vals,vec)
end

@inline function Arrays.evaluate!(cache,k::AttachDirichletMap,matvec::Tuple,vals,mask)
  if mask
    mat, vec = matvec
    vec_with_bcs = evaluate!(cache,k.muladd,mat,vals,vec)
    (mat, vec_with_bcs)
  else
    matvec
  end
end

function Arrays.return_cache(k::AttachDirichletMap,mat::AbstractMatrix,vals,mask)
  cm = return_cache(*,mat,vals)
  cv = CachedArray(mat*vals)
  fill!(cv.array,zero(eltype(cv)))
  (cm,cv)
end

@inline function Arrays.evaluate!(cache,k::AttachDirichletMap,mat::AbstractMatrix,vals,mask)
  cm, cv = cache
  if mask
    vec_with_bcs = evaluate!(cm,*,mat,vals)
    scale_entries!(vec_with_bcs,-1)
    (mat, vec_with_bcs)
  else
    if size(mat,1) != size(cv,1)
      m = axes(mat,1)
      setaxes!(cv,(m,))
      fill!(cv.array,zero(eltype(cv)))
    end
    (mat, cv.array)
  end
end


text:
@deprecate push_ctrl_gate! push_gate!
@deprecate pushfirst_ctrl_gate! pushfirst_gate!


text:
#=
This contains all of the pieces needed to add a new code gen target.
The following three functions must be provided.

1. get_external_language_elements() - file extensions, comment chars etc.
2. generate_external_code_layer(var, entities, terms, lorr, vors) - Turns symbolic expressions into code
3. generate_external_files(var, lhs_vol, lhs_surf, rhs_vol, rhs_surf) - Writes all files based on generated code

You will have access to anything in the Femshop module scope because this file will be included there.
=#

function get_external_language_elements()
    file_extension = ".jl";
    comment_char = "#";
    block_comment = ["#="; "=#"];
    # I imagine there could be more here in the future. Maybe not.

    return (file_extension, comment_char, block_comment);
end

#=
Translate the symbolic layer into code for the elemental assembly.
Input:
    var - unknown variable or array of variables.
    entities - leaf nodes of the computational graph. Type is SymEntity(see symexpression.jl) or number.
    terms - array of additive terms, separated for convenience. [t1, t2, t3] for expression t1+t2+t3 (see note below)
    lorr - LHS or RHS will have constant values LHS or RHS which are "lhs", "rhs"
    vors - volume or surface will have values "volume" or "surface"

Note on terms: Each term is a SymExpression(see symexpression.jl) which is essentially a Julia Expr and list of entities.
The Expr holds the computational graph for each term. LHS and RHS are determined by separating terms that involve
unknowns from those that don't. Also, in well formed FEM problems each term should be multiplied by a test function,
and thus quadrature can be done for each term separately.
=#
function generate_external_code_layer(var, entities, terms, lorr, vors)
    # The whole body of code is stored in a string that will be inserted in the
    # appropriate function during file writing.
    code = "";

    # ... code generation here ...

    return code;
end

# Create all of the code files.
# The input are code strings generated above. They will be inserted into their appropriate files.
function generate_external_files(var, lhs_vol, lhs_surf, rhs_vol, rhs_surf)
    # This can be split up into smaller functions as you wish.
    # generate_custom_utils_file();
    # generate_custom_main_file();
    # generate_custom_config_file();
    # etc.
end


text:
using DataStructures

module ImportSandbox

    export include_test_files!

    function include_test_files!(test_files, add_load_path=nothing)
        if !(add_load_path === nothing)
            push!(LOAD_PATH, add_load_path)
        end

        eval(quote
            using Jute
            $((:(include($test_file)) for test_file in test_files)...)
        end)
    end

end

using .ImportSandbox

function find_test_files(dir, test_file_postfix)
    fnames = String[]
    for (root, dirs, files) in walkdir(dir)
        for file in files
            if endswith(file, test_file_postfix)
                push!(fnames, joinpath(root, file))
            end
        end
    end
    fnames
end

function get_runtests_dir()
    runtests_path = abspath(PROGRAM_FILE)
    runtests_dir, _ = splitdir(runtests_path)
    runtests_dir
end

struct TestcaseInfo
    path :: Array{String, 1}
    name :: String
    tags :: Set{Symbol}
end

path_pair(tcinfo::TestcaseInfo) = tcinfo.path, tcinfo.name

path_string(tcinfo::TestcaseInfo) = join([tcinfo.path; tcinfo.name], "/")

function tag_string(tcinfo::TestcaseInfo, labels; full::Bool=false)
    fixtures_tag = isempty(labels) ? "" : join(labels, ",")
    tc_tag = full ? path_string(tcinfo) : tcinfo.name
    if length(labels) > 0
        tc_tag * " [" * fixtures_tag * "]"
    else
        tc_tag
    end
end

function _get_testcases(obj_dict, parent_path=String[])
    testcases = []

    for obj in obj_dict
        if isa(obj, TestGroup)
            group_testcases = _get_testcases(get_testcases(obj), [parent_path; obj.name])
            append!(testcases, group_testcases)
        elseif isa(obj, Testcase)
            tcinfo = TestcaseInfo(parent_path, obj.name, obj.tags)
            push!(testcases, tcinfo => obj)
        end
    end

    testcases
end

function get_testcases(obj_dict)
    _get_testcases(obj_dict)
end


text:
@testset "Flatten Function Tests" begin

  @test isdefined(Julz, :flatten) == true

end


text:
using Distributions
using Random

"""
    pnorm(x)

Julia equivalent of R's `pnorm`
"""
function pnorm(x)
    cdf(Normal(), x)
end

"""
    dnorm(x)

Julia equivalent of R's `dnorm`
"""
function dnorm(x)
    pdf(Normal(), x)
end

# """
#     cond_mean_on_LR(L::T, R::T, a1::T, a2::T, a3::T, μ1::T, μ2::T, μ3::T) where T <: Real

# Conditional mean on (L, R), defined as μ_{L,R} = (∫_L^R f(x) dx)^(-1) / (∫_L^R x f(x) dx)
# """
# function cond_mean_on_LR(L::T, R::T, a1::T, a2::T, a3::T, μ1::T, μ2::T, μ3::T) where T <: Real
#     pL = a1 * pnorm(L-μ1) + a2 * pnorm(L-μ2) + a3 * pnorm(L-μ3)
#     pR = a1 * pnorm(R-μ1) + a2 * pnorm(R-μ2) + a3 * pnorm(R-μ3)

#     ∫_L = -(a1 * dnorm(L-μ1) + a2 * dnorm(L-μ2) + a3 * dnorm(L-μ3)) +
#           a1 * μ1 * pnorm(L-μ1) + a2 * μ2 * pnorm(L-μ2) + a3 * μ3 * pnorm(L-μ3)
#     ∫_R = -(a1 * dnorm(R-μ1) + a2 * dnorm(R-μ2) + a3 * dnorm(R-μ3)) +
#           a1 * μ1 * pnorm(R-μ1) + a2 * μ2 * pnorm(R-μ2) + a3 * μ3 * pnorm(R-μ3)

#     μ_LR = (∫_R - ∫_L)/(pR - pL)

#     (pL, pR, ∫_L, ∫_R, μ_LR)
# end

"""
    cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)

Conditional mean on (L, R), defined as μ_{L,R} = (∫_L^R f(x) dx)^(-1) / (∫_L^R x f(x) dx)
"""
function cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)
    # (∫_l^r f(x) dx) = pR - pL
    pL = a1 * pnorm(L.-μ1) + a2 * pnorm(L.-μ2) #+ a3 * pnorm(L.-μ3)
    pR = a1 * pnorm(R.-μ1) + a2 * pnorm(R.-μ2) #+ a3 * pnorm(R.-μ3)

    # (∫_l^r x f(x) dx) = ∫_R - ∫_L
    ∫_L = -(a1 * dnorm(L.-μ1) + a2 * dnorm(L.-μ2)) + # + a3 * dnorm(L.-μ3)) +
          a1 * μ1 * pnorm(L.-μ1) + a2 * μ2 * pnorm(L.-μ2) #+ a3 * μ3 * pnorm(L.-μ3)
    ∫_R = -(a1 * dnorm(R-μ1) + a2 * dnorm(R-μ2)) + # + a3 * dnorm(R-μ3)) +
          a1 * μ1 * pnorm(R-μ1) + a2 * μ2 * pnorm(R-μ2) #+ a3 * μ3 * pnorm(R-μ3)
    ## adding a3/μ3 terms yields inaccurate result.

    # population conditional mean on (L, R)
    μ_LR = (∫_R .- ∫_L) ./ (pR .- pL)

    (pL, pR, ∫_L, ∫_R, μ_LR)
end

"""
    find_split(d, a1, a2, μ1, μ2, L, R;
               a3=0., μ3=0., find_deltas=false, find_split=!find_deltas)

Split the cluster (L, R) into subclusters that maximizes G_{L,R}.
Returns centered conditional mean (deltas) of the left and right subcluster
if `find_deltas` is true. Otherwise return the split point.
"""
function find_split(d, a1, a2, μ1, μ2, L, R;
                    a3=0., μ3=0., find_deltas=false, find_split=!find_deltas)
    # p, int, int_en evaluated at `d`
    ∫_d = -(a1 * dnorm(d.-μ1) + a2 * dnorm(d.-μ2)) + # + a3 * dnorm(d.-μ3)) +
          a1 * μ1 * pnorm(d.-μ1) + a2 * μ2 * pnorm(d.-μ2) #+ a3 * μ3 * pnorm(d.-μ3)
    p_d = a1 * pnorm(d.-μ1) + a2 * pnorm(d.-μ2) #+ a3 * pnorm(d.-μ3)

    pL, pR, ∫_L, ∫_R, μ_LR = cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)

    # find split point with normlized G'_{L,R} (=h)(supp. p.32)
    h = μ_LR .* (p_d .- pL).^2 + (pL + pR .- 2p_d) .* (∫_d .- ∫_L) - d .* (p_d .- pL) .* (pR .- p_d)
    split = d[2:end][h[2:end] .* h[1:end-1] .< 0]

    if find_deltas && (sum(μ1 .< split .< μ2) == 0)
        return(1., 1.)
    end

    # split point
    s = max(split[μ1 .< split .< μ2]...)

    if find_split
        return(s)
    else
        _, _, _, _, δ_1 = cond_mean_on_LR(L, s, a1, a2, a3, μ1, μ2, μ3)
        δ_1 -= (L + s)/2  # centering

        _, _, _, _, δ_2 = cond_mean_on_LR(s, R, a1, a2, a3, μ1, μ2, μ3)
        δ_2 -= (s + R)/2  # centering

        return(δ_1, δ_2)
    end

end

"""
    find_truncation(a1, a2, μ1, μ2,
                    R_ub=(μ2 + 8), L_lb=(μ1 - 8),
                    tuning_par1=0.01, tuning_par2=0.005,
                    a3=0, μ3=0, σ=1)

Find population right-truncation points for bimodal Normal mixtures.
Inputs:
    a1, a2, a3: weights for Normal distributions. a1 + a2 + a3 should be == 1
    μ1, μ2, μ3: means for normals μ1 = -μ2

Possible returns:
    0: negligible result (not a symmetric problem)
    -1: no truncation
"""
function find_truncation(a2, μ2;
                         a1=1-a2, μ1=-μ2, R_ub=μ2+8., L_lb=μ1-8.,
                         tuning_par1=0.01, tuning_par2=0.005,
                         a3=0., μ3=0., σ=1.)

    r = Array(R_ub:(-tuning_par1):μ1)

    # if the bimodal distn is symmetric, and a1 == a2,
    # then truncate both ends equally
    if a1 == a2 && μ2 == -μ1
        for i=1:length(r)

            # truncate both left and right with the same amount
            R = r[i]
            L = -R

            d = Array(L:0.01:R)
            δ_1, δ_2 = find_split(d, a1, a2, μ1, μ2, L, R; find_deltas=true)

            # L^* is the smallest L that satisfies δ_1 < 0 < δ_2. (supp. p.30)
            # Thus returns right-trunctation point
            if δ_1 < 0 < δ_2
                return(R)
            end
        end
    end

    # possible cases:
    #   a1 != a2 && μ2 != -μ1
    #   a1 != a2 && μ2 == -μ1
    if a1 != a2
        L_old = -1000
        for i=1:length(r)
            R = r[i]
            L = Array{Float64,1}(L_old:tuning_par2:(R-0.05))

            pL, pR, ∫_L, ∫_R, μ_LR = cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)

            # locate R_L using the equation μ_{L,R_L} - (L + R_L)/2.  (supp. p.30)
            temp = sign.(μ_LR .- (L .+ R)./2)

            if sum(temp[2:end] .* temp[1:end-1] .< 0) ≠ 0
                # this is L^* from (supp. p.28)
                L = L[2:end][temp[2:end] .* temp[1:end-1] .< 0][1]

                d = Array(L:0.01:R)
                δ_1, δ_2 = find_split(d, a1, a2, μ1, μ2, L, R; find_deltas=true)

                pL = a1 * pnorm(L-μ1) + a2 * pnorm(L-μ2) + a3 * pnorm(L-μ3)
                pR = a1 * pnorm(R-μ1) + a2 * pnorm(R-μ2) + a3 * pnorm(R-μ3)

                if sign((R-L)/(pR-pL) * (a1 * dnorm(L-μ1) + a2 * dnorm(L-μ2) + a3 * dnorm(L-μ3)) - 1) ≠ -1
                    println("🤔🤔")
                end

                # L^* is the smallest L that satisfies δ_1 < 0 < δ_2. (supp. p.30)
                # Thus returns right-trunctation point
                if δ_1 < 0 < δ_2 && L_old ≤ L

                    # if L^* is achieved at i==1, it means there is no truncation.
                    if i == 1
                        return(-1)
                    end

                    return(R)
                end

                # increment L
                if L_old ≤ L
                    L_old = L
                end

                # if L > 0, it means roughly that L > m_ (m_≤0)  (supp. p.29 table and p.28)
                if L > 0
                    return(-1)
                end

            end
        end
    end

    # If a1 == a2 and μ1 != -μ2, then return 0
    # : do not need to assess in this case, since it can WLOG be reduced to a symmetric case)
    return(0)
end

"""
    clusterpath_pop(a2::T, μ2::T; a1=1-a2, μ1=-μ1, a3=0., μ3=0, σ=1) where T <: Real

Inputs:
    a2: weight on 2nd normal distribution
    μ2: mean of 2nd normal distribution
Returns:
    A dictionary of weights (`a1`, `a2`), means (`μ1`, `μ2`) and
    description of truncation (`L*`, `R*`) and split point (`s`).
"""
function clusterpath_pop(a2::T, μ2::T; a1=1-a2, μ1=-μ2, a3=0., μ3=0, σ=1.) where T <: Real
    # search right truncation point
    R = find_truncation(a2, μ2; a1=a1, μ1=μ1, a3=a3, μ3=μ3, σ=σ)

    # R == 0 or R == -1
    if R == 0 || R == -1
        d = Array(μ1:1e-04:μ2)
        # d[argmin(abs.(a1 * dnorm(d.-μ1) .- a2 * dnorm(d.-μ2)))]  ## ???????
        # d[argmin(a1 * dnorm(d.-μ1) .+ a2 * dnorm(d.-μ2))]

        return(Dict("a1" => a1, "a2" => a2,
                    "mu1" => μ1, "mu2" => μ2,
                    "s" => NaN, "L*" => NaN, "R*" => NaN))
    end

    # If R is a proper right-truncation point, then find split point
    L_lb = -3 * max(R, 3)  # lower bound of searching grid of L
    L_ub = R - 0.01  # upper bound
    L = Array(L_lb:5e-04:L_ub)

    _, _, _, _, μ_LR = cond_mean_on_LR(L, R, a1, a2, a3, μ1, μ2, μ3)
    δs = sign.(μ_LR .- (L .+ R)/2)

    # L^* is the smallest L that satisfies δ_1 < 0 < δ_2
    L = L[2:end][δs[2:end] .* δs[1:end-1] .< 0][1]

    d = Array(L:0.001:R)
    split = find_split(d, a1, a2, μ1, μ2, L, R)

    return(Dict("a1" => a1, "a2" => a2,
                "mu1" => μ1, "mu2" => μ2,
                "s" => split, "L*" => L, "R*" => R))
end


text:
using BinaryBuilder

name = "NOMAD"
version = v"4.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gist.github.com/amontoison/06dac8b63424854f754264597af6b09e/raw/942aba554a82800d0c7438288d3da1a827ef2974/NOMAD.zip", "a6653af375be8006e742239af3914ba48034015d57a2be7a07fc14c6ff245d1b"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd "${WORKSPACE}/srcdir/NOMAD"
# atomic_patch -p1 "${WORKSPACE}/srcdir/patches/nomad_openmp.patch"
# atomic_patch -p1 "${WORKSPACE}/srcdir/patches/sgtelib_openmp.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/cache_corrections.patch"
if [[ "${target}" == *-musl* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/include_sys_time_missing_timeval_musl.patch"
elif [[ "${target}" == *-apple-* ]] || [[ "${target}" == *-freebsd* ]]; then
    CC=gcc
    CXX=g++
fi
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake -DCMAKE_BUILD_TYPE=Release ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libnomadInterface", :libnomadInterface),
    LibraryProduct("libnomadAlgos", :libnomadAlgos),
    LibraryProduct("libnomadEval", :libnomadEval),
    LibraryProduct("libnomadUtils", :libnomadUtils),
    LibraryProduct("libsgtelib", :libsgtelib),
    ExecutableProduct("nomad", :nomad),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8.1.0")


text:
using Pochoir
using Test

@testset "Pochoir.jl" begin
    # Write your tests here.
end


text:
# popup selection widget
# default behavior is a simple scrollable box of strings that is selectable
# quickselect enable a searchbox. Users can type in a string and the cursor will jump to the first item matching
# that as a prefix

# Additional modification of the behavior:
# if substr is enabled, it'd search for any substr instead of a faster startswith
# if hideunmatched is enabled, any choice that doesn't match will be hidden. deleting the search string will revert
# if sortmatched is enabled (usually in conjunction with substr and but not hideunmatched), a levenstein distance
#  score will be generated and the result sorted according to the match score
# if allownew is enabled, the search box is also an entry box to enter new text, as long as there is no match.
#   Use trailing space to disambiguate. They will be stripped afterwards.

# But if no additional modification is used, tab-completion is enabled.

POPUPQUICKSELECT    = 1
POPUPSUBSTR         = 2
POPUPHIDEUNMATCHED  = 4
POPUPSORTMATCHED    = 8
POPUPALLOWNEW       = 16

defaultPopupHelpText = utf8("""
arrows : move cursor
home   : move to top
end    : move to bottom
enter  : select
""")

defaultPopupQuickHelpText = utf8("""
arrows : move item cursor
home   : move to top
end    : move to bottom
enter  : select

Search box:
ctrl-a : move search cursor to start
ctrl-e : move search cursor to end
ctrl-k : empty search entry
ctrl-r : Toggle insertion/overwrite mode

ctrl-n : move to the next matched item
ctrl-p : move to the previous matched item
""")

type TwPopupData
    choices::Array{UTF8String,1}
    datalist::Array{Any, 1}
    maxchoicelength::Int
    searchbox::Any
    currentLine::Int
    currentLeft::Int
    currentTop::Int
    selectmode::Int
    helpText::UTF8String
    TwPopupData( arr::Array{UTF8String,1} ) = new( arr, Any[], maximum( map( z->length(z), arr ) ), nothing, 1, 1, 1, 0, utf8("") )
end
TwPopupData{ T<:AbstractString}( arr::Array{T, 1 } ) = TwPopupData( map( x->utf8( x ), arr ) )

# the ways to use it:
# standalone panel
# as a subwin as part of another widget (see next function)
# w include title width, if it's shown on the left
function newTwPopup( scr::TwObj, arr::Array{Symbol,1};
        posy::Any=:center,posx::Any=:center,
        title = utf8(""), maxwidth = 50, maxheight = 15, minwidth = 20,
        quickselect = false, substrsearch=false, hideunmatched=false, sortmatched=false, allownew=false )

    return( newTwPopup( scr, map(x->utf8(string(x)),arr),
        posy=posy,posx=posx,
        title=title,maxwidth=maxwidth,maxheight=maxheight,minwidth=minwidth,
        quickselect=quickselect,substrsearch=substrsearch,
        hideunmatched=hideunmatched,sortmatched=sortmatched,allownew=allownew ) )
end

function newTwPopup{T<:AbstractString}( scr::TwObj, arr::Array{T,1};
        posy::Any=:center,posx::Any=:center,
        title = utf8(""), maxwidth = 50, maxheight = 15, minwidth = 20,
        quickselect = false, substrsearch=false, hideunmatched=false, sortmatched=false, allownew=false )
    obj = TwObj( TwPopupData(arr), Val{ :Popup } )
    obj.box = true
    obj.title = title
    obj.borderSizeV= 1
    obj.borderSizeH= 1
    if quickselect
        obj.data.selectmode |= POPUPQUICKSELECT
    end
    if substrsearch
        obj.data.selectmode |= POPUPQUICKSELECT | POPUPSUBSTR
    end
    if hideunmatched
        obj.data.selectmode |= POPUPQUICKSELECT | POPUPHIDEUNMATCHED
    end
    if sortmatched
        obj.data.selectmode |= POPUPQUICKSELECT | POPUPSUBSTR | POPUPSORTMATCHED
    end
    if allownew
        obj.data.selectmode |= POPUPQUICKSELECT | POPUPSUBSTR | POPUPHIDEUNMATCHED | POPUPALLOWNEW
    end
    usedatalist = popup_use_datalist( obj )
    if usedatalist
        rebuild_popup_datalist( obj )
    end
    if obj.data.selectmode & POPUPQUICKSELECT != 0
        obj.data.helpText = defaultPopupQuickHelpText
    else
        obj.data.helpText = defaultPopupHelpText
    end

    h = 2 + min( length( arr ), maxheight )
    w = 2 + max( min( max( length( title ), obj.data.maxchoicelength ), maxwidth ), minwidth )

    link_parent_child( scr, obj, h, w, posy, posx )

    obj.data.searchbox = newTwEntry( obj, UTF8String; width=minwidth, posy=:bottom, posx = 1, box=false )
    obj.data.searchbox.title = utf8("?")
    obj.data.searchbox.hasFocus = false # so it looks dimmer than main cursor
    obj
end

function popup_use_datalist( o::TwObj )
    o.data.selectmode & POPUPHIDEUNMATCHED != 0 || o.data.selectmode & POPUPSORTMATCHED != 0
end

function rebuild_popup_datalist( o::TwObj{TwPopupData} )
    o.data.datalist = Any[]
    for (i, c) in enumerate( o.data.choices )
        searchstring = c
        push!( o.data.datalist, Any[lowercase( searchstring ), c, i, 0.0 ] )
    end
end

function draw( o::TwObj{TwPopupData} )
    werase( o.window )
    if o.box
        box( o.window, 0,0 )
    end
    if !isempty( o.title ) && o.box
        mvwprintw( o.window, 0, (@compat round(Int, ( o.width - length(o.title) )/2 )), "%s", o.title )
    end
    starty = o.borderSizeV
    viewContentHeight = o.height - o.borderSizeV * 2
    viewContentWidth  = o.width - o.borderSizeH * 2
    usedatalist = popup_use_datalist( o )
    if usedatalist
        n = length( o.data.datalist )
    else
        n = length( o.data.choices )
    end
    for r in o.data.currentTop:min( o.data.currentTop + viewContentHeight-1, n )
        flag = 0
        if r == o.data.currentLine
            flag = A_BOLD | COLOR_PAIR( o.hasFocus ? 15 : 30 )
        end
        if usedatalist
            s = o.data.datalist[r][2]
        else
            s = o.data.choices[r]
        end
        s = substr_by_width( s, o.data.currentLeft - 1, viewContentWidth )

        wattron( o.window, flag )
        mvwprintw( o.window, r - o.data.currentTop + starty, o.borderSizeH, "%s", s )
        wattroff( o.window, flag )
    end
    if o.data.selectmode & POPUPQUICKSELECT != 0
        draw( o.data.searchbox )
    end
end

function popup_search_next( o::TwObj{TwPopupData}, step::Int, trivialstop::Bool )
    st = o.data.currentLine
    tmpstr = lowercase(o.data.searchbox.data.inputText)
    if length(tmpstr) == 0
        TermWin.beep()
        return 0
    end

    usedatalist = popup_use_datalist( o )
    if usedatalist
        n = length( o.data.datalist )
    else
        n = length( o.data.choices )
    end

    i = trivialstop ? st : ( mod( st-1+step, n ) + 1 )
    while true
        if o.data.selectmode & POPUPSUBSTR != 0
            if usedatalist
                if contains( o.data.datalist[i][1], tmpstr )
                    o.data.currentLine = i
                    return i
                end
            else
                if contains( lowercase( o.data.choices[i] ), tmpstr )
                    o.data.currentLine = i
                    return i
                end
            end
        else
            if startswith( lowercase( o.data.choices[i] ), tmpstr )
                o.data.currentLine = i
                return i
            end
        end
        i = mod( i-1+step, n ) + 1
        if i == st
            TermWin.beep()
            return 0
        end
    end
end

function update_popup_score( o::TwObj{TwPopupData} )
    searchterm = o.data.searchbox.data.inputText
    needx = o.data.maxchoicelength

    l1 = length(searchterm)
    usedatalist = popup_use_datalist( o )

    if usedatalist
        prevchoice = ""
        if length( o.data.datalist ) >= o.data.currentLine >= 1
            prevchoice = o.data.datalist[ o.data.currentLine ][2]
        else
            o.data.currentLine = 1
        end
        if l1 == 0
            rebuild_popup_datalist( o )
            for (i, row) in enumerate( o.data.datalist)
                if row[2] == prevchoice
                    o.data.currentLine = i
                end
            end
        else
            if o.data.selectmode & POPUPHIDEUNMATCHED != 0
                o.data.datalist = Any[]
                if o.data.selectmode & POPUPSUBSTR != 0
                    for (i,c) in enumerate( o.data.choices )
                        if contains( lowercase( c ), lowercase( searchterm ) )
                            push!( o.data.datalist, Any[ lowercase( c ), c, i, search( lowercase( c ), lowercase( searchterm ) )[1] + length( c ) / needx ] )
                        end
                    end
                else
                    for (i,c) in enumerate( o.data.choices )
                        if startswith( lowercase( c ), lowercase( searchterm ) )
                            push!( o.data.datalist, Any[ lowercase( c ), c, i, length(c) ] )
                        end
                    end
                end
                if o.data.selectmode & POPUPSORTMATCHED != 0
                    sort!( o.data.datalist, lt=(x,y)-> x[4] < y[4] )
                    for (i,row) in o.data.datalist
                        if prevchoice == row[2]
                            o.data.currentLine = i
                        end
                    end
                end
            else # show everything
                # sort
                if o.data.selectmode & POPUPSORTMATCHED != 0
                    if o.data.selectmode & POPUPSUBSTR != 0
                        for row in o.data.datalist
                            ld = levenstein_distance( lowercase( searchterm ), row[1] )
                            l2 = length( row[1] )
                            minld = abs( l1 - l2 )
                            maxld = max( l1, l2 )
                            # ld closer to the theoretical minimum should be deemed almost as good as a full match
                            normld = ( ld - minld + 1) / (maxld - minld + 1) * (minld + 1 )
                            # finding the search term in the later part of a string should have a small penalty
                            substrpenalty = needx
                            r  = search( row[1], searchterm )
                            if length(r) != 0
                                substrpenalty = r.start
                            end
                            row[ 4] = ld + normld * 0.5 + substrpenalty
                        end
                    # prefix-based
                    else
                        for row in o.data.datalist
                            if startswith( row[1], lowercase( searchterm ) )
                                row[4] = length( row[1] )
                            else
                                row[4] = length( row[1] ) + needx
                            end
                        end
                    end
                    sort!( o.data.datalist, lt=(x,y)-> x[4] < y[4] )
                    o.data.currentLine = 1
                    o.data.currentTop = 1
                else # don't sort, but jump to the next match one
                    popup_search_next( o, 1, true )
                end
            end
        end
    else # just jump to the first term with the matched
        popup_search_next( o, 1, true )
    end
end

function inject( o::TwObj{TwPopupData}, token )
    dorefresh = false
    retcode = :got_it # default behavior is that we know what to do with it

    viewContentWidth = o.width - o.borderSizeH * 2
    viewContentHeight = o.height - 2 * o.borderSizeV

    usedatalist = popup_use_datalist( o )
    tabcomplete = ( o.data.selectmode & POPUPQUICKSELECT != 0 ) && (o.data.selectmode & POPUPSUBSTR == 0 )

    checkTop = () -> begin
        if o.data.currentTop > o.data.currentLine
            o.data.currentTop = o.data.currentLine
        elseif o.data.currentLine - o.data.currentTop > viewContentHeight - 1
            o.data.currentTop = o.data.currentLine - viewContentHeight + 1
        end
    end
    moveby = n -> begin
        if o.data.selectmode & POPUPHIDEUNMATCHED != 0
            sz = length( o.data.datalist )
        else
            sz = length( o.data.choices )
        end

        o.data.currentLine = max(1, min( sz, o.data.currentLine + n) )
        checkTop()
        return true
    end

    if o.data.selectmode & POPUPQUICKSELECT != 0 && token != :F1
        inputText = o.data.searchbox.data.inputText
        result = inject( o.data.searchbox, token )

        if result == :got_it
            if inputText != o.data.searchbox.data.inputText
                update_popup_score( o )
                checkTop()
            end
            refresh( o )
            return result
        end
    end

    if token == :esc
        retcode = :exit_nothing
    elseif token == :up
        dorefresh = moveby(-1)
    elseif token == :down
        dorefresh = moveby(1)
    elseif token == :left
        if o.data.currentLeft > 1
            o.data.currentLeft -= 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :right
        if o.data.currentLeft + viewContentWidth < o.data.maxchoicelength
            o.data.currentLeft += 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :shift_left # TODO ctrl-left
        if o.data.currentLeft > 1
            o.data.currentLeft -= 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :ctrlshift_left # TODO ctrl-left
        if o.data.currentLeft > 1
            o.data.currentLeft = 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :shift_right
        if o.data.currentLeft + viewContentWidth < o.data.maxchoicelength
            o.data.currentLeft += 1
            dorefresh = true
        else
            beep()
        end
    elseif token == :ctrlshift_right
        if o.data.currentLeft + viewContentWidth < o.data.maxchoicelength
            o.data.currentLeft = o.data.maxchoicelength - viewContentWidth
            dorefresh = true
        else
            beep()
        end
    elseif token == :pageup
        dorefresh = moveby( -viewContentHeight )
    elseif token == :pagedown
        dorefresh = moveby( viewContentHeight )
    elseif token == :ctrl_n
        popup_search_next( o, 1, false )
        checkTop()
        dorefresh = true
    elseif token == :ctrl_p
        popup_search_next( o, -1, false )
        checkTop()
        dorefresh = true
    elseif token == :KEY_MOUSE
        (mstate,x,y, bs ) = getmouse()
        if mstate == :scroll_up
            dorefresh = moveby( -(@compat round(Int, viewContentHeight/10 )) )
        elseif mstate == :scroll_down
            dorefresh = moveby( @compat round(Int, viewContentHeight/10 ) )
        elseif mstate == :button1_pressed && o.data.trackLine
            (rely, relx) = screen_to_relative( o.window, y, x )
            if 0<=relx<o.width && 0<=rely<o.height
                o.data.currentLine = o.data.currentTop + rely - o.borderSizeH + 1
                dorefresh = true
            else
                retcode = :pass
            end
        end
    elseif token == :tab && tabcomplete
        # auto-input the search text to be the longest common prefix of the current line and the next line,
        # as long as the current search text content can be appended.
        nextstr = ""
        currstr = ""
        if usedatalist
            if o.data.currentLine < length( o.data.datalist )
                currstr = o.data.datalist[ o.data.currentLine     ][2]
                nextstr = o.data.datalist[ o.data.currentLine + 1 ][2]
            end
        else
            if o.data.currentLine < length( o.data.choices )
                currstr = o.data.choices[ o.data.currentLine    ]
                nextstr = o.data.choices[ o.data.currentLine + 1]
            end
        end
        lcp = longest_common_prefix( currstr, nextstr )
        if startswith( lcp, o.data.searchbox.data.inputText )
            o.data.searchbox.data.inputText = lcp
            inject( o.data.searchbox, :ctrl_e ) # move the cursor to the end
            dorefresh = true
        else
            beep()
        end
    elseif  token == :home
        if o.data.currentTop != 1 || o.data.currentLeft != 1 || o.data.currentLine != 1
            o.data.currentTop = 1
            o.data.currentLeft = 1
            o.data.currentLine = 1
            dorefresh = true
        else
            beep()
        end
    elseif in( token, Any[ Symbol("end") ] )
        if usedatalist
            n = length( o.data.datalist )
        else
            n = length( o.data.choices )
        end
        if o.data.currentLine != n
            o.data.currentLine = n
            checkTop()
            dorefresh = true
        else
            beep()
        end
    elseif token == :enter || token == Symbol( "return" )
        if usedatalist
            if o.data.currentLine <= length( o.data.datalist )
                o.value = o.data.datalist[ o.data.currentLine ][ 2]
                retcode = :exit_ok
            elseif o.data.selectmode & POPUPALLOWNEW != 0
                o.value = strip( o.data.searchbox.data.inputText )
                retcode = :exit_ok
            end
        else
            if o.data.currentLine <= length( o.data.choices )
                o.value = o.data.choices[ o.data.currentLine ]
                retcode = :exit_ok
            end
        end
    else
        retcode = :pass # I don't know what to do with it
    end

    if dorefresh
        refresh(o)
    end

    return retcode
end

function helptext( o::TwObj{TwPopupData} )
    s = o.data.helpText
    tabcomplete = ( o.data.selectmode & POPUPQUICKSELECT != 0 ) && (o.data.selectmode & POPUPSUBSTR == 0 )
    if tabcomplete
        s *= "tab    : tab-completion"
    end
    s
end


text:
export regex_neg_lookahead, regex_tempered_greedy
regex_neg_lookahead(e,match=AnyChar()) =
    !Repeat(Sequence(NegativeLookahead(e),match))
regex_tempered_greedy(s,e, flags="s"; withend=true) =
    Sequence(
        2,s,
        !Repeat(NegativeLookahead(e),AnyChar()),
        ( withend ? e : Always()))
    # Regex("^"*regex_string(s)*"((?:(?!"*regex_string(e)*").)*)"*
    #       ( withend ? regex_string(e) : ""),flags)

export alternate, alternate_stop
alternate_stop(x,delim,stop;kw...) =
    alternate(seq(NegativeLookahead(stop), x; transform=2),
              seq(NegativeLookahead(stop), delim; transform=2);
              kw...)

alternate(x::Vector, delim; kw...) = alternate(Either(x...), delim; kw...)
"""
optimized repeated alternations of `x``delim`, optionally starting/ending with `delim`. `delim` `is agg`ed as right borders.
`delim` can be discarded in agg(result,missing,delim).

if `agg` is nothing, default is to aggregate delim after match is `result_type(delim) <: result_type(x)`, if not missing.
"""
function alternate(x, delim;
                   agg = nothing,
                   kw...)
    @show T, S = result_type(typeof(x)), result_type(typeof(delim))
    af = if agg === nothing
        if S <: T
            ( r, xmatch, delimmatch ) -> begin
                xmatch !== missing && push!(r,xmatch)
                delimmatch !== missing && push!(r,delimmatch)
                r::Vector{T}
            end
        else
            ( r, xmatch, delimmatch ) -> begin
                xmatch !== missing && push!(r,xmatch)
                r::Vector{T}
            end
        end
    else
        agg
    end
    function tf(v)::Vector{T}
        ## @show v,i
        r = T[]
        if isempty(v[2])
            af(r,v[1],v[3])
        else
            ms = v[2]
            af(r,v[1],ms[1][1])
            for i in 2:lastindex(ms)
                af(r, ms[i-1][2],ms[i][1])
            end
            af(r, ms[end][2],v[3])
        end
        r::Vector{T}
    end

    ## todo: factor out this transform condition!!
    r=Sequence(tf,
             Optional(x; default=missing),
             Repeat(seq(delim, x)),
               Optional(delim;default=missing))
    @show result_type(r)
    r
end

export Repeat_delim
Repeat_delim(x, delim; kw...) =
    Repeat_delim(
        promote_type(result_type(x),result_type(delim)),
        x, delim; kw...)
function Repeat_delim(
    T::Type, x, delim;
    log=false,repf=Repeat,
    transform=(v) -> v,
    transform_each=(v) -> v, kw...)
    x = parser(x)
    delim = parser(delim)
    function t(v)
        L = vcat(v...)
        transform(map(p -> transform_each(p,i),  L  ),i)
    end
    seq(Vector{T},
        Optional(delim; default=T[], log=log),
        repf(Vector{T},
             seq(x, delim; log=log); log=log,
             transform = (v) -> vcat([ [x...] for x in v ]...)),
        Optional(x; default=T[], log=log)
        ; log=log,
        ## todo: factor out this transform condition!!
        transform = (t)
        , kw...)
end

export Repeat_delim_par
function Repeat_delim_par(x, delim; repf=Repeat, transform=(v) -> v, transform_each=v -> v, kw...)
    x = parser(x)
    delim = parser(delim)
    T = result_type(typeof(x))
    D = result_type(typeof(delim))
    seq(Vector{T},
        Optional(v -> D[v], delim),
        repf(Sequence(1, x, delim);
             transform=(v) -> v),
        Optional(v -> T[v], x)
        ;
        ## todo: factor out this transform condition!!
        transform = (v)  -> transform(
            map(p -> transform_each(p),
                vcat(v[2:end]...)),i)
        , kw...)
end

"""
tokenize(x, str; delta=200, errorfile=nothing)

Tokenize string or iterator `str` with parser `x`.
"""
function tokenize(x, str; partial=:error)
    i=firstindex(str)
    till=lastindex(str)
    r, i_ = tryparsenext(x, str, i, till, TextParse.default_opts)
    if i_<=till
        if partial isa AbstractString ## remove?
            make_org(s) = replace(s, r"^\*"m => " *")
            open(partial, "a") do io
                println(io, "* incomplete parsing stopped at $i_ ")
                println(io, "error at")
                println(io, make_org(str[min(i_,end):min(end, nextind(str,i_,200))]))
                println(io, "** data")
                println(io, make_org(str))
            end
        elseif partial == :warn
            @warn "incomplete parsing stopped at $i_ " str[min(i_,end):min(end, nextind(str,i_,200))]
        elseif partial == :error
            throw(PartialMatchException(i_, str, x))
        elseif partial == :nothing
            return nothing
        end
    end
    if isnull(r)
        if partial == :error
            error("no match")
        elseif partial == :warn
            @warn "no match"
        else
            nothing
        end
    else
        get(r)
    end
end

export splitter
splitter(S, parse; transform_split = v -> tokenize(S, v), kw...) =
    splitter(Regex(regex_string(S)), parse;
             transform_split = transform_split, kw...)

function splitter(## R::Type,
                  split::Transformation{Regex,S},
                  parse::CombinedParsers.ParserTypes;
                  log=false,
                  transform = (v) -> v) where {S}
    @warn "todo: using old regex splitting..."
    transform_split = split.transform ## (v,i) -> v
    T = result_type(parse)
    R = promote_type(S,T)
    function tpn(str, i, n, opts) ## from util.jl:_split
        ## @show str
        ## @show R
        strs = Vector{R}(undef, 0)#[]
        lstr = str[i:min(end,n)]
        r = eachmatch(split.parser, lstr)
        j = 0
        for m in r
            if j <= m.match.offset
                ## m.match.offset  is indexed at 0!!
                ## @show lstr nextind(lstr,j) m.match.offset m.match
                before = SubString(lstr,nextind(lstr,j),prevind(lstr, m.match.offset + (transform_split===nothing ? sizeof(m.match) : 1)))
                log && @info "before" before
                push!(strs, (tokenize(parse, before))) # , i+nextind(lstr,j))) ## todo pass pos!
            end
            if transform_split!==nothing
                log && @info "split" before
                push!(strs, ( transform_split(m, i))) # , i+j) )
            end
            j = m.match.offset + sizeof(m.match) # =.ncodeunits
        end
        ## j = prevind(lstr,j)
        if j <= n-i
            after = SubString(str,i+j,min(lastindex(str),n))
            log && @info "after" after
            push!(strs,
                  (tokenize(parse, after))) ## , i+j)) ## todo pass pos!
        end
        result = transform(strs,i)
        ## error()
        log && @info "split" lstr strs result i j n
        return Nullable(result), nextind(str,n)
    end
    CustomParser(tpn, R)
end


text:
using BenchmarkTools
using ITensors

function main(; Nrange = 1:8)
  # Don't warn about large tensor orders
  @disable_warn_order begin

  for N in Nrange
    println("#################################################")
    println("# order = ", 2 * N)
    println("#################################################")
    println()

    d = 1
    i = Index(QN(0, 2) => d, QN(1, 2) => d)
    is = IndexSet(n -> settags(i, "i$n"), N)
    A = randomITensor(is'..., dag(is)...)
    B = randomITensor(is'..., dag(is)...)

    # Use standard algorithm, without combining first
    println("Contract:")
    C_contract = @btime $A' * $B samples = 5
    println()

    # Reshape the ITensors into matrices before contracting
    println("Combine then contract:")
    ITensors.enable_combine_contract()
    C_combine_contract = @btime $A' * $B samples = 5
    ITensors.disable_combine_contract()
    println()

    @show C_contract ≈ C_combine_contract
    println()
  end

  end # @disable_warn_order
end


text:
"""

Type for storing steady state response (SSR) data.

The minimum required information for SSR recording is stored in the type.
Additional information can be stored in the `processing` field.
To facilitate processing, specific names are used in the processing dictionary.

#### Fields

* `data`: contains the recorded data
* `trigers`: contains information about timing for creation of epochs
* `system_codes`: contains system information
* `samplingrate`: the sampling rate of the data
* `modulationrate`: the modulation rate of the stimulus
* `reference_channel`: the channel the data has been referenced to
* `file_path` and `file_name`: where the file was read in from
* `channel_names`: the names of the channels
* `processing`: dictionary type to store analysis
* `header`: additional information read from the file

#### Additional `processing` fields
The following standard names are used when saving data to the processing dictionary.

* `Name`: The identifier for the participant
* `Side`: Side of stimulation
* `Carrier_Frequency`: Carrier frequency of the stimulus
* `Amplitude`: Amplitude of the stimulus
* `epochs`: The epochs extracted from the recording
* `sweeps`: The extracted sweeps from the recording

#### Example

Put an example here

```julia
s = SSR("filename")
```

"""
mutable struct SSR
    data::Array
    sensors::Array{Sensor}
    triggers::Dict
    system_codes::Dict
    samplingrate::typeof(1.0u"Hz")
    modulationrate::typeof(1.0u"Hz")
    reference_channel::Array{AbstractString,1}
    file_path::AbstractString
    file_name::AbstractString
    processing::Dict
    header::Dict
end

#######################################
#
# SSR info
#
#######################################

"""
Return the sampling rate of a steady state type.
If no type is provided, the sampling rate is returned as a floating point.

#### Example

Return the sampling rate of a recording

```julia
s = read_SSR(filename)
samplingrate(s)
```
"""
samplingrate(t, s::SSR) = convert(t, ustrip(s.samplingrate))
samplingrate(s::SSR) = samplingrate(AbstractFloat, s)

"""
Return the modulation rate of a steady state type.
If no type is provided, the modulation rate is returned as a floating point.

#### Example

Return the modulation rate of a recording

```julia
s = read_SSR(filename)
modulationrate(s)
```
"""
modulationrate(t, s::SSR) = convert(t, ustrip(s.modulationrate))
modulationrate(s::SSR) = modulationrate(AbstractFloat, s)

"""
Return the names of sensors in SSR measurement.

#### Example

```julia
s = read_SSR(filename)
channelnames(s)
```
"""
channelnames(s::SSR) = labels(s.sensors)

"""
Change the names of sensors in SSR measurement.

#### Example

```julia
s = read_SSR(filename)
channelnames(s, 1, "Fp1")
```
"""
function channelnames(s::SSR, i::Int, l::S) where {S<:AbstractString}

    s.sensors[i].label = l
    return s
end
function channelnames(s::SSR, l::AbstractVector{S}) where {S<:AbstractString}
    for li = 1:length(l)
        s = channelnames(s, li, l[li])
    end
    return s
end

#######################################
#
# Show
#
#######################################

import Base.show
function Base.show(io::IO, a::SSR)
    time_length = round.(size(a.data, 1) / samplingrate(a) / 60)
    println(
        io,
        "SSR measurement of $time_length mins with $(size(a.data,2)) channels sampled at $(a.samplingrate)",
    )
    println(io, "  Modulation frequency: $(a.modulationrate )")

    if haskey(a.processing, "Amplitude")
        println(io, "  Stimulation amplitude: $(a.processing["Amplitude"]) dB")
    end
    if haskey(a.processing, "Name")
        println(io, "  Participant name: $(a.processing["Name"] )")
    end
    if haskey(a.processing, "Side")
        println(io, "  Stimulation side: $(a.processing["Side"] )")
    end
    if haskey(a.processing, "Carrier_Frequency")
        println(io, "  Carrier frequency: $(a.processing["Carrier_Frequency"] ) Hz")
    end

end

#######################################
#
# Operate on type
#
#######################################

import Base.hcat
"""
Append one SSR type to another, simulating a longer recording.

#### Example

```julia
hcat(a, b)
```
"""
function hcat(a::SSR, b::SSR)

    if channelnames(a) != channelnames(b)
        throw(
            ArgumentError(
                string("Channels do not match $(channelnames(a)) != $(channelnames(b))"),
            ),
        )
    end

    if haskey(a.processing, "epochs")
        @warn("Epochs have already been extracted and will no longer be valid")
    end
    if haskey(a.processing, "statistics")
        @warn("Statistics have already been calculated and will no longer be valid")
    end

    @debug(
        "Appending two SSRs with $(size(a.data, 2)) .& $(size(b.data, 2)) channels and lengths $(size(a.data, 1)) $(size(b.data, 1))"
    )

    join_triggers(a, b)
    a.data = [a.data; b.data]

    return a
end

"""
Append the trigger information of one SSR type to another.
Places the trigger information at the end of first file

#### Example

```julia
join_triggers(a, b)
```
"""
function join_triggers(a, b; offset = size(a.data, 1))

    a.triggers["Index"] = [a.triggers["Index"]; (b.triggers["Index"] .+ offset)]
    a.triggers["Code"] = [a.triggers["Code"]; b.triggers["Code"]]
    a.triggers["Duration"] = [a.triggers["Duration"]'; b.triggers["Duration"]']'

    a
end

#######################################
#
# Manipulate channels
#
#######################################

"""
Add a channel to the SSR type with specified channel names.

#### Example

Add a channel called `Merged`

```julia
s = read_SSR(filename)
new_channel = mean(s.data, 2)
s = add_channel(s, new_channel, "Merged")
```
"""
function add_channel(a::SSR, data::Vector, chanLabel::AbstractString; kwargs...)

    @info("Adding channel $chanLabel")

    a.data = hcat(a.data, data)
    push!(a.sensors, Electrode(chanLabel, Talairach(NaN, NaN, NaN), Dict()))

    return a
end

"""
Remove specified channels from SSR.

#### Example

Remove channel Cz and those in the set called `EEG_Vanvooren_2014_Right`

```julia
a = read_SSR(filename)
remove_channel!(a, [EEG_Vanvooren_2014_Right, "Cz"])
```
"""
function remove_channel!(
    a::SSR,
    channel_names::Array{S};
    kwargs...,
) where {S<:AbstractString}
    @debug("Removing channels $(join(channel_names, " "))")
    remove_channel!(
        a,
        Int[something(findfirst(isequal(c), channelnames(a)), 0) for c in channel_names],
    )
end

function remove_channel!(a::SSR, channel_name::S; kwargs...) where {S<:AbstractString}
    @debug("Removing channel $(channel_name)")
    remove_channel!(a, [channel_name])
end

remove_channel!(a::SSR, channel_names::Int; kwargs...) =
    remove_channel!(a, [channel_names]; kwargs...)

function remove_channel!(a::SSR, channel_idx::Array{Int}; kwargs...)

    channel_idx = channel_idx[channel_idx.!=0]

    @debug("Removing channel(s) $channel_idx")
    if any(channel_idx .== 0)
        @warn("Failed to remove a channel")
    end

    keep_idx = [1:size(a.data)[end];]
    for c in sort(channel_idx, rev = true)
        try
            splice!(keep_idx, c)
        catch
            # Nothing
        end
    end

    if haskey(a.processing, "epochs")
        if size(a.processing["epochs"], 3) == size(a.data, 2)
            @debug("Removing channel(s) from epoch data")
            a.processing["epochs"] = a.processing["epochs"][:, :, keep_idx]
        end
    end
    if haskey(a.processing, "sweeps")
        if size(a.processing["sweeps"], 3) == size(a.data, 2)
            @debug("Removing channel(s) from sweep data")
            a.processing["sweeps"] = a.processing["sweeps"][:, :, keep_idx]
        end
    end

    a.data = a.data[:, keep_idx]
    a.sensors = a.sensors[keep_idx]

    return a
end

"""
Remove all channels except those requested from SSR.

#### Example

Remove all channels except Cz and those in the set called `EEG_Vanvooren_2014_Right`

```julia
a = read_SSR(filename)
keep_channel!(a, [EEG_Vanvooren_2014_Right, "Cz"])
```
"""
function keep_channel!(a::SSR, channel_names::Array{S}; kwargs...) where {S<:AbstractString}
    @info("Keeping channel(s) $(join(channel_names, " "))")
    keep_channel!(
        a,
        vec(
            round.(
                Int,
                [
                    something(findfirst(isequal(c), channelnames(a)), 0) for
                    c in channel_names
                ],
            ),
        ),
    )
end

function keep_channel!(a::SSR, channel_name::AbstractString; kwargs...)
    keep_channel!(a, [channel_name]; kwargs...)
end

function keep_channel!(a::SSR, channel_idx::AbstractVector{Int}; kwargs...)

    remove_channels = [1:size(a.data, 2);]

    channel_idx = sort(channel_idx, rev = true)
    for c in channel_idx
        splice!(remove_channels, c)
    end

    remove_channel!(a, remove_channels; kwargs...)
end

#######################################
#
# Trim channels
#
#######################################

"""
Trim SSR recording by removing data after `stop` specifed samples.

#### Optional Parameters

* `start` Remove samples before this value

#### Example

Remove the first 8192 samples and everything after 8192*300 samples

```julia
s = trim_channel(s, 8192*300, start=8192)
```
"""
function trim_channel(a::SSR, stop::Int; start::Int = 1, kwargs...)

    @info("Trimming $(size(a.data)[end]) channels between $start and $stop")

    a.data = a.data[start:stop, :]

    to_keep = findall((a.triggers["Index"] .>= start) .& (a.triggers["Index"] .<= stop))
    a.triggers["Index"] = a.triggers["Index"][to_keep]
    a.triggers["Duration"] = a.triggers["Duration"][to_keep]
    a.triggers["Code"] = a.triggers["Code"][to_keep]
    a.triggers["Index"] .-= (start - 1)

    to_keep =
        findall((a.system_codes["Index"] .>= start) .& (a.system_codes["Index"] .<= stop))
    a.system_codes["Index"] = a.system_codes["Index"][to_keep]
    a.system_codes["Duration"] = a.system_codes["Duration"][to_keep]
    a.system_codes["Code"] = a.system_codes["Code"][to_keep]
    a.system_codes["Index"] .-= (start - 1)

    return a
end

#######################################
#
# Merge channels
#
#######################################

"""
Merge `SSR` channels listed in `merge_Chans` and label the averaged channel as `new_name`

If multiple channels are listed then the average of those channels will be added.

#### Example

```julia
s = merge_channels(s, ["P6", "P8"], "P68")
```
"""
function merge_channels(
    a::SSR,
    merge_Chans::Array{S},
    new_name::S;
    kwargs...,
) where {S<:AbstractString}

    @debug("Number of original channels: $(length(channelnames(a)))")

    keep_idxs =
        vec([something(findfirst(isequal(i), channelnames(a)), 0) for i in merge_Chans])

    if sum(keep_idxs .== 0) > 0
        @warn(
            "Could not merge as these channels don't exist: $(join(vec(merge_Chans[keep_idxs .== 0]), " "))"
        )
        keep_idxs = keep_idxs[keep_idxs.>0]
    end

    @info("Merging channels $(join(vec(channelnames(a)[keep_idxs,:]), " "))")
    @debug("Merging channels $keep_idxs")

    a = add_channel(a, vec(mean(a.data[:, keep_idxs], dims = 2)), new_name; kwargs...)
end

function merge_channels(
    a::SSR,
    merge_Chans::S,
    new_name::S;
    kwargs...,
) where {S<:AbstractString}
    a = merge_channels(a, [merge_Chans], new_name; kwargs...)
end

#######################################
#
# Helper functions
#
#######################################

function assr_frequency(
    rounded_freq::Number;
    stimulation_samplingrate::Number = 32000,
    stimulation_frames_per_epoch::Number = 32768,
)

    round.(rounded_freq / (stimulation_samplingrate / stimulation_frames_per_epoch)) *
    stimulation_samplingrate / stimulation_frames_per_epoch
end

function assr_frequency(rounded_freq::AbstractVector)

    [assr_frequency(f) for f in rounded_freq]
end


text:
"""
    calc_pco2(T,S,dic,ta,pt=0.0,sit=0.0)

Calculate pCO2 for given Temperature(T), Salinity(S), DIC, Alkalinity, etc..
Efficient solver following Follows et al (2005)

```
pco2eq = atmospheric reference pCO2 level (atmospheres)
               for which to find equilibrium dic, csat
      csat = equilibrium total inorganic carbon (mol/m³)
            where 1 T = 1 metric ton = 1000 kg
      ta  = total alkalinity (eq/m³)
      pt  = inorganic phosphate (mol/m³)
      sit = inorganic silicate (mol/m³)
      T   = temperature (degrees C)
      S   = salinity (PSU)
    hg  = first guess of [H+] (10e-8 is a good cold start value)
```
"""
function calc_pco2(T,S,dic,ta,pt=0.0,sit=0.0)
    # set first guess for [H⁺]
    hg = 1.0e-8
    # convert degrees C to degrees K
    TK = T + 273.0

    # estimated concentration of borate(bt) based on salinity
    scl=S/1.80655;
    bt=0.000232 * scl/10.811

    # Coefficient algorithms as used in OCMIP2 protocols K1, K2 Millero (1995) using Mehrbach data
    k1 = 10^(-1.0*(3670.7/TK - 62.008 + 9.7944*log(TK) - 1.18e-2*S + 1.16e-4*S*S))
    k2 = 10^(-1.0*(1394.7/TK + 4.777 - 1.84e-2*S + 1.18e-4*S*S))

    # K1p, K2p, K3p, DOE (1994)
    k1p = exp(-4576.752/TK + 115.525 - 18.453*log(TK) + (-106.736/TK+ 0.69171)*sqrt(S)+ (-0.65643/TK- 0.01844)*S)
    k2p = exp(-8814.715/TK + 172.0883- 27.927*log(TK) + (-160.34/TK + 1.3566)*sqrt(S) + (0.37335/TK - 0.05778)*S)
    k3p = exp(-3070.75/TK + 18.141 + (17.27039/TK + 2.81197)*sqrt(S) + (-44.99486/TK - 0.09984)*S)

    # Kb, Millero (1995) using data from Dickson
    kb = exp((-8966.90 - 2890.53*sqrt(S) - 77.942*S + 1.728*S^1.5 - 0.0996*S*S)/TK + (148.0248 + 137.1942*sqrt(S) + 1.62142*S) + (-24.4344 - 25.085*sqrt(S) - 0.2474*S) * log(TK) + 0.053105*TK*sqrt(S))

    # Kw, Millero (1995)
    kw = exp(-13847.26/TK + 148.9652 - 23.6521*log(TK) + (118.67/TK - 5.977 + 1.0495*log(TK))*sqrt(S) - 0.01615*S)

    # fugacity, Weiss and Price, Marine Chem, 8, 347 (1990)
    TK1 = TK/100.0
    ff = exp(-162.8301 + 218.2968/TK1 + 90.9241*log(TK1) - 1.47696*(TK1*TK1) + S*(0.025695 - 0.025225*TK1 + 0.0049867*(TK1*TK1)))

    # Ksi, Millero (1995) if needed
    I = (19.924*S)/(1000 - 1.005*S)
    ksi = exp(-8904.2/TK + 117.385 - 19.334*log(TK) + (-458.79/TK + 3.5913)*sqrt(I) + (188.74/TK - 1.5998)*I + (-12.1652/TK + 0.07871)*I*I + log(1.0 - 0.001005*S))

    # First guess of [H⁺]: from last timestep *OR* fixed for cold start here iterate for accurate solution
    for i in 1:1000
        # estimate contributions to total alk from borate, silicate, phosphate
        bohg = (bt*kb)/(hg+kb)
        siooh3g = (sit*ksi)/(ksi + hg)
        denom = (hg*hg*hg) + (k1p*hg*hg) + (k1p*k2p*hg) + (k1p*k2p*k3p)
        h3po4g = (pt*hg*hg*hg)/denom
        h2po4g = (pt*k1p*hg*hg)/denom
        hpo4g = (pt*k1p*k2p*hg)/denom
        po4g = (pt*k1p*k2p*k3p)/denom
        # estimate carbonate alkalinity
        cag = ta - bohg - (kw/hg) + hg - hpo4g - 2*po4g + h3po4g - siooh3g
        # improced estimate of hydrogen ion concentration
        gamm = dic/cag
        dummy = (1 - gamm)*(1 - gamm)*k1*k1 - 4*k1*k2*(1 - 2*gamm)
        Hnew = 0.5*((gamm - 1)*k1 + sqrt(dummy))
        hg = Hnew
    end
    co2s = dic/(1 + (k1/hg) + ((k1*k2)/(hg*hg)))
    pco2 = co2s/ff
    return pco2
end


text:
struct LocationIterator{w,T}
    range::UnitRange{Int}
    index::FMIndex{w,T}
end

Base.length(iter::LocationIterator) = length(iter.range)

function Base.iterate(iter::LocationIterator, i::Int=1)
    if i > length(iter)
        return nothing
    end
    return sa_value(iter.range[i], iter.index) + 1, i + 1
end


text:
parse(::Type{Quantity{T}}, s::AbstractString; kwargs...) where {T <: Real} =
    convert(Quantity{T}, tryparse_internal(Quantity{T}, s, firstindex(s), lastindex(s), 10, true; kwargs...))

function tryparse_internal(::Type{Quantity{T}}, sbuff::Union{String,SubString{String}},
    startpos::Int, endpos::Int, base::Integer, raise::Bool) where {T<:Real}
    if isempty(sbuff)
        raise && throw(ArgumentError("input string is empty"))
        return nothing
    end

    orig_start = startpos
    orig_end   = endpos

    # Ignore leading and trailing whitespace
    while isspace(sbuff[startpos]) && startpos <= endpos
        startpos = nextind(sbuff, startpos)
    end
    while isspace(sbuff[endpos]) && endpos >= startpos
        endpos = prevind(sbuff, endpos)
    end

    # Find first character of unit specification
    unitpos = startpos
    while sbuff[unitpos] ∈ "Ee+-0123456789.," && unitpos <= endpos
        unitpos = nextind(sbuff, unitpos)
    end

    numlen = unitpos - startpos + 1
    unitlen = endpos - unitpos + 1

    if numlen > 0 && unitlen > 0
        num = parse(T, sbuff[startpos:unitpos - 1])
        suni = strip(sbuff[unitpos:endpos], ['[', ']', ' '])
        if '/' ∉ suni
            suni = replace(suni, r"[∙· *]" => '∙')
            sunits = split(suni, '∙')
            uni =  lookup_units(MechanicalUnits, Symbol(sunits[1]))
            for string_unit in sunits[2:lastindex(sunits)]
                uni *= lookup_units(MechanicalUnits, Symbol(string_unit))
            end
            return num * uni
        end
    end

    if raise
        substr = SubString(sbuff, orig_start, orig_end) # show input string in the error to avoid confusion
        if all(isspace, substr)
            throw(ArgumentError("input string only contains whitespace"))
        elseif '/' ∉ sbuff[unitpos:endpos]
            throw(ArgumentError("input string contains '/' unit. Replace input with multiplied units: 'm/s' =>  'm∙s⁻¹'"))
        else
            throw(ArgumentError("invalid quantity representation: $(repr(substr))"))
        end
    end
    return nothing
end

max_stars_repo_name: hustf/MechanicalUnits

max_issues_repo_name: hustf/MechanicalUnits


text:
let
    dims = (2,2)
    # initialize a dense grid graph
    G = initialize_grid_graph_from_vtx_grid(initialize_dense_vtx_grid(dims...))
    # construct the corresponding vertex coordinate list
    vtxs = [(i,j) for i in 1:dims[1] for j in 1:dims[2]]
    vtx_map = construct_vtx_map(vtxs,dims)
    edge_cache = construct_edge_cache(vtxs,vtx_map)
    @test validate_edge_cache(G,vtxs,edge_cache)
end
let
    grid = initialize_regular_vtx_grid()
    dims = size(grid)
    vtxs = [(i,j) for i in 1:dims[1] for j in 1:dims[2] if grid[i,j] > 0]
    sort!(vtxs, by=vtx->grid[vtx...])
    zones = [v for (v,vtx) in enumerate(vtxs) if sum(grid[clip(vtx[1]-1,1,dims[1]):clip(vtx[1]+1,1,dims[1]),clip(vtx[2]-1,1,dims[2]):clip(vtx[2]+1,1,dims[2])] .> 0) == 7]
    vtx_map = construct_vtx_map(vtxs,dims)
    expanded_zones = construct_expanded_zones(vtxs,vtx_map,zones)
    @test TaskGraphs.validate_expanded_zones(vtx_map,expanded_zones)
end
let
    grid = initialize_regular_vtx_grid(;n_obstacles_x=1,n_obstacles_y=1)
    dims = size(grid)
    vtxs = [(i,j) for i in 1:dims[1] for j in 1:dims[2] if grid[i,j] > 0]
    # 1   2   3   4   5   6
    # 7   8   9  10  11  12
    # 13  14   0   0  15  16
    # 17  18   0   0  19  20
    # 21  22  23  24  25  26
    # 27  28  29  30  31  32
    dist_mtx_map = DistMatrixMap(grid,vtxs)
    dist_matrix = get_dist_matrix(initialize_grid_graph_from_vtx_grid(grid))
    v1 = 13
    v2 = 15
    @test get_distance(dist_mtx_map,v1,v2,(1,1)) == dist_matrix[v1,v2]
    @test get_distance(dist_mtx_map,v1,v2,(1,2)) == dist_matrix[v1,v2]
    @test get_distance(dist_mtx_map,v1,v2,(2,1)) == dist_matrix[v1,v2] + 2
    @test get_distance(dist_mtx_map,v1,v2,(2,2)) == dist_matrix[v1,v2] + 2
    v3 = 18
    v4 = 20
    config = 4
    @test get_distance(dist_mtx_map,v3,v4,(2,2),config) == dist_matrix[v1,v2] + 2
end
let
    env = construct_regular_factory_world()
    get_x_dim(env)
    get_y_dim(env)
    get_cell_width(env)
    get_transition_time(env)
    get_vtxs(env)
    get_pickup_zones(env)
    get_dropoff_zones(env)
    get_obstacles(env)
    get_pickup_vtxs(env)
    get_dropoff_vtxs(env)
    get_obstacle_vtxs(env)
    GridFactoryEnvironment(env,get_graph(env))
    Base.zero(env)
    LightGraphs.edges(env)
    LightGraphs.edgetype(env)
    LightGraphs.has_edge(env,1,2)
    LightGraphs.has_vertex(env,1)
    LightGraphs.inneighbors(env,1)
    LightGraphs.is_directed(env)
    LightGraphs.ne(env)
    LightGraphs.nv(env)
    LightGraphs.outneighbors(env,1)
    LightGraphs.vertices(env)
    get_num_free_vtxs(env)
    get_free_zones(env)

    filename = "env.toml"
    open(filename,"w") do io
        TOML.print(io, TOML.parse(env))
    end
    env2 = read_env(filename)

    @test get_vtxs(env) == get_vtxs(env2)
    @test nv(get_graph(env)) == nv(get_graph(env2))
    @test ne(get_graph(env)) == ne(get_graph(env2))

    construct_factory_env_from_vtx_grid(env.vtx_map)
    # construct_random_factory_world()
    # sample_random_robot_locations(env,10)
end
let
    xdim = 4
    ydim = 4
    vtx_grid = initialize_dense_vtx_grid(xdim,ydim)
    env_graph = initialize_grid_graph_from_vtx_grid(vtx_grid)
    vtxs = [(i,j) for i in 1:xdim for j in 1:ydim]
    factory_env = GridFactoryEnvironment(
        graph = env_graph,
        x_dim=xdim,
        y_dim=ydim,
        vtxs=vtxs
    )
end
# test custom envs for different-sized robot sizes
let
    vtx_grid = initialize_dense_vtx_grid(4,4)
    #  1   2   3   4
    #  5   6   7   8
    #  9  10  11  12
    # 13  14  15  16
    env = construct_factory_env_from_vtx_grid(vtx_grid)
    for i in 1:3
        for j in 1:3
            # @show env.expanded_zones[vtx_grid[i,j]]
        end
    end
end
let
    factory_env = construct_regular_factory_world(;n_obstacles_x=1,n_obstacles_y=1)
    # 1   2   3   4   5   6
    # 7   8   9  10  11  12
    # 13  14   0   0  15  16
    # 17  18   0   0  19  20
    # 21  22  23  24  25  26
    # 27  28  29  30  31  32
    dist_matrix = get_dist_matrix(factory_env)
    dist_mtx_map = factory_env.dist_function

    v1 = 13
    v2 = 15
    @test get_distance(dist_mtx_map,v1,v2) == dist_matrix[v1,v2]
    @test get_distance(dist_mtx_map,v1,v2,(1,1)) == dist_matrix[v1,v2]
    @test get_distance(dist_mtx_map,v1,v2,(1,2)) == dist_matrix[v1,v2]
    @test get_distance(dist_mtx_map,v1,v2,(2,1)) == dist_matrix[v1,v2] + 2
    @test get_distance(dist_mtx_map,v1,v2,(2,2)) == dist_matrix[v1,v2] + 2

    v3 = 18
    v4 = 20
    config = 4
    @test get_distance(dist_mtx_map,v3,v4,(2,2),config) == dist_matrix[v1,v2] + 2
end


text:
@testset "gainobj" begin

box_n = Cuboid(-[0.5,0.5,0.5],[1,1,1])
gbox_n = GainObject{3,3}(box_n, d::Real->-d)
@test (d = rand(); gbox_n.D₀fun(d) == -d)

box_p = Cuboid([0.5,0.5,0.5],[1,1,1])
gbox_p = GainObject{3,3}(box_p, d::Real->d)
@test (d = rand(); gbox_p.D₀fun(d) == d)

gobj_vec = [gbox_n, gbox_p]

# grid: from [-2,-2,-2] to [2,2,2] with ∆ = 1.  Exclude ghost points.
l = (([-2,-1,0,1],[-2,-1,0,1],[-2,-1,0,1]), ([-1.5,-0.5,0.5,1.5],[-1.5,-0.5,0.5,1.5],[-1.5,-0.5,0.5,1.5]))
N = SVector(4,4,4)

D₀ = zeros(3*prod(N))
d = 1.0
gt₀ = @SVector [PRIM,PRIM,PRIM]
assign_pumpstr!(D₀, gobj_vec, d, gt₀, l, false, true)

# Constrcut the expected D₀.
D₀exp = zeros(3, N...)

# x-components
D₀exp[1, 2, 2:3, 2:3] .= -1.0  # gbox_n
D₀exp[1, 3, 3:4, 3:4] .= 1.0  # gbox_p

# y-components
D₀exp[2, 2:3, 2, 2:3] .= -1.0  # gbox_n
D₀exp[2, 3:4, 3, 3:4] .= 1.0  # gbox_p

# z-components
D₀exp[3, 2:3, 2:3, 2] .= -1.0  # gbox_n
D₀exp[3, 3:4, 3:4, 3] .= 1.0  # gbox_p

@test D₀exp[:] == D₀

end  # @testset "gainobj"