text:
#function merge{T<:Any}(t1::Vector{Timestamp{T}}, t2::Vector{Timestamp{T}}; sorted=true)
#function merge(t1::Vector{Timestamp}, t2::Vector{Timestamp}; sorted=true)
function merge(t1::TimestampArray, t2::TimestampArray; sorted=true)
    ts = vcat(t1, t2)
    if sorted
        sdates = sort(Date[t.timestamp for t in ts]) # get dates and sort
        udates = unique(sdates)
        res    = ts[udates[1]]  # initialize the array to be pushed
        for i in 2:length(udates)
            res = vcat(res, ts[udates[i]])
        end
        return res
    else
        return ts
    end
end


text:
function btest(i::Int,pos::Int)
  return  (((i)>>(pos)) & 1)  == 1
end

function hilbert3d(x::Int,y::Int,z::Int,bit_length::Int,npoint::Int)

state_diagram = [  1, 2, 3, 2, 4, 5, 3, 5,
                   0, 1, 3, 2, 7, 6, 4, 5,
                   2, 6, 0, 7, 8, 8, 0, 7,
                   0, 7, 1, 6, 3, 4, 2, 5,
                   0, 9,10, 9, 1, 1,11,11,
                   0, 3, 7, 4, 1, 2, 6, 5,
                   6, 0, 6,11, 9, 0, 9, 8,
                   2, 3, 1, 0, 5, 4, 6, 7,
                   11,11, 0, 7, 5, 9, 0, 7,
                   4, 3, 5, 2, 7, 0, 6, 1,
                   4, 4, 8, 8, 0, 6,10, 6,
                   6, 5, 1, 2, 7, 4, 0, 3,
                   5, 7, 5, 3, 1, 1,11,11,
                   4, 7, 3, 0, 5, 6, 2, 1,
                   6, 1, 6,10, 9, 4, 9,10,
                   6, 7, 5, 4, 1, 0, 2, 3,
                   10, 3, 1, 1,10, 3, 5, 9,
                   2, 5, 3, 4, 1, 6, 0, 7,
                   4, 4, 8, 8, 2, 7, 2, 3,
                   2, 1, 5, 6, 3, 0, 4, 7,
                   7, 2,11, 2, 7, 5, 8, 5,
                   4, 5, 7, 6, 3, 2, 0, 1,
                   10, 3, 2, 6,10, 3, 4, 4,
                   6, 1, 7, 0, 5, 2, 4, 3 ]
state_diagram = reshape(state_diagram, 8 ,2, 12)

x_bit_mask = falses(bit_length)
y_bit_mask = falses(bit_length)
z_bit_mask = falses(bit_length)

i_bit_mask = falses(3*bit_length)

order=0.     #zeros(npoint); npoint=1, therefore no array

#for ip=1:npoint
  # convert to binary
  for i=1:bit_length
    x_bit_mask[i]=btest(x,i-1)
    y_bit_mask[i]=btest(y,i-1)
    z_bit_mask[i]=btest(z,i-1)
  end

  # interleave bits
  for i=0:(bit_length-1)
    i_bit_mask[3*i+2+1]=x_bit_mask[i+1]
    i_bit_mask[3*i+1+1]=y_bit_mask[i+1]
    i_bit_mask[3*i+1]=z_bit_mask[i+1]

  end

  # build Hilbert ordering using state diagram
  cstate=0
  staterange = range(bit_length-1, stop=0, step=-1)
  for i in staterange
    #println(i)

    if i_bit_mask[3*i+2+1] b2=1 else b2=0 end
    if i_bit_mask[3*i+1+1] b1=1 else b1=0 end
    if i_bit_mask[3*i+1]   b0=1 else b0=0 end

    sdigit=b2*4+b1*2+b0
    nstate=state_diagram[sdigit+1,1,cstate+1]
    hdigit=state_diagram[sdigit+1,2,cstate+1]
    #println("nstate: ", nstate)
    #println("hdigit: ", hdigit)
    i_bit_mask[3*i+2+1]=btest(hdigit,2)
    i_bit_mask[3*i+1+1]=btest(hdigit,1)
    i_bit_mask[3*i+1]  =btest(hdigit,0)
    #println("i_bit_mask: ", i_bit_mask[3*i+2+1])
    #println("i_bit_mask: ", i_bit_mask[3*i+1+1])
    #println("i_bit_mask: ", i_bit_mask[3*i+1])
    cstate=nstate
  end

  # save Hilbert key as double precision real
  #order[ip]=0.
  order=0.     #zeros(npoint); npoint=1, therefore no array
  for i=1:(3*bit_length)
    if i_bit_mask[i] b0=1 else b0=0 end
    #order[ip]=order[ip]+b0*2^i
    order=order+b0*2^(i-1)
    #println(i-1, " ", order[ip])
  end

#end
    return order
end

ext: jl


text:
@testset "promote 1" begin
  @test typeof(promote(one(Double64), one(Double32))[1]) === Double64
  @test typeof(promote(one(Double64), one(Double16))[1]) === Double64
  @test typeof(promote(one(Double32), one(Double16))[1]) === Double32

  @test typeof(promote(one(Double64), one(BigInt))[1]) === BigInt
  @test typeof(promote(one(Double32), one(BigInt))[1]) === BigInt
  @test typeof(promote(one(Double16), one(BigInt))[1]) === BigInt

  @test typeof(promote(one(Double64), one(BigFloat))[1]) === BigFloat
  @test typeof(promote(one(Double32), one(BigFloat))[1]) === BigFloat
  @test typeof(promote(one(Double16), one(BigFloat))[1]) === BigFloat

  @test promote_type(Double64, Double32) === Double64
  @test promote_type(Double64, Double16) === Double64
  @test promote_type(Double32, Double16) === Double32
end

@testset "promote 2" begin

 f2 = 2.0

  @test promote_type(DoubleFloat{Float64}, DoubleFloat{Float32}) == DoubleFloat{Float64}
  @test promote_type(DoubleFloat{Float64}, DoubleFloat{Float16}) == DoubleFloat{Float64}
  @test promote_type(DoubleFloat{Float32}, DoubleFloat{Float16}) == DoubleFloat{Float32}

  @test promote_type(DoubleFloat{Float64}, BigFloat) == BigFloat
  @test promote_type(DoubleFloat{Float32}, BigFloat) == BigFloat
  @test promote_type(DoubleFloat{Float16}, BigFloat) == BigFloat

  @test promote_type(DoubleFloat{Float64}, BigInt) == BigInt
  @test promote_type(DoubleFloat{Float32}, BigInt) == BigInt
  @test promote_type(DoubleFloat{Float16}, BigInt) == BigInt

  @test promote(DoubleFloat{Float64}(f2), DoubleFloat{Float32}(f2))[2] == DoubleFloat{Float64}(f2)
  @test promote(DoubleFloat{Float64}(f2), DoubleFloat{Float16}(f2))[2] == DoubleFloat{Float64}(f2)
  @test promote(DoubleFloat{Float32}(f2), DoubleFloat{Float16}(f2))[2] == DoubleFloat{Float32}(f2)

  @test promote(DoubleFloat{Float64}(f2), BigFloat(f2))[2] == DoubleFloat{Float64}(f2)
  @test promote(DoubleFloat{Float32}(f2), BigFloat(f2))[2] == DoubleFloat{Float32}(f2)
  @test promote(DoubleFloat{Float16}(f2), BigFloat(f2))[2] == DoubleFloat{Float16}(f2)

  @test promote(DoubleFloat{Float64}(f2), BigInt(f2))[2] == DoubleFloat{Float64}(f2)
  @test promote(DoubleFloat{Float32}(f2), BigInt(f2))[2] == DoubleFloat{Float32}(f2)
  @test promote(DoubleFloat{Float16}(f2), BigInt(f2))[2] == DoubleFloat{Float16}(f2)

end

@testset "promote_rule" begin

  @test promote_rule(DoubleFloat{Float64}, DoubleFloat{Float32}) == DoubleFloat{Float64}
  @test promote_rule(DoubleFloat{Float64}, DoubleFloat{Float16}) == DoubleFloat{Float64}
  @test promote_rule(DoubleFloat{Float32}, DoubleFloat{Float16}) == DoubleFloat{Float32}
  @test promote_rule(DoubleFloat{Float64}, Float64) == DoubleFloat{Float64}
  @test promote_rule(DoubleFloat{Float32}, Float32) == DoubleFloat{Float32}
  @test promote_rule(DoubleFloat{Float16}, Float16) == DoubleFloat{Float16}
  @test promote_rule(DoubleFloat{Float64}, Int64) == DoubleFloat{Float64}
  @test promote_rule(DoubleFloat{Float32}, Int32) == DoubleFloat{Float32}
  @test promote_rule(DoubleFloat{Float16}, Int16) == DoubleFloat{Float16}
  @test promote_rule(DoubleFloat{Float64}, BigInt) == BigInt
  @test promote_rule(DoubleFloat{Float32}, BigInt) == BigInt
  @test promote_rule(DoubleFloat{Float16}, BigInt) == BigInt
  @test promote_rule(BigInt, DoubleFloat{Float64}) == BigInt
  @test promote_rule(BigInt, DoubleFloat{Float32}) == BigInt
  @test promote_rule(BigInt, DoubleFloat{Float16}) == BigInt
  @test promote_rule(DoubleFloat{Float64}, BigFloat) == BigFloat
  @test promote_rule(DoubleFloat{Float32}, BigFloat) == BigFloat
  @test promote_rule(DoubleFloat{Float16}, BigFloat) == BigFloat
  @test promote_rule(BigFloat, DoubleFloat{Float64}) == BigFloat
  @test promote_rule(BigFloat, DoubleFloat{Float32}) == BigFloat
  @test promote_rule(BigFloat, DoubleFloat{Float16}) == BigFloat
  @test promote_rule(AbstractIrrational, DoubleFloat{Float16}) == DoubleFloat{Float16}
  @test promote_rule(AbstractIrrational, DoubleFloat{Float32}) == DoubleFloat{Float32}
  @test promote_rule(AbstractIrrational, DoubleFloat{Float64}) == DoubleFloat{Float64}
  @test widen(DoubleFloat{Float16}) == DoubleFloat{Float32}
  @test widen(DoubleFloat{Float32}) == DoubleFloat{Float64}
  @test widen(DoubleFloat{Float64}) == BigFloat

end


text:
using LinearAlgebra: sqrt
using Printf: Threads
using Plots: Threads
using ProgressMeter
using Printf
include("PolygonBase.jl")
include("ESBase.jl")
include("silhouette&pieces.jl")

# X[1つめのピースのx座標, 1つめのピースのy座標, １つめの回転角Θ, 2つめのx座標, ...]
function loss_poly(X::Array, loss_args::Array{Float64, 1})
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    tmp = LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))
    A = LibGEOS.area(tmp) / LibGEOS.area(MYPolygon2LibGEOS(silhouette))

    # 各2pieceの共通部分を求める
    inter = 1.0 - LibGEOS.area(unioned) / sum_pieces
    outer = 1.0 - LibGEOS.area(tmp)/sum_pieces

    # 頂点の誤差成分 (シルエットベース 緩い)
    Δv = 0.0
    for v_s in 1:silhouette.n
        Δmin = 100.0
        for j in 1:length(pieces)
            tmp_p = move(pieces[j], X[3j - 2], X[3j - 1])
            rotate!(tmp_p, X[3j] * 2π)
            for v_p in 1:tmp_p.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - tmp_p.vertexes[:, v_p]).^2))
            end
        end
        Δv += Δmin / silhouette.n
        # Δv = max(Δv, Δmin)
    end

    # 頂点の誤差成分 (ピースベース 厳しい)
    """
    Δv = 0.0
    for i in 1:length(pieces)
        p = move(pieces[i], X[3i - 2], X[3i - 1])
        rotate!(p, X[3i] * 2π)
        Δmax = 0.0
        # 頂点ごとに誤差を見る
        for v_p in 1:p.n
            Δmin = 100.0
            # シルエットについて見る
            for v_s in 1:silhouette.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - p.vertexes[:, v_p]).^2))
            end
            # その他のピースについて見る
            for j in 1:length(pieces)
                if j == i continue end
                q = move(pieces[j], X[3j - 2], X[3j - 1])
                rotate!(p, X[3j] * 2π)
                for v_q in 1:q.n
                    Δmin = min(Δmin, sum((q.vertexes[:,v_q] - p.vertexes[:, v_p]).^2))
                end
            end
            Δmax = max(Δmax, Δmin)
        end
        # Δv += 2tanh(2Δmin) / (π*silhouette.n)
        # Δv = max(Δv, Δmin)
        Δv += Δmax
    end
    """

    # 辺の誤差成分
    E = 0.0
    for i in 1:length(pieces)
        P = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(P, X[3i] * 2π)

        cosθ = 0.0

        # シルエットの各頂点と距離が最小な点を探す
        min_d = 100.0
        for p in 0:P.n-1, s in 0:silhouette.n-1
            d = sum((P.vertexes[:, p + 1] - silhouette.vertexes[:, s + 1]).^2)
            # 最小の大きさの角度を探す
            if d < min_d
                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ1 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ2 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ3 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ4 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                cosθ = max(cosθ1, cosθ2, cosθ3, cosθ4)
            end
        end
        E += cosθ / length(pieces)
    end

    return -loss_args[1] * A + loss_args[2] * Δv + loss_args[3] * outer + loss_args[4] * (E - 1)
end

# X[1つめのピースのx座標, 1つめのピースのy座標, １つめの回転角Θ, 2つめのx座標, ...]
function loss_poly_light(X::Array)
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    tmp = LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))
    A = LibGEOS.area(tmp) / LibGEOS.area(MYPolygon2LibGEOS(silhouette))

    return -A * 100
end

function loss_poly2(X::Array, loss_args::Array{Float64, 1})
    global pices, silhouette

    all_area = 0
    for P in pieces
        all_area += P.area
    end

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    F = LibGEOS.area(LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))) / silhouette.area

    # 辺の誤差成分
    E = 0.0
    for P in pieces
        sinθ_ave = 1
        for i in 1:P.n
            minimum_e = 1
            p1 = P.vertexes[:, i]
            p2 = i < P.n ? P.vertexes[:, i+1] : P.vertexes[:, 1]

            p21 = p2 - p1
            norm_p21 = norm(p21)

            # シルエットを見る
            for j in 1:silhouette.n
                s1 = silhouette.vertexes[:, j]
                s2 = j < silhouette.n ? silhouette.vertexes[:, j+1] : silhouette.vertexes[:, 1]
                sinθ = cross(p21, s2-s1)/(norm_p21*norm(s2-s1))
                minimum_e = min(sinθ, minimum_e)
            end
            # 他のピースを見る
            for Q in pieces
                if P == Q
                    continue
                end
                for j in 1:Q.n
                    q1 = Q.vertexes[:, j]
                    q2 = j < Q.n ? Q.vertexes[:, j+1] : Q.vertexes[:, 1]
                    sinθ = cross(p21, q2-q1)/(norm_p21*norm(q2-q1))
                    minimum_e = min(abs(sinθ), minimum_e)
                end
            end
           sinθ_ave += minimum_e / P.n
        end
        E += P.area *sinθ_ave / (length(pieces) * all_area)
    end

    # 頂点の誤差成分
    V = 0.0
    for P in pieces
        v = 0.0
        for i in 1:P.n
            p = P.vertexes[:, i]

            min_d = 100.0
            # シルエットについて見る
            for j in 1:silhouette.n
                ps = norm(p - silhouette.vertexes[:, j])
                min_d = min(min_d, ps)
            end

            # 他のピースについて見る
            for Q in pieces
                if Q == P
                    continue
                end
                for j in 1:Q.n
                    pq = norm(p - Q.vertexes[:, j])
                    min_d = min(min_d, pq)
                end
            end
            v += min_d / P.n
        end
        V += P.area * v / (length(pieces) * all_area)
    end

    return -loss_args[1] * F + loss_args[2] * E + loss_args[3] * V
end

silhouette = hexagon_m
pieces = [tri_s, tri_s, parallelogram, tri_m]

# 初期化
max_gen = 128
sample_num = 128

# c2 0~200, c3 0~200
for c2 in 40.0:20:100.0, c3 in 350.0:50:550.0
    loss_args = [100.0, c3, 0.0, c2]

    fpath = @sprintf "outputs\\outputs_c2-%03d_c3-%03d" Int(c2) Int(c3)
    rm(fpath, force=true, recursive=true)
    mkdir(fpath )

    best_fitnesses = Array{Float64, 1}()
    p = Progress(sample_num, desc="c2:$(Int(c2)) c3:$(Int(c3)): ")
    for t in 1:sample_num
        # local cmaes, rng
        cmaes = init_CMAES(zeros(3 * length(pieces)), 0.8, 18)
        fitnesses = zeros(cmaes.dim)
        X = zeros((cmaes.dim, cmaes.λ))
        rng = MersenneTwister(t)

        for gen in 1:max_gen
            # 個体生成
            X = samplePopulation(cmaes, rng=rng)

            # 回転角のパラメータを -π ~ π までに正規化
            for j in 1:size(X)[2]
                for i in 3:3:size(X)[1]
                    X[i, j] = mod(X[i, j], 1.0)
                end
            end

            # 評価値
            loss(x) = loss_poly(x, loss_args)
            fitnesses = get_fitness(X, loss)

            # 更新
            update!(cmaes, X, fitnesses, gen)
        end
        push!(best_fitnesses, maximum(fitnesses))

        # 表示用
        best_arg = argmin(fitnesses)
        best_pieces = Array{MYPolygon, 1}()
        for i in 1:length(pieces)
            todisplay = move(pieces[i], X[3i-2, best_arg], X[3i-1, best_arg])
            rotate!(todisplay, X[3i, best_arg] * 2π)
            push!(best_pieces, todisplay)
        end
        display(silhouette, best_pieces...)
        savefig(@sprintf "outputs\\outputs_c2-%03d_c3-%03d\\%02.0f_%d_results" Int(c2) Int(c3) -best_fitnesses[end] t)

        next!(p)
    end
end

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
@testset "Devices" begin
  include("Amplifier/AmplifierTests.jl")
  include("DAQ/DAQTests.jl")
  include("GaussMeter/GaussMeterTests.jl")
  include("Virtual/VirtualDeviceTests.jl")
  include("Robots/RobotTests.jl")
end


text:
# a convenient way of interacting with a simulation
# maintained by @zsunberg

"""
    sim(polfunc::Function, mdp::MDP[, initial_state]; [kwargs...])
    sim(polfunc::Function, pomdp::POMDP[, initial_state])

Alternative way of running a simulation with a function specifying how to calculate the action at each timestep.

# Usage

    sim(mdp) do s
        # code that calculates action `a` based on `s` - this is the policy
        # you can also do other things like display something
        return a
    end

for an MDP or

    sim(pomdp) do o
        # code that does belief updates with observation `o` and calculates `a`
        # you can also do other things like display something
        return a
    end

for a POMDP.

# Keyword Arguments

Use the `simulator` keyword argument to specify any simulator to run the simulation. If nothing is specified for the simulator, a HistoryRecorder will be used as the simulator, with all keyword arguments forwarded to it, e.g.

    sim(mdp, max_steps=100, show_progress=true) do s
        # ...
    end

will limit the simulation to 100 steps.

The POMDP version also has two additional keyword arguments:
- `initialobs`: this will control the initial observation given to the policy function.
- `updater`: if provided, this updater will be used to update the belief, and the belief will be used as the argument to the policy function. If a custom updater is provided, the `initialobs` keyword argument should be used to specify the initial belief.
"""
function sim end

function sim(polfunc::Function, mdp::MDP,
             initialstate=nothing;
             simulator=nothing,
             kwargs...
            )

    kwargd = Dict(kwargs)
    if initialstate==nothing && statetype(mdp) != Nothing
        if haskey(kwargd, :initialstate)
            initialstate = pop!(kwargd, :initialstate)
        else
            initialstate = default_init_state(mdp)
        end
    end
    delete!(kwargd, :initialstate)
    if simulator==nothing
        simulator = HistoryRecorder(;kwargd...)
    end
    policy = FunctionPolicy(polfunc)
    simulate(simulator, mdp, policy, initialstate)
end

function sim(polfunc::Function, pomdp::POMDP,
             initialstate=nothing;
             simulator=nothing,
             initialobs=nothing,
             updater=nothing,
             kwargs...
            )

    kwargd = Dict(kwargs)
    if initialstate==nothing && statetype(pomdp) != Nothing
        if haskey(kwargd, :initialstate)
            initialstate = pop!(kwargd, :initialstate)
        else
            initialstate = default_init_state(pomdp)
        end
    end
    delete!(kwargd, :initialstate)
    if simulator==nothing
        simulator = HistoryRecorder(;kwargd...)
    end
    if initialobs==nothing && obstype(pomdp) != Nothing
        initialobs = default_init_obs(pomdp, initialstate)
    end
    if updater==nothing
        updater = PreviousObservationUpdater()
    end
    policy = FunctionPolicy(polfunc)
    simulate(simulator, pomdp, policy, updater, initialobs, initialstate)
end

function default_init_obs(p::POMDP, s)
    if implemented(generate_o, Tuple{typeof(p), typeof(s), typeof(Random.GLOBAL_RNG)})
        return generate_o(p, s, Random.GLOBAL_RNG)
    else
        return missing
    end
end

@generated function default_init_state(p::Union{MDP,POMDP})
    if implemented(initialstate, Tuple{p, typeof(Random.GLOBAL_RNG)})
        return :(initialstate(p, Random.GLOBAL_RNG))
    else
        return quote
            error("""
                  Error in sim(::$(typeof(p))): No initial state specified.

                  Please supply it as an argument after the mdp or define the method POMDPs.initialstate(::$(typeof(p)), ::$(typeof(Random.GLOBAL_RNG))) or define the method POMDPs.initialstate_distribution(::$(typeof(p))).

                  """)
        end
    end
end


text:
@testset "minimal_intersection_test" begin
    @test intersect(interval(1.0,3.0), interval(2.1,4.0)) == interval(2.1,3.0)
    @test intersect(interval(1.0,3.0), interval(3.0,4.0)) == interval(3.0,3.0)
    @test intersect(interval(1.0,3.0), emptyinterval()) == emptyinterval()
    @test intersect(entireinterval(), emptyinterval()) == emptyinterval()
    @test intersect(interval(1.0,3.0), entireinterval()) == interval(1.0,3.0)
end

@testset "minimal_convex_hull_test" begin
    @test hull(interval(1.0,3.0), interval(2.1,4.0)) == interval(1.0,4.0)
    @test hull(interval(1.0,1.0), interval(2.1,4.0)) == interval(1.0,4.0)
    @test hull(interval(1.0,3.0), emptyinterval()) == interval(1.0,3.0)
    @test hull(emptyinterval(), emptyinterval()) == emptyinterval()
    @test hull(interval(1.0,3.0), entireinterval()) == entireinterval()
end


text:
using FastBEAST

function test_fullrank(N)
    store = Vector{MatrixView}(undef, N)
    n = 100
    for i=1:N
        store[i] = FullMatrixView(rand(n,n), [1; 2], [1; 2], 10, 20)
    end

    return store
end

function test_fullrank2(N)
    store = MatrixView[]
    n = 100
    for i=1:N
        push!(store, FullMatrixView(rand(n,n), [1; 2], [1; 2], 10, 20))
    end

    return store
end

function test_fullrank3(N)
    store = Vector{FullMatrixView{Float32,Int64}}(undef, N)
    n = 100
    for i=1:N
        store[i] = FullMatrixView(rand(n,n), [1; 2], [1; 2], 10, 20)
    end

    return store
end

function test_lowrank(N)
    store = Vector{MatrixView}(undef, N)
    n = 100
    for i=1:N
        store[i] = LowRankMatrixView(rand(n,n), rand(n,n), [1; 2], [1; 2], 10, 20)
    end

    return store
end

function test_lowrank2(N)
    store = MatrixView[]
    n = 100
    for i=1:N
        push!(store, LowRankMatrixView(rand(n,n), rand(n,n), [1; 2], [1; 2], 10, 20))
    end

    return store
end

function test_lowrank3(N)
    store = Vector{LowRankMatrixView}(undef, N)
    n = 100
    for i=1:N
        store[i] = LowRankMatrixView(rand(n,n), rand(n,n), [1; 2], [1; 2], 10, 20)
    end

    return store
end

function test_all(N)
@time    A = test_fullrank(N)
@time    A = test_fullrank2(N)
@time    A = test_fullrank3(N)
@time    A = test_lowrank(N)
@time    A = test_lowrank2(N)
@time    A = test_lowrank3(N)
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
# missing methods (in Gtk)
function define_missing_gtk_methods()
    @eval ComplexVisualGtk begin
        close_window(win::Gtk.GtkWindow) = ccall(
            (:gtk_window_close, Gtk.libgtk),
            Nothing, (Ptr{Gtk.GObject},), win)
    end
end
define_missing_gtk_methods()

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
# Band-pass wavelet
# Author: Gabrio Rizzuti; Ziyi Yin, ziyi.yin@gatech.edu
# Date: August 2020

export cfreqfilt, cfreq_wavelet

###############################################################################
# Generate band-pass source wavelet

function cfreqfilt(nt::Int64, dt::Float32, cfreqs::NTuple{4, Float32}; type::String = "bandpass", edge = hamming, flag_fft::Bool = true)
	# Initialization
	Filter = zeros(ComplexF32, nt)
	# Index corresponding to corner frequencies
	i1 = Int64(ceil(cfreqs[1]*nt*dt)+1)
	i2 = Int64(ceil(cfreqs[2]*nt*dt)+1)
	i3 = Int64(floor(cfreqs[3]*nt*dt)+1)
	i4 = Int64(floor(cfreqs[4]*nt*dt)+1)
	if type == "lowpass"
		i1 = 1
		i2 = 1
	elseif type == "highpass"
		i3 = nt
		i4 = nt
	end
	# Populating filter in the frequency domain
	Filter[1:i1-1] .= 0
	k = collect(i1:i2-1); Filter[k] .= (k.-i1)./(i2-i1)
	Filter[i2:i3] .= 1
	k = collect(i3+1:i4); Filter[k] .= Float32(1).-(k.-i3)./(i4-i3)
	Filter[i4+1:end] .= 0
	SL = Int64((nt-1)/2)
	Filter[SL+2:end] = Filter[SL+1:-1:2]
	# Zeroing edge values in time domain
	filter = ifftshift(Float32.(edge(length(Filter))).*fftshift(real.(ifft(Filter))))
	# Final Fourier transform
	if flag_fft
		return fft(filter)
	else
		return filter
	end
end

function cfreq_wavelet(SL::Int64, nt::Int64, dt::Float32, cfreqs::NTuple{4, Float32}; type::String = "bandpass", edge = hamming, flag_fft::Bool = true)
	return [ifftshift(cfreqfilt(2*SL+1, dt, cfreqs; type = type, edge = edge, flag_fft = false)); zeros(Float32, nt-(2*SL+1))]
end


text:
@testset "GPS L1" begin

    gpsl1 = GPSL1()
    @test @inferred(get_center_frequency(gpsl1)) == 1.57542e9Hz
    @test @inferred(get_code_length(gpsl1)) == 1023
    @test @inferred(get_secondary_code_length(gpsl1)) == 1
    @test @inferred(get_code(gpsl1, 0, 1)) == 1
    @test @inferred(get_code(gpsl1, 0.0, 1)) == 1
    @test @inferred(GNSSSignals.get_code_unsafe(gpsl1, 0.0, 1)) == 1
    @test @inferred(get_data_frequency(gpsl1)) == 50Hz
    @test @inferred(get_code_frequency(gpsl1)) == 1023e3Hz
    @test get_code.(gpsl1, 0:1022, 1) == L1_SAT1_CODE

end


text:
"""
    OCEstimator <: Estimator

Compatible estimator minimizing expected MSE for response rate `p` see also:

> Kunzmann K, Kieser M. Point estimation and p‐values in phase II adaptive two‐stage designs with a binary endpoint. Statistics in medicine. 2017 Mar 15;36(6):971-84.
"""
struct OCEstimator{TD<:Design} <: Estimator

  design::TD
  estimates::DataFrames.DataFrame

  function OCEstimator{TD}(
    design::TD,
    solver::TS,
    prior::Function,
    k::Integer
  ) where {TD<:Design,TS<:MathProgBase.AbstractMathProgSolver}

    epsilon = 10.0^(-6.0)
    n1   = interimsamplesize(design)
    nmax = maximum(samplesize(design))
    supp = support(design)
    model = JuMP.Model(solver = solver)
    JuMP.@variable(model, 0.0 <= est[x1 = 0:n1, x2 = 0:(nmax - n1); ispossible(design, x1, x2)] <= 1.0)
    h = 1 / k
    vrho = collect(linspace(0, 1, k + 1))

    function coefs(x1, x2, rho)

      if !ispossible(design, x1, x2)
        return 0.0
      else
        pdf = prior(rho)
        pdf = (pdf > 10.0) | !isfinite(pdf) ? 10 : pdf
        return pdf * rho^(x1 + x2) * (1.0 - rho)^(samplesize(design, x1) - x1 - x2)
      end

    end # coefs

    JuMP.@objective(model, Min,
      sum(binomial(BigInt(n1), BigInt(x1))*binomial(BigInt(samplesize(design, x1) - n1), BigInt(x2))*h/3*(
              coefs(x1, x2, vrho[2*i - 1])*(est[x1, x2] - vrho[2*i - 1])^2 +
            4*coefs(x1, x2,     vrho[2*i])*(est[x1, x2] -     vrho[2*i])^2 +
              coefs(x1, x2, vrho[2*i + 1])*(est[x1, x2] - vrho[2*i + 1])^2) for
        x1 in 0:n1, x2 in 0:(nmax - n1), i in 1:convert(Int64, floor(k/2)) if
        ispossible(design, x1, x2)
      )
    )
    JuMP.@variable(model, 0 <= minREJECT <= 1)
    JuMP.@variable(model, 0 <= maxACCEPT <= 1)
    for i in 1:size(supp, 1)
      x1, x2 = supp[i, :]
      if test(design, x1, x2) == true
        JuMP.@constraint(model, minREJECT - est[x1, x2] <= 0)
      else
        JuMP.@constraint(model, maxACCEPT - est[x1, x2] >= 0)
      end
    end
    JuMP.@constraint(model, minREJECT - maxACCEPT - epsilon >= 0) # smallest estimate of reject must be larger than largest estimate of accept
    for i in size(supp, 1)
      for j in size(supp, 1)
        if i < j
          x1a, x2a = supp[i, :]
          x1b, x2b = supp[j, :]
          if (x1a == x1b) & (x2a >= x2b) # same first stage, larger or equal second stage ~~> larger or equal estimate
            JuMP.@constraint(model, est[x1a, x2a] - est[x1b, x2b] >= 0)
          end
          n2a = samplesize(design, x1a)
          n2a = samplesize(design, x1a)
          if (n2a == interimsamplesize(design)) & (n2b == interimsamplesize(design)) & (x1a >= x1b) # no second stage, larger or equal first stage ~~> larger or equal estimate
            JuMP.@constraint(model, est[x1a, x2a] - est[x1b, x2b] >= 0)
          end
        end
      end
    end
    JuMP.solve(model)
    estimates = Array{Float64}(n1 + 1, nmax - n1 + 1)
    for x1 in 0:n1
      for x2 in 0:(nmax - n1) # brute force, not all possible
        if ispossible(design, x1, x2)
          estimates[x1 + 1, x2 + 1] = JuMP.getvalue(est[x1, x2])
        else
          estimates[x1 + 1, x2 + 1] = convert(Float64, NaN)
        end
      end
    end
    new(design, estimates)

  end # inner constructor

end

"""
    OCEstimator{TS<:MathProgBase.AbstractMathProgSolver}(
      design::Design,
      solver::TS;
      prior::Function = jeffreysprior(design),
      k = 100
    )

Create compatible estimator minimizing expected MSE for response rate `p` see also:

> Kunzmann K, Kieser M. Point estimation and p‐values in phase II adaptive two‐stage designs with a binary endpoint. Statistics in medicine. 2017 Mar 15;36(6):971-84.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| design       | Design |
| solver       | MathProgBase solver, must support quadratic expressions |
| prior        | weight function for MSE values at different `p`, must be of form `f(p::Real)::Real` |
| k            | number of equally spaced grid-points for evaluation of MSE and prior |
"""
function OCEstimator(
  design::TD, solver::TS;
  prior::Function = jeffreysprior(design), k::Integer = 100
) where {TD<:Design,TS<:MathProgBase.AbstractMathProgSolver}

  isodd(k) ?
    error("k must be even") : nothing
  return OCEstimator{TD}(design, solver, prior, k)

end # outer constructor

Base.show(io::IO, estimator::OCEstimator) = print("OCEstimator")

function estimate(estimator::OCEstimator, x1::T, x2::T) where {T<:Integer}

  !ispossible(design(estimator), x1, x2) ?
      error("x1/x2 combination not possible under design") : nothing
    return estimator.estimates[x1 + 1, x2 + 1]

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
#=
    Os can be either
        - a nonzero scalar
        - a tensormap
        - zero

    if it is a scalar; implicitly we have the identity operator at that site - so we can optimize the contraction
    if it is zero then there is nothing to optimize, the result is zero anyway
    if it's a tensormap, then we do the contraction

    In a later stage we can perhaps support nonzero scalar -> scalar*isometry - though the usecases for that seem very limited

    The principles are that you can write code without knowing anything about the senvse structure, and it should just work (potentially a bit slower)

    Unhappy about this design because :
        - the constructor is a mess
        - need to know about the internals to use efficiently
        - it's annoying to initialize

    An alternative is to store 2 fields;
        - the dense blocks (so store the zero opperators)
        - an informative (isscal,value)
=#
"
    MPOHamiltonian

    represents a general periodic quantum hamiltonian
"
struct MPOHamiltonian{S,T<:MPOTensor,E<:Number}<:Hamiltonian
    Os::PeriodicArray{Union{E,T},3}

    domspaces::PeriodicArray{S,2}
    pspaces::PeriodicArray{S,1}
end

function Base.getproperty(h::MPOHamiltonian,f::Symbol)
    if f==:odim
        return size(h.domspaces,2)
    elseif f==:period
        return size(h.pspaces,1)
    elseif f==:imspaces
        return PeriodicArray(circshift(adjoint.(h.domspaces),(-1,0)))
    else
        return getfield(h,f)
    end
end

#=
allow passing in
        - non strictly typed matrices
        - missing fields
        - 2leg tensors
        - only mpo tensors
=#

# bit of a helper - accept non strict typed data
MPOHamiltonian(x::AbstractArray{Any,3}) = MPOHamiltonian(union_split(x));
#another helper - artificially create a union and reuse next constructor
MPOHamiltonian(x::AbstractArray{T,3}) where T<: TensorMap = MPOHamiltonian(convert(AbstractArray{Union{T,eltype(T)},3},x));
function MPOHamiltonian(x::AbstractArray{T,3}) where T<:Union{A} where A
    (Sp,M,E) = _envsetypes(union_types(T));

    nx = similar(x,Union{E,M});

    for (i,t) in enumerate(x)
        if t isa MPSBondTensor
            nx[i] = permute(add_util_leg(t),(1,2),(4,3))
        elseif ismissing(t)
            nx[i] = zero(E)
        elseif t isa Number
            nx[i] = convert(E,t);
        else
            nx[i] = t;
        end
    end

    MPOHamiltonian{Sp,M,E}(nx);
end

#default constructor
function MPOHamiltonian{Sp,M,E}(x::AbstractArray{Union{E,M},3}) where {Sp,M<:MPOTensor{Sp},E<:Number}
    (period,numrows,numcols) = size(x);

    E == eltype(M) || throw(ArgumentError("scalar type should match mpo eltype $E ≠ $(eltype(M))"))
    numrows == numcols || throw(ArgumentError("mpos have to be square"))

    domspaces = PeriodicArray{Union{Missing,Sp}}(missing,period,numrows);
    pspaces = PeriodicArray{Union{Missing,Sp}}(missing,period)

    isused = fill(false,period,numrows,numcols);
    isstopped = false
    while !isstopped
        isstopped = true;

        for i = 1:period, j in 1:numrows, k in 1:numcols
            isused[i,j,k] && continue;

            if x[i,j,k] isa MPOTensor
                isused[i,j,k] = true;
                isstopped = false;

                #asign spaces when possible
                dom = space(x[i,j,k],1);im = space(x[i,j,k],3);p = space(x[i,j,k],2)

                ismissing(pspaces[i]) && (pspaces[i] = p);
                pspaces[i] != p && throw(ArgumentError("physical space for $((i,j,k)) incompatible : $(pspaces[i]) ≠ $(p)"))

                ismissing(domspaces[i,j]) && (domspaces[i,j] = dom)
                domspaces[i,j] != dom && throw(ArgumentError("Domspace for $((i,j,k)) incompatible : $(domspaces[i,j]) ≠ $(dom)"))

                ismissing(domspaces[i+1,k]) && (domspaces[i+1,k] = im')
                domspaces[i+1,k] != im' && throw(ArgumentError("Imspace for $((i,j,k)) incompatible : $(domspaces[i+1,k]) ≠ $(im')"))

                #if it's zero -> store zero
                #if it's the identity -> store identity
                if x[i,j,k] ≈ zero(x[i,j,k])
                    x[i,j,k] = zero(E) #the element is zero/missing
                else
                    ii,sc = isid(x[i,j,k])

                    if ii #the tensor is actually proportional to the identity operator -> store this knowledge
                        x[i,j,k] = sc
                    end
                end
            elseif x[i,j,k] != zero(E)
                if !ismissing(domspaces[i,j])
                    isused[i,j,k] = true;
                    isstopped = false;

                    ismissing(domspaces[i+1,k]) && (domspaces[i+1,k] = domspaces[i,j])
                    domspaces[i+1,k] != domspaces[i,j] && throw(ArgumentError("Identity incompatible at $((i,j,k)) : $(domspaces[i+1,k]) ≠ $(domspaces[i,j])"))
                elseif !ismissing(domspaces[i+1,k])
                    isused[i,j,k] = true;
                    isstopped = false;

                    ismissing(domspaces[i,j]) && (domspaces[i,j] = domspaces[i+1,k])
                    domspaces[i+1,k] != domspaces[i,j] && throw(ArgumentError("Identity incompatible at $((i,j,k)) : $(domspaces[i+1,k]) ≠ $(domspaces[i,j])"))
                end

            else
                isused[i,j,k] = true;
            end
        end
    end

    sum(ismissing.(pspaces)) == 0 || throw(ArgumentError("Not all physical spaces were assigned"))
    sum(ismissing.(domspaces)) == 0 || @warn "faied to deduce all domspaces"
    f_domspaces = map(x-> ismissing(x) ? oneunit(Sp) : x,domspaces) #missing domspaces => oneunit ; should also not happen

    ndomspaces = PeriodicArray{Sp}(f_domspaces)
    npspaces = PeriodicArray{Sp}(pspaces)

    return MPOHamiltonian{Sp,M,E}(PeriodicArray(x[:,:,:]),ndomspaces,npspaces)
end

#allow passing in regular tensormaps
MPOHamiltonian(t::TensorMap) = MPOHamiltonian(decompose_localmpo(add_util_leg(t)));

#a very simple utility constructor; given our "localmpo", constructs a mpohamiltonian
function MPOHamiltonian(x::Array{T,1}) where T<:MPOTensor{Sp} where Sp
    domspaces = Sp[space(y,1) for y in x]
    push!(domspaces,space(x[end],3)')

    pspaces=[space(x[1],2)]

    nOs = PeriodicArray{Union{eltype(T),T}}(fill(zero(eltype(T)),1,length(x)+1,length(x)+1))

    for (i,t) in enumerate(x)
        nOs[1,i,i+1]=t
    end

    nOs[1,1,1] = one(eltype(T));
    nOs[1,end,end] = one(eltype(T));

    ndomspace = PeriodicArray{Sp,2}(undef,1,length(x)+1);
    ndomspace[1,:] = domspaces[:]

    return MPOHamiltonian{Sp,T,eltype(T)}(nOs,ndomspace,PeriodicArray(pspaces))
end

function _envsetypes(d::Tuple)
    a = Base.first(d);
    b = Base.tail(d);

    if a <: MPOTensor
        return spacetype(a),a,eltype(a);
    elseif a <: MPSBondTensor
        return spacetype(a),tensormaptype(spacetype(a),2,2,eltype(a)),eltype(a)
    else
        @assert !isempty(b)
        return _envsetypes(b);
    end
end

#utility functions for finite mpo
function Base.getindex(x::MPOHamiltonian{S,T,E},a::Int,b::Int,c::Int)::T where {S,T,E}
    b <= x.odim && c <= x.odim || throw(BoundsError(x,[a,b,c]))
    if x.Os[a,b,c] isa E
        if x.Os[a,b,c] == zero(E)
            return TensorMap(zeros,E,x.domspaces[a,b]*x.pspaces[a],x.imspaces[a,c]'*x.pspaces[a])
        else
            return x.Os[a,b,c]*isomorphism(Matrix{E},x.domspaces[a,b]*x.pspaces[a],x.imspaces[a,c]'*x.pspaces[a])
        end
    else
        return x.Os[a,b,c]
    end
end

function Base.setindex!(x::MPOHamiltonian{S,T,E},v::T,a::Int,b::Int,c::Int)  where {S,T,E}
    b <= x.odim && c <= x.odim || throw(BoundsError(x,[a,b,c]))

    (ii,scal) = isid(v);

    if ii
        x.Os[a,b,c] = scal
    elseif v ≈ zero(v)
        x.Os[a,b,c] = zero(E)
    else
        x.Os[a,b,c] = v;
    end

    return x
end
Base.getindex(x::MPOHamiltonian,a::Colon,b::Int,c::Int) = [x[i,b,c] for i in 1:x.period];

Base.eltype(x::MPOHamiltonian) = typeof(x[1,1,1])
Base.size(x::MPOHamiltonian) = (x.period,x.odim,x.odim)
Base.size(x::MPOHamiltonian,i) = size(x)[i]

Base.keys(x::MPOHamiltonian) = Iterators.filter(a->contains(x,a[1],a[2],a[3]),Iterators.product(1:x.period,1:x.odim,1:x.odim))
Base.keys(x::MPOHamiltonian,i::Int) = Iterators.filter(a->contains(x,i,a[1],a[2]),Iterators.product(1:x.odim,1:x.odim))

opkeys(x::MPOHamiltonian) = Iterators.filter(a-> !isscal(x,a[1],a[2],a[3]),keys(x));
opkeys(x::MPOHamiltonian,i::Int) = Iterators.filter(a-> !isscal(x,i,a[1],a[2]),keys(x,i));

scalkeys(x::MPOHamiltonian) = Iterators.filter(a-> isscal(x,a[1],a[2],a[3]),keys(x));
scalkeys(x::MPOHamiltonian,i::Int) = Iterators.filter(a-> isscal(x,i,a[1],a[2]),keys(x,i));

Base.contains(x::MPOHamiltonian{S,T,E},a::Int,b::Int,c::Int) where {S,T,E} = !(x.Os[a,b,c] == zero(E))
isscal(x::MPOHamiltonian{S,T,E},a::Int,b::Int,c::Int) where {S,T,E} = x.Os[a,b,c] isa E && contains(x,a,b,c)

"
checks if ham[:,i,i] = 1 for every i
"
isid(ham::MPOHamiltonian{S,T,E},i::Int) where {S,T,E}= reduce((a,b) -> a && isscal(ham,b,i,i) && abs(ham.Os[b,i,i]-one(E))<1e-14,1:ham.period,init=true)

"
to be valid in the thermodynamic limit, these hamiltonians need to have a peculiar structure
"
function sanitycheck(ham::MPOHamiltonian)
    for i in 1:ham.period

        @assert isid(ham[i,1,1])[1]
        @assert isid(ham[i,ham.odim,ham.odim])[1]

        for j in 1:ham.odim
            for k in 1:(j-1)
                if contains(ham,i,j,k)
                    return false
                end
            end
        end
    end

    return true
end

"
checks if the given 4leg tensor is the identity (needed for infinite mpo hamiltonians)
"
function isid(x::MPOTensor)
    cod = space(x,1)*space(x,2);
    dom = space(x,3)'*space(x,4)';

    #would like to have an 'isisomorphic'
    for c in union(blocksectors(cod), blocksectors(dom))
        blockdim(cod, c) == blockdim(dom, c) || return false,0.0;
    end

    id = isomorphism(Matrix{eltype(x)},cod,dom)
    scal = dot(id,x)/dot(id,id)
    diff = x-scal*id

    scal = (scal ≈ 0.0) ? 0.0 : scal #shouldn't be necessary (and I don't think it is)

    return norm(diff)<1e-14,scal
end

include("utility.jl")
include("hamslice.jl")


text:
using FluxPrune
using Test

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


text:
################################################################################
# Spring Impulses
################################################################################

function spring_force(relative::Symbol, joint::Translational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion;
    unitary::Bool=false) where T

    spring = unitary ? 1.0 : joint.spring
    distance = joint.spring_offset .- minimal_coordinates(joint, xa, qa, xb, qb)
    force = spring * zerodimstaticadjoint(nullspace_mask(joint)) * distance

    return [force; szeros(T, 3)]
end

function spring_impulses(relative::Symbol, joint::Translational, pbody::Node, cbody::Node, timestep; unitary::Bool=false)
    spring_impulses(relative, joint,
        current_configuration(pbody.state)...,
        current_configuration(cbody.state)...,
        timestep, unitary=unitary)
end

function spring_impulses(relative::Symbol, joint::Translational,
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    timestep; unitary::Bool=false)
    timestep * impulse_transform(relative, joint, xa, qa, xb, qb) * spring_force(relative, joint, xa, qa, xb, qb; unitary=unitary)[SVector{3,Int}(1,2,3)]
end

spring_impulses(relative::Symbol, joint::Translational{T,3}, pbody::Node, cbody::Node, timestep; unitary::Bool=false) where T = szeros(T, 6)

################################################################################
# Spring Jacobians
################################################################################

function spring_force_jacobian_configuration(jacobian_relative::Symbol,
        joint::Translational{T},
        xa::AbstractVector, qa::Quaternion,
        xb::AbstractVector, qb::Quaternion;
        unitary::Bool=false,
        attjac=true) where T
    spring = unitary ? 1.0 : joint.spring
    return -spring * zerodimstaticadjoint(nullspace_mask(joint)) * minimal_coordinates_jacobian_configuration(jacobian_relative, joint, xa, qa, xb, qb, attjac=attjac)
end

function spring_jacobian_configuration(relative::Symbol, jacobian_relative::Symbol,
    joint::Translational{T},
    xa::AbstractVector, qa::Quaternion,
    xb::AbstractVector, qb::Quaternion,
    timestep::T;
    unitary::Bool=false,
    attjac=true) where T

    force = spring_force(relative, joint, xa, qa, xb, qb, unitary=unitary)[SVector{3,Int}(1,2,3)]

    J = impulse_transform(relative, joint, xa, qa, xb, qb) *
        spring_force_jacobian_configuration(jacobian_relative, joint, xa, qa, xb, qb, unitary=unitary, attjac=attjac)
    J += impulse_transform_jacobian(relative, jacobian_relative, joint, xa, qa, xb, qb, force, attjac=attjac)

    return timestep * J
end

function spring_jacobian_configuration(relative::Symbol, jacobian::Symbol,
    joint::Translational,
    pbody::Node, cbody::Node,
    timestep::T;
    attjac::Bool=true,
    unitary=false) where T
    return spring_jacobian_configuration(relative, jacobian,
        joint,
        current_configuration(pbody.state)...,
        current_configuration(cbody.state)...,
        timestep; unitary=false, attjac=attjac)
end

spring_jacobian_configuration(relative::Symbol, jacobian::Symbol, joint::Translational{T,3}, pbody::Node, cbody::Node, timestep::T; attjac::Bool = true, unitary::Bool=false) where T = attjac ? szeros(T, 6, 6) : szeros(T, 6, 7)
spring_jacobian_velocity(relative::Symbol, jacobian::Symbol, joint::Translational, pbody::Node, cbody::Node, timestep::T, unitary::Bool=false) where T = szeros(T, 6, 6)


text:
############################################################################################
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenStrict,
    x::AbstractArray{R}
) where {T<:AbstractFloat,F<:FlattenTypes,R<:Real}
    len = length(x)
    buffer_flat = zeros(T, size(vec(x)))
    function flatten_to_Real_Arr(v::AbstractArray{R}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return fill_array!(buffer_flat, v)
    end
    buffer_unflat = zeros(R, size(x))
    function unflatten_to_Real_Arr(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return fill_array!(buffer_unflat, v)
    end
    return flatten_to_Real_Arr, unflatten_to_Real_Arr
end
function construct_flatten(
    output::Type{T},
    flattentype::F,
    unflattentype::UnflattenFlexible,
    x::AbstractArray{R}
) where {T<:AbstractFloat,F<:FlattenTypes,R<:Real}
    len = length(x)
    function flatten_to_Real_Arr_AD(v::AbstractArray{R}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return fill_array!(zeros(R, length(x)), v)
    end
    function unflatten_to_Real_Arr_AD(v::Union{R,AbstractVector{R}}) where {R<:Real}
        ArgCheck.@argcheck length(v) == len
        return fill_array!(zeros(eltype(v), size(x)), v)
    end
    return flatten_to_Real_Arr_AD, unflatten_to_Real_Arr_AD
end

############################################################################################
#Export
export
    construct_flatten


text:
"""
    peakproms(peaks, x;
        strict=true,
        minprom=nothing,
        maxprom=nothing
    ) -> (peaks, proms)

Calculate the prominences of `peaks` in `x`, filtering peaks with prominences less than
`minprom` and greater than `maxprom`, if either are given.

Peak prominence is the absolute height difference between the current peak and the larger of
the two adjacent smallest magnitude points between the current peak and adjacent larger
peaks or signal ends.

The prominence for a peak with a `NaN` or `missing` between the current peak and either
adjacent larger peaks will be `NaN` or `missing` if `strict == true`, or it will be
the larger of the smallest non-`NaN` or `missing` values between the current peak and
adjacent larger peaks for `strict == false`.

See also: [`findminima`](@ref), [`findmaxima`](@ref)

# Examples
```jldoctest
julia> x = [0,5,2,3,3,1,4,0];

julia> xpks = argmaxima(x)
3-element Vector{Int64}:
 2
 4
 7

julia> peakproms(xpks, x)
([2, 4, 7], Union{Missing, Int64}[5, 1, 3])

julia> x = [missing,5,2,3,3,1,4,0];

julia> peakproms(xpks, x)
([2, 4, 7], Union{Missing, Int64}[missing, 1, 3])

julia> peakproms(xpks, x; strict=false)
([2, 4, 7], Union{Missing, Int64}[5, 1, 3])
```
"""
function peakproms(peaks::AbstractVector{Int}, x::AbstractVector{T};
    strict=true, minprom=nothing, maxprom=nothing, strictbounds=nothing
) where T
    if !isnothing(strictbounds)
        Base.depwarn("Keyword `strictbounds` has been renamed to `strict`", :peakproms!)
        strict=strictbounds
    end

    if !isnothing(minprom) || !isnothing(maxprom)
        _peaks = copy(peaks)
    else
        # peaks will not be modified
        _peaks = peaks
    end
    return peakproms!(_peaks, x; strict=strict, minprom=minprom, maxprom=maxprom)
end

"""
    peakproms!(peaks, x;
        strict=true,
        minprom=nothing,
        maxprom=nothing
    ) -> (peaks, proms)

Calculate the prominences of `peaks` in `x`, removing `peaks` with prominences less than
`minprom` or greater than `maxprom`, if either are given. Returns the modified peaks and
their prominences.

See also: [`peakproms`](@ref), [`findminima`](@ref), [`findmaxima`](@ref)
"""
function peakproms!(peaks::AbstractVector{Int}, x::AbstractVector{T};
    strict=true, minprom=nothing, maxprom=nothing
) where T
    if !isnothing(minprom) && !isnothing(maxprom)
        minprom < maxprom || throw(ArgumentError("minprom must be less than maxprom"))
    end
    all(∈(eachindex(x)), peaks) ||
        throw(ArgumentError("peaks contains invalid indices to x"))
    isempty(peaks) && return peaks, T[]

    # if peaks was calculated with strict=false, first(peaks) could be minima at firstindex
    fp = length(peaks) > 1 ? peaks[2] : first(peaks)
    if fp > 1 && ((x[fp] < x[fp-1]) === true)
        pktype = :minima
    else
        pktype = :maxima
    end
    cmp = pktype === :maxima ? (≥) : (≤)
    exm = pktype === :maxima ? minimum : maximum
    exa = pktype === :maxima ? max : min

    _ref = Missing <: T ? missing :
           Float64 <: T ? NaN :
           Float32 <: T ? NaN32 :
           Float16 <: T ? NaN16 :
                          missing

    proms = similar(peaks,promote_type(T,typeof(_ref)))

    if strict
        lbegin, lend = firstindex(x), lastindex(x)

        @inbounds for i in eachindex(peaks, proms)
            # Find left and right bound (self-intersections)
            lb = something(findprev(y -> cmp(y, x[peaks[i]]) === true, x, peaks[i] - 2),
                lbegin)
            rb = something(findnext(y -> cmp(y, x[peaks[i]]) === true, x, peaks[i] + 2),
                lend)

            # Find extremum of left and right bounds
            if isempty(lb:(peaks[i] - 1))
                lref = _ref
            else
                lref = exm(view(x, lb:(peaks[i] - 1)))
            end

            if isempty((peaks[i] + 1):rb)
                rref = _ref
            else
                rref = exm(view(x, (peaks[i] + 1):rb))
            end

            proms[i] = abs(x[peaks[i]] - exa(lref, rref))
        end
    else
        # The extremum search space in the bounding intervals can be reduced by
        # restricting the search space to known peaks/reverse peaks. The cost of
        # finding all peaks/reverse peaks should be mitigated by the fact that
        # the same peaks/reverse peaks will be the pivotal elements for
        # numerous peaks.
        if pktype === :maxima
            peaks′ = argmaxima(x, 1; strict=false)
            notm = argminima(x, 1; strict=false)
        else
            peaks′ = argminima(x, 1; strict=false)
            notm = argmaxima(x, 1; strict=false)
        end

        notmval = x[notm]

        for i in eachindex(peaks, proms)
            j = searchsorted(peaks′, peaks[i])

            # Find left and right bounding peaks
            _lb = findprev(y -> cmp(x[y], x[peaks[i]]) === true, peaks′, first(j)-1)
            peaks′[j] === peaks[i] && (j += 1)
            _rb = findnext(y -> cmp(x[y], x[peaks[i]]) === true, peaks′, last(j)+1)

            # Find left and right reverse peaks just inside the bounding peaks
            lb = isnothing(_lb) ? firstindex(notm) :
                                  searchsortedfirst(notm, peaks′[_lb])
            rb = isnothing(_rb) ? lastindex(notm) :
                                  searchsortedlast(notm, peaks′[_rb])

            k = searchsortedfirst(notm, peaks[i])

            if isempty(lb:(k-1))
                lref = missing
            else
                lref = exm(view(notmval, lb:(k-1)))
            end

            rb > lastindex(notm) && (rb -= 1)
            if isempty(k:rb)
                rref = missing
            else
                rref = exm(view(notmval, k:rb)) # Slice corollary upper side
            end

            proms[i] = abs(x[peaks[i]] - exa(coalesce(lref, rref), coalesce(rref, lref)))
        end
    end

    if !isnothing(minprom) || !isnothing(maxprom)
        lo = something(minprom, zero(eltype(x)))
        up = something(maxprom, typemax(Base.nonmissingtype(eltype(x))))
        matched = findall(x -> !ismissing(x) && !(lo ≤ x ≤ up), proms)
        deleteat!(peaks, matched)
        deleteat!(proms, matched)
    end

    return peaks, proms
end

struct Maxima; end
struct Minima; end

@deprecate peakprom(x::AbstractVector, w::Int=1; strictbounds=true, minprom=nothing) peakproms(argmaxima(x, w; strict=strictbounds), x; strict=strictbounds, minprom=minprom)
@deprecate peakprom(m::Minima, x::AbstractVector, w::Int=1; strictbounds=true, minprom=nothing) peakproms(argminima(x, w; strict=strictbounds), x; strict=strictbounds, minprom=minprom)
@deprecate peakprom(m::Maxima, x::AbstractVector, w::Int=1; strictbounds=true, minprom=nothing) peakproms(argmaxima(x, w; strict=strictbounds), x; strict=strictbounds, minprom=minprom)


text:
#### Entrainment-Detrainment kernels
using Flux

function compute_turbulent_entrainment(param_set, εδ_model_vars)
    FT = eltype(εδ_model_vars)
    c_t = FT(CPEDMF.c_t(param_set))
    K_ε = εδ_model_vars.a_up * c_t * sqrt(max(εδ_model_vars.tke, 0)) * εδ_model_vars.R_up

    ε_turb = if εδ_model_vars.w_up * εδ_model_vars.a_up > 0
        (2 / εδ_model_vars.R_up^2) * K_ε / (εδ_model_vars.w_up * εδ_model_vars.a_up)
    else
        FT(0)
    end

    return (ε_turb, K_ε)
end

function compute_inverse_timescale(param_set, εδ_model_vars)
    FT = eltype(εδ_model_vars)
    Δb = εδ_model_vars.b_up - εδ_model_vars.b_en
    Δw = get_Δw(param_set, εδ_model_vars)
    c_λ = FT(CPEDMF.c_λ(param_set))

    l_1 = c_λ * abs(Δb / sqrt(εδ_model_vars.tke + 1e-8))
    l_2 = abs(Δb / Δw)
    l = SA.SVector(l_1, l_2)
    return lamb_smooth_minimum(l, FT(0.1), FT(0.0005))
end

function get_Δw(param_set, εδ_model_vars)
    FT = eltype(εδ_model_vars)
    Δw = εδ_model_vars.w_up - εδ_model_vars.w_en
    Δw += copysign(FT(CPEDMF.w_min(param_set)), Δw)
    return Δw
end

function get_c_δ(param_set, εδ_model_vars)
    FT = eltype(εδ_model_vars)
    c_δ = if !TD.has_condensate(εδ_model_vars.q_cond_up + εδ_model_vars.q_cond_en)
        FT(0)
    else
        FT(CPEDMF.c_δ(param_set))
    end
    return c_δ
end

function get_MdMdz(εδ_model_vars)
    FT = eltype(εδ_model_vars)
    MdMdz = max(εδ_model_vars.dMdz / max(εδ_model_vars.M, eps(FT)), 0)
    MdMdz = max(-εδ_model_vars.dMdz / max(εδ_model_vars.M, eps(FT)), 0)
    return MdMdz
end

function dimensional_part(param_set, εδ_model_vars)
    Δw = get_Δw(param_set, εδ_model_vars)
    λ = compute_inverse_timescale(param_set, εδ_model_vars)
    return (λ / Δw)
end

function max_area_limiter(param_set, εδ_model_vars)
    FT = eltype(εδ_model_vars)
    γ_lim = FT(ICP.area_limiter_scale(param_set))
    β_lim = FT(ICP.area_limiter_power(param_set))
    logistic_term = (2 - 1 / (1 + exp(-γ_lim * (εδ_model_vars.max_area - εδ_model_vars.a_up))))
    area_limiter = dimensional_part(param_set, εδ_model_vars) * (logistic_term^β_lim - 1)
    return area_limiter
end

function non_dimensional_groups(param_set, εδ_model_vars)
    Δw = get_Δw(param_set, εδ_model_vars)
    Δb = εδ_model_vars.b_up - εδ_model_vars.b_en
    Π_1 = Δw^2 / (εδ_model_vars.zc_i * Δb)
    Π_2 = Δw^2 / εδ_model_vars.tke
    Π_3 = √(εδ_model_vars.a_up)
    Π_4 = εδ_model_vars.RH_up - εδ_model_vars.RH_en
    return SA.SVector(Π_1, Π_2, Π_3, Π_4)
end

"""
    entr_detr(param_set, εδ_model_vars, εδ_model_type::MDEntr)

Returns the dynamic entrainment and detrainment rates,
as well as the turbulent entrainment rate, following
Cohen et al. (JAMES, 2020), given:
 - `param_set`      :: parameter set
 - `εδ_model_vars`  :: structure containing variables
 - `εδ_model_type`  :: MDEntr - Moisture deficit entrainment closure
"""
function entr_detr(param_set, εδ_model_vars, εδ_model_type::MDEntr)
    FT = eltype(εδ_model_vars)
    dim_scale = dimensional_part(param_set, εδ_model_vars)
    area_limiter = max_area_limiter(param_set, εδ_model_vars)

    # Moisture deficit closure
    Δw = get_Δw(param_set, εδ_model_vars)
    D_ε, D_δ, M_δ, M_ε = nondimensional_exchange_functions(param_set, Δw, εδ_model_vars)
    c_ε = FT(CPEDMF.c_ε(param_set))
    c_δ = get_c_δ(param_set, εδ_model_vars)
    nondim_ε = (c_ε * D_ε + c_δ * M_ε)
    nondim_δ = (c_ε * D_δ + c_δ * M_δ)
    c_div = FT(ICP.entrainment_massflux_div_factor(param_set))
    MdMdz = get_MdMdz(εδ_model_vars)
    drybubble_adj = MdMdz * c_div

    # dynamic entrainment / detrainment
    ε_dyn = dim_scale * nondim_ε + drybubble_adj
    δ_dyn = dim_scale * nondim_δ + drybubble_adj + area_limiter

    # turbulent entrainment
    ε_turb, K_ε = compute_turbulent_entrainment(param_set, εδ_model_vars)

    return EntrDetr(ε_dyn, δ_dyn, ε_turb, K_ε)
end

"""
    entr_detr(param_set, εδ_model_vars, εδ_model_type::NNEntr)

Returns the dynamic entrainment and detrainment rates,
as well as the turbulent entrainment rate, using a neural
network to predict the non-dimensional component of dynamical entrainment:
 - `param_set`      :: parameter set
 - `εδ_model_vars`  :: structure containing variables
 - `εδ_model_type`  :: NNEntr - Neural network entrainment closure
"""
function entr_detr(param_set, εδ_model_vars, εδ_model_type::NNEntr)
    c_gen = ICP.c_gen(param_set)
    dim_scale = dimensional_part(param_set, εδ_model_vars)
    area_limiter = max_area_limiter(param_set, εδ_model_vars)

    # Neural network closure
    nn_arc = (4, 2, 2)  # (#inputs, #neurons, #outputs)
    nn_model = Chain(
        Dense(reshape(c_gen[1:8], nn_arc[2], nn_arc[1]), c_gen[9:10], sigmoid),
        Dense(reshape(c_gen[11:14], nn_arc[3], nn_arc[2]), c_gen[15:16], softplus),
    )

    nondim_groups = non_dimensional_groups(param_set, εδ_model_vars)
    nondim_ε, nondim_δ = nn_model(nondim_groups)

    # dynamic entrainment / detrainment
    ε_dyn = dim_scale * nondim_ε
    δ_dyn = dim_scale * nondim_δ + area_limiter

    # turbulent entrainment
    ε_turb, K_ε = compute_turbulent_entrainment(param_set, εδ_model_vars)

    return EntrDetr(ε_dyn, δ_dyn, ε_turb, K_ε)
end


text:
#######################################################################################################################################################################################################
#
# Changes to this struct
# General
#     2021-Aug-10: refactor the structure with renamed fields
#     2021-Aug-10: add a constructor within the structure to avoid external initialization
#     2021-Oct-19: sort variable to prognostic and dignostic catergories
#     2022-Jan-24: fix documentation
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure that stores wave length information.

# Fields

$(TYPEDFIELDS)

"""
mutable struct WaveLengthSet{FT<:AbstractFloat}
    # parameters that do not change with time
    "Indicies of Λ_NIR in Λ"
    IΛ_NIR::Vector{Int}
    "Indicies of Λ_PAR in Λ"
    IΛ_PAR::Vector{Int}
    "Indicies of Λ_SIF in Λ"
    IΛ_SIF::Vector{Int}
    "Indicies of Λ_SIFE in Λ"
    IΛ_SIFE::Vector{Int}
    "Number of wavelength bins"
    NΛ::Int
    "Number of wavelength bins for PAR"
    NΛ_PAR::Int
    "Number of wavelength bins for SIF"
    NΛ_SIF::Int
    "Number of wavelength bins for SIF excitation"
    NΛ_SIFE::Int
    "Wavelength limits for NIR `[nm]`"
    WL_NIR::Vector{FT}
    "Wavelength limits for PAR `[nm]`"
    WL_PAR::Vector{FT}
    "Wavelength limits for SIF emission `[nm]`"
    WL_SIF::Vector{FT}
    "Wavelength limits for SIF excitation `[nm]`"
    WL_SIFE::Vector{FT}
    "Standard wavelength (boundaries) `[nm]`"
    SΛ::Vector{FT}
    "Differential wavelength `[nm]`"
    ΔΛ::Vector{FT}
    "Differential wavelength for PAR `[nm]`"
    ΔΛ_PAR::Vector{FT}
    "Differential wavelength for SIF excitation `[nm]`"
    ΔΛ_SIFE::Vector{FT}
    "Wavelength (bins) `[nm]`"
    Λ::Vector{FT}
    "Wavelength bins for PAR `[nm]`"
    Λ_PAR::Vector{FT}
    "Wavelength bins for SIF `[nm]`"
    Λ_SIF::Vector{FT}
    "Wavelength bins for SIF excitation `[nm]`"
    Λ_SIFE::Vector{FT}
end

#######################################################################################################################################################################################################
#
# Changes to this constructor
# General
#     2021-Nov-24: isolate the constructor from structure
#     2022-Jan-24: fix documentation
#
#######################################################################################################################################################################################################
"""

    WaveLengthSet{FT}(swl::Vector=WAVELENGTHS; opti::String=OPTI_2021) where {FT<:AbstractFloat}

Constructor for [`WaveLengthSet`](@ref), given
- `swl` Standard wave length boundaries `[nm]`, default is `ClimaCache.WAVELENGTHS`
- `opti` Optical properties file path, default is `ClimaCache.OPTI_2021`

---
# Examples
```julia
wls = WaveLengthSet{Float64}();
wls = WaveLengthSet{Float64}(collect(400:5:2500));
wls = WaveLengthSet{Float64}(collect(400:5:2500); opti=ClimaCache.OPTI_2017);
```
"""
WaveLengthSet{FT}(swl::Vector=WAVELENGTHS; opti::String=OPTI_2021) where {FT<:AbstractFloat} = (
    _dwl    = diff(swl);
    _λ      = zeros(FT, length(swl)-1);
    _opti   = matread(opti)["optipar"];
    _λ_opti = _opti["wl"];
    @inbounds for _i in 1:length(swl)-1
        _wo    = findall( swl[_i] .<= _λ_opti .< swl[_i+1] );
        _λ[_i] = mean(_λ_opti[_wo]);
    end;
    _iλ_nir  = findall( 700 .<= _λ .<= 2500 );
    _iλ_par  = findall( 400 .<= _λ .<= 700  );
    _iλ_sif  = findall( 640 .<= _λ .<= 850  );
    _iλ_sife = findall( 400 .<= _λ .<= 750  );

    return WaveLengthSet{FT}(_iλ_nir, _iλ_par, _iλ_sif, _iλ_sife, length(_λ), length(_iλ_par), length(_iλ_sif), length(_iλ_sife), FT[700,2500], FT[400,700], FT[640,850], FT[400,750], swl, _dwl,
                             _dwl[_iλ_par], _dwl[_iλ_sife], _λ, _λ[_iλ_par], _λ[_iλ_sif], _λ[_iλ_sife])
);


text:
# utilities for x ↔ scale factor ↔ redshift
a2z(a::T) where T = one(T)/a - one(T)
z2a(z::T) where T = one(T)/(z + one(T))
a2x(a) = log(a)
x2a(x) = exp(x)
z2x(z) = a2x(z2a(z))
x2z(x) = a2z(x2a(x))

# utility function for constructing an interpolator
spline(x, y) = scale(interpolate(y, BSpline(Cubic(Line(OnGrid())))), x)
spline_∂ₓ(f, x_grid) = spline(x_grid, [Interpolations.gradient(f, x)[1] for x in x_grid])
spline_∂ₓ²(f, x_grid) = spline(x_grid, [Interpolations.hessian(f, x)[1] for x in x_grid])

δ_kron(i, j) = (i == j) ? 1 : 0


text:
# Command: Move To
"""
    moveto!(session::Session; x::Integer = 0, y::Integer = 0)
    moveto!(element::Element; x::Integer = 0, y::Integer = 0)

Move the mouse by an offset of the specificed element.

!!! note

    This command uses the JSON Wire Protocol instead of the current W3c WebDriver API.
"""
function moveto!(session::Session; x::Integer = 0, y::Integer = 0)
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/moveto!",
        [("Content-Type" => "application/json")],
        JSON3.write(Dict("x" => x, "y" => y)),
    )
    @assert response.status == 200
    nothing
end
function moveto!(element::Element; x::Integer = 0, y::Integer = 0)
    @unpack addr, id = element.session
    response = HTTP.post(
        "$addr/session/$id/moveto",
        [("Content-Type" => "application/json")],
        JSON3.write(Dict("element" => element.id, "x" => x, "y" => y)),
    )
    @assert response.status == 200
    nothing
end


text:
"""
    change_basis(basis,changeofbasis::AbstractMatrix)

# Examples

Compute the Lagrangian basis associated with a set of nodes

```jldoctests
using LinearAlgebra
using InplaceArrays.Fields
using InplaceArrays.Polynomials

D = 2
order = 1
f = MonomialBasis{D}(Float64,order)

nodes = Point{2,Int}[(0,0),(1,0),(0,1),(1,1)]
change = inv(evaluate(f,nodes))

g = change_basis(f,change)
println(evaluate(g,nodes))

# output
[1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]

```
"""
function change_basis(basis,changeofbasis::AbstractMatrix)
  BasisFromChangeOfBasis(basis,changeofbasis)
end

struct BasisFromChangeOfBasis{B,M} <: Field
  basis::B
  change::M
  function BasisFromChangeOfBasis(basis,change::AbstractMatrix)
    B = typeof(basis)
    M = typeof(change)
    new{B,M}(basis,change)
  end
end

function field_cache(b::BasisFromChangeOfBasis,x)
  cb = field_cache(b.basis,x)
  bx = evaluate_field!(cb,b.basis,x)
  c = CachedArray(bx*b.change)
  (c,cb)
end

function evaluate_field!(cache,b::BasisFromChangeOfBasis,x)
  c, cb = cache
  bx = evaluate_field!(cb,b.basis,x)
  setsize!(c,size(bx))
  mul!(c.array,bx,b.change)
  c.array
end

function gradient_cache(b::BasisFromChangeOfBasis,x)
  cb = gradient_cache(b.basis,x)
  bx = evaluate_gradient!(cb,b.basis,x)
  c = CachedArray(bx*b.change)
  (c,cb)
end

function evaluate_gradient!(cache,b::BasisFromChangeOfBasis,x)
  c, cb = cache
  bx = evaluate_gradient!(cb,b.basis,x)
  setsize!(c,size(bx))
  mul!(c.array,bx,b.change)
  c.array
end

function hessian_cache(b::BasisFromChangeOfBasis,x)
  cb = hessian_cache(b.basis,x)
  bx = evaluate_hessian!(cb,b.basis,x)
  c = CachedArray(bx*b.change)
  (c,cb)
end

function evaluate_hessian!(cache,b::BasisFromChangeOfBasis,x)
  c, cb = cache
  bx = evaluate_hessian!(cb,b.basis,x)
  setsize!(c,size(bx))
  mul!(c.array,bx,b.change)
  c.array
end


text:
module InputInts

Base.Experimental.@optlevel 1

using Redux
using CImGui
import ..ReduxImGui: get_label, get_values, get_value

# actions
abstract type AbstractInputIntAction <: AbstractSyncAction end

"""
    Rename(label, new_label)
Change label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractInputIntAction
    label::String
    new_label::String
end

"""
    ChangeSize(label, size)
Change widget size to `size`.
Note changing size also changes widget appearance. If the old size is bigger,
those unused states are left intact.
"""
struct ChangeSize <: AbstractInputIntAction
    label::String
    size::Int
end

"""
    ChangeStep(label, step)
Change value's moving step to `step`.
This only takes effect when the widget size is 1.
"""
struct ChangeStep <: AbstractInputIntAction
    label::String
    step::Int
end

"""
    ChangeSpeed(label, speed)
Change value's moving speed to `speed`.
This only takes effect when the widget size is 1.
"""
struct ChangeSpeed <: AbstractInputIntAction
    label::String
    speed::Int
end

# state
struct State <: AbstractImmutableState
    label::String
    vals::Vector{Cint}
    size::Int
    step::Int
    speed::Int
    flags::CImGui.ImGuiInputTextFlags_
end
State(
    label::AbstractString;
    vals::Vector{Cint} = Cint[0, 0, 0, 0],
    size = 1,
    step = 1,
    speed = 100,
    flags = CImGui.ImGuiInputTextFlags_None,
) = State(label, vals, size, step, speed, flags)

State(label::AbstractString, val::Integer; args...) =
    State(label; vals = Cint[val, 0, 0, 0], args...)

# reducers
input_int(state::AbstractState, action::AbstractAction) = state
input_int(state::Vector{<:AbstractState}, action::AbstractAction) = state
input_int(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

function input_int(state::Dict{String,State}, action::AbstractInputIntAction)
    s = Dict{String,State}()
    for (k, v) in state
        s[k] = get_label(v) == action.label ? input_int(v, action) : v
    end
    return s
end

input_int(s::State, a::Rename) = State(a.new_label, s.vals, s.size, s.step, s.speed, s.flags)
input_int(s::State, a::ChangeSize) = State(s.label, s.vals, a.size, s.step, s.speed, s.flags)
input_int(s::State, a::ChangeStep) = State(s.label, s.vals, s.size, a.step, s.speed, s.flags)
input_int(s::State, a::ChangeSpeed) = State(s.label, s.vals, s.size, s.step, a.speed, s.flags)

# helper
"""
    InputInt(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
"""
function InputInt(store::AbstractStore, get_state = Redux.get_state)
    s = get_state(store)
    if s.size == 1
        is_triggered = CImGui.InputInt(get_label(s), s.vals, s.step, s.speed, s.flags)
    elseif s.size == 2
        is_triggered = CImGui.InputInt2(get_label(s), s.vals, s.flags)
    elseif s.size == 3
        is_triggered = CImGui.InputInt3(get_label(s), s.vals, s.flags)
    elseif s.size == 4
        is_triggered = CImGui.InputInt4(get_label(s), s.vals, s.flags)
    else
        is_triggered = false
    end
    return is_triggered
end

get_label(s::State) = s.label
get_values(s::State) = s.vals
get_value(s::State) = first(s.vals)

end # module


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ ea80a138-0847-11ec-3f65-ed1ea62e1f73
html"""

Macro 318

 Julia for aspiring economists

"""

# ╔═╡ ab0e1d94-5f52-4f33-bcc3-51a60cb0cec5
md""" # Introduction """

# ╔═╡ f3d6382d-418e-472a-bfc7-f7d0a67fc842
md""" In this notebook we will give a broad overview of the Julia programming language and how it can be useful for you as an aspiring economist. Do not be too concerned about the language and syntax, most of you are new to the idea of programming and will find this too be potentially scary at first. """

# ╔═╡ Cell order:
# ╟─ea80a138-0847-11ec-3f65-ed1ea62e1f73
# ╟─ab0e1d94-5f52-4f33-bcc3-51a60cb0cec5
# ╟─f3d6382d-418e-472a-bfc7-f7d0a67fc842


text:
#module A6

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\A6.jl"),
  informat="script",
  out_path = "$(pwd())\\A6.html",
  doctype = "md2html")
=#
if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt, HCubature, StaticArrays, CTNumerical

const WORKERS = 8 #seems like logical cores * (3/4) is the magic value
const CLEAN_WORKERS = true

if nworkers() < WORKERS || (nworkers()==1 && WORKERS == 1)
  addprocs(nworkers()==1?WORKERS-nworkers()+1:WORKERS-nworkers())
end

@everywhere using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTModCopy, JuMP, NLopt, HCubature, StaticArrays, CTNumerical

@everywhere if !isdefined(:constDef)
  const constDef = true
  const DATA_PATH = pwd() * "\\data" #path where we store the data
  const OUTPUT_PATH = pwd() * "\\output"
  const ALT_DATE_FORMAT_STR = "yyyymmdd"
  const DATE_FORMAT_STR = "m/d/y"
  const SP500_NAME = "GSPC"
  const TBILL_NAME = "t-bill"
  const TERM_NAME = "termStructure"
  const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
  const HEADER_NAME = "header.tex"
  const OPTION_STD_NAME = "VIXoptionsStd"
  const OPTION_NAME = "VIXoptions"

  const PlotContainer = Union{Plot,Gadfly.Compose.Context}
  const CSRegSpec = Tuple{Symbol, Symbol, Int}
  const CPRegSpec = Tuple{Symbol, Int}

  const MIN_VALS_FOR_EST = 10
  const Ι = Complex(0.0,1.0)
  const STANDARD_DAYS = [30, 60, 91, 122]
end

#pre-process the s&p500 data
function preProcessSP500()::Void
  dateFormat::String = DATE_FORMAT_STR
  sp500DF::DataFrame = readtable("$DATA_PATH\\$SP500_NAME.csv")

  #fix the dates
  rename!(sp500DF, :Date, :DATE)
  sp500DF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(sp500DF[:,:DATE])::DataVector{Date}
  sort!(sp500DF, cols = [:DATE])
  sp500DF[:,:returns] = similar(sp500DF[:,:Adj_Close])
  sp500DF[2:end,:returns] = log.(sp500DF[2:end,:Adj_Close] ./  sp500DF[1:(end-1),:Adj_Close])

  sp500DF[:,:dayOfQuarter] = zeros(Int,size(sp500DF,1))
  for i::Int ∈ 2:size(sp500DF,1)
    if Dates.firstdayofquarter(sp500DF[i,:DATE]) ≠ Dates.firstdayofquarter(sp500DF[i-1,:DATE])
      sp500DF[i,:dayOfQuarter] = 1
    else
      sp500DF[i,:dayOfQuarter] = sp500DF[i-1,:dayOfQuarter] + 1
    end
  end
  #showcols(sp500DF)

  stream::IOStream = open("$DATA_PATH\\$SP500_NAME.jls", "w")
  serialize(stream, sp500DF)
  close(stream)

  return nothing
end

#a simple routine for making an easy to read dataframe binary
function preProcessTBills()::Void
  dateFormat::String = DATE_FORMAT_STR
  tBillDF::DataFrame = readtable("$DATA_PATH\\$TBILL_NAME.csv")

  #fix the dates
  tBillDF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(tBillDF[:,:DATE])::DataVector{Date}
  sort!(tBillDF, cols = [:DATE])
  tBillDF[:,:TB3MS] .*= .01
  tBillDF[:,:rates]= tBillDF[:,:TB3MS] ./12.0 #(1.0+tBillDF[:,:TB3MS]).^(1.0/12.0) .- 1.0

  stream::IOStream = open("$DATA_PATH\\$TBILL_NAME.jls", "w")
  serialize(stream, tBillDF)
  close(stream)

  return nothing
end

#pre-process standard options data
function preProcessVIXOptionsStd(;verbose::Bool = false)::Void
  dateFormat::String = ALT_DATE_FORMAT_STR
  stdDF::DataFrame = readtable("$DATA_PATH\\$OPTION_STD_NAME.csv")

  stdDays::SVector{length(STANDARD_DAYS), Int} = SVector{length(STANDARD_DAYS), Int}(STANDARD_DAYS)

  #adjust date type
  stdDF[:, :date] = ((i::Int)->Date("$i", dateFormat)::Date).(stdDF[:,:date])::DataVector{Date}
  rename!(stdDF, :date, :DATE)

  stdDF = stdDF[completecases(stdDF[:,[:DATE, :days, :forward_price, :delta, :cp_flag]]),:]

  if verbose
    println("total record")
    showcols(stdDF)

    for d::Int ∈ unique(stdDF[:,:days])
      println("records for $d")
      showcols(stdDF[stdDF[:,:days].==d,:])
    end
  end

  maximumTau::Int = maximum(STANDARD_DAYS)
  stdDF = stdDF[stdDF[:,:days].≤maximumTau,:]

  #make this of type char
  delete!(stdDF, [:secid, :ticker, :index_flag])
  stdDF[:, :cp_flag] = ((s::String)->s[1]).(stdDF[:,:cp_flag])::DataVector{Char}
    sort!(stdDF, cols = [:DATE, :cp_flag, :days])

  #showcols(stdDF)

  stream::IOStream = open("$DATA_PATH\\$OPTION_STD_NAME.jls", "w")
  serialize(stream, stdDF)
  close(stream)

  return nothing

end

#pre-process standard options data
function preProcessVIXOptions()::Void
  #constants that will be used later
  dateFormat::String = ALT_DATE_FORMAT_STR
  optionDF::DataFrame = readtable("$DATA_PATH\\$OPTION_NAME.csv")

  #adjust date type
  optionDF[:, :date] = ((i::Int)->Date("$i", dateFormat)::Date).(optionDF[:,:date])::DataVector{Date}
  rename!(optionDF, :date, :DATE)

  optionDF[:, :exdate] = ((i::Int)->Date("$i", dateFormat)::Date).(optionDF[:,:exdate])::DataVector{Date}
  rename!(optionDF, :exdate, :exDATE)

  maximumTau::Int = maximum(STANDARD_DAYS)
  optionDF[:,:days] = (Dates.value).(optionDF[:, :exDATE])  .- (Dates.value).(optionDF[:, :DATE])
  optionDF = optionDF[(!isna).(optionDF[:,:days]), :]
  optionDF = optionDF[optionDF[:,:days] .≤ maximumTau, :]
  optionDF[:, :τ] = optionDF[:,:days] ./ 365.0

  #sort and re-scale
  optionDF[:, :strike_price] = optionDF[:, :strike_price] ./ 1000.0
  sort!(optionDF, cols = [:DATE, :exDATE, :cp_flag, :strike_price])

  optionDF[:, :price] = (optionDF[:, :best_offer] .+ optionDF[:, :best_bid]) ./ 2.0

  #filter out bad data
  optionDF[:, :cp_flag] = ((s::String)->s[1]).(optionDF[:,:cp_flag])::DataVector{Char}
  optionDF = optionDF[completecases(optionDF[:,[:price, :delta, :cp_flag]]),:]
  delete!(optionDF, [:best_offer, :best_bid, :secid, :ticker, :index_flag, :issuer, :exercise_style])

  #showcols(optionDF)

  stream::IOStream = open("$DATA_PATH\\$OPTION_NAME.jls", "w")
  serialize(stream, optionDF)
  close(stream)

  return nothing

end

#this gets a mapping of a given day to a standard day, and the weight for interpolation
function getInterpolationDayIndex()::Tuple{Vector{Int}, Vector{Float64}}
  standardDays::Vector{Int} = STANDARD_DAYS
  maxStandardDays = maximum(standardDays)

  #allocate the return containers, which include the mapping and the weights
  standardDayIndex::Vector{Int} = Vector{Int}(maximum(standardDays))
  weights::Vector{Float64} = Vector{Float64}(maxStandardDays)

  #iterate through all of the standardized days and map them to a standard day
  standardCtr::Int = 1
  for i::Int = 1:maxStandardDays
    if i>standardDays[standardCtr]
      standardCtr += 1
    end
    standardDayIndex[i] = standardCtr

    #calculate the weight
    weights[i] = standardCtr==1? 0.0 :
      (i-standardDays[standardCtr-1]) / (standardDays[standardCtr] - standardDays[standardCtr-1])

  end

  return (standardDayIndex, weights)
end

#merge the optionDF and stdDF
function mergeOptionDFs!(optionDF::DataFrame, stdDF::DataFrame)::DataFrame

  #calculate information about the standard time periods
  standardDays::Vector{Int} = STANDARD_DAYS
  minStandardDays = minimum(standardDays)
  maxStandardDays = maximum(standardDays)
  standardDayIndex::Vector{Int}, weights::Vector{Float64} = getInterpolationDayIndex();

  #index up the standard maturities (will help with merge)
  Dateτ::DataVector{Tuple{Date, Int, Char}} =
    ((d::Date, f::Int, c::Char)->(d,f,c)).(stdDF[:,:DATE],stdDF[:,:days], stdDF[:,:cp_flag])

  stdDict::Dict = Dict(Dateτ[i]=>i for i::Int ∈ 1:size(stdDF,1))
  FStd::DataVector{Float64} =  stdDF[:, :forward_price]

  #for each date and ex date pair
  N::Int =size(optionDF,1)

  optionDF[:,:F] = DataVector{Float64}(Vector{Float64}(N))
  optionDF[:,:F] .= NA
  optionDF[:,:invalid] = false

  #interpolate the forward price
  by(optionDF, :DATE) do dateDF::SubDataFrame
    by(dateDF, :exDATE) do exDateDF::SubDataFrame
      by(exDateDF, :cp_flag) do subDF::SubDataFrame

        period::Int = min(subDF[1,:days], maxStandardDays)
        current::Date = subDF[1,:DATE]
        cp::Char = subDF[1,:cp_flag]

        #now calculate the interpolated value
        if !haskey(stdDict, (current, standardDays[standardDayIndex[period]], cp))
          subDF[:,:invalid] = true
        elseif period ∈ standardDays #if no need for interpolation
          subDF[:,:F] = FStd[stdDict[(current, period, cp)]]
        elseif standardDayIndex[period] == 1
          subDF[:,:F] = FStd[stdDict[(current, standardDays[1], cp)]]
        else
          #take the weighted average of the nearest values
          subDF[:,:F] =
            FStd[stdDict[(current, standardDays[standardDayIndex[period]], cp)]] * weights[period] +
            FStd[stdDict[(current, standardDays[standardDayIndex[period]-1], cp)]] * (1.0-weights[period])
        end
      end
    end
  end

  optionDF = optionDF[(!isna).(optionDF[:,:F]),:]
  optionDF = optionDF[!optionDF[:,:invalid], :]
  delete!(optionDF, :invalid)
  N = size(optionDF,1)

  optionDF[:, :M] = optionDF[:, :F] ./ optionDF[:, :strike_price]

  #now get the high/low metric for the VIX

  medianST::Float64 = median(stdDF[stdDF[:,:days].==minStandardDays,:forward_price])
  medianLT::Float64 = median(stdDF[stdDF[:,:days].==maxStandardDays,:forward_price])

  optionDF[:,:LT_VIX] = Vector{Symbol}(N)
  optionDF[:,:LT_VIX] = NA

  optionDF[:,:ST_VIX] = Vector{Symbol}(N)
  optionDF[:,:ST_VIX] = NA

  by(optionDF, :DATE) do dateDF::SubDataFrame
    by(dateDF, :cp_flag) do subDF::SubDataFrame
      #appropriate standard indices

      #create a view for calculating if the VIX is high or low
      stdSubDF =
        view(stdDF, (stdDF[:,:DATE] .== subDF[1,:DATE]) .& (stdDF[:,:cp_flag] .== subDF[1,:cp_flag]))

      #compute if high or low
      maxDays::Int = maximum(stdSubDF[:,:days])
      minDays::Int = minimum(stdSubDF[:,:days])
      LTVIX::Float64 = (stdSubDF[stdSubDF[:,:days].==maxDays,:forward_price])[1]
      STVIX::Float64 = (stdSubDF[stdSubDF[:,:days].==minDays,:forward_price])[1]

      #update the dataframe
      subDF[:,:LT_VIX] = LTVIX>medianLT?(:H):(:L)
      subDF[:,:ST_VIX] = STVIX>medianLT?(:H):(:L)

    end
  end

  return optionDF
end

#this function merges in the risk-free rate
function mergeOptionDFρ!(optionDF::DataFrame, tBillDF::DataFrame)::DataFrame

  rates::Vector{Float64} = tBillDF[:,:TB3MS]
  rateDates::DataVector{Date} = tBillDF[:,:DATE]
  rateDict::Dict = Dict(rateDates[i]=>rates[i] for i::Int ∈ 1:length(rates))

  #allocate the interest rate column
  optionDF[:, :ρ] = DataVector{Float64}(Vector{Float64}(size(optionDF,1)))
  optionDF[:, :ρ] .= NA
  #for each date, look up the interest rate
  by(optionDF, :DATE) do subDF
    dt::Date = Dates.firstdayofmonth(subDF[1, :DATE])
    if haskey(rateDict, dt)
      subDF[:, :ρ] = rateDict[dt]
    end
  end

  #get rid of invalid interest rate data
  optionDF = optionDF[(!isna).(optionDF[:,:ρ]), :]
  optionDF[:,:ρ] = (log).(1.0 .+ optionDF[:,:ρ])

  return optionDF

end

#convert all puts to calls
function conformToCalls!(optionDF::DataFrame)::DataFrame
  putDF::SubDataFrame = view(optionDF, optionDF[:,:cp_flag].=='P')

  #use put call parity to convert all puts to calls
  putDF[:, :price] += ((exp).(-1.0 .* putDF[:,:τ] .* putDF[:,:ρ])) .*
    (putDF[:,:F] .- putDF[:,:strike_price])

  delete!(optionDF, :cp_flag)

  return optionDF
end

#black model with no dividends
@everywhere function getImpliedBlackσ(price::Float64, F::Float64, τ::Float64, X::Float64,
  r::Float64, startσ::Float64 = 0.15)::Float64
  #Algorithms: :LD_LBFGS, :LN_COBYLA, :LN_BOBYQA
  #  :LD_SLSQP, :LD_TNEWTON, :LD_VAR2/:LD_VAR1,
  #  :LD_MMA, :GN_ESCH, :GN_DIRECT, GN_DIRECT_L, :GD_STOGO
  # GN_ISRES, :LD_LBFGS, :LN_SBPLX,
  #GN_CRS2_LM
  #:GN_CRS2_LM***, :GN_ISRES*** (slow), :LD_MMA**, LN_SBPLX**, LD_TNEWTON** best:62,848

  B::Float64 = exp(-1.0*r*τ)

  function obj(σ::Real)::Real
    d1 = (log(F/X) + τ*(σ^2/2.0))/ (σ*τ^0.5)

    calcPrice::Real = (F*cdf(Normal(),d1) - X*cdf(Normal(),d1-σ*τ^0.5))*B
    return (price-calcPrice)*(price-calcPrice)
  end

  #println("test obj: $(obj(0.3))")

  #works well:LN_NELDERMEAD , LN_SBPLX, LD_SLSQP, LN_BOBYQA
  mod = Model(solver=NLoptSolver(algorithm=:LD_SLSQP, maxtime=1))
  JuMP.register(mod, :obj, 1, obj, autodiff=true)

  #constrain vol to 0 and 1000%
  @variable(mod, 3.0 >= σ >= 10.^-6.0, start = startσ)

  #set the objective function
  @NLobjective(mod, Min, obj(σ))

  σBest::Float64  = -99.0
  objBest::Float64 = -9999.0

  try
    output = solve(mod)
    σBest  = getvalue(σ)
  catch err
    σBest = -99.0
  end

  if σBest ≤   2.0*10.0 ^ -6.0
    σBest = -99.0
  end

  return σBest
end

#calcualte the IV for all options in the dataframe
function calculateIVs!(optionDF)::DataFrame
  N::Int = size(optionDF,1)

  Cs::Array{Float64} = optionDF[:, :price]
  Fs::Array{Float64} = optionDF[:, :F]
  startσs::Array{Float64} = optionDF[:, :impl_volatility]
  τs::Array{Float64} = optionDF[:, :τ]
  Xs::Array{Float64} = optionDF[:, :strike_price]
  ρs::Array{Float64} = optionDF[:, :ρ]

  #pre-allocate the answer space
  σs::Array{Float64} = Vector{Float64}(N)
  #  σs .= NA

  #=@sync @parallel for i  ∈  1:N
    σs[i] = getImpliedBlackσ(Cs[i], Fs[i], τs[i], Xs[i], ρs[i], startσs[i])
  end=#

  σs .= pmap(getImpliedBlackσ, Cs, Fs, τs, Xs, ρs, startσs)

  optionDF[:,:σ] = DataVector{Float64}(Vector{Float64}(N))
  optionDF[:,:σ] .= ((f::Float64)->f<0.0?NA:f).(σs)
  println("Solution failures: $(length(optionDF[isna.(optionDF[:,:σ]),:σ]))
    out of $(length(optionDF[:,:σ]))")

  return optionDF
end

#this encapsulates the processing of the option DF
function processOptionData!(optionDF::DataFrame, stdDF::DataFrame, tBillDF::DataFrame)
  optionDF = mergeOptionDFs!(optionDF, stdDF) #get the merged data
  optionDF = mergeOptionDFρ!(optionDF, tBillDF) #merge in interest rates
  optionDF = conformToCalls!(optionDF) #convert puts to calls

  println("beginning implied volatility calculation")
  optionDF =  calculateIVs!(optionDF) #calculate the implied volatilities

  #write this out to a file
  stream::IOStream = open("$DATA_PATH\\$OPTION_NAME-merged.jls", "w")
  serialize(stream, optionDF)
  close(stream)
end

#makes the volatility smile graph
function makeSmileGraphs(optionDF::DataFrame; label::Symbol = :full,
    binScaleM::Float64 = 25.0, binScaleDays::Float64 = 0.1)::Void

  ####first smile plot
  plotDF::DataFrame = deepcopy(optionDF)
  delete!(plotDF, [:DATE, :exDATE, :strike_price, :F, :ρ, :impl_volatility, :τ, :price, :ST_VIX, :LT_VIX])
  completecases!(plotDF)

  plotDFNoBin::DataFrame = deepcopy(plotDF)

  plotDF[:,:M] = ((f::Float64)->round(f*binScaleM)/binScaleM).(plotDF[:,:M])
  plotDF[:,:days] = ((i::Int)->Int(round(i*binScaleDays)/binScaleDays)).(plotDF[:,:days])
  plotDF=aggregate(plotDF, [:days, :M], mean)

  #get the graphs
  Gadfly.push_theme(:default)
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  push!(plotNames, "multiSmile")
  push!(plots, plot(plotDF, x=:M, y=:σ_mean, color=:days,
    Geom.line,
    Guide.ylabel("σ"), Guide.xlabel("F/X"),
    Coord.Cartesian(ymin=0.0,ymax=2.5, xmin=0.5, xmax=2.0),
    Guide.title("Volatility Smile  by Time and Moneyness")))

  push!(plotNames, "multiSmileDot")
  push!(plots, plot(plotDF, x=:M, y=:σ_mean, color=:days, shape=:days,
    Geom.point,
    Guide.ylabel("σ"), Guide.xlabel("F/X"),
    Coord.Cartesian(ymin=0.0,ymax=2.5, xmin=0.5, xmax=2.0),
    Guide.title("Volatility Smile  by Time and Moneyness")))

  for i ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i])_$label.svg", 4inch, 4inch),plots[i])
  end

  return nothing
end

function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

function makeAllSmileGraphs(optionDF::DataFrame)
  makeSmileGraphs(optionDF, label = :full_smooth, binScaleM = 10.0, binScaleDays=0.05)
  makeSmileGraphs(optionDF, label = :full_fine, binScaleM = 1000.0, binScaleDays=0.025)

  makeSmileGraphs(optionDF[(optionDF[:,:ST_VIX] .== :H) .& (optionDF[:,:LT_VIX] .== :H), :],
    label = :HH, binScaleM = 10.0, binScaleDays=0.05)
  makeSmileGraphs(optionDF[(optionDF[:,:ST_VIX] .== :L) .& (optionDF[:,:LT_VIX] .== :H), :],
    label = :LH, binScaleM = 10.0, binScaleDays=0.05)
  makeSmileGraphs(optionDF[(optionDF[:,:ST_VIX] .== :H) .& (optionDF[:,:LT_VIX] .== :L), :],
    label = :HL, binScaleM = 10.0, binScaleDays=0.05)
  makeSmileGraphs(optionDF[(optionDF[:,:ST_VIX] .== :L) .& (optionDF[:,:LT_VIX] .== :L), :],
    label = :LL, binScaleM = 10.0, binScaleDays=0.05)
end

#script for getting and processing the data
function A6Script(;refreshData::Bool=true, reprocessOptions=true)::Void

  #read the binary file
  tBillDF::DataFrame = getDF("$DATA_PATH\\$TBILL_NAME.jls", refreshData, preProcessTBills)
  sp500DF::DataFrame = getDF("$DATA_PATH\\$SP500_NAME.jls", refreshData, preProcessSP500)
  stdDF::DataFrame = getDF("$DATA_PATH\\$OPTION_STD_NAME.jls", refreshData, preProcessVIXOptionsStd)
  optionDF::DataFrame = getDF("$DATA_PATH\\$OPTION_NAME.jls", refreshData, preProcessVIXOptions)

  #we need to create an anonymous function to use the getDF interface
  #mergeOptionDFs!(deepcopy(optionDF),deepcopy(stdDF))
  optionDF = getDF("$DATA_PATH\\$OPTION_NAME-merged.jls", reprocessOptions,
    ()->processOptionData!(optionDF, stdDF,tBillDF))

  makeAllSmileGraphs(optionDF)

  #print the summary stats
  #preProcessVIXOptionsStd(verbose = true)

  #write out the processed option file for testing purposes
  writetable("$OUTPUT_PATH\\$OPTION_NAME.csv", optionDF)
  return nothing
end

function cleanworkers()
  myProcs = procs()
  if CLEAN_WORKERS
    for i in 2:length(myProcs)
      rmprocs(myProcs[i])
    end
  end
end

@time begin
  A6Script(refreshData=false, reprocessOptions=true)
  cleanworkers()
end

#end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# bessel_minus_struve.jl

module bessel_minus_struve
export BesselI_minus_StruveL

using ArbNumerics

"""
Implementation of I[n, z] - L[-n, z] where I[n, z] is the modified Bessel function of the first kind, and L[-n, z] is the modified Struve function. n is an integer ≥ 0 and the order of the functions, and z∈ℜ>0 a positive real argument. Notice the order are opposite signs.

The expansion of the BesselI and StruveL functions are:

    I[n, z] = ∑_{k=0}^∞ (z/2)^{2k+n} / (Γ(k+1)⋅Γ(k+n+1))
    L[-n, z] = ∑_{k=0}^∞ (z/2)^{2k-n+1} / (Γ(k+3/2)⋅Γ(k-n+3/2))

with Γ(x) the Gamma function. The values of I[n, z] and L[-n, z] can be incredibly close, so ArbReal types are used as they allow for accurate arbitrary precision floats so that the difference, I[n, z] - L[-n, z], can be distinguished. ArbReal types are faster than BigFloat types, and the error of a float can be tracked using ball(::ArbReal)[2].

Unfortunately, as ArbReal(bits = 64) * ArbReal(bits = 128) = ArbReal(bits = 64), if the starting precision is too small to accurately produce I[n, z] - L[-n, z], then the entire summation has to be restarted at a higher precision otherwise rounding errors will be too large, which slows the algorithm sometimes. Therefore, if the precision is found to be too low, the precision is just doubled (rather than increased by some set amount) to reach a required precision quickly.

The general structure of this arbitrary precision summation algorithm is used elsewhere too since you just change the arguments and term appropriately.
"""

function BesselI_minus_StruveL(n, z, a; prec = 64) # z > 0 & n >= 0

    # Initialise precision of ArbReal to prec.
    p = prec
    setextrabits(0)
    setprecision(ArbReal, p) # ArbReal(bit = 64 + 8) has same precision as Float64 (which is why we add 8 bits).

    n = ArbReal("$n")
    z = ArbReal("$z")
    a = ArbReal("$a")

    k = ArbReal("0.0")
    result = ArbReal("0.0")
    term = ArbReal("1.0")
    err = eps(result)  # Machine accuracy of specified precision = prec.

    while abs(midpoint(term)) > err * abs(midpoint(result))
        term = ArbReal((abs(z) * a / 2)^(2 * k + n + 1) / (ArbNumerics.gamma(k + 1) * ArbNumerics.gamma(k + n + 2)) - (abs(z) * a / 2)^(2 * k - n) / (ArbNumerics.gamma(k + 3//2) * ArbNumerics.gamma(k - n + 1//2)))
        result += term
        # println("term: k = ", k, "\nterm value: ", ball(term), "\ncumulant result: ", ball(result), "\n")
        k += 1
        # Double precision if rounding error in result exceeds accuracy specified by prec.
        while radius(result) > err * abs(midpoint(result))
            p *= 2
            setprecision(ArbReal, p)
            # variables have to be re-parsed into the higher precision ArbReal type.
            # println("Not precise enough. Error = ", radius(result), " > ", ArbReal("$err"), ". Increasing precision to ", p, " bits.\n")
            n = ArbReal("$n")
            z = ArbReal("$z")
            a = ArbReal("$a")
            k = ArbReal("0")
            result = ArbReal("0.0")
            term = ArbReal("1.0")
        end
    end
    # println("z: ", ArbReal(z, bits = prec), ". Final result: ", ArbReal(result, bits = prec))
    ArbReal(result * √π * ArbNumerics.gamma(-n - 1/2) * (abs(z) / 2)^n * z / 2, bits = prec) # return to specified precision.
end
# @time bsi = BesselI_minus_StruveL(3, 55, 1)
# @show(bsi)

# β = ArbReal("2.0")
# α = ArbReal("7.0")
# v = ArbReal("5.8")
# w = ArbReal("1.6")
# R = ArbReal((v^2 - w^2) / (w^2 * v))
# a = ArbReal(sqrt(β^2 / 4 + R * β * coth(β * v / 2)))
# z = ArbReal("90.61")
# n = ArbReal("12.0")
# c = BesselI_minus_StruveL(n, z, a; prec = 24)

end # end of module


text:
# surface(::Matrix, ::Matrix, ::Matrix)
function _default(main::Tuple{MatTypes{T}, MatTypes{T}, MatTypes{T}}, s::Style{:surface}, data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        position_x = main[1] => (Texture, "x position, must be a `Matrix{Float}`")
        position_y = main[2] => (Texture, "y position, must be a `Matrix{Float}`")
        position_z = main[3] => (Texture, "z position, must be a `Matrix{Float}`")
        scale = Vec3f0(0) => "scale must be 0, for a surfacemesh"
    end
    surface(position_z, s, data)
end

# surface(Vector or Range, Vector or Range, ::Matrix)
function _default(main::Tuple{VectorTypes{T}, VectorTypes{T}, MatTypes{T}}, s::Style{:surface}, data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        position_x = main[1] => (Texture, "x position, must be a `Vector{Float}`")
        position_y = main[2] => (Texture, "y position, must be a `Vector{Float}`")
        position_z = main[3] => (Texture, "z position, must be a `Matrix{Float}`")
        scale = Vec3f0(0) => "scale must be 0, for a surfacemesh"
    end
    surface(position_z, s, data)
end

# surface(::Matrix)
function _default(main::MatTypes{T}, s::Style{:surface}, data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        ranges = ((-1f0, 1f0), (-1f0,1f0)) => "x, and y position given as `(start, endvalue)` or any `Range`"
    end
    delete!(data, :ranges) # no need to have them in the OpenGL data
    _default((Grid(to_value(main), to_value(ranges)), main), s, data)
end

# surface(::Matrix)
function _default(main::Tuple{G, MatTypes{T}}, s::Style{:surface}, data::Dict) where {G <: Grid{2}, T <: AbstractFloat}
    xrange = main[1].dims[1]; yrange = main[1].dims[2]
    xscale = (maximum(xrange) - minimum(xrange)) / (length(xrange)-1)
    yscale = (maximum(yrange) - minimum(yrange)) / (length(yrange)-1)
    @gen_defaults! data begin
        position    = main[1] =>" Position given as a `Grid{2}`.
        Can be constructed e.g. `Grid(LinRange(0,2,N1), LinRange(0,3, N2))`"
        position_z  = main[2] => (Texture, "height offset for the surface, must be `Matrix{Float}`")
        scale       = Vec3f0(xscale, yscale, 1) => "scale of the grid planes forming the surface. Can be made smaller, to let the grid show"
    end
    surface(position_z, s, data)
end

_extrema(x::FRect3D) = Vec2f0(minimum(x)[3], maximum(x)[3])
nothing_or_vec(x) = x
nothing_or_vec(x::Array) = vec(x)

function normal_calc(x::Bool, invert_normals::Bool = false)
    i = invert_normals ? "-" : ""
    if x
        "$(i)getnormal(position, position_x, position_y, position_z, o_uv);"
        # "getnormal_fast(position_z, ind2sub(dims, index1D));"
    else
        "vec3(0, 0, $(i)1);"
    end
end

function light_calc(x::Bool)
    if x
        """
        vec3 L      = normalize(o_lightdir);
        vec3 N      = normalize(o_normal);
        vec3 light1 = blinnphong(N, o_camdir, L, color.rgb);
        color       = vec4(light1, color.a); // + light2 * 0.4
        """
    else
        ""
    end
end

function native_triangle_mesh(mesh)
    return gl_convert(triangle_mesh(mesh))
end

function surface(main, s::Style{:surface}, data::Dict)
    @gen_defaults! data begin
        primitive = Rect2D(0f0,0f0,1f0,1f0) => native_triangle_mesh
        scale = nothing
        position = nothing
        position_x = nothing => Texture
        position_y = nothing => Texture
        position_z = nothing => Texture
        wireframe = false
        glow_color = RGBA{Float32}(0,0,0,0) => GLBuffer
        stroke_color = RGBA{Float32}(0,0,0,1) => GLBuffer
        stroke_width = wireframe ? 0.03f0 : 0f0
        glow_width = 0f0
        uv_offset_width = Vec4f0(0) => GLBuffer
        shape = RECTANGLE
        image = nothing => Texture
        distancefield = nothing => Texture
        shading = true
        normal = shading
        invert_normals = false
    end
    @gen_defaults! data begin
        color = nothing => Texture
        color_map = nothing => Texture
        color_norm = nothing
        fetch_pixel = false
        matcap = nothing => Texture
        uv_scale = Vec2f0(1)
        instances = const_lift(x->(size(x,1)-1) * (size(x,2)-1), main) => "number of planes used to render the surface"
        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "surface.vert",
            to_value(wireframe) ? "distance_shape.frag" : "standard.frag",
            view = Dict(
                "position_calc" => position_calc(position, position_x, position_y, position_z, Texture),
                "normal_calc" => normal_calc(normal, to_value(invert_normals)),
                "light_calc" => light_calc(shading),
            )
        )
    end
end

function position_calc(x...)
    _position_calc(Iterators.filter(x->!isa(x, Nothing), x)...)
end

function _position_calc(
        position_x::MatTypes{T}, position_y::MatTypes{T}, position_z::MatTypes{T}, target::Type{Texture}
    ) where T<:AbstractFloat
    """
    int index1D = index + offseti.x + offseti.y * dims.x + (index/(dims.x-1));
    ivec2 index2D = ind2sub(dims, index1D);
    vec2 index01 = vec2(index2D) / (vec2(dims)-1.0);
    pos = vec3(
        texelFetch(position_x, index2D, 0).x,
        texelFetch(position_y, index2D, 0).x,
        texelFetch(position_z, index2D, 0).x
    );
    """
end

function _position_calc(
        position_x::VectorTypes{T}, position_y::VectorTypes{T}, position_z::MatTypes{T},
        target::Type{Texture}
    ) where T<:AbstractFloat
    """
    int index1D = index + offseti.x + offseti.y * dims.x + (index/(dims.x-1));
    ivec2 index2D = ind2sub(dims, index1D);
    vec2 index01 = vec2(index2D) / (vec2(dims)-1.0);
    pos = vec3(
        texelFetch(position_x, index2D.x, 0).x,
        texelFetch(position_y, index2D.y, 0).x,
        texelFetch(position_z, index2D, 0).x
    );
    """
end

function _position_calc(
        position_xyz::VectorTypes{T}, target::Type{TextureBuffer}
    ) where T <: StaticVector
    "pos = texelFetch(position, index).xyz;"
end

function _position_calc(
        position_xyz::VectorTypes{T}, target::Type{GLBuffer}
    ) where T <: StaticVector
    len = length(T)
    filler = join(ntuple(x->0, 3-len), ", ")
    needs_comma = len != 3 ? ", " : ""
    "pos = vec3(position $needs_comma $filler);"
end

function _position_calc(
        grid::Grid{2}, position_z::MatTypes{T}, target::Type{Texture}
    ) where T<:AbstractFloat
    """
    int index1D = index + offseti.x + offseti.y * dims.x + (index/(dims.x-1));
    ivec2 index2D = ind2sub(dims, index1D);
    vec2 index01 = vec2(index2D) / (vec2(dims)-1.0);
    float height = texture(position_z, index01).x;
    pos = vec3(grid_pos(position, index01), height);
    """
end


text:
function __init__()
    @require TaylorModels = "314ce334-5f6e-57ae-acf6-00b6e903104a" load_taylormodels()
end

function load_taylormodels()
    eval(load_taylormodels_overapproximation())
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Point light sources used by point-light tracer

@doc raw"""
    PointLight

A point light (used by [`PointLightRenderer`](@ref)).

This type holds information about a point light.

# Fields

- `position::Point`: a [`Point`](@ref) object holding the position of the point light in 3D space.
- `color::RGB{Float32}`: the color of the point light.
- `linear_radius::Float32`: radius of the source, used to compute solid angle subtended by the light.

If `linear_radius` is non-zero, it is used to compute the solid angle subtended by the light at a given
distance `d` through the formula:

```math
\left(\frac{\mathrm{linear\_radius}}{d}\right)^2
```
"""
Base.@kwdef struct PointLight
    position::Point = ORIGIN
    color::RGB{Float32} = WHITE
    linear_radius::Float32 = 0f0
end

@doc """
    PointLight(position::Point, color::RGB{Float32}, linear_radius::Float32)

Constructor for a [`PointLight`](@ref) instance.
""" PointLight(::Point, ::RGB{Float32}, ::Float32)

@doc """
    PointLight(; position::Point = ORIGIN,
                 color::RGB{Float32} = WHITE,
                 linear_radius::Float32 = 0f0)

Constructor for a [`PointLight`](@ref) instance.

If no parameter is specified, it return a white point light in the origin with no radius.
""" PointLight(; ::Point, ::RGB{Float32}, ::Float32)

"""
    Lights

Alias of `Vector{PointLight}`, to store a list of [`PointLight`](@ref) sources.
"""
const Lights = Vector{PointLight}