text:
using DataFrames
import CSV
using StatsPlots
using Statistics

tuning = DataFrame(CSV.File("hpc/outputs/tuning.csv"))

logistic = (x) -> 1.0 / (1.0 + exp(-x))
logit = (p) -> log(p/(1.0-p))

tuning.evidence = (tuning.updated ./ tuning.initial) .- 1.0
tuning.P = logistic.(tuning.evidence)

@df tuning density(:P, group=(:value), fill=(0, 0.2), legend=:topleft, frame=:box, dpi=300)
xaxis!((0,1), "Post-imputation probability")
yaxis!((0, 8))
savefig("model_performance/probabilities.png")

# Thresholding
thresholding_results = DataFrame(model=String[], rank=Int64[], AUC=Float64[], cutoff=Float64[],
    TPR=Float64[], TNR=Float64[], PPV=Float64[], NPV=Float64[],
    FNR=Float64[], FPR=Float64[], FDR=Float64[], FOR=Float64[],
    CSI=Float64[], ACC=Float64[], J=Float64[]
    )

for mod in unique(tuning.model)
    for rank in unique(tuning.rank)
        comb = tuning[(tuning.model .== mod) .& (tuning.rank .== rank), :]
        S = LinRange(minimum(comb.P), maximum(comb.P), 1000)
        TP = zeros(Int64, length(S))
        FP = zeros(Int64, length(S))
        TN = zeros(Int64, length(S))
        FN = zeros(Int64, length(S))
        for (i, s) in enumerate(S)
            pred = comb.P .>= s
            TP[i] = sum((comb.value .== pred) .& comb.value)
            FP[i] = sum(pred .> comb.value)
            TN[i] = sum((comb.value .== pred) .& .!(comb.value))
            FN[i] = sum(pred .< comb.value)
        end

        TPR = TP ./ (TP .+ FN)
        TNR = TN ./ (TN .+ FP)
        PPV = TP ./ (TP .+ FP)
        NPV = TN ./ (TN .+ FN)
        FNR = FN ./ (FN .+ TP)
        FPR = FP ./ (FP .+ TN)
        FDR = FP ./ (FP .+ TP)
        FOR = FN ./ (FN .+ TN)
        CSI = TP ./ (TP .+ FN .+ FP)
        ACC = (TP .+ TN) ./ (TP .+ TN .+ FP .+ FN)

        J = (TP ./ (TP .+ FN)) .+ (TN ./ (TN .+ FP)) .- 1.0
        best_J = last(findmax(J))
        p_cutoff = S[best_J]

        dx = [reverse(FPR)[i] - reverse(FPR)[i - 1] for i in 2:length(FPR)]
        dy = [reverse(TPR)[i] + reverse(TPR)[i - 1] for i in 2:length(TPR)]
        AUC = sum(dx .* (dy ./ 2.0))

        push!(thresholding_results,
                (mod, rank, AUC, p_cutoff,
                    TPR[best_J],
                    TNR[best_J],
                    PPV[best_J],
                    NPV[best_J],
                    FNR[best_J],
                    FPR[best_J],
                    FDR[best_J],
                    FOR[best_J],
                    CSI[best_J],
                    ACC[best_J],
                    J[best_J])
            )

        plot([0,1], [0,1], lab="", aspectratio=1, frame=:box, c=:grey, ls=:dash, dpi=300)
        yaxis!((0, 1), "True Positive Rate")
        xaxis!((0, 1), "False Positive Rate")
        plot!(FPR, TPR, c=:orange, lw=1, lab="", ls=:dash, fill=(:orange, [FPR, FPR], 0.1))
        plot!([FPR[best_J],FPR[best_J]], [FPR[best_J],TPR[best_J]], c=:orange, lab="")
        scatter!([FPR[best_J]], [TPR[best_J]], c=:orange, lab="", msw=0.0)
        p_cutoff = round(S[best_J]; digits=2)
        cutoff_text = text("P ≈ $(round(p_cutoff; digits=2))\nAUC ≈ $(round(AUC; digits=2))", :black, :left, 7)
        annotate!(([FPR[best_J] + 0.05], [FPR[best_J]], cutoff_text))
        savefig(joinpath("model_performance", "roc", "rank-$(rank)-model-$(mod).png"))
    end
end

sort!(thresholding_results, :AUC, rev=true)
CSV.write("artifacts/modelselection.csv", thresholding_results)

diags = ["AUC" => :AUC,
"Probability cutoff" => :cutoff,
"True positive rate" => :TPR,
"True negative rate" => :TNR,
"Positive pred. val." => :PPV,
"Negative pred. val." => :NPV,
"False negative rate" => :FNR,
"False positive rate" => :FPR,
"False discov. rate" => :FDR,
"False omission rate" => :FOR,
"Crit. succ. index" => :CSI,
"Accuracy" => :ACC,
"Youden's index" => :J]

thr_plot = sort(thresholding_results, :rank)

for (t, p) in diags
    plot(thr_plot[!, :rank], thr_plot[!, p], group=(thr_plot[!,:model]), legend=:outertopright, m=:circle, msw=0.0, dpi=300)
    xaxis!("Approximation rank", (1, 20))
    yaxis!(t, (0, 1))
    savefig("model_performance/metrics/svd-tuning-output-$(p).png")
end

max_stars_repo_name: viralemergence/trefle

max_issues_repo_name: viralemergence/trefle

max_forks_repo_name: viralemergence/trefle


text:
using ImageUtils, HDF5, LinearAlgebra, RegularizedLeastSquares

imagesBart = Any[]
imagesMRIReco = Any[]
imagesMRIReco2 = Any[]

rf = [1,2,3,4]

f_img  = "./reco/imgCG_bart.h5"
for i = 1:length(rf)
  h5open(f_img, "r") do file
    A = read(file, "/rf$(rf[i])_re")
    B = read(file, "/rf$(rf[i])_im")
    im_ = A.+im.*B
    push!(imagesBart, im_)
    filename = "./reco/imgCG_bart_rf$(rf[i]).png"
    exportImage(filename, abs.(im_), colormap="viridis")
  end
end

f_img  = "./reco/imgCG_mrireco_toeplitz1_oversamp2.0.h5"
for i = 1:length(rf)
  h5open(f_img, "r") do file
    A = read(file, "/rf$(rf[i])_re")
    B = read(file, "/rf$(rf[i])_im")
    im_ = A.+im.*B
    push!(imagesMRIReco, im_)
    filename = "./reco/imgCG_mrireco_toeplitz1_oversamp2.0_rf$(rf[i]).png"
    exportImage(filename, abs.(im_),colormap="viridis" )
  end
end

f_img  = "./reco/imgCG_mrireco_toeplitz0_oversamp1.25.h5"
for i = 1:length(rf)
  h5open(f_img, "r") do file
    A = read(file, "/rf$(rf[i])_re")
    B = read(file, "/rf$(rf[i])_im")
    im_ = A.+im.*B
    push!(imagesMRIReco2, im_)
    filename = "./reco/imgCG_mrireco_toeplitz0_oversamp1.25_rf$(rf[i]).png"
    exportImage(filename, abs.(im_), colormap="viridis" )
  end
end

function optimalScaling(I,Ireco)
  N = length(I)

  # This is a little trick. We usually are not interested in simple scalings
  # and therefore "calibrate" them away
  alpha = norm(Ireco)>0 ? (dot(vec(I),vec(Ireco))+dot(vec(Ireco),vec(I))) /
          (2*dot(vec(Ireco),vec(Ireco))) : 1.0
  I2 = Ireco.*alpha

  return I2
end

A1 = abs.(imagesMRIReco[1])
A1 .= A1 ./maximum(A1)

B1 = abs.(imagesBart[1])
B1 .= B1 ./maximum(B1)

for i = 1:length(rf)
  A = abs.(imagesMRIReco[i])
  A .= A ./maximum(A)

  B = abs.(imagesBart[i])
  B .= optimalScaling(A,B)

  D = A.-B
  filename = "./reco/imgCG_diff_rf$(rf[i]).png"
  exportImage(filename, D, colormap="viridis", vmin=-0.1, vmax=0.1, normalize=false )

  ###
  A .= optimalScaling(A1,A)
  D = A.-A1
  filename = "./reco/imgCG_diff_mrireco_rf$(rf[i]).png"
  exportImage(filename, D, colormap="viridis", vmin=-0.1, vmax=0.1, normalize=false )

  B .= optimalScaling(B1,B)
  D = B.-B1
  filename = "./reco/imgCG_diff_bart_rf$(rf[i]).png"
  exportImage(filename, D, colormap="viridis", vmin=-0.1, vmax=0.1, normalize=false )

  @info nrmsd(A,B)   nrmsd(abs.(imagesMRIReco[1]),A)  nrmsd(abs.(imagesBart[1]),B)

end

ext: jl


text:
#================================#
#=== 3D expansion simulations ===#
#================================#
using Pkg
Pkg.activate(".")
using FourierGPE, Test

## density with Makie
using Makie, AbstractPlotting
function dense(phi)
    ψm = xspace(phi,sim)
    density = abs2.(ψm)
    pmax = maximum(density)
    return density/pmax
end

function densityfilm(sol,Nt;file="3dquenchtrap.gif")
    dir = "media"
    !isdir(dir) && mkdir(dir)
    saveto=joinpath("media",file)

    scene = Scene()
    tindex = Makie.Node(1)
    scene = volume(lift(i -> dense(sol[i]), tindex),
    algorithm = :iso,
    color = (:mediumseagreen,0.15),
    show_axis=false,
    isovalue=3f0(.15))

    R = 70
    eyeat = Vec3f0(R,R,R)
    lookat = Vec3f0(0,0,0)

    record(scene, saveto, 1:Nt) do i
        update_cam!(scene, eyeat, lookat)
        rotate_cam!(scene, 0., -0.1, 0.)
        tindex[] = i
    end
    p = scene[end];
    return
end

#==================================#

## 3D expansion
using Plots, Interact, LaTeXStrings, FourierGPE
gr(titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

## Scaling GPE
import FourierGPE: nlin!, Lgp!

function k2_scaling(K,λ,t)
    Ks = K./λ(t)
    k2(Ks)
end

#TODO we need to learn a bit about iterators, do, map

# const yt = y'
# const zt = reshape(z,1,1,length(z))

function nlin!(dϕ,ϕ,sim::Sim{3},t)
    @unpack g,X,V0 = sim; x,y,z = X
    yt = y'
    zt = reshape(z,1,1,length(z))
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + g*abs2(dϕ)/λx(t)/λy(t)/λz(t) + V(λx(t)*x,λy(t)*yt,λz(t)*zt,t)
    kspace!(dϕ,sim)
end

function Lgp!(dϕ,ϕ,sim::Sim{3},t)
    @unpack γ,μ,espec,K = sim
    nlin!(dϕ,ϕ,sim,t)
    espec = 0.5*k2_scaling(K,λ,t)
    @. dϕ = -im*(1.0 - im*γ)*(dϕ + (espec - μ)*ϕ)
end

function r2_scaling(X,λ,λ̇,t)
    @assert eltype(sqrt.(λ(t).*λ̇(t))) <: Real
    rs = X.*sqrt.(λ(t).*λ̇(t))
    ρs = Iterators.product(rs...)
    map(x->sum(abs2.(x)),ρs)
end

function remove_phase!(ψ,t)
    ψ .*= exp.(-im*r2_scaling(X,λ,λ̇,t)/2)
end

function restore_phase!(ψ,t)
    ψ .*= exp.(im*r2_scaling(X,λ,λ̇,t)/2)
end

init_phase!(ψ) = remove_phase!(ψ,0)

## First, we verify Castin-Dum scaling inversion
# using parameters of Mewes et al 1996 (or smaller N)

# initialize default sim
L = (50,50,3)
N = (32,32,32)
sim = Sim(L,N)
@unpack_Sim sim;

# convenient units
mum = 1e-6
nm = 1e-9
amu = 1.661e-27
ħ = 1.055e-34
a0 = 0.05nm

# Rb87 parameters
m = 87amu
as = 5.8nm

# simulation params in units of wr, ar
wr = 18*2pi
wz = 320*2pi
az = sqrt(ħ/m/wz)
ar = sqrt(ħ/m/wr)
gdim = 4*pi*ħ^2*as/m
g = gdim/(ħ*wr)/ar^3

ar = 1.
az = wr/wz |> sqrt
wz /= wr
wr = 1.

γ = 0.1
μ = 50.0
Rr_tf = sqrt(2*μ)
Rz_tf = sqrt(2*μ/wz^2)

# make grids
X,K,dX,dK = makearrays(L,N)
x,y,z = X
yt = y'
zt = reshape(z,1,1,length(z))

# time info
tf = .2/γ
Nt = 100
t = LinRange(0.,tf,Nt)

# potential for initial state (prolate)
import FourierGPE:V
V(x,y,z,t) = 0.5*(x^2 + y^2 + wz^2*z^2)
λx(t) = 1.
λy(t) = 1.
λz(t) = 1.

λẋ(t) = 0.
λẏ(t) = 0.
λż(t) = 0.

λ(t) = (λx(t),λy(t),λz(t))
λ̇(t) = (λẋ(t),λẏ(t),λż(t))

# random initial state
# ψi = randn(N)+im*randn(N)

#Thomas-Fermi initial state
zmid = N[3]/2 |> Int
ψ0(x,y,z,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,z,0.0)/μ,0.0)+im*0.0)
ψi = ψ0.(x,yt,zt,μ,g)
showpsi(x,y,ψi[:,:,zmid])

ϕi = kspace(ψi,sim)

# check atom number
Na = sum(abs2.(ψi))*prod(dX)

#save files?
nfiles = true
path = joinpath(@__DIR__,"oblate_ground_state")

## Evolve in k space
@pack_Sim! sim
@time sol = runsim(sim)

## slice through z=0:
ϕg = sol[end]

# loadfile = joinpath(path,"save200.jld2")
# @load loadfile ψ
# ϕg = ψ
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg[:,:,zmid])

## one D plot
Vplot(x,y,z,t) = 0.5*(x^2 + y^2 + wz^2*z^2)
ψtfplot(x,y,z,μ,g)=sqrt(μ/g)*sqrt(max(1.0-Vplot(x,y,z,0.0)/μ,0.0)+im*0.0)
Plots.plot(x,g*abs2.(ψi[:,16,16]),label=L"|\psi(x)|^2")
Plots.plot!(x,g*abs2.(ψtfplot.(x,0,0,μ,g)),label=L"n_{TF}(x)")

## file for easy load
function loadpsi(i::Int,sim)
    padi = lpad(string(i),ndigits(length(sim.t)),"0")
    fromfile = joinpath(sim.path,sim.filename*padi*".jld2")
    ϕt = load(fromfile,"ψ")
end

## initial state for expansion
# tind = 200
# ϕt = sol[tind]
# ψt = xspace(ϕt,sim)

tind = 100
ϕt = loadpsi(tind,sim)
ψt = xspace(ϕt,sim)

@manipulate for i in 1:zmid
    showpsi(x,y,ψt[:,:,i])
end

## save initial state
# using JLD2, FileIO

tofile = "initial_state_oblate.jld2"
save(tofile,"ϕt",ϕt)

## Free expansion
fromfile = tofile
load(fromfile,"ϕt")

#save files?
nfiles = true
path = joinpath(@__DIR__,"oblate_expansion")

# potential
V(x,y,z,t) = 0.0

# Scaling GPE: free expansion scaling law (Castin-Dum)
λx(t) = sqrt(1+(wr*t)^2)
λy(t) = sqrt(1+(wr*t)^2)
λz(t) = sqrt(1+(wz*t)^2)

λẋ(t) = wr^2*t/λx(t)
λẏ(t) = wr^2*t/λy(t)
λż(t) = wz^2*t/λz(t)

λ(t) = (λx(t),λy(t),λz(t))
λ̇(t) = (λẋ(t),λẏ(t),λż(t))

# Parameters
γ = 0.0
ti = 0.0
tf = 10
Nt = 200
t = LinRange(ti,tf,Nt)
ϕi .= ϕt

ψi .= xspace(ϕi,sim)
init_phase!(ψi)
ϕi .= kspace(ψi,sim)

simes = Sim(L,N)
@pack_Sim! simes

simes
## Evolve
@time soles = runsim(simes)

## Plot
@manipulate for i in 1:zmid, tind in 1:length(t)
    ϕf = loadpsi(tind,simes)
    ψf = xspace(ϕf,simes)
    s = t[tind]
    # showpsi(x*λx(s),y*λy(s),ψf[:,:,i])
    # showpsi(x*λy(s),y*λz(s),ψf[i,:,:])
    p1 = Plots.heatmap(y*λy(s),z*λz(s),abs2.(ψf[i,:,:]),c=:inferno,colorbar=true)
    Plots.xlabel!(L"y/a_z")
    Plots.ylabel!(L"z/a_z")
end

## aspect ratio
aspect = zero.(t)

x = X[1]
z = X[3]

for i in eachindex(t)
    ϕf = loadpsi(i,simes)
    ψf = xspace(ϕf,simes)
    s = t[i]
    zt = reshape(z,1,1,length(z))
    sx = sum(abs2.(ψf).*x.^2) |> sqrt
    sz = sum(abs2.(ψf).*zt.^2) |> sqrt
    aspect[i] = (λz(s)/λx(s))*sz/sx
end

# aspect_tinf = 2/π/ϵ
# aspect_analytic(t) = (Rr_tf/Rz_tf)*λx(t)/λz(t)

p1 = Plots.plot(t,aspect,size=(500,250),legend=:bottomright,label="sGPE")
# Plots.plot!(t,aspect_analytic.(t),c=:red,label="Castin-Dum")
# Plots.ylims!(1e-1,1e1)
# Plots.plot!(t,one.(t)*aspect_tinf,ls=:dash,label=L"2\omega_\perp/(\pi\omega_z)")
Plots.xlabel!(L"\omega_z t")
Plots.ylabel!(L"R_z(t)/R_\perp(t)")

# <x^2>/<z^2> and Rx^2/Rz^2 are equivalent for Thomas-Fermi approximation

# savefig(p1,"expand_oblate_comparison.pdf")
## make movie
# @time p = densityfilm(soles,Nt,file="3dfree_expansion.gif")

## plot column density
@manipulate for i in 1:4:lastindex(t)
    ϕf = soles[i]
    ψf = xspace(ϕf,simes)
    ψcd = dropdims(sum(ψf,dims=3),dims=3)
    s = t[i]
    showpsi(x*λ(s),y*λ(s),ψcd/λ(s))
end

## compute aspect ratio from standard deviations in r,z
dx = diff(x)[1]; dy = diff(y)[1];dz = diff(z)[1]
Na(ψ) = sum(abs2.(ψ))*dx*dy*dz
Rz(ψ) = sum(abs2.(ψ).*zt.^2)*dx*dy*dz/Na(ψ) |> abs |> sqrt
Rr(ψ) = sum(abs2.(ψ).*(x.^2 .+ yt.^2))*dx*dy*dz/Na(ψ) |> abs |> sqrt
aspect_ratio(ψ) = Rr(ψ)/Rz(ψ)

α0 = aspect_ratio(xspace(soles[1],simes))
α = zero(t)

for j in eachindex(t)
    α[j] = aspect_ratio(xspace(soles[j],simes))
end

plot(t,α,legend=false)
ylims!(0,3.2)
plot!(t,one.(t)*2/α0/pi,ls=:dash)

## check against TF radii with fitting
Rz_tf0 = sqrt(2*μ)
Rperp_tf0 = sqrt(2*μ/64)

a0 = Rperp_tf0/Rz_tf0

2/pi/a0


text:
using SLEEF
using BenchmarkTools
using JLD, DataStructures
using Printf

const RETUNE  = false
const VERBOSE = true
const DETAILS = false
const test_types = (Float64, Float32) # Which types do you want to bench?
const bench = ("Base", "SLEEF")
const suite = BenchmarkGroup()
for n in bench
    suite[n] = BenchmarkGroup([n])
end

bench_reduce(f::Function, X) = mapreduce(x -> reinterpret(Unsigned,x), |, f(x) for x in X)

using Base.Math.IEEEFloat

MRANGE(::Type{Float64}) = 10000000
MRANGE(::Type{Float32}) = 10000
IntF(::Type{Float64}) = Int64
IntF(::Type{Float32}) = Int32
x_trig(::Type{T}) where {T<:IEEEFloat} = begin
    x_trig = T[]
    for i = 1:10000
        s = reinterpret(T, reinterpret(IntF(T), T(pi)/4 * i) - IntF(T)(20))
        e = reinterpret(T, reinterpret(IntF(T), T(pi)/4 * i) + IntF(T)(20))
        d = s
        while d <= e
            append!(x_trig, d)
            d = reinterpret(T, reinterpret(IntF(T), d) + IntF(T)(1))
        end
    end
    x_trig = append!(x_trig, -10:0.0002:10)
    x_trig = append!(x_trig, -MRANGE(T):200.1:MRANGE(T))
end
x_exp(::Type{T}) where {T<:IEEEFloat}        = map(T, vcat(-10:0.0002:10, -1000:0.1:1000))
x_exp2(::Type{T}) where {T<:IEEEFloat}       = map(T, vcat(-10:0.0002:10, -120:0.023:1000, -1000:0.02:2000))
x_exp10(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(-10:0.0002:10, -35:0.023:1000, -300:0.01:300))
x_expm1(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(-10:0.0002:10, -1000:0.021:1000, -1000:0.023:1000, 10.0.^-(0:0.02:300), -10.0.^-(0:0.02:300), 10.0.^(0:0.021:300), -10.0.^-(0:0.021:300)))
x_log(::Type{T}) where {T<:IEEEFloat}        = map(T, vcat(0.0001:0.0001:10, 0.001:0.1:10000, 1.1.^(-1000:1000), 2.1.^(-1000:1000)))
x_log10(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(0.0001:0.0001:10, 0.0001:0.1:10000))
x_log1p(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(0.0001:0.0001:10, 0.0001:0.1:10000, 10.0.^-(0:0.02:300), -10.0.^-(0:0.02:300)))
x_atrig(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(-1:0.00002:1))
x_atan(::Type{T}) where {T<:IEEEFloat}       = map(T, vcat(-10:0.0002:10, -10000:0.2:10000, -10000:0.201:10000))
x_cbrt(::Type{T}) where {T<:IEEEFloat}       = map(T, vcat(-10000:0.2:10000, 1.1.^(-1000:1000), 2.1.^(-1000:1000)))
x_trigh(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(-10:0.0002:10, -1000:0.02:1000))
x_asinhatanh(::Type{T}) where {T<:IEEEFloat} = map(T, vcat(-10:0.0002:10, -1000:0.02:1000))
x_acosh(::Type{T}) where {T<:IEEEFloat}      = map(T, vcat(1:0.0002:10, 1:0.02:1000))
x_pow(::Type{T}) where {T<:IEEEFloat} = begin
    xx1 = map(Tuple{T,T}, [(x,y) for x = -100:0.20:100, y = 0.1:0.20:100])[:]
    xx2 = map(Tuple{T,T}, [(x,y) for x = -100:0.21:100, y = 0.1:0.22:100])[:]
    xx3 = map(Tuple{T,T}, [(x,y) for x = 2.1, y = -1000:0.1:1000])
    xx = vcat(xx1, xx2, xx2)
end

import Base.atanh
for f in (:atanh,)
    @eval begin
        ($f)(x::Float64) = ccall($(string(f)), Float64, (Float64,), x)
        ($f)(x::Float32) = ccall($(string(f,"f")), Float32, (Float32,), x)
    end
end

const micros = OrderedDict(
    "sin"   => x_trig,
    "cos"   => x_trig,
    "tan"   => x_trig,
    "asin"  => x_atrig,
    "acos"  => x_atrig,
    "atan"  => x_atan,
    "exp"   => x_exp,
    "exp2"  => x_exp2,
    "exp10" => x_exp10,
    "expm1" => x_expm1,
    "log"   => x_log,
    "log2"  => x_log10,
    "log10" => x_log10,
    "log1p" => x_log1p,
    "sinh"  => x_trigh,
    "cosh"  => x_trigh,
    "tanh"  => x_trigh,
    "asinh" => x_asinhatanh,
    "acosh" => x_acosh,
    "atanh" => x_asinhatanh,
    "cbrt"  => x_cbrt
    )

for n in bench
    for (f,x) in micros
        suite[n][f] = BenchmarkGroup([f])
        for T in test_types
            fex = Expr(:., Symbol(n), QuoteNode(Symbol(f)))
            suite[n][f][string(T)] = @benchmarkable bench_reduce($fex, $(x(T)))
        end
    end
end

tune_params = joinpath(@__DIR__, "params.jld")
if !isfile(tune_params) || RETUNE
    tune!(suite; verbose=VERBOSE, seconds = 2)
    save(tune_params, "suite", params(suite))
    println("Saving tuned parameters.")
else
    println("Loading pretuned parameters.")
    loadparams!(suite, load(tune_params, "suite"), :evals, :samples)
end

println("Running micro benchmarks...")
results = run(suite; verbose=VERBOSE, seconds = 2)

printstyled("Benchmarks: median ratio SLEEF/Base\n", color = :blue)
for f in keys(micros)
    printstyled(string(f) color = :magenta)
    for T in test_types
        println()
        print("time: ", )
        tratio = ratio(median(results["SLEEF"][f][string(T)]), median(results["Base"][f][string(T)])).time
        tcolor = tratio > 3 ? :red : tratio < 1.5 ? :green : :blue
        printstyled(@sprintf("%.2f",tratio), " ", string(T), color = tcolor)
        if DETAILS
            printstyled("details SLEEF/Base\n", color=:blue)
            println(results["SLEEF"][f][string(T)])
            println(results["Base"][f][string(T)])
            println()
        end
    end
    println("\n")
end


text:
function remove_tips!(sdg::Graphs.SDG, min_size::Integer)
    @info "Beginning tip removal process"
    #sdg = graph(ws)
    pass = 1
    tips = Graphs.find_tip_nodes(sdg, min_size)
    ntips = length(tips)
    utgs = Vector{Graphs.SequenceDistanceGraphPath{typeof(sdg)}}()
    newnodes = Vector{Graphs.NodeID}()
    while true
        @info string("Pass number: ", pass)
        if isempty(tips)
            @info "No more tips in graph"
            break
        end
        @info string("Found ", ntips, " tips")
        for tip in tips
            Graphs.remove_node!(sdg, tip)
        end
        @info string("Removed ", ntips, " tips")
        @info "Collapsing any resulting transient paths"
        Graphs.collapse_all_unitigs!(utgs, newnodes, sdg, 2, true)
        pass = pass + 1
        tips = Graphs.find_tip_nodes(sdg, min_size)
        ntips = length(tips)
    end
    @info string("Finished tip removal process in ", pass, " passes")
    #return ws
    return sdg
end

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
using GeoJSON, FactCheck, Turf

point1 = GeoJSON.Feature(GeoJSON.Point(Float64[-75.343, 39.984]))
point2 = GeoJSON.Feature(GeoJSON.Point(Float64[-75.534, 39.123]))

@fact Turf.distance(point1, point2, :miles) --> roughly(60.37218405837491)
@fact Turf.distance(point1, point2, :kilometers) --> roughly(97.15957803131901)
@fact Turf.distance(point1, point2, :radians) --> roughly(0.015245501024842149)
@fact Turf.distance(point1, point2, :degrees) --> roughly(0.8735028650863799)
@fact Turf.distance(point1, point2) --> distance(point1, point2, :kilometers) # default

polygon = GeoJSON.Polygon(Vector{Vector{Float64}}[Vector{Float64}[[0,0],[10,0],[10,10],[0,10],[0,0]]])
@fact Turf.area(polygon) --> roughly(1232921098571.292)
point = GeoJSON.Point(Float64[0, 0])
@fact Turf.area(point) --> roughly(0)
line = GeoJSON.LineString(Vector{Float64}[[0,0],[1,2]])
@fact Turf.area(line) --> roughly(0)
polygonFC = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(polygon, Dict{AbstractString,Any}())])
@fact Turf.area(polygonFC) --> roughly(1232921098571.292)
polygonFC = GeoJSON.FeatureCollection(GeoJSON.Feature[GeoJSON.Feature(polygon, Dict{AbstractString,Any}()),
                                                      GeoJSON.Feature(polygon, Dict{AbstractString,Any}())])
@fact Turf.area(polygonFC) --> roughly(1232921098571.292 * 2)
feature = GeoJSON.Feature(polygon, Dict{AbstractString,Any}())
@fact Turf.area(feature) --> roughly(1232921098571.292)

pt = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/pt.geojson"))
@fact typeof(pt) --> GeoJSON.Feature
pts = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/pts.geojson"))
@fact typeof(pts) --> GeoJSON.FeatureCollection
closest = Turf.nearest(pt, pts)
@fact typeof(closest) --> GeoJSON.Feature
@fact typeof(closest.geometry) --> GeoJSON.Point
@fact closest.geometry.coordinates --> roughly([-75.33, 39.44])

polygon = Turf.bboxpolygon([0.0,0.0,10.0,10.0])
@fact typeof(polygon) --> GeoJSON.Polygon
@fact length(polygon.coordinates) --> 5
@fact polygon.coordinates[1] --> roughly(polygon.coordinates[end])

fc = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/fc.geojson"))
@fact typeof(fc) --> GeoJSON.FeatureCollection
bbox = Turf.envelope(fc)
@fact typeof(bbox) --> GeoJSON.Polygon

pt = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/Point.geojson"))
line = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/LineString.geojson"))
poly = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/Polygon.geojson"))
multiline = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/MultiLineString.geojson"))
multipoly = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/MultiPolygon.geojson"))
@fact Turf.extent(fc) --> roughly([20, -10, 130, 4])
@fact Turf.extent(pt) --> roughly([102, 0.5, 102, 0.5])
@fact Turf.extent(line) --> roughly([102, -10, 130, 4])
@fact Turf.extent(poly) --> roughly([100, 0, 101, 1])
@fact Turf.extent(multiline) --> roughly([100, 0, 103, 3])
@fact Turf.extent(multipoly) --> roughly([100, 0, 103, 3])

@fact Turf.square(Float64[0,0,5,10]) --> roughly([-2.5,0,7.5,10])
@fact Turf.square(Float64[0,0,10,5]) --> roughly([0,-2.5,10,7.5])

@fact Turf.size(Float64[0,0,10,10], 2) --> roughly([-5,-5,15,15])
@fact Turf.size(Float64[0,0,4,4], 1) --> roughly([0,0,4,4])
@fact Turf.size(Float64[0,0,4,4], 2) --> roughly([-2,-2,6,6])
@fact Turf.size(Float64[0,0,4,4], 0.5) --> roughly([1,1,3,3])
@fact Turf.size(Float64[-10,-10,0,0], 2) --> roughly([-15,-15,5,5])
@fact Turf.size(Float64[0,0,10,10], 1.5) --> roughly([-2.5,-2.5,12.5,12.5])
@fact Turf.size(Float64[0,0,10,10], 0.5) --> roughly([2.5,2.5,7.5,7.5])

feature = Turf.center(fc)
@fact typeof(feature) --> GeoJSON.Feature
@fact typeof(feature.geometry) --> GeoJSON.Point
@fact feature.geometry.coordinates --> roughly([75, -3])

pt1 = GeoJSON.Point(Float64[0, 0])
pt2 = GeoJSON.Point(Float64[10, 0])
midpt = Turf.midpoint(pt1, pt2)
@fact typeof(midpt) --> GeoJSON.Point
@fact midpt.coordinates --> roughly([0, 5])
pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[1, 1]))
pt2 = GeoJSON.Feature(GeoJSON.Point(Float64[11, 11]))
@fact typeof(midpt) --> GeoJSON.Feature
@fact midpt.geometry.coordinates --> roughly([6, 6])

pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[-75.4, 39.4]))
pt2 = GeoJSON.Feature(GeoJSON.Point(Float64[-75.534, 39.123]))
bearing = Turf.bearing(pt1, pt2) # what's the expected value?

pt1 = GeoJSON.Feature(GeoJSON.Point(Float64[-75.0, 39.0]))
dist = 100.0
bearing = 180.0
pt2 = Turf.destination(pt1, dist, bearing) # what's the expected value?

route1 = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/route1.geojson"))
route2 = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/route2.geojson"))
Turf.linedistance(route1, :miles)
# t.equal(Math.round(lineDistance(route1, 'miles')), 202);
#     t.true((lineDistance(route2, 'kilometers') - 742) < 1 &&
#             (lineDistance(route2, 'kilometers') - 742) > (-1) );

line = GeoJSON.parsefile(joinpath(dirname(@__FILE__),"geojson/dc-line.geojson"))
pt1 = GeoJSON.Feature(Turf.along(line, 1, :miles))
pt2 = GeoJSON.Feature(Turf.along(line, 1.2, :miles))
pt3 = GeoJSON.Feature(Turf.along(line, 1.4, :miles))
pt4 = GeoJSON.Feature(Turf.along(line, 1.6, :miles))
pt5 = GeoJSON.Feature(Turf.along(line, 1.8, :miles))
pt6 = GeoJSON.Feature(Turf.along(line, 2, :miles))
pt7 = GeoJSON.Feature(Turf.along(line, 100, :miles))
pt8 = GeoJSON.Feature(Turf.along(line, 0, :miles))
fc = GeoJSON.FeatureCollection(GeoJSON.Feature[pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8])
for feature in fc.features
    @fact typeof(feature.geometry) --> GeoJSON.Point
end
@fact length(fc.features) --> 8
@fact fc.features[end].geometry.coordinates --> roughly(pt8.geometry.coordinates)


text:
using Lifting
using Test

@testset "Maths" begin
    @test isapprox(calcRPE(3, calcIntensity(3, 8)), 8)
    @test isapprox(calcReps(calcIntensity(3, 8), calcRPE(3, calcIntensity(3, 8))), 3)

    @test calcIntensity(1, 3) == 0.7927699381639448
    @test calcIntensity(1, 5) == 0.8369601606963508
    @test calcIntensity(1, 6) == 0.8609556607834696
    @test calcIntensity(2, 1) == 0.6837373081262179
    @test calcIntensity(2, 10) == 0.9320967516428205
    @test calcIntensity(3, 5) == 0.7674204440806303
    @test calcIntensity(3, 8) == 0.8416327675690841
    @test calcIntensity(4, 1) == 0.5781098697807519
    @test calcIntensity(4, 7) == 0.7856495497947491
    @test calcIntensity(4, 8) == 0.810257864565398
    @test calcIntensity(5, 1) == 0.536711034778875
    @test calcIntensity(5, 2) == 0.6135345726731701
    @test calcIntensity(5, 4) == 0.6833868653044488
    @test calcIntensity(5, 5) == 0.7092198581560283
    @test calcIntensity(6, 3) == 0.6263113393667992
    @test calcIntensity(6, 5) == 0.6833479493866953
    @test calcIntensity(7, 5) == 0.6593075387107712
    @test calcIntensity(7, 10) == 0.7750478868872971
    @test calcIntensity(8, 2) == 0.5319184303087121
    @test calcIntensity(8, 5) == 0.6369071787400383
    @test calcIntensity(9, 2) == 0.5093407432979247
    @test calcIntensity(9, 9) == 0.704666457876605
    @test calcIntensity(9, 10) == 0.7262574340517902
    @test calcIntensity(10, 3) == 0.5366678294469638

    @test calcRPE(1, 0.8) == 3.342342342342342
    @test calcRPE(2, 0.3) == 0.05156764251678613
    @test calcRPE(2, 0.7) == 1.3099066269984132
    @test calcRPE(2, 0.9) == 8.88856743957622
    @test calcRPE(2, 1.0) == 10.0
    @test calcRPE(3, 0.9) == 10.0
    @test calcRPE(4, 0.8) == 7.586200791837665
    @test calcRPE(5, 0.4) == 0.3934674694286991
    @test calcRPE(6, 0.1) == 0.05901737061084574
    @test calcRPE(7, 0.4) == 0.6265524536794053
    @test calcRPE(7, 0.8) == 10.0
    @test calcRPE(7, 1.0) == 10.0
    @test calcRPE(9, 0.2) == 0.23682903555711834
    @test calcRPE(9, 0.6) == 4.380666149879444
    @test calcRPE(10, 1.0) == 10.0

    @test calcReps(0.2, 3) == 57
    @test calcReps(0.2, 5) == 72
    @test calcReps(0.2, 8) == 86
    @test calcReps(0.3, 5) == 41
    @test calcReps(0.4, 8) == 32
    @test calcReps(0.5, 4) == 14
    @test calcReps(0.6, 1) == 4
    @test calcReps(0.6, 6) == 11
    @test calcReps(0.7, 1) == 2
    @test calcReps(0.7, 8) == 8
    @test calcReps(0.8, 3) == 1
    @test calcReps(0.8, 6) == 3
    @test calcReps(0.8, 10) == 6
    @test calcReps(1.0, 9) == 0

    @test calcIntensityRatio(1, 2, 8, 6) == 0.8554603294340484
    @test calcIntensityRatio(1, 8, 10, 1) == 0.4327753512915275
    @test calcIntensityRatio(2, 2, 10, 3) == 0.7401991037647249
    @test calcIntensityRatio(2, 7, 10, 3) == 0.6316772019672708
    @test calcIntensityRatio(2, 9, 7, 3) == 0.6656842312192123
    @test calcIntensityRatio(3, 8, 7, 8) == 0.8661098904150482
    @test calcIntensityRatio(3, 10, 3, 10) == 1.0
    @test calcIntensityRatio(4, 4, 7, 7) == 0.9922434668550097
    @test calcIntensityRatio(4, 5, 4, 7) == 1.0657925380128117
    @test calcIntensityRatio(4, 6, 10, 7) == 0.842533486694681
    @test calcIntensityRatio(5, 7, 7, 8) == 0.9626050463922357
    @test calcIntensityRatio(5, 8, 4, 10) == 1.1034631381591704
    @test calcIntensityRatio(5, 10, 4, 1) == 0.6957552282811349
    @test calcIntensityRatio(6, 8, 8, 3) == 0.7664551245770161
    @test calcIntensityRatio(6, 10, 7, 2) == 0.694007341735247
    @test calcIntensityRatio(8, 4, 5, 1) == 0.8793406504937742
    @test calcIntensityRatio(8, 4, 6, 1) == 0.8205956008180643
    @test calcIntensityRatio(8, 8, 9, 2) == 0.7220860050141765
    @test calcIntensityRatio(9, 7, 8, 8) == 1.0658963533980068
    @test calcIntensityRatio(9, 8, 1, 8) == 1.3366720451385716

    @test calcRPERatio(1, 0.1, 10, 0.8) == 2.9993973518287506e14
    @test calcRPERatio(2, 0.4, 1, 0.4) == -4.614083979793281e-14
    @test calcRPERatio(2, 0.8, 10, 0.8) == 2.0069922576450354
    @test calcRPERatio(4, 0.6, 2, 0.4) == 0.07359916350276016
    @test calcRPERatio(4, 0.7, 3, 1.0) == 2.805871332615087
    @test calcRPERatio(4, 0.9, 5, 0.9) == 1.0
    @test calcRPERatio(4, 0.9, 8, 0.7) == 0.7755866011749226
    @test calcRPERatio(4, 0.9, 9, 0.7) == 0.8782128594994774
    @test calcRPERatio(5, 0.2, 6, 0.3) == 2.4186452245040173
    @test calcRPERatio(5, 0.5, 1, 0.5) == 8.813339622784159e-15
    @test calcRPERatio(5, 0.9, 6, 0.2) == 0.014388869488003317
    @test calcRPERatio(6, 0.2, 7, 0.6) == 20.620074858037235
    @test calcRPERatio(6, 0.4, 6, 1.0) == 19.73972406773003
    @test calcRPERatio(6, 0.5, 4, 0.7) == 3.5867973963742124
    @test calcRPERatio(6, 0.6, 9, 0.2) == 0.10156686550096347
    @test calcRPERatio(6, 1.0, 2, 0.7) == 0.13099066269984133
    @test calcRPERatio(6, 1.0, 3, 0.7) == 0.246570448138738
    @test calcRPERatio(8, 0.1, 9, 0.1) == 1.147351193973504
    @test calcRPERatio(8, 0.5, 8, 0.9) == 6.495419905922825
    @test calcRPERatio(8, 1.0, 3, 0.7) == 0.246570448138738
    @test calcRPERatio(9, 0.2, 5, 1.0) == 42.22455230827558
    @test calcRPERatio(9, 0.3, 5, 1.0) == 21.486366386831257

    @test calcRepRatio(0.1, 7, 0.1, 7) == 1.0
    @test calcRepRatio(0.1, 7, 1.0, 3) == 0.0
    @test calcRepRatio(0.1, 8, 0.2, 8) == 0.441025641025641
    @test calcRepRatio(0.2, 8, 0.7, 7) == 0.08139534883720931
    @test calcRepRatio(0.3, 5, 0.6, 9) == 0.36585365853658536
    @test calcRepRatio(0.3, 7, 0.3, 8) == 1.0638297872340425
    @test calcRepRatio(0.3, 9, 0.9, 2) == 0.0
    @test calcRepRatio(0.4, 1, 0.2, 5) == 7.2
    @test calcRepRatio(0.5, 2, 1.0, 7) == 0.0
    @test calcRepRatio(0.5, 3, 0.6, 5) == 0.8333333333333334
    @test calcRepRatio(0.6, 2, 0.9, 8) == 0.2
    @test calcRepRatio(0.6, 10, 0.8, 3) == 0.0625
    @test calcRepRatio(0.7, 1, 0.5, 6) == 9.0
    @test calcRepRatio(0.7, 6, 0.6, 6) == 1.8333333333333333
    @test calcRepRatio(0.8, 3, 0.9, 8) == 1.0
    @test calcRepRatio(0.8, 7, 0.1, 3) == 33.0
    @test calcRepRatio(0.8, 10, 0.1, 10) == 34.666666666666664
    @test calcRepRatio(0.9, 3, 0.5, 8) == Inf
    @test calcRepRatio(0.9, 7, 0.1, 9) == 202.0
    @test calcRepRatio(0.9, 8, 0.2, 9) == 90.0
    @test isnan(calcRepRatio(1.0, 8, 1.0, 3))

    weight = rand(10:30, 20)
    actualReps = rand(1:10, 20)
    actualRPE = rand(1:10, 20)
    targetReps = rand(1:10, 20)
    targetRPE = rand(1:10, 20)

    @test calcRepMax.(weight, actualReps, actualRPE, targetReps, targetRPE) ==
          weight .* calcIntensityRatio.(actualReps, actualRPE, targetReps, targetRPE)

    @test intensityArb(1) == 0.9724788485850433
    @test intensityArb(2) == 0.9419743782969102
    @test intensityArb(3) == 0.9133254178463787
    @test intensityArb(4) == 0.8863676653075695
    @test intensityArb(5) == 0.8609556607834696
    @test intensityArb(6) == 0.8369601606963508
    @test intensityArb(7) == 0.814265939255761
    @test intensityArb(8) == 0.7927699381639448
    @test intensityArb(9) == 0.7723797018614351
    @test intensityArb(10) == 0.753012048192771
    @test intensityArb(11) == 0.7345919341805627
    @test intensityArb(12) == 0.7170514842965725
    @test intensityArb(13) == 0.7003291547027102
    @test intensityArb(14) == 0.684369011771147
    @test intensityArb(15) == 0.6691201070592172
    @test intensityArb(16) == 0.6545359340227779
    @test intensityArb(17) == 0.6405739542630196
    @test intensityArb(18) == 0.6271951831409934
    @test intensityArb(19) == 0.61436382625791
    @test intensityArb(20) == 0.6020469596628537
end

max_stars_repo_name: dcelisgarza/Lifting


text:
# remove deps.jl if it exists, in case build.jl fails
isfile("deps.jl") && rm("deps.jl")

using Libdl

println("Pardiso library")
println("===============")

const LIBPARDISONAMES =
if Sys.iswindows()
[
    "libpardiso600-WIN-X86-64.dll",
    "libpardiso500-WIN-X86-64.dll",
]
elseif Sys.isapple()
[
    "libpardiso600-MACOS-X86-64.dylib",
    "libpardiso500-MACOS-X86-64.dylib",
]
elseif Sys.islinux()
[
    "libpardiso600-GNU800-X86-64.so",
    "libpardiso600-GNU720-X86-64.so",
    "libpardiso500-GNU461-X86-64.so",
    "libpardiso500-GNU472-X86-64.so",
    "libpardiso500-GNU481-X86-64.so",
]
else
    error("unhandled OS")
end

println("Looking for libraries with name: ", join(LIBPARDISONAMES, ", "), ".")

PATH_PREFIXES = [@__DIR__; get(ENV, "JULIA_PARDISO", [])]

if !haskey(ENV, "JULIA_PARDISO")
    println("INFO: use the `JULIA_PARDISO` environment variable to set a path to " *
            "the folder where the Pardiso library is located")
end

pardiso_version = 0
function find_paradisolib()
    found_lib = false
    for prefix in PATH_PREFIXES
        println("Looking in \"$(abspath(prefix))\" for libraries")
        for libname in LIBPARDISONAMES
            local path
            try
                path = joinpath(prefix, libname)
                if isfile(path)
                    println("    found \"$(abspath(path))\", attempting to load it...")
                    Libdl.dlopen(path, Libdl.RTLD_GLOBAL)
                    println("    loaded successfully!")
                    global PARDISO_LIB_FOUND = true
                    if occursin("600", libname)
                        global pardiso_version = 6
                    else
                        global pardiso_version = 5
                    end
                    return path, true
                end
            catch e
                println("    failed to load due to:")
                Base.showerror(stderr, e)
            end
        end
    end
    println("did not find libpardiso, assuming PARDISO 5/6 is not installed")
    return "", false
end

pardisopath, found_pardisolib = find_paradisolib()

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

println("\nMKL Pardiso")
println("=============")
function find_mklparadiso()
    if haskey(ENV, "MKLROOT")
        println("found MKLROOT environment variable, using it")
        return ENV["MKLROOT"], true
    end
    println("did not find MKLROOT environment variable, assuming MKL is not installed")
    return "", false
end

mklroot, found_mklpardiso = find_mklparadiso()

if !(found_mklpardiso || found_pardisolib)
    println("WARNING: no Pardiso library managed to load")
end

open("deps.jl", "w") do f
    print(f,
"""
const MKL_PARDISO_LIB_FOUND = $found_mklpardiso
const PARDISO_LIB_FOUND = $found_pardisolib
const PARDISO_VERSION = $pardiso_version
const MKLROOT = $(repr(mklroot))
const PARDISO_PATH = raw"$pardisopath"
"""
)

end


text:
# accuracies: (0.7895149180648964,0.7598746490827197,0.37677090814128095)
include("CenterOfMass.jl")
include("setup_parameters.jl")
include("plotting_tools.jl")
using PyCall
@pyimport skimage.measure as measure

plt[:close]("all")
pygui(false)

feasibility_constraint = 0.9

folder = "feasibility"

actuator_limit = 6.0
max_robots = 4
n_trial = 20
n_measurement = 20

mass_resolution = 0.05
masses = collect(min_mass:mass_resolution:max_mass)
resolution = 0.05

function analyze_feasibility(robots)
  # initialize feasibility histogram
  prior = initialize_prior(circle_ps, resolution, interior_q, masses)
  feasibility = Histogram(get_range(prior))

  all_robots = collect(1:length(attachment_ws))
  chosen_ws = attachment_ws[robots]
  remaining_ws = attachment_ws[setdiff(all_robots, robots)]

  ranges = get_range(feasibility)
  for (ii, x) = enumerate(ranges[1])
    for (jj, y) = enumerate(ranges[2])
      for (kk, m) = enumerate(ranges[3])
        Convex.clearmemory()
        gravity_w = g * m * point_to_wrench([x;y])

        feasible = check_feasible_configuration(chosen_ws, remaining_ws,
          gravity_w, actuator_limit, max_robots)
        get_data(feasibility)[ii,jj,kk] = feasible
      end
    end
  end

  f = figure()
  f[:add_subplot](111,projection="3d")

  fill3d(hcat(circle_ps...), alpha=0.2)
  plot_attachment_points(attachment_ps)
  plot_occupied_points(attachment_ps[robots])
  plot_quadrotor([0;0;0.025], color="k", scale=0.15, alpha=0.0)

  cloud = to_cloud(feasibility, interior_q)
  scaled_ps = (15*cloud[4,:]'/maximum(cloud[4,:]))
  #scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', "z", scaled_ps.^2, "gray", alpha = 0.5)
  verts, faces = measure.marching_cubes(get_data(feasibility), level=0.5,
  spacing=(resolution,resolution,mass_resolution))
  verts[:,[1,2]] -= 1.0
  verts[:,3] += min_mass
  plot_trisurf(verts[:,1],verts[:,2],faces,verts[:,3], cmap="BuPu", lw=0.5)

  xlabel("X")
  ylabel("Y")
  zlabel("M")

  xlim(-1.0,1.0)
  ylim(-1.0,1.0)
  zlim(0.0,1.405)

  if ~isdir(folder)
    mkdir(folder)
  end
  savefig("$(folder)/feasibility_$(robots).png", pad_inches=0.01, bbox_inches="tight")
  feasibility_probability = sum(get_data(prior).*get_data(feasibility))
  feasibility_probability
end

f1 = analyze_feasibility(Int64[])
f2 = analyze_feasibility([5,10])
f3 = analyze_feasibility([1,3,5,8])
pygui(true)
@show (f1,f2,f3)
Void

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
using Test
using TDD
@testset "graph_test" begin
    graph1 = [
          [2, 3],     # list of nodes that can be reached from node 1
          [1],        # list of nodes that can be reached from node 2
          [1],        # ", node 3
          [5],        # ", node 4
          [4]]        # ", node 5

    graph2 = [
          [],     # node 1
          [1,3],      # node 2
          [],        # ", node 3
          [],        # ", node 4
          [4],        # ", node 5
          [3,4]]        # ", node 6

    graph3 = Bool[
        1 1 1 0 0;
        1 1 0 0 0;
        1 0 1 0 0;
        0 0 0 0 1;
        0 0 0 1 1;
      ]
    nodes = neighbors(graph1, 1)
    @test nodes == [1,2,3]
    nodes = neighbors(graph2, 4)
    @test nodes == [4]
    nodes = neighbors(graph3, 4)
    @test nodes == [4,5]

    nodes = connectednodes(graph1, 1)
    @test nodes == [1,2,3]
    nodes = connectednodes(graph2, 5)
    @test nodes == [4,5]
    nodes = connectednodes(graph3, 4)
    @test nodes == [4,5]

    parts=connectedparts(graph1)
    @test parts==[[1,2,3], [4,5]]
    parts=connectedparts(graph2)
    @test parts==[[1,2,3,4,5,6]]
    nodes = connectedparts(graph3)
    @test nodes == [[1,2,3], [4,5]]

end;


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

println("Hello World")


text:
module InlineStrings

import Base: ==

using Parsers

export InlineString, InlineStringType, inlinestrings

"""
    InlineString

A set of custom string types of various fixed sizes. Each inline string
is a custom primitive type and can benefit from being stack friendly by
avoiding allocations/heap tracking in the GC. When used in an array,
the elements are able to be stored inline since each one has a fixed
size. Currently support inline strings from 1 byte up to 255 bytes.
See more details by looking at individual docs for [`String1`](@ref),
[`String3`](@ref), [`String7`](@ref), [`String15`](@ref), [`String31`](@ref),
[`String63`](@ref), [`String127`](@ref), or [`String255`](@ref).
"""
abstract type InlineString <: AbstractString end

for sz in (1, 4, 8, 16, 32, 64, 128, 256)
    nm = Symbol(:String, max(1, sz - 1))
    nma = Symbol(:InlineString, max(1, sz - 1))
    @eval begin
        """
            $($nm)(str::AbstractString)
            $($nm)(bytes::AbstractVector{UInt8}, pos, len)
            $($nm)(ptr::Ptr{UInt8}, [len])

        Custom fixed-size string with a fixed size of $($sz) bytes.
        1 byte is used to store the length of the string. If an
        inline string is shorter than $($(max(1, sz - 1))) bytes, the entire
        string still occupies the full $($sz) bytes since they are,
        by definition, fixed size. Otherwise, they can be treated
        just like normal `String` values. Note that `sizeof(x)` will
        return the # of _codeunits_ in an $($nm) like `String`, not
        the total fixed size. For the fixed size, call `sizeof($($nm))`.
        $($nm) can be constructed from an existing `String` (`$($nm)(x::AbstractString)`),
        from a byte buffer with position and length (`$($nm)(buf, pos, len)`),
        from a pointer with optional length (`$($nm)(ptr, len)`)
        or built iteratively by starting with `x = $($nm)()` and calling
        `x, overflowed = InlineStrings.addcodeunit(x, b::UInt8)` which returns a
        new $($nm) with the new codeunit `b` appended and an `overflowed` `Bool`
        value indicating whether too many codeunits have been appended for the
        fixed size. When constructed from a pointer, note that the `ptr` must
        point to valid memory or program data may become corrupt. If the `len`
        argument is specified with the pointer, it must fit within the fixed size
        of $($nm); if no length is provided, the C-string is assumed to be
        NUL-terminated. If the NUL-terminated string ends up longer than can
        fit in $($nm), an ArgumentError will be thrown.
        """
        primitive type $nm <: InlineString $(sz * 8) end
        const $nma = $nm
        export $nm
        export $nma
    end
end

_bswap(x::T) where {T <: InlineString} = T === InlineString1 ? x : Base.bswap_int(x)

const InlineStringTypes = Union{InlineString1,
                            InlineString3,
                            InlineString7,
                            InlineString15,
                            InlineString31,
                            InlineString63,
                            InlineString127,
                            InlineString255}

function Base.promote_rule(::Type{T}, ::Type{S}) where {T <: InlineString, S <: InlineString}
    T === InlineString1 && return S
    S === InlineString1 && return T
    T === InlineString3 && return S
    S === InlineString3 && return T
    T === InlineString7 && return S
    S === InlineString7 && return T
    T === InlineString15 && return S
    S === InlineString15 && return T
    T === InlineString31 && return S
    S === InlineString31 && return T
    T === InlineString63 && return S
    S === InlineString63 && return T
    T === InlineString127 && return S
    S === InlineString127 && return T
    return InlineString255
end

Base.promote_rule(::Type{T}, ::Type{String}) where {T <: InlineString} = String

Base.widen(::Type{InlineString1}) = InlineString3
Base.widen(::Type{InlineString3}) = InlineString7
Base.widen(::Type{InlineString7}) = InlineString15
Base.widen(::Type{InlineString15}) = InlineString31
Base.widen(::Type{InlineString31}) = InlineString63
Base.widen(::Type{InlineString63}) = InlineString127
Base.widen(::Type{InlineString127}) = InlineString255
Base.widen(::Type{InlineString255}) = String

Base.ncodeunits(::InlineString1) = 1
Base.ncodeunits(x::InlineString) = Int(Base.trunc_int(UInt8, x))
Base.codeunit(::InlineString) = UInt8

Base.@propagate_inbounds function Base.codeunit(x::T, i::Int) where {T <: InlineString}
    @boundscheck checkbounds(Bool, x, i) || throw(BoundsError(x, i))
    if T === InlineString1
        return Base.bitcast(UInt8, x)
    else
        return Base.trunc_int(UInt8, Base.lshr_int(x, 8 * (sizeof(T) - i)))
    end
end

function Base.String(x::T) where {T <: InlineString}
    len = ncodeunits(x)
    out = Base._string_n(len)
    if T === InlineString1
        GC.@preserve out unsafe_store!(pointer(out), codeunit(x, 1))
        return out
    end
    ref = Ref{T}(_bswap(x))
    GC.@preserve ref out begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        unsafe_copyto!(pointer(out), ptr, len)
    end
    return out
end

function Base.Symbol(x::T) where {T <: InlineString}
    ref = Ref{T}(_bswap(x))
    return ccall(:jl_symbol_n, Ref{Symbol},
        (Ref{T}, Int), ref, sizeof(x))
end

# add a codeunit to end of string method
function addcodeunit(x::T, b::UInt8) where {T <: InlineString}
    if T === InlineString1
        return Base.bitcast(InlineString1, b), false
    end
    len = Base.trunc_int(UInt8, x)
    sz = Base.trunc_int(UInt8, sizeof(T))
    shf = Base.zext_int(Int16, max(0x01, sz - len - 0x01)) << 3
    x = Base.or_int(x, Base.shl_int(Base.zext_int(T, b), shf))
    return Base.add_int(x, Base.zext_int(T, 0x01)), (len + 0x01) >= sz
end

# from String
InlineString1(byte::UInt8=0x00) = Base.bitcast(InlineString1, byte)

function InlineString1(x::AbstractString)
    sizeof(x) == 1 || stringtoolong(InlineString1, sizeof(x))
    return Base.bitcast(InlineString1, codeunit(x, 1))
end

function InlineString1(buf::AbstractVector{UInt8}, pos, len)
    len == 1 || stringtoolong(InlineString1, len)
    return Base.bitcast(InlineString1, buf[pos])
end

function InlineString1(ptr::Ptr{UInt8}, len=nothing)
    ptr == Ptr{UInt8}(0) && nullptr(InlineString1)
    if len === nothing
        y, _ = addcodeunit(InlineString1(), unsafe_load(ptr))
        unsafe_load(ptr, 2) === 0x00 || stringtoolong(InlineString1, 2)
        return y
    else
        len == 1 || stringtoolong(InlineString1, len)
        return Base.bitcast(InlineString1, unsafe_load(ptr))
    end
end

function InlineString1(x::S) where {S <: InlineString}
    sizeof(x) == 1 || stringtoolong(InlineString1, sizeof(x))
    return Base.bitcast(InlineString1, codeunit(x, 1))
end

for T in (:InlineString3, :InlineString7, :InlineString15, :InlineString31, :InlineString63, :InlineString127, :InlineString255)
    @eval $T() = Base.zext_int($T, 0x00)

    @eval function $T(x::AbstractString)
        if typeof(x) === String && sizeof($T) <= sizeof(UInt)
            len = sizeof(x)
            len < sizeof($T) || stringtoolong($T, len)
            y = GC.@preserve x unsafe_load(convert(Ptr{$T}, pointer(x)))
            sz = 8 * (sizeof($T) - len)
            return Base.or_int(Base.shl_int(Base.lshr_int(_bswap(y), sz), sz), Base.zext_int($T, UInt8(len)))
        else
            len = ncodeunits(x)
            len < sizeof($T) || stringtoolong($T, len)
            y = $T()
            for i = 1:len
                @inbounds y, _ = addcodeunit(y, codeunit(x, i))
            end
            return y
        end
    end

    @eval function $T(buf::AbstractVector{UInt8}, pos, len)
        blen = length(buf)
        blen < len && buftoosmall(len)
        len < sizeof($T) || stringtoolong($T, len)
        if (blen - pos + 1) < sizeof($T)
            # if our buffer isn't long enough to hold a full $T,
            # then we can't do our unsafe_load trick below because we'd be
            # unsafe_load-ing memory from beyond the end of buf
            # we need to build the InlineString byte-by-byte instead
            y = $T()
            for i = pos:(pos + len - 1)
                @inbounds y, _ = addcodeunit(y, buf[i])
            end
            return y
        else
            y = GC.@preserve buf unsafe_load(convert(Ptr{$T}, pointer(buf, pos)))
            sz = 8 * (sizeof($T) - len)
            return Base.or_int(Base.shl_int(Base.lshr_int(_bswap(y), sz), sz), Base.zext_int($T, UInt8(len)))
        end
    end

    @eval function $T(ptr::Ptr{UInt8}, len=nothing)
        ptr == Ptr{UInt8}(0) && nullptr($T)
        y = $T()
        if len === nothing
            i = 1
            while true
                b = unsafe_load(ptr, i)
                b == 0x00 && break
                @inbounds y, overflowed = addcodeunit(y, b)
                overflowed && stringtoolong($T, i)
                i += 1
            end
        else
            len < sizeof($T) || stringtoolong($T, len)
            for i = 1:len
                @inbounds y, _ = addcodeunit(y, unsafe_load(ptr, i))
            end
        end
        return y
    end

    # between InlineStringTypes
    @eval function $T(x::S) where {S <: InlineString}
        if $T === S
            return x
        elseif sizeof($T) < sizeof(S)
            # trying to compress
            len = sizeof(x)
            len > (sizeof($T) - 1) && stringtoolong($T, len)
            y = Base.trunc_int($T, Base.lshr_int(x, 8 * (sizeof(S) - sizeof($T))))
            return Base.add_int(y, Base.zext_int($T, UInt8(len)))
        else
            # promoting smaller InlineString to larger
            if S === InlineString1
                y = Base.shl_int(Base.zext_int($T, x), 8 * (sizeof($T) - sizeof(S)))
            else
                y = Base.shl_int(Base.zext_int($T, Base.lshr_int(x, 8)), 8 * (sizeof($T) - sizeof(S) + 1))
            end
            return Base.add_int(y, Base.zext_int($T, UInt8(sizeof(x))))
        end
    end
end

@noinline nullptr(T) = throw(ArgumentError("cannot convert NULL to $T"))
@noinline buftoosmall(n) = throw(ArgumentError("input buffer too short for requested length: $n"))
@noinline stringtoolong(T, n) = throw(ArgumentError("string too large ($n) to convert to $T"))

function InlineStringType(n::Integer)
    n > 255 && stringtoolong(InlineString, n)
    return n == 1  ? InlineString1   : n < 4  ? InlineString3  :
           n < 8   ? InlineString7   : n < 16 ? InlineString15 :
           n < 32  ? InlineString31  : n < 64 ? InlineString63 :
           n < 128 ? InlineString127 : InlineString255
end

InlineString(x::InlineString) = x
InlineString(x::AbstractString)::InlineStringTypes = (InlineStringType(ncodeunits(x)))(x)

(==)(x::T, y::T) where {T <: InlineString} = Base.eq_int(x, y)
function ==(x::String, y::T) where {T <: InlineString}
    sizeof(x) == sizeof(y) || return false
    ref = Ref{T}(_bswap(y))
    return ccall(:memcmp, Cint, (Ptr{UInt8}, Ref{T}, Csize_t),
            pointer(x), ref, sizeof(x)) == 0
end
==(y::InlineString, x::String) = x == y

Base.cmp(a::T, b::T) where {T <: InlineString} =
    Base.eq_int(a, b) ? 0 : Base.ult_int(a, b) ? -1 : 1

function Base.hash(x::T, h::UInt) where {T <: InlineString}
    h += Base.memhash_seed
    ref = Ref{T}(_bswap(x))
    return ccall(Base.memhash, UInt,
        (Ref{T}, Csize_t, UInt32),
        ref, sizeof(x), h % UInt32) + h
end

function Base.write(io::IO, x::T) where {T <: InlineString}
    ref = Ref{T}(x)
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        Int(unsafe_write(io, ptr, reinterpret(UInt, sizeof(T))))::Int
    end
end

function Base.read(s::IO, ::Type{T}) where {T <: InlineString}
    return read!(s, Ref{T}())[]::T
end

function Base.print(io::IO, x::T) where {T <: InlineString}
    x isa InlineString1 && return print(io, Char(Base.bitcast(UInt8, x)))
    ref = Ref{T}(_bswap(x))
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        unsafe_write(io, ptr, sizeof(x))
        return
    end
end

function Base.isascii(x::T) where {T <: InlineString}
    if T === InlineString1
        return codeunit(x, 1) < 0x80
    end
    len = ncodeunits(x)
    x = Base.lshr_int(x, 8 * (sizeof(T) - len))
    for _ = 1:(len >> 2)
        y = Base.trunc_int(UInt32, x)
        (y & 0xff000000) >= 0x80000000 && return false
        (y & 0x00ff0000) >= 0x00800000 && return false
        (y & 0x0000ff00) >= 0x00008000 && return false
        (y & 0x000000ff) >= 0x00000080 && return false
        x = Base.lshr_int(x, 32)
    end
    return true
end

Base.chop(s::InlineString1; kw...) = chop(String3(s); kw...)
function Base.chop(s::InlineString; head::Integer = 0, tail::Integer = 1)
    if isempty(s)
        return s
    end
    n = ncodeunits(s)
    i = min(n + 1, max(nextind(s, firstindex(s), head), 1))
    j = max(0, min(n, prevind(s, lastindex(s), tail)))
    newlen = max(0, n - ((i - 1) + (n - j)))
    s = Base.shl_int(Base.lshr_int(s, 8), 8)
    return Base.or_int(Base.shl_int(s, (i - 1) * 8), _oftype(typeof(s), newlen))
end

Base.chomp(s::InlineString1) = chomp(String3(s))
function Base.chomp(s::InlineString)
    i = lastindex(s)
    if i < 1 || codeunit(s, i) != 0x0a
        return s
    elseif i < 2 || codeunit(s, i - 1) != 0x0d
        return Base.sub_int(s, _oftype(typeof(s), 1))
    else
        return Base.sub_int(s, _oftype(typeof(s), 2))
    end
end

Base.first(s::InlineString1, n::Integer) = first(String3(s), n)
function Base.first(s::InlineString, n::Integer)
    i = min(lastindex(s), nextind(s, 0, n))
    newlen = nextind(s, i) - 1
    s = Base.shl_int(Base.lshr_int(s, 8), 8)
    return Base.or_int(s, _oftype(typeof(s), newlen))
end

Base.last(s::InlineString1, n::Integer) = last(String3(s), n)
function Base.last(s::InlineString, n::Integer)
    nc = ncodeunits(s) + 1
    i = max(1, prevind(s, nc, n))
    i == 1 && return s
    newlen = nc - i
    return Base.or_int(Base.shl_int(s, (i - 1) * 8), _oftype(typeof(s), newlen))
end

Base.reverse(x::String1) = x
function Base.reverse(s::T) where {T <: InlineString}
    nc = ncodeunits(s)
    if isascii(s)
        len = Base.zext_int(T, Base.trunc_int(UInt8, s))
        x = Base.or_int(Base.shl_int(_bswap(s), 8 * (sizeof(T) - nc)), len)
        return x
    end
    x = Base.zext_int(T, Base.trunc_int(UInt8, s))
    i = 1
    while i <= nc
        j = nextind(s, i)
        _x = Base.lshr_int(s, 8 * (sizeof(T) - (j - 1)))
        n = j - i
        _x = Base.and_int(_x, n == 1 ? Base.zext_int(T, 0xff) :
            n == 2 ? Base.zext_int(T, 0xffff) :
            n == 3 ? Base.zext_int(T, 0xffffff) :
                     Base.zext_int(T, 0xffffffff))
        _x = Base.shl_int(_x, 8 * (sizeof(T) - (nc - (i - 1))))
        x = Base.or_int(x, _x)
        i = j
    end
    return x
end

@inline function Base.__unsafe_string!(out, x::T, offs::Integer) where {T <: InlineString}
    n = sizeof(x)
    ref = Ref{T}(_bswap(x))
    GC.@preserve ref out begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        unsafe_copyto!(pointer(out, offs), ptr, n)
    end
    return n
end

const BaseStrs = Union{Char, String, SubString{String}}
Base.string(a::InlineString) = String(a)
Base.string(a::InlineString...) = _string(a...)
Base.string(a::BaseStrs, b::InlineString) = _string(a, b)
Base.string(a::BaseStrs, b::BaseStrs, c::InlineString) = _string(a, b, c)
@inline function _string(a::Union{BaseStrs, InlineString}...)
    n = 0
    for v in a
        if v isa Char
            n += ncodeunits(v)
        else
            n += sizeof(v)
        end
    end
    out = Base._string_n(n)
    offs = 1
    for v in a
        offs += Base.__unsafe_string!(out, v, offs)
    end
    return out
end

function Base.repeat(x::T, r::Integer) where {T <: InlineString}
    r < 0 && throw(ArgumentError("can't repeat a string $r times"))
    r == 0 && return ""
    r == 1 && return s
    n = sizeof(x)
    out = Base._string_n(n * r)
    if n == 1 # common case: repeating a single-byte string
        @inbounds b = codeunit(x, 1)
        ccall(:memset, Ptr{Cvoid}, (Ptr{UInt8}, Cint, Csize_t), out, b, r)
    else
        for i = 0:r-1
            ref = Ref{T}(_bswap(x))
            GC.@preserve ref out begin
                ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
                unsafe_copyto!(pointer(out, i * n + 1), ptr, n)
            end
        end
    end
    return out
end

# copy/pasted from strings/util.jl
function Base.startswith(a::T, b::Union{String, SubString{String}, InlineString}) where {T <: InlineString}
    cub = ncodeunits(b)
    ncodeunits(a) < cub && return false
    ref = Ref{T}(_bswap(a))
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        if Base._memcmp(ptr, b, sizeof(b)) == 0
            nextind(a, cub) == cub + 1
        else
            false
        end
    end
end

function Base.endswith(a::T, b::Union{String, SubString{String}, InlineString}) where {T <: InlineString}
    cub = ncodeunits(b)
    astart = ncodeunits(a) - ncodeunits(b) + 1
    astart < 1 && return false
    ref = Ref{T}(_bswap(a))
    return GC.@preserve ref begin
        ptr = convert(Ptr{UInt8}, Base.unsafe_convert(Ptr{T}, ref))
        if Base._memcmp(ptr + (astart - 1), b, sizeof(b)) == 0
            thisind(a, astart) == astart
        else
            false
        end
    end
end

Base.match(r::Regex, s::InlineString, i::Integer) = match(r, String(s), i)
Base.findnext(r::Regex, s::InlineString, i::Integer) = findnext(r, String(s), i)

# the rest of these methods are copy/pasted from Base strings/string.jl file
# for efficiency
Base.@propagate_inbounds function Base.isvalid(x::InlineString, i::Int)
    @boundscheck checkbounds(Bool, x, i) || throw(BoundsError(x, i))
    return @inbounds thisind(x, i) == i
end

Base.@propagate_inbounds function Base.thisind(s::InlineString, i::Int)
    i == 0 && return 0
    n = ncodeunits(s)
    i == n + 1 && return i
    @boundscheck Base.between(i, 1, n) || throw(BoundsError(s, i))
    @inbounds b = codeunit(s, i)
    (b & 0xc0 == 0x80) & (i-1 > 0) || return i
    @inbounds b = codeunit(s, i-1)
    Base.between(b, 0b11000000, 0b11110111) && return i-1
    (b & 0xc0 == 0x80) & (i-2 > 0) || return i
    @inbounds b = codeunit(s, i-2)
    Base.between(b, 0b11100000, 0b11110111) && return i-2
    (b & 0xc0 == 0x80) & (i-3 > 0) || return i
    @inbounds b = codeunit(s, i-3)
    Base.between(b, 0b11110000, 0b11110111) && return i-3
    return i
end

Base.@propagate_inbounds function Base.nextind(s::InlineString, i::Int)
    i == 0 && return 1
    n = ncodeunits(s)
    @boundscheck Base.between(i, 1, n) || throw(BoundsError(s, i))
    @inbounds l = codeunit(s, i)
    (l < 0x80) | (0xf8 ≤ l) && return i+1
    if l < 0xc0
        i′ = @inbounds thisind(s, i)
        return i′ < i ? @inbounds(nextind(s, i′)) : i+1
    end
    # first continuation byte
    (i += 1) > n && return i
    @inbounds b = codeunit(s, i)
    b & 0xc0 ≠ 0x80 && return i
    ((i += 1) > n) | (l < 0xe0) && return i
    # second continuation byte
    @inbounds b = codeunit(s, i)
    b & 0xc0 ≠ 0x80 && return i
    ((i += 1) > n) | (l < 0xf0) && return i
    # third continuation byte
    @inbounds b = codeunit(s, i)
    ifelse(b & 0xc0 ≠ 0x80, i, i+1)
end

Base.@propagate_inbounds function Base.iterate(s::InlineString, i::Int=firstindex(s))
    (i % UInt) - 1 < ncodeunits(s) || return nothing
    b = @inbounds codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
    return iterate_continued(s, i, u)
end

function iterate_continued(s::InlineString, i::Int, u::UInt32)
    u < 0xc0000000 && (i += 1; @goto ret)
    n = ncodeunits(s)
    # first continuation byte
    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16
    # second continuation byte
    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8
    # third continuation byte
    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i)
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b); i += 1
@label ret
    return reinterpret(Char, u), i
end

Base.@propagate_inbounds function Base.getindex(s::InlineString, i::Int)
    b = codeunit(s, i)
    u = UInt32(b) << 24
    Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u)
    return getindex_continued(s, i, u)
end

function getindex_continued(s::InlineString, i::Int, u::UInt32)
    if u < 0xc0000000
        # called from `getindex` which checks bounds
        @inbounds isvalid(s, i) && @goto ret
        Base.string_index_err(s, i)
    end
    n = ncodeunits(s)

    (i += 1) > n && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 1
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 16

    ((i += 1) > n) | (u < 0xe0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 2
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b) << 8

    ((i += 1) > n) | (u < 0xf0000000) && @goto ret
    @inbounds b = codeunit(s, i) # cont byte 3
    b & 0xc0 == 0x80 || @goto ret
    u |= UInt32(b)
@label ret
    return reinterpret(Char, u)
end

Base.length(s::InlineString) = length_continued(s, 1, ncodeunits(s), ncodeunits(s))

Base.@propagate_inbounds function Base.length(s::InlineString, i::Int, j::Int)
    @boundscheck begin
        0 < i ≤ ncodeunits(s)+1 || throw(BoundsError(s, i))
        0 ≤ j < ncodeunits(s)+1 || throw(BoundsError(s, j))
    end
    j < i && return 0
    @inbounds i, k = thisind(s, i), i
    c = j - i + (i == k)
    length_continued(s, i, j, c)
end

@inline function length_continued(s::InlineString, i::Int, n::Int, c::Int)
    i < n || return c
    @inbounds b = codeunit(s, i)
    @inbounds while true
        while true
            (i += 1) ≤ n || return c
            0xc0 ≤ b ≤ 0xf7 && break
            b = codeunit(s, i)
        end
        l = b
        b = codeunit(s, i) # cont byte 1
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xe0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 2
        c -= (x = b & 0xc0 == 0x80)
        x & (l ≥ 0xf0) || continue

        (i += 1) ≤ n || return c
        b = codeunit(s, i) # cont byte 3
        c -= (b & 0xc0 == 0x80)
    end
end

# Parsers.xparse
function Parsers.xparse(::Type{T}, source::Union{AbstractVector{UInt8}, IO}, pos, len, options::Parsers.Options, ::Type{S}=T)::Parsers.Result{S} where {T <: InlineString, S}
    res = Parsers.xparse(String, source, pos, len, options)
    code = res.code
    overflowed = false
    poslen = res.val
    if !Parsers.valueok(code) || Parsers.sentinel(code)
        x = T()
    else
        poslen = res.val
        if T === InlineString1
            if poslen.len != 1
                overflowed = true
                x = T()
            else
                Parsers.fastseek!(source, poslen.pos)
                x = InlineString1(Parsers.peekbyte(source, poslen.pos))
                Parsers.fastseek!(source, pos + res.tlen - 1)
            end
        elseif Parsers.escapedstring(code) || !(source isa AbstractVector{UInt8})
            if poslen.len > (sizeof(T) - 1)
                overflowed = true
                x = T()
            else
                # manually build up InlineString
                i = poslen.pos
                maxi = i + poslen.len
                x = T()
                Parsers.fastseek!(source, i - 1)
                while i < maxi
                    b = Parsers.peekbyte(source, i)
                    if b == options.e
                        i += 1
                        Parsers.incr!(source)
                        b = Parsers.peekbyte(source, i)
                    end
                    x, overflowed = addcodeunit(x, b)
                    i += 1
                    Parsers.incr!(source)
                end
                Parsers.fastseek!(source, maxi)
            end
        else
            vlen = poslen.len
            if vlen > (sizeof(T) - 1)
                # @show T, vlen, sizeof(T)
                overflowed = true
                x = T()
            else
                # @show poslen.pos, vlen
                x = T(source, poslen.pos, vlen)
            end
        end
    end
    if overflowed
        code |= Parsers.OVERFLOW
    end
    return Parsers.Result{S}(code, res.tlen, x)
end

## InlineString sorting
using Base.Sort, Base.Order

# Only small-ish InlineStrings benefit from RadixSort algorithm
const SmallInlineStrings = Union{String1, String3, String7, String15}

# And under certain thresholds, MergeSort is faster than RadixSort, even for small InlineStrings
const MergeSortThresholds = Dict(
    1 => 2^5,
    4 => 2^7,
    8 => 2^9,
    16 => 2^23
)

struct InlineStringSortAlg <: Algorithm end
const InlineStringSort = InlineStringSortAlg()

Base.Sort.defalg(::AbstractArray{<:Union{SmallInlineStrings, Missing}}) = InlineStringSort

struct Radix
    size::Int
    pow::Int
    mask::UInt16
end

Radix(size) = Radix(size, 2^size, typemax(UInt16) >> (16 - size))

sortvalue(o::By,   x     ) = sortvalue(Forward, o.by(x))
sortvalue(o::Perm, i::Int) = sortvalue(o.order, o.data[i])
sortvalue(o::Lt,   x     ) = error("sortvalue does not work with general Lt Orderings")
sortvalue(rev::ReverseOrdering, x) = Base.not_int(sortvalue(rev.fwd, x))
sortvalue(::Base.ForwardOrdering, x) = x

_oftype(::Type{T}, x::S) where {T, S} = sizeof(T) == sizeof(S) ? Base.bitcast(T, x) : sizeof(T) > sizeof(S) ? Base.zext_int(T, x) : Base.trunc_int(T, x)

radix(v::T, j, radix_size, radix_mask) where {T} = _oftype(Int64, Base.and_int(Base.lshr_int(v, (j - 1) * radix_size), _oftype(T, radix_mask))) + 1

@noinline requireprimitivetype(T) = throw(ArgumentError("InlineStringSort requires isprimitivetype input: `$T` invalid"))

function Base.sort!(vs::AbstractVector, lo::Int, hi::Int, ::InlineStringSortAlg, o::Ordering)
    # Input checking
    lo >= hi && return vs

    # Make sure we're sorting a primitive type
    T = Base.Order.ordtype(o, vs)
    isprimitivetype(Base.nonmissingtype(T)) || requireprimitivetype(T)

    if hi - lo < MergeSortThresholds[sizeof(T)]
        return sort!(vs, lo, hi, MergeSort, o)
    end

    # setup
    ts = similar(vs)
    rdx = Radix(sizeof(T) == 1 ? 8 : 11)
    radix_size = rdx.size
    radix_mask = rdx.mask
    radix_size_pow = rdx.pow
    # iters is the # of 11-bit chunks we split each element up into
    # they each represent a "significant digit" we'll be sorting on
    iters = cld(sizeof(T) * 8, radix_size)
    # bin has a row for each unique 11-bit pattern
    # and a column for each 11-bit chunk we'll split each element up into
    bin = zeros(UInt32, radix_size_pow, iters)
    # if for some reason our lo isn't 1, we want to start our
    # 1st row bin values as the 1st index we'll start at in the output
    # i.e. we're assuming firstindex(vs):(lo - 1) is already sorted
    if lo > 1;  bin[1, :] .= lo-1; end

    # for each element, split into 11-bit chunks (radix)
    # and accumulate counts per unique pattern in bin
    for i = lo:hi
        v = sortvalue(o, vs[i])
        for j = 1:iters
            idx = radix(v, j, radix_size, radix_mask)
            @inbounds bin[idx, j] += 1
        end
    end

    # now we sort elements by sorting each radix using counting sort
    swaps = 0
    len = hi - lo + 1
    @inbounds for j = 1:iters
        # we first check if the radix for each element happened to be
        # the exact same bit pattern; if so, they're "already sorted"
        # for this radix and we can skip to the next. This would be common
        # if we, for example, had many small integer values stored in Int64
        # which would result in many "wasted" zero bits in most elements
        v = sortvalue(o, vs[hi])
        idx = radix(v, j, radix_size, radix_mask)

        # if every element was counted at this bit pattern
        # we can skip to the next radix chunk
        bin[idx, j] == len && continue

        # otherwise, we perform the counting sort for this radix
        # by doing a cumulative sum for this radix column in bin
        x = bin[1, j]
        for i = 2:radix_size_pow
            x += bin[i, j]
            bin[i, j] = x
        end
        # now we extract the output index for our 1st element (vs[hi])
        ci = bin[idx, j]
        # and decrement the count for that bit pattern which
        # will result in a subsequent identical bit pattern being
        # placed one index ahead of the current one
        bin[idx, j] -= 1
        ts[ci] = vs[hi]

        # now we sort the rest of the elements' radix similarly
        for i in (hi - 1):-1:lo
            v = sortvalue(o, vs[i])
            idx = radix(v, j, radix_size, radix_mask)
            ci = bin[idx, j]
            bin[idx, j] -= 1
            ts[ci] = vs[i]
        end
        # we keep 2 arrays, vs and ts
        # because we can't overwrite where the current
        # element will go in the output before we've sorted
        # the element already there
        vs, ts = ts, vs
        swaps += 1
    end

    if isodd(swaps)
        vs, ts = ts, vs
        @inbounds for i = lo:hi
            vs[i] = ts[i]
        end
    end
    return vs
end

# collections of InlineStrings
"""
    inlinestrings(itr) => Vector

    Utility function that takes any iterator of `AbstractString` values
and attempts to produce a `Vector` with a single promoted `InlineString` type. That is,
all iterated elements will be promoted to the smallest `InlineString` subtype
that can fit all elements. If any value is larger than the current largest InlineString
type (256 bytes), the entire collection will be promoted to `String` instead.
`missing` values are also allowed and will result in a result eltype of `Union{Missing, X}`
where `X` is an `InlineString` subtype or `String`.
"""
function inlinestrings(itr::T) where {T}
    # x must be iterable
    IS = Base.IteratorSize(T)
    state = iterate(itr)
    state === nothing && return []
    y, st = state
    x = y === missing ? missing : sizeof(y) < 256 ? InlineString(y) : String(y)
    eT = typeof(x)
    # allocate res, which will either be same length as `itr` if
    # IS <: HasLength, or length of 0 if Base.SizeUnknown
    res = allocate(eT, IS, itr)
    i = 1
    # set! push!-es for Base.SizeUnknown, or setindex! for HasLength
    set!(IS, res, x, i)
    i += 1
    # dispatch to separate function for type stability
    return _inlinestrings(itr, st, eT, IS, res, i)
end

const HasLength = Union{Base.HasShape, Base.HasLength}
allocate(::Type{T}, ::HasLength, itr) where {T} = Vector{T}(undef, length(itr))
allocate(::Type{T}, IS, itr) where {T} = Vector{T}(undef, 0)
set!(::HasLength, res, x, i) = setindex!(res, x, i)
set!(IS, res, x, i) = push!(res, x)

function _inlinestrings(itr, st, ::Type{eT}, IS, res, i) where {eT}
    while true
        state = iterate(itr, st)
        state === nothing && break
        y, st = state
        if y === missing && eT >: Missing
            set!(IS, res, missing, i)
        elseif y !== missing && sizeof(y) < sizeof(eT) || sizeof(y) == 1
            set!(IS, res, Base.nonmissingtype(eT)(y), i)
        else
            # need to promote and widen res,
            # then re-dispatch on _inlinestrings for new eltype
            x = y === missing ? missing : sizeof(y) < 256 ? InlineString(y) : String(y)
            new_eT = promote_type(typeof(x), eT)
            newres = allocate(new_eT, Base.HasLength(), res)
            copyto!(newres, 1, res, 1, i - 1)
            set!(IS, newres, x, i)
            return _inlinestrings(itr, st, new_eT, IS, newres, i + 1)
        end
        i += 1
    end
    return res
end

Base.Broadcast.broadcasted(::Type{InlineString}, A::AbstractArray) = inlinestrings(A)
Base.map(::Type{InlineString}, A::AbstractArray) = inlinestrings(A)
Base.collect(::Type{InlineString}, A::AbstractArray) = inlinestrings(A)

end # module


text:
Base.:*(::C, x::Real) where {C<:Currency} = Cash{C}(x)
Base.:*(x::Real, ::C) where {C<:Currency} = Cash{C}(x)

Base.:-(x::Cash{C}) where {C<:Currency} = Cash{C}(-x.value)

Base.:+(x::Cash{C}, y::Cash{C}) where {C<:Currency} = Cash{C}(x.value + y.value)
Base.:-(x::Cash{C}, y::Cash{C}) where {C<:Currency} = Cash{C}(x.value - y.value)

Base.:*(c::Cash{C}, x::Real) where {C<:Currency} = Cash{C}(c.value * x)
Base.:*(x::Real, c::Cash{C}) where {C<:Currency} = Cash{C}(c.value * x)
Base.:*(x::Cash{C}, y::Cash{C}) where {C<:Currency} = Cash{C}(x.value * y.value)

Base.:/(c::Cash{C}, x::Real) where {C<:Currency} = Cash{C}(c.value / x)
Base.:/(c1::Cash{C}, c2::Cash{C}) where {C<:Currency} = c1.value / c2.value

Base.:(==)(x::Cash{C}, y::Cash{C}) where {C<:Currency} = x.value == y.value
Base.:(<)(x::Cash{C}, y::Cash{C}) where {C<:Currency} = x.value < y.value
Base.:(<=)(x::Cash{C}, y::Cash{C}) where {C<:Currency} = x.value <= y.value

Base.isless(x::Cash{C}, y::Cash{C}) where {C<:Currency} = isless(x.value, y.value)
Base.isequal(x::Cash{C}, y::Cash{C}) where {C<:Currency} = isequal(x.value, y.value)

Base.abs(x::Cash{C}) where {C<:Currency} = Cash{C}(abs(x.value))

Base.isapprox(x::Cash{C}, y::Cash{C}; kwargs...) where {C<:Currency} = isapprox(x.value, y.value; kwargs...)

Base.zero(::Cash{C,T}) where {C,T} = Cash{C,T}(zero(T))
Base.zero(::Type{Cash{C,T}}) where {C,T} = Cash{C,T}(zero(T))
Base.one(::Cash{C,T}) where {C,T} = Cash{C,T}(one(T))
Base.one(::Type{Cash{C,T}}) where {C,T} = Cash{U,T}(one(T))

Base.trunc(x::Cash{C}; digits=0) where {C<:Currency} = Cash{C}(trunc(x.value, digits=digits))
Base.round(x::Cash{C}; digits=0) where {C<:Currency} = Cash{C}(round(x.value, digits=digits))


text:
function \(cur_matrix::FourBlockMatrix, cur_vector::Vector)

  A, B, C, D = get_sub_matrices(cur_matrix)

  iszero(cur_matrix.n) && return D \ cur_vector
  iszero(cur_matrix.m) && return A \ cur_vector

  top_b = view(cur_vector, 1:cur_matrix.n)
  bot_b = view(cur_vector, cur_matrix.n .+ (1:cur_matrix.m) )

  cur_lu = lu(A)
  F = SharedArray{Float64}(cur_matrix.m,cur_matrix.m)

  @sync @distributed for j = 1:size(B,2)
    tmp_B = B[:,j]
    isempty(tmp_B) && continue
    F[:,j] = C * ( cur_lu \ Vector(tmp_B) )
  end

  G = D - F

  L_1 = A \ top_b

  L_2 = G \ ( C * L_1 )

  top_left = L_1 + A \ ( B * L_2 )
  bot_left = - L_2

  R_1 = G \ bot_b

  top_right = - A \ ( B * R_1 )
  bot_right = R_1

  vcat(
    ( top_left .+ top_right ),
    ( bot_left .+ bot_right )
  )

end


text:
"""
    ObjectId

Complete, unique ID of a `ModelObject`

Contains own id and a vector of parent ids, so one knows exactly where the object
is placed in the model tree.
"""
struct ObjectId
    # Store IDs as vector, not tuple (b/c empty tuples are tricky)
    # "Youngest" member is positioned last in vector
    ids :: Vector{SingleId}
end

## ---------  Constructors

"""
	$(SIGNATURES)

Construct an ObjectId from a vector of `SingleId`.
"""
ObjectId() = ObjectId(Vector{SingleId}());

# Without a parent or index
function ObjectId(ownId :: SingleId)
    return ObjectId([ownId]);
end

# With parent; no index
function ObjectId(name :: Symbol, parentIds :: ObjectId = ObjectId())
    return ObjectId(vcat(parentIds.ids, SingleId(name)))
end

# With everything
function ObjectId(name :: Symbol, index :: Vector{T1},
    parentIds :: ObjectId = ObjectId()) where T1 <: Integer

    return ObjectId(vcat(parentIds.ids, SingleId(name, index)))
end

function ObjectId(name :: Symbol, idx :: T1,
    parentIds :: ObjectId = ObjectId()) where T1 <: Integer

    return ObjectId(vcat(parentIds.ids,  SingleId(name, [idx])))
end

ObjectId(name :: Symbol, descr :: String,
    parentIds :: ObjectId = ObjectId()) =
    ObjectId(vcat(parentIds.ids,  SingleId(name, descr)));

## ------  Parent info

function has_parent(oId :: ObjectId)
    return length(oId.ids) > 1
end

function get_parent_id(oId :: ObjectId)
    if has_parent(oId)
        return ObjectId(oId.ids[1 : (end-1)])
    else
        return ObjectId()
    end
end

is_parent_of(pId :: ObjectId,  oId :: ObjectId) = isequal(pId, get_parent_id(oId))

"""
	$(SIGNATURES)

Number of parents of an `ObjectId`.
"""
n_parents(pId :: ObjectId) = length(pId.ids) - 1;

"""
	$(SIGNATURES)

Make child ObjectId from parent ObjectId.
"""
function make_child_id(obj :: T1, name :: Symbol,
    index :: Vector{T2} = Vector{Int}()) where {T1, T2 <: Integer}

    return ObjectId(name, index, obj.objId)
end

# Make child ID from parent's ID
function make_child_id(parentId :: ObjectId, name :: Symbol,
    index :: Vector{T2} = Vector{Int}()) where {T2 <: Integer}

    return ObjectId(name, index, parentId)
end

make_child_id(parentId :: ObjectId, name :: Symbol, descr :: String) =
    ObjectId(name, descr, parentId);

## -----------  Basic properties

Base.broadcastable(oId :: ObjectId) = Ref(oId);

"""
	$(SIGNATURES)

Checks whether two `ObjectId`s are the same. Does not consider descriptions.
"""
function Base.:(==)(id1 :: ObjectId,  id2 :: ObjectId)
    outVal = (length(id1.ids) == length(id2.ids))  &&  all(isequal.(id1.ids, id2.ids))
    return outVal
end

Base.isequal(id1 :: ObjectId,  id2 :: ObjectId) = id1 == id2;

# Must be implmented for storage in `Dict`.
Base.hash(id1 :: ObjectId, h :: UInt) =
    hash(id1.ids, hash(:ObjectId, h));

"""
	$(SIGNATURES)

Description of an `ObjectId`.
"""
description(oId :: ObjectId) = description(own_id(oId));

own_index(oId :: ObjectId) = oId.ids[end].index
# Return own SingleId
own_id(oId :: ObjectId) = oId.ids[end];

"""
	$(SIGNATURES)

Return object's own name as `Symbol`.
"""
function own_name(oId :: ObjectId)
    return name(own_id(oId))
end

"""
	$(SIGNATURES)

Make string from ObjectId. Such as "p > q > r[4, 2]".
"""
function make_string(id :: ObjectId)
    outStr = "";
    for i1 = 1 : length(id.ids)
        if i1 > 1
            outStr = outStr  * ObjIdSeparator;
        end
        outStr = outStr * ModelObjectsLH.make_string(id.ids[i1]);
    end
    return outStr
end

"""
	$(SIGNATURES)

The inverse of `make_string`.
"""
function make_object_id(s :: T1) where T1 <: AbstractString
    if occursin(ObjIdSeparator, s)
        strV = split(s, ObjIdSeparator);
        singleIdV = similar(strV, SingleId);
        for (j, str) in enumerate(strV)
            singleIdV[j] = make_single_id(str);
        end
        return ObjectId(singleIdV)
    else
        return ObjectId(Symbol(s));
    end
end

"""
	$(SIGNATURES)

Show an object id.
"""
function Base.show(io :: IO,  id :: ObjectId)
    print(io,  "ObjectId: " * make_string(id));
end

# ----------------


text:
# --- PRIVATE METHODS BELOW HERE ----------------------------------------------------------------- #
function _compute_crr_index_array(numberOfLevels::Int64)

    # ok, so lets build an index array -
    number_items_per_level = [i for i = 1:numberOfLevels]
    tmp_array = Array{Int64,1}()
    theta = 0
    for value in number_items_per_level
        for index = 1:value
            push!(tmp_array, theta)
        end
        theta = theta + 1
    end

    N = sum(number_items_per_level[1:(numberOfLevels - 1)])
    index_array = Array{Int64,2}(undef, N, 3)
    for row_index = 1:N

        # index_array[row_index,1] = tmp_array[row_index]
        index_array[row_index,1] = row_index
        index_array[row_index,2] = row_index + 1 + tmp_array[row_index]
        index_array[row_index,3] = row_index + 2 + tmp_array[row_index]
    end

    return index_array
end

function _american_decision_logic(currentValue::Float64, futureValue::Float64)::Float64
    return max(currentValue, futureValue)
end

function _crr_movement_logic(latticeModel::VLBinomialLattice)::NamedTuple

    # initialize - get stuff from the lattice model -
    µ = latticeModel.µ
    σ = latticeModel.σ
    𝝙t = latticeModel.𝝙t

    # compute the movement up, down and probability -
    u = exp(σ * √𝝙t)
    d = exp(-σ * √𝝙t)
    p = (exp(µ * 𝝙t) - d) / (u - d)
    dfactor = exp(-µ * 𝝙t)

    # setup the tuple -
    movement_tuple = (up=u, down=d, probability=p, discount=dfactor)

    # return -
    return movement_tuple
end

# --- PRIVATE METHODS ABOVE HERE ----------------------------------------------------------------- #

# --- PUBLIC METHODS BELOW HERE ------------------------------------------------------------------ #
"""
    binomial_price(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPrice::Float64;
        decisionLogic::Function = _american_decision_logic, movementLogic::Function = _crr_movement_logic) -> VLResult

Magical description goes here!
"""
function binomial_price(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPrice::Float64;
    decisionLogic::Function=_american_decision_logic, movementLogic::Function = _crr_movement_logic)::VLResult

    try

        # initialize -
        numberOfLevels = latticeModel.numberOfLevels
        crr_index_array = _compute_crr_index_array(numberOfLevels)
        total_number_of_lattice_nodes = crr_index_array[end,end]    # the number of nodes will be the end, end element of the index array -
        number_of_nodes_to_evaluate = crr_index_array[end,1]        # the number of nodes that we need to evaluate is goven by the end x 1 element

        # compute the movement using the movement function -
        movement_tuple = movementLogic(latticeModel)
        u = movement_tuple.up
        d = movement_tuple.down
        p = movement_tuple.probability
        dfactor = movement_tuple.discount

        # init the tree -
        tree_value_array = Array{Float64,2}(undef, total_number_of_lattice_nodes, 3) # nodes x 3 = col1: underlying price, col2: intrinsic value, col3: option price

        # First: let's compute the underlying price on the lattice -
        tree_value_array[1,1] = underlyingPrice
        for node_index = 1:number_of_nodes_to_evaluate

            # get index -
            parent_node_index = crr_index_array[node_index,1]
            up_node_index = crr_index_array[node_index,2]
            down_node_index = crr_index_array[node_index,3]

            # compute prices -
            parent_price = tree_value_array[parent_node_index,1]
            up_price = parent_price * u
            down_price = parent_price * d

            # store prices -
            tree_value_array[up_node_index,1] = up_price
            tree_value_array[down_node_index,1] = down_price
        end

        # Next: compute the intrinsice value for each node -
        for node_index = 1:total_number_of_lattice_nodes

             # ok, get the underlying price -
            underlying_price_value = tree_value_array[node_index,1]

            # compute the intrinsic value -
            iv_value = intrinsic_value(contractSet, underlying_price_value) |> check

            # capture -
            tree_value_array[node_index,2] = iv_value
            tree_value_array[node_index,3] = iv_value
        end

        # Last: compute the option price -
        reverse_node_index_array = range(number_of_nodes_to_evaluate, stop=1, step=-1) |> collect
        for (_, parent_node_index) in enumerate(reverse_node_index_array)

            # ok, get the connected node indexes -
            up_node_index = crr_index_array[parent_node_index,2]
            down_node_index = crr_index_array[parent_node_index,3]

            # ok, let's compute the payback *if* we continue -
            future_payback = dfactor * (p * tree_value_array[up_node_index,3] + (1 - p) * tree_value_array[down_node_index,3])
            current_payback = tree_value_array[parent_node_index,2]

            # use the decision logic to compute price -
            node_price = decisionLogic(current_payback, future_payback)

            # capture -
            tree_value_array[parent_node_index,3] = node_price
        end

        # return -
        return VLResult(tree_value_array)
    catch error
        return VLResult(error)
    end
end

"""
    binomial_price(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPriceArray::Array{Float64,1};
        decisionLogic::Function = _american_decision_logic, movementLogic::Function = _crr_movement_logic) -> VLResult

Magical description goes here!
"""
function binomial_price(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPriceArray::Array{Float64,1};
    decisionLogic::Function=_american_decision_logic, movementLogic::Function = _crr_movement_logic)::VLResult

    try

        # initialize -
        number_of_underlying_prices = length(underlyingPriceArray)
        contract_price_array = Array{Float64,2}(undef, number_of_underlying_prices, 2)

        # ok, so let's simulate updating the strike price -
        for (index, underlying_price) in enumerate(underlyingPriceArray)

            # compute -
            price_tree = binomial_price(contractSet, latticeModel, underlying_price; decisionLogic=decisionLogic, movementLogic=movementLogic) |> check
            price_value = price_tree[1,3]   # the value of the contract is *always* the 1,3 element -

            # grab -
            contract_price_array[index,1] = underlying_price
            contract_price_array[index,2] = price_value
        end

        # return -
        return VLResult(contract_price_array)
    catch error
        return VLResult(error)
    end
end

"""
    binomial_price(contract::Union{VLCallOptionContract,VLPutOptionContract}, latticeModel::VLBinomialLattice,
        underlyingPrice::Float64, strikPriceArray::Array{Float64,1};
            decisionLogic::Function = _american_decision_logic, movementLogic::Function = _crr_movement_logic) -> VLResult

Magical description goes here!
"""
function binomial_price(contract::Union{VLCallOptionContract,VLPutOptionContract}, latticeModel::VLBinomialLattice,
    underlyingPrice::Float64, strikPriceArray::Array{Float64,1};
        decisionLogic::Function=_american_decision_logic, movementLogic::Function = _crr_movement_logic)::VLResult

    try

        # initialize -
        number_of_strike_prices = length(strikPriceArray)
        contract_price_array = Array{Float64,2}(undef, number_of_strike_prices, 2)

        # ok, so lets set the strike price and process the contract =
        for (index, strike_price) in enumerate(strikPriceArray)

            # set the price -
            contract.strikePrice = strike_price

            # package -
            contract_set = Set{VLAbstractAsset}()
            push!(contract_set, contract)

            # calculate the contract price -
            price_tree = binomial_price(contract_set, latticeModel, underlyingPrice; decisionLogic=decisionLogic, movementLogic=movementLogic) |> check
            price_value = price_tree[1,3]   # the value of the contract is *always* the 1,3 element -

            # grab -
            contract_price_array[index,1] = strike_price
            contract_price_array[index,2] = price_value
        end

        # return -
        return VLResult(contract_price_array)
    catch error
        return VLResult(error)
    end
end

# --- PUBLIC METHODS ABOVE HERE ------------------------------------------------------------------ #


text:
using MicroLogging

function checkerror(err)
    if err[] != C_NULL
        funcname = unsafe_string(rs2_get_failed_function(err[]))
        funcargs = unsafe_string(rs2_get_failed_args(err[]))
        message = unsafe_string(rs2_get_error_message(err[]))
        @error "rs_error was raised when calling $funcname($funcargs):\n    $message"
    end
end

function deviceinfo(dev)
    err = Ref{Ptr{rs2_error}}(0)
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_NAME, err)
    checkerror(err)
    @info "This device is an $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_SERIAL_NUMBER, err)
    checkerror(err)
    @info "  Serial number: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_FIRMWARE_VERSION, err)
    checkerror(err)
    @info "  Firmware version: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_PHYSICAL_PORT, err)
    checkerror(err)
    @info "  Physical port: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_ADVANCED_MODE, err)
    checkerror(err)
    @info "  Advanced mode: $(unsafe_string(device_info))"
    device_info = rs2_get_device_info(dev, RS2_CAMERA_INFO_PRODUCT_ID, err)
    checkerror(err)
    @info "  Product ID: $(unsafe_string(device_info))"
end


text:
using Revise
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters
using Test
using SparseArrays
using BenchmarkTools
using Profile
# using ProfileView
using InteractiveUtils

using Base.Iterators: product, OneTo

using AxisAlgorithms: A_mul_B_md!
using ShaleDrillingLikelihood: beta_1minusbeta, update_IminusTVp!, IminusTEVp
using LinearAlgebra
using IterativeSolvers
using IncompleteLU

using ShaleDrillingLikelihood: DCDPTmpVars,
    q,
    lse,
    tmp,
    ubV,
    tmp_cart,
    Πψtmp,
    IminusTEVp,
    ubVfull,
    reward,
    dp1space,
    sprimes,
    vfit!,
    dEV,
    EV,
    anticipate_t1ev,
    _dmax,
    solve_inf_vfit!,
    solve_inf_pfit!,
    solve_inf_vfit_pfit!,
    gradinf!,
    update_static_payoffs!

nz = 51
nψ = 51
num_zt = 15

royalty_rates = RoyaltyRatesSmall()
rwrd    = DefaultDrillReward()
θ_drill   = Theta(DefaultDrillReward())
wp = FullLeasedProblem()
anticipate = true
seed=1234

z = TestPriceProcess(;nsim=num_zt, ngrid=nz)
_zchars = ShaleDrillingLikelihood.SDLParameters.series(z)

ddm  = TestDynamicDrillModel(z, anticipate; reward=rwrd)
tmpv = DCDPTmpVars(ddm)
evs = ValueFunction(ddm)

itype = (4.67, 0.25,)
dograd = false

i = length(wp) - 1 # state_index

idxd  = dp1space(wp,i)
idxs  = collect(sprimes(statespace(ddm),i))

EV0   = view(EV(evs) ,    :, :, i)
EV1   = view(EV(evs) ,    :, :, idxs)
dEV0  = view(dEV(evs), :, :, :, i)
dEV1  = view(dEV(evs), :, :, :, idxs)

fill!(evs, 0)

tmp_vw = view(tmpv, idxd)
update_static_payoffs!(tmp_vw, ddm, θ_drill, i, itype, dograd)

ubV(tmp_vw)  .+= discount(ddm) .* EV1

converged, iter, bnds = solve_inf_vfit!(EV0, tmp_vw, ddm; maxit=10, vftol=1e-5)

let t = tmp_vw
    bnds = ShaleDrillingLikelihood.pfit!(EV0, t, ddm; vftol=1e-10)
    ubVfull(t)[:,:,1] .= discount(ddm) .* EV0
end

t = tmp_vw
vftol = 1e-16

ΔEV = lse(t)
q0 = ubV(t)

if anticipate_t1ev(ddm)
    logsumexp_and_softmax!(ΔEV, q0, tmp(t), ubV(t), 1)
else
    findmax!(add_1_dim(ΔEV), add_1_dim(tmp_cart(t)), ubV(t))
    q0[:,:,1] .= last.(getfield.(tmp_cart(t), :I)) .== 1         # update q0 as Pr(d=0|x)
end
A_mul_B_md!(tmp(t), ztransition(ddm), ΔEV, 1)

# compute difference & check bnds
bnds = extrema(ΔEV .= EV0 .- tmp(t)) .* -beta_1minusbeta(ddm)
if all(abs.(bnds) .< vftol)
    EV0 .= tmp(t)
    return bnds
end

j = 1
q0j  = view(q0, :, j, 1)
ΔEVj_orig = view(ΔEV, :, j)
ΔEVj = similar(ΔEVj_orig)

# Consider https://juliamath.github.io/IterativeSolvers.jl/dev/preconditioning/#Preconditioning-1
# with https://github.com/haampie/IncompleteLU.jl as preconditioner
# https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method
# https://en.wikipedia.org/wiki/Generalized_minimal_residual_method
# Econ paper - https://doi.org/10.1016/S1474-6670(17)33089-6
# Mrkaic and Pauletto (2001) Preconditioning in Economic Stochastic Growth Models
update_IminusTVp!(t, ddm, q0j)

function solvedirect(t, ΔEVj, ΔEVj_orig)
    ΔEVj .= IminusTEVp(t) \ ΔEVj_orig
end

using Preconditioners

function solveindirect(t, ΔEVj, ΔEVj_orig; kwargs...)
    # bicgstabl!(ΔEVj, IminusTEVp(t), ΔEVj_orig; kwargs...)
    gmres!(ΔEVj, IminusTEVp(t), ΔEVj_orig; kwargs...)
end

function solveindirect_plu(t, ΔEVj, ΔEVj_orig, tau = 0.02; kwargs...)
    LU = ilu(IminusTEVp(t); τ = tau)
    # bicgstabl!(ΔEVj, IminusTEVp(t), ΔEVj_orig; Pl=LU, kwargs...)
    gmres!(ΔEVj, IminusTEVp(t), ΔEVj_orig; Pl=LU, kwargs...)
end

function solveindirect_diag(t, ΔEVj, ΔEVj_orig, dp::DiagonalPreconditioner; kwargs...)
    A = IminusTEVp(t)
    # dp = DiagonalPreconditioner(A)
    UpdatePreconditioner!(dp, A)
    # bicgstabl!(ΔEVj, A, ΔEVj_orig; Pl=LU, kwargs...)
    gmres!(ΔEVj, A, ΔEVj_orig; Pl=dp, kwargs...)
end

dp = DiagonalPreconditioner(IminusTEVp(t))

tol = 1e-13
@btime solvedirect(  $t, $ΔEVj, $ΔEVj_orig)
@btime solveindirect($t, $ΔEVj, $ΔEVj_orig; tol=tol)
@btime solveindirect_plu($t, $ΔEVj, $ΔEVj_orig; tol=tol)
@btime solveindirect_diag($t, $ΔEVj, $ΔEVj_orig, $dp; tol=tol)

A = IminusTEVp(t)
ILU = ilu(A; τ = 1e-4)
LU = lu(A)
@btime ilu($A; τ = 1e-4)
@btime lu($A)
@btime gmres!($ΔEVj, $A, $ΔEVj_orig; Pl=$ILU)
@btime solvedirect(  $t, $ΔEVj, $ΔEVj_orig)

# for tau in (2.0, 1.0, 0.1, 0.01, 0.001)
#     @btime solveindirect_plu($t, $ΔEVj, $ΔEVj_orig, $tau)
# end

bicgstabl!(ΔEVj, IminusTEVp(t), ΔEVj_orig; log=true)

DiagonalPreconditioner(IminusTEVp(t))

for tau in (1.0, 0.1, 0.02, 0.01, 0.001)
    LU = ilu(IminusTEVp(t); τ = tau)
    xplu, hist = gmres!(ΔEVj, IminusTEVp(t), ΔEVj_orig; Pl=LU, log=true)
    @show hist
end

x, hist = gmres!(ΔEVj, IminusTEVp(t), ΔEVj_orig; log=true)
@show hist

UpdatePreconditioner!(dp, IminusTEVp(t))
x, hist = gmres!(ΔEVj, IminusTEVp(t), ΔEVj_orig; Pl=dp, log=true)
@show hist

# converged, iterp, bnds = solve_inf_vfit_pfit!(EV0, tmp_vw, ddm; vftol=1e-10, maxit0=20, maxit1=40)


text:
export @intrinsic

# NOTE:
# we only define basic staff for now
# but in principal we should generate simulation
# instruction and QASM directly from here
# since these are not real "intrinsic" APIs

macro intrinsic(ex)
    return esc(intrinsic_m(ex))
end

function intrinsic_m(name::Symbol)
    tname = Symbol(name, :Gate)
    return quote
        Core.@__doc__ struct $tname <: $IntrinsicRoutine end
        Core.@__doc__ const $name = $tname()
    end
end

is_one_qubit_gate(x) = false

# function intrinsic_m(ex::Expr)
#     ex.head === :call || error("expect a function call or a symbol")
#     name = ex.args[1]::Symbol
#     args = rm_annotations.(ex.args[2:end])

#     body = Expr(:block)
#     def = Expr(:struct, false, name, body)

#     return quote
#         Core.@__doc__ const $name = $IntrinsicRoutine{$(QuoteNode(name))}()

#         function (self::$IntrinsicRoutine{$(QuoteNode(name))})($(ex.args[2:end]...))
#             return $IntrinsicSpec(self, $(rm_annotations.(ex.args[2:end])...))
#         end
#     end
# end

module Intrinsics

using MLStyle
import YaoCompiler
using ..YaoCompiler: @intrinsic, IntrinsicRoutine
export X, Y, Z, H, S, T, shift, Rx, Ry, Rz

@intrinsic X
@intrinsic Y
@intrinsic Z
@intrinsic H
@intrinsic S
@intrinsic T

YaoCompiler.is_one_qubit_gate(::Type{<:XGate}) = true
YaoCompiler.is_one_qubit_gate(::Type{<:YGate}) = true
YaoCompiler.is_one_qubit_gate(::Type{<:ZGate}) = true
YaoCompiler.is_one_qubit_gate(::Type{<:HGate}) = true
YaoCompiler.is_one_qubit_gate(::Type{<:SGate}) = true
YaoCompiler.is_one_qubit_gate(::Type{<:TGate}) = true

# TODO: implement intrinsic macro for function calls
# @intrinsic shift(θ::Real)
# @intrinsic Rx(θ::Real)
# @intrinsic Ry(θ::Real)
# @intrinsic Rz(θ::Real)

for name in [:shift, :Rx, :Ry, :Rz]
    @eval struct $name{T<:Real} <: IntrinsicRoutine
        θ::T

        $name(θ) = new{typeof(θ)}(θ)
    end

    @eval YaoCompiler.is_one_qubit_gate(::Type{<:$name}) = true
end

struct UGate{T} <: IntrinsicRoutine
    α::T
    β::T
    γ::T
end

YaoCompiler.is_one_qubit_gate(::Type{<:UGate}) = true
# this is a workaround since MLStyle doesn't support a.b pattern

@as_record XGate
@as_record YGate
@as_record ZGate
@as_record HGate
@as_record SGate
@as_record TGate

@as_record shift
@as_record Rx
@as_record Ry
@as_record Rz
@as_record UGate

end


text:
__precompile__()

module AsymptoticPosteriors

using   LinearAlgebra, # special functions and LinearAlgebra ought to be clear enough.
        SpecialFunctions,
        PaddedMatrices, VectorizationBase, SIMDPirates,
        DifferentiableObjects,
        Statistics,
        ForwardDiff, DiffResults

# import

using PaddedMatrices: AbstractFixedSizeVector,
                    AbstractMutableFixedSizeVector,
                    AbstractMutableFixedSizeMatrix,
                    MutableFixedSizeVector
using DifferentiableObjects: AbstractDifferentiableObject

export AsymptoticPosteriorFD, mode

# const SizedArray{P,T} = Union{SizedSIMDVector{P,T}, }

# debug() = true
# debug_rootsearch() = true
debug() = false
debug_rootsearch() = false

include("function_wrappers.jl")
include("brent.jl")
include("MAP.jl")
include("AsymptoticPosteriorInterface.jl")
include("AsymptoticPosteriorConstructors.jl")
# include("nested.jl")
include("initial_root_search.jl")
# include("plot_recipes.jl")

end # module


text:
@testset "FullSearch CUDA" begin
    original = Float32.(imresize(testimage("mri-stack"), (64, 64)));
    ref = original[:, :, 1];
    img = original[:, :, 2];

    @testset "best_match" begin
        for T in (Float32, N0f8, Gray{Float32}, Gray{N0f8})
            cu_ref = CuArray(T.(ref))
            cu_img = CuArray(T.(img))
            # Although supported, not all metrics are useful for blockmatching purposes
            for d in [
                Euclidean(),
                SqEuclidean(),
                Cityblock(),
                TotalVariation(),
                Minkowski(1.5f0),
            ]
                S = FullSearch(d, 2, patch_radius=2, search_radius=5)
                cu_matches = OffsetArrays.parent_call(Array, best_match(S, cu_ref, cu_img))
                matches = best_match(S, ref, img)

                # Because the inner loop orders are different, the output results of
                # CPU and CUDA version are different in constant fields.
                @test_broken cu_matches == matches
                @test sum(cu_matches .!= matches)/length(matches) < 0.5

                cu_matches = OffsetArrays.parent_call(Array, best_match(S, cu_ref, cu_img; offset=true))
                matches = best_match(S, ref, img; offset=true)
                @test_broken cu_matches == matches
                @test sum(cu_matches .!= matches)/length(matches) < 0.5
            end
        end
    end

    @testset "multi_match" begin
        for T in (Float32, N0f8, Gray{Float32}, Gray{N0f8})
            cu_ref = CuArray(T.(ref))
            cu_img = CuArray(T.(img))

            # Although supported, not all metrics are useful for blockmatching purposes
            for d in [
                Euclidean(),
                SqEuclidean(),
                Cityblock(),
                TotalVariation(),
                Minkowski(1.5f0),
            ]
                S = FullSearch(d, 2, patch_radius=2, search_radius=5)
                num_patches = 20
                cu_matches = OffsetArrays.parent_call(Array, multi_match(S, cu_ref, cu_img; num_patches))
                matches = multi_match(S, ref, img; num_patches)

                # Because the inner loop orders are different, the output results of
                # CPU and CUDA version are different in constant fields.
                @test_broken cu_matches == matches
                diffcount = mapreduce(+, CartesianIndices(axes(matches)[1:2])) do p
                    cu_v = Set(view(cu_matches, p, :))
                    v = Set(view(matches, p, :))
                    length(cu_v) - length(intersect(cu_v, v))
                end
                @test diffcount/length(matches) < 0.02

                cu_matches = OffsetArrays.parent_call(Array, multi_match(S, cu_ref, cu_img; num_patches, offset=true))
                matches = multi_match(S, ref, img; num_patches, offset=true)
                @test_broken cu_matches == matches
                diffcount = mapreduce(+, CartesianIndices(axes(matches)[1:2])) do p
                    cu_v = Set(view(cu_matches, p, :))
                    v = Set(view(matches, p, :))
                    length(cu_v) - length(intersect(cu_v, v))
                end
                @test diffcount/length(matches) < 0.02
            end
        end
    end
end