text:
module KohnLeiboviciSzkup3 #introduces a new global scope (see Julia documentation for more info.)

## ##################################################
# This code solves the model in
# No Credit, No Gain: Trade Liberalization Dynamics, Production Inputs, and Financial Development
# Kohn, Leibovici, Szkup 2021
####################################################

# Dependencies (packages used)
using LinearAlgebra, Distributions
using Plots # used to plot graphs
using Parameters # used to create named tuples
using NLsolve # (multivariable) root-finder [main solver]
using LeastSquaresOptim # alternative solver [transition only]
using QuantEcon #used to generate a markov chain (tauchen discretization)
using Dates #used to save the date of different workspaces
using JLD2 #allows to save the workspace or individual variables in .jld2 format (akin to matlab's .mat)
using MAT #allows to save the workspace or individual variables in .mat format (useful for matlab comparisons and the welfare graphs)

### Benchmarking packages
#using Profile
#using TimerOutputs.jl # allows to print nicely formatted tables with the timing of different sections of the code

# Model files
include("parameters_settings.jl") #contains the default setup of settings and parameters (+tauchen function)
include("functions.jl") #contains all functions used in the computation of the SS
include("functions_trans.jl") #contains (additional) functions used in the transition

p=parameters_default() #instantiating the parameters with default values
s=settings_default() #instantiating the settings with default values

say(what) = run(`osascript -e "say \"$(what)\""`, wait=false)
#function that lets you know when code finished running

#Step 1: Solve initial steady state

    if s.tariffsincome == 1
        if s.load_prices ==0
            const g1=(w= 0.021645899110592, # In Julia, it's better perfomance-wise not to use global variables such as this. If they are used and they are constants (no not change through the code), labeling them as such using "const" helps mitigate the effect on performance.
            Yc = 0.075394893630123,
            ξ = 0.280906535178547,
            Pk = 0.969892818482190,
            Yk = 0.065988921968899)
        elseif s.load_prices ==1
            @load "KLS3_Prices.jld2" Prices
            const g1=(w = Prices[2,1],
            Yc = Prices[1,1],
            ξ = Prices[3,1],
            Pk = Prices[4,1],
            Yk = Prices[5,1])
        end
        x0 = [g1.w g1.Yc g1.ξ g1.Pk g1.Yk]

    else

        if s.load_prices ==0
            const g1=(w = 0.035690209553739,
            ϕ_h = 0.16403,
            ξ =0.369060692554216,
            Pk = 0.975043532756500)
        elseif s.load_prices ==1
            @load "KLS3_Prices.jld2" Prices
            const g1=(w = Prices[2,1],
            ϕ_h = Prices[1,1],
            ξ = Prices[3,1],
            Pk = Prices[4,1])
        end
        x0 = [g1.w g1.ϕ_h g1.ξ g1.Pk]

    end

@time begin
    if s.GE == 1
        results_GE =
         nlsolve((F,x) ->KLS3_GE_par!(F,x,p,s),log.(x0),autodiff = :forward, method=s.method_GE, xtol=s.xtol_GE, ftol=s.ftol_GE, show_trace=s.show_trace_GE)
         mcc_0, p_o_0, sim_0 = KLS3_GE_par(results_GE.zero,p,s)
    elseif s.GE ==0
        mcc_0, p_o_0, sim_0 =KLS3_GE_par(log.(x0),p,s)
    end
    say("Initial Steady State Finished - you'd better come and take a look....")
end

#clearing previously defined global vars (so they don't retain their values when computing the 2nd SS)
guessM=nothing
guessV=nothing

    if s.tariffsincome_PE_flag==1
        s=settings_default(tariffsincome_PE = 1) #note that here what you do is instantiate "setting_default" as before but changing the default value of tariffsincome_PE (instead of changing the value in the previous instantiation, which cannot be done)
    end

# Step 2: Solve final steady state
    if s.transition == 1

        #Shocks
        if s.tariffsincome_PE_flag==1
            p=parameters_default(tariffsincome=p_o_0.tariffsincome,
            #The change in p.tarrifsincome was originally in the previous block of code. However, the fact that what we do here is instantiate p again, it's convenient to do all the changes to p together (otherwise, we'd have to add the change in tariffsincome again)
            r = p.rv[end],
            Pf = p.Pfv[end],
            Yf = p.Yfv[end],
            β = p.β_v[end],
            θ = p.θ_v[end],
            τ = p.τ_v[end],
            δ = p.δ_v[end],
            τ_x = p.τ_x_v[end],
            τ_m_c = p.τ_m_c_v[end],
            τ_m_k = p.τ_m_k_v[end])
        else
            p=parameters_default(r = p.rv[end],
            Pf = p.Pfv[end],
            Yf = p.Yfv[end],
            β = p.β_v[end],
            θ = p.θ_v[end],
            τ = p.τ_v[end],
            δ = p.δ_v[end],
            τ_x = p.τ_x_v[end],
            τ_m_c = p.τ_m_c_v[end],
            τ_m_k = p.τ_m_k_v[end])
        end

         if s.tariffsincome == 1

            if s.load_prices ==0  &&  s.PE==0
                 const g2=(w = 0.021645899110592,
                Yc = 0.075394893630123,
                ξ = 0.280906535178547,
                Pk = 0.969892818482190,
                Yk = 0.065988921968899)
            elseif s.load_prices ==1
                @load "KLS3_Prices.jld2" Prices
                const g2=(w = Prices[2,end],
                Yc = Prices[1,end],
                ξ = Prices[3,end],
                Pk = Prices[4,end],
                Yk = Prices[5,end])
            elseif s.load_prices ==0  &&  s.PE==1
                const g2=(w = p_o_0.w,
                Yc = p_o_0.Yc,
                ξ = p_o_0.ξ,
                Pk = p_o_0.Pk,
                Yk = p_o_0.Yk)
            end
            x0 = [g2.w g2.Yc g2.ξ g2.Pk g2.Yk]

        else

            if s.load_prices ==0  &&  s.PE==0
                const g2=(w = 0.035690209553739,
                ϕ_h = 0.16403,
                ξ =0.369060692554216,
                Pk = 0.975043532756500)
            elseif s.load_prices ==1
                @load "KLS3_Prices.jld2" Prices
                const g2=(w = Prices[2,1],
                ϕ_h = Prices[1,1],
                ξ = Prices[3,1],
                Pk = Prices[4,1])
            elseif s.load_prices ==0 && s.PE==0
                const g2=(w = p_o_0.w,
                ϕ_h = p_o_0.ϕ_h,
                ξ = p_o_0.ξ,
                Pk = p_o_0.Pk)
            end

            x0 = [g2.w g2.ϕ_h g2.ξ g2.Pk]
        end

        @time begin
            if s.GE_end == 1  && s.PE == 0
                results_GE = nlsolve((F,x) ->KLS3_GE_par!(F,x,p,s),log.(x0), autodiff = :forward, method=s.method_GE, xtol=s.xtol_GE, ftol=s.ftol_GE, show_trace=s.show_trace_GE)
                mcc_end, p_o_end, sim_end = KLS3_GE_par(results_GE.zero,p,s)
            else
                mcc_end, p_o_end, sim_end =KLS3_GE_par(log.(x0),p,s)
            end
            say("Final Steady State Finished - you'd better come and take a look....")
        end

        #clearing previously defined global vars
        guessV=nothing
        guessM=nothing
        # Real variables
        sim_end=merge(sim_end,(X_Laspeyres = sum(sim_end.measure.*p_o_0.ξ.*p_o_0.pf.*p_o_end.yf),
        D_Laspeyres = sum(sim_end.measure.* p_o_0.pd.*p_o_end.yd),
        Sales_Laspeyres = sum(sim_end.measure.*(p_o_0.pd.*p_o_end.yd + p_o_0.ξ.*p_o_0.pf.*p_o_end.yf)),
        GDP_Laspeyres = sum(sim_end.measure.*(p_o_0.pd.*p_o_end.yd + p_o_0.ξ*p_o_0.pf.*p_o_end.yf -p_o_end.m*p_o_0.Pk))))

 # STEP 3: guess sequence of aggregate prices for N period, with N sufficiently large [this was in a separate script in matlab, KLS3_trans_obj, you could that here but there isn't much gain and it's useful to mantain a low number of scripts]

 # Store results from initial and final steady states

rt=Vector{Any}(undef,p.N) #this creates a vector that accept every type in its entries (hence the "any")

rt[1] = merge(p_o_0,(measure=sim_0.measure,))
rt[p.N] = merge(p_o_end,(measure=sim_end.measure,))
 #r_0=nothing #[Julia's equivalent to clear]
 #r_end=nothing

 # Guess sequence of aggregate prices for N period, with N sufficiently large
 #ϕht = zeros(p.N,1)
 p_o_t=(Pkt = zeros(p.N,1),
 wt = zeros(p.N,1),
 ξt = zeros(p.N,1),
 Ykt = zeros(p.N,1),
 Yct = zeros(p.N,1))

    if s.transition_AD == 1 # If AutoDiff is enabled, tariffsincomet must be pre-allocated in a way that accepts Real numbers (and not just Float64 as default). This is because when the solver uses AD, it uses a different number type (DualNumbers).
    p_o_t=merge(p_o_t,(tariffsincomet = zeros(Real,p.N,1),))
    else
    p_o_t=merge(p_o_t,(tariffsincomet = zeros(p.N,1),))
    end

    if s.tariffsincome == 1
     # Guess
     #p_o_t.ϕht[1] = p_o_0.ϕ_h
     p_o_t.wt[1] = p_o_0.w
     p_o_t.ξt[1] = p_o_0.ξ
     p_o_t.Pkt[1] = p_o_0.Pk
     p_o_t.Yct[1] = p_o_0.Yc
     p_o_t.Ykt[1] = p_o_0.Yk

        if s.PE == 0
         #p_o_t.ϕht[p.N] = p_o_end.ϕ_h
         p_o_t.wt[p.N] = p_o_end.w
         p_o_t.ξt[p.N] = p_o_end.ξ
         p_o_t.Pkt[p.N] = p_o_end.Pk
         p_o_t.Yct[p.N] = p_o_end.Yc
         p_o_t.Ykt[p.N] = p_o_end.Yk

        elseif s.PE ==1

         #p_o_t.ϕht[p.N]=p_o_0.ϕ_h
         p_o_t.wt[p.N] = p_o_0.w
         p_o_t.ξt[p.N] = p_o_0.ξ
         p_o_t.Pkt[p.N]=p_o_end.Pk #This was in the MATLAB code, is it a typo and should it be p_o_0.Pk as the rest of the prices?
         p_o_t.Yct[p.N] = p_o_0.Yc
         p_o_t.Ykt[p.N] = p_o_0.Yk

        end

     const period2_change = 0.5
     #I labeled every init/expo differently because as I am using "const" (see lines above for an explanation), using the same name for all would be problematic (because it would make them not constant after all)
     const init1 =p_o_t.wt[1]+ (p_o_t.wt[p.N]-p_o_t.wt[1])*period2_change
     const expo1 =log.(p_o_t.wt[p.N]./init1)./log.(p.N)   #p_o_t.wt[1]*(p.N)^expo1=p_o_t.wt[p.N] =>expo1=log.(p_o_t.wt[p.N]/p_o_t.wt[1])/log.(p.N)
     p_o_t.wt[2:p.N-1].= @. init1*(2:p.N-1).^expo1

     const init2 =p_o_t.Yct[1]+ (p_o_t.Yct[p.N]-p_o_t.Yct[1])*period2_change
     const expo2 =log.(p_o_t.Yct[p.N]./init2)./log.(p.N) #p_o_t.wt[1]*(p.N)^expo=p_o_t.wt[p.N] =>expo2=log.(p_o_t.wt[p.N]/p_o_t.wt[1])/log.(p.N)
     p_o_t.Yct[2:p.N-1].= @. init2*(2:p.N-1).^expo2

     const init3 =p_o_t.Ykt[1]+ (p_o_t.Ykt[p.N]-p_o_t.Ykt[1])*period2_change
     const expo3 =log.(p_o_t.Ykt[p.N]./init3)./log.(p.N)   #p_o_t.wt[1]*(p.N).^expo3=p_o_t.wt[p.N] =>expo3=log.(p_o_t.wt[p.N]./p_o_t.wt[1])/log.(p.N)
     p_o_t.Ykt[2:p.N-1].= @. init3*(2:p.N-1).^expo3

     p_o_t.Pkt[2:p.N-1].= p_o_end.Pk
     p_o_t.ξt[2:p.N-1].=p_o_end.ξ

     Guess = [p_o_t.Yct[2:p.N-1] p_o_t.wt[2:p.N-1] p_o_t.ξt[2:p.N-1] p_o_t.Pkt[2:p.N-1] p_o_t.Ykt[2:p.N-1]]

        if s.load_prices == 1
        @load "KLS3_Prices.jld2" Prices
         Guess = [Prices[1,2:p.N-1] Prices[2,2:p.N-1] Prices[3,2:p.N-1] Prices[4,2:p.N-1] Prices[5,2:p.N-1]]
        end
    end

## STEP 4: Solving for the GE prices and quantities

optim_trans(F,x)=KLS3_transition_vec2!(F,x,p,p_o_t,rt) #"instantiating" the objective functions with the previously specified values of the parameters/pre-allocations (in order to use it as an argument for the solvers below)

@time begin
        if s.transition_GE == 1 && s.solver_LM==0
            if s.transition_AD == 1
                #using NLsolve (Trust Region Dogleg Algorithm)
                results_trans = nlsolve(optim_trans,log.(Guess), autodiff = :forward, method=s.method_trans, factor=1, autoscale=false, xtol=s.xtol_trans, ftol=s.ftol_trans,iterations=s.MaxIter_trans, show_trace=s.show_trace_trans,store_trace=true)
            else
                results_trans = nlsolve(optim_trans,log.(Guess), method=s.method_trans, factor=1, autoscale=true, xtol=s.xtol_trans, ftol=s.ftol_trans,iterations=s.MaxIter_trans, show_trace=s.show_trace_trans,store_trace=true)
            end
            Prices_sol=results_trans.zero
            mc, p_o_t, sim_fun, rt = KLS3_transition_vec2(Prices_sol,p,p_o_t,rt)
        elseif s.transition_GE == 1 && s.solver_LM==1
            #using LeastSquaresOptim #(L-M algorithm)
            # #[This solver allows the use of the Levenberg-Marquardt algorithm, which seems to be better at computing the transition (it throws out almost the same solution as matlab's, while NLsolve throws a much noiser one; this changes if the guess is closer to the solution, though). Unfortunately, it takes around 10x the time NLsolve takes] IMPORTANT: If using this solver, uncomment first lines of function KLS3_transition_vec2!

            if s.transition_AD == 1
                solver_LSO=optimize!(LeastSquaresProblem(x = vec(log.(Guess)), f! = optim_trans, output_length = length(vec(log.(Guess))), autodiff = :forward), x_tol=s.xtol_trans, f_tol=s.ftol_trans,iterations=s.MaxIter_trans,LevenbergMarquardt())
            else
                solver_LSO=optimize!(LeastSquaresProblem(x = vec(log.(Guess)), f! = optim_trans, output_length = length(vec(log.(Guess))), autodiff = :forward), x_tol=s.xtol_trans, f_tol=s.ftol_trans,iterations=s.MaxIter_trans,LevenbergMarquardt())
            end

            Prices_sol=hcat(solver_LSO.minimizer[1:p.N-2],solver_LSO.minimizer[p.N-1:2*p.N-4],solver_LSO.minimizer[2*p.N-3:3*p.N-6],solver_LSO.minimizer[3*p.N-5:4*p.N-8],solver_LSO.minimizer[4*p.N-7:5*p.N-10])
            mc_lso, p_o_t_lso, sim_fun_lso, rt_lso = KLS3_transition_vec2(Prices_sol,p,p_o_t,rt)
        else
            Prices_sol = log.(Guess)
            mc, p_o_t, sim_fun, rt = KLS3_transition_vec2(Prices_sol,p,p_o_t,rt)
        end
    say("Transition Finished - you'd better come and take a look....")
end

end

 ## Welfare analysis #INCOMPLETE [3D FIGURES PENDING]
 if s.welfare==1
     include("welfare.jl")
 end

 ## End

 if s.save_prices==1
     if s.tariffsincome == 0
         Prices[1,:] = [m.ϕht[1] exp.(Prices_sol[1:p.N-2])' m.ϕht[p.N]]
         Prices[2,:] = [m.wt[1] exp.(Prices_sol[p.N-1:2*p.N-4])' m.wt[p.N]]
         Prices[3,:] = [m.ξt[1] exp.(Prices_sol[2*p.N-3:3*p.N-6])' m.ξt[p.N]]
         Prices[4,:] = [m.Pkt[1] exp.(Prices_sol[3*p.N-5:4*p.N-8])' m.Pkt[p.N]]
     elseif s.tariffsincome == 1
         Prices[1,:] = [m.Yct[1] exp.(Prices_sol[1:p.N-2])' m.Yct[p.N]]
         Prices[2,:] = [m.wt[1] exp.(Prices_sol[p.N-1:2*p.N-4])' m.wt[p.N]]
         Prices[3,:] = [m.ξt[1] exp.(Prices_sol[2*p.N-3:3*p.N-6])' m.ξt[p.N]]
         Prices[4,:] = [m.Pkt[1] exp.(Prices_sol[3*p.N-5:4*p.N-8])' m.Pkt[p.N]]
         Prices[5,:] = [m.Ykt[1] exp.(Prices_sol[4*p.N-7:end])' m.Ykt[p.N]]
     end
    @save "KLS3_prices.jld2" Prices
 end

 if s.save_workspace==1
     rt = nothing
     r_0 = nothing
     r_end = nothing
     @save "workspace_$(Dates.year(now()))_$(Dates.monthname(now()))_$(Dates.day(now()))_$(Dates.hour(now()))_$(Dates.minute(now()))"

 end

end


text:
cd(@__DIR__); include("setups/rgc100.jl")
plotlyjs(dpi = 200)

gplot(W, X3; grid = true, shape = :circle, mcolor = :yellow)
    plot!(xlab = "x", ylab = "y", zlab = "z", size = (500, 500))

ext: jl


text:
# Example 2.1 of
# A. Ahmadi, and P. Parrilo
# "Joint spectral radius of rank one matrices and the maximum cycle mean problem."
# CDC, 731-733, 2012.
# The JSR is 1

@testset "[AP12] Example 2.1" begin
    for factory in sdp_factories
        tol = 1e-4
        for d in 1:3
            # Get a fresh system to discard lyapunovs, smp and sets lb to 0
            s = discreteswitchedsystem([[0 1; 0 0], [0 0; 1 0]])
            smp = periodicswitching(s, [1, 2])
            @test smp.growthrate == 1
            @test periodicswitching(s, [1, 1]) != periodicswitching(s, [2, 2])

            lb, ub = soslyapb(s, d, factory=factory, tol=tol)
            @test lb ≈ 1 / 2^(1/(2*d)) rtol=tol
            @test ub ≈ 1 rtol=tol

            @test getlb(s) == 1.
            @test hassmp(s)
            @test getsmp(s) == smp

            cyc = sosextractcycle(s, d, tol=tol)
            @test cyc !== nothing
            @test cyc == smp

            seq = sosbuildsequence(s, d, p_0=:Primal)
            psw = findsmp(seq)
            @test psw !== nothing
            @test psw.growthrate == 1
        end
    end
end


text:
using RandomBasedArrays
using Random, Test

@testset "RandomBasedArrays.jl" begin
    Random.seed!(1)
    p = reshape(collect(1:6), 2, 3)
    A = RandomBasedArray(p)
    @test parent(A) === p
    @test parent(A) == [1 3 5; 2 4 6]
    if Sys.WORD_SIZE == 64
        @test size(A) === size(p)
        @test A[42] == 1
        @test A[-35] == 1
        @test A[4,0] == 6
        A[-314] = 0
        @test parent(A) == [1 3 5; 2 0 6]
        A[0,0] = 42
        @test parent(A) == [1 3 5; 42 0 6]
        @test repr("text/plain", A) == "2×3 Array{Int64,2}:\n 3  42  6\n 3   1  5"
    elseif Sys.WORD_SIZE == 32
        @test size(A) === size(p)
        @test A[42] == 3
        @test A[-35] == 3
        @test A[4,0] == 6
        A[-314] = 0
        @test parent(A) == [1 3 5; 2 4 0]
        A[0,0] = 42
        @test parent(A) == [1 3 5; 2 42 0]
        @test repr("text/plain", A) == "2×3 Array{Int32,2}:\n 1  2  42\n 1  1   5"
    else
        @info "We are not on a 32- or 64-bit system"
    end
end


text:
@diffusion_process LotkaVolterra begin
    :dimensions
    process --> 2
    wiener --> 2

    :parameters
    (α, β, γ, δ, σ1, σ2) --> Float64

    :conjugate
    phi(t, x) --> (
        (0.0, 0.0),
        (x[1], 0.0),
        (-x[1]*x[2], 0.0),
        (0.0, x[1]*x[2]),
        (0.0, -x[2]),
        (0.0, 0.0),
        (0.0, 0.0)
    )
    nonhypo(x) --> x
    num_non_hypo --> 2

    :additional
    constdiff --> true
    statespace --> LowerBoundedStateSpace((1,2),(0.0, 0.0))
end

function b(t, x, P::LotkaVolterra)
    ℝ{2}(
        P.α*x[1] - P.β*x[1]*x[2],
        P.δ*x[1]*x[2] - P.γ*x[2],
    )
end

σ(t, x, P::LotkaVolterra) = SDiagonal(P.σ1, P.σ2)


text:
using CalibrateEDMF, Documenter
using DocumenterCitations

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"))

#! format: off
api = [
    "Diagnostics" => "API/Diagnostics.md",
    "ReferenceModels" => "API/ReferenceModels.md",
    "ReferenceStats" => "API/ReferenceStats.md",
    "Pipeline" => "API/Pipeline.md",
    "TurbulenceConvectionUtils" => "API/TurbulenceConvectionUtils.md",
    "KalmanProcessUtils" => "API/KalmanProcessUtils.md",
    "DistributionUtils" => "API/DistributionUtils.md",
    "LESUtils" => "API/LESUtils.md",
    "HelperFuncs" => "API/HelperFuncs.md",
]

pages = Any[
    "Home" => "index.md",
    "Installation instructions" => "installation.md",
    "Getting started" => "getting_started.md",
    "API" => api,
    "Contributing" => "contributing.md",
    "References" => "References.md",
]

mathengine = MathJax(Dict(
    :TeX => Dict(
        :equationNumbers => Dict(:autoNumber => "AMS"),
        :Macros => Dict(),
    ),
))

format = Documenter.HTML(
    prettyurls = get(ENV, "CI", nothing) == "true",
    mathengine = mathengine,
    collapselevel = 1,
)
#! format: on

makedocs(
    bib,
    sitename = "CalibrateEDMF.jl",
    # strict = true, # TODO: make docs strict
    format = format,
    # checkdocs = :exports, # TODO: remove comment when doc strings are added to API
    clean = true,
    doctest = true,
    modules = [CalibrateEDMF],
    pages = pages,
)

deploydocs(
    repo = "github.com/CliMA/CalibrateEDMF.jl.git",
    target = "build",
    push_preview = true,
    devbranch = "main",
    forcepush = true,
)


text:
using MultilayerPerceptron
using Random: shuffle

using Makie

function sine_example()
    N = 2000

    train_x = (0.0:0.1:(2π + 0.1))' |> collect |> shuffle
    train_y = sin.(train_x)

    hidden_layer1 = Layer((5, 1); activation_fn=logistic, grad_fn=∇logistic)
    hidden_layer2 = Layer((3, 5); activation_fn=logistic, grad_fn=∇logistic)
    output_layer = Layer((1, 3))

    mlp = MLP([hidden_layer1, hidden_layer2, output_layer])
    fit!(mlp, train_x, train_y; learning_rate=0.5, batch_size=5, n_epochs=N)

    test_x = (0.0:0.1:(2π + 0.1))' |> collect
    test_y = sin.(test_x)

    # Make PackageCompiler-precompiled Makie show the window.
    # See https://github.com/JuliaPlots/Makie.jl/issues/300.
    AbstractPlotting.__init__()

    scene = Scene()
    lines!(scene, test_x[:], test_y[:]; color=:blue)
    lines!(scene, test_x[:], predict(mlp, Float64.(test_x))[:]; color=:red)
    scene_legend = legend(scene.plots[2:end], ["sin(x)", "ŷ"])
    vbox(scene, scene_legend)
end

sine_example()

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
"""
Start a notebook server using:

```julia
julia> Pluto.run()
```

Have a look at the FAQ:
https://github.com/fonsp/Pluto.jl/wiki
"""
module Pluto

import Pkg

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))
include_dependency(joinpath(PKG_ROOT_DIR, "Project.toml"))
const PLUTO_VERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(PKG_ROOT_DIR, "Project.toml"))["version"])
const PLUTO_VERSION_STR = 'v' * string(PLUTO_VERSION)
const JULIA_VERSION_STR = 'v' * string(VERSION)
const ENV_DEFAULTS = Dict(
    "PLUTO_WORKSPACE_USE_DISTRIBUTED" => "true",
    "PLUTO_RUN_NOTEBOOK_ON_LOAD" => "true",
    "PLUTO_WORKING_DIRECTORY" => let
        preferred_dir = startswith(Sys.BINDIR, pwd()) ? homedir() : pwd()
        joinpath(preferred_dir, "") # must end with / or \
    end,
)
get_pl_env(key::String) = haskey(ENV, key) ? ENV[key] : ENV_DEFAULTS[key]

if get(ENV, "PLUTO_SHOW_BANNER", "true") == "true"
@info """\n
    Welcome to Pluto $(PLUTO_VERSION_STR) 🎈
    Start a notebook server using:

  julia> Pluto.run()

    Have a look at the FAQ:
    https://github.com/fonsp/Pluto.jl/wiki
\n"""
end

include("./evaluation/Tokens.jl")
include("./runner/PlutoRunner.jl")
include("./analysis/ExpressionExplorer.jl")

include("./notebook/PathHelpers.jl")
include("./notebook/Cell.jl")
include("./notebook/Notebook.jl")
include("./webserver/Session.jl")

include("./analysis/Errors.jl")
include("./analysis/Parse.jl")
include("./analysis/Topology.jl")

include("./evaluation/WorkspaceManager.jl")
include("./evaluation/Update.jl")
include("./evaluation/Run.jl")

include("./webserver/Configuration.jl")
include("./webserver/MsgPack.jl")
include("./webserver/PutUpdates.jl")
include("./webserver/SessionActions.jl")
include("./webserver/Static.jl")
include("./webserver/Dynamic.jl")
include("./webserver/REPLTools.jl")
include("./webserver/WebServer.jl")

end


text:
#=
The following is a Julia 1.4.2 translation of the core functionality of R lanague packages
https://cran.r-project.org/web/packages/InvariantCausalPrediction/InvariantCausalPrediction.pdf
https://cran.r-project.org/web/packages/nonlinearICP/nonlinearICP.pdf

"In practice, first apply ICP with linear models.  Apply a nonlinear version if all linear
models are rejected by linear IPC."

Two Jupyter labs showcases the R language version of ICP.
1. Linear Invariant Causal Prediction Using Employment Data From The Work Bank
https://notes.quantecon.org/submission/5e851bfecc00b7001acde469

2.  Nonlinear Invariant Causal Prediction Using Unemployment Data and Inflation Adjusted Prices from the USA Bureau of Labor
https://notes.quantecon.org/submission/5e8e2a6cd079ab001915ca09

New contributions here in Julia are improvements to code readability and support ability.
Also, improvements to program speed such parallelism of random forest computations, p-value linear computations.
There are two version of the InvariantCausalPrediction main functions. One version is sequential and the other parallel.
In addition, there are new VegaLite plots of the InvariantCausalPrediction results.

Cheers,
Clarman Cruz
June 2020 =#

println("Loading Libraries... This can take a few minutes")
using Distributed;
using  Hwloc;                               println("Hwloc")
if length(workers()) == 1
    addprocs(Hwloc.num_physical_cores())  # get ready for parallelism
end
println("Distributed with ", length(workers()), " workers")
using DataFrames;                           println("DataFrames")
using Distributions: Normal, quantile;      println("Distributions")
using Query: @rename, @orderby;             println("Query")
using VegaLite;                             println("VegaLite")
include("GetPValue.jl");                    println("GetPValue.jl")
include("GetBlankets.jl");                  println("GetBlankets.jl")
using .GetBlankets:  getBlanketAll, getBlanketBoosting, getBlanketRandomForest;
include("InvariantEnvironmentPrediction.jl");
using .InvariantEnvironmentPredictionModule; println("InvariantEnvironmentPrediction.jl")
include("ICPlibrary.jl");                    println("ICPlibrary.jl")
println("Libraries Loaded")

"""
boxPlotICPconfidenceIntervals(X, ConfInt)
Creates a box plot of the `ConfInt` created from 'X'
# Arguments
* `X`:                 as returned by InvariantCausalPrediction function
* `ConfInt`:           as returned by InvariantCausalPrediction function
* `width`:             optial plot width
* `height`:            optial plot height
"""
function boxPlotICPconfidenceIntervals(
    X::DataFrame,
    ConfInt::Array{Float64,2},
    width=600,
    height=500,
)
    if ConfInt === missing
        println("boxPlotICPconfidenceIntervals:  There are not any Confidence Intervals")
        return (missing)
    end

    dfGraph = DataFrame(ConfInt, names(X))
    dfTranspose = DataFrames.stack(dfGraph, 1:ncol(dfGraph))  |> @rename(:value => :ConfidenceInterval) |> DataFrame

    plot = dfTranspose |>
        @vlplot(
            mark = {:boxplot, extent = "min-max"},
            x = {"variable:o", axis = {title = "Predictors"}},
            y = {:ConfidenceInterval, axis = {title = "Invariant Causal Prediction"}},
            size = {value = 25},
            width  = width,
            height = height)

    VegaLite.display(plot)
end # boxPlotICPconfidenceIntervals

"""
plotICPacceptedSets(X, acceptedSets)
Creates a circle plot of the predictor count in the 'acceptedSets' created from 'X'
The larger the circle, the strong the predictor
# Arguments
* `X`:                 as returned by InvariantCausalPrediction function
* `acceptedSets`:      as returned by InvariantCausalPrediction function
* `width`:             optial plot width
* `height`:            optial plot height
"""
function plotICPacceptedSets(
    X::DataFrame,
    acceptedSets::Array{Array{Int64,1},1},
    width=600,
    height=500,
)
    hash = Dict{Int64,Int64}()
    for s in acceptedSets
        for e in collect(s)
            if haskey(hash, Int(e))
                hash[Int(e)] += 1
            else
                hash[Int(e)] = 1
            end
        end
    end

    dfGraph = DataFrames.stack(DataFrame(hash), 1:length(keys(hash))) |> @orderby(_.variable) |> DataFrame
    predictorFromIndex = getPredictorFromIndex(X)
    dfGraph.predictor = [predictorFromIndex[i] for i in keys(hash)]

    plot = dfGraph |>
    @vlplot(
        :circle,
        width = width,
        height = height,
        color = { "predictor:n", legend = nothing},
        y = {"value:o", axis = {title = "Count in Accepted Sets"}},
        x = "predictor:o",
        size = :value
    )
    VegaLite.display(plot)
end # plotICPacceptedSets

"""
barPlotICPCasualPredictors(X, casualPredictors)
Creates a bar plot of the 'casualPredictors' created from 'X'
# Arguments
* `X`:                 as returned by InvariantCausalPrediction function
* `casualPredictors`:  as returned by InvariantCausalPrediction function
* `width`:             optial plot width
* `height`:            optial plot height
"""
function barPlotICPCasualPredictors(
    X::DataFrame,
    casualPredictors::Union{Dict{Int64,String},Missing},
    width=600,
    height=500,
)
    if casualPredictors === missing
        println("barPlotICPCasualPredictors:  There are not any Casual Predictors")
        return (missing)
    end

    xAxis = [String(n) for n in names(X)]
    yAxis = [n in keys(casualPredictors) ? 1 : 0 for n in 1:ncol(X)]
    plot = @vlplot(:bar,
        width = width,
        height = height,
        y = {yAxis, axis = {title = "In All Accepted Sets"}},
        x = {xAxis, axis = {title = "Predictor"}}
    )
    VegaLite.display(plot)
end # barPlotICPCasualPredictors

"""
LinearInvariantCausalPrediction!(X, Y, ExpInd; α=0.01, selection=booster, verbose=true)
Searching in parallel over subsets in `X[,S]` for direct causes of `Y`
# Arguments
* `X`:                 DataFrame of the predictors. No missing values are allowed
* `Y`:                 DataFrame of the target.  Binary for classification or floats for regression. No missing values are allowed
* `ExpInd`:            DataFrame of environment indicators (e.g. 1, 2, ...)
* `α`:                 significance level (e.g. 0.01 or 0.05)
* `selection`:         if all, will include all the predictors else XGBooster selects the predictors
* `verbose`:           if true, will print each subset tested and some other output
"""
function LinearInvariantCausalPrediction!(
    X::DataFrame,
    Y::DataFrame,
    ExpInd::DataFrame;
    α=0.01,
    selection="booster",
    verbose=false,
)
    @assert nrow(X) == nrow(Y) == nrow(ExpInd)
    @assert ncol(X) >= 2
    @assert ncol(Y) == 1
    @assert ncol(ExpInd) == 1

    X = hotEncoder(X)
    numberOfTargets = length(unique(Y[:, 1]))
    numberOfEnvironments = nrow(unique(ExpInd))
    @assert numberOfEnvironments >= 2
    variableUniverse = Set(1:ncol(X))

    println("Experiment Indicator Integer Mapper")
    categoryIntegerMapper!(ExpInd)

    acceptedSets = [Int64[]]
    ConfInt = zeros(2, ncol(X))
    ConfInt[1, :] .= -Inf
    ConfInt[2, :] .= Inf
    Coeff = Dict{Int64,Array{Float64}}()
    CoeffVar = Dict{Int64,Array{Float64}}()
    PValuesAccepted = Dict{Int64,Array{Float64}}()
    PValuesNotAccepted = Dict{Int64,Array{Float64}}()
    Pall = Float64[]
    gof = max(0.01, α)

    if numberOfTargets > 2
        println("More than two targests in Y thus doing Linear Regression")
    else
        println("Two targests in Y thus doing Logistic Classification")
        println("Y Integer Mapper")
        categoryIntegerMapper!(Y,0)
        @assert maximum(Y[:,1]) == 1
    end
    fY = float.(Y[:, 1])

    (blankets, featureImportances) = getBlanketBoosting(X, Y, selection)
	pValues = pmap(testSet -> getPValueLinear(
            DataFrame(X[:, testSet]),
            fY,
            ExpInd,
            numberOfTargets,
            numberOfEnvironments
			), blankets, retry_delays=zeros(5))
    @assert length(pValues) == length(blankets)

	testSet = 1
    for r in pValues
		notUsedVariables = sort(collect(setdiff(variableUniverse, blankets[testSet])))
        push!(Pall, r.pvalue)

        if r.pvalue > α
            push!(acceptedSets, blankets[testSet])
            v = getQuantileNormalStd(α, r.coefsStdError)
            ConfInt[1, blankets[testSet]] = max(ConfInt[1, blankets[testSet]], r.coefs .+ v)
            ConfInt[2, blankets[testSet]] = min(ConfInt[2, blankets[testSet]], r.coefs .- v)

            if length(notUsedVariables) >= 1
                ConfInt[1, notUsedVariables] = max(
                    ConfInt[1, notUsedVariables],
                    zeros(length(ConfInt[1, notUsedVariables])),
                )
                ConfInt[2, notUsedVariables] = min(
                    ConfInt[2, notUsedVariables],
                    zeros(length(ConfInt[2, notUsedVariables])),
                )
            end # of notUsedVariables

            StoreAcceptedOuput!(r, blankets[testSet], PValuesAccepted, Coeff, CoeffVar)
        else
            StoreNotAcceptedOuput!(r, blankets[testSet], PValuesNotAccepted)
        end  # of acceptedSets

        if verbose == true
            println("TestSet:  ", blankets[testSet], "\nResults:  ", r)
            println("ConfInt ", ConfInt)
            println("Pvalues All: ", Pall)
        end

		testSet += 1
	end # of testSets

    casualPredictors = doOutput( X, featureImportances, acceptedSets, PValuesAccepted,
                                 PValuesNotAccepted, Pall, gof, verbose, ConfInt, Coeff, CoeffVar,)

    return((X = X, acceptedSets = acceptedSets, ConfInt = ConfInt, Coeff = Coeff, CoeffVar = CoeffVar,
                PValuesAccepted = PValuesAccepted, PValuesNotAccepted = PValuesNotAccepted, Pall = Pall, casualPredictors = casualPredictors))
end # LinearInvariantCausalPrediction

"""
ForestInvariantCausalPrediction!(X, Y, ExpInd; α=0.01, selection=forest, verbose=true)
Searching in parallel over subsets in `X[,S]` for direct causes of `Y`
# Arguments
* `X`:                 DataFrame of the predictors. No missing values are allowed
* `Y`:                 DataFrame of the target.  Binary for classification or floats for regression. No missing values are allowed
* `ExpInd`:            DataFrame of environment indicators (e.g. 1, 2, ...)
* `α`:                 significance level (e.g. 0.01 or 0.05)
* `selection`:         if all, will include all the predictors.  If booster, XGBooster selects the predictors.  Else random forest with Shapley selects the predictors
* `verbose`:           if true, will print each subset tested and some other output
"""
function ForestInvariantCausalPrediction!(
    X::DataFrame,
    Y::DataFrame,
    ExpInd::DataFrame;
    α=0.05,
    selection="forest",
    verbose=false,
)
    @assert nrow(X) == nrow(Y) == nrow(ExpInd)
    @assert ncol(X) >= 1
    @assert ncol(Y) == 1
    @assert ncol(ExpInd) == 1

    X = labelEncoder(X)
    numberOfTargets = length(unique(Y[:, 1]))
    numberOfEnvironments = nrow(unique(ExpInd))
    @assert numberOfEnvironments >= 2
    variableUniverse = Set(1:ncol(X))

    println("Experiment Indicator Integer Mapper")
    categoryIntegerMapper!(ExpInd, 0)
    cExpInd = ExpInd[:,1] |> categorical

    acceptedSets = [Int64[]]
    PValuesAccepted = Dict{Int64,Array{Float64}}()
    PValuesNotAccepted = Dict{Int64,Array{Float64}}()
    Pall = Float64[]
    gof = max(0.01, α)

    if numberOfTargets > 2
        println("More than two targests in Y thus doing Random Forest Regression")
    else
        println("Two targests in Y thus doing Random Forest Classification")
        println("Y Integer Mapper")
        categoryIntegerMapper!(Y,0)
        @assert maximum(Y[:,1]) == 1
    end
    fY = float.(Y[:, 1])

    (blankets, featureImportances) = getBlanketRandomForest(X, Y, selection)
    pValues = pmap(testSet -> InvariantEnvironmentPrediction(
				DataFrame(X[:, testSet]),
				fY,
				cExpInd,
				numberOfTargets,
				numberOfEnvironments
			), blankets, retry_delays=zeros(5))
    @assert length(pValues) == length(blankets)

    testSet = 1
    for r in pValues
        push!(Pall, r.pvalue)
        if r.pvalue > α
            push!(acceptedSets, blankets[testSet])
            StoreAcceptedOuput!(r, blankets[testSet], PValuesAccepted)
        else
            StoreNotAcceptedOuput!(r, blankets[testSet], PValuesNotAccepted)
        end # of acceptedSets
        testSet += 1
    end # of testSets

    casualPredictors = doOutput(X, featureImportances, acceptedSets, PValuesAccepted, PValuesNotAccepted, Pall, gof, verbose)

    return((X = X, acceptedSets = acceptedSets, PValuesAccepted = PValuesAccepted, PValuesNotAccepted = PValuesNotAccepted,
			Pall = Pall, casualPredictors = casualPredictors, featureImportances = featureImportances))
end # ForestInvariantCausalPrediction

"""
LinearInvariantCausalPredictionSequential!(X, Y, ExpInd; α=0.01, selection=booster, verbose=true)
Searching over subsets in `X[,S]` for direct causes of `Y`
# Arguments
* `X`:                 DataFrame of the predictors. No missing values are allowed
* `Y`:                 DataFrame of the target.  Binary for classification or floats for regression. No missing values are allowed
* `ExpInd`:            DataFrame of environment indicators (e.g. 1, 2, ...)
* `α`:                 significance level (e.g. 0.01 or 0.05)
* `selection`:         if all, will include all the predictors else XGBooster selects the predictors
* `verbose`:           if true, will print each subset tested and some other output
"""
function LinearInvariantCausalPredictionSequential!(
    X::DataFrame,
    Y::DataFrame,
    ExpInd::DataFrame;
    α=0.01,
    selection="booster",
    verbose=false,
)
    @assert nrow(X) == nrow(Y) == nrow(ExpInd)
    @assert ncol(X) >= 2
    @assert ncol(Y) == 1
    @assert ncol(ExpInd) == 1

    X = hotEncoder(X)
    numberOfTargets = length(unique(Y[:, 1]))
    numberOfEnvironments = nrow(unique(ExpInd))
    @assert numberOfEnvironments >= 2
    variableUniverse = Set(1:ncol(X))

    println("Experiment Indicator Integer Mapper")
    categoryIntegerMapper!(ExpInd)

    acceptedSets = [Int64[]]
    ConfInt = zeros(2, ncol(X))
    ConfInt[1, :] .= -Inf
    ConfInt[2, :] .= Inf
    Coeff = Dict{Int64,Array{Float64}}()
    CoeffVar = Dict{Int64,Array{Float64}}()
    PValuesAccepted = Dict{Int64,Array{Float64}}()
    PValuesNotAccepted = Dict{Int64,Array{Float64}}()
    Pall = Float64[]
    gof = max(0.01, α)

    if numberOfTargets > 2
        println("More than two targests in Y thus doing Linear Regression")
    else
        println("Two targests in Y thus doing Logistic Classification")
        println("Y Integer Mapper")
        categoryIntegerMapper!(Y,0)
        @assert maximum(Y[:,1]) == 1
    end
    fY = float.(Y[:, 1])

    (blankets, featureImportances) = getBlanketBoosting(X, Y, selection)
    for testSet in blankets
        notUsedVariables = sort(collect(setdiff(variableUniverse, testSet)))

        r = getPValueLinear(
            DataFrame(X[:, testSet]),
            fY,
            ExpInd,
            numberOfTargets,
            numberOfEnvironments
        )
        push!(Pall, r.pvalue)
        if r.pvalue > α
            push!(acceptedSets, testSet)
            v = getQuantileNormalStd(α, r.coefsStdError)
            ConfInt[1, testSet] = max(ConfInt[1, testSet], r.coefs .+ v)
            ConfInt[2, testSet] = min(ConfInt[2, testSet], r.coefs .- v)

            if length(notUsedVariables) >= 1
                ConfInt[1, notUsedVariables] = max(
                    ConfInt[1, notUsedVariables],
                    zeros(length(ConfInt[1, notUsedVariables])),
                )
                ConfInt[2, notUsedVariables] = min(
                    ConfInt[2, notUsedVariables],
                    zeros(length(ConfInt[2, notUsedVariables])),
                )
            end # of notUsedVariables

            StoreAcceptedOuput!(r, testSet, PValuesAccepted, Coeff, CoeffVar)
        else
            StoreNotAcceptedOuput!(r, testSet, PValuesNotAccepted)
        end # of testSets

        if verbose == true
            println("TestSet:  ", testSet, "\nResults:  ", r)
            println("ConfInt ", ConfInt)
            println("Pvalues All: ", Pall)
        end
    end # of getBlanketAll

    casualPredictors = doOutput( X, featureImportances, acceptedSets, PValuesAccepted,
                                 PValuesNotAccepted, Pall, gof, verbose, ConfInt, Coeff, CoeffVar,)

    return((X = X, acceptedSets = acceptedSets, ConfInt = ConfInt, Coeff = Coeff, CoeffVar = CoeffVar,
                PValuesAccepted = PValuesAccepted, PValuesNotAccepted = PValuesNotAccepted, Pall = Pall, casualPredictors = casualPredictors))
end # LinearInvariantCausalPredictionSequential

"""
ForestInvariantCausalPredictionSequential!(X, Y, ExpInd; α=0.01, selection=forest, verbose=true)
Searching over subsets in `X[,S]` for direct causes of `Y`
# Arguments
* `X`:                 DataFrame of the predictors. No missing values are allowed
* `Y`:                 DataFrame of the target.  Binary for classification or floats for regression. No missing values are allowed
* `ExpInd`:            DataFrame of environment indicators (e.g. 1, 2, ...)
* `α`:                 significance level (e.g. 0.01 or 0.05)
* `selection`:         if all, will include all the predictors.  If booster, XGBooster selects the predictors.  Else random forest with Shapley selects the predictors
* `verbose`:           if true, will print each subset tested and some other output
"""
function ForestInvariantCausalPredictionSequential!(
    X::DataFrame,
    Y::DataFrame,
    ExpInd::DataFrame;
    α=0.05,
    selection="forest",
    verbose=false,
)
    @assert nrow(X) == nrow(Y) == nrow(ExpInd)
    @assert ncol(X) >= 1
    @assert ncol(Y) == 1
    @assert ncol(ExpInd) == 1

    X = labelEncoder(X)
    numberOfTargets = length(unique(Y[:, 1]))
    numberOfEnvironments = nrow(unique(ExpInd))
    @assert numberOfEnvironments >= 2
    variableUniverse = Set(1:ncol(X))

    println("Experiment Indicator Integer Mapper")
    categoryIntegerMapper!(ExpInd, 0)
    cExpInd = ExpInd[:,1] |> categorical

    acceptedSets = [Int64[]]
    PValuesAccepted = Dict{Int64,Array{Float64}}()
    PValuesNotAccepted = Dict{Int64,Array{Float64}}()
    Pall = Float64[]
    gof = max(0.01, α)

    if numberOfTargets > 2
        println("More than two targests in Y thus doing Random Forest Regression")
    else
        println("Two targests in Y thus doing Random Forest Classification")
        println("Y Integer Mapper")
        categoryIntegerMapper!(Y,0)
        @assert maximum(Y[:,1]) == 1
    end
    fY = float.(Y[:, 1])

    (blankets, featureImportances) = getBlanketRandomForest(X, Y, selection)
    for testSet in blankets
        r = InvariantEnvironmentPrediction(
            DataFrame(X[:, testSet]),
            fY,
            cExpInd,
            numberOfTargets,
            numberOfEnvironments
        )
        push!(Pall, r.pvalue)
        if r.pvalue > α
            push!(acceptedSets, testSet)
            StoreAcceptedOuput!(r, testSet, PValuesAccepted)
        else
            StoreNotAcceptedOuput!(r, testSet, PValuesNotAccepted)
        end # of acceptedSets

        if verbose == true
            println("TestSet:  ", testSet, "\nResults:  ", r)
            println("Pvalues All: ", Pall)
        end
    end # of testSets

    casualPredictors = doOutput(X, featureImportances, acceptedSets, PValuesAccepted, PValuesNotAccepted, Pall, gof, verbose)

    return((X = X, acceptedSets = acceptedSets, PValuesAccepted = PValuesAccepted, PValuesNotAccepted = PValuesNotAccepted,
			Pall = Pall, casualPredictors = casualPredictors, featureImportances = featureImportances))
end # ForestInvariantCausalPredictionSequential


text:
# Miscellaneous functions

"""
    logmvgamma(p::Int, a::Real)

Return the logarithm of [multivariate gamma function](https://en.wikipedia.org/wiki/Multivariate_gamma_function) ([DLMF 35.3.1](https://dlmf.nist.gov/35.3.1)).
"""
function logmvgamma(p::Int, a::Real)
    # NOTE: one(a) factors are here to prevent unnecessary promotion of Float32
    res = p * (p - 1) * log(pi * one(a)) / 4
    for ii in 1:p
        res += lgamma(a + (1 - ii) * one(a)/ 2)
    end
    return res
end

"""
    lstirling_asym(x)

The remainder term after
[Stirling's approximation](https://en.wikipedia.org/wiki/Stirling%27s_approximation)
to [`lgamma`](@ref):

```math
\\log \\Gamma(x) \\approx x \\log(x) - x + \\log(2\\pi/x)/2 = \\log(x)*(x-1/2) + \\log(2\\pi)/2 - x
```

In Julia syntax, this means:

    lstirling_asym(x) = lgamma(x) + x - (x-0.5)*log(x) - 0.5*log(2π)

For sufficiently large `x`, this can be approximated using the asymptotic
_Stirling's series_ ([DLMF 5.11.1](https://dlmf.nist.gov/5.11.1)):

```math
\\frac{1}{12x} - \\frac{1}{360x^3} + \\frac{1}{1260x^5} - \\frac{1}{1680x^7} + \\ldots
```

The truncation error is bounded by the first omitted term, and is of the same sign.

Relative error of approximation is bounded by
    (174611/125400 x^-19) / (1/12 x^-1 - 1/360 x^-3)
which is < 1/2 ulp for x >= 10.0, and total numeric error appears to be < 2 ulps

# References

* Temme, N. (1996) Special functions: An introduction to the classical functions of
   mathematical physics, Wiley, New York, ISBN: 0-471-11313-1, Chapter 3.6, pp 61-65.
* Weisstein, Eric W. ["Stirling's Series."](http://mathworld.wolfram.com/StirlingsSeries.html).
  MathWorld.
* [OEIS A046968](http://oeis.org/A046968) and [OEIS A046969](http://oeis.org/A046969)
  for the series coefficients
"""
function lstirling_asym end

lstirling_asym(x::BigFloat) = lgamma(x) + x - log(x)*(x - big(0.5)) - log2π/big(2)

lstirling_asym(x::Integer) = lstirling_asym(float(x))

const lstirlingF64 = Float64[lstirling_asym(k) for k in big(1):big(64)]
const lstirlingF32 = Float64[lstirling_asym(k) for k in big(1):big(40)]

function lstirling_asym(x::Float64)
    isinteger(x) && (0 < x ≤ length(lstirlingF64)) && return lstirlingF64[Int(x)]
    t = inv(abs2(x))
    @horner(t,
             8.33333333333333333e-2, #  1/12 x^-1
            -2.77777777777777778e-3, # -1/360 x^-3
             7.93650793650793651e-4, #  1/1260 x^-5
            -5.95238095238095238e-4, # -1/1680 x^-7
             8.41750841750841751e-4, #  1/1188 x^-9
            -1.91752691752691753e-3, # -691/360360 x^-11
             6.41025641025641026e-3, #  1/156 x^-13
            -2.95506535947712418e-2, # -3617/122400 x^-15
             1.79644372368830573e-1)/x #  43867/244188 x^-17
end

function lstirling_asym(x::Float32)
    isinteger(x) && (0 < x ≤ length(lstirlingF32)) && return lstirlingF32[Int(x)]
    t = inv(abs2(x))
    @horner(t,
             8.333333333333f-2, #  1/12 x^-1
            -2.777777777777f-3, # -1/360 x^-3
             7.936507936508f-4, #  1/1260 x^-5
            -5.952380952381f-4, # -1/1680 x^-7
             8.417508417508f-4)/x #  1/1188 x^-9
end


text:
############################################################################################
#=
!NOTE 1: flatten functions taken and adjusted from excellent package ParameterHandling.jl, see:
    https://github.com/invenia/ParameterHandling.jl/blob/8e998949e7fcf24d5c3f8bed5018ec300542151b/src/flatten.jl#LL1-L17
!NOTE 2: Functions redefined from discussion in:
    https://github.com/invenia/ParameterHandling.jl/issues/27
=#

############################################################################################
#!NOTE: I am very open to change type/struct names in this section :)

"""
$(TYPEDEF)

Supertype for dispatching different types of flatten. Determines if all inputs are flattened (FlattenAll) or only continuous values (FlattenContinuous).

# Fields
$(TYPEDFIELDS)
"""
abstract type FlattenTypes end
struct FlattenAll <: FlattenTypes end
struct FlattenContinuous <: FlattenTypes end

"""
$(TYPEDEF)

Determines if unflatten returns original type or if type may change (AD friendly).

# Fields
$(TYPEDFIELDS)
"""
abstract type UnflattenTypes end
struct UnflattenStrict <: UnflattenTypes end
struct UnflattenFlexible <: UnflattenTypes end

############################################################################################
"""
$(TYPEDEF)

Default arguments for flatten function.

# Fields
$(TYPEDFIELDS)
"""
struct FlattenDefault{T<:AbstractFloat,F<:FlattenTypes}
    "Type of flatten output"
    output::Type{T}
    "Determines if all inputs are flattened (FlattenAll) or only continuous values (FlattenContinuous)."
    flattentype::F
end
function FlattenDefault(;
    output::Type{T}=Float64,
    flattentype::F=FlattenContinuous(),
) where {T,F<:FlattenTypes}
    return FlattenDefault(output, flattentype)
end

############################################################################################
"""
    $(FUNCTIONNAME)(x )
Construct a flatten function for 'x' given specifications in 'df'.

# Examples
```julia
```
"""
function construct_flatten end
construct_flatten(x) = construct_flatten(FlattenDefault(), UnflattenStrict(), x)
construct_flatten(constraint, x) = construct_flatten(FlattenDefault(), UnflattenStrict(), constraint, x)

construct_flatten(df::F, x) where {F<:FlattenDefault} = construct_flatten(df.output, df.flattentype, UnflattenStrict(), x)
construct_flatten(df::F, constrained, x) where {F<:FlattenDefault} = construct_flatten(df.output, df.flattentype, UnflattenStrict(), constrained, x)
construct_flatten(df::F, unflattentype::U, x) where {F<:FlattenDefault, U<:UnflattenTypes} = construct_flatten(df.output, df.flattentype, unflattentype, x)
construct_flatten(df::F, unflattentype::U, constraint, x) where {F<:FlattenDefault, U<:UnflattenTypes} = construct_flatten(df.output, df.flattentype, unflattentype, constraint, x)

############################################################################################
"""
$(TYPEDEF)

Contains information for flatten construct.

# Fields
$(TYPEDFIELDS)
"""
struct FlattenConstructor{S<:Function, T<:Function}
    strict::S
    flexible::T
    function FlattenConstructor(strict::S, flexible::T) where {S<:Function, T<:Function}
        return new{S,T}(strict, flexible)
    end
end

############################################################################################
"""
$(TYPEDEF)

Contains information for unflatten construct.

# Fields
$(TYPEDFIELDS)
"""
struct UnflattenConstructor{S<:Function, T<:Function}
    strict::S
    flexible::T
    function UnflattenConstructor(strict::S, flexible::T) where {S<:Function, T<:Function}
        return new{S,T}(strict, flexible)
    end
end

############################################################################################
"""
$(TYPEDEF)

Contains information for flatten/unflatten construct.

# Fields
$(TYPEDFIELDS)
"""
struct ReConstructor{F<:FlattenDefault, S<:FlattenConstructor, T<:UnflattenConstructor}
    default::F
    flatten::S
    unflatten::T
    function ReConstructor(default::F, flatten::S, unflatten::T) where {F<:FlattenDefault, S<:FlattenConstructor, T<:UnflattenConstructor}
        return new{F,S,T}(default, flatten, unflatten)
    end
end
function ReConstructor(flattendefault::FlattenDefault, x)
# Assign flatten constructors
    flatten, unflatten = construct_flatten(flattendefault, UnflattenStrict(), x)
    flattenAD, unflattenAD = construct_flatten(flattendefault, UnflattenFlexible(), x)
    flatten_constructor = FlattenConstructor(flatten, flattenAD)
    unflatten_constructor = UnflattenConstructor(unflatten, unflattenAD)
# Return structs
    return ReConstructor(flattendefault, flatten_constructor, unflatten_constructor)
end
function ReConstructor(flattendefault::FlattenDefault, constraint, x)
# Assign flatten constructors
    flatten, unflatten = construct_flatten(flattendefault, UnflattenStrict(), constraint, x)
    flattenAD, unflattenAD = construct_flatten(flattendefault, UnflattenFlexible(), constraint, x)
    flatten_constructor = FlattenConstructor(flatten, flattenAD)
    unflatten_constructor = UnflattenConstructor(unflatten, unflattenAD)
# Return structs
    return ReConstructor(flattendefault, flatten_constructor, unflatten_constructor)
end
function ReConstructor(x)
    return ReConstructor(FlattenDefault(), x)
end
function ReConstructor(constraint, x)
    return ReConstructor(FlattenDefault(), constraint, x)
end

############################################################################################
"""
    $(FUNCTIONNAME)(x )
Convert 'x' into a Vector.

# Examples
```julia
```
"""
function flatten end
function flatten(constructor::ReConstructor, x)
    return constructor.flatten.strict(x)
end

"""
    $(FUNCTIONNAME)(x )
Convert 'x' into a Vector that is AD compatible.

# Examples
```julia
```
"""
function flattenAD end
function flattenAD(constructor::ReConstructor, x)
    return constructor.flatten.flexible(x)
end

"""
    $(FUNCTIONNAME)(x )
Unflatten 'x' into original shape.

# Examples
```julia
```
"""
function unflatten end
function unflatten(constructor::ReConstructor, x)
    return constructor.unflatten.strict(x)
end

"""
    $(FUNCTIONNAME)(x )
Unflatten 'x' into original shape but keep type information of 'x' for AD compatibility.

# Examples
```julia
```
"""
function unflattenAD end
function unflattenAD(constructor::ReConstructor, x)
    return constructor.unflatten.flexible(x)
end

############################################################################################
"""
    $(FUNCTIONNAME)(x )
Contains information to constrain and unconstrain parameter.

# Examples
```julia
```
"""
struct TransformConstructor{S, T}
    constrain::S
    unconstrain::T
    function TransformConstructor(constraint, x)
        #!NOTE: Transform is used to unconstrain, and inverse-transform to constrain parameter back.
        transform, inverse_transform = construct_transform(constraint, x)
        return new{typeof(inverse_transform), typeof(transform)}(inverse_transform, transform)
    end
end

function constrain(transform::T, val) where {T<:TransformConstructor}
    return constrain(transform.constrain, val)
end
function unconstrain(transform::T, val) where {T<:TransformConstructor}
    return unconstrain(transform.unconstrain, val)
end

############################################################################################
# Export
export FlattenTypes,
    FlattenAll,
    FlattenContinuous,
    UnflattenTypes,
    UnflattenStrict,
    UnflattenFlexible,
    FlattenDefault,
    FlattenConstructor,
    UnflattenConstructor,
    ReConstructor,
    flatten,
    flattenAD,
    unflatten,
    unflattenAD,
    TransformConstructor


text:
"""
`module Electrostatics`

This module implements various functions related to (long range) electrostatic interactions
"""

module Electrostatics

using LinearAlgebra

export get_dipole, electrostatic_energy, soft_coulomb, soft_q_μ, soft_μ_μ

const μ_0 = 4.0e-7 * π
const c = 299792458.0
const ϵ_0 = 1 / μ_0 / c^2
const e = 1.602176634e-19
const c_light = 299792458

"""
Get the overall dipole moment of a set of oint charges and point dipoles
Arguments:
  pos::Array natoms x 3 Array
  charges::Array natoms x 1 Array
  dipoles::Array natoms x 3 Array
Returns:
  total_dipole::Array 1 x 3 Array
"""
function get_dipole(pos::Array, charges::Array, dipoles::Array, pbc::Bool=false)
  @assert pbc == false "Periodic boundary condition not yet supported"
  return sum(pos .* charges .* (1e-11/c_light/e) .+ dipoles, dims=1)
end

"""
Total electrostatic enenrgy of a set of point charges and point dipoles
calculated using the soft core potentials.
λ = 1.0 returns the non-soft core version.
"""
function electrostatic_energy(pos::Array, charges::Array, dipoles::Array, λ::Any, pbc::Bool=false)
  @assert pbc == false "Periodic boundary condition not yet supported"
  qq = 0
  qμ = 0
  μμ = 0
  for (i, R) in enumerate(pos)
    for j = (i+1):length(charges)
      qq += soft_coulomb(R, charges[i], pos[j], charges[j], λ)
      qμ += soft_q_μ(R, charges[i], pos[j], dipoles[j], λ)
      μμ += soft_μ_μ(R, dipoles[i], pos[j], dipoles[j], λ)
    end
  end
  return qq + qμ + μμ
end

"""
Soft core Coulomb interaction between two point charges
If λ=1.0 the normal Coulomb is returned
Adapted from LAMMPS `pair_style coul/long/soft`
"""
function soft_coulomb(pos1, q1, pos2, q2, λ=0.9, α=10.0)
  r = norm(pos1 - pos2)
  return λ * q1 * q2 / (4*π*ϵ_0 * (α * (1 - λ)^2 + r^2)^0.5) * e * 1e10
end

"""
Soft core Coulomb interaction between point charge and point dipole
If λ=1.0 the normal Coulomb is returned
Analogously defined to LAMMPS `pair_style coul/long/soft`
"""
function soft_q_μ(pos1, q1, pos2, μ, λ=0.9, α=10.0)
  r12 = pos1 - pos2
  r = norm(r12)
  return λ * q1 * r12 ⋅ μ / (4*π*ϵ_0 * (α * (1 - λ)^2 + r^2)^1.5) *1e10^2 * (1e-21/c_light)
end

"""
Soft core Coulomb interaction between two point dipoles
If λ=1.0 the normal Coulomb is returned
Analogously defined to LAMMPS `pair_style coul/long/soft`
"""
function soft_μ_μ(pos1, μ1, pos2, μ2, λ=0.9, α=10.0)
  r12 = pos1 - pos2
  r = norm(r12)
  return λ / (4*π*ϵ_0 * (α * (1 - λ)^2 + r^2)^1.5) * (μ1 ⋅ μ2 - 3 * (μ1 ⋅ r12) * (r12 ⋅ μ2) / (α * (1 - λ)^2 + r^2)) *1e10^3 * (1e-21/c_light)^2 / e
end

end  # end of module Electro


text:
#=
function _unpack_local_solve!(m::GlobalOptimizer, opt::T) where T

    m._maximum_node_id = 0
    m._termination_status_code = MOI.get(opt, MOI.TerminationStatus())
    m._result_status_code = MOI.get(opt, MOI.PrimalStatus())
    m._feasible_solution_found = m._result_status_code == MOI.FEASIBLE_POINT

    if MOI.get(opt, MOI.ResultCount()) > 0

        objective_value = MOI.get(opt, MOI.ObjectiveValue())

        # corrects for standard printing multiplier
        if _optimization_sense(m._input_problem) == MOI.MAX_SENSE
            objective_value *= -1.0
        end

        m._global_lower_bound = objective_value
        m._global_upper_bound = objective_value
        m._objective_value = objective_value
        m._best_upper_value = objective_value
        m._solution_value = objective_value
    end

    map!(x -> MOI.get(opt, MOI.VariablePrimal(), x), m._solution, m._relaxed_variable_index)

    m._run_time = time() - m._start_time

    return nothing
end

function add_sv_or_aff_obj!(m::GlobalOptimizer, opt::T) where T
    if _objective_type(_input_problem(m)) == SINGLE_VARIABLE
        MOI.set(opt, MOI.ObjectiveFunction{SV}(), m._input_problem._objective_sv)
    elseif _objective_type(_input_problem(m)) == SCALAR_AFFINE
        MOI.set(opt, MOI.ObjectiveFunction{SAF}(), m._input_problem._objective)
    end
    return nothing
end
=#

add_nl_functions!(m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S} = add_nonlinear_functions!(m, m._input_problem._nlp_data)
add_nl_functions!(m::GlobalOptimizer{N,T,S}, x::Nothing) where {N,T<:Real,S} = nothing
add_nl_functions!(m::GlobalOptimizer{N,T,S}, x::MOI.NLPBlockData) where {N,T<:Real,S} = add_nl_functions!(m, x.evaluator)
function add_nonlinear_functions!(m::GlobalOptimizer{N,T,S}, x::JuMP.NLPEvaluator) where {N,T<:Real,S}

    nlp_data = m._input_problem._nlp_data
    MOI.initialize(x, Symbol[:Grad, :ExprGraph])

    # set nlp data structure
    m._working_problem._nlp_data = nlp_data

    # add subexpressions (assumes they are already ordered by JuMP)
    # creates a dictionary that lists the subexpression sparsity
    # by search each node for variables dict[2] = [2,3] indicates
    # that subexpression 2 depends on variables 2 and 3
    # this is referenced when subexpressions are called by other
    # subexpressions or functions to determine overall sparsity
    # the sparsity of a function is the collection of indices
    # in all participating subexpressions and the function itself
    # it is necessary to define this as such to enable reverse
    # McCormick constraint propagation
    relax_evaluator = m._working_problem._relaxed_evaluator
    has_subexpressions = length(x.m.nlp_data.nlexpr) > 0
    dict_sparsity = Dict{Int64,Vector{Int64}}()
    if has_subexpressions
        for i = 1:length(x.subexpressions)
            subexpr = x.subexpressions[i]
            push!(relax_evaluator.subexpressions, NonlinearExpression!(subexpr, dict_sparsity, i,
                                                                      x.subexpression_linearity,
                                                                      m.relax_tag))
        end
    end

    # scrubs udf functions using Cassette to remove odd data structures...
    # alternatively convert udfs to JuMP scripts...
    m.presolve_scrubber_flag && Script.scrub!(m._working_problem._nlp_data)
    if m.presolve_to_JuMP_flag
        Script.udf_loader!(m)
    end

    parameter_values = copy(x.parameter_values)

    # add nonlinear objective
    if x.has_nlobj
        # TODO: Should add a constraint to nlp_data... or evaluator not objective
        m._working_problem._objective_nl = BufferedNonlinearFunction(x.objective, MOI.NLPBoundsPair(-Inf, Inf),
                                                                     dict_sparsity, x.subexpression_linearity,
                                                                     m.relax_tag)
    end

    # add nonlinear constraints
    constraint_bounds = m._working_problem._nlp_data.constraint_bounds
    for i = 1:length(x.constraints)
        constraint = x.constraints[i]
        bnds = constraint_bounds[i]
        push!(m._working_problem._nonlinear_constr, BufferedNonlinearFunction(constraint, bnds, dict_sparsity,
                                                                              x.subexpression_linearity,
                                                                              m.relax_tag))
    end

    m._input_problem._nonlinear_count = length(m._working_problem._nonlinear_constr)
    m._working_problem._nonlinear_count = length(m._working_problem._nonlinear_constr)

    return nothing
end

function add_nl_evaluator!(m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S}
    add_nl_evaluator!(m, m._input_problem._nlp_data)
end
add_nl_evaluator!(m::GlobalOptimizer{N,T,S}, x::Nothing) where {N,T<:Real,S} = nothing
function add_nonlinear_evaluator!(m::GlobalOptimizer{N,T,S}, x::JuMP.NLPEvaluator) where {N,T<:Real,S}
    m._working_problem._relaxed_evaluator = Evaluator{N,T}()

    relax_evaluator = m._working_problem._relaxed_evaluator
    relax_evaluator.variable_count = length(m._working_problem._variable_info)
    relax_evaluator.user_operators = x.m.nlp_data.user_operators

    relax_evaluator.lower_variable_bounds = zeros(relax_evaluator.variable_count)
    relax_evaluator.upper_variable_bounds = zeros(relax_evaluator.variable_count)
    relax_evaluator.x                     = zeros(relax_evaluator.variable_count)
    relax_evaluator.num_mv_buffer         = zeros(relax_evaluator.variable_count)
    relax_evaluator.treat_x_as_number     = fill(false, relax_evaluator.variable_count)
    relax_evaluator.ctx                   = GuardCtx(metadata = GuardTracker(m.domain_violation_ϵ,
                                                                             m.domain_violation_guard_on))
    relax_evaluator.subgrad_tol           = m.subgrad_tol

    m._nonlinear_evaluator_created = true

    return nothing
end

function add_subexpr_buffers!(m::GlobalOptimizer{N,T,S}) where {N,T<:Real,S}
    relax_evaluator = m._working_problem._relaxed_evaluator
    relax_evaluator.subexpressions_eval = fill(false, length(relax_evaluator.subexpressions))

    return nothing
end

function add_variables(m::GlobalOptimizer, optimizer::T, variable_number::Int) where T

    variable_index = fill(VI(1), variable_number)
    for i = 1:variable_number
        @inbounds variable_index[i] = MOI.add_variable(optimizer)
        single_variable = SV(@inbounds variable_index[i])
        v_info = @inbounds m._working_problem._variable_info[i]
        if _is_fixed(v_info)
            MOI.add_constraint(optimizer, single_variable, _ET(vi))
        elseif _is_real_interval(vi)
            MOI.add_constraint(optimizer, single_variable, _IT(vi))
        elseif _is_greater_than(vi)
            MOI.add_constraint(optimizer, single_variable, _GT(vi))
        elseif _is_less_than(vi)
            MOI.add_constraint(optimizer, single_variable, _LT(vi))
        elseif _is_zero_one(vi)
            MOI.add_constraint(optimizer, single_variable, _ZO(vi))
        elseif _is_int_interval(vi)
            MOI.add_constraint(optimizer, single_variable, _Int(vi))
        end
    end

    return variable_index
end

"""
$(TYPEDSIGNATURES)

Loads variables, linear constraints, and empty storage for first nlp and
quadratic cut.
"""
function load_relaxed_problem!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    opt = m.relaxed_optimizer

    # add variables and indices and constraints
    wp = m._working_problem
    branch_num = 0

    variable_count = wp._variable_num
    for i = 1:variable_count

        vinfo = wp._variable_info[i]
        is_branch = m._branch_variables[i]
        is_branch && (branch_num += 1)

        var_indx = MOI.add_variable(opt)
        var = SV(var_indx)
        push!(m._relaxed_variable_index, var_indx)

        if vinfo.is_integer
        elseif vinfo.is_fixed
            ci_et = MOI.add_constraint(opt, var, ET(vinfo.lower_bound))
            is_branch && push!(m._relaxed_variable_eq, (ci_sv_et, branch_num))
        else
            if vinfo.has_lower_bound
                ci_gt = MOI.add_constraint(opt, var, GT(vinfo.lower_bound))
                is_branch && push!(m._relaxed_variable_gt, (ci_gt, branch_num))
            end
            if vinfo.has_upper_bound
                ci_lt = MOI.add_constraint(opt, var, LT(vinfo.upper_bound))
                is_branch && push!(m._relaxed_variable_lt, (ci_lt, branch_num))
            end
        end
    end

    # set node index to single variable constraint index maps
    m._branch_variable_num = branch_num
    m._node_to_sv_leq_ci = fill(CI{SV,LT}(-1), branch_num)
    m._node_to_sv_geq_ci = fill(CI{SV,GT}(-1), branch_num)
    for v in m._relaxed_variable_lt
        m._node_to_sv_leq_ci[v[2]] = v[1]
    end
    for v in m._relaxed_variable_gt
        m._node_to_sv_geq_ci[v[2]] = v[1]
    end

    # add linear constraints (ASSUME NOT BRIDGED... TODO: FIX LATER)
    _add_linear_constraints!(opt, m._input_problem)
    MOI.set(opt, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    return
end

function _mod_vinfo!(wp, ip, ci::CI{SV,S}) where S
    sv = MOI.get(ip, MOI.ConstraintFunction(), ci)
    sv_set = MOI.get(ip, MOI.ConstraintSet(), ci)
    wp._variable_info[sv.value] = VariableInfo(wp._variable_info[sv.value], sv_set)
    return nothing
end

function _add_decision_variables!(wp, ip::InputModel{T}) where {T<:Real}
    imodel = ip._input_model
    ip_variable_num = MOI.get(imodel, MOI.NumberOfVariables())
    wp._variable_num = ip_variable_num
    wp._variable_info = fill(VariableInfo{T}(), ip_variable_num)
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, LT{T}}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, GT{T}}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, ET{T}}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, ZO}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, MOI.Integer}()))
    any(isempty, wp._variable_info) && error("Variable bounds lead to infeasible problem.") # TODO: Message + Model Issue Return Code
    return nothing
end

function _add_constraint!(wp, ip::InputModel{T}) where {T<:Real}
    imodel = ip._input_model
    ip_variable_num = MOI.get(imodel, MOI.NumberOfVariables())
    wp._variable_num = ip_variable_num
    wp._variable_info = fill(VariableInfo{T}(), ip_variable_num)
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, LT{T}}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, GT{T}}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, ET{T}}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, ZO}()))
    foreach(x -> _mod_vinfo!(wp, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SV, MOI.Integer}()))
    any(isempty, wp._variable_info) && error("Variable bounds lead to infeasible problem.") # TODO: Message + Model Issue Return Code
    return nothing
end

for Q in (LT,GT)
    @eval function _add_constraint!(wp::ParsedProblem, ip::InputModel{T}, ci::CI{SAF{T},$Q{T}}) where {T<:Real}
        imodel = ip._input_model
        f = MOI.get(imodel, MOI.ConstraintFunction(), ci)
        s = MOI.get(imodel, MOI.ConstraintSet(), ci)
        push!(wp._saf_leq, AffineFunctionIneq(f, s))
        return nothing
    end
    @eval function _add_constraint!(wp::ParsedProblem, ip::InputModel{T}, ci::CI{SQF{T},$Q{T}}) where {T<:Real}
        imodel = ip._input_model
        f = MOI.get(imodel, MOI.ConstraintFunction(), ci)
        s = MOI.get(imodel, MOI.ConstraintSet(), ci)
        push!(wp._sqf_leq, BufferedQuadraticIneq(f, s))
        return nothing
    end
end
function _add_constraint!(wp::ParsedProblem, ip::InputModel{T}, ci::CI{SAF{T},ET{T}}) where {T<:Real}
    imodel = ip._input_model
    f = MOI.get(imodel, MOI.ConstraintFunction(), ci)
    s = MOI.get(imodel, MOI.ConstraintSet(), ci)
    push!(wp._saf_eq, AffineFunctionEq(f, s))
    return nothing
end
function _add_constraint!(wp::ParsedProblem, ip::InputModel{T}, ci::CI{SQF{T},ET{T}}) where {T<:Real}
    imodel = ip._input_model
    f = MOI.get(imodel, MOI.ConstraintFunction(), ci)
    s = MOI.get(imodel, MOI.ConstraintSet(), ci)
    push!(wp._sqf_eq, BufferedQuadraticEq(f, s))
    return nothing
end

function _add_constraint!(wp::ParsedProblem, ip::InputModel{T}, ci::CI{VECVAR,SOC_CONE}) where {T<:Real}
    imodel = ip._input_model
    f = MOI.get(imodel, MOI.ConstraintFunction(), ci)
    s = MOI.get(imodel, MOI.ConstraintSet(), ci)
    first_variable_loc = f.variables[1].value
    prior_lbnd = wp._variable_info[first_variable_loc].lower_bound
    wp._variable_info[first_variable_loc].lower_bound = max(prior_lbnd, zero(T))
    push!(wp._conic_second_order, BufferedSOC(f, s))
    return nothing
end

function _add_linear_constraints!(m::S, ip::InputModel{T}) where {S,T<:Real}
    imodel = ip._input_model
    foreach(x -> _add_constraint!(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SAF{T}, LT{T}}()))
    foreach(x -> _add_constraint!(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SAF{T}, GT{T}}()))
    foreach(x -> _add_constraint!(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SAF{T}, ET{T}}()))
    return nothing
end

function _add_quadratic_constraints!(m::S, ip::InputModel{T}) where {S,T<:Real}
    imodel = ip._input_model
    foreach(x -> _add_constraint!(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SQF{T}, LT{T}}()))
    foreach(x -> _add_constraint!(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SQF{T}, GT{T}}()))
    foreach(x -> _add_constraint!(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{SQF{T}, ET{T}}()))
    return nothing
end

function _add_conic_constraints!(m::S, ip::InputModel{T}) where {S,T<:Real}
    imodel = ip._input_model
    foreach(x -> MOI.add_constraint(m, imodel, x), MOI.get(imodel, MOI.ListOfConstraintIndices{VECVAR,SOC_CONE}()))
    return nothing
end

#=
Adds objective, converting problem to MinSensae
=#
_add_objective(d::Nothing, m::S, ip::InputModel{T}) where {S,T<:Real} = false
function _add_objective(d, m::S, ip::InputModel{T}) where {S,T<:Real}
    if d.has_objective
        # TODO: other manipulations of nonlinear objective...
        imodel = ip._input_model
        vi = VI(MOI.get(imodel, MOI.NumberOfVariables()) + 1)
        push!(m._variable_info, VariableInfo{T}())
        m._objective = SAF{T}(SAT{T}[SAT{T}(one(T), vi)], zero(T))
        return true
    end
    return false
end
function _add_objective!(m::S, ip::InputModel{T}) where {S,T<:Real}
    nl_obj_set = _add_objective(ip._nlp_data, m, ip)
    if !nl_obj_set
        saf = MOI.get(ip._input_model, MOI.ObjectiveFunction{SAF{T}}())
        m._objective = copy(saf)
        #if MOI.get(ip._input_model, MOI.OptimizationSense(1)) == MOI.MAX_SENSE # TODO: Fix me!
        #    m._objective = MOIU.operate(-, T, m._objective)
        #end
    end
    m._objective_parsed = AffineFunction(copy(m._objective), m._objective.constant)
    m._optimization_sense = MOI.MIN_SENSE
    return nothing
end

"""
$(TYPEDSIGNATURES)

Creates an initial node with initial box constraints and adds it to the stack.
"""
function create_initial_node!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}
    push!(m._stack, NodeBB{N,T}(ntuple(i -> _lower_bound(BranchVar, m, i), N),
                                ntuple(i -> _upper_bound(BranchVar, m, i), N),
                                -Inf, Inf, 1, 1, 0, BD_NONE))
    m._node_count = 1
    m._maximum_node_id += 1

    return nothing
end
function presolve_global!(t::ExtensionType, m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    wp = _working_problem(m)
    ip = _input_problem(m)

    _add_decision_variables!(wp, ip)      # copy constraints with storage
    _add_linear_constraints!(wp, ip)
    _add_quadratic_constraints!(wp, ip)
    _add_conic_constraints!(wp, ip)

    _add_objective!(wp, ip)               # set objective function
    _label_branch_variables!(m)           # labels default variables to branch on

    add_nl_evaluator!(m)
    add_nl_functions!(m)
    add_subexpression_buffers!(m)

    load_relaxed_problem!(m)
    create_initial_node!(m)

    branch_variable_num = m._branch_variable_num

    m._current_xref             = fill(zero(T), branch_variable_num)
    m._candidate_xref           = fill(zero(T), branch_variable_num)
    m._current_objective_xref   = fill(zero(T), branch_variable_num)
    m._prior_objective_xref     = fill(zero(T), branch_variable_num)
    m._lower_lvd                = fill(zero(T), branch_variable_num)
    m._lower_uvd                = fill(zero(T), branch_variable_num)

    # populate in full space until local MOI nlp solves support constraint deletion
    # uses input model for local nlp solves... may adjust this if a convincing reason
    # to use a reformulated upper problem presents itself
    m._lower_result_count_max  = 100
    m._lower_result_count      = 0
    m._lower_solution          = Vector{T}[]
    for j = 1:_lower_result_count_max
        push!(m._lower_solution, zeros(m._working_problem._variable_num))
    end

    m._cut_solution        = zeros(T, m._working_problem._variable_num)
    m._solution            = zeros(T, m._working_problem._variable_num)
    m._upper_solution      = zeros(T, m._working_problem._variable_num)
    m._upper_variables     = fill(VI(-1), m._working_problem._variable_num)

    # add storage for fbbt
    m._lower_fbbt_buffer   = zeros(T, m._working_problem._variable_num)
    m._upper_fbbt_buffer   = zeros(T, m._working_problem._variable_num)

    # add storage for obbt ( perform obbt on all relaxed variables, potentially)
    m._obbt_working_lower_index = fill(false, branch_variable_num)
    m._obbt_working_upper_index = fill(false, branch_variable_num)
    m._old_low_index            = fill(false, branch_variable_num)
    m._old_upp_index            = fill(false, branch_variable_num)
    m._new_low_index            = fill(false, branch_variable_num)
    m._new_upp_index            = fill(false, branch_variable_num)
    m._lower_indx_diff          = fill(false, branch_variable_num)
    m._upper_indx_diff          = fill(false, branch_variable_num)
    m._obbt_variable_num        = branch_variable_num

    # set subgradient refinement flag
    wp._relaxed_evaluator.is_post = m.subgrad_tighten
    wp._relaxed_evaluator.subgrad_tighten = m.subgrad_tighten
    wp._relaxed_evaluator.reverse_subgrad_tighten =  m.reverse_subgrad_tighten

    m._presolve_time = time() - m._parse_time

    return nothing
end


text:
module SimpleQuantum
using LinearAlgebra, SparseArrays, SimpleRandom

import Base: show, ==, length, kron

export measure!, (⊗), lg

RR = Float64            # shortcut for real
CC = Complex{RR}        # shortcut for complex

"""
`A⊗B` is the Kronecker (tensor) product of `A` and `B`.
"""
⊗(A) = A
⊗(A, B) = kron(A, B)
⊗(A, B, C...) = ⊗(A, ⊗(B, C...))

"""
`lg(n)` is a fast, base-2 log of `n` assuming `n` is a power of 2.
"""
function lg(n::Int)::Int
    count = 0
    while n>1
        n = n>>1
        count += 1
    end
    return count
end

include("quantum_state.jl")
include("qubit.jl")
include("register.jl")
include("gate.jl")

end  # end of module


text:
using MINE
using Base.Test

function test()
    para = Parameter(0.6, 15)
    x = linspace(0,1,1000)
    y = sin(10π*x) +x
    prob = Problem(x,y)
    ppara = pointer(para)
    pprob = pointer(prob)
    score = compute_score(pprob, ppara)
    mic(pointer(score))
    free_score(pointer(score))
end

test()


text:
model = @reaction_network SSTest begin
    (d_x, d_y), (x, y) --> 0
    mm(y, 2, 1), 0 --> x
    1., 0 --> y
end d_x d_y

u0 = [0.1, 0.1]
p = [1., 1.]

prob = SteadyStateProblem{true}(model, u0, p)
sol = solve(prob, SSRootfind())
@test sol.u ≈ [1, 1]

u0 = [0.1, 0.1]
prob = SteadyStateProblem(model, u0, p)
sol = solve(prob, SSRootfind())
@test sol.u ≈ [1, 1]


text:
module ABCDBeamTrace

export FreeSpace, Interface, ThinLens, Mirror, Tan, Sag, RTM, Beam
export beamtrace, spotsize, location, discretize
export waistradiusfunc
export WithBeam

using RecipesBase
import Colors, Interpolations

include("elements.jl")
include("beamtrace.jl")
include("comparisons.jl")
include("plot-recipes.jl")

end # module


text:
"""
    mutable struct Diagnostic{T, N} <: AbstractDiagnostic

A diagnostic that includes `N` elements of type `T`.

$(TYPEDFIELDS)
"""
mutable struct Diagnostic{T, N} <: AbstractDiagnostic
    "function that returns the diagnostic via `calc(prob)`"
   calc :: Function
    "the relevant problem for this diagnostic"
   prob :: Problem
    "vector where the diagnostic time-series is saved"
   data :: Vector{T}
    "vector with the times for which the diagnostic was saved"
      t :: Vector{Float64}
    "vector with the problem's `step` for which the diagnostic was saved"
  steps :: Vector{Int}
    "integer denoting how often (every how many `problem.step`s) to save the diagnostic"
   freq :: Int
    "integer denoting how many times the `diagnostic.data` was updated"
      i :: Int
end

"""
    Diagnostic(calc, prob; freq=1, nsteps=100, ndata=ceil(Int, (nsteps+1)/freq))

Construct a diagnostic that stores the result of `calc(prob)` with frequency `freq`.

Keywords
========
  - `freq`: Diagnostic is saved every `freq` steps.
  - `nsteps`: The total number of steps in problem.
  - `ndata`: The number of diagnostics to be saved.
"""
function Diagnostic(calc, prob; freq=1, nsteps=100, ndata=ceil(Int, (nsteps+1)/freq))
  firstvalue = calc(prob)
  T = typeof(firstvalue)

  data = Vector{T}(undef, ndata)
  t = Vector{Float64}(undef, ndata)
  steps = Vector{Int}(undef, ndata)

  data[1] = firstvalue
  t[1] = prob.clock.t
  steps[1] = prob.clock.step

  return Diagnostic{T,ndata}(calc, prob, data, t, steps, freq, 1)
end

"""
    extend!(diag, n)

Extend the `data`, `time`, and `steps` vectors of the diagnostic `diag` by `n`.
"""
function extend!(diag, n)
  resize!(diag.data, length(diag.steps)+n)
  resize!(diag.t, length(diag.steps)+n)
  resize!(diag.steps, length(diag.steps)+n)

  return nothing
end

"""
    extend!(diag::Diagnostic{T,N})

Double the extend of the `data`, `time`, and `steps` vectors of the diagnostic `diag`.
"""
extend!(diag::Diagnostic{T,N}) where {T,N} = extend!(diag, N)

"""
    update!(diag)

Update `diag`.
"""
function update!(diag, i)
  i <= length(diag.steps) || extend!(diag)

   diag.data[i] = diag() #diag.calc(prob)
      diag.t[i] = diag.prob.clock.t
  diag.steps[i] = diag.prob.clock.step
         diag.i = i

  return nothing
end

"""
    increment!(diag)

Increment the Diagnostic `diag`, or an array of Diagnostics `diags`.
"""
function increment!(diag)
  diag.prob.clock.step % diag.freq == 0 && update!(diag, diag.i+1)

  return nothing
end

function increment!(diags::AbstractVector)
  for diag in diags
    increment!(diag)
  end

  return nothing
end

"""
e.g. `plot(energydiag[:t], energydiag[:data])`.
"""
getindex(d::Diagnostic, s::Union{AbstractString,Symbol}) = getfield(d, Symbol(s))[1:d.i]
getindex(d::Diagnostic, idx...) = getindex(d.data, idx...)

# d() = current value
# e = Diagnostic(energy, ...)
# e() returns energy at the current time.
# sensible?
(d::Diagnostic)() = d.calc(d.prob)

show(io::IO, d::Diagnostic{T, N}) where {T, N} =
     print(io, "Diagnostic\n",
               "  ├─── calc: ", d.calc, '\n',
               "  ├─── prob: ", summary(d.prob), '\n',
               "  ├─── data: ", summary(d.data), '\n',
               "  ├────── t: ", summary(d.t), '\n',
               "  ├── steps: ", summary(d.steps), '\n',
               "  ├─── freq: ", d.freq, '\n',
               "  └────── i: ", d.i)


text:
using Test

# for convenience, when debugging, the build step can be disabled when running the tests with
#   julia test/workdir/tests.jl skipbuild
if !("skipbuild" in ARGS)
    rm(joinpath(@__DIR__, "builds"), recursive=true, force=true) # cleanup of previous test run
    include(joinpath(@__DIR__, "make.jl"))
end

@testset "makedocs workdir" begin
    # test for the default build
    @test isfile(joinpath(@__DIR__, "builds", "default", "root_index.txt"))
    @test isfile(joinpath(@__DIR__, "builds", "default", "root_file.txt"))
    @test isfile(joinpath(@__DIR__, "builds", "default", "subdir", "subdir_index.txt"))
    @test isfile(joinpath(@__DIR__, "builds", "default", "subdir", "subdir_file.txt"))

    # absolute path
    @test !isfile(joinpath(@__DIR__, "builds", "absolute", "root_index.txt"))
    @test !isfile(joinpath(@__DIR__, "builds", "absolute", "root_file.txt"))
    @test !isfile(joinpath(@__DIR__, "builds", "absolute", "subdir", "subdir_index.txt"))
    @test !isfile(joinpath(@__DIR__, "builds", "absolute", "subdir", "subdir_file.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "absolute-workdir", "root_index.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "absolute-workdir", "root_file.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "absolute-workdir", "subdir_index.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "absolute-workdir", "subdir_file.txt"))

    # relative path
    @test !isfile(joinpath(@__DIR__, "builds", "relative", "root_index.txt"))
    @test !isfile(joinpath(@__DIR__, "builds", "relative", "root_file.txt"))
    @test !isfile(joinpath(@__DIR__, "builds", "relative", "subdir", "subdir_index.txt"))
    @test !isfile(joinpath(@__DIR__, "builds", "relative", "subdir", "subdir_file.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "relative-workdir", "root_index.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "relative-workdir", "root_file.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "relative-workdir", "subdir_index.txt"))
    @test  isfile(joinpath(@__DIR__, "builds", "relative-workdir", "subdir_file.txt"))
end