text:
# Kernel Op
# Input: x[nd,nx]
# Support vector matrix: s[nd,ns]
# Kernel function: K[ns,nx] = kernel(s,x)  -- a new representation for x
# Output: y[nc,nx] = w[nc,ns] * K[ns,nx]  -- like the mmul layer
# Initially: s[nd,0], w[nc,0]
#
# A subtype of kernel needs to define:
# Fields: s, w, v, x, y, dy
# Functions: kernel
# Initialize: w[nc,0]

abstract Kernel <: Op

# The forw function is just mmul applied to a transformed input:

function forw(l::Kernel, x; train=true, o...)
    initforw(l, x)
    # kernel fn is defined by subtypes of kernel and returns the kernel matrix K
    l.y = (train ? l.v : l.w) * kernel(l)
end

function back(l::Kernel, dy; returndx=false, o...)
    @assert issimilar(dy, l.y)
    l.dy = dy  # we assume this is from PercLoss: -1 for correct answer, +1 for predicted answer
    returndx && error("Kernel layers do not know how to return dx")
end

function update!(l::Kernel; o...)
    (yrows, ycols) = size(l.dy)
    (xrows, xcols) = size(l.x)
    @assert ycols == xcols
    snew = similar(l.x, (xrows, 0))
    wnew = similar(l.w, (yrows, 0))
    for j=1:ycols
        for i=1:yrows
            if l.dy[i,j] != 0
                snew=[snew l.x[:,j]] ### BIG ALLOC
                wnew=[wnew (-l.dy[:,j])]
                break
            end
        end
    end
    l.w += l.v
    l.w = [l.w wnew]
    l.v = [l.v wnew]
    l.s = [l.s snew]
end

function initforw(l::Kernel, x)
    (xrows, xcols) = size2(x)
    l.x = (size(x)==(xrows,xcols) ? x : reshape(x, xrows, xcols))
    isdefined(l,:s) || (l.s=similar(x, (xrows, 0)))
    (srows, scols) = size(l.s)
    @assert srows == xrows
    isdefined(l,:w) || (l.w=Array(eltype(x), l.n, 0))
    (wrows, wcols) = size(l.w)
    @assert scols == wcols
    # l.v is the perceptron weights, l.w is the averaged (summed) weights
    isdefined(l,:v) || (@assert isempty(l.w); l.v = l.w)
end

# TODO: explain back and update
# TODO: what if we froze the SV and used regular updates (copy sv into mmul)?

# using Base.LinAlg.BLAS: gemm!

# OK forget about cuda, just do Array and SparseArray for now.
# Use generic code, forget about efficiency.

# DONE: forw will not work for the first batch of x when sv is empty!
# Can we multiply zero dim matrices? yes.

# No need to alloc l.y, no in-place ops for sparse
# similar!(l, :y, l.sv, (scols, xcols))

# Careful: poly may turn a sparse array into dense if b!=0.
# If b==0 we can just go over nonzero entries with sparse.

# if GPU
# poly(y::AbstractCudaArray{Float32}, a, b, d)=(ccall((poly32,libkunet),Void,(Ptr{Float32},Float32,Float32,Float32),y,a,b,d);y)
# poly(y::AbstractCudaArray{Float64}, a, b, d)=(ccall((poly64,libkunet),Void,(Ptr{Float64},Float64,Float64,Float64),y,a,b,d);y)
# # TODO: test sparse on cpu/gpu
# # TODO: implement poly32/64 in cuda
# end

    # back: Kmul passes back its own dy, nonzero columns indicate mistakes
    # We do nothing here except record it
    # Actual creation of support vectors done in update!()
    # TODO: we could compute dx like we do in mmul?

# TODO: back
# q: what is dy? - use the same dy as Kmul
# a: we need indices of x vectors with wrong answers to turn into sv.
# a: that means dy not the same size as y?
# a: argue for integrating beta here? pro:more like other layers, con:repeat code for other kernels
# q: what does back do? nothing? (store incoming info?)
# q: what does update do? (add new sv?)
# q: what do we pass back? (nothing?)
# a: does primal perceptron compute dx?

# TODO: worry about efficiency of concat later (consider making sv a list?)
# See this for efficiency: https://groups.google.com/forum/#!topic/julia-users/B4OUYPFM5L8
# TODO: check if concat works for cpu/gpu, sparse/full
# TODO: gpu version

# TODO: fix the starting empty sv problem
# TODO: implement Kmul
# TODO: implement Rbfk
# TODO: see about code repetition, Only forward will be different, back and update the same.
# we can define abstract Kernel to avoid repetition.
# if we define Kernel, we can also have the beta mult in there!

# TODO: uniq sv, cotter13


text:
using Base.Test
using Lint
using Mocking

include("../src/Playground.jl")
using Playground

TEST_DIR = pwd()

#TEST_TMP_DIR = mktempdir("../test/")
TEST_TMP_DIR = "$TEST_DIR/tmp/"
TEST_PLAYGROUND_DIR = joinpath(TEST_TMP_DIR, "playground")
TEST_CONFIG = load_config(DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)

mkpath(TEST_TMP_DIR)
mkpath(TEST_PLAYGROUND_DIR)

Playground.init(TEST_CONFIG)

# Order matters.
tests = [
    "lint",
    "utils",
    "list",
    "parsing",
    "install",
    "create",
    "activate",
    "execute",
    "clean"
]

println("Running tests:")

for t in tests
    tfile = string(t, ".jl")
    println(" * $(tfile) ...")
    include(joinpath(TEST_DIR, tfile))
end

main(["install", "download", "0.4", "--labels", "julia-0.3"], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)
main(["install", "link", joinpath(TEST_CONFIG.dir.bin, "julia-0.3"), "--labels", "julia-stable-dir"], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)
main(["create"], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)
main(["exec", "ls -al", joinpath(TEST_DIR, ".playground")], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)
main(["list"], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)
main(["clean"], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)
main(["rm", "--dir", joinpath(TEST_DIR, ".playground")], DEFAULT_CONFIG, TEST_PLAYGROUND_DIR)

ext: jl


text:
module RebuggerTesting

const cbdata1 = Ref{Any}(nothing)
const cbdata2 = Ref{Any}(nothing)

# Do not alter the line number at which `foo` occurs
foo(x, y) = nothing

snoop0()             = snoop1("Spy")
snoop1(word)         = snoop2(word, "on")
snoop2(word1, word2) = snoop3(word1, word2, "arguments")
snoop3(word1, word2, word3::T; adv="simply", morekws...) where T = error("oops")

kwvarargs(x; kw1=1, kwargs...)  = kwvarargs2(x; kw1=kw1, kwargs...)
kwvarargs2(x; kw1=0, passthrough=true) = (x, kw1, passthrough)

destruct(x, (a, b), y) = a

struct HasValue
    x::Float64
end
const hv_test = HasValue(11.1)

(hv::HasValue)(str::String) = hv.x

@inline kwfuncerr(y) = error("stop")
@noinline kwfuncmiddle(x::T, y::Integer=1; kw1="hello", kwargs...) where T = kwfuncerr(y)
@inline kwfunctop(x; kwargs...) = kwfuncmiddle(x, 2; kwargs...)

end

module RBT2

using ..RebuggerTesting

bar(::Int) = 5
RebuggerTesting.foo() = bar(1)

end


text:
function print_and_run_cmd(cmd::Base.AbstractCmd)

    println(cmd)

    run(cmd)

end

export print_and_run_cmd


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, true, false),),
    ((2, false, true),),
    ((2, false, true), :SOCExpPSD, relaxed_tols),
    ]
insts["fast"] = [
    ((10, true, false),),
    ((10, false, true), nothing, relaxed_tols),
    ((10, false, true), :SOCExpPSD),
    ((100, true, false),),
    ((100, false, true), nothing, relaxed_tols),
    ((100, false, true), :SOCExpPSD, relaxed_tols),
    ((1000, true, false),),
    ]
insts["various"] = [
    ((500, true, false),),
    ((500, false, true),),
    ((500, true, false), :SOCExpPSD),
    ((500, false, true), :SOCExpPSD),
    ((1000, true, false),),
    ((1000, false, true),),
    ((1000, true, false), :SOCExpPSD),
    ((1000, false, true), :SOCExpPSD, relaxed_tols),
    ((2000, true, false),),
    ((2000, false, true),),
    ((2000, true, false), :SOCExpPSD),
    ]
return (MaxVolumeJuMP, insts)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.331138821460632 <= q <= 11.287379993959185)

                     add_NL_constraint(m, :(swish(0.16893740747119512 + 0.361142465903193*swish(0.3759145594542668 + 0.23452173305395663*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.21914525995597334*swish(-0.017859443401365915 + -0.7776638835685663*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.7287835676612944*swish(0.6615910141517634 + 0.18731902906517428*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + swish(-0.9123133183860701 + 0.273513749523993*swish(0.3759145594542668 + 0.23452173305395663*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.052599029592083735*swish(-0.017859443401365915 + -0.7776638835685663*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.2642412570560455*swish(0.6615910141517634 + 0.18731902906517428*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + swish(0.5514504616591211 + 0.27732859795053555*swish(0.3759145594542668 + 0.23452173305395663*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + -0.08702067567597149*swish(-0.017859443401365915 + -0.7776638835685663*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.40713382739201665*swish(0.6615910141517634 + 0.18731902906517428*swish(-0.8645963938709578 + -0.1467408759040536*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*swish(-0.9609056399980664 + 0.06060702618594416*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*swish(-0.9973666264873402 + 0.1126936245130743*swish(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*swish(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*swish(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@with_kw struct NKTable
    wavelength::Vector{Float64}
    n::Vector{Float64}
    k::Vector{Float64}
    wavelength_unit::Float64 = 1 * 10^-6
end

epsilon_real(model::NKTable) = model.n .^ 2 .- model.k .^ 2
epsilon_imag(model::NKTable) = 2 * model.n .* model.k

epsilon(model::NKTable) = epsilon_real(model) + im * epsilon_imag(model)

wavelength_range(model::NKTable) = (first(model.wavelength),
    last(model.wavelength))

@recipe function plot(model::NKTable;
    xlabel = "Wavelength μm", color = [:blue :red], label = ["n" "k"], linewidth = [3, 3])
    xlabel --> xlabel
    color --> color
    model.wavelength, [model.n, model.k]
end

n_from_epsilon(eps) = real(sqrt(eps))
k_from_epsilon(eps) = imag(sqrt(eps))
epsilon_from_nk(n, k) = n^2 - k^2 + 2 * im * n * k

function load_nk(filepath::String)
    data = readdlm(filepath)
    return NKTable(wavelength = data[:,1],n = data[:,2], k = data[:,3])
end


text:
function taxicabdistanceXY(d::DisplacementPositions)
    return TaxiCabDistanceXY(
        abs(d.p.x - d.q.x),
        abs(d.p.y - d.q.y))
end

function taxicabdistance(d::DisplacementPositions)
    TaxiCabDistance(abs(d.p.x - d.q.x) + abs(d.p.y - d.q.y))
end

function periodicdomainshift(distance,domain)
     return (distance > domain/2 ? domain - distance : distance)
end

function taxicabdistance(d::DisplacementPositions,domain)
     TaxiCabDistance(
          periodicdomainshift(abs(d.p.x - d.q.x),domain) +
               periodicdomainshift(abs(d.p.y - d.q.y),domain)
     )
end

function taxicabdistanceXY(d::TaxiCabDistanceXY,domain)
   dx = abs(periodicdomainshift(d.dₓ,domain))
   dy = abs(periodicdomainshift(d.dy,domain))
   return dx + dy
end

walker_distance(w) = taxicabdistance(DisplacementPositions(w[1],w[end]),domain)

function dₜ(w,domain)
     d = DisplacementPositions(
          Walker2D(w[1][1],w[1][2]),
          Walker2D(w[end][1],w[end][2])
     )
     D = taxicabdistanceXY(
          taxicabdistanceXY(d),domain
          )
     return D
end


text:
##########
##########
# Functions to perform perturbations on the food webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

# set up =======================================================================
import Pkg
Pkg.instantiate()

include("./04_eqn_component_functions.jl")
include("./01_species_list_creation.jl")
#include("./05_network_metrics.jl")
include("./08_perturbation_sims_functions.jl")

#@info{"All prereqs done..."}

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# pulse perturbations ==========================================================

all_species_lists = readdlm("./data/intermediate-files/species_lists_for_simulations.csv", ',', Any)

# little line to set the network id (this happens from the shell line)
network_id = parse(Int64, ARGS[1])

species_list = all_species_lists[all_species_lists[:,network_id_pos] .== network_id,:]
network = create_interaction_matrix(species_list)

# set parameters - pulse
c = 2 # constant amount
x = 4000 # duration
f = 1 # frequency
z = c/(x*(1/f))

# set vector callback
extinction_callback =
    VectorContinuousCallback(extinction_threshold,
                            extinction_affect!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))
steady_callback = TerminateSteadyState(1e-12,
                                        1e-9,
                                        DiffEqCallbacks.allDerivPass)

# set callback set
callback_set_perturb = extinction_callback
callback_set = CallbackSet(extinction_callback, steady_callback)

#  allocate x0
u0 = Float64[]
for i in 1:size(species_list, 1)
    append!(u0, species_list[i, biomass_pos])
end

# set perturbation tspan
tspan = (1000.0, 5000.0)

# set constant stop values
const stop1 = [1000.0]
const stop2 = [5000.0]
const tstop = [1000.;5000.]

# loop through random selections of species to perturb
#@info{"starting loop..."}

for j in 1:10

    # pick set of random species to perturb
    all = Array(1:size(species_list, 1))
    perturb_sp = rand(Array(1:size(species_list, 1)),
                        Int(round(size(species_list, 1)*(1/3))))
    non_perturb = setdiff(Set(all), Set(perturb_sp))

    # set times to save from the simulation
    save_times = [collect(range(1000.0, 1250.0, step = 0.1));
                collect(range(1250, 4800, step = sample(200:400, 1)[1]));
                collect(range(4800.0, 5000.0, step = 0.1))]

    # set p
    p = [species_list, network, 0.3, z, perturb_sp, non_perturb]

    prob_press = ODEProblem(biomass_perturb_random,
                                        u0,
                                        tspan,
                                        p,
                                        tstops=tstop)
    @time sol_prob_press = solve(prob_press,
                                        Tsit5(),
                                        callback = callback_set,
                                        # save every part for the perturbation and
                                        # the first 150 time steps after the
                                        # perturbation, but then just grab a few
                                        # from there to the end
                                        saveat = save_times)
    # save output of the solver
    output = hcat(sol_prob_press.t, sol_prob_press.u)
    # write out the biomass outputs
    name = string("press-", network_id, "-", j, "-", "sims.csv")
    path = "./data/simulation-outputs/press-sims/"
    writedlm(string(path,name), output, ",")
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
"""
# Description

Convert a value `x`, into an expression that evaluates to that value. Useful
for taking a value constructed during macro evaluation and placing it in the
resulting output expression. Only handles a very limited set of types needed
to construct `Expression` objects.

# Arguments

1. `x::Any`: A value to be converted into an expression.

# Return Values

1. `e′::Any`: An expression that evaluates to `x`.

# Examples

```
julia> column_names_tuple_expr((ColumnName(:a, false), ColumnName(:b, true)))
:((:a, $(Expr(:escape, :b))))
```
"""
function column_names_tuple_expr(@nospecialize(x::Any))
    if isa(x, NTuple{N, ColumnName} where N)
        Expr(
            :tuple,
            map(column_names_tuple_expr, x)...,
        )
    elseif isa(x, ColumnName)
        if x.is_dynamic === false
            QuoteNode(x.name)
        else
            esc(x.name)
        end
    else
        error("Input type cannot be converted by column_names_tuple_expr")
    end
end


text:
abstract type Method end

abstract type BracketingMethod <: Method end
abstract type LinearMethod <: Method end
abstract type NonlinearMethod <: Method end

abstract type DirectMethod <: LinearMethod end
abstract type IterativeMethod <: LinearMethod end

abstract type NewtonMethod <: NonlinearMethod end
abstract type LinesearchMethod <: NonlinearMethod end

struct Newton <: NewtonMethod end
struct DFP <: NewtonMethod end
struct BFGS <: NewtonMethod end

struct Backtracking <: LinesearchMethod end
struct Bisection <: LinesearchMethod end
struct Quadratic <: LinesearchMethod end

struct Static{T} <: LinesearchMethod
    alpha::T
    Static(alpha::T = 1.0) where {T} = new{T}(alpha)
end

Base.show(io::IO, alg::Newton) = print(io, "Newton")
Base.show(io::IO, alg::DFP) = print(io, "DFP")
Base.show(io::IO, alg::BFGS) = print(io, "BFGS")

Base.show(io::IO, alg::Static) = print(io, "Static")
Base.show(io::IO, alg::Backtracking) = print(io, "Backtracking")
Base.show(io::IO, alg::Bisection) = print(io, "Bisection")
Base.show(io::IO, alg::Quadratic) = print(io, "Quadratic Polynomial")


text:
using Test
using BetweenFlags

@testset "Scaling" begin
  flag_set = FlagSet([
        FlagPair{GreedyType}(
        StartFlag("|", [" "], [" "]),
        StopFlag( "|", [" "], [" "])
        )
        ])
  base_text = "baz | foo | baz | bar | foobaz "
  timings = Dict()
  for n in 10 .^ (1:3)
    text = repeat(base_text,n)
    Δt = @elapsed TokenStream(text, flag_set)
    @show n, Δt
  end
end

# using Profile
# @testset "Profile" begin
#   flag_set = FlagSet([
#         FlagPair{GreedyType}(
#         StartFlag("|", [" "], [" "]),
#         StopFlag( "|", [" "], [" "])
#         )
#         ])
#   base_text = "baz | foo | baz | bar | foobaz "
#   base_text = "foo | θ | bar | ψ | foobar "

#   text = repeat(base_text,1)
#   ts = TokenStream(text, flag_set)

#   text = repeat(base_text,10^3)
#   ts = @profile TokenStream(text, flag_set)
#   Profile.print()
# end


text:
@system Mass begin
    pheno: phenology ~ hold
    NU: nodal_units ~ hold

    initial_leaf_ratio ~ hold
    potential_leaf_area_increase ~ hold

    shoot_carbon ~ hold
    root_carbon ~ hold
    leaf_carbon ~ hold
    # sheath_carbon ~ hold
    # scape_carbon ~ hold
    # bulb_carbon ~ hold
    stem_carbon ~ hold
    nonstructural_carbon_mass ~ hold

    agefn ~ hold

    # seed weight g/seed
    initial_seed_mass => 0.275 ~ preserve(u"g", parameter)

    # 10% available
    seed_mass_export_limit => 10 ~ preserve(u"percent/d")
    seed_mass_export_rate(seed_mass, seed_mass_export_limit, β=pheno.BF.ΔT) => begin
        # reserved in the propagule (e.g., starch in endosperm of seeds)
        #HACK ratio should depend on growth stage, but fix it for now
        T_effect = β
        seed_mass * T_effect * seed_mass_export_limit
    end ~ track(u"g/d")

    #HACK carbon mass of seed is pulled in the reserve
    seed_mass(seed_mass_export_rate) => begin
        -seed_mass_export_rate
    end ~ accumulate(init=initial_seed_mass, u"g")

    #stem(x=NU["*"].stem.mass) => begin # for maize
    # total_sheath_mass(x=NU["*"].sheath.mass) => begin # for garlic
    #     # dt the addition of C_reserve here only serves to maintain a total for the mass. It could have just as easily been added to total mass.
    #     # C_reserve is added to stem here to represent soluble TNC, SK
    #     #sum(typeof(0.0u"g")[nu.stem.mass' for nu in NU]) + self.p.carbon.reserve
    #     #sum(typeof(0.0u"g")[nu.sheath.mass' for nu in NU]) + self.p.carbon.reserve
    #     #FIXME carbon not ready yet
    #     sum(x)
    # end ~ track(u"g")

    initial_leaf_mass(initial_seed_mass, initial_leaf_ratio) => begin
        initial_seed_mass * initial_leaf_ratio
    end ~ track(u"g")

    # this is the total mass of active leaves that are not entirely dead (e.g., dropped).
    # It would be slightly greather than the green leaf mass because some senesced leaf area is included until they are complely aged (dead), SK
    active_leaf_mass(NU #=, x=NU["*"].leaf.mass =#) => begin
        sum(typeof(0.0u"g")[nu.leaf.mass' for nu in NU if !nu.leaf.dropped'])
    end ~ track(u"g")
    #TODO: support complex composition (i.e. `!`(leaf.dropped)) in condition syntax?
    #active_leaf_mass(x=NU["*/!leaf.dropped"].leaf.mass) => sum(x) ~ track(u"g")

    dropped_leaf_mass(NU #=, x=NU["*"].leaf.mass =#) => begin
        sum(typeof(0.0u"g")[nu.leaf.mass' for nu in NU if nu.leaf.dropped'])
    end ~ track(u"g")
    #TODO: support more referencing options (i.e. "leaf.dropped") in condition syntax?
    #dropped_leaf(x=NU["*/leaf.dropped"].leaf.mass) => sum(x) ~ track(u"g")

    total_leaf_mass(x=NU["*"].leaf.mass) => begin
        # this should equal to activeLeafMass + droppedLeafMass
        sum(x)
    end ~ track(u"g")

    # for maize

    # ear_mass(ear.mass) ~ track(u"g")

    # for garlic

    root_mass(root_carbon) ~ accumulate(u"g")
    leaf_mass(leaf_carbon) ~ accumulate(u"g")
    stem_mass(stem_carbon) ~ accumulate(u"g")
    # sheath_mass(sheath_carbon) ~ accumulate(u"g")
    # scape_mass(scape_carbon) ~ accumulate(u"g", when=!pheno.scape_removed)
    # visible_scape_mass(scape_carbon) ~ accumulate(u"g", when=pheno.scape_appeared & !pheno.scape_removed)

    # bulb_structural_mass(bulb_carbon) ~ accumulate(u"g")
    # bulb_mass(bulb_structural_mass, nonstructural_carbon_mass) => begin
    #     bulb_structural_mass + nonstructural_carbon_mass
    # end ~ track(u"g")

    # stalk_mass(sheath_mass, scape_mass) => begin
    #     #FIXME inconsistency: stem vs. sheath
    #     sheath_mass + scape_mass
    # end ~ track(u"g")

    #shoot_mass(shoot_carbon) => begin
    # shoot_mass(seed_mass, stalk_mass, leaf_mass, bulb_mass) => begin
    #     # for maize
    #     #seed_mass + stem_mass + leaf_mass + ear_mass
    #     # for garlic
    #     seed_mass + stalk_mass + leaf_mass + bulb_mass
    # end ~ track(u"g")

    #Jim's lettuce revision
    shoot_mass(seed_mass, stem_mass, leaf_mass) => begin
        # for maize
        #seed_mass + stem_mass + leaf_mass + ear_mass
        # for lettuce
        seed_mass + stem_mass + leaf_mass
    end ~ track(u"g")

    total_mass(shoot_mass, root_mass) => begin
        #HACK include mobilized carbon pool (for more accurate mass under germination)
        #shoot_mass + root_mass + carbon.pool
        shoot_mass + root_mass
    end ~ track(u"g")

    living_root_mass(agefn, root_mass) => begin
        agefn * root_mass
    end ~ track(u"g")

    living_leaf_mass(agefn, leaf_mass) => begin
        agefn * leaf_mass
    end ~ track(u"g")

    living_stem_mass(agefn, stem_mass) => begin
        agefn * stem_mass
    end ~ track(u"g")

    # living_sheath_mass(agefn, sheath_mass) => begin
    #     agefn * sheath_mass
    # end ~ track(u"g")

    # living_stalk_mass(living_sheath_mass, scape_mass) => begin
    #     living_sheath_mass + scape_mass
    # end ~ track(u"g")

    # living_shoot_mass(seed_mass, living_stalk_mass, living_leaf_mass) => begin
    #     seed_mass + living_stalk_mass + living_leaf_mass
    # end ~ track(u"g")

    #Jim's lettuce revision
    living_shoot_mass(seed_mass, living_stem_mass,  living_leaf_mass) => begin
        seed_mass + living_stem_mass + living_leaf_mass
    end ~ track(u"g")

    living_total_mass(living_shoot_mass, living_root_mass, bulb_mass) => begin
        living_shoot_mass + living_root_mass + bulb_mass
    end ~ track(u"g")

    #Jim's lettuce revision
    living_total_mass(living_shoot_mass, living_root_mass) => begin
        living_shoot_mass + living_root_mass
    end ~ track(u"g")

    # this will only be used for total leaf area adjustment.
    # If the individual leaf thing works out this will be deleted.
    potential_carbon_demand(potential_leaf_area_increase, SLA=200u"cm^2/g") => begin
        # Just a mocking value for now. Need to find a more mechanistic way to simulate change in SLA YY
        # SK 8/20/10: changed it to 200 cm2/g based on data from Kim et al. (2007) EEB
        # units are biomass not carbon
        leaf_mass_demand = potential_leaf_area_increase / SLA
        # potential_carbon_demand = carbon_demand # for now only carbon demand for leaf is calculated.
        leaf_mass_demand
    end ~ track(u"g")
end


text:
@testset "Links" begin

  #Logit link
  @test isequal(Unmarked.invlink(0.5, Unmarked.LogitLink()), logistic(0.5))
  @test isequal(Unmarked.invlink([0.3,0.5], Unmarked.LogitLink()),
                logistic.([0.3,0.5]))

  #Log link
  @test isequal(Unmarked.invlink(0.5, Unmarked.LogLink()), exp(0.5))
  @test isequal(Unmarked.invlink([0.3,0.5], Unmarked.LogLink()),
                exp.([0.3,0.5]))

end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


text:
using StaticArrays, LinearAlgebra

"""
    AbstractAxis

An axis of a [`AxisTensor`](@ref).
"""
abstract type AbstractAxis{I} end

Base.Broadcast.broadcastable(a::AbstractAxis) = a

"""
    dual(ax::AbstractAxis)

The dual axis to `ax`.

```
julia> using ClimaCore.Geometry

julia> Geometry.dual(Geometry.Covariant12Axis())
ClimaCore.Geometry.ContravariantAxis{(1, 2)}()

julia> Geometry.dual(Geometry.Cartesian123Axis())
ClimaCore.Geometry.CartesianAxis{(1, 2, 3)}()
```
"""
function dual end

struct CovariantAxis{I} <: AbstractAxis{I} end
symbols(::CovariantAxis) = (:u₁, :u₂, :u₃)

struct ContravariantAxis{I} <: AbstractAxis{I} end
symbols(::ContravariantAxis) = (:u¹, :u², :u³)
dual(::CovariantAxis{I}) where {I} = ContravariantAxis{I}()
dual(::ContravariantAxis{I}) where {I} = CovariantAxis{I}()

struct LocalAxis{I} <: AbstractAxis{I} end
symbols(::LocalAxis) = (:u, :v, :w)
dual(::LocalAxis{I}) where {I} = LocalAxis{I}()

struct CartesianAxis{I} <: AbstractAxis{I} end
symbols(::CartesianAxis) = (:u1, :u2, :u3)
dual(::CartesianAxis{I}) where {I} = CartesianAxis{I}()

coordinate_axis(::Type{<:XPoint}) = (1,)
coordinate_axis(::Type{<:YPoint}) = (2,)
coordinate_axis(::Type{<:ZPoint}) = (3,)
coordinate_axis(::Type{<:XYPoint}) = (1, 2)
coordinate_axis(::Type{<:XZPoint}) = (1, 3)
coordinate_axis(::Type{<:XYZPoint}) = (1, 2, 3)

coordinate_axis(::Type{<:Cartesian1Point}) = (1,)
coordinate_axis(::Type{<:Cartesian2Point}) = (2,)
coordinate_axis(::Type{<:Cartesian3Point}) = (3,)

coordinate_axis(::Type{<:Cartesian123Point}) = (1, 2, 3)

coordinate_axis(::Type{<:LatLongPoint}) = (1, 2)

coordinate_axis(coord::AbstractPoint) = coordinate_axis(typeof(coord))

@inline function idxin(I::Tuple{Vararg{Int}}, i::Int)
    N = length(I)
    for n in 1:N
        if I[n] == i
            return n
        end
    end
    return nothing
end

@inline function symidx(ax::AbstractAxis{I}, name::Symbol) where {I}
    S = symbols(ax)
    name == S[1] ? idxin(I, 1) :
    name == S[2] ? idxin(I, 2) :
    name == S[3] ? idxin(I, 3) : error("$ax has no symbol $name")
end

# most of these are required for printing
Base.length(ax::AbstractAxis{I}) where {I} = length(I)
Base.unitrange(ax::AbstractAxis) = StaticArrays.SOneTo(length(ax))
Base.LinearIndices(axes::Tuple{Vararg{AbstractAxis}}) =
    1:prod(map(length, axes))
Base.checkindex(::Type{Bool}, ax::AbstractAxis, i) =
    Base.checkindex(Bool, Base.unitrange(ax), i)
Base.lastindex(ax::AbstractAxis) = length(ax)
Base.getindex(m::AbstractAxis, i::Int) = i

"""
    AxisTensor(axes, components)

An `AxisTensor` is a wrapper around a `StaticArray`, where each dimension is
labelled with an [`AbstractAxis`](@ref). These axes must be consistent for
operations such as addition or subtraction, or be dual to each other for
operations such as multiplication.

# See also
[`components`](@ref) to obtain the underlying array.
"""
struct AxisTensor{
    T,
    N,
    A <: NTuple{N, AbstractAxis},
    S <: StaticArray{<:Tuple, T, N},
} <: AbstractArray{T, N}
    axes::A
    components::S
end

AxisTensor(
    axes::A,
    components::S,
) where {
    A <: Tuple{Vararg{AbstractAxis}},
    S <: StaticArray{<:Tuple, T, N},
} where {T, N} = AxisTensor{T, N, A, S}(axes, components)

AxisTensor(axes::Tuple{Vararg{AbstractAxis}}, components) =
    AxisTensor(axes, SArray{Tuple{map(length, axes)...}}(components))

# if the axes are already defined
(AxisTensor{T, N, A, S} where {S})(
    components::AbstractArray{T, N},
) where {T, N, A} = AxisTensor(A.instance, components)
(AxisTensor{T, N, A, S} where {T, S})(
    components::AbstractArray{<:Any, N},
) where {N, A} = AxisTensor(A.instance, components)

Base.axes(a::AxisTensor) = getfield(a, :axes)
Base.axes(::Type{AxisTensor{T, N, A, S}}) where {T, N, A, S} = A.instance
Base.size(a::AxisTensor) = map(length, axes(a))

"""
    components(a::AxisTensor)

Returns a `StaticArray` containing the components of `a` in its stored basis.
"""
components(a::AxisTensor) = getfield(a, :components)

Base.getindex(v::AxisTensor, i...) = getindex(components(v), i...)

function Base.getindex(
    v::AxisTensor{<:Any, 2, Tuple{A1, A2}},
    ::Colon,
    i::Integer,
) where {A1, A2}
    AxisVector(axes(v, 1), getindex(components(v), :, i))
end
function Base.getindex(
    v::AxisTensor{<:Any, 2, Tuple{A1, A2}},
    i::Integer,
    ::Colon,
) where {A1, A2}
    AxisVector(axes(v, 2), getindex(components(v), i, :))
end

Base.map(f::F, a::AxisTensor) where {F} =
    AxisTensor(axes(a), map(f, components(a)))
Base.map(
    f::F,
    a::AxisTensor{Ta, N, A},
    b::AxisTensor{Tb, N, A},
) where {F, Ta, Tb, N, A} =
    AxisTensor(axes(a), map(f, components(a), components(b)))
#Base.map(f, a::AxisTensor{Ta,N}, b::AxisTensor{Tb,N}) where {Ta,Tb,N} =
#    map(f, promote(a,b)...)

Base.zero(a::AxisTensor{T, N, A, S}) where {T, N, A, S} = zero(typeof(a))
Base.zero(::Type{AxisTensor{T, N, A, S}}) where {T, N, A, S} =
    AxisTensor(axes(AxisTensor{T, N, A, S}), zero(S))

import Base: +, -, *, /, \, ==

# Unary ops
@inline +(a::AxisTensor) = map(+, a)
@inline -(a::AxisTensor) = map(-, a)

# Binary ops
# Between arrays
@inline +(a::AxisTensor, b::AxisTensor) = map(+, a, b)
@inline -(a::AxisTensor, b::AxisTensor) = map(-, a, b)
# Scalar-array
@inline *(a::Number, b::AxisTensor) = map(c -> a * c, b)
@inline *(a::AxisTensor, b::Number) = map(c -> c * b, a)
@inline /(a::AxisTensor, b::Number) = map(c -> c / b, a)
@inline \(a, b::AxisTensor) = map(c -> a \ c, b)

@inline (==)(a::AxisTensor, b::AxisTensor) =
    axes(a) == axes(b) && components(a) == components(b)

# vectors
const AxisVector{T, A1, S} = AxisTensor{T, 1, Tuple{A1}, S}

AxisVector(ax::A1, v::SVector{N, T}) where {A1 <: AbstractAxis, N, T} =
    AxisVector{T, A1, SVector{N, T}}((ax,), v)

(AxisVector{T, A, SVector{1, T}} where {T})(arg1::Real) where {A} =
    AxisVector(A.instance, SVector(arg1))
(AxisVector{T, A, SVector{2, T}} where {T})(arg1::Real, arg2::Real) where {A} =
    AxisVector(A.instance, SVector(arg1, arg2))
(AxisVector{T, A, SVector{3, T}} where {T})(
    arg1::Real,
    arg2::Real,
    arg3::Real,
) where {A} = AxisVector(A.instance, SVector(arg1, arg2, arg3))

const CovariantVector{T, I, S} = AxisVector{T, CovariantAxis{I}, S}
const ContravariantVector{T, I, S} = AxisVector{T, ContravariantAxis{I}, S}
const CartesianVector{T, I, S} = AxisVector{T, CartesianAxis{I}, S}
const LocalVector{T, I, S} = AxisVector{T, LocalAxis{I}, S}

Base.propertynames(x::AxisVector) = symbols(axes(x, 1))
function Base.getproperty(x::AxisVector, name::Symbol)
    n = symidx(axes(x, 1), name)
    isnothing(n) ? zero(eltype(x)) : components(x)[n]
end

const AdjointAxisVector{T, A1, S} = Adjoint{T, AxisVector{T, A1, S}}

components(va::AdjointAxisVector) = components(parent(va))'
Base.getindex(va::AdjointAxisVector, i::Int) = getindex(components(va), i)
Base.getindex(va::AdjointAxisVector, i::Int, j::Int) =
    getindex(components(va), i, j)

# 2-tensors
const Axis2Tensor{T, A, S} = AxisTensor{T, 2, A, S}
Axis2Tensor(
    axes::Tuple{AbstractAxis, AbstractAxis},
    components::AbstractMatrix,
) = AxisTensor(axes, components)

const AdjointAxis2Tensor{T, A, S} = Adjoint{T, Axis2Tensor{T, A, S}}
components(va::AdjointAxis2Tensor) = components(parent(va))'

const Axis2TensorOrAdj{T, A, S} =
    Union{Axis2Tensor{T, A, S}, AdjointAxis2Tensor{T, A, S}}

# based on 1st dimension
const Covariant2Tensor{T, A, S} =
    Axis2Tensor{T, A, S} where {A <: Tuple{CovariantAxis, AbstractAxis}}
const Contravariant2Tensor{T, A, S} =
    Axis2Tensor{T, A, S} where {A <: Tuple{ContravariantAxis, AbstractAxis}}
const Cartesian2Tensor{T, A, S} =
    Axis2Tensor{T, A, S} where {A <: Tuple{CartesianAxis, AbstractAxis}}
const Local2Tensor{T, A, S} =
    Axis2Tensor{T, A, S} where {A <: Tuple{LocalAxis, AbstractAxis}}

const CovariantTensor = Union{CovariantVector, Covariant2Tensor}
const ContravariantTensor = Union{ContravariantVector, Contravariant2Tensor}
const CartesianTensor = Union{CartesianVector, Cartesian2Tensor}
const LocalTensor = Union{LocalVector, Local2Tensor}

for I in [(1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
    strI = join(I)
    N = length(I)
    strUVW = join(map(i -> [:U, :V, :W][i], I))
    @eval begin
        const $(Symbol(:Covariant, strI, :Axis)) = CovariantAxis{$I}
        const $(Symbol(:Covariant, strI, :Vector)){T} =
            CovariantVector{T, $I, SVector{$N, T}}

        const $(Symbol(:Contravariant, strI, :Axis)) = ContravariantAxis{$I}
        const $(Symbol(:Contravariant, strI, :Vector)){T} =
            ContravariantVector{T, $I, SVector{$N, T}}

        const $(Symbol(:Cartesian, strI, :Axis)) = CartesianAxis{$I}
        const $(Symbol(:Cartesian, strI, :Vector)){T} =
            CartesianVector{T, $I, SVector{$N, T}}

        const $(Symbol(strUVW, :Axis)) = LocalAxis{$I}
        const $(Symbol(strUVW, :Vector)){T} = LocalVector{T, $I, SVector{$N, T}}
    end
end

# LinearAlgebra

check_axes(::A, ::A) where {A} = nothing
check_axes(ax1, ax2) where {A} =
    throw(DimensionMismatch("$ax1 and $ax2 do not match"))

check_dual(ax1, ax2) = _check_dual(ax1, ax2, dual(ax2))
_check_dual(::A, _, ::A) where {A} = nothing
_check_dual(ax1, ax2, _) =
    throw(DimensionMismatch("$ax1 is not dual with $ax2"))

function LinearAlgebra.dot(x::AxisVector, y::AxisVector)
    check_dual(axes(x, 1), axes(y, 1))
    return LinearAlgebra.dot(components(x), components(y))
end

function Base.:*(x::AxisVector, y::AdjointAxisVector)
    AxisTensor((axes(x, 1), axes(y, 2)), components(x) * components(y))
end
function Base.:*(A::Axis2TensorOrAdj, x::AxisVector)
    check_dual(axes(A, 2), axes(x, 1))
    return AxisVector(axes(A, 1), components(A) * components(x))
end
function Base.:*(A::Axis2TensorOrAdj, B::Axis2TensorOrAdj)
    check_dual(axes(A, 2), axes(B, 1))
    return AxisTensor((axes(A, 1), axes(B, 2)), components(A) * components(B))
end

function Base.inv(A::Axis2TensorOrAdj)
    return AxisTensor((dual(axes(A, 2)), dual(axes(A, 1))), inv(components(A)))
end
function Base.:\(A::Axis2TensorOrAdj, x::AxisVector)
    check_axes(axes(A, 1), axes(x, 1))
    return AxisVector(dual(axes(A, 2)), components(A) \ components(x))
end

# can only compute norm if self-dual
function LinearAlgebra.norm(x::AxisVector)
    check_dual(axes(x, 1), axes(x, 1))
    LinearAlgebra.norm(components(x))
end
function LinearAlgebra.norm_sqr(x::AxisVector)
    check_dual(axes(x, 1), axes(x, 1))
    LinearAlgebra.norm_sqr(components(x))
end

LinearAlgebra.cross(x::Cartesian12Vector, y::Cartesian3Vector) =
    Cartesian12Vector(x.u2 * y.u3, -x.u1 * y.u3)
LinearAlgebra.cross(y::Cartesian3Vector, x::Cartesian12Vector) =
    Cartesian12Vector(-x.u2 * y.u3, x.u1 * y.u3)
LinearAlgebra.cross(x::UVVector, y::WVector) = UVVector(x.v * y.w, -x.u * y.w)
LinearAlgebra.cross(y::WVector, x::UVVector) = UVVector(-x.v * y.w, x.u * y.w)

function Base.:(+)(A::Axis2Tensor, b::LinearAlgebra.UniformScaling)
    check_dual(axes(A)...)
    AxisTensor(axes(A), components(A) + b)
end
function Base.:(-)(A::Axis2Tensor, b::LinearAlgebra.UniformScaling)
    check_dual(axes(A)...)
    AxisTensor(axes(A), components(A) - b)
end

function _transform(
    ato::Ato,
    x::AxisVector{T, Afrom, SVector{N, T}},
) where {Ato <: AbstractAxis{I}, Afrom <: AbstractAxis{I}} where {I, T, N}
    x
end

function _project(
    ato::Ato,
    x::AxisVector{T, Afrom, SVector{N, T}},
) where {Ato <: AbstractAxis{I}, Afrom <: AbstractAxis{I}} where {I, T, N}
    x
end

@generated function _transform(
    ato::Ato,
    x::AxisVector{T, Afrom, SVector{N, T}},
) where {
    Ato <: AbstractAxis{Ito},
    Afrom <: AbstractAxis{Ifrom},
} where {Ito, Ifrom, T, N}
    errcond = false
    for n in 1:N
        i = Ifrom[n]
        if i ∉ Ito
            errcond = :($errcond || x[$n] != zero(T))
        end
    end
    vals = []
    for i in Ito
        val = :(zero(T))
        for n in 1:N
            if i == Ifrom[n]
                val = :(x[$n])
                break
            end
        end
        push!(vals, val)
    end
    quote
        if $errcond
            throw(InexactError(:transform, Ato, x))
        end
        AxisVector(ato, SVector($(vals...)))
    end
end

@generated function _project(
    ato::Ato,
    x::AxisVector{T, Afrom, SVector{N, T}},
) where {
    Ato <: AbstractAxis{Ito},
    Afrom <: AbstractAxis{Ifrom},
} where {Ito, Ifrom, T, N}
    vals = []
    for i in Ito
        val = :(zero(T))
        for n in 1:N
            if i == Ifrom[n]
                val = :(x[$n])
                break
            end
        end
        push!(vals, val)
    end
    return :(AxisVector(ato, SVector($(vals...))))
end

function _transform(
    ato::Ato,
    x::Axis2Tensor{T, Tuple{Afrom, A2}},
) where {
    Ato <: AbstractAxis{I},
    Afrom <: AbstractAxis{I},
    A2 <: AbstractAxis{J},
} where {I, J, T}
    x
end

function _project(
    ato::Ato,
    x::Axis2Tensor{T, Tuple{Afrom, A2}},
) where {
    Ato <: AbstractAxis{I},
    Afrom <: AbstractAxis{I},
    A2 <: AbstractAxis{J},
} where {I, J, T}
    x
end

@generated function _transform(
    ato::Ato,
    x::Axis2Tensor{T, Tuple{Afrom, A2}},
) where {
    Ato <: AbstractAxis{Ito},
    Afrom <: AbstractAxis{Ifrom},
    A2 <: AbstractAxis{J},
} where {Ito, Ifrom, J, T}
    N = length(Ifrom)
    M = length(J)
    errcond = false
    for n in 1:N
        i = Ifrom[n]
        if i ∉ Ito
            for m in 1:M
                errcond = :($errcond || x[$n, $m] != zero(T))
            end
        end
    end
    vals = []
    for m in 1:M
        for i in Ito
            val = :(zero(T))
            for n in 1:N
                if i == Ifrom[n]
                    val = :(x[$n, $m])
                    break
                end
            end
            push!(vals, val)
        end
    end
    quote
        if $errcond
            throw(InexactError(:transform, Ato, x))
        end
        Axis2Tensor((ato, axes(x, 2)), SMatrix{$(length(Ito)), $M}($(vals...)))
    end
end

@generated function _project(
    ato::Ato,
    x::Axis2Tensor{T, Tuple{Afrom, A2}},
) where {
    Ato <: AbstractAxis{Ito},
    Afrom <: AbstractAxis{Ifrom},
    A2 <: AbstractAxis{J},
} where {Ito, Ifrom, J, T}
    N = length(Ifrom)
    M = length(J)
    vals = []
    for m in 1:M
        for i in Ito
            val = :(zero(T))
            for n in 1:N
                if i == Ifrom[n]
                    val = :(x[$n, $m])
                    break
                end
            end
            push!(vals, val)
        end
    end
    return :(Axis2Tensor(
        (ato, axes(x, 2)),
        SMatrix{$(length(Ito)), $M}($(vals...)),
    ))
end

transform(ato::CovariantAxis, v::CovariantTensor) = _transform(ato, v)
transform(ato::ContravariantAxis, v::ContravariantTensor) = _transform(ato, v)
transform(ato::CartesianAxis, v::CartesianTensor) = _transform(ato, v)
transform(ato::LocalAxis, v::LocalTensor) = _transform(ato, v)

project(ato::CovariantAxis, v::CovariantTensor) = _project(ato, v)
project(ato::ContravariantAxis, v::ContravariantTensor) = _project(ato, v)
project(ato::CartesianAxis, v::CartesianTensor) = _project(ato, v)
project(ato::LocalAxis, v::LocalTensor) = _project(ato, v)

"""
    outer(x, y)
    x ⊗ y

Compute the outer product of `x` and `y`. Typically `x` will be a vector, and
`y` can be either a number, vector or tuple/named tuple.

```julia
# vector ⊗ scalar = vector
julia> [1.0,2.0] ⊗ 2.0
2-element Vector{Float64}:
 2.0
 4.0

# vector ⊗ vector = matrix
julia> [1.0,2.0] ⊗ [1.0,3.0]
2×2 Matrix{Float64}:
 1.0  3.0
 2.0  6.0

# vector ⊗ tuple = recursion
julia> [1.0,2.0] ⊗ (1.0, (a=2.0, b=3.0))
([1.0, 2.0], (a = [2.0, 4.0], b = [3.0, 6.0]))
```
"""
function outer end
const ⊗ = outer

function outer(x::AbstractVector, y::AbstractVector)
    x * y'
end
function outer(x::AbstractVector, y::Number)
    x * y
end
function outer(x::AbstractVector, y)
    RecursiveApply.rmap(y -> x ⊗ y, y)
end


text:
"""
    _variable_dbbt!

Tightens the bounds of the `_current_node` using the current global upper bound
and the duality information obtained from the relaxation.
"""
function _variable_dbbt!(m::GlobalOptimizer{N,T,S}) where {N,T<:AbstractFloat,S}

    LBD = m._lower_objective_value
    UBD = m._global_upper_bound
    delta = UBD - LBD
    if LBD <= UBD
        @inbounds for i = 1:m._branch_variable_num
            ml = m._lower_lvd[i]
            if ml > zero(T)
                cut = _lower_bound(BranchVar,m,i) + delta/ml
                if cut < _upper_bound(BranchVar,m,i)
                    _set_upper_bound!(BranchVar,m,i,cut)
                end
            else
                mh = m._lower_uvd[i]
                if mh > zero(T)
                    cut = _upper_bound(BranchVar,m,i) - delta/mh
                    if cut > _lower_bound(BranchVar,m,i)
                        _set_lower_bound!(BranchVar,m,i,cut)
                    end
                end
            end
        end
    end
    return nothing
end

"""
$(SIGNATURES)

Default postprocess perfoms duality-based bound tightening only.
"""
function postprocess!(t::ExtensionType, m::GlobalOptimizer)
    (m.dbbt_depth > m._iteration_count) && _variable_dbbt!(m)
    return nothing
end


text:
using GDF_test
using Test

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


text:
# This file includes the template plugin

struct TemplatePlugin <: AbstractPlugin
    process(plg::TemplatePlugin, x) = println("In the template plugin. Doing nothing")
end


text:
using Plots
using Random

Random.seed!(123)

# Solver options
tf = 5.0
N = 101
integration = :rk4
opts = SolverOptions()
opts.verbose = true
opts.cost_tolerance = 1e-8
opts.cost_tolerance_intermediate = 1e-8
opts.constraint_tolerance = 1e-8
opts.constraint_tolerance_intermediate = 1e-8
opts.gradient_norm_tolerance = 1e-12
opts.gradient_norm_tolerance_intermediate = 1e-12
opts.use_gradient_aula = true
opts.active_constraint_tolerance = 0.0
opts.penalty_scaling = 10.0
opts.penalty_initial = 1e-3
opts.penalty_update_frequency = 1
opts.constraint_decrease_ratio = .25
opts.iterations = 1000
opts.iterations_outerloop = 50
opts.iterations_innerloop = 300
opts.outer_loop_update_type = :individual

# Set up model, objective, solver
model, = TrajectoryOptimization.Dynamics.quadrotor
n = model.n
m = model.m

# -initial state
x0 = zeros(n)
x0[1:3] = [0.; 0.; 0.]
q0 = [1.;0.;0.;0.]
x0[4:7] = q0

# -final state
xf = copy(x0)
xf[1:3] = [0.;40.;0.] # xyz position
xf[4:7] = q0

# -control limits
u_min = 0.0
u_max = 10.0

# Q = (1e-1)*Matrix(I,n,n)
# Q[4,4] = 1.0
# Q[5,5] = 1.0
# Q[6,6] = 1.0
# Q[7,7] = 1.0
# R = (1e-1)*Matrix(I,m,m)
Q = (1e-1)*Matrix(I,n,n)
R = (1e-1)*Matrix(I,m,m)
Qf = (1000.0)*Matrix(I,n,n)

# obstacles constraint
r_quad = 3.0
r_sphere = 3.0
spheres = ((0.,10.,0.,r_sphere),(0.,20.,0.,r_sphere),(0.,30.,0.,r_sphere))
n_spheres = 3

function cI(c,x,u)
    for i = 1:n_spheres
        c[i] = sphere_constraint(x,spheres[i][1],spheres[i][2],spheres[i][3],spheres[i][4]+r_quad)
    end
    c
end

# unit quaternion constraint
function cE(c,x,u)
    c = sqrt(x[4]^2 + x[5]^2 + x[6]^2 + x[7]^2) - 1.0
end

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
obj_con = TrajectoryOptimization.ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max,cI=cI,cE=cE)

solver_uncon = Solver(model,obj_uncon,integration=integration,N=N,opts=opts)
solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)

U0 = zeros(solver_uncon.model.m, solver_uncon.N-1)
@time results_uncon, stats_uncon = solve(solver_uncon,U0)

plot(to_array(results_uncon.X)[1:3,:]')

solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)

results_con.μ

@time results_con, stats_con = solve(solver_con,U0)

plot(to_array(results_con.U)[:,1:solver_con.N-1]')
plot(to_array(results_con.X)[1:3,:]')
plot(to_array(results_con.λ[1:N-1])')

max_violation(results_con)
total_time(solver_con,results_con)

plot!(log.(stats_con["max_condition_number"]))
plot(log.(stats_con["c_max"]).+15.5)
plot(log.(stats_con["cost"]))


text:
module Const

const kByteCountMask  = Int64(0x40000000)
const kByteCountVMask = Int64(0x4000)
const kClassMask      = Int64(0x80000000)
const kNewClassTag    = Int64(0xFFFFFFFF)

const kIsOnHeap       = UInt32(0x01000000)
const kIsReferenced   = UInt32(16)

const kMapOffset      = 2

const kNullTag              = 0
const kNotDeleted           = UInt32(0x02000000)
const kZombie               = UInt32(0x04000000)
const kBitMask              = UInt32(0x00FFFFFF)
const kDisplacementMask     = UInt32(0xFF000000)

# core/zip/inc/Compression.h
const kZLIB                 = 1
const kLZMA                 = 2
const kOldCompressionAlgo   = 3
const kLZ4                  = 4
const kZSTD                 = 5
const kUndefinedCompressionAlgorithm = 6

# constants for streamers
const kBase                 = 0
const kChar                 = 1
const kShort                = 2
const kInt                  = 3
const kLong                 = 4
const kFloat                = 5
const kCounter              = 6
const kCharStar             = 7
const kDouble               = 8
const kDouble32             = 9
const kLegacyChar           = 10
const kUChar                = 11
const kUShort               = 12
const kUInt                 = 13
const kULong                = 14
const kBits                 = 15
const kLong64               = 16
const kULong64              = 17
const kBool                 = 18
const kFloat16              = 19
const kOffsetL              = 20
const kOffsetP              = 40
const kObject               = 61
const kAny                  = 62
const kObjectp              = 63
const kObjectP              = 64
const kTString              = 65
const kTObject              = 66
const kTNamed               = 67
const kAnyp                 = 68
const kAnyP                 = 69
const kAnyPnoVT             = 70
const kSTLp                 = 71

const kSkip                 = 100
const kSkipL                = 120
const kSkipP                = 140

const kConv                 = 200
const kConvL                = 220
const kConvP                = 240

const kSTL                  = 300
const kSTLstring            = 365

const kStreamer             = 500
const kStreamLoop           = 501

# TBranchElement fTypes
# https://groups.google.com/d/msg/polyglot-root-io/yeC0mAizQcA/zuUHOFBABwAJ

const kTopLevelTClonesArray   = 3
const kSubbranchTClonesArray  = 31
const kTopLevelSTLCollection  = 4
const kSubbranchSTLCollection = 41

# constants from core/foundation/inc/ESTLType.h

const kNotSTL               = 0
const kSTLvector            = 1
const kSTLlist              = 2
const kSTLdeque             = 3
const kSTLmap               = 4
const kSTLmultimap          = 5
const kSTLset               = 6
const kSTLmultiset          = 7
const kSTLbitset            = 8
const kSTLforwardlist       = 9
const kSTLunorderedset      = 10
const kSTLunorderedmultiset = 11
const kSTLunorderedmap      = 12
const kSTLunorderedmultimap = 13
const kSTLend               = 14
const kSTLany               = 300

# IOFeatures

const kGenerateOffsetMap    = 1

end


text:
"""
    NoOp <: Augmentor.AffineOperation

Identity transformation that does not do anything with the given
image, but instead passes it along unchanged (without copying).

Usually used in combination with [`Either`](@ref) to denote a
"branch" that does not perform any computation.
"""
struct NoOp <: AffineOperation end

@inline supports_eager(::Type{NoOp}) = false
@inline supports_stepview(::Type{NoOp}) = true
@inline supports_view(::Type{NoOp}) = true

# TODO: implement method for n-dim arrays
toaffinemap(::NoOp, img::AbstractMatrix) = AffineMap(@SMatrix([1. 0; 0 1.]), @SVector([0.,0.]))
applyeager(::NoOp, img::AbstractArray, param) = plain_array(img)
applylazy(::NoOp, img::AbstractArray, param) = img

function applyview(::NoOp, img::AbstractArray, param)
    idx = map(i->1:length(i), axes(img))
    indirect_view(img, idx)
end

function applystepview(::NoOp, img::AbstractArray, param)
    idx = map(i->1:1:length(i), axes(img))
    indirect_view(img, idx)
end

function showconstruction(io::IO, op::NoOp)
    print(io, typeof(op).name.name, "()")
end

function Base.show(io::IO, op::NoOp)
    if get(io, :compact, false)
        print(io, "No operation")
    else
        print(io, "Augmentor.")
        showconstruction(io, op)
    end
end


text:
include("definitions.jl")

function multifilter(;conditions...)
    f = (x) -> all(map((k, v) -> x[k] == v, keys(conditions), values(conditions)))
    return f
end

function timeify(results::AbstractDataFrame, var; id = :ID)
    gdf = groupby(results, [id])
    df = combine(gdf, var => x -> [x])
    return df
end

function bins(var::Number, edges, values; minimum = missing)
    @assert (length(edges) == length(values)) "Mismatch in bin dimensions"
    bin = sum(var .>= edges)
    return bin == 0 ? minimum : values[bin]
end

function bins(vector::AbstractArray, edges, values; minimum = missing)
    return bins.(vector, Ref(edges), Ref(values), minimum = minimum)
end

function ss!(data, metric::Symbol, x0::Pair, xF::Pair; col::Union{Symbol, AbstractString} = :ss)
    initial = filter(PKPD.multifilter(;x0), data)[metric]
    final = filter(PKPD.multifilter(;xF), data)[metric]

    @assert !(col in names(data)) "Column already exists"
    @assert length(initial) == length(final) "Initial and finals vectors are of different lengths"
    data[!, col] = final - initial
    return nothing
end

function read(path::String)
    open(path) do file
        file_array = readlines(file)
        colnames = filter(x -> x != "", split(file_array[2], ' '))
        colnames = [Symbol(col) for col in colnames]
        data = DataFrame(fill(Float64[], length(colnames)), colnames; makeunique=true)
        for i in 3:length(file_array)
            row = file_array[i]
            if (row == file_array[1]) | (row == file_array[2])
                continue
            else
                row_vector = split(row, ' ')
                row_vector = filter(x -> x != "", row_vector)
                row_vector = map(x -> parse(Float64, x), row_vector)
                push!(data, row_vector)
            end
        end
        return data
    end
end


text:
"""
`module  JAC.InteractionStrengthQED`
    ... a submodel of JAC that contains methods for computing local single-electron QED corrections.
"""
module InteractionStrengthQED

    using Printf, ..Basics, ..Bsplines, ..Defaults, ..HydrogenicIon, ..ManyElectron, ..Radial, ..RadialIntegrals, ..Nuclear

    #                   Z/10       1s_1/2      2p_1/2      2p_3/2      3d_3/2      3d_5/2        Zi
    const   fez_qed =         [    4.6542     -0.1148      0.1304     -0.0427      0.0408   ;   #  10
                                   3.2462     -0.0925      0.1438     -0.0420      0.0417   ;   #  20
                                   2.5518     -0.0643      0.1606     -0.0410      0.0432   ;   #  30
                                   2.1347     -0.0310      0.1796     -0.0396      0.0452   ;   #  40
                                   1.8633      0.0080      0.2001     -0.0378      0.0475   ;   #  50
                                   1.6820      0.0547      0.2216     -0.0353      0.0503   ;   #  60
                                   1.5637      0.1126      0.2441     -0.0321      0.0536   ;   #  70
                                   1.4955      0.1877      0.2671     -0.0279      0.0572   ;   #  80
                                   1.4721      0.2912      0.2904     -0.0229      0.0612   ;   #  90
                                   1.4961      0.4450      0.3135     -0.0154      0.0654   ;   # 100
                                   1.5771      0.6961      0.3356     -0.0063      0.0699   ;   # 110
                                   1.7335      1.1559      0.3548      0.0051      0.0745   ]   # 120

    """
    `InteractionStrengthQED.qedLocal(a::Orbital, b::Orbital, nm::Nuclear.Model, qed::ManyElectron.AbstractQedModel,
                                     pot::Radial.Potential, grid::Radial.Grid)`
        ... to calculate the local, single-electron QED correction due to a chosen QedModel.  The basis::Basis is needed to
            compute the electron density for the Uehling potential. The function also determines whether
            the 'stored' hydrogenic values for the lambda_C-dampled overlap integrals belong the given Z-value, and re-calculates
            these overlap integrals whenever necessary.  A single-electron amplitud wa::Float64 is returned.
    """
    function qedLocal(a::Orbital, b::Orbital, nm::Nuclear.Model, qed::ManyElectron.AbstractQedModel,
                      pot::Radial.Potential, grid::Radial.Grid)
        # Define a grid for the t-integration
        qgrid = Radial.GridGL("QED",  7)
        if  a.subshell == b.subshell == Subshell("1s_1/2")    println(" ")    end
        # Re-define the hydrogenic values of the
        if  nm.Z != Defaults.GBL_QED_NUCLEAR_CHARGE
            println(">> Calculate hydrogenic values only if necessary.")
            alpha = Defaults.getDefaults("alpha")

            wp           = Bsplines.generatePrimitives(grid)
            subshellList = [Subshell("1s_1/2"), Subshell("2p_1/2"), Subshell("2p_3/2"), Subshell("3d_3/2"), Subshell("3d_5/2")]
            orbitals     = Bsplines.generateOrbitalsHydrogenic(wp, nm, subshellList; printout=true)
            wc1          = wc2 = wc3 = wc4 = wc5 = 1.0
            for  subsh in subshellList
                orb = orbitals[subsh]
                if      orb.subshell == Subshell("1s_1/2")      wc1 = RadialIntegrals.qedDampedOverlap(alpha, orb, orb, grid)
                elseif  orb.subshell == Subshell("2p_1/2")      wc2 = RadialIntegrals.qedDampedOverlap(alpha, orb, orb, grid)
                elseif  orb.subshell == Subshell("2p_3/2")      wc3 = RadialIntegrals.qedDampedOverlap(alpha, orb, orb, grid)
                elseif  orb.subshell == Subshell("3d_3/2")      wc4 = RadialIntegrals.qedDampedOverlap(alpha, orb, orb, grid)
                elseif  orb.subshell == Subshell("3d_5/2")      wc5 = RadialIntegrals.qedDampedOverlap(alpha, orb, orb, grid)
                end
            end
            Defaults.setDefaults("QED: damped-hydrogenic", nm.Z, [wc1, wc2, wc3, wc4, wc5] )
            println("Redefined damped radial integrals GBL_QED_HYDROGENIC_LAMBDAC = $(Defaults.GBL_QED_HYDROGENIC_LAMBDAC)")
        end

        if      qed == QedSydney()
            wa = RadialIntegrals.qedUehlingSimple(a, b, pot, grid, qgrid) +
                 ## RadialIntegrals.qedWichmannKrollSimple(a, b, pot, grid, qgrid) +
                 ## RadialIntegrals.qedElectricFormFactor(a, b, pot, grid, qgrid) +
                 ## RadialIntegrals.qedMagneticFormFactor(a, b, pot, grid, qgrid) +
                 RadialIntegrals.qedLowFrequency(a, b, nm, grid, qgrid)
        elseif   qed == QedPetersburg()
            ## RadialIntegrals.qedUehling(a, b, nm, grid, qgrid)
            wa = RadialIntegrals.qedUehlingSimple(a, b, pot, grid, qgrid) +
                 InteractionStrengthQED.selfEnergyVolotka(a, b, nm, grid, qgrid)
        else    error("stop a")
        end
        return( wa )
    end

    """
    `InteractionStrengthQED.selfEnergyVolotka(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid,
                                              qgrid::Radial.GridGL)`
        ... to calculate the local, single-electron self-energy contribution due to the Petersburg model (PRA, 2013), further
            simplified by Andrey Volotka. A non-zero self-energy contribution is returned only for orbitals with 1 <= n <= 4 and
            kappa = -1, 1, -2, 2, -3, since the method appears rather inaccurate for higher (n, kappa).
            A single-electron self-energy amplitude wa::Float64 is returned.
    """
    function selfEnergyVolotka(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid, qgrid::Radial.GridGL)
        # Non-zero local amplitudes only for kappa_a = kappa_b
        if      a.subshell != b.subshell  ||  a.subshell.n > 4         return( 0. )    end
        if      a.subshell.kappa == 3  ||  abs(a.subshell.kappa) > 3   return( 0. )    end

        if      a.subshell.kappa == -1   nn = 1;    elseif  a.subshell.kappa ==  1   nn = 2;    elseif  a.subshell.kappa == -2   nn = 2
        elseif  a.subshell.kappa ==  2   nn = 3;    elseif  a.subshell.kappa == -3   nn = 3
        else    error("stop a")
        end

        alpha       = Defaults.getDefaults("alpha");
        whydrogenic = InteractionStrengthQED.tabulateFzeOverHydrogenic(nm.Z, a.subshell) * alpha / pi * (alpha*nm.Z)^4 / nn^3 / alpha^2
        wb          = RadialIntegrals.qedDampedOverlap(alpha, a, a, grid)
        wa          = whydrogenic * wb

        ## println("<e^-r>, Coulomb   = $(wb)")
        ## println("<e^-r>, DH        = $(wb)")
        ## println("local SE, DHF     = $(wa)  [a.u.]")

        # Printout the corresponding  F(alpha Z) value for the given orbital if a == b
        ## FalphaZ = alpha/pi * (alpha*nm.Z)^4 / a.subshell.n^3 / alpha^2
        ## FalphaZ = wa / FalphaZ
        ## println("Self-energy function for $(a.subshell) and $(a.subshell) is F(alpha Z) = $FalphaZ,  wdamped = $wb,  wa = $wa ")

        println("QED single-electron strength <$(a.subshell)| h^(SE, Volotka) | $(b.subshell)> = $wa    with <$(a.subshell)| e^-alpha r | $(b.subshell)> = $wb")

        return( wa )
    end

    """
    `InteractionStrengthQED.tabulateFzeOverHydrogenic(Z::Float64, sh::Subshell)`
        ... to return the value of the FZE() table as shown by Shabaev et al., PRA 88, 012513  (2013) but just for the lowest
            possible subshell of a given kappa-value; this tabulated value is devided by the lambda_C-dampled overlap integral,
            calculated with hydrogenic functions; a Float64 is returned. The table below only includes the data for extended
            nuclei, if available in the tabulations above, while a point-nucleus was assumed otherwise. All these data
            are calculated with hydrogenic function and, hence, need to be 're-scaled' for many-electron ions and atoms
            (as realized for these local, single-electron amplitudes).
    """
    function tabulateFzeOverHydrogenic(Z::Float64, sh::Subshell)
        global  fez_qed
        # Select the proper column from the table
        if      sh.kappa == -1   k = 1;    elseif  sh.kappa ==  1   k = 2;    elseif  sh.kappa == -2   k = 3
        elseif  sh.kappa ==  2   k = 4;    elseif  sh.kappa == -3   k = 5
        else    error("stop a")
        end
        #  Select the proper raw from the table
        if        1.0 <= Z <  10.0    return( 0. )
        elseif   10.0 <= Z <  20.0    wa =  fez_qed[1,k] + (Z -  10.) / 10. *( fez_qed[2,k] -  fez_qed[1,k])
        elseif   20.0 <= Z <  30.0    wa =  fez_qed[2,k] + (Z -  20.) / 10. *( fez_qed[3,k] -  fez_qed[2,k])
        elseif   30.0 <= Z <  40.0    wa =  fez_qed[3,k] + (Z -  30.) / 10. *( fez_qed[4,k] -  fez_qed[3,k])
        elseif   40.0 <= Z <  50.0    wa =  fez_qed[4,k] + (Z -  40.) / 10. *( fez_qed[5,k] -  fez_qed[4,k])
        elseif   50.0 <= Z <  60.0    wa =  fez_qed[5,k] + (Z -  50.) / 10. *( fez_qed[6,k] -  fez_qed[5,k])
        elseif   60.0 <= Z <  70.0    wa =  fez_qed[6,k] + (Z -  60.) / 10. *( fez_qed[7,k] -  fez_qed[6,k])
        elseif   70.0 <= Z <  80.0    wa =  fez_qed[7,k] + (Z -  70.) / 10. *( fez_qed[8,k] -  fez_qed[7,k])
        elseif   80.0 <= Z <  90.0    wa =  fez_qed[8,k] + (Z -  80.) / 10. *( fez_qed[9,k] -  fez_qed[8,k])
        elseif   90.0 <= Z < 100.0    wa =  fez_qed[9,k] + (Z -  90.) / 10. *(fez_qed[10,k] -  fez_qed[9,k])
        elseif  100.0 <= Z < 110.0    wa = fez_qed[10,k] + (Z - 100.) / 10. *(fez_qed[11,k] - fez_qed[10,k])
        elseif  110.0 <= Z < 120.0    wa = fez_qed[11,k] + (Z - 110.) / 10. *(fez_qed[12,k] - fez_qed[11,k])
        else    error("stop c")
        end

        wb = Defaults.GBL_QED_HYDROGENIC_LAMBDAC[k]
        ## println("QED  <$(sh)| h^(SE, hydrogenic) | $(sh)> = $wa      e^-alpha r, hydrogenic = $wb")
        ## println("local SE, Coulomb = $(wa)    [F]")
        ## println("<e^-r>, Coulomb   = $(wb)    from storage")

        return( wa / Defaults.GBL_QED_HYDROGENIC_LAMBDAC[k] )
    end

end # module


text:
using Knet
libknet8handle = Libdl.dlopen(Knet.libknet8)

SIZE = 100000
ITER = 100000
x32 = KnetArray(rand(Float32,SIZE))
y32 = similar(x32)
s32 = rand(Float32)
x64 = KnetArray(rand(Float64,SIZE))
y64 = similar(x64)
s64 = rand(Float64)

function cuda01test(fname, jname=fname, o...)
    println(fname)
    fcpu = eval(parse(jname))
    f32 = Libdl.dlsym(libknet8handle, fname*"_32_01")
    @time cuda01rep(f32,s32,x32,y32)
    isapprox(Array(y32),fcpu(s32,Array(x32))) || warn("$fname 32")
    f64 = Libdl.dlsym(libknet8handle, fname*"_64_01")
    @time cuda01rep(f64,s64,x64,y64)
    isapprox(Array(y64),fcpu(s64,Array(x64))) || warn("$fname 64")
end

function cuda01rep{T}(f,s::T,x::KnetArray{T},y::KnetArray{T})
    n = Cint(length(y))
    for i=1:ITER
        ccall(f,Void,(Cint,T,Ptr{T},Ptr{T}),n,s,x,y)
    end
    Knet.@cuda(cudart,cudaDeviceSynchronize,())
    Knet.@cuda(cudart,cudaGetLastError,())
end

for f in Knet.cuda01
    isa(f,Tuple) || (f=(f,))
    cuda01test(f...)
end


text:
using Test

@testset "GraphicalLinearAlgebra" begin
#   include("GLA.jl")
#   include("StructuredGLA.jl")
  include("bayes.jl")
end


text:
"""
    struct Parity

Represents the parity of a quantum system, taking two possible values: `even` or `odd`.

The integer values that correspond to `even` and `odd` parity are `+1` and `-1`, respectively.
`Base.convert` can be used to convert integers into `Parity` values.

```jldoctest
julia> convert(Parity, 1)
even

julia> convert(Parity, -1)
odd
```
"""
struct Parity
    p::Bool
end

"""
    @p_str

A string macro to easily construct [`Parity`](@ref) values.

```jldoctest
julia> p"even"
even

julia> p"odd"
odd
```
"""
macro p_str(ps)
    if ps == "even"
        Parity(true)
    elseif ps == "odd"
        Parity(false)
    else
        throw(ArgumentError("Invalid parity string $(ps)"))
    end
end

function Base.convert(::Type{Parity}, i::I) where {I<:Integer}
    i == 1 && return p"even"
    i == -1 && return p"odd"
    throw(ArgumentError("Don't know how to convert $(i) to parity"))
end

Base.iseven(p::Parity) = p.p
Base.isodd(p::Parity) = !p.p
Base.isless(a::Parity, b::Parity) = isodd(a) && iseven(b)

Base.:*(a::Parity, b::Parity) = Parity(a == b)
Base.:^(p::Parity, i::I) where {I<:Integer} =
    p.p || iseven(i) ? Parity(true) : Parity(false)
Base.:-(p::Parity) = Parity(!p.p)

Base.show(io::IO, p::Parity) =
    write(io, iseven(p) ? "even" : "odd")

UnicodeFun.to_superscript(p::Parity) = iseven(p) ? "" : "ᵒ"

"""
    parity(object) -> Parity

Returns the parity of `object`.
"""
function parity end

export Parity, @p_str, parity


text:
using SRKGenerator

dx = 1/50
imin,imax,jmin,jmax = -12,1,-10,10
x = [-0.00575177,3.48225,-2.59676,-0.364729,1.14808,-0.783625,0.13798,-0.397231,0.668615,-0.299582,1.19174,0.107113,-0.0913619,1.58563,0.061024,-4.21259,0.605065,3.15894,-0.372408,2.31455,-1.22656,1.19005,-0.694235,-0.681141,0.97351,-0.454075,0.560372,-0.0796419,-1.82528,1.5744,0.64214,0.608841,1.6336,-2.09727,0.239234,0.224391,2.82355,-1.57304,-0.64245,-0.608839,-4.40931,4.06816,1.23893,-0.897945]
# 13.18955568047541

@time f1 = SRKGenerator.cpu_f(x,dx,imin,imax,jmin,jmax)
@time f2 = SRKGenerator.parallel_f(x,dx,imin,imax,jmin,jmax)
@time f2 = SRKGenerator.threaded_f(x,dx,imin,imax,jmin,jmax)

f1 == f2

using CUDArt

NLoptRandSeed = 0
parameter_minmax=5
max_eval=Int(1e8)
initCon = ones(44)
tol = 1e-2
ftol = 1e-15
tol2 = 1e-5
counterSteps=Int(1e5)
counterSteps2=Int(1e6)
initStepSize=[]
gpuEnabled=true
ptx_str  = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationWin.ptx")
cudaCores = 1664
initStepSize2=1e-6
outfile=""
constrain_c = true

N = 26
N2= 16
M = 44
count = [0]

## Script Start
x_L = -parameter_minmax*ones(M)
x_U =  parameter_minmax*ones(M)
g_L = -tol*ones(N)
g_U =  tol*ones(N)
sizei = length(imin:dx:imax)
sizej = length(jmin:dx:jmax)
totArea = (imax-imin)*(jmax-jmin)/(sizei*sizej)
if gpuEnabled
  iarr = convert(Vector{Float32},collect(imin:dx:imax))
  jarr = convert(Vector{Float32},collect(jmin:dx:jmax))
  numCards = length(cudaCores)
  totCores = sum(cudaCores)
  equalDiv = sizei*sizej÷totCores + 1
  portions = equalDiv*cudaCores
  if numCards > 1
    startIdx = [0 cumsum(portions)[1:end-1]]
  else
    startIdx = [0]
  end
  g_iarr = Vector{CUDArt.CudaArray{Float32,1}}(numCards)
  g_jarr = Vector{CUDArt.CudaArray{Float32,1}}(numCards)
  g_tmp = Vector{CUDArt.CudaArray{Int32,1}}(numCards)
  g_coefs = Vector{CUDArt.CudaArray{Float32,1}}(numCards)
  ans = Vector{Int32}(numCards)
  integrationFuncs = Vector{CUDArt.CuFunction}(numCards)
  CUDArt.init(devices(dev->true))
  devlist = devices(dev->true)
  streams = [(device(dev); Stream()) for dev in devlist]
  for i = 1:numCards
    device(i-1)
    g_iarr[i]  = CudaArray(iarr)
    g_jarr[i]  = CudaArray(jarr)
    g_tmp[i]   = CudaArray(Int32,cudaCores[i])
    md = CuModule(ptx_str,false)
    integrationFuncs[i] = CuFunction(md,"integration")
  end
end

eV = [1;1;1;1]

eval_f = (x,grad) -> SRKGenerator.f_maker(x,ans,integrationFuncs,cudaCores,numCards,g_coefs,g_iarr,g_jarr,sizei,sizej,equalDiv,startIdx,g_tmp,totArea,counterSteps,counterSteps2,outfile,gpuEnabled,count)
eval_g = (tmp,x,grad) -> SRKGenerator.g_maker(x,tmp,eV,counterSteps,counterSteps2,outfile,count)
eval_g_ineq = (tmp,x,grad) -> SRKGenerator.g_ineq_maker(x,tmp,eV,counterSteps,counterSteps2,outfile,count)

f3 = eval_f(x,0)

#=
### Now test some gradients

f = (x) -> SRKGenerator.cpu_f(x,dx,imin,imax,jmin,jmax)

using ForwardDiff

ForwardDiff.gradient(f,x)

using Calculus

Calculus.gradient(f,x)

gpu_f = (x) -> eval_f(x,0)

Calculus.gradient(gpu_f,x)
=#


text:
using WebSockets
import WebSockets:Response,
                  Request
using Dates
using Sockets

global LASTREQ = 0
global LASTWS = 0
global LASTMSG = 0
global LASTSERVER = 0

const CLOSEAFTER = Dates.Second(1800)
const HTTPPORT = 8080
const LOCALIP = string(Sockets.getipaddr())
const USERNAMES = Dict{String, WebSocket}()
const HTMLSTRING = read(joinpath(@__DIR__, "chat_explore.html"), String)

@info """
A chat server application. For each browser (tab) that connects,
an 'asyncronous function' aka 'coroutine' aka 'task' is started.

To use:
    - include("chat_explore.jl") in REPL
    - start a browser on the local ip address, e.g.: http://192.168.0.4:8080
    - inspect global variables starting with 'LAST' while the chat is running asyncronously

"""

# Since we are to access a websocket from outside
# it's own websocket handler coroutine, we need some kind of
# mutable container for storing references:
const WEBSOCKETS = Dict{WebSocket, Int}()

"""
Called by 'gatekeeper', this function will be running in a task while the
particular websocket is open. The argument is an open websocket.
Other instances of the function run in other tasks.
"""
function coroutine(thisws)
    global lastws = thisws
    push!(WEBSOCKETS, thisws => length(WEBSOCKETS) +1 )
    t1 = now() + CLOSEAFTER
    username = ""
    while now() < t1
        # This next call waits for a message to
        # appear on the socket. If there is none,
        # this task yields to other tasks.
        data, success = readguarded(thisws)
        !success && break
        global LASTMSG = msg = String(data)
        print("Received: $msg ")
        if username == ""
            username = approvedusername(msg, thisws)
            if username != ""
                println("from new user $username ")
                !writeguarded(thisws, username) && break
                println("Tell everybody about $username")
                foreach(keys(WEBSOCKETS)) do ws
                    writeguarded(ws, username * " enters chat")
                end
            else
                println(", username taken!")
                !writeguarded(thisws, "Username taken!") && break
            end
        else
            println("from $username ")
            distributemsg(username * ": " * msg, thisws)
            startswith(msg, "exit") && break
        end
    end
    exitmsg = username == "" ? "unknown" : username * " has left"
    distributemsg(exitmsg, thisws)
    println(exitmsg)
    # No need to close the websocket. Just clean up external references:
    removereferences(thisws)
    nothing
end

function removereferences(ws)
    haskey(WEBSOCKETS, ws) && pop!(WEBSOCKETS, ws)
    for (discardname, wsref) in USERNAMES
        if wsref === ws
            pop!(USERNAMES, discardname)
            break
        end
    end
    nothing
end

function approvedusername(msg, ws)
    !startswith(msg, "userName:") && return ""
    newname = msg[length("userName:") + 1:end]
    newname =="" && return ""
    haskey(USERNAMES, newname) && return ""
    push!(USERNAMES, newname => ws)
    newname
end

function distributemsg(msgout, not_to_ws)
    foreach(keys(WEBSOCKETS)) do ws
        if ws !== not_to_ws
            writeguarded(ws, msgout)
        end
    end
    nothing
end

"""
`Server => gatekeeper(Request, WebSocket) => coroutine(WebSocket)`

The gatekeeper makes it a little harder to connect with
malicious code. It inspects the request that was upgraded
to a a websocket.
"""
function gatekeeper(req, ws)
    global LASTREQ = req
    global LASTWS = ws
    orig = WebSockets.origin(req)
    if occursin(LOCALIP, orig)
        coroutine(ws)
    else
        @warn("Unauthorized websocket connection, $orig not approved by gatekeeper, expected $LOCALIP")
    end
    nothing
end

"Request to response. Response is the predefined HTML page with some javascript"
req2resp(req::Request) = HTMLSTRING |> Response

# The following lines disblle detail messages from spilling into the
# REPL. Remove the it to gain insight.
using Logging
import Logging.shouldlog
function shouldlog(::ConsoleLogger, level, _module, group, id)
    if _module == WebSockets.HTTP.Servers
        if level == Logging.Warn || level == Logging.Info
            return false
        else
            return true
        end
    else
        return true
    end
end

# ServerWS takes two functions; the first a http request handler function for page requests,
# one for opening websockets (which javascript in the HTML page will try to do)
global LASTSERVER = WebSockets.ServerWS(req2resp, gatekeeper)

# Start the server asyncronously, and stop it later
@async WebSockets.serve(LASTSERVER, LOCALIP, HTTPPORT)
@async begin
    println("HTTP server listening on $LOCALIP:$HTTPPORT for $CLOSEAFTER")
    sleep(CLOSEAFTER.value)
    println("Time out, closing down $HTTPPORT")
    put!(LASTSERVER.in, "I can send anything, you close")
    nothing
end
# for inspecting in REPL or Atom / Juno - update after starting some clients.
LASTWS
LASTSERVER.out
WEBSOCKETS
#take!(LASTSERVER.out)|>string
nothing


text:
using VLOptionsModeling
using Plots

# create an emprt contract set for now -
contract_set = Set{VLAbstractAsset}()
assetSymbol = "XYZ"
strikePrice = 50.0
underlyingPrice = 50.0
expirationDate = nothing
premimumValue = 0.0
numberOfContracts = 1
sense = :buy
contractMultiplier = 1.0
put_contract = VLPutOptionContract(assetSymbol, strikePrice, premimumValue; sense=sense, contractMultiplier=contractMultiplier)
push!(contract_set, put_contract)

# setup the lattice (Example 19.1/Fig. 19.3 Hull)
𝝙t = 30.42 / (365.0) # years
μ = log(1.0084) * (1 / 𝝙t)
σ = (1 / √𝝙t) * log(1.1224)
numberOfLevels = 6

# Create the lattice model -
latticeModel = VLBinomialLattice(μ, σ, 𝝙t, numberOfLevels)

# compute delta -
delta_value = 𝝙(contract_set, latticeModel, underlyingPrice) |> check


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/qr_utilities.jl
# Functions used for preconditioning with QR factorization.
#############################################################################

#=
"""
QRDenseStorage

Provides preallocated storage for the QR factorization, Q, and the inverse of Q.

$(TYPEDFIELDS)
"""
mutable struct QRDenseStorage
    "QR Factorization"
    factorization::LinearAlgebra.QR{Float64,Array{Float64,2}}
    "Orthogonal matrix Q"
    Q::Array{Float64,2}
    "Inverse of Q"
    inv::Array{Float64,2}
end

"""
QRDenseStorage(nx::Int)

A constructor for QRDenseStorage assumes `Q` is of size `nx`-by-`nx` and of
type `Float64`.
"""
function QRDenseStorage(nx::Int)
    A = Float64.(Matrix(I, nx, nx))
    factorization = LinearAlgebra.qrfactUnblocked!(A)
    Q = similar(A)
    inverse = similar(A)
    QRDenseStorage(factorization, Q, inverse)
end

function Base.copy(q::QRDenseStorage)
    factors_c = q.factorization.factors
    τ_c = q.factorization.τ
    f_copy = LinearAlgebra.QR{Float64,Array{Float64,2}}(factors_c,τ_c)
    QRDenseStorage(f_copy, copy(q.Q), copy(q.inv))
end
=#

#=
"""
calculateQ!

Computes the QR factorization of `A` of size `(nx,nx)` and then stores it to
fields in `qst`.
"""
function calculateQ!(qst::QRDenseStorage, A::Matrix{Float64}, nx::Int)
    qst.factorization = LinearAlgebra.qrfactUnblocked!(A)
    qst.Q .= qst.factorization.Q*Matrix(I,nx,nx)
    nothing
end

"""
calculateQinv!

Computes `inv(Q)` via transpose! and stores this to `qst.inverse`.
"""
function calculateQinv!(qst::QRDenseStorage)
    transpose!(qst.inv, qst.Q)
    nothing
end
=#

function calculateQ!(Q, A::Matrix{Float64}, nx::Int)
    if nx == 1
        Q[1,1] = 1.0
    else
        F = LinearAlgebra.qr(A)
        Q = view(F.Q, 1:nx, 1:nx)
    end
    nothing
end
function calculateQinv!(Qinv, Q, nx::Int)
    if nx == 1
        Qinv[1,1] = 1.0
    else
        transpose!(Qinv, Q)
    end
    nothing
end
#=
"""
An circular buffer of fixed capacity and length which allows
for access via getindex and copying of an element to the last then cycling
the last element to the first and shifting all other elements. See
[DataStructures.jl](https://github.com/JuliaCollections/DataStructures.jl).
"""
function DataStructures.CircularBuffer(a::T, length::Int) where T
    cb = CircularBuffer{T}(length)
    append!(cb, [zero.(a) for i=1:length])
    cb
end

"""
reinitialize!(x::CircularBuffer{QRDenseStorage})

Sets the first QR storage to the identity matrix.
"""
function reinitialize!(x::CircularBuffer{QRDenseStorage})
    fill!(x[1].Q, 0.0)
    for i in 1:size(x[1].Q, 1)
        x[1].Q[i,i] = 1.0
    end
    nothing
end
=#


text:
# Optimal package install / uninstall
# http://www.cs.ucsd.edu/~lerner/papers/opium.pdf

import PicoSAT

index = Dict(
    "a" => Dict(:depends => ["b", "c", "z"]),
    "b" => Dict(:depends => ["d"]),
    "c" => Dict(:depends => ["d|e", "f|g"]),
    "d" => Dict(:conflicts => ["e"]),
    "e" => Dict(:conflicts => ["d"]),
    "f" => Dict(:conflicts => ["g"]),
    "g" => Dict(:conflicts => ["f"]),
    "y" => Dict(:depends => ["z"]),
    "z" => Dict{Any,Any}()
)

name2var = Dict{Any,Any}()
var2name = Dict{Any,Any}()

for (i,name) in enumerate(keys(index))
    name2var[name] = i
    var2name[i] = name
end

clauses = Any[]

function tocnf(name)
    # dependencies
    depends = get(index[name], :depends, Any[])
    for r in depends
        clause = Any[-(name2var[name])]
        append!(clause, Any[name2var[n] for n in split(r, '|')])
        push!(clauses, clause)
    end
    # conflicts
    conflicts = get(index[name], :conflicts, Any[])
    for c in conflicts
        push!(clauses, Any[-(name2var[name]), -(name2var[c])])
    end
end

foreach(tocnf, keys(index))

println(name2var)
push!(clauses, [name2var["a"]])

for sol in PicoSAT.itersolve(clauses)
    s = [var2name[v] for v in sol[sol .> 0]]
    println(sort(s))
end


text:
#=
BellmanFordTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BellmanFord.jl")
matrix =  GBMatrix([[0, Inf, Inf, 0.2, Inf, Inf, Inf] [0.3, 0, Inf, Inf, Inf, Inf, Inf] [Inf, Inf, 0, 0.4, Inf, 0.5, 0.1] [0.8, Inf, Inf, 0, Inf, Inf, 0.5] [Inf, 0.1, Inf, Inf, 0, Inf, 0.9] [Inf, Inf, 0.5, Inf, 0.1, 0, Inf] [Inf, 0.7, Inf, Inf, Inf, Inf, 0]])
@test bellmanford(matrix, 1, 7) == GBVector([0, 0.3, 1, 0.8, 0.4, 0.5, 1])


text:
# Autogenerated wrapper script for libcxxwrap_julia_jll for armv6l-linux-musleabihf-cxx11-julia_version+1.7.0
export libcxxwrap_julia, libcxxwrap_julia_stl

JLLWrappers.@generate_wrapper_header("libcxxwrap_julia")
JLLWrappers.@declare_library_product(libcxxwrap_julia, "libcxxwrap_julia.so.0")
JLLWrappers.@declare_library_product(libcxxwrap_julia_stl, "libcxxwrap_julia_stl.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcxxwrap_julia,
        "lib/libcxxwrap_julia.so",
        RTLD_GLOBAL,
    )

    JLLWrappers.@init_library_product(
        libcxxwrap_julia_stl,
        "lib/libcxxwrap_julia_stl.so",
        RTLD_GLOBAL,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function EnumeratedIntegerDistribution(arg0::Vector{jint})
    return EnumeratedIntegerDistribution((Vector{jint},), arg0)
end

function EnumeratedIntegerDistribution(arg0::Vector{jint}, arg1::Vector{jdouble})
    return EnumeratedIntegerDistribution((Vector{jint}, Vector{jdouble}), arg0, arg1)
end

function cumulative_probability(obj::EnumeratedIntegerDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_numerical_mean(obj::EnumeratedIntegerDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::EnumeratedIntegerDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_pmf(obj::EnumeratedIntegerDistribution)
    return jcall(obj, "getPmf", List, ())
end

function get_support_lower_bound(obj::EnumeratedIntegerDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::EnumeratedIntegerDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::EnumeratedIntegerDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function probability(obj::EnumeratedIntegerDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end


text:
include("getdigit_n.jl")
include("myprintf.jl")
include("simpson.jl")
using Printf
using .GetDigit_n
using .MyPrintf
using .Simpson

const DIGIT = 17
const MAXN = 10000

function atan(n)
    return BigFloat(4) * Simpson.simpson(x::BigFloat -> BigFloat(1) / (BigFloat(1) + x * x), BigFloat(0), BigFloat(1), n)
end

function main()
    n = GetDigit_n.getn_even(atan, MAXN, DIGIT, false)
    precision = GetDigit_n.getdigit(atan, n, DIGIT, false)
    precision = GetDigit_n.getdigit_2(atan, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, atan(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
using MixedSubdivisions
const MS = MixedSubdivisions
import MultivariatePolynomials
const MP = MultivariatePolynomials
import PolynomialTestSystems: equations, cyclic, ipp2, cyclooctane
using Test

@testset "MixedSubdivisions" begin
	@testset "Basic" begin
	    A₁ = [0 0 1 1; 0 2 0 1]
	    A₂ = [0 0 1 2; 0 1 1 0]

	    A = MS.cayley(A₁, A₂)

	    @test A == [0  0  1  1  0  0  1  2
	                0  2  0  1  0  1  1  0
	                1  1  1  1  0  0  0  0
	                0  0  0  0  1  1  1  1]

	    @test_throws ErrorException MS.cayley([1 0; 0 1], [1 0; 0 0; 0 2])

	    w₁ = [0, 0, 0, -2]
	    w₂ = [0, -3, -4, -8]
	    w = [w₁; w₂]

	    v₁ = [0, 0, 0, -1]
	    v₂ = copy(w₂)
	    v = [v₁; v₂]

	    mixed_cell_indices = [(2, 3), (1, 3)]
	    indexing = MS.CayleyIndexing(size.((A₁, A₂), 2))
	    ord = MS.DotOrdering(Int32.(w))
	    cell = MS.MixedCellTable(mixed_cell_indices, A, indexing)
	    @test cell.volume == 3
	    @test cell.circuit_table == [1 2; 3 0; 0 0; 1 -1; 0 3; 1 2; 0 0; -2 -1]
	    ineq = MS.first_violated_inequality(cell, v, ord)
	    @test ineq.config_index == 1
	    @test ineq.col_index == 4

	    @test MS.exchange_column(cell, MS.exchange_first, ineq) == MS.MixedCellTable([(4, 3), (1, 3)], A, indexing)
	    @test MS.exchange_column(cell, MS.exchange_second, ineq) == MS.MixedCellTable([(2, 4), (1, 3)], A, indexing)

	    ind_back = MS.reverse_index(ineq, cell, MS.exchange_second)
	    cell2 = MS.exchange_column(cell, MS.exchange_second, ineq)
	    @test cell2.volume == 2
	    @test cell == MS.exchange_column(cell2, MS.exchange_second, ind_back)

	    ind_back = MS.reverse_index(ineq, cell, MS.exchange_first)
	    cell2 = MS.exchange_column(cell, MS.exchange_first, ineq)
	    @test cell2.volume == 1
	    @test cell == MS.exchange_column(cell2, MS.exchange_first, ind_back)
	end

	@testset "Mixed Volume" begin
		@test mixed_volume(equations(cyclic(5))) == 70
		@test mixed_volume(equations(cyclic(7))) == 924
		@test mixed_volume(equations(cyclic(10))) == 35940
		@test mixed_volume(equations(cyclic(11))) == 184756
		@test mixed_volume(equations(ipp2())) == 288

		@test mixed_volume(equations(cyclic(5)), algorithm=:total_degree) == 70
		@test mixed_volume(equations(ipp2()), algorithm=:total_degree) == 288
	end

	@testset "Mixed Cells" begin
		A₁ = [0 0 1 1; 0 2 0 1]
		A₂ = [0 0 1 2; 0 1 1 0]

		A = [A₁, A₂]

		w₁ = [0, 0, 0, 2]
		w₂ = [0, 3, 4, 8]
		w = [w₁, w₂]

		v₁ = [0, 0, 0, 1]
		v₂ = copy(w₂)
		v = [v₁, v₂]

		cells_v = mixed_cells(A, v)
		@test length(cells_v) == 3
		@test sum(volume, cells_v) == 4
		@test sort(volume.(cells_v)) == [1, 1, 2]
		@test all(cells_v) do c
		    all(sort.(map(A_i -> A_i' * normal(c), A) .+ v)) do r
		        isapprox(r[1], r[2], atol=1e-12)
		    end
		end

		cells_w = mixed_cells(A, w)
		@test length(cells_w) == 2
		@test sum(volume, cells_w) == 4
		@test sort(volume.(cells_w)) == [1, 3]
		@test all(cells_w) do c
		    all(sort.(map(A_i -> A_i' * normal(c), A) .+ w)) do r
		        isapprox(r[1], r[2], atol=1e-12)
		    end
		end
	end

	@testset "Fine mixed cells" begin
		f = equations(cyclic(7))
		cells, lift = fine_mixed_cells(f)
		@test sum(c -> c.volume, cells) == 924
		@test lift isa Vector{Vector{Int32}}
	end

	@testset "Overflow error messages" begin
		f = equations(cyclooctane())
		@test_throws ArgumentError MS.mixed_volume(f)
		F = [f; randn(2, 18) * [MP.variables(f);1]]
		A = support(F)
		lifting = map(Ai -> MS.gaussian_lifting_sampler(size(Ai,2)), A)
		@test_throws OverflowError mixed_cells(A, lifting)
		@test fine_mixed_cells(F) === nothing
	end
end


text:
## Import packages
using StatsBase
using Plots

## Define Functions
function gillespie(transitions, initial, N, select_hyp="none", t0=0.0)

    """
    Implementation of the Gillespie algorithm for Lotka-Volterra model

    Args:
    transitions(dict): A dictionary with keys "prey" and "pred", and each element
        is a list containing the row of the transition matrix.
    initial(dict): A dictionary with keys 'prey' and 'predator' giving the initial
        counts for each species at the start of the sim.
    N(int): Number of desired simulation iterations.
    select_hyp(str): Default is "none", wherein a uniform prior is assigned to the
        three hypotheses.  Otherwise, "H1", "H2", or "H3" concentrates all prior
        probability on one hypothesis (useful in forward simulation).
    t0(float): initial time point.
    """

    trajectory = Dict("prey" => [initial["prey"]], "pred" => [initial["pred"]])
    times = [t0]
    ecology = initial

    theta_hypotheses = Dict(
        "H1" => Dict("spawn_prey" => 1, "prey2pred" => 0.005, "pred_dies" => 0.6),
        "H2" => Dict("spawn_prey" => .9, "prey2pred" => 0.004, "pred_dies" => 0.4),
        "H3" => Dict("spawn_prey" => 1.5, "prey2pred" => 0.01, "pred_dies" => 0.8))

    if select_hyp != "none"
        theta = theta_hypotheses[select_hyp]
    else
        hypoth_list = ["H1", "H2", "H3"]
        weights = [1 - 0.6666, 0.3333, 0.3333]
        theta = theta_hypotheses[sample(hypoth_list, Weights(weights))]
    end

    t = times[1]

    for i = 1:N
        hazards = get_hazards(ecology, theta)
        transition = transitions[sample(collect(keys(hazards)), Weights(collect(values(hazards))))]
        t = t + exp(sum(values(hazards)))

        ecology["prey"] += transition[1]
        ecology["pred"] += transition[2]
        # Enforce only positive integers
        ecology["prey"] = max(0, ecology["prey"])
        ecology["pred"] = max(0, ecology["pred"])

        append!(trajectory["prey"], ecology["prey"])
        append!(trajectory["pred"], ecology["pred"])
        append!(times, t)
    end

    return Dict(
        "prey" => trajectory["prey"],
        "pred" => trajectory["pred"],
        "times" => times
        )
end

function get_hazards(ecology, theta)
    """
    Compute the hazard function given the current states.  Note that in this
        model the function depends only on state, not time, though this is not
        the case in general. "spawn_prey" represents the event of a prey being born,
        "prey2pred" represents a predator consuming a new prey and consequently spawning
        a new predator, "pred_dies" represents the death of a predator.

    args:
        ecology(dict): a dictionary containing species counts.
        theta(dict): A dictionary where keys represent events (reactions) and
            values are the location hyperparameters for rate parameter constants
            corresponding to those events.
    """

    return Dict(
        "spawn_prey" => theta["spawn_prey"] * ecology["prey"],
        "prey2pred" => theta["prey2pred"] * ecology["prey"] * ecology["pred"],
        "pred_dies" => theta["pred_dies"] * ecology["pred"]
        )
end

## Initialize Function
initial = Dict("prey" => 50, "pred" => 100)
Pre = [[1, 0], [1, 1], [0, 1]]
Post = [[2, 0], [0, 2], [0, 0]]
transition_mat = Post - Pre
transitions = Dict("spawn_prey" => transition_mat[1,],
                    "prey2pred" => transition_mat[2,],
                    "pred_dies" => transition_mat[3,])

gillespie_sim = gillespie(transitions, initial, 10000, "H2", 0.0)

plot_vals = hcat(values(gillespie_sim["prey"]), values(gillespie_sim["pred"]))
plot(1:10001, plot_vals,
title = "Non-Omega Simulation",
xlabel = "Time",
ylabel = "Quantity",
label = ["Prey" "Predators"],
lw = 1.25)