text:
"""
	HighShelf(f, dBGain, BWoct, fs=96e3)

Return Biquad of sampling frequency fs as DSP.SecondOrderSections that resembles Lake High Shelf EQ.

The implementation is according to:

https://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html (12.6.2019)

where the analog Q formula is utilized instead of the prewarped version to
fit the output of the Lake processor.
"""
function HighShelf(f, dBGain, BWoct, fs=96e3)
	#intermediate variables
	A = 10^(dBGain/40)
	w0 = 2*pi*f/fs
	#alpha = sin(w0)*sinh(log(2)/2*BWoct*w0/sin(w0))
	alpha = sin(w0)*sinh(log(2)/2*BWoct)
	#coefficients
	b0 = A*((A+1)+(A-1)*cos(w0)+2*sqrt(A)*alpha)
	b1 = -2*A*((A-1)+(A+1)*cos(w0))
	b2 = A*((A+1)+(A-1)*cos(w0)-2*sqrt(A)*alpha)
	a0 = (A+1)-(A-1)*cos(w0)+2*sqrt(A)*alpha
	a1 = 2*((A-1)-(A+1)*cos(w0))
	a2 = (A+1)-(A-1)*cos(w0)-2*sqrt(A)*alpha
	#normalize
	result = [b0 b1 b2 a1 a2]./a0
	return SecondOrderSections(Biquad(result...))
end

"""
	LowShelf(f, dBGain, BWoct, fs=96e3)

Return Biquad of sampling frequency fs as DSP.SecondOrderSections that resembles Lake Low Shelf EQ.

The implementation is according to:

https://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html (12.6.2019)

where the analog Q formula is utilized instead of the prewarped version to
fit the output of the Lake processor.
"""
function LowShelf(f, dBGain, BWoct, fs=96e3)
	#intermediate variables
	A = 10^(dBGain/40)
	w0 = 2*pi*f/fs
	#alpha = sin(w0)*sinh(log(2)/2*BWoct*w0/sin(w0))
	alpha = sin(w0)*sinh(log(2)/2*BWoct)
	#coefficients
	b0 = A*((A+1)-(A-1)*cos(w0)+2*sqrt(A)*alpha)
	b1 = 2*A*((A-1)-(A+1)*cos(w0))
	b2 = A*((A+1)-(A-1)*cos(w0)-2*sqrt(A)*alpha)
	a0 = (A+1)+(A-1)*cos(w0)+2*sqrt(A)*alpha
	a1 = -2*((A-1)+(A+1)*cos(w0))
	a2 = (A+1)+(A-1)*cos(w0)-2*sqrt(A)*alpha
	#normalize
	result = [b0 b1 b2 a1 a2]./a0
	return SecondOrderSections(Biquad(result...))
end


text:
abstract type FitBackend end
abstract type OptimFitBackend <: FitBackend end

backend_type(::Val{:optim}) = OptimFitBackend

fit(args...; backend::Symbol = :optim, kwargs...) = fit( backend_type(Val(backend)), args...; kwargs...)
fit(::Type{OptimFitBackend}, args...; kwargs...) = opt_fit(args...; kwargs...)

ext: jl


text:
function nontrivial_intersection(dim,N)
    intersecting_volumes = zeros(Float64, N)
    for i = 1:N
        S1,S2 = nontrivially_intersecting_simplices(dim)
        intersecting_volumes[i] = simplexintersection(copy(transpose(S1)), copy(transpose(S2)))
    end
    return intersecting_volumes
end

function shared_vertex_intersection(dim,N)
    intersecting_volumes = zeros(Float64, N)
    for i = 1:N
        S1,S2 = simplices_sharing_vertices(dim)
        intersecting_volumes[i] = simplexintersection(copy(transpose(S1)), copy(transpose(S2)))
    end
    return intersecting_volumes

end


text:
# Constraint template: Power balance constraint including candidate storage
#=function constraint_power_balance_acne_dcne_strg(pm::_PM.AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = PowerModels.ref(pm, nw, :bus, i)
    bus_arcs = PowerModels.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = PowerModels.ref(pm, nw, :ne_bus_arcs, i)
    bus_arcs_dc = PowerModels.ref(pm, nw, :bus_arcs_dc, i)
    bus_gens = PowerModels.ref(pm, nw, :bus_gens, i)
    bus_convs_ac = PowerModels.ref(pm, nw, :bus_convs_ac, i)
    bus_convs_ac_ne = PowerModels.ref(pm, nw, :bus_convs_ac_ne, i)
    bus_loads = PowerModels.ref(pm, nw, :bus_loads, i)
    bus_shunts = PowerModels.ref(pm, nw, :bus_shunts, i)
    bus_storage = PowerModels.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = PowerModels.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => PowerModels.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => PowerModels.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_dcne_strg(pm, nw, i, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end
#Power balance constraint including candidate storage
function constraint_power_balance_acne_dcne_strg(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p = _PM.var(pm, n, :p)
    pg = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne = _PM.var(pm, n, :pconv_ac_ne)
    p_ne = _PM.var(pm, n, :p_ne)
    ps   = _PM.var(pm, n, :ps)
    ps_ne   = _PM.var(pm, n, :ps_ne)
    v = 1

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) -sum(ps_ne[s] for s in bus_storage_ne) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end

# Constraint template: Power balance constraint (without DC equipment) including candidate storage
function constraint_power_balance_acne_strg(pm::_PM.AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = PowerModels.ref(pm, nw, :bus, i)
    bus_arcs = PowerModels.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = PowerModels.ref(pm, nw, :ne_bus_arcs, i)
    bus_gens = PowerModels.ref(pm, nw, :bus_gens, i)
    bus_loads = PowerModels.ref(pm, nw, :bus_loads, i)
    bus_shunts = PowerModels.ref(pm, nw, :bus_shunts, i)
    bus_storage = PowerModels.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = PowerModels.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => PowerModels.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => PowerModels.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_strg(pm, nw, i, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end
#Power balance constraint (without DC equipment) including candidate storage
function constraint_power_balance_acne_strg(pm::_PM.AbstractWModels, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p     = _PM.var(pm, n, :p)
    q     = _PM.var(pm, n, :q)
    p_ne  = _PM.var(pm, n, :p_ne)
    q_ne  = _PM.var(pm, n, :q_ne)
    pg    = _PM.var(pm, n, :pg)
    qg    = _PM.var(pm, n, :qg)
    ps    = _PM.var(pm, n, :ps)
    qs    = _PM.var(pm, n, :qs)
    ps_ne = _PM.var(pm, n, :ps_ne)
    qs_ne = _PM.var(pm, n, :qs_ne)
    w     = _PM.var(pm, n, :w, i)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) - sum(ps_ne[s] for s in bus_storage_ne) - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*w)
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) == sum(qg[g] for g in bus_gens) - sum(qs[s] for s in bus_storage) - sum(qs_ne[s] for s in bus_storage_ne) - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*w)
end=#

#Constraint template:  Power Balance constraint including candidate storage & flexible demand
#=function constraint_power_balance_acne_dcne_flex(pm::_PM.AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = PowerModels.ref(pm, nw, :bus, i)
    bus_arcs = PowerModels.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = PowerModels.ref(pm, nw, :ne_bus_arcs, i)
    bus_arcs_dc = PowerModels.ref(pm, nw, :bus_arcs_dc, i)
    bus_gens = PowerModels.ref(pm, nw, :bus_gens, i)
    bus_convs_ac = PowerModels.ref(pm, nw, :bus_convs_ac, i)
    bus_convs_ac_ne = PowerModels.ref(pm, nw, :bus_convs_ac_ne, i)
    bus_loads = PowerModels.ref(pm, nw, :bus_loads, i)
    bus_shunts = PowerModels.ref(pm, nw, :bus_shunts, i)
    bus_storage = PowerModels.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = PowerModels.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => PowerModels.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => PowerModels.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_dcne_flex(pm, nw, i, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end=#
# Power Balance constraint including candidate storage & flexible demand
#=function constraint_power_balance_acne_dcne_flex(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p = _PM.var(pm, n, :p)
    pg = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne = _PM.var(pm, n, :pconv_ac_ne)
    p_ne = _PM.var(pm, n, :p_ne)
    ps   = _PM.var(pm, n, :ps)
    ps_ne   = _PM.var(pm, n, :ps_ne)
    pflex = _PM.var(pm, n, :pflex)
    v = 1

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) -sum(ps_ne[s] for s in bus_storage_ne) - sum(pflex[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end=#

#Constraint template:  Power Balance constraint (without DC equipment) including candidate storage & flexible demand
#=function constraint_power_balance_acne_flex(pm::_PM.AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = PowerModels.ref(pm, nw, :bus, i)
    bus_arcs = PowerModels.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = PowerModels.ref(pm, nw, :ne_bus_arcs, i)
    bus_gens = PowerModels.ref(pm, nw, :bus_gens, i)
    bus_loads = PowerModels.ref(pm, nw, :bus_loads, i)
    bus_shunts = PowerModels.ref(pm, nw, :bus_shunts, i)
    bus_storage = PowerModels.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = PowerModels.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => PowerModels.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => PowerModels.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)
    constraint_power_balance_acne_flex(pm, nw, i, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end=#
# Power Balance constraint (without DC equipment) including candidate storage & flexible demand
#=function constraint_power_balance_acne_flex(pm::_PM.AbstractWModels, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_gens, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p                = _PM.var(pm, n, :p)
    q                = _PM.var(pm, n, :q)
    p_ne             = _PM.var(pm, n, :p_ne)
    q_ne             = _PM.var(pm, n, :q_ne)
    pg               = _PM.var(pm, n, :pg)
    qg               = _PM.var(pm, n, :qg)
    ps               = _PM.var(pm, n, :ps)
    qs               = _PM.var(pm, n, :qs)
    ps_ne            = _PM.var(pm, n, :ps_ne)
    qs_ne            = _PM.var(pm, n, :qs_ne)
    pflex            = _PM.var(pm, n, :pflex)
    qflex            = _PM.var(pm, n, :qflex)
    w                = _PM.var(pm, n, :w, i)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) - sum(ps_ne[s] for s in bus_storage_ne) - sum(pflex[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*w)
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) == sum(qg[g] for g in bus_gens) - sum(qs[s] for s in bus_storage) - sum(qs_ne[s] for s in bus_storage_ne) - sum(qflex[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*w)
end=#

#Constraint template:  Power Balance constraint with power interrupted by contingency - pinter (including candidate storage & flexible demand)
#=function constraint_power_balance_reliability(pm::_PM.AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = PowerModels.ref(pm, nw, :bus, i)
    bus_arcs = PowerModels.ref(pm, nw, :bus_arcs, i)
    bus_arcs_ne = PowerModels.ref(pm, nw, :ne_bus_arcs, i)
    bus_arcs_dc = PowerModels.ref(pm, nw, :bus_arcs_dc, i)
    bus_gens = PowerModels.ref(pm, nw, :bus_gens, i)
    bus_convs_ac = PowerModels.ref(pm, nw, :bus_convs_ac, i)
    bus_convs_ac_ne = PowerModels.ref(pm, nw, :bus_convs_ac_ne, i)
    bus_loads = PowerModels.ref(pm, nw, :bus_loads, i)
    bus_shunts = PowerModels.ref(pm, nw, :bus_shunts, i)
    bus_storage = PowerModels.ref(pm, nw, :bus_storage, i)
    bus_storage_ne = PowerModels.ref(pm, nw, :bus_storage_ne, i)

    pd = Dict(k => PowerModels.ref(pm, nw, :load, k, "pd") for k in bus_loads)
    qd = Dict(k => PowerModels.ref(pm, nw, :load, k, "qd") for k in bus_loads)

    gs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bs = Dict(k => PowerModels.ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)

    constraint_power_balance_reliability(pm, nw, i, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
end=#

# Power Balance constraint  with power interrupted by contingency - pinter (including candidate storage & flexible demand)
#=function constraint_power_balance_reliability(pm::_PM.AbstractDCPModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, bus_storage, bus_storage_ne, pd, qd, gs, bs)
    p = _PM.var(pm, n, :p)
    pg = _PM.var(pm, n, :pg)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    pconv_ac = _PM.var(pm, n, :pconv_ac)
    pconv_ac_ne = _PM.var(pm, n, :pconv_ac_ne)
    p_ne = _PM.var(pm, n, :p_ne)
    ps   = _PM.var(pm, n, :ps)
    ps_ne   = _PM.var(pm, n, :ps_ne)
    pflex = _PM.var(pm, n, :pflex)
    pinter = _PM.var(pm, n, :pinter)
    v = 1

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens) - sum(ps[s] for s in bus_storage) -sum(ps_ne[s] for s in bus_storage_ne) - sum(pflex[d] for d in bus_loads) + sum(pinter[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*v^2)
end=#


text:
using ImputedGenotypes
using Base.Test

bases = [ "A", "G", "C", "T" ]
function random_base_pair()
	pair=[ bases[rand(1:4)] for i in 1:2 ]
	while first(pair)==last(pair)
		pair=[ bases[rand(1:4)] for i in 1:2 ]
	end
	return pair
end

random_genotype() = join(map(x->@sprintf("%0.3f",x),normalize([rand(),rand(),rand()],1)),":")

generate_snp(position) = vcat(position,random_base_pair() )
generate_individual(i) = "indv$i"
generate_n(func,n) = [ func(i) for i in 1:n ]

vcf_header = "##fileformat=VCFv4.1\n##FORMAT=<ID=GP,Number=3,Type=Float,Description=\"...\">"
vcf_header_line(individuals) = join(vcat(["#CHROM","POS","ID","REF","ALT","QUAL","FILTER","INFO","FORMAT"],individuals),'\t')
vcf_line(snp,individuals) = join(vcat("42",snp,join(snp,'_'),".",".","GP",[random_genotype() for i in individuals]),'\t')

generate_imputed_vcf(snps, individuals) = join(vcat(vcf_header, vcf_header_line(individuals), [ vcf_line(snp,individuals) for snp in snps ]),'\n')

@test dosage_to_likelihoods(0.) == (1.,0.,0.)
@test dosage_to_likelihoods(1.) == (0.,1.,0.)
@test dosage_to_likelihoods(2.) == (0.,0.,1.)

@test likelihoods_to_dosage(1.,0.,0.) == (0.)
@test likelihoods_to_dosage(0.,1.,0.) == (1.)
@test likelihoods_to_dosage(1.,0.,1.) == (2.)

for i in 0:0.01:2
	@test i ≈ likelihoods_to_dosage(dosage_to_likelihoods(i))
end

@test proper_info([(0.3333333, 0.3333333, 0.3333333), (1.0, 0.0, 0.0)])≈0.428571367347

max_stars_repo_name: cth/merge-imputations

max_issues_repo_name: cth/merge-imputations

max_forks_repo_name: cth/merge-imputations


text:
#=
Author : Guy Farmer, Michael Hilborn
Runs simulations for all models (generated in Ambulance_Deployment_experiments.jl) and returns the data as a dictionary
=#

import JuMP, Gurobi
import DataStructures: PriorityQueue, enqueue!, dequeue!
import DataFrames: DataFrame, nrow
import Distributions
import Pkg
import CSV, Query
using Gurobi, CSV, Query, JuMP, Dates, JLD, GLPK
using DataFrames, CSV, Random, Plots, JSON, Distributions

using AmbulanceDeployment

function generate_simulation(model_name::String, namb::Int, ncalls::Int)

    turnaround = Distributions.LogNormal(3.65, 0.3)
    lambda = 0
    local_path = ""

    solverstats = JSON.parsefile(PROJECT_ROOT * "/src/outputs/solver_stats.json")
    hospitals = CSV.File(string(local_path, PROJECT_ROOT * "/test/austin-data/hospitals.csv")) |> DataFrame
    stations = CSV.File(string(local_path, PROJECT_ROOT * "/test/austin-data/stations.csv")) |> DataFrame
    hourly_calls = CSV.File(PROJECT_ROOT * "/test/austin-data/Full_WeekdayCalls.csv") |> DataFrame
    adjacent_nbhd = CSV.File(PROJECT_ROOT * "/test/austin-data/adjacent_nbhd.csv") |> DataFrame
    coverage = CSV.read(PROJECT_ROOT * "/test/austin-data/coverage_real.csv", DataFrame, header=false)
    coverage = convert(Array{Bool, 2}, coverage[:, :])
    incidents = CSV.File(PROJECT_ROOT * "/test/austin-data/austin_incidents.csv") |> DataFrame
    amb_deployment = solverstats["amb_deployment"]
    model_dict = Dict{String, Symbol}("Stochastic"=>:Stochastic, "Robust01"=>:Robust01, "Robust005"=>:Robust005, "Robust001"=>:Robust001, "Robust0001"=>:Robust0001,
    "Robust00001"=>:Robust00001, "MEXCLP"=>:MEXCLP, "MALP"=>:MALP)
    test_calls = CSV.File(PROJECT_ROOT *"/test/austin-data/austin_test_calls.csv")|> DataFrame
    # remember to reset number of calls
    test_calls = test_calls[1:ncalls,:]
    #filtering out invalid neighborhoods

    test_calls = filter(x->x[:neighborhood]>0,test_calls)

    p = DeploymentProblem(
    hourly_calls,
    adjacent_nbhd,
    coverage,
    namb = namb,
    train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 3)
    )

    # We focus on emergency calls during the "peak period" (8AM - 8PM),
    # with the emergency calls from the first 3 month as our training set,
    # and the subsequent emergency calls from the remaining months as our test set

    println("running $namb ambuances & $ncalls calls")
    x = amb_deployment[model_name][string(namb)]
    x = convert(Array{Int64}, x)
    problem = DispatchProblem(test_calls, hospitals, stations, p.coverage, x, turnaround=turnaround)
    dispatch = ClosestDispatch(p, problem)

    # id 145 dispatch to nbhd 88
    Random.seed!(1234); # reset seed
    @time df, guiArray = simulate_events!(problem, dispatch);
    println("wait time : ", df[!,:waittime])
    println("response time: ", df[!,:responsetime])

    d = df[!,:responsetime]
    r = d
    d = filter(x->x!=Inf,d)
    m =  mean(d)
    println("mean response time = $m")
    json_string = JSON.json(guiArray)
    open(PROJECT_ROOT *"/src/outputs/guiArray.json","w") do f
        write(f, json_string)
    end
    events = filter(x->typeof(x)==gui_event,guiArray)
    responded_array = filter(x->x.event_type == "call responded", events)
    arrived_array = filter(x->x.event_type == "ambulance arrived", events)
    return guiArray, responded_array, arrived_array
    #results[j,i] = mean(df[!,:waittime] + df[!,:responsetime])
end

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
"""
    ilc_weights(cij)

This function returns weights (a vector of the number of frequency channels) of the ILC method.

*Reference*: Equation (12) of Tegmark et al., Phys. Rev. D, 68, 123523 (2003)

# Arguments
- `cij::Array{<:AbstractFloat,2}`: symmetric covariance matrix with the dimention of `(nν, nν)` where `nν` is the number of frequency bands.
"""
function ilc_weights(cij::Array{<:AbstractFloat,2})
    if size(cij)[1] ≠ size(cij)[2]
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    else
        M = Symmetric(cij)
    end
    nν = size(cij)[1]
    e = ones(nν) # e = [1,...,1]
    x = M \ e
    w = x / e'x  # ILC weights
end

"""
    ilc_weights(cij, e)

This function returns weights (a vector of the number of frequency channels) of the ILC method.

*Reference*: Equation (12) of Tegmark et al., Phys. Rev. D, 68, 123523 (2003)

# Arguments
- `cij::Array{<:AbstractFloat,2}`: symmetric covariance matrix with the dimention of `(nν, nν)` where `nν` is the number of frequency bands.
- `e::Array{<:AbstractFloat,1}`: vector of the frequency response. E.g., `e = [1,...,1]` for CMB.
"""
function ilc_weights(cij::Array{T,2}, e::Array{T,1}) where {T<:AbstractFloat}
    if size(cij)[1] ≠ size(cij)[2]
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    elseif size(cij)[1] ≠ length(e)
        throw(DimensionMismatch("dimensions of the covariance matrix and the frequency response vector do not match"))
    else
        M = Symmetric(cij)
    end
    x = M \ e
    w = x / e'x  # ILC weights
end

"""
    ilc_weights(cijℓ[, ℓid=3])

This function returns weights (a `nν`-by-`nℓ` matrix) of the ILC method.

Here, `nν` is the number of frequency channels and `nℓ` is the number of elements in the relevant domain, e.g., multipoles, band-power bins, pixels, etc.

*Reference*: Equation (12) of Tegmark et al., Phys. Rev. D, 68, 123523 (2003)

# Arguments
- `cijℓ::Array{<:AbstractFloat,3}`: symmetric covariance matrix with the dimention of `(nℓ, nν, nν)`, `(nν, nℓ, nν)` or `(nν, nν, nℓ)` (default).

# Optional Arguments
- `ℓid::Integer=3`: location of the index for the `nℓ` domain. `ℓid=1` if `cijℓ[nℓ,nν,nν]`, `ℓid=2` if `cijℓ[nν,nℓ,nν]`, and `ℓid=3` (the default value) if `cijℓ[nν,nν,nℓ]`.
"""
function ilc_weights(
    cijℓ::Array{<:AbstractFloat,3},
    ℓid::Integer = 3,
)
    if ℓid > 3 || ℓid < 1
        throw(DomainError(ℓid, "ℓid must be 1, 2, or 3"))
    end
    if (ℓid == 3 && size(cijℓ)[1] ≠ size(cijℓ)[2]) ||
       (ℓid == 2 && size(cijℓ)[1] ≠ size(cijℓ)[3]) ||
       (ℓid == 1 && size(cijℓ)[2] ≠ size(cijℓ)[3])
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    end
    nℓ = size(cijℓ)[ℓid]
    nν = ifelse(ℓid == 3, size(cijℓ)[1], size(cijℓ)[3])
    e = ones(nν)      # e = [1,...,1]
    wℓ = zeros(nν, nℓ) # ILC weights
    for iℓ = 1:nℓ
        if ℓid == 3
            M = Symmetric(cijℓ[:, :, iℓ])
        elseif ℓid == 2
            M = Symmetric(cijℓ[:, iℓ, :])
        else
            M = Symmetric(cijℓ[iℓ, :, :])
        end
        x = M \ e
        wℓ[:, iℓ] = x / e'x
    end
    return wℓ
end

"""
    ilc_weights(cijℓ, e[, ℓid=3])

This function returns weights (a `nν`-by-`nℓ` matrix) of the ILC method.

Here, `nν` is the number of frequency channels and `nℓ` is the number of elements in the relevant domain, e.g., multipoles, band-power bins, pixels, etc.

*Reference*: Equation (12) of Tegmark et al., Phys. Rev. D, 68, 123523 (2003)

# Arguments
- `cijℓ::Array{<:AbstractFloat,3}`: symmetric covariance matrix with the dimention of `(nℓ, nν, nν)`, `(nν, nℓ, nν)` or `(nν, nν, nℓ)` (default).
- `e::Array{<:AbstractFloat,1}`: vector of the frequency response. E.g., `e = [1,...,1]` for CMB.

# Optional Arguments
- `ℓid::Integer=3`: location of the index for the `nℓ` domain. `ℓid=1` if `cijℓ[nℓ,nν,nν]`, `ℓid=2` if `cijℓ[nν,nℓ,nν]`, and `ℓid=3` (the default value) if `cijℓ[nν,nν,nℓ]`.
"""
function ilc_weights(
    cijℓ::Array{T,3},
    e::Array{T,1},
    ℓid::Integer = 3,
) where {T<:AbstractFloat}
    if ℓid > 3 || ℓid < 1
        throw(DomainError(ℓid, "ℓid must be 1, 2, or 3"))
    end
    if (ℓid == 3 && size(cijℓ)[1] ≠ size(cijℓ)[2]) ||
       (ℓid == 2 && size(cijℓ)[1] ≠ size(cijℓ)[3]) ||
       (ℓid == 1 && size(cijℓ)[2] ≠ size(cijℓ)[3])
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    end
    nℓ = size(cijℓ)[ℓid]
    nν = ifelse(ℓid == 3, size(cijℓ)[1], size(cijℓ)[3])
    if length(e) ≠ nν
        throw(DimensionMismatch("dimensions of the covariance matrix and the frequency response vector do not match"))
    end
    wℓ = zeros(nν, nℓ) # ILC weights
    for iℓ = 1:nℓ
        if ℓid == 3
            M = Symmetric(cijℓ[:, :, iℓ])
        elseif ℓid == 2
            M = Symmetric(cijℓ[:, iℓ, :])
        else
            M = Symmetric(cijℓ[iℓ, :, :])
        end
        x = M \ e
        wℓ[:, iℓ] = x / e'x
    end
    return wℓ
end

"""
    ilc_clean_cij(cij, w)

This function returns power of the extracted component for a given element of the relevant domain, e.g., multipole, band-power bin, pixel, etc.

*Reference*: Tegmark et al., Phys. Rev. D, 68, 123523 (2003)

# Arguments
- `cij::Array{<:AbstractFloat,2}`: symmetric covariance matrix with the dimention of `(nν, nν)` where `nν` is the number of frequency bands.
- `w::Array{<:AbstractFloat,1}`: ILC weights.
"""
function ilc_clean_cij(cij::Array{T,2}, w::Array{T,1}) where {T<:AbstractFloat}
    if size(cij)[1] ≠ size(cij)[2]
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    elseif size(cij)[1] ≠ length(w)
        throw(DimensionMismatch("dimensions of the covariance matrix and the ILC weight do not match"))
    else
        M = Symmetric(cij)
    end
    cl = w' * M * w
end

"""
    ilc_clean_cij(cijℓ, wℓ[, ℓid::Integer=3])

This function returns a vector of the power of the extracted component, with elements in the relevant domain, e.g., multipole, band-power bin, pixel, etc.

*Reference*: Tegmark et al., Phys. Rev. D, 68, 123523 (2003)

# Arguments
- `cijℓ::Array{<:AbstractFloat,3}`: symmetric covariance matrix with the dimention of `(nℓ, nν, nν)`, `(nν, nℓ, nν)` or `(nν, nν, nℓ)` (default).
- `wℓ::Array{<:AbstractFloat,2}`: ILC weights (a `nν`-by-`nℓ` matrix).
    - Here, `nν` is the number of frequency bands and `nℓ` is the number of elements in the relevant domain.

# Optional Arguments
- `ℓid::Integer=3`: location of the index for the `nℓ` domain. `ℓid=1` if `cijℓ[nℓ,nν,nν]`, `ℓid=2` if `cijℓ[nν,nℓ,nν]`, and `ℓid=3` (the default value) if `cijℓ[nν,nν,nℓ]`.
"""
function ilc_clean_cij(
    cijℓ::Array{T,3},
    wℓ::Array{T,2},
    ℓid::Integer = 3,
) where {T<:AbstractFloat}
    if ℓid > 3 || ℓid < 1
        throw(DomainError(ℓid, "ℓid must be 1, 2, or 3"))
    end
    if (ℓid == 3 && size(cijℓ)[1] ≠ size(cijℓ)[2]) ||
       (ℓid == 2 && size(cijℓ)[1] ≠ size(cijℓ)[3]) ||
       (ℓid == 1 && size(cijℓ)[2] ≠ size(cijℓ)[3])
        throw(DimensionMismatch("covariance matrix must be a square matrix"))
    end
    nℓ = size(cijℓ)[ℓid]
    nν = ifelse(ℓid == 3, size(cijℓ)[1], size(cijℓ)[3])
    if nℓ ≠ size(wℓ)[2] || nν ≠ size(wℓ)[1]
        throw(DimensionMismatch("dimensions of the covariance matrix and the ILC weight do not match"))
    end
    cl = zeros(nℓ)
    for iℓ = 1:nℓ
        if ℓid == 3
            M = Symmetric(cijℓ[:, :, iℓ])
        elseif ℓid == 2
            M = Symmetric(cijℓ[:, iℓ, :])
        else
            M = Symmetric(cijℓ[iℓ, :, :])
        end
        cl[iℓ] = wℓ[:, iℓ]' * M * wℓ[:, iℓ]
    end
    return cl
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    _SOCtoPSDaff{T}(f::MOI.VectorAffineFunction{T}, g::MOI.ScalarAffineFunction{T})

Builds a VectorAffineFunction representing the upper (or lower) triangular part of the matrix
[ f[1]     f[2:end]' ]
[ f[2:end] g * I     ]
"""
function _SOCtoPSDaff(
    ::Type{T},
    f::Union{MOI.AbstractVectorFunction,AbstractVector{T}},
    g::Union{MOI.AbstractScalarFunction,T},
) where {T}
    F = MOIU.promote_operation(vcat, T, typeof(g), T)
    f_scalars = MOIU.eachscalar(f)
    dim = length(f_scalars)
    n = div(dim * (dim + 1), 2)
    h = MOIU.zero_with_output_dimension(F, n)
    MOIU.operate_output_index!(+, T, MOIU.trimap(1, 1), h, f_scalars[1])
    for i in 2:dim
        MOIU.operate_output_index!(+, T, MOIU.trimap(1, i), h, f_scalars[i])
        MOIU.operate_output_index!(+, T, MOIU.trimap(i, i), h, g)
    end
    return h
end

"""
The `SOCtoPSDBridge` transforms the second order cone constraint ``\\lVert x \\rVert \\le t`` into the semidefinite cone constraints
```math
\\begin{pmatrix}
  t & x^\\top\\\\
  x & tI
\\end{pmatrix} \\succeq 0
```
Indeed by the Schur Complement, it is positive definite iff
```math
\\begin{align*}
  tI & \\succ 0\\\\
  t - x^\\top (tI)^{-1} x & \\succ 0
\\end{align*}
```
which is equivalent to
```math
\\begin{align*}
  t & > 0\\\\
  t^2 & > x^\\top x
\\end{align*}
```

!!! warning
    This bridge is not added by default by [`MOI.Bridges.full_bridge_optimizer`](@ref)
    as bridging second order cone constraints to semidefinite constraints can be
    achieved by the [`SOCtoRSOCBridge`](@ref) followed by the [`RSOCtoPSDBridge`](@ref)
    while creating a smaller semidefinite constraint.
"""
struct SOCtoPSDBridge{T,F,G} <: SetMapBridge{
    T,
    MOI.PositiveSemidefiniteConeTriangle,
    MOI.SecondOrderCone,
    F,
    G,
}
    constraint::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle}
end

const SOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoPSDBridge{T},OT}

function concrete_bridge_type(
    ::Type{<:SOCtoPSDBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.SecondOrderCone},
) where {T}
    F = MOIU.promote_operation(vcat, T, MOIU.scalar_type(G), T)
    return SOCtoPSDBridge{T,F,G}
end

function MOIB.map_set(::Type{<:SOCtoPSDBridge}, set::MOI.SecondOrderCone)
    return MOI.PositiveSemidefiniteConeTriangle(MOI.dimension(set))
end

function MOIB.inverse_map_set(
    ::Type{<:SOCtoPSDBridge},
    set::MOI.PositiveSemidefiniteConeTriangle,
)
    return MOI.SecondOrderCone(MOI.side_dimension(set))
end

function MOIB.map_function(::Type{<:SOCtoPSDBridge{T}}, func) where {T}
    return _SOCtoPSDaff(T, func, MOIU.eachscalar(func)[1])
end

function MOIB.inverse_map_function(::Type{<:SOCtoPSDBridge}, func)
    scalars = MOIU.eachscalar(func)
    dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars))
    return scalars[MOIU.trimap.(1, 1:dim)]
end

function MOIB.adjoint_map_function(::Type{<:SOCtoPSDBridge{T}}, func) where {T}
    scalars = MOIU.eachscalar(func)
    dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars))
    tdual = sum(i -> func[MOIU.trimap(i, i)], 1:dim)
    return MOIU.operate(vcat, T, tdual, func[MOIU.trimap.(2:dim, 1)] * 2)
end

function MOIB.inverse_adjoint_map_function(
    ::Type{<:SOCtoPSDBridge{T}},
    func,
) where {T}
    # func is (t, x) such that x'x ≤ t^2 and we need to find a PSD matrix
    # [a   x'/2]
    # [x/2 Q   ]
    # such that a + tr(Q) == t
    # By choosing a = t/2 and Q = t/(2x'x) * xx', we get
    # a + tr(Q) = t/2 + t/(2x'x) * tr(xx') = t/2 + t/(2x'x) * tr(x'x) = t
    # Moreover, by the Schur complement, the matrix is PSD iff
    # xx'/(2t) ⪯ t/(2x'x) * xx'
    # 1/(2t) ≤ t/(2x'x)
    # x'x ≤ t^2
    # which is the SOC inequality
    t = func[1]
    x = func[2:end]
    return MOIB.inverse_adjoint_map_function(RSOCtoPSDBridge{T}, [t / 2; t; x])
end

"""
The `RSOCtoPSDBridge` transforms the second order cone constraint
``\\lVert x \\rVert \\le 2tu`` with ``u \\ge 0`` into the semidefinite cone
constraints
```math
\\begin{pmatrix}
  t & x^\\top\\\\
  x & 2uI
\\end{pmatrix} \\succeq 0
```

Indeed by the Schur Complement, it is positive definite iff
```math
\\begin{align*}
  uI & \\succ 0\\\\
  t - x^\\top (2uI)^{-1} x & \\succ 0
\\end{align*}
```
which is equivalent to
```math
\\begin{align*}
  u & > 0\\\\
  2tu & > x^\\top x
\\end{align*}
```
"""
struct RSOCtoPSDBridge{T,F,G} <: SetMapBridge{
    T,
    MOI.PositiveSemidefiniteConeTriangle,
    MOI.RotatedSecondOrderCone,
    F,
    G,
}
    constraint::MOI.ConstraintIndex{F,MOI.PositiveSemidefiniteConeTriangle}
end

const RSOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoPSDBridge{T},OT}

function concrete_bridge_type(
    ::Type{<:RSOCtoPSDBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.RotatedSecondOrderCone},
) where {T}
    S = MOIU.scalar_type(G)
    H = MOIU.promote_operation(*, T, T, S)
    F = MOIU.promote_operation(vcat, T, S, H, T)
    return RSOCtoPSDBridge{T,F,G}
end

function MOIB.map_set(
    ::Type{<:RSOCtoPSDBridge},
    set::MOI.RotatedSecondOrderCone,
)
    return MOI.PositiveSemidefiniteConeTriangle(MOI.dimension(set) - 1)
end

function MOIB.inverse_map_set(
    ::Type{<:RSOCtoPSDBridge},
    set::MOI.PositiveSemidefiniteConeTriangle,
)
    return MOI.RotatedSecondOrderCone(MOI.side_dimension(set) + 1)
end

function MOIB.map_function(::Type{<:RSOCtoPSDBridge{T}}, func) where {T}
    scalars = MOIU.eachscalar(func)
    if length(scalars) < 3
        error(
            "Unable to bridge RotatedSecondOrderCone to PSD because the ",
            "dimension is too small: got $(length(scalars)), expected >= 3.",
        )
    end
    h = MOIU.operate!(*, T, scalars[2], convert(T, 2))
    return _SOCtoPSDaff(T, scalars[[1; 3:length(scalars)]], h)
end

function MOIB.inverse_map_function(::Type{<:RSOCtoPSDBridge{T}}, func) where {T}
    scalars = MOIU.eachscalar(func)
    dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars))
    t = scalars[1]
    # It is (2u*I)[1,1] so it needs to be divided by 2 to get u
    u = MOIU.operate!(/, T, scalars[3], convert(T, 2))
    return MOIU.operate(
        vcat,
        T,
        t,
        u,
        scalars[[MOIU.trimap(1, i) for i in 2:dim]],
    )
end

function MOIB.adjoint_map_function(::Type{<:RSOCtoPSDBridge{T}}, func) where {T}
    scalars = MOIU.eachscalar(func)
    dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars))
    udual = sum(i -> func[MOIU.trimap(i, i)], 2:dim)
    return MOIU.operate(
        vcat,
        T,
        func[1],
        2udual,
        func[MOIU.trimap.(2:dim, 1)] * 2,
    )
end

function MOIB.inverse_adjoint_map_function(
    ::Type{<:RSOCtoPSDBridge{T}},
    func,
) where {T}
    # func is (t, u, x) such that x'x ≤ 2tu and we need to find a PSD matrix
    # [t   x'/2]
    # [x/2 Q   ]
    # such that 2tr(Q) == u
    # By choosing Q = u/(2x'x) * xx', we get
    # 2tr(Q) = u/(x'x) * tr(xx') = u/(x'x) * tr(x'x) = u
    # Moreover, by the Schur complement, the matrix is PSD iff
    # xx'/(4t) ⪯ u/(2x'x) * xx'
    # 1/(4t) ≤ u/(2x'x)
    # x'x ≤ 2tu
    # which is the RSOC inequality
    t = func[1]
    u = func[2]
    x = func[3:end]
    Q = (x * x') * (u / (2 * (x' * x)))
    M = [
        t x'/2
        x/2 Q
    ]
    return [M[i, j] for j in 1:size(M, 2) for i in 1:j]
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
#!/usr/bin/env julia

Fract_PT = RheoModelClass(
          # Model name
          name="fPT",
          # Model parameters,
          p = (:cₐ, :a, :cᵦ, :β, :cᵧ, :γ),
          # Relaxation modulus
          G = quote
                  G(s) = (1/s)*(cᵧ*s^γ *(cₐ * s^a + cᵦ * s^β))/(cᵧ*s^γ + cₐ*s^a + cᵦ*s^β)
                  InverseLaplace.talbot(s -> G(s), t)
              end,
          # Creep modulus
          J = quote
                  t^a/cₐ* mittleff(a - β, 1 + a, -cᵦ*t^(a-β)/cₐ) + t^γ / (cᵧ * gamma(1 + γ))
              end,
         # Storage modulus
          Gp = quote
                  denominator = (cₐ * ω^a)^2 + (cᵦ * ω^β)^2 + (cᵧ * ω^γ)^2 + 2* (cₐ * ω^a) * (cᵦ * ω^β) * cos((a-β)*π/2) + 2* (cₐ * ω^a) * (cᵧ * ω^γ) * cos((a-γ)*π/2) + 2* (cᵦ * ω^β) * (cᵧ * ω^γ)* cos((β-γ)*π/2)
                  numerator = (cᵧ * ω^γ) * cos(γ*π/2) * ((cₐ * ω^a)^2 + (cᵦ * ω^β)^2) + (cᵧ * ω^γ)^2 * ((cₐ * ω^a)*cos(a*π/2)+(cᵦ * ω^β)*cos(β*π/2)) + (cₐ * ω^a) * (cᵦ * ω^β) * (cᵧ * ω^γ) * (cos((a-β-γ)*π/2)+cos((β-a-γ)*π/2))
                  numerator/denominator
               end,
         # Loss modulus
          Gpp = quote
                  denominator = (cₐ * ω^a)^2 + (cᵦ * ω^β)^2 + (cᵧ * ω^γ)^2 + 2* (cₐ * ω^a) * (cᵦ * ω^β) * cos((a-β)*π/2) + 2* (cₐ * ω^a) * (cᵧ * ω^γ) * cos((a-γ)*π/2) + 2* (cᵦ * ω^β) * (cᵧ * ω^γ)* cos((β-γ)*π/2)
                  numerator = (cᵧ * ω^γ) * sin(γ*π/2) * ((cₐ * ω^a)^2 + (cᵦ * ω^β)^2) + (cᵧ * ω^γ)^2 * ((cₐ * ω^a)*sin(a*π/2)+(cᵦ * ω^β)*sin(β*π/2)) + (cₐ * ω^a) * (cᵦ * ω^β) * (cᵧ * ω^γ) * (sin((a-β-γ)*π/2)+sin((β-a-γ)*π/2))
                  numerator/denominator
                end,
        # Constraints
        constraint = quote
                 all([   (a<1) & (a>0)
                         (β<1) & (β>0)
                          -a+β < 0
                         (γ<1) & (γ>0)] )
                end,
          # Network
          info= "
                     _________╱╲_________
                    |         ╲╱ cₐ, a   |
                ____|                    |______╱╲____
                    |                    |      ╲╱
                    |_________╱╲_________|        cᵧ, γ
                              ╲╱
                                 cᵦ, β
                                 "
          )

FractSLS_PT = RheoModelClass(
          # Model name
          name="fSLS_PT",
          # Model parameters,
          p = (:cₐ, :a, :kᵦ, :kᵧ),
          # Relaxation modulus
          G = quote
                Zkᵦ = (kᵧ)^2/(kᵦ + kᵧ)
                Zcₐ = cₐ * (kᵧ)^2 / (kᵦ + kᵧ)^2
                Zkᵧ = kᵦ * kᵧ / (kᵦ + kᵧ)
                Zkᵦ*mittleff(a, -(Zkᵦ/Zcₐ)*t^a) + Zkᵧ
              end,
          # Creep modulus
          J = quote
                  1/kᵧ + t^(a)/cₐ * mittleff(a, 1 + a, -kᵦ*t^(a)/cₐ)
              end,
          # Storage modulus
          Gp = quote
                  Zkᵦ = (kᵧ)^2/(kᵦ + kᵧ)
                  Zcₐ = cₐ * (kᵧ)^2 / (kᵦ + kᵧ)^2
                  Zkᵧ = kᵦ * kᵧ / (kᵦ + kᵧ)
                  denominator = (Zcₐ*ω^a)^2 + Zkᵦ^2 + 2*(Zcₐ*ω^a)*Zkᵦ*cos(a*π/2)
                  numerator = Zkᵦ^2*(Zcₐ*ω^a)*cos(a*π/2) + ((Zcₐ*ω^a)^2)*Zkᵦ
                  numerator/denominator + Zkᵧ
               end,
         # Loss modulus
          Gpp = quote
                  Zkᵦ = (kᵧ)^2/(kᵦ + kᵧ)
                  Zcₐ = cₐ * (kᵧ)^2 / (kᵦ + kᵧ)^2
                  Zkᵧ = kᵦ * kᵧ / (kᵦ + kᵧ)
                  denominator = (Zcₐ*ω^a)^2 + Zkᵦ^2 + 2*(Zcₐ*ω^a)*Zkᵦ*cos(a*π/2)
                  numerator = Zkᵦ^2*(Zcₐ*ω^a)*sin(a*π/2)
                  numerator/denominator
                end,
        # Constraints
          constraint = quote
                 (a<1) & (a>0)
                end,
          # Network
          info= "
                     _________╱╲_________
                    |         ╲╱ cₐ, a   |
                ____|                    |______╱╲  ╱╲  ╱╲  ____
                    |                    |        ╲╱  ╲╱  ╲╱
                    |____╱╲  ╱╲  ╱╲  ____|                   kᵧ
                           ╲╱  ╲╱  ╲╱
                                     kᵦ
                                 "
          )

SLS_PT = RheoModelClass(
          # Model name
          name="SLS_PT",
          # Model parameters,
          p = (:η, :kᵦ, :kᵧ),
          # Relaxation modulus
          G = quote
                Zkᵦ = (kᵧ)^2/(kᵦ + kᵧ)
                Zη = η * (kᵧ)^2 / (kᵦ + kᵧ)^2
                Zkᵧ = kᵦ * kᵧ / (kᵦ + kᵧ)
                Zkᵧ + Zkᵦ*exp(-t*Zkᵦ/Zη)
              end,
          # Creep modulus
          J = quote
                  1/kᵧ + (1 - exp(-kᵦ*t/η))/kᵦ
              end,
          # Storage modulus
          Gp = quote
                  Zkᵦ = (kᵧ)^2/(kᵦ + kᵧ)
                  Zη = η * (kᵧ)^2 / (kᵦ + kᵧ)^2
                  Zkᵧ = kᵦ * kᵧ / (kᵦ + kᵧ)
                  τ = Zη/Zkᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω^2*τ^2*Zkᵦ
                  numerator/denominator + Zkᵧ
               end,
         # Loss modulus
          Gpp = quote
                  Zkᵦ = (kᵧ)^2/(kᵦ + kᵧ)
                  Zη = η * (kᵧ)^2 / (kᵦ + kᵧ)^2
                  Zkᵧ = kᵦ * kᵧ / (kᵦ + kᵧ)
                  τ = Zη/Zkᵦ
                  denominator = 1 + τ^2*ω^2
                  numerator = ω*τ*Zkᵦ
                  numerator/denominator
                end,
          # Network
          info= "
                             ___
                     _________| |________
                    |        _|_| η      |
                ____|                    |______╱╲  ╱╲  ╱╲  ____
                    |                    |        ╲╱  ╲╱  ╲╱
                    |____╱╲  ╱╲  ╱╲  ____|                   kᵧ
                           ╲╱  ╲╱  ╲╱
                                     kᵦ
                                 "
          )

FractJeffreys_PT = RheoModelClass(
              # Model name
              name="fjeff_PT",
              # Model parameters,
              p = (:ηₐ, :cᵦ, :β, :ηᵧ),
              # Relaxation modulus
              G = quote
                      Zηₐ = (ηᵧ)^2/(ηₐ + ηᵧ)
                      Zcᵦ = cᵦ * (ηᵧ)^2 / (ηₐ + ηᵧ)^2
                      Zηᵧ = ηₐ * ηᵧ / (ηₐ + ηᵧ)
                      diracterm = t!=0.0 ? 0.0 : Inf
                      Zcᵦ*t^(-β)*mittleff(1-β, 1-β, -Zcᵦ*t^(1-β)/Zηₐ) + Zηᵧ*diracterm
                  end,
              # Creep modulus
              J = quote
                     t/ηₐ * mittleff(1-β, 2, -cᵦ*t^(1-β)/ηₐ) + t/ηᵧ
                  end,
              # Storage modulus
              Gp = quote
                      Zηₐ = (ηᵧ)^2/(ηₐ + ηᵧ)
                      Zcᵦ = cᵦ * (ηᵧ)^2 / (ηₐ + ηᵧ)^2
                      Zηᵧ = ηₐ * ηᵧ / (ηₐ + ηᵧ)
                      denominator = (Zηₐ*ω)^2 + (Zcᵦ*ω^β)^2 + 2*(Zηₐ*ω)*(Zcᵦ*ω^β)*cos((1 - β)*π/2)
                      numerator = ((Zηₐ*ω)^2)*(Zcᵦ*ω^β)*cos(β*π/2)
                      # cosine floating point error work-around
                      if β==1.0
                          0.0
                      else
                          numerator/denominator
                      end
                   end,
             # Loss modulus
              Gpp = quote
                      Zηₐ = (ηᵧ)^2/(ηₐ + ηᵧ)
                      Zcᵦ = cᵦ * (ηᵧ)^2 / (ηₐ + ηᵧ)^2
                      Zηᵧ = ηₐ * ηᵧ / (ηₐ + ηᵧ)
                      denominator = (Zηₐ*ω)^2 + (Zcᵦ*ω^β)^2 + 2*(Zηₐ*ω)*(Zcᵦ*ω^β)*cos((1 - β)*π/2)
                      numerator = ((Zcᵦ*ω^β)^2)*(Zηₐ*ω) + ((Zηₐ*ω)^2)*(Zcᵦ*ω^β)*sin(β*π/2)
                      numerator/denominator + Zηᵧ*ω
                    end,
              # Constraints
              constraint = quote
                       (β<1) & (β>0)
                      end,
              # Network
              info= "
                                      ___
                              _________| |________
                             |        _|_| ηₐ     |        ___
                         ____|                    |_________| |_____
                             |                    |        _|_| ηᵧ
                             |_________╱╲_________|
                                       ╲╱
                                          cᵦ, β
                                          "
              )

Jeffreys_PT = RheoModelClass(
          # Model name
          name="jeffreys_PT",
          # Model parameters,
          p = (:ηₐ, :k, :ηᵧ),
          # Relaxation modulus
          G = quote
              Zηₐ = (ηᵧ)^2/(ηₐ + ηᵧ)
              Zk = k * (ηᵧ)^2 / (ηₐ + ηᵧ)^2
              Zηᵧ = ηₐ * ηᵧ / (ηₐ + ηᵧ)
              diracterm = t!=0.0 ? 0.0 : Inf
              Zk*exp(-Zk*t/Zηₐ) + Zηᵧ*diracterm
          end,
          # Creep modulus
          J = quote
                (- exp(-k*t/ηₐ) + 1 )/k + t/ηᵧ
              end,
          Gp = quote
                  Zηₐ = (ηᵧ)^2/(ηₐ + ηᵧ)
                  Zk = k * (ηᵧ)^2 / (ηₐ + ηᵧ)^2
                  Zηᵧ = ηₐ * ηᵧ / (ηₐ + ηᵧ)
                  denominator = (Zηₐ*ω)^2 + Zk^2
                  numerator = ((Zηₐ*ω)^2)*Zk
                  numerator/denominator
               end,
         # Loss modulus
          Gpp = quote
                  Zηₐ = (ηᵧ)^2/(ηₐ + ηᵧ)
                  Zk = k * (ηᵧ)^2 / (ηₐ + ηᵧ)^2
                  Zηᵧ = ηₐ * ηᵧ / (ηₐ + ηᵧ)
                  denominator = (Zηₐ*ω)^2 + Zk^2
                  numerator = (Zk^2)*(Zηₐ*ω)
                  numerator/denominator + Zηᵧ*ω
                end,
          # Network
          info= "
                             ___
                     _________| |________
                    |        _|_| ηₐ     |        ___
                ____|                    |_________| |_____
                    |                    |        _|_| ηᵧ
                    |____╱╲  ╱╲  ╱╲  ____|
                           ╲╱  ╲╱  ╲╱
                                     k
                                 "
          )


text: # # Training interface


text:
"""
    savegif(filename::String, o::Output; kw...)

Write the output array to a gif.

# Arguments

- `filename`: File path to save the gif file to.
- `output`: An [`Output`](@ref) object. Note that to make a gif, the output should stores
    frames, and run with `store=true`, and `@assert DynamicGrids.istored(o)` should pass.

# Keywords

[`ImageConfig`](@ref) keywords:

$IMAGECONFIG_KEYWORDS
"""
function savegif(filename::String, o::Output, ruleset=Ruleset();
    minval=minval(o), maxval=maxval(o),
    scheme=ObjectScheme(), renderer=autorenderer(init(o), scheme),
    font=autofont(), text=TextConfig(font=font), textconfig=text, kw...
)
    im_o = NoDisplayImageOutput(o;
        imageconfig=ImageConfig(init(o);
            minval=minval, maxval=maxval, renderer=renderer, textconfig=textconfig
        )
    )
    savegif(filename, im_o, ruleset; kw...)
end
function savegif(filename::String, o::ImageOutput, ruleset=Ruleset(); fps=fps(o), kw...)
    length(o) == 1 && @warn "The output has length 1: the saved gif will be a single image"
    simdata = SimData(o, ruleset)
    images = map(firstindex(o):lastindex(o)) do f
        @set! simdata.currentframe = f
        Array(render!(o, simdata, o[f]))
    end
    array = cat(images..., dims=3)
    FileIO.save(filename, array; fps=fps, kw...)
    array
end

"""
    GifOutput <: ImageOutput

    GifOutput(init; filename, tspan, kw...)

Output that stores the simulation as images and saves a Gif file on completion.

# Arguments:

- `init`: initialisation `AbstractArrayArray` or `NamedTuple` of `AbstractArrayArray`.

# Keywords

Storing the gif:
- `filename`: File path to save the gif file to.

$IMAGEOUTPUT_KEYWORDS
"""
mutable struct GifOutput{T,F<:AbstractVector{T},E,GC,IC,G,N} <: ImageOutput{T,F}
    frames::F
    running::Bool
    extent::E
    graphicconfig::GC
    imageconfig::IC
    gif::G
    filename::N
end
function GifOutput(; frames, running, extent, graphicconfig, imageconfig, filename, kw...)
    gif = _allocgif(imageconfig, extent)
    GifOutput(frames, running, extent, graphicconfig, imageconfig, gif, filename)
end

# Getters
filename(o::GifOutput) = o.filename
gif(o::GifOutput) = o.gif

# Output/ImageOutput interface methods
maybesleep(output::GifOutput, f) = nothing
showimage(image, o::GifOutput, data::AbstractSimData) = gif(o)[:, :, currentframe(data)] .= image
finalisegraphics(o::GifOutput, data::AbstractSimData) = savegif(o)

# The gif is already generated, just save it again if neccessary
savegif(o::GifOutput) = savegif(filename(o), o)
function savegif(filename::String, o::GifOutput, fps=fps(o); kw...)
    FileIO.save(filename, gif(o); fps=fps, kw...)
end

# Preallocate the 3 dimensional gif array
_allocgif(i::ImageConfig, e::Extent) = _allocgif(renderer(i), i::ImageConfig, e::Extent)
function _allocgif(::Renderer, i::ImageConfig, e::Extent)
    zeros(ARGB32, gridsize(e)..., length(tspan(e)))
end
function _allocgif(p::Layout, i::ImageConfig, e::Extent)
    zeros(ARGB32, (gridsize(e) .* size(p.layout))..., length(tspan(e)))
end


text:
using VarStructs
using Test

@testset "VarStructs.jl" begin

    @testset "parse" begin
        include("parse.jl")
    end

    @testset "@var" begin
        include("@var.jl")
    end

    @testset "@shared_var" begin
        include("@shared_var.jl")
    end

    include("../example/schema.jl")

end


text:
#these tests use 64-bit floating point values as the "ground truth".
#sometimes, the input values may need to be trimmed from the exact 64-bit value.
FloatFamily = Union{Posit, Float32, Float16, Float64}
trim_float(arr::Array{Float64}, T::Type{<:FloatFamily}) = Float64.(T.(arr))
trim_float(arr::Array{Complex{Float64}}, T::Type{<:FloatFamily}) = Complex{Float64}.(Complex{T}.(arr))

@testset "blas-level-1-asum-real" begin
  arr64 = rand(10) - 0.5
  arr16_0 = Posit{16,0}.(arr64)
  arr32_1 = Posit{32,1}.(arr64)

  #test exact asum with no striding against the julia native asum function.
  @test Posit{16,0}(BLAS.asum(4, trim_float(arr64, Posit{16,0}), 1)) == BLAS.asum(4, arr16_0, 1)
  @test Posit{32,1}(BLAS.asum(4, trim_float(arr64, Posit{32,1}), 1)) == BLAS.asum(4, arr32_1, 1)

  #test exact asum with striding
  @test Posit{16,0}(BLAS.asum(4, trim_float(arr64, Posit{16,0}), 2)) == BLAS.asum(4, arr16_0, 2)
  @test Posit{32,1}(BLAS.asum(4, trim_float(arr64, Posit{32,1}), 2)) == BLAS.asum(4, arr32_1, 2)
end

@testset "blas-level-1-asum-cplx" begin
  arr_c64 = (rand(10) - 0.5) + (rand(10) - 0.5) * im
  arr_c32_1 = Complex{Posit{32,1}}.(arr_c64)

  #currently, asum_cmplx does not meet exactitude requirements.  Instead, we'll show that it's better
  #than the equivalent basic floating point
  arr_c32 = Complex{Float32}.(arr_c64)

  delta_p32 = abs(
    BLAS.asum(4, trim_float(arr_c64, Posit{32,1}), 1) -
    Float64(BLAS.asum(4, arr_c32_1, 1)))

  delta_f32 = abs(
    BLAS.asum(4, trim_float(arr_c64, Float32), 1) -
    Float64(BLAS.asum(4, arr_c32_1, 1)))

  @test delta_p32 < delta_f32

  #test exact asum with no striding against the julia native asum function.
  #@test Posit{16,0}(BLAS.asum(4, trim_float(arr_c64, Posit{16,0}), 1)) == BLAS.asum(4, arr_c16_0, 1)
  #@test Posit{32,1}(BLAS.asum(4, trim_float(arr_c64, Posit{32,1}), 1)) == BLAS.asum(4, arr_c32_1, 1)

  #test exact asum with striding
  #@test Posit{16,0}(BLAS.asum(4, trim_float(arr_c64, Posit{16,0}), 2)) == BLAS.asum(4, arr_c16_0, 2)
  #@test Posit{32,1}(BLAS.asum(4, trim_float(arr_c64, Posit{32,1}), 2)) == BLAS.asum(4, arr_c32_1, 2)
end

@testset "blas-level-1-dot" begin
  arra_64 = rand(10) - 0.5
  arra_16_0 = Posit{16,0}.(arra_64)
  arra_32_1 = Posit{32,1}.(arra_64)

  arrb_64 = rand(10) - 0.5
  arrb_16_0 = Posit{16,0}.(arrb_64)
  arrb_32_1 = Posit{32,1}.(arrb_64)

  #test exact dot product as the basic dot product.
  @test Posit{16,0}(trim_float(arra_64, Posit{16,0}) ⋅ trim_float(arrb_64, Posit{16,0})) == arra_16_0 ⋅ arrb_16_0
  @test Posit{16,0}(trim_float(arra_64, Posit{16,0}) ⋅ trim_float(arrb_64, Posit{16,0})) == dot(arra_16_0, arrb_16_0)
  @test Posit{32,1}(trim_float(arra_64, Posit{32,1}) ⋅ trim_float(arrb_64, Posit{32,1})) == arra_32_1 ⋅ arrb_32_1

  #test exact asum with striding
  #@test Posit{16,0}(BLAS.dot(4, trim_float(arr64, Posit{16,0}), 2)) == BLAS.asum(4, arr16_0, 2)
  #@test Posit{32,1}(BLAS.dot(4, trim_float(arr64, Posit{32,1}), 2)) == BLAS.asum(4, arr32_1, 2)
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers