text:
"""
    declaration(cinfo, func, argtypes, Rtype)

Get a wat string with function declaration.
"""
function declaration(cinfo, func, argtypes, Rtype)
    Nparams = length(argtypes)
    slotnames = cinfo.slotnames[2:end]
    slottypes = cinfo.slottypes[2:end]
    slottypes = get_slottypes(slottypes)

    decl = ["func \$$func (export \"$func\")"]
    locals = []
    for i = 1:length(slotnames)
        if i <= Nparams
            push!(decl, "(param \$$(slotnames[i]) $(slottypes[i]))")
        else
            push!(locals, "(local \$$(slotnames[i]) $(slottypes[i]))")
        end
    end

    ret = []
    if length(Rtype.parameters) > 0
        #isarray = Rtype == Array{Float64,1}
        isarray = Rtype <: Array
        if isarray
            push!(ret, "(result i32)")
        else
            # DataType Tuple
            for rtype in Rtype.parameters
                if !(rtype <: Nothing)
                    rt = rtype <: AbstractFloat ? "f32" : "i32"
                    push!(ret, "(result $rt)")
                end
            end
        end
    else
        # DataType
        if !(Rtype <: Nothing)
            rt = Rtype <: AbstractFloat ? "f32" : "i32"
            push!(ret, "(result $rt)")
        end
    end

    if length(locals)>0
        return join([decl; ret; "\n"; locals], " ")
    else
        return join([decl; ret], " ")
    end
end

function get_slottypes(slottypes)
    st = []
    for slottype in slottypes
        if slottype <: AbstractFloat
            push!(st, "f32")
        elseif length(slottype.parameters) > 1 && istuple(slottype)
            t = slottype.parameters[1] <: AbstractFloat ? "f32" : "i32"
            push!(st, t)
        else
            push!(st, "i32")
        end
    end
    return st
end

"""
    translate(i::Integer, cinfo::CodeInfo, item)

Get a wat string from item, specialized on item type.
"""
translate(i::Integer, ci::CodeInfo, item) = item
translate(i::Integer, ci::CodeInfo, item::AbstractFloat) = "(f32.const $item)"
translate(i::Integer, ci::CodeInfo, item::Number) = "(i32.const $item)"
# translate(i::Integer, ci::CodeInfo, item::Nothing) = "(i32.const 0)"
translate(i::Integer, ci::CodeInfo, item::Bool) = item ? "(i32.const 1)" : "(i32.const 0)"
translate(i::Integer, ci::CodeInfo, item::SlotNumber) = "(local.get \$$(ci.slotnames[item.id]))"
translate(i::Integer ,ci::CodeInfo, item::GlobalRef) = "call \$$(item.name)"
translate(i::Integer ,ci::CodeInfo, item::NewvarNode) = nothing
translate(i::Integer, ci::CodeInfo, item::Compiler.Const) = translate(i, ci, item.val)
function translate(i::Integer, ci::CodeInfo, item::TypedSlot)
    if isa(item.typ, Compiler.Const)
        return translate(i, ci, item.typ.val[1])
    else
        return item
    end
end

"""
    translate(ci::CodeInfo, items::Array)

Get a wat string from items, branching to special cases based on items[1].
"""
function translate(i::Integer, ci::CodeInfo, items::Array)
    # dont translate?
    if isnothing(items)
        return nothing

    # builtin to webassembly?
    elseif is_floatop(ci, items)
        return ["$(floatops[items[1].name])"; translate(i, ci, items[2:end])]
    elseif is_intop(ci, items)
        return ["$(intops[items[1].name])"; translate(i, ci, items[2:end])]

    # special expression?

    elseif hasname(items[1], :(size)) && length(items)==3
        return ["call \$size$(items[3])"; translate(i, ci, items[2])]
    elseif hasname(items[1], :(tuple))
        #ignore construction of tuples...
        return [translate.((i,), (ci,), items[2:end])]
    elseif hasname(items[1], :(getindex)) && length(items)==3
        if isa(items[2],SlotNumber) && istuple(ci.slottypes[items[2].id])
            #allow indexing into tuples
            id = items[2].id
            fieldnumber = items[3]
            name = fieldnumber == 1 ? "$(ci.slotnames[id])" : "$(ci.slotnames[id])$(fieldnumber)"
            return "(local.get \$$name)"
        else
            return ["call \$getlinearindex"; translate(i, ci, items[2:end])]
        end
    elseif hasname(items[1], :(setindex!)) && length(items)==4
        return ["call \$setlinearindex"; translate(i, ci, items[2:end])]
    elseif hasname(items[1], :(ifelse))
        return ["select"; translate(i, ci, items[2:end])]
    elseif hasname(items[1], Symbol("return"))
        if hasname(items[2], Symbol("nothing"))
            return ["return"]
        elseif isa(items[2], Array) && hasname(items[2][1], :(tuple))
            return ["return"; translate(i, ci, items[2][2:end])]
        else
            return ["return"; translate(i, ci, items[2:end])]
        end
    elseif hasname(items[1], :(not_int))
        return ["i32.eqz", translate(i, ci, items[2])]
    elseif hasname(items[1], Symbol("==="))
        if typeof(items[3]) <: Nothing
            if isa(items[2], SlotNumber) && string(ci.slotnames[items[2].id])[1] == '_' # iterator variable begins with "_"
                return ["i32.eqz", "(local.get \$$(ci.slotnames[items[2].id])i)"] # the iteratorindex
            elseif itemtype(ci, items[2]) <: AbstractFloat
                return ["f32.eq", translate(i, ci, items[2]), "(f32.const 0.0)"]
            else
                return ["i32.eqz", translate(i, ci, items[2])]
            end
        else
            return ["i32.eq", translate(i, ci, items[2:3])]
        end
    # elseif hasname(items[1], :(=))
    #    slotname = ci.slotnames[items[2].id]
    #    return ["local.set \$$(slotname)", translate(i, ci, items[3])]

    elseif hasname(items[1], :(=))
        slotname = ci.slotnames[items[2].id]
        if isa(items[3], Array) && is_iterate(items[3][1])
            return ["local.set \$$(slotname) ( local.set \$$(slotname)i", translate(i, ci, items[3]), ")"]
        elseif isa(items[3], Array) && istuple(ci.slottypes[items[2].id])
            params = ci.slottypes[items[2].id].parameters
            settuples = []
            closetuples = []
            for n = 2:length(params)
                push!(settuples,"(local.set \$$(slotname)$(n)")
                push!(closetuples, ")")
            end
            settuples = spacedjoin(settuples)
            closetuples = spacedjoin(closetuples)
            return ["local.set \$$(slotname)",settuples, translate(i, ci, items[3]), closetuples]
        else
            return ["local.set \$$(slotname)", translate(i, ci, items[3])]
        end
    # default individual translation
    else
        return translate.((i,), (ci,), items)
    end
end

is_iterate(item) = isa(item, GlobalRef) && length(string(item.name)) >= 7 && string(item.name)[1:7] == "iterate"
is_floatop(ci::CodeInfo, items) = isa(items[1], GlobalRef) && items[1].name in keys(floatops) && (hasitemtype(ci, items[2], AbstractFloat) || (length(items) > 2 && hasitemtype(ci, items[3], AbstractFloat)))
is_intop(ci::CodeInfo, items) = isa(items[1], GlobalRef) && items[1].name in keys(intops) && (hasitemtype(ci, items[2], [Integer, Bool, Nothing]) || (length(items) > 2 && hasitemtype(ci, items[3], [Integer, Bool, Nothing])))

floatops = Dict(
:(+) => "f32.add", # these consume 2 args
:(-) => "f32.sub",
:(*) => "f32.mul",
:(/) => "f32.div",
:(==) => "f32.eq",
:(!=) => "f32.ne",
:(<) => "f32.lt",
:(>) => "f32.gt",
:(<=) => "f32.le",
:(>=) => "f32.ge",
:(min) => "f32.min",
:(max) => "f32.max",
:(copysign) => "f32.copysign",
:(abs) => "f32.abs", # these consume 1 arg
:(ceil) => "f32.ceil",
:(floor) => "f32.floor",
:(trunc) => "f32.trunc",
:(round) => "f32.nearest",
:(sqrt) => "f32.sqrt",
:(float) => "f32.convert_i32_s",
:(Int) => "i32.trunc_f32_s",
# :(^) => ["call \$pow",2],
# ""=>["f32.neg",1],
# ""=>["f32.load",1],
# ""=>["f32.store",2],
# ""=>["f32.const",1],
)

intops = Dict(
:(+) => "i32.add", # these consume 2 args
:(-) => "i32.sub",
:(*) => "i32.mul",
:(div) => "i32.div_s", # julia 4/2 will return float. but div(4,2) wont
:(%) => "i32.rem_s",
:(%) => "i32.rem_s",
:(rem) => "i32.rem_s",
:(&&) => "i32.and",
:(||) => "i32.or",
:(⊻) => "i32.xor",
:(xor) => "i32.xor",
:(<<) => "i32.shl",
:(>>) => "i32.shr_s",
:(==) => "i32.eq",
# :(===) => "i32.eq", #handle this in a special way.
:(!=) => "i32.ne",
:(<) => "i32.lt_s",
:(>) => "i32.gt_s",
:(<=) => "i32.le_s",
:(>=) => "i32.ge_s",
:(!) => "i32.eqz", # these consume 1 arg
:(not_int) => "i32.eqz",
:(leading_zeros) => "i32.clz",
:(trailing_zeros) => "i32.ctz",
:(count_ones) => "i32.popcnt",
:(float) => "f32.convert_i32_s",
:(Int) => "i32.trunc_f32_s",
# :(^) => ["call \$powi",2],
# ""=>["i32.load",1]
# ""=>["i32.store",2]
# ""=>["i32.const",1]
)


text:
using Revise
src_dir = "C:/Users/stravaglino3/Downloads/julia_code/src/"
# includet(src_dir*"model_fit_global.jl")
includet(src_dir*"model_fit_global_v2.jl")
includet(src_dir*"models.jl")
includet(src_dir*"helper_functions.jl")
includet(src_dir*"global_plotting.jl")
using Plots
plot(0:10)
# --------------------------create dummy data-----------------------------------
# f = (du,u,p,t) -> bimolecular!(du,u,p,t,dens,cons) # enclose constants
#
# # u0 = [0.0;0.0]
# u₀ = [0.0]          # Initial condition
# tspan = (0.0,10.0)  # Simulation interval
# p = [1e-4,0.5]      # equation parameter. p = [k₁, k₋₁]
# dens = [100., 10., 1.]
# cons = []
#
# # Setup the ODE problem, then solve
# prob = ODEProblem(f,u₀,tspan,p)
# sol = solve(prob,Tsit5())
# # sol = solve(prob,Tsit5())
# plot(sol)
# dataset1 = generate_data(sol)
#
# dens = [50., 35., 1.]
# prob = ODEProblem(f,u₀,tspan,p)
# sol = solve(prob,Tsit5())
# # sol = solve(prob,Tsit5())
# plot!(sol)
# dataset2 = generate_data(sol)
#
# model_opt2 = (1.0e-6,10.0)
# # pinit = [2.0e-4,0.6]
# pinit = [-4.0,0.0]
# t = collect(range(0,stop=10,length=101))
# dens_comb = [[100., 10., 1.],[50., 35., 1.]]
# rates,hes,l = optimization(bimolecular!,[vec(dataset1),vec(dataset2)],
#                             dens_comb,[[],[]],model_opt2,[t,t])

#----------------------------------------------------------------------------
# try TCR-pMHC only data from Muaz

using CSV

data_dir = "C:/Users/stravaglino3/Downloads/julia_code/TCR_CD4_project/data/"
tcr_data = CSV.read(data_dir*"tcr_only_data.csv", DataFrame)

@time t, Pa, n_se_tcr = clean_data(tcr_data)
n = Pa2n(Pa)

# scatter(t,Pa)
# using Plots
scatter(t,n)
scatter!(n_se_tcr[:,1],n_se_tcr[:,2],yerror=n_se_tcr[:,3],markercolor=:black)

densities_tcr = [[25.0, 38.0, 1.0]]
tspan = (0,maximum(t)+1)
u₀ = [0.0]
bounds = (1.0e-6,100.0)
pinit = [-4.0,0.0]
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_tcr,[]) # enclose constants

pinit = [-4.5,-0.7]

tcr = Data([n],[t])
M_tcr = M(bimolecular!,u₀,tspan,densities_tcr,[[]],0.0,Rodas5)

# rates_tcr,hes_tcr,l_tcr = optimization(bimolecular!,[n],[densities_tcr],[[]],model_opt2,[t])
using BlackBoxOptim
rates_tcr,hes_tcr,l_tcr = optimization(M_tcr,tcr,bounds)
# σ_tcr = hessian2σ(hes_tcr)
σ_tcr = hessian2σ(hes_tcr,l_tcr,tcr.n[1],rates_tcr)
# u₀ = [0.0,0.]
# rates_tcr,hes_tcr,l_tcr = optimization(trimolecular!,n,densities_tc1',rates_tc,model_opt2,t)
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_tcr[1],[])
plot_sum(rates_tcr)

# unique_t = [unique(times) for times in tcr.t]
# loss_in = (p) -> loss(p,tcr,M_tcr,unique_t)
# hes_for =ForwardDiff.hessian(loss_in,rates_tcr)
# σ_tcr = hessian2σ(hes_for)
# plot_curve_sds(hes_tcr,l_tcr,n,rates_tcr)
# calculate_sds(hes_tcr,l_tcr,n,rates_tcr)
# prob = ODEProblem(f,u₀,tspan,10 .^ rates_tcr)
# # prob = remake(prob, p=10 .^ rates)
# sol = solve(prob,Rodas4())
# plot!(sol)
# # fig = plot!(sol,linecolor=:black)
# summ = reduce(vcat,sum(sol,dims=1))
# h = plot!(sol.t,summ,linecolor=:black)

# prob = ODEProblem(f,u₀,tspan,[1e-4,0.5])

#------------------------------------------------------------------------------
# try TCR-pMHC-CD4 fit with 2-path, 2-step

tcr_cd4_1 = CSV.read(data_dir*"tcr_cd4_1.csv", DataFrame)
tcr_cd4_2 = CSV.read(data_dir*"tcr_cd4_2.csv", DataFrame)
tcr_cd4_3 = CSV.read(data_dir*"tcr_cd4_3.csv", DataFrame)
t_tc1, Pa_tc1, n_se_1 = clean_data(tcr_cd4_1)
t_tc2, Pa_tc2,  n_se_2  = clean_data(tcr_cd4_2)
t_tc3, Pa_tc3,  n_se_3 = clean_data(tcr_cd4_3)

n_tc1 = Pa2n(Pa_tc1)
n_tc2 = Pa2n(Pa_tc2)
n_tc3 = Pa2n(Pa_tc3)

# scatter(t,Pa)
scatter(t_tc1,n_tc1)
scatter!(t_tc2,n_tc2)
scatter!(t_tc3,n_tc3)

densities_tc1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
densities_tc2 = [8.0,  38.0, 17.0, 1.0]
densities_tc3 = [3.0,  38.0, 22.0, 1.0]

u₀ = [0.0, 0.0]
# cons_tcr = 10 .^ rates_tcr
cons_tcr = rates_tcr

# CD4 rates from Muaz, need to get full CD4 curve
kon2 = 7.01979340915551e-06 / 6
koff2 = 0.4000
ka = 8.09e-7 / 6
kon2 = ka*0.4000
koff2 = 0.4000
cons_tc = vcat(cons_tcr, kon2, koff2)
tspan = (0,maximum(t)+1)

pinit = [0.0,0.0]

tc1 = Data([n_tc1],[t_tc1])
tc2 = Data([n_tc2],[t_tc2])
tc3 = Data([n_tc3],[t_tc3])
M_tc1 = M(trimolecular!,u₀,tspan,[densities_tc1],[cons_tc],0.0,Rodas5)
M_tc2 = M(trimolecular!,u₀,tspan,[densities_tc2],[cons_tc],0.0,Rodas5)
M_tc3 = M(trimolecular!,u₀,tspan,[densities_tc3],[cons_tc],0.0,Rodas5)

rates_tc1,hes_tc1,l_tc1 = optimization(M_tc1,tc1,bounds)
σ_tc1 = hessian2σ(hes_tc1,l_tc1,tc1.n[1],rates_tc1)
rates_tc2,hes_tc2,l_tc2 = optimization(M_tc2,tc2,bounds)
σ_tc2 = hessian2σ(hes_tc2,l_tc2,tc2.n[1],rates_tc2)
rates_tc3,hes_tc3,l_tc3 = optimization(M_tc3,tc3,(1e-6,1e3))
σ_tc3 = hessian2σ(hes_tc3,l_tc3,tc3.n[1],rates_tc3)

# rates_tc1,hes_tc1,l_tc1,boot_tc1 = optimization(trimolecular!,n_tc1,densities_tc1',
#                                     10 .^ rates_tcr,model_opt2,t_tc1)

# boot = bootstrapping(rates_tc1,n_tc1,trimolecular!,densities_tc1,10 .^rates_tcr,t_tc1)

f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc1, rates_tcr)
plot_sum(rates_tc1)
plot_sum([0., 0.])
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc2, rates_tcr)
plot_sum(rates_tc2)
plot_sum([0., 0.])
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc3, rates_tcr)
plot_sum(rates_tc3)
plot_sum([0., 0.])

display_rates([rates_tc1,rates_tc2,rates_tc3])
# plot_curve_sds(hes_tc1,l_tc1,n_tc1,rates_tc1)
# rates_tc1_sds = calculate_sds(hes_tc1,l_tc1,n_tc1,rates_tc1)

# rates_tc1,hes_tc2,l_tc2  = optimization(trimolecular!,n_tc2,densities_tc2',cons_tc,model_opt2,t_tc2)
# rates_tc1,hes_tc2,l_tc2  = optimization(trimolecular!,n_tc3,densities_tc3',cons_tc,model_opt2,t_tc3)

scatter(n_se_1[:,1],n_se_1[:,2],yerror=n_se_1[:,3],markercolor=:red)
scatter!(n_se_2[:,1],n_se_2[:,2],yerror=n_se_2[:,3],markercolor=:blue)
scatter!(n_se_3[:,1],n_se_3[:,2],yerror=n_se_3[:,3],markercolor=:cyan)

f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc1,10 .^rates_tcr)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_tc1_sds)
sol = solve(prob,Vern7(),saveat=1)
summ1 = reduce(vcat,sum(sol,dims=1))
time = sol.t
plot!(summ1,linecolor=:red)

f = (du,u,p,t) -> two_path!(du,u,p,t,densities_tc2,cons_tc)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_tc2_2paths)
sol = solve(prob,Rodas4())
summ = reduce(vcat,sum(sol,dims=1))
time = sol.t
plot!(sol.t,summ,linecolor=:blue)

f = (du,u,p,t) -> two_path!(du,u,p,t,densities_tc3,cons_tc)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_tc3_2paths)
sol = solve(prob,Rodas4())
summ = reduce(vcat,sum(sol,dims=1))
plot!(sol.t,summ,linecolor=:cyan)

function plot_sum(rates)
    prob = ODEProblem(f,u₀,tspan,rates)
    # prob = remake(prob, p=10 .^ rates)
    sol = solve(prob,Rosenbrock23())
    summ = reduce(vcat,sum(sol,dims=1))
    plot!(sol.t,summ)
end

function display_rates(rates)
    for rate in rates
        display(10 .^ rate')
    end
end

display_rates([rates_tc1_2paths,rates_tc2_2paths,rates_tc3_2paths])

# function find_stds(model,rates,data)
# end

# rates = optimization(bimolecular!,dataset1,dens',[],model_opt2,t)
# loss(p,data,model,dens,cons,t)
using ForwardDiff
using Zygote
using ReverseDiff
# loss_in = (p) -> loss(p,data,model,dens,cons,t)
loss_in = (p) -> loss_2(p,dataset1,bimolecular!,dens',[],t)
loss(p,dataset1,bimolecular!,dens',[],t)
grad_zyg = Zygote.gradient(loss_in,[-4,-0.5])
grad_for = ForwardDiff.gradient(loss_in,[-4,-0.5])
grad_rev = ReverseDiff.gradient(loss_in,[-4,-0.5])
hes_for = ForwardDiff.hessian(loss_in,[-4,-0.297])
hes_rev = ReverseDiff.hessian(loss_in,[-4,-0.297])
grad_zyg = Zygote.gradient(loss_in,rates)
grad_for = ForwardDiff.gradient(loss_in,rates)
grad_rev = ReverseDiff.gradient(loss_in,rates)
hes_for = ForwardDiff.hessian(loss_in,rates)
hes_rev = ReverseDiff.hessian(loss_in,rates)
L = loss_in(rates)
covrp = inv(hes_rev)

m = size(dataset1,2)
n = size(rates,1)
ν = m - n
MSE = L/ν
stds = sqrt.(diag(covrp).*MSE)
using LinearAlgebra
# stds = sqrt.(diag(covrp).*L)

using Measurements
rates_std = rates .± stds

prob = ODEProblem(f,u₀,tspan,10 .^ rates_std)
# prob = remake(prob, p=10 .^ rates)
sol = solve(prob,Tsit5())
plot(sol,errorbars = :ribbon, alpha=0.2)
# fig = plot!(sol,linecolor=:black)

loss_in = (p) -> loss(p,n',bimolecular!,densities_tcr',[],t)
loss(rates_tcr,n',bimolecular!,densities_tcr',[],t)

grad_zyg = Zygote.gradient(loss_in,rates_tcr)
grad_for = ForwardDiff.gradient(loss_in,rates_tcr)
grad_rev = ReverseDiff.gradient(loss_in,rates_tcr)
hes_for = ForwardDiff.hessian(loss_in,rates_tcr)
hes_rev = ReverseDiff.hessian(loss_in,rates_tcr)
hes_zyg = Zygote.hessian(loss_in,rates_tcr)

L = loss_in(rates_tcr)
covrp = inv(hes_for)
using LinearAlgebra
function plot_curve_sds(hes,l,data,rates)
    covrp = inv(hes)
    m = size(data,1)
    n = size(rates,1)
    ν = m - n
    MSE = l/ν
    stds = sqrt.(diag(covrp).*MSE)
    rates_std = rates .± stds

    prob = ODEProblem(f,u₀,tspan,10 .^ rates_std)
    # prob = remake(prob, p=10 .^ rates)
    sol = solve(prob,Vern7(),abs_tol=1e-10,rel_tol=1e-10)
    summ = reduce(vcat,sum(sol,dims=1))
    # time = sol.t
    plot!(sol,errorbars = :ribbon, alpha=0.2)
    # plot!(time,summ,errorbars = :ribbon, alpha=0.2)
    # plot!(time,summ)
end

function calculate_sds(hes,l,data,rates)
    covrp = inv(hes)
    m = size(data,1)
    n = size(rates,1)
    ν = m - n
    MSE = l/ν
    stds = sqrt.(diag(covrp).*MSE)
    rates_std = rates .± stds
end

ext: jl


text:
module RandomTweaks

using ForwardDiff
using QuadGK

include("load.jl")
export load_file_line_by_line

include("curve.jl")
export curve_length,curve_fraction

end


text:
include("blackpaper.jl")
using DelimitedFiles

# predefined hyperparameters for solution

# center location

const x_00 = 0.0
const y_00 = 0.0

# velocity field radius

const r_1 = 1.0

# center movement hyperparameters

const u_00 = 1.0 # transfer speed for center
const r_00 = 0.1 # radius of the circle where the center of the initial condition travels
const T_0  = 1.0  # time of the single rotation of the initial condition

# hyperparameters for transfer speeds

const T_1 = 0.5 # time of the single rotation around the center
const T_2 = 0.5 # Decaying factor

# hyperparameters for initial conditions

const H_0   = 1.0 # height of the cone
const r_0   = 0.5 # radius of the cone
#angles defining the slice of the cone
const phi_0 = pi / 4.0
const phi_1 = pi / 2.0

# define new fuctions from templates
# to build the solution, we need a single sample of
# initial conditions, transfer speeds, and center point movement functions

# radius-vector from central point defined above
function r(x, y)
	return r_template(x, y, x_00, y_00)
end

# x coordinate for center point
function x0_linear(t)
	return x0_linear_template(t, u_00, x_00)
end

# y coordinate for center point
function y0_linear(t)
	return x0_linear_template(t, u_00, y_00)
end

function a1_linear(r)
	return r
end

function a2_constant(t)
	return a2_constant_template(t, T_1)
end

# initial conditions in polar coordinates
function cone_polar(r, phi)
	return cone_template(r, phi, r_0, H_0)
end

function cone_slice_polar(r, phi)
	return cone_slice_template(r, phi, r_0, H_0, phi_0, phi_1)
end

# initial conditions in cartesian coordinates
function cone_cartesian(xs, ys)

	grid = [(x, y) for x in xs, y in ys]

	radius = inflate(
			(x, y) -> ((x .- x_00) ^ 2 + (y .- y_00) ^ 2) ^ 0.5,
			xs, ys
	)

	rlh = [
		r != 0.0 ?
		acos((x .- x_00) / r) :
		0.5 * pi
		for (r, (x, y)) in zip(radius, grid)
	]

	rrh = [
		r != 0.0 ?
		(2.0 * pi .- acos((x .- x_00) / r)) :
		1.5 * pi
		for  (r, (x, y)) in zip(radius, grid)
	]

	phi     = zeros((length(xs), length(ys)))
	mask    = [y .<  y_00 for (x, y) in grid]
	notmask	= [y .>= y_00 for (x, y) in grid]
	phi[mask]    = rrh[mask]
	phi[notmask] = rlh[notmask]

	periods = floor.(phi / (2.0 * pi))
	phi     = phi - periods * 2.0 * pi

	return cone_template(radius, phi, r_0, H_0)

function cone_slice_cartesian(xs, ys)

	grid = [(x, y) for x in xs, y in ys]

	radius = inflate(
			(x, y) -> ((x .- x_00) ^ 2 + (y .- y_00) ^ 2) ^ 0.5,
			xs, ys
	)

	rlh = [
		r != 0.0 ?
		acos((x .- x_00) / r) :
		0.5 * pi
		for (r, (x, y)) in zip(radius, grid)
	]

	rrh = [
		r != 0.0 ?
		(2.0 * pi .- acos((x .- x_00) / r)) :
		1.5 * pi
		for  (r, (x, y)) in zip(radius, grid)
	]

	phi     = zeros((length(xs), length(ys)))
	mask    = [y .<  y_00 for (x, y) in grid]
	notmask	= [y .>= y_00 for (x, y) in grid]
	phi[mask]    = rrh[mask]
	phi[notmask] = rlh[notmask]

	periods = floor.(phi / (2.0 * pi))
	phi     = phi - periods * 2.0 * pi

	return cone_slice_template(r, phi, r_0, H_0, phi_0, phi_1)

# transfer velocities
#function v1(x, y, t)
#	return v1_template(x, y, t, )

function solution(x, y, t)

	return solution_template(
		x, y, t,
		cone_slice,
		a1_linear,
		a2_constant,
		x0_linear,
		y0_linear
	)

end

# setup the domain

x = -1.0:0.01:1.0
y = -1.0:0.01:1.0
t = 0.0:0.1:10

res = solution(x, y, t)

# save as .csv file

writedlm("result.csv", res, ',')

max_stars_repo_name: HSE-LAMBDA/RheologyReconstruction

max_issues_repo_name: HSE-LAMBDA/RheologyReconstruction

max_forks_repo_name: HSE-LAMBDA/RheologyReconstruction


text:
#!/usr/bin/env julia
# 1.7: Roots of a polynomial
using Roots

function polynomialGenerator(a...)
  n = length(a) - 1
  poly = function (x)
    return sum([a[i + 1] * x^i for i in 0:n])
  end
  return poly
end

polynomial = polynomialGenerator(1, 3, -10)
zeroVals = find_zeros(polynomial, -10, 10)
println("Zeros of the function f(x): ", zeroVals)


text:
__precompile__()
module Mayday

using JuMP
import JuMP: getvalue
using MultiPoly
import MultiPoly: print, evaluate
import DataStructures: OrderedDict
export monomials,
       chebyshev_basis_first_kind,
       polynomial_basis,
       generator,
       generators,
       evaluate,
       deg,
       grad,
       jacobian,
       diff,
       MPoly,
       defPolynomial,
       defSoSPolynomial,
       addPolynomialEqualityConstraint,
       addSoSConstraint

generators(s::Symbol...) = MultiPoly.generators(MPoly{Float64}, s...)
generator(s::Symbol) = generators(s)[1]

include("basis.jl")

function print{T<:JuMP.GenericAffExpr}(io::IO, p::MPoly{T})
    first = true
    for (m, c) in p
        if first
            print(io, c)
            MultiPoly.printmonomial(io, m, vars(p))
            first = false
        else
            print(io, " + ($c)")
            MultiPoly.printmonomial(io, m, vars(p))
        end
    end
    if first
        print(io, zero(T))
    end
end

grad(poly::MPoly, symbols::Symbol...) = [diff(poly, s) for s in symbols]
jacobian{T}(polys::Vector{MPoly{T}}, symbols::Symbol...) = permutedims(hcat([grad(poly, symbols...) for poly in polys]...), [2 1])

function evaluate{PolyType, ArgType}(polys::Array{MPoly{PolyType}}, args::ArgType...)
    result_type = promote_type(PolyType, ArgType)
    result = Array(result_type, size(polys)...)
    for i = 1:length(polys)
        result[i] = evaluate(polys[i], args...)
    end
    return result
end

function defPolynomial{T}(model::JuMP.Model, variables::Vector{Symbol}, basis::Vector{MPoly{T}})
    @variable(model, coeffs[1:length(basis)])
    return sum(basis .* coeffs)
end

defPolynomial(model::JuMP.Model, variables::Vector{Symbol}, max_monomial_degree::Integer, basis_generator::Function=monomial) = defPolynomial(model, variables, polynomial_basis(variables, 0:max_monomial_degree, basis_generator))

function quad{Tx, TQ, Ty}(x::AbstractVector{Tx}, Q::AbstractMatrix{TQ}, y::AbstractVector{Ty})
    result = zero(promote_type(Tx, promote_type(TQ, Ty)))
    for (i, xi) in enumerate(x)
        for (j, yj) in enumerate(y)
            result += xi * Q[i,j] * yj
        end
    end
    result
end
quad(x::AbstractVector, Q::AbstractMatrix) = quad(x, Q, x)

function defSoSPolynomial{T}(model::JuMP.Model, variables::Vector{Symbol}, basis::Vector{MPoly{T}})
    @variable(model, Q[1:length(basis), 1:length(basis)], SDP)
    quad(basis, Q)
end

defSoSPolynomial(model::JuMP.Model, variables::Vector{Symbol}, max_monomial_degree::Real, basis_generator::Function=monomial) = defSoSPolynomial(model, variables, polynomial_basis(variables, 0:max_monomial_degree, basis_generator))

function addPolynomialEqualityConstraint(model::JuMP.Model, poly1::MPoly, poly2::MPoly)
    @assert poly1.vars == poly2.vars
    constrained_powers = []
    for (powers, coefficient) in poly1
        @constraint(model, poly2[powers] == coefficient)
        push!(constrained_powers, powers)
    end
    for (powers, coefficient) in poly2
        if !in(powers, constrained_powers)
            @constraint(model, poly1[powers] == coefficient)
        end
    end
end

function addSoSConstraint(model, poly, basis_generator::Function=monomial)
    max_degree = maximum([maximum(x) for x in keys(poly.terms)])
    basis = polynomial_basis(poly.vars, 0:Int(ceil(max_degree / 2)), basis_generator)
    addSoSConstraint(model, poly, basis)
end

function addSoSConstraint{T}(model, poly, basis::AbstractVector{MPoly{T}})
    sos_poly = defSoSPolynomial(model, poly.vars, basis)
    addPolynomialEqualityConstraint(model, sos_poly, poly)
end

function getvalue{T<:JuMP.GenericAffExpr}(poly::MPoly{T})
    MPoly{Float64}(OrderedDict{Vector{Int}, Float64}(zip(keys(poly.terms), map(getvalue, values(poly.terms)))), poly.vars)
end

end


text:
##  SCRIPT NAME:   mUtilH_cH.jl
##  DATE:          March 2015
#
#   This script calculates the husband's marginal utility with
#   respect to his private consumption.
#
#
##  ALEXANDROS THELOUDIS, UCL

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

function mUtilH_cH(cH, K)

  marginal_utility = alphaCH * (cH^(-gammaCH))

  return marginal_utility

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using SampledSignals, PortAudio, LibSndFile, DSP
using Logging, Observables, Unitful, Makie

include("./visualizer.jl")
include("./types.jl")

"A dictionary of different representations of the input audio signal that update when the signal updates."
struct SignalRepresentations
  map::Dict{String, Union{SignalRepresentations,Observable}}
end

# function TrackSpectrumSlice(TrackSlice)::SpectrumSlice
#   fft_periodogram = periodogram(TrackSlice, fs=TrackSlice.samplerate)
#   return SpectrumBuf(power(fft_periodogram, TrackSlice.samplerate)
# end

function SignalRepresentations(audio_sample_buffer::Observable{SampleBuf{T, 1}})::SignalRepresentations where T <: AbstractFloat

  smoothSignal(sample_buffer, width=100, σ=.5) = conv(gaussian(width, σ), sample_buffer)
  firstDerivative(sample_buffer) = conv(sample_buffer, [-1, 0, 1,])
  secondDerivative(sample_buffer) = conv(sample_buffer, [-1, 2, -1,])
  generateSpectrogram(mono_sample_buffer) = spectrogram(mono_sample_buffer; fs=mono_sample_buffer.samplerate)

  rMap = SignalRepresentations(Dict(
    "audio" => audio_sample_buffer,
    "smooth audio" => map(smoothSignal, audio_sample_buffer),
    "da/dt" => SignalRepresentations(Dict(
      "audio" => map(firstDerivative, audio_sample_buffer),
      "smooth audio" => map(firstDerivative, rMap["smooth audio"]),
    )),
    "(da/dt)²" => SignalRepresentations(Dict(
      "audio" => map(secondDerivative, audio_sample_buffer),
      "smooth audio" => map(secondDerivative, rmap["smooth audio"]),
    )),
    "spectrogram" => map(generateSpectrogram, audio_sample_buffer),
  ))

  return rMap
end

t = Template(;
             user="mostlysimilar",
             license="MIT",
             authors=["Michael Muszynski"],
             dir="~/Soft/projects",
             julia_version=v"1.3",
             plugins=[
               TravisCI(),
               Coveralls(),
               GitHubPages(),
             ],
             )

generate("QualitativeAudioFeatures", t)


text:
# This is the generic fallback
# function instrument(filename::String,
#                     expr::Expr,
#                     ::Val,
#                     lnn::LineNumberNode;
#                     prepend_tracker::Bool = true)
#     return insert_prepended_tracker(filename,
#                                     expr,
#                                     lnn;
#                                     prepend_tracker = prepend_tracker)
# end

function instrument(filename::String,
                    lnn_1::LineNumberNode,
                    lnn_2::LineNumberNode;
                    prepend_tracker::Bool = true)
    always_assert(lnn_1 == lnn_2, "lnn_1 == lnn_2")
    return lnn_1
end

function instrument(filename::String,
                    expr::Expr;
                    prepend_tracker::Bool = true)
    return instrument(filename,
                      expr,
                      Val(expr.head);
                      prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    return instrument(filename,
                      expr,
                      Val(expr.head),
                      lnn;
                      prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:block};
                    prepend_tracker::Bool = true)
    original_args = expr.args
    always_assert(original_args isa Vector, "original_args isa Vector")
    if length(original_args) == 0
        return expr
    end
    always_assert(length(original_args) >= 1, "length(original_args) >= 1")
    always_assert( original_args[1] isa LineNumberNode,
                  "original_args[1] isa LineNumberNode")
    return instrument(filename,
                      expr,
                      Val(:block),
                      original_args[1];
                      prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:block},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    arg_to_lnn = Vector{LineNumberNode}(undef, num_args)
    current_lnn = lnn
    arg_to_lnn[1] = current_lnn
    for i = 1:num_args
        if original_args[i] isa LineNumberNode
            current_lnn = original_args[i]
        end
        arg_to_lnn[i] = current_lnn
    end
    new_args = Vector{Any}(undef, num_args)
    for i = 1:num_args
        new_args[i] = instrument(filename, original_args[i], arg_to_lnn[i])
    end
    return insert_prepended_tracker(filename,
                                    Expr(:block, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:function},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    new_args = Vector{Any}(undef, num_args)
    new_args[1] = original_args[1]
    for i = 2:num_args
        new_args[i] = instrument(filename, original_args[i], lnn)
    end
    return insert_prepended_tracker(filename,
                                    Expr(:function, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:if},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    new_args = Vector{Any}(undef, num_args)
    new_args[1] = instrument(filename,
                             original_args[1],
                             lnn;
                             prepend_tracker = false)
    for i = 2:num_args
        new_args[i] = instrument(filename,
                                 original_args[i],
                                 lnn;
                                 prepend_tracker = prepend_tracker)
    end
    return insert_prepended_tracker(filename,
                                    Expr(:if, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:||},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    new_args = Vector{Any}(undef, num_args)
    for i = 1:num_args
        new_args[i] = instrument(filename, original_args[i], lnn)
    end
    return insert_prepended_tracker(filename,
                                    Expr(:||, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:&&},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    new_args = Vector{Any}(undef, num_args)
    for i = 1:num_args
        new_args[i] = instrument(filename, original_args[i], lnn)
    end
    return insert_prepended_tracker(filename,
                                    Expr(:&&, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:call},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    new_args = Vector{Any}(undef, num_args)
    new_args[1] = original_args[1]
    for i = 2:num_args
        new_args[i] = instrument(filename, original_args[i], lnn)
    end
    return insert_prepended_tracker(filename,
                                    Expr(:call, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end

function instrument(filename::String,
                    sym::Symbol,
                    lnn::LineNumberNode)
    return sym
end

function instrument(filename::String,
                    str::AbstractString,
                    lnn::LineNumberNode)
    return str
end

function instrument(filename::String,
                    n::Number,
                    lnn::LineNumberNode)
    return n
end

function instrument(filename::String,
                    expr::Expr,
                    ::Val{:return},
                    lnn::LineNumberNode;
                    prepend_tracker::Bool = true)
    original_args = expr.args
    num_args = length(original_args)
    new_args = Vector{Any}(undef, num_args)
    for i = 1:num_args
        new_args[i] = instrument(filename, original_args[i], lnn)
    end
    return insert_prepended_tracker(filename,
                                    Expr(:return, new_args...),
                                    lnn;
                                    prepend_tracker = prepend_tracker)
end


text:
export Abstract_Model2Fit
export parameter_size, accumulate_y!, eval_y, alloc_y

@doc raw"""

Abstract type, base of all model to fit.

# Interface

The interface is very simple

- [`parameter_size`](@ref)
- [`accumulate_y!`](@ref)

# Convenience functions

- [`eval_y`](@ref)

"""
abstract type Abstract_Model2Fit end

# Visit  ================
#
visit_submodel_size(model::Abstract_Model2Fit) = @assert(false,"To implement for $(typeof(model)) type!")
visit_get_submodel(model::Abstract_Model2Fit,submodel_idx::Int) = @assert(false,"To implement for $(typeof(model)) type!")
visit_get_Y(model::Abstract_Model2Fit,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = @assert(false,"To implement for $(typeof(model)) type!")
visit_get_X(model::Abstract_Model2Fit,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = @assert(false,"To implement for $(typeof(model)) type!")
visit_get_θ(model::Abstract_Model2Fit,submodel_idx::Int,Y::AbstractVector,X::AbstractVector,θ::AbstractVector) = @assert(false,"To implement for $(typeof(model)) type!")

# Public functions ================
#
@doc raw"""
```julia
parameter_size(::Abstract_Model2Fit)::Int
```

Return length of the expected parameter vector `θ`

# Also see
- [`Abstract_Model2Fit`](@ref)
"""
parameter_size(model::Abstract_Model2Fit) = @assert(false,"To implement for $(typeof(model)) type!")

@doc raw"""
```julia
accumulate_y!(::Abstract_Model2Fit,Y::AbstractVector,X::AbstractVector,θ::AbstractVector)::AbstractVector
```

Accumulate model contribution into vector `Y`.

# Also see
- [`Abstract_Model2Fit`](@ref)
- [`eval_y`](@ref)
"""
accumulate_y!(model::Abstract_Model2Fit,Y::AbstractVector,X::AbstractVector,θ::AbstractVector)  = @assert(false,"To implement for $(typeof(model)) type!")

# convenience

function alloc_y(m::Abstract_Model2Fit,X::AbstractVector,θ::AbstractVector)
    zeros(promote_type(eltype(X),eltype(θ)),length(X))
end

@doc raw"""
```julia
eval_y(m::Abstract_Model2Fit,X::AbstractVector,θ::AbstractVector)::AbstractVector
```

A convenience function that call [`accumulate_y!`](@ref) using a zero
initialized `Y` vector. This returned vector contains model values.

# Also see
- [`Abstract_Model2Fit`](@ref)
"""
eval_y(m::Abstract_Model2Fit,X::AbstractVector,θ::AbstractVector) = accumulate_y!(m,alloc_y(m,X,θ),X,θ)


text:
using AppliAR

using Test
using AppliSales
using AppliGeneralLedger
using Dates

# TEST MODEL
@testset "Orders" begin
    orders = AppliSales.process()
    @test length(orders) == 3
    @test orders[1].org.name == "Scrooge Investment Bank"
    @test orders[1].training.name == "Learn Smiling"
end

@testset "Retrieve UnpaidInvoices" begin
    orders = AppliSales.process()
    AppliAR.process(orders)
    unpaid_invoices = retrieve_unpaid_invoices()
    unpaid_invoice = unpaid_invoices[1]

    @test id(unpaid_invoice) == "A1001"
    @test currency_ratio(meta(unpaid_invoice)) == 1.0
    @test name(header(unpaid_invoice)) == "Scrooge Investment Bank"
    @test length(students(body(unpaid_invoice))) == 1
    @test price_per_student(body(unpaid_invoice)) == 1000.0
    @test first(students(body(unpaid_invoice))) == "Scrooge McDuck"
    @test vat_perc(body(unpaid_invoice)) == 0.21

    cmd = `rm test_invoicing.txt invoicenbr.txt`
    run(cmd)
end

@testset "Retrieve BankStatement from CSV" begin
    #stms = read_bank_statements("./bank.csv")
    stms = [BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)]
    @test length(stms) == 1
    @test amount(stms[1]) == 2420.0
end

@testset "JournalEntry's" begin
    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)
    stms = [stm1]

    orders = AppliSales.process()
    AppliAR.process(orders)
    invoices = retrieve_unpaid_invoices()

    AppliAR.process(invoices, stms)
    paid_invoices = retrieve_paid_invoices()

    @test length(paid_invoices) == 1
    @test id(paid_invoices[1]) == "A1002"
    @test amount(stm((paid_invoices[1]))) == 2420.0

    cmd = `rm test_invoicing.txt invoicenbr.txt`
    run(cmd)
end

@testset "process(orders)" begin
    orders = AppliSales.process()
    entries = AppliAR.process(orders)
    @test length(entries) == 3
    @test entries[1].from == 1300
    @test entries[1].to == 8000
    @test entries[1].debit == 1000.0
    @test entries[1].vat == 210.0
    @test entries[1].descr == "Learn Smiling"

    cmd = `rm test_invoicing.txt invoicenbr.txt`
    run(cmd)
end

@testset "retrieve unpaid invoices" begin
    orders = AppliSales.process()
    entries = AppliAR.process(orders)
    unpaid_invoices = retrieve_unpaid_invoices()

    @test length(unpaid_invoices) == 3
    @test id(unpaid_invoices[1]) == "A1001"

    cmd = `rm test_invoicing.txt invoicenbr.txt`
    run(cmd)
end

@testset "process(unpaid_invoices)" begin
    orders = AppliSales.process()
    AppliAR.process(orders)
    #unpaid_invoices = UnpaidInvoice[invoice for invoice in read_from_file("./test_invoicing.txt")]
    unpaid_invoices = retrieve_unpaid_invoices()

    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)
    stms = [stm1]
    entries = AppliAR.process(unpaid_invoices, stms)
    @test length(entries) == 1
    @test entries[1].from == 1150
    @test entries[1].to == 1300
    @test entries[1].debit == 2420.0

    cmd = `rm test_invoicing.txt test_invoicing_paid.txt invoicenbr.txt`
    run(cmd)
end

@testset "report" begin
    orders = AppliSales.process()
    AppliAR.process(orders)
    unpaid_invoices = retrieve_unpaid_invoices()

    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)
    stms = [stm1]

    AppliAR.process(unpaid_invoices, stms)

    #r = Reporting.aging(path)
    r = report()
    @test r[1].csm == "Scrooge Investment Bank"
    @test r[1].days == Day(0)

    cmd = `rm test_invoicing.txt test_invoicing_paid.txt invoicenbr.txt`
    run(cmd)
end

@testset "report with file names" begin
    orders = AppliSales.process()
    AppliAR.process(orders; path="./test3_invoicing.txt")
    unpaid_invoices = retrieve_unpaid_invoices(; path="./test3_invoicing.txt")

    stm1 = BankStatement(Date(2020-01-15), "Duck City Chronicals Invoice A1002", "NL93INGB", 2420.0)
    stms = [stm1]

    AppliAR.process(unpaid_invoices, stms; path="./test3_invoicing_paid.txt")

    #r = Reporting.aging(path)
    r = report(path_unpaid="./test3_invoicing.txt", path_paid="./test3_invoicing_paid.txt")
    @test r[1].csm == "Scrooge Investment Bank"
    @test r[1].days == Day(0)

    cmd = `rm test3_invoicing.txt test3_invoicing_paid.txt invoicenbr.txt`
    run(cmd)
end


text:
"""
`optimalPLD_runDNCHeuristicObjFn(P, obj_fn, total_num_products, maxIter, heuristic_seed)`

Run the divide and conquer heuristic (see Green and Krieger 1993, Belloni et al. 2008).

`P` = number of products in the product line.

`obj_fn` = arbitrary function that will compute expected revenue of a product line (= an array of product indices). Should only include actual products (not no-purchase option).

`total_num_products` = total number of candidate products from which the product line may be chosen.

`maxIter` = integer specifying number of random starting points to use.

`heuristic_seed` = integer to set the seed of the random number generator (to be used for generating the random starting points).

# References

P. E. Green and A. M. Krieger. Conjoint analysis with product-positioning applications. In J. Eliashberg and G. L. Lilien, editors, Handbooks in Operations Research and Management Science, volume 5, pages 467–515. Elsevier, 1993.

A. Belloni, R. Freund, M. Selove, and D. Simester. Optimizing product line designs: Efficient methods and comparisons. Management Science, 54(9):1544–1552, 2008.

"""
function optimalPLD_runDNCHeuristicObjFn(P, obj_fn, total_num_products, maxIter, heuristic_seed)

tic()

heuristicRevenues = zeros(maxIter,1);
heuristicSolns = zeros(maxIter,P);

srand(heuristic_seed)

for iter=1:maxIter
	currentProduct = 1;
	chosenProducts = randperm(total_num_products)[1:P]

	currentRevenue = obj_fn(chosenProducts)
	initialRevenue = currentRevenue;

	numExhaustedProducts = 0;

	while (numExhaustedProducts < P)
		# Keep going

		# Evaluate
		bestExpectedCandRevenue = 0;
		bestCandProduct = -1;
		for p2=1:total_num_products
			if (!in(p2, chosenProducts))
				cand_chosenProducts = copy(chosenProducts);
				cand_chosenProducts[currentProduct] = p2;
                expectedCandRevenue = obj_fn(cand_chosenProducts); #sum(lambdaFull .* revenues_sorted); #sum(sum(lambda .* revenues_sorted));
                if (bestExpectedCandRevenue < expectedCandRevenue)
                	bestCandProduct = p2;
                	bestExpectedCandRevenue = expectedCandRevenue;
                end
			end
		end

		if (bestExpectedCandRevenue > currentRevenue)
			numExhaustedProducts = 1;
			chosenProducts[currentProduct] = bestCandProduct;
			currentRevenue = bestExpectedCandRevenue;
		else
			numExhaustedProducts = numExhaustedProducts + 1;
		end

		currentProduct = currentProduct + 1;
		if (currentProduct > P)
			currentProduct = 1;
		end
	end

	print(currentRevenue)
	print("\t")
	println(chosenProducts)
	heuristicRevenues[iter] = currentRevenue
	heuristicSolns[iter,:] = chosenProducts

end

bestSolnInd = indmax(heuristicRevenues)
bestHeuristicSoln = heuristicSolns[bestSolnInd,:][:]
bestHeuristicRevenue = maximum(heuristicRevenues)

heuristicTime = toc()

return bestHeuristicRevenue, bestHeuristicSoln, heuristicTime

end

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
__precompile__(false)

module ASE

using Reexport

@reexport using JuLIP
@reexport using PyCall

# the functions to be implemented
import JuLIP:
      positions, set_positions!,
      cell, set_cell!,             # ✓
      pbc, set_pbc!,               # ✓
      calculator, set_calculator!, # ✓
      neighbourlist,                # ✓
      energy, forces, virial, stress,
      momenta, set_momenta!,
      masses, set_masses!,
      set_transient!,
      atomic_numbers,
      Atoms, chemical_symbols,
      get_data, has_data, set_data!,
      bulk

import Base.length, Base.deleteat!, Base.deepcopy

# from arrayconversions:
using JuLIP: mat, vecs, JVecF, JMatF,
      AbstractAtoms,
      AbstractCalculator, maxdist, SVec,
      Dofs, set_dofs!

using LinearAlgebra: det

# extra ASE functionality:
import Base: repeat         # ✓

export ASEAtoms,      # ✓
      AbstractASECalculator, ASECalculator,
      extend!, get_info, set_info!, get_array, set_array!, has_array, has_info,
      get_transient, has_transient,
      velocities, set_velocities!,
      static_neighbourlist,
      read_xyz, write_xyz

ase_build = pyimport("ase.build")
ase_atoms = pyimport("ase.atoms")
ase_io    = pyimport("ase.io")

"""
`type ASEAtoms <: AbstractAtoms`

Julia wrapper for the ASE `Atoms` class.

### Constructors:
```julia
ASEAtoms(po::PyObject)         # from a given ASE object
ASEAtoms(s::AbstractString)    # e.g., "Si2" for a Si cluster containing 2 atoms
ASEAtoms(at::Atoms)            # from a JuLIP Atoms object
```

"""
mutable struct ASEAtoms <: AbstractAtoms{Float64}
   po::PyObject       # ase.Atoms instance
   calc::Union{AbstractCalculator, Nothing}
end

Base.eltype(::ASEAtoms) = Float64

ASEAtoms(po::PyObject) = ASEAtoms(po, nothing)
pyobject(a::ASEAtoms) = a.po

# this one is needed e.g. for JuLIP conversions
ASEAtoms(s::AbstractString) = ASEAtoms(ase_atoms.Atoms(s))

ASEAtoms(syms::AbstractVector{<:String}, X::AbstractVector{<: JVec}) =
      ASEAtoms(ase_atoms.Atoms(collect(syms), collect(mat(X)')))

set_calculator!(at::ASEAtoms, calc::Union{AbstractCalculator, Nothing}) = (at.calc = calc; at)
calculator(at::ASEAtoms) = at.calc

positions(at::ASEAtoms) = at.po.get_positions()' |> vecs |> collect
set_positions!(at::ASEAtoms, X::Matrix) = (at.po.set_positions(convert(Matrix, X')); at)
set_positions!(at::ASEAtoms, X::AbstractVector{JVecF}) = set_positions!(at, Matrix(mat(X)))
# TODO: revert to "clever" set_positions!
# TODO: how to get rid of the `convert`? There is an awful lot of copying going on here!
# function set_positions!(a::ASEAtoms, p::JVecsF)
#    pold = positions(a)
#    r = maxdist(pold, p)
#    p_py = PyReverseDims(mat(p))
#    a.po.set_positions(p_py)
#    update_transient_data!(a, r)
#    return a
# end

length(at::ASEAtoms) = length(positions(at))
# TODO: could be made efficient by staying in Python
#       e.g. size(at.po.get_positions(), 1)

set_pbc!(at::ASEAtoms, val::Bool) = set_pbc!(at, (val,val,val))
set_pbc!(at::ASEAtoms, val::NTuple{3,Bool}) = (at.po.pbc = val; at)
set_pbc!(at::ASEAtoms, val::AbstractVector{Bool}) = set_pbc!(at, tuple(val...))
pbc(a::ASEAtoms) = a.po.pbc

cell(at::ASEAtoms) = _convert_cell(at.po.get_cell())
_convert_cell(C::Array) = C
_convert_cell(C::PyObject) = C.__array__()
set_cell!(at::ASEAtoms, p::Matrix) = (at.po.set_cell(p); at)
# TODO: if we want to return to transient data
# A = pinv(cell(a)) * p
# update_transient_data!(a, vecnorm(A))

# special arrays: momenta, velocities, masses, chemical_symbols
"get the momenta array"
momenta(at::ASEAtoms) = at.po.get_momenta()' |> vecs
"set the momenta array"
set_momenta!(at::ASEAtoms, p::AbstractVector{<:JVec}) = at.po.set_momenta(p |> mat |> PyReverseDims)
"get the velocities array (convert from momenta)"
velocities(at::ASEAtoms) = at.po.get_velocities()' |> vecs
"convert to momenta, then set the momenta array"
set_velocities!(at::ASEAtoms, v::AbstractVector{<:JVec}) = at.po.set_velocities(v |> mat |> PyReverseDims)
"get Vector of atom masses"
masses(at::ASEAtoms) = at.po.get_masses()
"set atom mass array as Vector{Float64}"
set_masses!(at::ASEAtoms, m::Vector{Float64}) = at.po.set_masses(m)
"return vector of chemical symbols as strings"
chemical_symbols(at::ASEAtoms) = pyobject(at).get_chemical_symbols()
"set the chemical symbols"
set_chemical_symbols!(at::ASEAtoms, s::Vector{T}) where {T <: AbstractString} =
   pyobject(at).set_chemical_symbols(s)
"return vector of atomic numbers"
atomic_numbers(at::ASEAtoms) = pyobject(at).get_atomic_numbers()

# ===================== Translations / Conversion ======================

JuLIP.Chemistry.rnn(s::AbstractString) = rnn(Symbol(s))

#
# TODO: write tests for consistency of these conversions
#       for the new test_ase code  \\  X, P, M, Z, cell, pbc
#       also the |> collect should not be needed, probably too strict
#       typing in JuLIP
Atoms(at_ase::ASE.ASEAtoms) =
   Atoms( positions(at_ase) |> collect,
          momenta(at_ase) |> collect,
          masses(at_ase) |> collect,
          atomic_numbers(at_ase) |> collect,
          JMat{Float64}(cell(at_ase)),
          pbc(at_ase),
          calculator(at_ase) )

function ASEAtoms(at::Atoms)
   # simplify by assuming there is only one species
   syms = string.(chemical_symbols(at))
   at_ase = ASEAtoms(syms, positions(at))
   set_momenta!(at_ase, momenta(at))
   set_masses!(at_ase, masses(at))
   set_cell!(at_ase, Matrix(cell(at)))
   set_pbc!(at_ase, tuple(pbc(at)...))
   set_calculator!(at_ase, calculator(at))
end

# ===================== BUILD ATOMS =================================

# @doc ase_build.bulk.__doc__ ->
bulk(s::AbstractString, args...; pbc=true, kwargs...) =
   set_pbc!(ASEAtoms(ase_build.bulk(s, args...; kwargs...)), pbc)

build(sym::Symbol, args...; kwargs...) = ase_build[sym](args...; kwargs...)

repeat(a::ASEAtoms, n::NTuple{3, Int64}) = ASEAtoms(a.po.repeat(n))
repeat(a::ASEAtoms, n::AbstractArray) = repeat(a, tuple(n...))

import Base.*
*(at::ASEAtoms, n) = repeat(at, n)
*(n::NTuple{3, Int64}, at::ASEAtoms) = repeat(at, n)
*(at::ASEAtoms, n::Integer) = repeat(at, (n,n,n))
*(n::Integer, at::ASEAtoms) = repeat(at, (n,n,n))

function deleteat!(at::ASEAtoms, n::Integer)
   at.po.__delitem__(n-1) # delete in the actual array
   # update_transient_data!(at, Inf)
   return at
end

# function deleteat!(at::ASEAtoms, nn)
#    deleteat!(at, collect(nn) - 1)
#    return at
# end

# ======================== CALCULATORS ===============================

######################################################
#    Attaching an ASE-style calculator
#    ASE-style aliases
######################################################

"""
abstract type for all calculators that interface to ASE
"""
abstract type AbstractASECalculator <: AbstractCalculator end

"""
Concrete subtype of ASECalculator for classical potentials
"""
mutable struct ASECalculator <: AbstractASECalculator
   po::PyObject
end

function set_calculator!(at::ASEAtoms, calc::AbstractASECalculator)
   at.po.set_calculator(calc.po)
   at.calc = calc
   return at
end

set_calculator!(at::ASEAtoms, po::PyObject) =
      set_calculator!(at, AbstractASECalculator(po))

forces(calc::AbstractASECalculator, at::ASEAtoms) = calc.po.get_forces(at.po)' |> vecs
energy(calc::AbstractASECalculator, at::ASEAtoms) = calc.po.get_potential_energy(at.po)

function virial(calc::AbstractASECalculator, at::ASEAtoms)
    s = calc.po.get_stress(at.po)
    vol = det(cell(at))
    if size(s) == (6,)
      # unpack stress from compressed Voigt vector form
      s11, s22, s33, s23, s13, s12 = s
      return -JMatF([s11 s12 s13;
                     s12 s22 s23;
                     s13 s23 s33]) * vol
    elseif size(s) == (3,3)
      return -JMatF(s) * vol
    else
      error("got unxpected size(stress) $(size(stress)) from ASE")
    end
end

################### extra ASE functionality  ###################
# TODO: we probably want more of this
#       and a more structured way to translate
#

"""
* `extend!(at::ASEAtoms, atadd::ASEAtoms)::ASEAtoms`
* `extend!(at::ASEAtoms, S::AbstractString, x::JVecF)::ASEAtoms`

add `atadd` atoms to `at` and returns `at`; only `at` is modified.

A short variant is
```julia
#  extend!(at, (s, x))  <<<< deprecated
extend!(at, s, x)
```
where `s` is a string, `x::JVecF` a position
"""
function extend!(at::ASEAtoms, atadd::ASEAtoms)::ASEAtoms
   at.po.extend(atadd.po)
   return at
end

function extend!(at::ASEAtoms, atnew::Tuple{S,JVecF}) where S <: AbstractString
   @warn("`extend!(at, (s,x))` is deprecated; use `extend!(at, s, x)` instead")
   extend!(at, atnew[1], atnew[2])
end

extend!(at::ASEAtoms, S::AbstractString, x::JVecF) = extend!(at, ASEAtoms(S, [x]))

"""
* `write_xyz(filename, at, mode=:write)` : write atoms object to `filename`
* `write_xyz(filehandle, at)` : write atoms object as xyz file
* `write_xyz(filename, ats::Vector{ASEAtoms}, mode=:write)` : write a time series to a file
* `write_xyz(filename, at, x::Vector{Dofs}, mode=:write)` : write a time series to a file

to append to an existing file, use `:append` or `"a"` instead of `:write`.
"""
write_xyz(filename::AbstractString, at::ASEAtoms, mode=:write) =
   mode == :write ? ase_io.write(filename, at.po) : write_xyz(filename, [at], mode)

write_xyz(filehandle::PyObject, at::ASEAtoms) = ase_io.write(filehandle, at.po, format="xyz")

# open and close files from Python (to get a python filehandle)
pyopenf(filename::AbstractString, mode::AbstractString) = py"open($(filename), $(mode))"
pyclosef(filehandle) = filehandle.close()

function write_xyz(filename::AbstractString, at::ASEAtoms, xs::AbstractVector{Dofs}, mode=:write)
   x0 = dofs(at) # save the dofs
   filehandle = pyopenf(filename, string(mode)[1:1])
   for x in xs
     write_xyz(filehandle, set_dofs!(at, x))
   end
   pyclosef(filehandle)
   set_dofs!(at, x0)   # revert to original configuration
end

function write_xyz(filename::AbstractString, ats::AbstractVector{ASEAtoms}, mode=:write)
   filehandle = pyopenf(filename, string(mode)[1:1])
   for at in ats
      write_xyz(filehandle, at)
   end
   pyclosef(filehandle)
end

read_xyz(filename::AbstractString) = ASEAtoms(ase_io.read(filename))

# ====================== TODO REWRITE =================================
#
# * deepcopy
# *

# ASEAtoms(s::AbstractString, X::JVecsF) = ASEAtoms(s, mat(X))
# ASEAtoms(s::AbstractString, X::Matrix) = ASEAtoms(ase_atoms.Atoms(s, X'))
# ASEAtoms(s::AbstractString) = ASEAtoms(ase_atoms.Atoms(s))

# function deleteat!(at::ASEAtoms, n::Integer)
#    at.po.__delitem__(n-1) # delete in the actual array
#    update_transient_data!(at, Inf)
#    return at
# end
#
# function deleteat!(at::ASEAtoms, nn)
#    deleteat!(at, collect(nn) - 1)
#    return at
# end

# "nanotube(n, m, length=1, bond=1.42, symbol=\"C\", verbose=False)"
# function nanotube(args...; kwargs...)
#    at = ASEAtoms(ase_build.nanotube(args...; kwargs...))
#    C = Matrix(cell(at))
#    if det(C) ≈ 0.0
#       C[1,1] = 1.0
#       C[2,2] = 1.0
#       @assert !(det(C) ≈ 0.0)
#    end
#    set_cell!(at, C)
#    return at
# end

include("nlist.jl")

include("models.jl")

# -------------- Calling a JuLIP Calculator on ASEAtoms --------------

# energy(V::AbstractCalculator, at::ASEAtoms) = energy(V, Atoms(at))
# forces(V::AbstractCalculator, at::ASEAtoms) = forces(V, Atoms(at))
# virial(V::AbstractCalculator, at::ASEAtoms) = virial(V, Atoms(at))
# stress(V::AbstractCalculator, at::ASEAtoms) = stress(V, Atoms(at))

# -------------- and ASE Calculator on JuLIP Atoms --------------

energy(V::AbstractASECalculator, at::Atoms) = energy(V, ASEAtoms(at))
forces(V::AbstractASECalculator, at::Atoms) = forces(V, ASEAtoms(at))
virial(V::AbstractASECalculator, at::Atoms) = virial(V, ASEAtoms(at))
# stress(V::AbstractASECalculator, at::Atoms) = stress(V, ASEAtoms(at))

end


text:
function simplify(t::Expr)
    postwalk(t) do x
        if !isa(x, Expr)
            return x
        end
        if x.head == :call
            # remove ^0 or ^1
            if x.args[1] == :(^)
                if x.args[3] == 0
                    return 1
                end
                if x.args[3] == 1
                    return x.args[2]
                end
                # replace ^1/2 with sqrt
                if x.args[3] == :(1/2)
                    return :(sqrt($(x.args[2])))
                end
            end
            # remove /1
            if x.args[1] == :(/)
                if x.args[3] == 1
                    return x.args[2]
                end
            end
        end
        return x
    end
end

"""    symbolic_symplify(ex::Expr)

a more thorouugh simplification of an expression,
removes -1/1, *(a,*(b,c)) and *(a,/(b,c)).
"""
function symbolic_symplify(ex::Expr)
    iscall(x) = false
    iscall(x, name::Symbol) = false
    iscall(x::Expr) = x.head == :call
    iscall(x::Expr, name::Symbol) = iscall(x) && x.args[1] == name

    MacroTools.postwalk(ex) do x
        if x == :param
            return :T
        end
        if x == :state
            return :u
        end
        # -1/1 => -1
        if x == :(-1/1)
            return :(-1)
        end
        # +(x) => x
        if iscall(x, :+) && length(x.args)==2
            return x.args[2]
        end
        # 1*x => x
        if iscall(x, :*) && length(x.args)==3 && x.args[2] == 1
            return x.args[3]
        end
        # *(a, *(b,c)) => *(a,b,c)
        if iscall(x, :*) && length(x.args)==3 && iscall(x.args[end], :*)
            return :(*($(x.args[2]), $(x.args[3].args[2:end]...)))
        end
        # *(a, /(b,c)) => /(*(a,b),c)
        if iscall(x, :*) && length(x.args)==3 && iscall(x.args[end], :/)
            a = x.args[2]
            b = x.args[3].args[2]
            c = x.args[3].args[3]
            num = :(*($(a), $(b)))
            # apply *(a, *(b,c)) => *(a,b,c) again
            if iscall(num, :*) && length(num.args)==3 && iscall(num.args[end], :*)
                num = :(*($(num.args[2]), $(num.args[3].args[2:end]...)))
            end
            f = :($num / $c)
            return f
        end
        return x
    end
end

stripnullterms(e) = begin
    newex = MacroTools.postwalk(e) do ex
        if typeof(ex) != Expr; return ex end
        if ex.args[1] != :(*); return ex end
        if ex.args[2] != 0; return ex end
        if ex.args[2] == 0; return 0 end
        return ex
    end
    newex = MacroTools.postwalk(newex) do ex
        if typeof(ex) != Expr; return ex end
        if ex.args[1] != :(+); return ex end
        if ex.args[2] == 0
            return ex.args[3]
        end
        if ex.args[3] == 0
            return ex.args[2]
        end
        return ex
    end
    return newex
end

function funcbody(ex::Equation, ctx=:state)
  return ex.lhs.op.name => funcbody(ex.rhs, ctx)
end

function funcbody(ex::Operation, ctx=:state)
  args = Symbol[]
  body = postwalk(convert(Expr, ex)) do x
    # @show x, typeof(x);
    if typeof(x) == Expr && x.head == :call
      if length(x.args) == 1
        var = x.args[1]
        push!(args, var)
        return :($ctx.$var)
      end
    end
    return x
  end
  return body, Set(args)
end

funckit(fname, args, body) = quote $fname($(collect(args)...)) = $body end
funckit(fname::Symbol, arg::Symbol, body) = quote $fname($arg) = $body end

""" funckit(p::Problem, ctx=:state)

Compile petri net problem to native Julia expressions for faster solving
"""
function funckit(p::Problem, ctx=:state)
  # @show "Λs"
  λf = map(p.m.Λ) do λ
    body, args = funcbody(λ, ctx)
    fname = gensym("λ")
    q = funckit(fname, ctx, body)
    return q
  end
  # @show "Δs"
  δf = map(p.m.Δ) do δ
    q = quote end
    map(δ) do f
      vname, vfunc = funcbody(f, ctx)
      body, args = vfunc
      qi = :(state.$vname = $body)
      push!(q.args, qi)
    end
    sym = gensym("δ")
    :($sym(state) = $(q) )
  end

  # @show "Φs"
  ϕf = map(p.m.Φ) do ϕ
    body, args = funcbody(ϕ, ctx)
    fname = gensym("ϕ")
    q = funckit(fname, ctx, body)
  end
  return Model(p.m.S, δf, λf, ϕf)
end

quotesplat(x) = quote $(x...) end


text:
using FromFile
using SymbolicUtils
@from "Core.jl" import CONST_TYPE, Node, Options
@from "Utils.jl" import isgood, isbad, @return_on_false

const SYMBOLIC_UTILS_TYPES = Union{<:Number,SymbolicUtils.Symbolic{<:Number}}

const SUPPORTED_OPS = (cos, sin, exp, cot, tan, csc, sec, +, -, *, /)

isgood(x::SymbolicUtils.Symbolic) = SymbolicUtils.istree(x) ? all(isgood.([SymbolicUtils.operation(x);SymbolicUtils.arguments(x)])) : true
subs_bad(x) = isgood(x) ? x : Inf

function parse_tree_to_eqs(tree::Node, options::Options, index_functions::Bool=false, evaluate_functions::Bool=false)
    if tree.degree == 0
        # Return constant if needed
        tree.constant && return subs_bad(tree.val)
        return SymbolicUtils.Sym{Number}(Symbol("x$(tree.feature)"))
    end
    # Collect the next children
    children = tree.degree >= 2 ? (tree.l, tree.r) : (tree.l,)
    # Get the operation
    op = tree.degree > 1 ? options.binops[tree.op] : options.unaops[tree.op]
    # Create an N tuple of Numbers for each argument
    dtypes = map(x->Number, 1:tree.degree)
    #
    if !(op ∈ SUPPORTED_OPS) && index_functions
        op = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple{dtypes...}, Number}}(Symbol(op))
    end

    return subs_bad(op(map(x->parse_tree_to_eqs(x, options, index_functions, evaluate_functions), children)...))
end

## Convert symbolic function back
convert_to_function(x, args...) = x

function convert_to_function(x::SymbolicUtils.Sym{SymbolicUtils.FnType{T, Number}}, options::Options) where {T <: Tuple}
    degree = length(T.types)
    if degree == 1
        ind = findoperation(x.name, options.unaops)
        return options.unaops[ind]
    elseif degree == 2
        ind = findoperation(x.name, options.binops)
        return options.binops[ind]
    else
        throw(AssertionError("Function $(String(x.name)) has degree > 2 !"))
    end
end

# Split equation
function split_eq(op, args, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
    !(op ∈ (sum, prod, +, *)) && throw(error("Unsupported operation $op in expression!"))
    if Symbol(op) == Symbol(sum)
        ind = findoperation(+, options.binops)
    elseif Symbol(op) === Symbol(prod)
        ind = findoperation(*, options.binops)
    else
        ind = findoperation(op, options.binops)
    end
    return Node(ind, convert(Node, args[1], options; varMap=varMap), convert(Node, op(args[2:end]...), options; varMap=varMap))
end

function findoperation(op, ops)
    for (i,oi) in enumerate(ops)
        Symbol(oi) == Symbol(op) && return i
    end
    throw(error("Operation $op in expression not found in operationrs $ops!"))
end

function Base.convert(::typeof(Node), x::Number, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
    return Node(CONST_TYPE(x))
end

function Base.convert(::typeof(Node), x::Symbol, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
    varMap === nothing && return Node(String(x))
    return Node(String(x), varMap)
end

function Base.convert(::typeof(Node), x::SymbolicUtils.Symbolic, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
    if !SymbolicUtils.istree(x)
        varMap === nothing && return Node(String(x.name))
        return Node(String(x.name), varMap)
    end

    op = convert_to_function(SymbolicUtils.operation(x))
    args = SymbolicUtils.arguments(x)

    length(args) > 2 && return split_eq(op, args, options; varMap=varMap)
    ind = length(args) == 2 ? findoperation(op, options.binops) : findoperation(op, options.unaops)

    return Node(ind, map(x->convert(Node, x, options, varMap=varMap), args)...)
end

"""
    node_to_symbolic(tree::Node, options::Options;
                varMap::Union{Array{String, 1}, Nothing}=nothing,
                evaluate_functions::Bool=false,
                index_functions::Bool=false)

The interface to SymbolicUtils.jl. Passing a tree to this function
will generate a symbolic equation in SymbolicUtils.jl format.

## Arguments

- `tree::Node`: The equation to convert.
- `options::Options`: Options, which contains the operators used in the equation.
- `varMap::Union{Array{String, 1}, Nothing}=nothing`: What variable names to use for
    each feature. Default is [x1, x2, x3, ...].
- `evaluate_functions::Bool=false`: Whether to evaluate the operators, or
    leave them as symbolic.
- `index_functions::Bool=false`: Whether to generate special names for the
    operators, which then allows one to convert back to a `Node` format
    using `symbolic_to_node`.
"""
function node_to_symbolic(tree::Node, options::Options;
                     varMap::Union{Array{String, 1}, Nothing}=nothing,
                     evaluate_functions::Bool=false,
                     index_functions::Bool=true
                     )
    expr = subs_bad(parse_tree_to_eqs(tree, options, index_functions, evaluate_functions))
    # Check for NaN and Inf
    @assert isgood(expr) "The recovered equation contains NaN or Inf."
    # Return if no varMap is given
    varMap === nothing && return expr
    # Create a substitution tuple
    subs = Dict(
        [SymbolicUtils.Sym{Number}(Symbol("x$(i)")) => SymbolicUtils.Sym{Number}(Symbol(varMap[i])) for i in 1:length(varMap)]...
    )
    return substitute(expr, subs)
end

function symbolic_to_node(eqn::T, options::Options;
                          varMap::Union{Array{String, 1}, Nothing}=nothing)::Node where {T<:SymbolicUtils.Symbolic}

    convert(Node, eqn, options; varMap=varMap)
end


text:
import DSGE: n_model_states, n_backward_looking_states

path = dirname(@__FILE__)

### Model
m = KrusellSmith()

### Steady State
@load "$path/reference/steady_state.jld2" saved_l saved_c saved_mu saved_K saved_KF

@testset "Steady State" begin
    @test saved_l ≈ m[:lstar].value
    @test saved_c ≈ m[:cstar].value
    @test saved_mu ≈ m[:μstar].value
    @test saved_K ≈ m[:Kstar].value
    @test @test_matrix_approx_eq saved_KF m[:KFstar].value
end

### Jacobian
m.testing = true
JJ       = DSGE.jacobian(m)
@load "$path/reference/jacobian.jld2" saved_JJ nw

# we will always order things XP YP X Y
# convention is that capital letters generally refer to indices
# XP
LMP   = 1     :nw
LELLP = nw+1  :2*nw
KKP   = 2*nw+1

# YP
ZP    = 2*nw+2
MP    = 2*nw+3:3*nw+2
ELLP  = 3*nw+3:4*nw+2

# X
LM    = 4*nw+3:5*nw+2
LELL  = 5*nw+3:6*nw+2
KK    = 6*nw+3

# Y
Z     = 6*nw+4
M     = 6*nw+5:7*nw+4
ELL   = 7*nw+5:8*nw+4

# create objects needed for solve.jl
F1 = 1     :nw
F2 = nw+1  :2*nw
F3 = 2*nw+1:3*nw
F4 = 3*nw+1:4*nw
F5 = 4*nw+1:4*nw+1
F6 = 4*nw+2:4*nw+2

@testset "Jacobian" begin
    @testset "Euler Equation" begin
        @test @test_matrix_approx_eq saved_JJ[F1, KKP]  JJ[F1, KKP]
        @test @test_matrix_approx_eq saved_JJ[F1, ZP]   JJ[F1, ZP]
        @test @test_matrix_approx_eq saved_JJ[F1, ELLP] JJ[F1, ELLP]
        @test @test_matrix_approx_eq saved_JJ[F1, ELL]  JJ[F1, ELL]
    end

    @testset "Kolmogorov Forward Equation" begin
        @test @test_matrix_approx_eq saved_JJ[F2, LM]   JJ[F2, LM]
        @test @test_matrix_approx_eq saved_JJ[F2, LELL] JJ[F2, LELL]
        @test @test_matrix_approx_eq saved_JJ[F2, KK]   JJ[F2, KK]
        @test @test_matrix_approx_eq saved_JJ[F2, Z]    JJ[F2, Z]
        @test @test_matrix_approx_eq saved_JJ[F2, M]    JJ[F2, M]
    end

    @testset "LM(t+1) = M(t)" begin
        @test @test_matrix_approx_eq saved_JJ[F3, LMP]   JJ[F3, LMP]
        @test @test_matrix_approx_eq saved_JJ[F3, M]     JJ[F3, M]
    end

    @testset "LELL(t+1) = ELL(t)" begin
        @test @test_matrix_approx_eq saved_JJ[F4, LELLP]   JJ[F4, LELLP]
        @test @test_matrix_approx_eq saved_JJ[F4, ELL]     JJ[F4, ELL]
    end

    @testset "LOM K" begin
        @test @test_matrix_approx_eq saved_JJ[F5, ELL]   JJ[F5, ELL]
        @test @test_matrix_approx_eq saved_JJ[F5, M]     JJ[F5, M]
        @test @test_matrix_approx_eq saved_JJ[F5, KKP]   JJ[F5, KKP]
    end

    @testset "TFP" begin
        @test @test_matrix_approx_eq saved_JJ[F6, ZP]    JJ[F6, ZP]
        @test @test_matrix_approx_eq saved_JJ[F6, Z]     JJ[F6, Z]
    end
end

### Klein
m.testing = false
TTT_jump, TTT_state = klein(m)
@load "$path/reference/klein.jld2" gx hx

@testset "Solve" begin
    @test @test_matrix_approx_eq gx TTT_jump
    @test @test_matrix_approx_eq hx TTT_state
end

### Filtering

# Setup
N    = 200
m    = KrusellSmith()
endo = m.endogenous_states

# Model Solution/Transition Equation
TTT, RRR, CCC = solve(m)

# Measurement Equation
meas = measurement(m, TTT, RRR, CCC)
ZZ  = meas[:ZZ]
DD  = meas[:DD]
EE  = fill(0.1, (1,1))
QQ  = meas[:QQ]

# Generate measurement errors and shocks
srand(42)
u_t     = EE*randn(n_observables(m), N)
ε_t     = QQ*randn(n_shocks_exogenous(m), N)

####################################
# Simulate states forward N periods
####################################
s_t = zeros(n_model_states(m), N)
for t in 1:N
    # Initialize with shock around steady state
    if t == 1
        s_t[:, t] = RRR*ε_t[:, t]
    else
        s_t[:, t] = TTT*s_t[:, t-1] + RRR*ε_t[:, t]
    end
end

# Simulation in terms of grid values
simulated_log_gdp = vec(s_t[endo[:z′_t], :] + (m[:α]/m[:Kstar])*s_t[endo[:K′_t], :])
meas_log_gdp      = simulated_log_gdp + u_t[1, :]

##########
# Filter!
##########

# Measurement equation
data = ZZ*s_t + u_t   # Should be the same as meas_log_gdp

# Initialize with 0 initial condition
s_0     = zeros(n_model_states(m))
P_0     = RRR*QQ*RRR'

# Testing the simulation and filter setup
file = jldopen("$path/reference/simulate_and_filter.jld2", "r")
saved_u_t  = read(file, "u_t")
saved_ε_t  = read(file, "eps_t")
saved_s_t  = read(file, "s_t")
saved_s_jump_t = read(file, "s_jump_t")
saved_simulated_log_gdp = read(file, "simulated_log_gdp")
saved_meas_log_gdp = read(file, "meas_log_gdp")
saved_data = read(file, "data")
saved_s_0  = read(file, "s_0")
saved_P_0  = read(file, "P_0")
close(file)

@testset "Simulate and Filter" begin
    @test saved_u_t ≈ u_t
    @test saved_ε_t ≈ ε_t
    @test saved_s_t ≈ s_t[1:n_backward_looking_states(m), :]
    @test saved_s_jump_t ≈ s_t[n_backward_looking_states(m)+1:end, :]
    @test saved_simulated_log_gdp ≈ simulated_log_gdp
    @test saved_meas_log_gdp ≈ meas_log_gdp
    @test saved_data ≈ data
    @test saved_s_0 ≈ s_0[1:n_backward_looking_states(m)]
    @test saved_P_0 ≈ P_0[1:n_backward_looking_states(m), 1:n_backward_looking_states(m)]
end

# Testing the filter inputs against the outputs
EE_zero = zeros(1, 1)
# reference_output = load("$path/reference/filter_inputs_output.jld2")
@load "$path/reference/filter_inputs_output.jld2" y loglik
test = kalman_filter(y, TTT, RRR, CCC, QQ, ZZ,
                     DD, EE_zero, s_0, P_0)

@testset "Check Filter Output" begin
    @test test[1] ≈ loglik
end


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"
using UUIDs
using JSON

function send(n)
    # 1. Create Connection
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        println("Created Connection")
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            println("Created Channel")
            # 3. Create queue for response
            success, queue_name, message_count, consumer_count = queue_declare(chan, "", exclusive=true)
            corr_id = string(uuid4())

            # 4. Receive the response
            on_receive = (msg) -> begin
                try
                    response = String(msg.data)
                    msg_corr_id = convert(String, msg.properties[:correlation_id])
                    println(msg_corr_id)
                    if msg_corr_id == corr_id
                        println("Received RPC response: $response")
                    end
                    basic_ack(chan, msg.delivery_tag)
                catch err
                    println(err)
                end
            end

            success, consumer_tag = basic_consume(chan, queue_name, on_receive)

            # 5. Send RPC Request in JSON encoding since integer is disallowed
            data = convert(Vector{UInt8}, codeunits(JSON.json(n)))
            msg = Message(data, content_type="text/plain", delivery_mode=PERSISTENT, correlation_id=corr_id,
                          reply_to=queue_name)
            routing_key = "rpc_queue"
            basic_publish(chan, msg; exchange="", routing_key=routing_key)
            println("RPC Request sent: $n, routing_key: $queue_name")

            # Listen for ever
            while true
                sleep(1)
            end

        end
    end
end

# Don't exit on Ctrl-C
Base.exit_on_sigint(false)
try
    send(30)
catch ex
    if ex isa InterruptException
        println("Interrupted")
    else
        println("Exception: $ex")
    end
end

max_forks_repo_name: tomashq/rabbitmq-tutorials


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    Plane{T,N} <: ParametricSurface{T,N}

Infinite planar surface where the positive normal side is outside the surface.

By default this will not create any geometry for visualization, the optional `vishalfsizeu` and `vishalfsizev` arguments can be used to draw the plane as a rectangle for visualization **note that this does not fully represent the surface**.
In this case, the rotation of the rectangle around the normal to the plane is defined by `visvec` - `surfacenormal×visvec` is taken as the vector along the u axis.

```julia
Plane(surfacenormal::SVector{N,T}, pointonplane::SVector{N,T}; interface::NullOrFresnel{T} = nullinterface(T), vishalfsizeu::T = 0.0, vishalfsizev::T = 0.0, visvec::SVector{N,T} = [0.0, 1.0, 0.0])
Plane(nx::T, ny::T, nz::T, x::T, y::T, z::T; interface::NullOrFresnel{T} = nullinterface(T), vishalfsizeu::T = 0.0, vishalfsizev::T = 0.0, visvec::SVector{N,T} = [0.0, 1.0, 0.0])
```
"""
struct Plane{T, N} <: ParametricSurface{T, N}
    normal::SVector{N, T}
    d::T
    pointonplane::SVector{N, T}
    interface::NullOrFresnel{T}
    # below only for visualization purposes
    vishalfsizeu::T
    vishalfsizev::T
    visuvec::SVector{N, T}
    visvvec::SVector{N, T}

    function Plane(
        surfacenormal::AbstractArray{T, 1},
        pointonplane::AbstractArray{T, 1};
        interface::NullOrFresnel{T} = NullInterface(T),
        vishalfsizeu::T = zero(T),
        vishalfsizev::T = zero(T),
        visvec::AbstractArray{T, 1} = [0.0, 1.0, 0.0],
    ) where {T}
        @assert length(surfacenormal) == length(pointonplane) == length(visvec)
        N = length(surfacenormal)
        return Plane(
            SVector{N, T}(surfacenormal),
            SVector{N, T}(pointonplane),
            interface = interface,
            vishalfsizeu = vishalfsizeu,
            vishalfsizev = vishalfsizev,
            visvec = SVector{N, T}(visvec),
        )
    end

    function Plane(
        surfacenormal::SVector{N, T},
        pointonplane::SVector{N, T};
        interface::NullOrFresnel{T} = NullInterface(T),
        vishalfsizeu::T = zero(T),
        vishalfsizev::T = zero(T),
        visvec::SVector{N, T} = SVector{3, T}(0.0, 1.0, 0.0),
    ) where {T <: Real, N}
        norml = normalize(surfacenormal)
        d = dot(norml, pointonplane)
        if abs(dot(visvec, norml)) == one(T)
            visvec = SVector{3, T}(1.0, 0.0, 0.0)
        end
        uvec = normalize(cross(normalize(visvec), norml))
        vvec = normalize(cross(norml, uvec))
        return new{T, N}(
            norml,
            d,
            pointonplane,
            interface,
            vishalfsizeu,
            vishalfsizev,
            uvec,
            vvec,
        )
    end

    function Plane(
        nx::T,
        ny::T,
        nz::T,
        x::T,
        y::T,
        z::T;
        interface::NullOrFresnel{T} = NullInterface(T),
        vishalfsizeu::T = zero(T),
        vishalfsizev::T = zero(T),
        visvec = SVector{3, T}(0.0, 1.0, 0.0),
    ) where {T <: Real}
        return Plane(
            SVector{3, T}(nx, ny, nz),
            SVector{3, T}(x, y, z),
            interface = interface,
            vishalfsizeu = vishalfsizeu,
            vishalfsizev = vishalfsizev,
            visvec = visvec,
        )
    end
end

function Base.show(io::IO, a::Plane{T}) where {T <: Real}
    return print(
        io,
        "Plane{$T}($(a.pointonplane), $(normal(a)), $(interface(a)))",
    )
end

interface(a::Plane{T, N}) where {T <: Real, N} = a.interface
normal(pln::Plane{T, N}) where {T <: Real, N} = pln.normal

function inside(pln::Plane{T, 3}, p::SVector{3, T}) where {T <: Real}
    return dot(normal(pln), p) - pln.d < zero(T)
end

function onsurface(pln::Plane{T, 3}, p::SVector{3, T}) where {T <: Real}
    return samepoint(dot(normal(pln), p), pln.d)
end

function distancefromplane(
    p::Plane{T, N},
    point::SVector{N, T},
) where {N, T <: Real}
    return dot(normal(p), point) - p.d
end

function uvrange(::Type{Plane{T, N}}) where {T <: Real, N}
    return ((-one(T), one(T)), (-one(T), one(T)))
end
function point(p::Plane{T}, u::T, v::T) where {T <: Real}
    return p.pointonplane +
           p.vishalfsizeu * u * p.visuvec +
           p.vishalfsizev * v * p.visvvec
end
normal(p::Plane{T}, ::T, ::T) where {T <: Real} = normal(p)
pointonplane(p::Plane{T}) where {T <: Real} = p.pointonplane

function surfaceintersection(
    pln::Plane{T, N},
    r::AbstractRay{T, N},
) where {T <: Real, N}
    n̂ = normal(pln)
    d = direction(r)
    o = origin(r)
    nd = dot(n̂, d)
    if samepoint(nd, zero(T))
        # ray and plane are parallel
        if inside(pln, o) || onsurface(pln, o)
            return rayorigininterval(Infinity(T))
        else
            # no intersection only if the ray is strictly outside of the palne
            return EmptyInterval(T)
        end
    end
    t = (pln.d - dot(n̂, o)) / nd
    if t < zero(T)
        if inside(pln, o)
            # if the ray starts 'inside' the surface then we want to return a ray so intersection works
            return rayorigininterval(Infinity(T))
        else
            return EmptyInterval(T) # no ray plane intersection
        end
    end
    temp = Intersection(t, point(r, t), n̂, zero(T), zero(T), interface(pln))
    if nd < zero(T)
        return positivehalfspace(temp)
    else
        return rayorigininterval(temp)
    end
end

function BoundingBox(pln::Plane{T, 3}) where {T <: Real}
    # TODO! this is far from ideal, we should try and do something better for intersection with non-axis-algined planes
    # valid for axis aligned planes, otherwise we have to assume an infinite bounding box
    if normal(pln) === SVector{3, T}(0, 0, 1)
        return BoundingBox(
            typemin(T),
            typemax(T),
            typemin(T),
            typemax(T),
            typemin(T),
            pln.pointonplane[3],
        )
    elseif normal(pln) === SVector{3, T}(0, 0, -1)
        return BoundingBox(
            typemin(T),
            typemax(T),
            typemin(T),
            typemax(T),
            pln.pointonplane[3],
            typemax(T),
        )
    elseif normal(pln) === SVector{3, T}(0, 1, 0)
        return BoundingBox(
            typemin(T),
            typemax(T),
            typemin(T),
            pln.pointonplane[2],
            typemin(T),
            typemax(T),
        )
    elseif normal(pln) === SVector{3, T}(0, -1, 0)
        return BoundingBox(
            typemin(T),
            typemax(T),
            pln.pointonplane[2],
            typemax(T),
            typemin(T),
            typemax(T),
        )
    elseif normal(pln) === SVector{3, T}(1, 0, 0)
        return BoundingBox(
            typemin(T),
            pln.pointonplane[1],
            typemin(T),
            typemax(T),
            typemin(T),
            typemax(T),
        )
    elseif normal(pln) === SVector{3, T}(-1, 0, 0)
        return BoundingBox(
            pln.pointonplane[1],
            typemax(T),
            typemin(T),
            typemax(T),
            typemin(T),
            typemax(T),
        )
    else
        return BoundingBox(
            typemin(T),
            typemax(T),
            typemin(T),
            typemax(T),
            typemin(T),
            typemax(T),
        )
    end
end


text:
using JGE:main

JGE.main()


text:
import Glucose
using Glucose: optim_all_vars_exo
using SuParameters
using RunTools
using Omega: normal, uniform, ciid

function runparams()
  φ = Params()
  φ.simname = "train"
  φ.runname = ciid(randrunname)
  φ.tags = ["odeoptim", "firsttry"]
  φ.logdir = ciid(ω -> logdir(runname = φ.runname(ω), tags = φ.tags))
  φ.runfile = @__FILE__
  φ.gitinfo = current_commit(@__FILE__)
  φ
end

"Optimization Parameters"
function optparams()
  Params((binsize = uniform([10]), # 1, 5, 10
          batchsize = uniform([10]), # 6, 8, 10, 12, 15
          maxiters = uniform([300]), # [150:150:1050...]
          lr = uniform([0.01]), # 0.001, 0.005, 0.01, 0.05, 0.1
          hypoid = uniform([1,2,3]))) # 1,2,3
end

function train(params)
  samples, losses, best_pred = optim_all_vars_no_exo(params.binsize,
                                                  params.batchsize,
                                                  params.maxiters,
                                                  params.lr,
                                                  hypo_id = params.hypoid,
                                                  log_dir = params.logdir)
end

function allparams()
  φ = Params()
  merge(φ, runparams(), optparams()) # netparams()
end

"Run with julia -L params.jl -E 'hyper(;)' -- --queue"
function hyper(; params = Params(), n = 100)
  params_ = allparams()
  paramsamples = rand(params_, n)
  display.(paramsamples)
  control(train, paramsamples)
end

function testrun()
  params_ = allparams()
  train(rand(params_))
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
in_channel = 10
out_channel = 5
num_V = 6
num_E = 7
adj = [0. 1. 0. 0. 0. 0.;
       1. 0. 0. 1. 1. 1.;
       0. 0. 0. 0. 0. 1.;
       0. 1. 0. 0. 1. 0.;
       0. 1. 0. 1. 0. 1.;
       0. 1. 1. 0. 1. 0.]

struct NewGNLayer <: GraphNet
end

V = rand(in_channel, num_V)
E = rand(in_channel, 2*num_E)
u = rand(in_channel)

@testset "gn" begin
    l = NewGNLayer()

    @testset "without aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg)

        fg = FeaturedGraph(adj, V)
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (0, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    @testset "with neighbor aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg, :add)

        fg = FeaturedGraph(adj, V, E, zeros(0))
        l = NewGNLayer()
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (in_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.update_edge(l::NewGNLayer, e, vi, vj, u) = rand(out_channel)
    @testset "update edge with neighbor aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg, :add)

        fg = FeaturedGraph(adj, V, E, zeros(0))
        l = NewGNLayer()
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (in_channel, num_V)
        @test size(edge_feature(fg_)) == (out_channel, 2*num_E)
        @test size(global_feature(fg_)) == (0,)
    end

    GeometricFlux.update_vertex(l::NewGNLayer, ē, vi, u) = rand(out_channel)
    @testset "update edge/vertex with all aggregation" begin
        (l::NewGNLayer)(fg) = propagate(l, fg, :add, :add, :add)

        fg = FeaturedGraph(adj, V, E, u)
        l = NewGNLayer()
        fg_ = l(fg)

        @test graph(fg_) === adj
        @test size(node_feature(fg_)) == (out_channel, num_V)
        @test size(edge_feature(fg_)) == (out_channel, 2*num_E)
        @test size(global_feature(fg_)) == (in_channel,)
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _offer(filter::_ViolationFilter, v::_Violation)::Nothing
    if v.monitored_line.offset ∉ keys(filter.queues)
        filter.queues[v.monitored_line.offset] =
            PriorityQueue{_Violation,Float64}()
    end
    q::PriorityQueue{_Violation,Float64} =
        filter.queues[v.monitored_line.offset]
    if length(q) < filter.max_per_line
        enqueue!(q, v => v.amount)
    else
        if v.amount > peek(q)[1].amount
            dequeue!(q)
            enqueue!(q, v => v.amount)
        end
    end
    return nothing
end

function _query(filter::_ViolationFilter)::Array{_Violation,1}
    violations = Array{_Violation,1}()
    time_queue = PriorityQueue{_Violation,Float64}()
    for l in keys(filter.queues)
        line_queue = filter.queues[l]
        while length(line_queue) > 0
            v = dequeue!(line_queue)
            if length(time_queue) < filter.max_total
                enqueue!(time_queue, v => v.amount)
            else
                if v.amount > peek(time_queue)[1].amount
                    dequeue!(time_queue)
                    enqueue!(time_queue, v => v.amount)
                end
            end
        end
    end
    while length(time_queue) > 0
        violations = [violations; dequeue!(time_queue)]
    end
    return violations
end


text:
using Test
using ITensors
using Random

Random.seed!(1234)

@testset "Test contractions with scalar-like ITensors" begin
  i = Index(2, "i")
  j = Index(2, "j")
  k = Index(2, "k")
  α = Index(1, "α")

  is = (i, j, k)

  A = randomITensor(is..., dag(α))
  B = ITensor(2, α, α', α'')

  C = A * B
  @test C ≈ B[1, 1, 1] * A * ITensor(1, inds(B))

  C = emptyITensor(is..., α', α'')
  C .= A .* B
  @test C ≈ B[1, 1, 1] * A * ITensor(1, inds(B))

  C = emptyITensor(shuffle([(is..., α', α'')...])...)
  C .= A .* B
  @test C ≈ B[1, 1, 1] * A * ITensor(1, inds(B))
end

@testset "NaN in-place contraction bug regression test" begin
  BlasFloats = (Float32, Float64, ComplexF32, ComplexF64)
  @testset "Scalar contract, no permutation" for ElA in BlasFloats, ElB in BlasFloats
    i = Index(2, "i")
    j = Index(3, "j")
    α = Index(1, "α")
    A = randomITensor(ElA, i, j, α')
    B = randomITensor(ElB, dag(α)', α)
    R = ITensor(promote_type(ElA, ElB), i, j, α)

    R .= NaN
    @test any(isnan, R)

    R .= A .* B
    @test !any(isnan, R)
    @test array(R) ≈ array(A) * array(B)[]

    R .= NaN
    @test any(isnan, R)

    R .= B .* A
    @test !any(isnan, R)
    @test array(R) ≈ array(A) * array(B)[]
  end

  @testset "Scalar contraction, permutation" for ElA in BlasFloats, ElB in BlasFloats
    i = Index(2, "i")
    j = Index(3, "j")
    α = Index(1, "α")
    A = randomITensor(ElA, i, j, α')
    B = randomITensor(ElB, dag(α)', α)
    R = ITensor(promote_type(ElA, ElB), j, i, α)

    R .= NaN
    @test any(isnan, R)

    R .= A .* B
    @test !any(isnan, R)
    @test array(R) ≈ permutedims(array(A), (2, 1, 3)) * array(B)[]

    R .= NaN
    @test any(isnan, R)

    R .= B .* A
    @test !any(isnan, R)
    @test array(R) ≈ permutedims(array(A), (2, 1, 3)) * array(B)[]
  end

  @testset "General contraction, permutation" for ElA in BlasFloats, ElB in BlasFloats
    i = Index(2, "i")
    j = Index(3, "j")
    α = Index(2, "α")
    A = randomITensor(ElA, i, j, α')
    B = randomITensor(ElB, dag(α)', α)
    R = ITensor(promote_type(ElA, ElB), j, i, α)

    R .= NaN
    @test any(isnan, R)

    R .= A .* B
    @test !any(isnan, R)
    @test reshape(array(R), 6, 2) ≈
      reshape(permutedims(array(A), (2, 1, 3)), 6, 2) * array(B)

    R .= NaN
    @test any(isnan, R)

    R .= B .* A
    @test !any(isnan, R)
    @test reshape(array(R), 6, 2) ≈
      reshape(permutedims(array(A), (2, 1, 3)), 6, 2) * array(B)
  end
end