text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{String, Any}(
"physics" => "Advection",
"run_type" => 1,
"order" => 1,
"IC_name" => "ICsinwave",
#"IC_name" => "ICFile",
#"ICfname" => "q_vec1.dat_old",
"numBC" => 1,
#"BC1" => [4, 10, 7, 13],
"BC1" => [ 0 ],
"BC1_name" => "sinwaveBC",
#"BC2" => [7, 13],
#"BC2_name" => "noPenetrationBC",
#"BC2_name" => "Rho1E2U3BC",
"delta_t" => 0.005,
"t_max" => 500.000,
#"smb_name" => "../../mesh_files/channel.smb",
#"dmg_name" => "../../mesh_files/channel.dmg",
"smb_name" => "SRCMESHES/channel_8el_large.smb",
#"dmg_name" => ".null",
"res_abstol" => 1e-12,
"res_reltol" => -1.0,
"step_tol" => 1e-9,
"itermax" => 20,
"write_face_vertnums" => true,
"solve" => false
)


text:
using Lexicon, TargetedLearning

fmnames = [("tmles", TargetedLearning.TMLEs),
           ("ctmles", TargetedLearning.CTMLEs),
           #("strategies", TargetedLearning.CTMLEs.Strategies),
           ("common", TargetedLearning.Common),
           ("lreg", TargetedLearning.LReg),
           ("parameters", TargetedLearning.Parameters),
           #("qmodels", TargetedLearning.Qmodels)
           ]

for (fname, mod) in fmnames
	save(joinpath(Pkg.dir("TargetedLearning"), "docs", "api", "$fname.md"), mod, mathjax=true, include_internal = false)
end

nbpath = joinpath(Pkg.dir("TargetedLearning"), "docs", "user-guide", "lalonde_example.ipynb")
mdpath = joinpath(Pkg.dir("TargetedLearning"), "docs", "user-guide", "lalonde_example.md")

run(`ipython nbconvert --to markdown $nbpath --output=$mdpath`)

ext: jl


text: const SYMATA_VERSION = v"0.4.6"


text:
const ALLOWED_OPTIONS = [:printlevel, :axtol, :atytol, :objtol, :pinftol, :dinftol, :maxiter, :minstepfrac, :maxstepfrac, :minstepp, :minstepd, :usexzgap, :tweakgap, :affine, :perturbobj, :fastmode, :write_prob]

# The :write_prob option is for the following function
function write_prob(m)
    let wrt = string(get(m.options, :write_prob, ""))
        if length(wrt) > 0
            k = 1
            wrtf = "$wrt.$k"
            while isfile(wrtf)
                wrtf = "$wrt.$k"
                k += 1
            end
            @info "Writing problem to $(pwd())/$(wrtf)"
            write_prob(wrtf, m.C, m.b, map(A->A.csdp, m.As))
        end
    end
end

options(params::Dict{Symbol}) = get(params, :printlevel, 1), paramstruc(params)


text:
module Keys
### This file contains all functions which are necessary for accessing .peacevote/keys folder. In future it would be great to have a password prtotection so one could do easier backups.

using Pkg.TOML
using Base: UUID
using ..Types: ID, DemeID, keydir
using ..Plugins: Notary

struct Signer
    uuid::UUID
    id::ID
    sign::Function
end

function _save(s,fname)
    mkpath(dirname(fname))
    dict = Dict(s)
    open(fname, "w") do io
        TOML.print(io, dict)
    end
end

function Signer(uuid::UUID,notary::Notary,account::AbstractString)
    fname = keydir(uuid) * account

    if !isfile(fname)
        @info "Creating a new signer for the community"
        s = notary.Signer()
        _save(s,fname)
    end

    dict = TOML.parsefile(fname)
    signer = notary.Signer(dict)
    sign(data) = notary.Signature(data,signer)

    test = "A test message."
    signature = sign(test)
    id = notary.verify(test,signature)

    @assert id!=nothing

    return Signer(uuid,id,sign)
end

#Signer(uuid::UUID,account::AbstractString) = Signer(uuid,Notary(DemeSpec(uuid)),account)

### Perhaps I actually do not need this method !!!
#Signer(deme::Deme,account::AbstractString) = Signer(deme.spec.uuid,deme.notary,account)

### In ususal computer science
### certified authorithy is called also a certificate.
### However the therm voucher, verifier, issuer (tust anchor) would be more appropriate in such setting.
### On the other hand one might look on it as a certificate from the chain of trust.
### Anyway there should be a difference in the terminolgy.
struct Certificate{T}
    document::T
    signature::Dict
end

Certificate(x,signer::Signer) = Certificate(x,Dict(signer.sign("$x")))

function Dict(config::Certificate{T}) where T
    sdict = Dict(config.signature)
    dict = Dict(config.document)
    return Dict(dict...,"signature"=>sdict)
end

function Certificate{T}(dict::Dict) where T
    document = T(dict)
    signature = dict["signature"]
    return Certificate(document,signature)
end

### This is the type which requires a dynamic dyspatch
struct Envelope{T}
    uuid::UUID
    contract::Certificate{T} # Perhaps Union{Contract{T},Certificate{T}}
end

Envelope(data,signer::Signer) = Envelope(signer.uuid,Certificate(data,signer))

verify(cert::Certificate,notary::Notary) = notary.verify("$(cert.document)",notary.Signature(cert.signature))

function verify(envelope::Envelope)
    demespec = DemeSpec(envelope.uuid)
    notary = Notary(demespec)
    id = verify(envelope.contract,notary)
    return DemeID(envelope.uuid,id)
end

struct Intent{T}
    document::T
    reference::Union{ID,DemeID}
end

function Intent(envelope::Envelope)
    demespec = DemeSpec(envelope.uuid)
    notary = Notary(demespec)

    contract = envelope.contract
    id = notary.verify("$(contract.document)",notary.Signature(contract.signature))

    ref = DemeID(envelope.uuid,id)
    return Intent(document,ref)
end

function Intent(contract::Certificate,notary::Notary)
    signature = notary.Signature(contract.signature)
    id = notary.verify("$(contract.document)",signature)
    return Intent(contract.document,id)
end

struct Contract{T}
    document::T
    signatures::Vector{Dict{String,Any}}
end

function Dict(contract::Contract)
    dict = Dict(contract.document)

    signatures = Dict[]
    for s in contract.signatures
        push!(signatures,Dict(s))
    end
    dict["signatures"] = signatures

    return dict
end

function Contract{T}(dict::Dict) where T
    braid = T(dict)
    signatures = dict["signatures"]

    return Contract(braid,signatures)
end

### This one is associated to the contract
# Gives permission to do something, like add a new key to the braidchain synchronically.
# Common consent
struct Consensus{T}
    document::T
    references::Vector{ID} ### The contract is allowed to be made only between equals. That is necessary to resolve disputes as the Deme who made the contract shall be held responsable for it.
end

function Consensus(contract::Contract,notary::Notary)
    doc = "$(contract.document)"
    refs = ID[]
    for s in contract.signatures
        signature = notary.Signature(s)
        id = notary.verify(doc,signature)
        push!(refs,id)
    end
    return Consensus(contract.document,refs)
end

end


text:
# Autogenerated wrapper script for DirectQhullHelper_jll for x86_64-unknown-freebsd
export libDirectQhullHelper

using Qhull_jll
JLLWrappers.@generate_wrapper_header("DirectQhullHelper")
JLLWrappers.@declare_library_product(libDirectQhullHelper, "libDirectQhullHelper.so")
function __init__()
    JLLWrappers.@generate_init_header(Qhull_jll)
    JLLWrappers.@init_library_product(
        libDirectQhullHelper,
        "lib/libDirectQhullHelper.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Flase, Test

world = World(
    v0 = 1.,
    n_dogs = 120,
    boxsize = 10.0,
    motion = BrownianMotion(
        noise = 0.5,
        friction = 1.0
        ),
    sheeps = DenseSheeps(
        10,
        n_sheeps = 10,
        )
    )
@testset "Test iteration" begin
    @test iterate( world.dogs )[2][2] == 1
    @test iterate( world.dogs, iterate(world.dogs)[2] )[2][2] == 2
end # testset


text:
@testset "Test `-`" begin
    @test typeof(-EngineeringStrain(1:6)) == EngineeringStrain{Int}
    @test typeof(-EngineeringStrain(1.0:6.0)) == EngineeringStrain{Float64}
    @test -EngineeringStrain(1:6) == EngineeringStrain(-1, -2, -3, -4, -5, -6)
    @test -EngineeringStress(1.0:6.0) == EngineeringStress(-1.0, -2, -3, -4, -5, -6)
    @test -(-EngineeringStress(1.0:6.0)) == EngineeringStress(1.0:6.0)
end

@testset "Test `*` and `/`" begin
    a = EngineeringStrain(1:6)
    @test typeof(2a) == EngineeringStrain{Int}
    @test typeof(2.0a) == EngineeringStrain{Float64}
    @test 2a == 2 * a == a * 2 == 2.0a == EngineeringStress(2:2:12)
    @test a / 2 == EngineeringStress(0.5, 1, 1.5, 2, 2.5, 3)
    @test 2a / 2 == a
end

@testset "Test addition and subtraction" begin
    a = EngineeringStrain(1:6)
    b = -EngineeringStrain(1:1:6)
    @test a - b == 2a
    @test b - a == -2a
    @test a + b == EngineeringStrain(zeros(6))
    a = TensorStress(EngineeringStress(1:6))
    b = -TensorStress(EngineeringStress(1:6))
    @test a - b == 2a
    @test b - a == -2a
    @test a + b == TensorStress(zeros(3, 3))
end


text:
include("calcmetrics.jl")
include("boolmetrics.jl")
include("metricWrapper.jl")
include("fairness.jl")


text:
# calculate probability by pdf
# sensor_data_200.txt
# plot as 2d-contour

module ContourPdf200
    using DataFrames, CSV, Statistics, LinearAlgebra
    using Base.Iterators, Plots
    pyplot()

    function pdf(x, mu, cov_mat)
        # determinant
        det_mat = det(cov_mat)

        diff_x = x - mu
        return exp(-((diff_x'/cov_mat)*diff_x)/2)/(2*pi*sqrt(det_mat))
    end

    function main(is_test=false)
        # input
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df = CSV.read(data_path, DataFrame,
                      header=["date", "time", "ir", "lidar"],
                      delim=' ')

        # variance
        var_ir = var(df.ir)
        var_lidar = var(df.lidar)

        #covariance
        diff_ir = df.ir .- mean(df.ir)
        diff_lidar = df.lidar .- mean(df.lidar)
        a = diff_ir .* diff_lidar
        covar = cov(df.ir, df.lidar, corrected=false)

        # covariance matrix
        cov_mat = [var_ir covar; covar var_lidar]

        # mean(mu)
        mu = [mean(df.ir); mean(df.lidar)]

        # plot contour
        vx = 280:340
        vy = 190:230
        z = [pdf([x; y], mu, cov_mat) for x in vx, y in vy]
        contour(vx, vy, z', label="contour", c=:haline, aspect_ratio=:equal)

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/multi_dim_gauss_dist/contour_pdf/contour_pdf_200.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# THIS FILE ADAPTED FROM GraphLayout.jl
# see https://github.com/IainNZ/GraphLayout.jl

"""
    Use the spring/repulsion model of Fruchterman and Reingold (1991):
        Attractive force:  f_a(d) =  d^2 / k
        Repulsive force:  f_r(d) = -k^2 / d
    where d is distance between two vertices and the optimal distance
    between vertices k is defined as C * sqrt( area / num_vertices )
    where C is a parameter we can adjust

    Arguments:
    adj_matrix Adjacency matrix of some type. Non-zero of the eltype
               of the matrix is used to determine if a link exists,
               but currently no sense of magnitude
    C          Constant to fiddle with density of resulting layout
    MAXITER    Number of iterations we apply the forces
    INITTEMP   Initial "temperature", controls movement per iteration
"""
function layout_spring_adj(
    adj_matrix::Array{T,2};
    C = 2.0,
    MAXITER = 100,
    INITTEMP = 2.0,
) where {T}

    size(adj_matrix, 1) != size(adj_matrix, 2) && error("Adj. matrix must be square.")
    N = size(adj_matrix, 1)

    # Initial layout is random on the square [-1,+1]^2
    locs_x = 2 * rand(N) .- 1.0
    locs_y = 2 * rand(N) .- 1.0

    # The optimal distance bewteen vertices
    K = C * sqrt(4.0 / N)

    # Store forces and apply at end of iteration all at once
    force_x = zeros(N)
    force_y = zeros(N)

    # Iterate MAXITER times
    @inbounds for iter = 1:MAXITER
        # Calculate forces
        for i = 1:N
            force_vec_x = 0.0
            force_vec_y = 0.0
            for j = 1:N
                i == j && continue
                d_x = locs_x[j] - locs_x[i]
                d_y = locs_y[j] - locs_y[i]
                d = sqrt(d_x^2 + d_y^2)
                if adj_matrix[i, j] != zero(eltype(adj_matrix)) ||
                   adj_matrix[j, i] != zero(eltype(adj_matrix))
                    # F = d^2 / K - K^2 / d
                    F_d = d / K - K^2 / d^2
                else
                    # Just repulsive
                    # F = -K^2 / d^
                    F_d = -K^2 / d^2
                end
                # d  /          sin θ = d_y/d = fy/F
                # F /| dy fy    -> fy = F*d_y/d
                #  / |          cos θ = d_x/d = fx/F
                # /---          -> fx = F*d_x/d
                # dx fx
                force_vec_x += F_d * d_x
                force_vec_y += F_d * d_y
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
        end
        # Cool down
        TEMP = INITTEMP / iter
        # Now apply them, but limit to temperature
        for i = 1:N
            force_mag = sqrt(force_x[i]^2 + force_y[i]^2)
            scale = min(force_mag, TEMP) / force_mag
            locs_x[i] += force_x[i] * scale
            #locs_x[i]  = max(-1.0, min(locs_x[i], +1.0))
            locs_y[i] += force_y[i] * scale
            #locs_y[i]  = max(-1.0, min(locs_y[i], +1.0))
        end
    end

    # Scale to unit square
    min_x, max_x = minimum(locs_x), maximum(locs_x)
    min_y, max_y = minimum(locs_y), maximum(locs_y)
    function scaler(z, a, b)
        2.0 * ((z - a) / (b - a)) - 1.0
    end
    locs_x = map(z -> scaler(z, min_x, max_x), locs_x)
    locs_y = map(z -> scaler(z, min_y, max_y), locs_y)

    return locs_x, locs_y
end

function _spring(G::SimpleGraph{T}, nits::Int = 1000) where {T}
    n = NV(G)
    A, vv = private_adj(G)

    d = Dict{T,Vector{Float64}}()

    if n == 0
        return d
    end

    if n == 1
        v = first(vv)
        d[v] = [0.0, 0.0]
        return d
    end

    x, y = layout_spring_adj(A, MAXITER = nits)

    for i = 1:n
        v = vv[i]
        d[v] = [x[i], y[i]]
    end
    return d
end


text:
# Julia wrapper for header: /usr/include/scip/scip_cons.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlr(scip, name, desc, sepapriority, enfopriority, chckpriority, sepafreq, propfreq, eagerfreq, maxprerounds, delaysepa, delayprop, needscons, proptiming, presoltiming, conshdlrcopy, consfree, consinit, consexit, consinitpre, consexitpre, consinitsol, consexitsol, consdelete, constrans, consinitlp, conssepalp, conssepasol, consenfolp, consenforelax, consenfops, conscheck, consprop, conspresol, consresprop, conslock, consactive, consdeactive, consenable, consdisable, consdelvars, consprint, conscopy, consparse, consgetvars, consgetnvars, consgetdivebdchgs, conshdlrdata)
    ccall((:SCIPincludeConshdlr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, Cint, Cint, Cint, Cint, Cint, UInt32, UInt32, UInt32, SCIP_PROPTIMING, SCIP_PRESOLTIMING, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_CONSHDLRDATA}), scip, name, desc, sepapriority, enfopriority, chckpriority, sepafreq, propfreq, eagerfreq, maxprerounds, delaysepa, delayprop, needscons, proptiming, presoltiming, conshdlrcopy, consfree, consinit, consexit, consinitpre, consexitpre, consinitsol, consexitsol, consdelete, constrans, consinitlp, conssepalp, conssepasol, consenfolp, consenforelax, consenfops, conscheck, consprop, conspresol, consresprop, conslock, consactive, consdeactive, consenable, consdisable, consdelvars, consprint, conscopy, consparse, consgetvars, consgetnvars, consgetdivebdchgs, conshdlrdata)
end

function SCIPincludeConshdlrBasic(scip, conshdlrptr, name, desc, enfopriority, chckpriority, eagerfreq, needscons, consenfolp, consenfops, conscheck, conslock, conshdlrdata)
    ccall((:SCIPincludeConshdlrBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONSHDLR}}, Cstring, Cstring, Cint, Cint, Cint, UInt32, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_CONSHDLRDATA}), scip, conshdlrptr, name, desc, enfopriority, chckpriority, eagerfreq, needscons, consenfolp, consenfops, conscheck, conslock, conshdlrdata)
end

function SCIPsetConshdlrSepa(scip, conshdlr, conssepalp, conssepasol, sepafreq, sepapriority, delaysepa)
    ccall((:SCIPsetConshdlrSepa, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}, Ptr{Cvoid}, Cint, Cint, UInt32), scip, conshdlr, conssepalp, conssepasol, sepafreq, sepapriority, delaysepa)
end

function SCIPsetConshdlrProp(scip, conshdlr, consprop, propfreq, delayprop, proptiming)
    ccall((:SCIPsetConshdlrProp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}, Cint, UInt32, SCIP_PROPTIMING), scip, conshdlr, consprop, propfreq, delayprop, proptiming)
end

function SCIPsetConshdlrEnforelax(scip, conshdlr, consenforelax)
    ccall((:SCIPsetConshdlrEnforelax, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consenforelax)
end

function SCIPsetConshdlrCopy(scip, conshdlr, conshdlrcopy, conscopy)
    ccall((:SCIPsetConshdlrCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}, Ptr{Cvoid}), scip, conshdlr, conshdlrcopy, conscopy)
end

function SCIPsetConshdlrFree(scip, conshdlr, consfree)
    ccall((:SCIPsetConshdlrFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consfree)
end

function SCIPsetConshdlrInit(scip, conshdlr, consinit)
    ccall((:SCIPsetConshdlrInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consinit)
end

function SCIPsetConshdlrExit(scip, conshdlr, consexit)
    ccall((:SCIPsetConshdlrExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consexit)
end

function SCIPsetConshdlrInitsol(scip, conshdlr, consinitsol)
    ccall((:SCIPsetConshdlrInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consinitsol)
end

function SCIPsetConshdlrExitsol(scip, conshdlr, consexitsol)
    ccall((:SCIPsetConshdlrExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consexitsol)
end

function SCIPsetConshdlrInitpre(scip, conshdlr, consinitpre)
    ccall((:SCIPsetConshdlrInitpre, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consinitpre)
end

function SCIPsetConshdlrExitpre(scip, conshdlr, consexitpre)
    ccall((:SCIPsetConshdlrExitpre, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consexitpre)
end

function SCIPsetConshdlrPresol(scip, conshdlr, conspresol, maxprerounds, presoltiming)
    ccall((:SCIPsetConshdlrPresol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}, Cint, SCIP_PRESOLTIMING), scip, conshdlr, conspresol, maxprerounds, presoltiming)
end

function SCIPsetConshdlrDelete(scip, conshdlr, consdelete)
    ccall((:SCIPsetConshdlrDelete, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consdelete)
end

function SCIPsetConshdlrTrans(scip, conshdlr, constrans)
    ccall((:SCIPsetConshdlrTrans, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, constrans)
end

function SCIPsetConshdlrInitlp(scip, conshdlr, consinitlp)
    ccall((:SCIPsetConshdlrInitlp, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consinitlp)
end

function SCIPsetConshdlrResprop(scip, conshdlr, consresprop)
    ccall((:SCIPsetConshdlrResprop, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consresprop)
end

function SCIPsetConshdlrActive(scip, conshdlr, consactive)
    ccall((:SCIPsetConshdlrActive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consactive)
end

function SCIPsetConshdlrDeactive(scip, conshdlr, consdeactive)
    ccall((:SCIPsetConshdlrDeactive, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consdeactive)
end

function SCIPsetConshdlrEnable(scip, conshdlr, consenable)
    ccall((:SCIPsetConshdlrEnable, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consenable)
end

function SCIPsetConshdlrDisable(scip, conshdlr, consdisable)
    ccall((:SCIPsetConshdlrDisable, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consdisable)
end

function SCIPsetConshdlrDelvars(scip, conshdlr, consdelvars)
    ccall((:SCIPsetConshdlrDelvars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consdelvars)
end

function SCIPsetConshdlrPrint(scip, conshdlr, consprint)
    ccall((:SCIPsetConshdlrPrint, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consprint)
end

function SCIPsetConshdlrParse(scip, conshdlr, consparse)
    ccall((:SCIPsetConshdlrParse, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consparse)
end

function SCIPsetConshdlrGetVars(scip, conshdlr, consgetvars)
    ccall((:SCIPsetConshdlrGetVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consgetvars)
end

function SCIPsetConshdlrGetNVars(scip, conshdlr, consgetnvars)
    ccall((:SCIPsetConshdlrGetNVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consgetnvars)
end

function SCIPsetConshdlrGetDiveBdChgs(scip, conshdlr, consgetdivebdchgs)
    ccall((:SCIPsetConshdlrGetDiveBdChgs, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONSHDLR}, Ptr{Cvoid}), scip, conshdlr, consgetdivebdchgs)
end

function SCIPfindConshdlr(scip, name)
    ccall((:SCIPfindConshdlr, libscip), Ptr{SCIP_CONSHDLR}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetConshdlrs(scip)
    ccall((:SCIPgetConshdlrs, libscip), Ptr{Ptr{SCIP_CONSHDLR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNConshdlrs(scip)
    ccall((:SCIPgetNConshdlrs, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{SCIP_CONSHDLR}, Ptr{SCIP_CONSDATA}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPparseCons(scip, cons, str, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode, success)
    ccall((:SCIPparseCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, Ptr{UInt32}), scip, cons, str, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode, success)
end

function SCIPcaptureCons(scip, cons)
    ccall((:SCIPcaptureCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPreleaseCons(scip, cons)
    ccall((:SCIPreleaseCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}), scip, cons)
end

function SCIPchgConsName(scip, cons, name)
    ccall((:SCIPchgConsName, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cstring), scip, cons, name)
end

function SCIPsetConsInitial(scip, cons, initial)
    ccall((:SCIPsetConsInitial, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, initial)
end

function SCIPsetConsSeparated(scip, cons, separate)
    ccall((:SCIPsetConsSeparated, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, separate)
end

function SCIPsetConsEnforced(scip, cons, enforce)
    ccall((:SCIPsetConsEnforced, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, enforce)
end

function SCIPsetConsChecked(scip, cons, check)
    ccall((:SCIPsetConsChecked, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, check)
end

function SCIPsetConsPropagated(scip, cons, propagate)
    ccall((:SCIPsetConsPropagated, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, propagate)
end

function SCIPsetConsLocal(scip, cons, _local)
    ccall((:SCIPsetConsLocal, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, _local)
end

function SCIPsetConsModifiable(scip, cons, modifiable)
    ccall((:SCIPsetConsModifiable, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, modifiable)
end

function SCIPsetConsDynamic(scip, cons, dynamic)
    ccall((:SCIPsetConsDynamic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, dynamic)
end

function SCIPsetConsRemovable(scip, cons, removable)
    ccall((:SCIPsetConsRemovable, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, removable)
end

function SCIPsetConsStickingAtNode(scip, cons, stickingatnode)
    ccall((:SCIPsetConsStickingAtNode, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32), scip, cons, stickingatnode)
end

function SCIPupdateConsFlags(scip, cons0, cons1)
    ccall((:SCIPupdateConsFlags, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_CONS}), scip, cons0, cons1)
end

function SCIPtransformCons(scip, cons, transcons)
    ccall((:SCIPtransformCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_CONS}}), scip, cons, transcons)
end

function SCIPtransformConss(scip, nconss, conss, transconss)
    ccall((:SCIPtransformConss, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_CONS}}, Ptr{Ptr{SCIP_CONS}}), scip, nconss, conss, transconss)
end

function SCIPgetTransformedCons(scip, cons, transcons)
    ccall((:SCIPgetTransformedCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_CONS}}), scip, cons, transcons)
end

function SCIPgetTransformedConss(scip, nconss, conss, transconss)
    ccall((:SCIPgetTransformedConss, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{Ptr{SCIP_CONS}}, Ptr{Ptr{SCIP_CONS}}), scip, nconss, conss, transconss)
end

function SCIPaddConsAge(scip, cons, deltaage)
    ccall((:SCIPaddConsAge, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cdouble), scip, cons, deltaage)
end

function SCIPincConsAge(scip, cons)
    ccall((:SCIPincConsAge, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPresetConsAge(scip, cons)
    ccall((:SCIPresetConsAge, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPenableCons(scip, cons)
    ccall((:SCIPenableCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPdisableCons(scip, cons)
    ccall((:SCIPdisableCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPenableConsSeparation(scip, cons)
    ccall((:SCIPenableConsSeparation, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPdisableConsSeparation(scip, cons)
    ccall((:SCIPdisableConsSeparation, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPenableConsPropagation(scip, cons)
    ccall((:SCIPenableConsPropagation, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPdisableConsPropagation(scip, cons)
    ccall((:SCIPdisableConsPropagation, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPmarkConsPropagate(scip, cons)
    ccall((:SCIPmarkConsPropagate, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPunmarkConsPropagate(scip, cons)
    ccall((:SCIPunmarkConsPropagate, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPaddConsLocksType(scip, cons, locktype, nlockspos, nlocksneg)
    ccall((:SCIPaddConsLocksType, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, SCIP_LOCKTYPE, Cint, Cint), scip, cons, locktype, nlockspos, nlocksneg)
end

function SCIPaddConsLocks(scip, cons, nlockspos, nlocksneg)
    ccall((:SCIPaddConsLocks, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cint, Cint), scip, cons, nlockspos, nlocksneg)
end

function SCIPcheckCons(scip, cons, sol, checkintegrality, checklprows, printreason, result)
    ccall((:SCIPcheckCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, UInt32, UInt32, UInt32, Ptr{SCIP_RESULT}), scip, cons, sol, checkintegrality, checklprows, printreason, result)
end

function SCIPenfopsCons(scip, cons, solinfeasible, objinfeasible, result)
    ccall((:SCIPenfopsCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32, UInt32, Ptr{SCIP_RESULT}), scip, cons, solinfeasible, objinfeasible, result)
end

function SCIPenfolpCons(scip, cons, solinfeasible, result)
    ccall((:SCIPenfolpCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, UInt32, Ptr{SCIP_RESULT}), scip, cons, solinfeasible, result)
end

function SCIPenforelaxCons(scip, cons, sol, solinfeasible, result)
    ccall((:SCIPenforelaxCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, UInt32, Ptr{SCIP_RESULT}), scip, cons, sol, solinfeasible, result)
end

function SCIPinitlpCons(scip, cons, infeasible)
    ccall((:SCIPinitlpCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{UInt32}), scip, cons, infeasible)
end

function SCIPsepalpCons(scip, cons, result)
    ccall((:SCIPsepalpCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_RESULT}), scip, cons, result)
end

function SCIPsepasolCons(scip, cons, sol, result)
    ccall((:SCIPsepasolCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_SOL}, Ptr{SCIP_RESULT}), scip, cons, sol, result)
end

function SCIPpropCons(scip, cons, proptiming, result)
    ccall((:SCIPpropCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, SCIP_PROPTIMING, Ptr{SCIP_RESULT}), scip, cons, proptiming, result)
end

function SCIPrespropCons(scip, cons, infervar, inferinfo, boundtype, bdchgidx, relaxedbd, result)
    ccall((:SCIPrespropCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}, Cint, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}, Cdouble, Ptr{SCIP_RESULT}), scip, cons, infervar, inferinfo, boundtype, bdchgidx, relaxedbd, result)
end

function SCIPpresolCons(scip, cons, nrounds, presoltiming, nnewfixedvars, nnewaggrvars, nnewchgvartypes, nnewchgbds, nnewholes, nnewdelconss, nnewaddconss, nnewupgdconss, nnewchgcoefs, nnewchgsides, nfixedvars, naggrvars, nchgvartypes, nchgbds, naddholes, ndelconss, naddconss, nupgdconss, nchgcoefs, nchgsides, result)
    ccall((:SCIPpresolCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Cint, SCIP_PRESOLTIMING, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{SCIP_RESULT}), scip, cons, nrounds, presoltiming, nnewfixedvars, nnewaggrvars, nnewchgvartypes, nnewchgbds, nnewholes, nnewdelconss, nnewaddconss, nnewupgdconss, nnewchgcoefs, nnewchgsides, nfixedvars, naggrvars, nchgvartypes, nchgbds, naddholes, ndelconss, naddconss, nupgdconss, nchgcoefs, nchgsides, result)
end

function SCIPactiveCons(scip, cons)
    ccall((:SCIPactiveCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPdeactiveCons(scip, cons)
    ccall((:SCIPdeactiveCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPprintCons(scip, cons, file)
    ccall((:SCIPprintCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{FILE}), scip, cons, file)
end

function SCIPgetConsVars(scip, cons, vars, varssize, success)
    ccall((:SCIPgetConsVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Ptr{SCIP_VAR}}, Cint, Ptr{UInt32}), scip, cons, vars, varssize, success)
end

function SCIPgetConsNVars(scip, cons, nvars, success)
    ccall((:SCIPgetConsNVars, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{Cint}, Ptr{UInt32}), scip, cons, nvars, success)
end


text:
#example of creating custom types

importall Equations

type Exp<:Component
	x
end
simplify(c::Exp)=isa(simplify(c.x),Number)?exp(simplify(c.x)):c
type Cos<:Component
	x
end
simplify(c::Cos)=isa(simplify(c.x),Number)?cos(simplify(c.x)):c
type Sin<:Component
	x
end
simplify(c::Sin)=isa(simplify(c.x),Number)?sin(simplify(c.x)):c
e1=Vec([1,0,0]);e2=Vec([0,1,0]);e3=Vec([0,0,1])

type Differentiable <: Component
	fun
	der::Dict
end
#getargs(d::Differentiable)=(d.fun,)
Equations.maketype(d::Differentiable,fun)=Differentiable(fun(d.fun),d.der)
Differentiable(a,b)=Differentiable(a,[:x=>b])
ed=Vec(Differentiable(1,0),Differentiable(0,1),0)
type DerOp #<: Operator
	dx
end
DerOp()=DerOp(:x)
*(d::DerOp,f::Differentiable)=haskey(f.der,d.dx)?f.der[d.dx]:Der(f,d.dx)
*(d::DerOp,ex::Factor)=Der(ex,d.dx)
type Grad <: Operator
	v
end
Grad()=Grad(Vec([DerOp(:x),DerOp(:y),DerOp(:z)]))
Base.print(io::IO,g::Grad)=print(io,"∇")
function dot(v1::Grad,v2::Differentiable)
	ex=Expression(Term[Factor[(v1.v[1]*v2)[1]]])
	for t in 2:length(v1.v)
		push!(ex,Factor[(v1.v[t]*v2)[t]])
	end
	return simplify(ex)
end

eq1=Cross(:j,:B)≖Grad()*:p
eq2=Cross(Grad(),:B)≖:μ0*:j
eq2p=Dot(Grad(),Cross(Grad(),:B))≖:μ0*Dot(Grad(),:j)
eq3=Dot(Grad(),:B)≖0

eqc=eq1≖0


text:
module TimeSeriesUtils
#
# This module contains some stuff for modelling financial years, quarters and the like. Used for calculating
# pension ages and in the uprating routines.
#
using Base: Integer
using TimeSeries, Dates

export fy, fy_array, get_birthdate, fyear, fy_from_date, fy_from_bits

"""
For Unit testing. A date age_in_years + 1 month before today.
"""
function get_birthdate( age_in_years :: Integer, from_date :: DateTime = now() ) :: Date
    return Date(from_date - Year( age_in_years ) - Month(1))
end

"""
financial year as an interval
"""
function fy(year::Integer)::StepRange
    return Date(year,04,06):Day(1):Date(year+1,04,05)
end

"""

"""
function fyear( year :: Integer  ) :: Date
    return Date( year, 04, 06 )
end

function fy_from_date( d :: Date ) :: Integer
    y = year( d )
    if (month( d ) >= 4) && (day( d ) >= 6)
        y += 1
    end
    return y
end

function fy_from_bits( year :: Integer, month :: Integer, day :: Integer = 1) :: Integer
    if month >= 4 && day >= 6
        return year + 1
    end
    return year
end

function fy_array( years::UnitRange ) :: Vector{Date}
    n = size(years)[1]
    d = Vector{Date}(undef, n)
    i = 0
    for y in years
        i+=1
        d[i] = Date(y,4,6)
    end
    return d
end
# FIXME the day needn't be the 6th; see for example IFS fiscal facts
const FY_2000 = fyear(2000)
const FY_2001 = fyear(2001)
const FY_2002 = fyear(2002)
const FY_2003 = fyear(2003)
const FY_2004 = fyear(2004)
const FY_2005 = fyear(2005)
const FY_2006 = fyear(2006)
const FY_2007 = fyear(2007)
const FY_2008 = fyear(2008)
const FY_2009 = fyear(2009)
const FY_2010 = fyear(2010)
const FY_2011 = fyear(2011)
const FY_2012 = fyear(2012)
const FY_2013 = fyear(2013)
const FY_2014 = fyear(2014)
const FY_2015 = fyear(2015)
const FY_2016 = fyear(2016)
const FY_2017 = fyear(2017)
const FY_2018 = fyear(2018)
const FY_2019 = fyear(2019)
const FY_2020 = fyear(2020)
const FY_2021 = fyear(2021)
const FY_2022 = fyear(2022)
const FY_2023 = fyear(2023)
const FY_2024 = fyear(2024)
const FY_2025 = fyear(2025)
const FY_2026 = fyear(2026)
const FY_2027 = fyear(2027)
const FY_2028 = fyear(2028)
const FY_2029 = fyear(2029)
const FY_2030 = fyear(2030)

# financial years as intervals - y in FYI_2021 and so on
export FY_2000
export FY_2001
export FY_2002
export FY_2003
export FY_2004
export FY_2005
export FY_2006
export FY_2007
export FY_2008
export FY_2009
export FY_2010
export FY_2011
export FY_2012
export FY_2013
export FY_2014
export FY_2015
export FY_2016
export FY_2017
export FY_2018
export FY_2019
export FY_2020
export FY_2021
export FY_2022
export FY_2023
export FY_2024
export FY_2025
export FY_2026
export FY_2027
export FY_2028
export FY_2029
export FY_2030

# financial years as dates d < FY_2012 and so on
const FYI_2000 = fy(2000)
const FYI_2001 = fy(2001)
const FYI_2002 = fy(2002)
const FYI_2003 = fy(2003)
const FYI_2004 = fy(2004)
const FYI_2005 = fy(2005)
const FYI_2006 = fy(2006)
const FYI_2007 = fy(2007)
const FYI_2008 = fy(2008)
const FYI_2009 = fy(2009)
const FYI_2010 = fy(2010)
const FYI_2011 = fy(2011)
const FYI_2012 = fy(2012)
const FYI_2013 = fy(2013)
const FYI_2014 = fy(2014)
const FYI_2015 = fy(2015)
const FYI_2016 = fy(2016)
const FYI_2017 = fy(2017)
const FYI_2018 = fy(2018)
const FYI_2019 = fy(2019)
const FYI_2020 = fy(2020)
const FYI_2021 = fy(2021)
const FYI_2022 = fy(2022)
const FYI_2023 = fy(2023)
const FYI_2024 = fy(2024)
const FYI_2025 = fy(2025)
const FYI_2026 = fy(2026)
const FYI_2027 = fy(2027)
const FYI_2028 = fy(2028)
const FYI_2029 = fy(2029)
const FYI_2030 = fy(2030)

# financial years as intervals - y in FYI_2021 and so on
export FYI_2000
export FYI_2001
export FYI_2002
export FYI_2003
export FYI_2004
export FYI_2005
export FYI_2006
export FYI_2007
export FYI_2008
export FYI_2009
export FYI_2010
export FYI_2011
export FYI_2012
export FYI_2013
export FYI_2014
export FYI_2015
export FYI_2016
export FYI_2017
export FYI_2018
export FYI_2019
export FYI_2020
export FYI_2021
export FYI_2022
export FYI_2023
export FYI_2024
export FYI_2025
export FYI_2026
export FYI_2027
export FYI_2028
export FYI_2029
export FYI_2030

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module test_h3_lib_unidirectional_edges

using H3.Lib

Lib.h3IndexesAreNeighbors
Lib.getH3UnidirectionalEdge
Lib.h3UnidirectionalEdgeIsValid
Lib.getOriginH3IndexFromUnidirectionalEdge
Lib.getDestinationH3IndexFromUnidirectionalEdge
Lib.getH3IndexesFromUnidirectionalEdge
Lib.getH3UnidirectionalEdgesFromHexagon
Lib.getH3UnidirectionalEdgeBoundary

end # module test_h3_lib_unidirectional_edges


text:
@testset "minimal.absRevBin_test" begin

    @test abs_rev(emptyinterval(), entireinterval())[2] === emptyinterval()

    @test abs_rev(interval(0.0, 1.0), emptyinterval())[2] === emptyinterval()

    @test abs_rev(interval(0.0, 1.0), interval(7.0, 9.0))[2] === emptyinterval()

    @test abs_rev(emptyinterval(), interval(0.0, 1.0))[2] === emptyinterval()

    @test abs_rev(interval(-2.0, -1.0), entireinterval())[2] === emptyinterval()

    @test abs_rev(interval(1.0, 1.0), entireinterval())[2] === Interval(-1.0, 1.0)

    @test abs_rev(interval(0.0, 0.0), entireinterval())[2] === Interval(0.0, 0.0)

    @test abs_rev(interval(-1.0, -1.0), entireinterval())[2] === emptyinterval()

    @test abs_rev(interval(0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023), entireinterval())[2] === Interval(-0x1.FFFFFFFFFFFFFp1023, 0x1.FFFFFFFFFFFFFp1023)

    @test abs_rev(interval(0x1p-1022, 0x1p-1022), entireinterval())[2] === Interval(-0x1p-1022, 0x1p-1022)

    @test abs_rev(interval(-0x1p-1022, -0x1p-1022), entireinterval())[2] === emptyinterval()

    @test abs_rev(interval(-0x1.FFFFFFFFFFFFFp1023, -0x1.FFFFFFFFFFFFFp1023), entireinterval())[2] === emptyinterval()

    @test abs_rev(interval(1.0, 2.0), entireinterval())[2] === Interval(-2.0, 2.0)

    @test abs_rev(interval(1.0, 2.0), interval(0.0, 2.0))[2] === Interval(1.0, 2.0)

    @test abs_rev(interval(0.0, 1.0), interval(-0.5, 2.0))[2] === Interval(-0.5, 1.0)

    @test abs_rev(interval(-1.0, 1.0), entireinterval())[2] === Interval(-1.0, 1.0)

    @test abs_rev(interval(-1.0, 0.0), entireinterval())[2] === Interval(0.0, 0.0)

    @test abs_rev(interval(0.0, Inf), entireinterval())[2] === entireinterval()

    @test abs_rev(entireinterval(), entireinterval())[2] === entireinterval()

    @test abs_rev(interval(-Inf, 0.0), entireinterval())[2] === Interval(0.0, 0.0)

    @test abs_rev(interval(1.0, Inf), interval(-Inf, 0.0))[2] === Interval(-Inf, -1.0)

    @test abs_rev(interval(-1.0, Inf), entireinterval())[2] === entireinterval()

    @test abs_rev(interval(-Inf, -1.0), entireinterval())[2] === emptyinterval()

    @test abs_rev(interval(-Inf, 1.0), entireinterval())[2] === Interval(-1.0, 1.0)

end


text: Dict{String, Any}("1" => Dict{String, Any}("value" => "1", "type" => "integer"))


text:
module Encoder

import ...AbstractCondition

# Splat fields
splat(x, idx=1:fieldcount(typeof(x))) = (getfield(x, i) for i ∈ idx)

"""
    Enc()

Wrap an `IOBuffer` to hold the outbound message
"""
struct Enc
  buf::IOBuffer
end
Enc() = Enc(IOBuffer(sizehint=64))

"""
    (e::Enc)(::T)

Define various encodings for known types
"""
(e::Enc)(::T) where T = error("Unknown Type: $T")

(e::Enc)(x::Union{AbstractString,Int,Float64,Symbol}) = print(e.buf, x, '\0')

(e::Enc)(::Nothing) = e("")

(e::Enc)(x::Bool) = e(x ? "1" : "0")

# Enums
(e::Enc)(x::Enum{Int32}) = e(Int(x))

# NamedTuples
function (e::Enc)(x::NamedTuple)

  for (n, v) ∈ pairs(x)
    v isa Union{AbstractString,Int,Float64} ||
      @warn "Unsupported Type in NamedTuple" N=n V=v T=typeof(v)

    print(e.buf, n, '=', v, ';')
  end

  print(e.buf, '\0')
end

# Condition Types
(e::Enc)(x::AbstractCondition{E}) where E = e(E, splat(x))

# Generators, as returned by splat()
(e::Enc)(x::Base.Generator) = foreach(e, x)

# Multiple arguments
function (e::Enc)(x, y...)
  e(x)
  foreach(e, y)
end

end


text:
# Copyright 2017,2018 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test

include("loaddata.jl")

@testset "encode by farthest points" begin
    using KernelMethods
    using SimilaritySearch: l2_distance, normalize!

    X, y = loadiris()
    dist = l2_distance
    # criterion = change_criterion(0.01)
    refs = Vector{typeof(X[1])}()
    dmax = 0.0
    function callback(c, _dmax)
        push!(refs, X[c])
        dmax = _dmax
    end

    fftraversal(callback, dist, X, sqrt_criterion())
    g = cauchy_kernel(dist, dmax/2)
    M = kmap(X, g, refs)
    nnc = NearNeighborClassifier(l2_distance, M, y)

    @test optimize!(nnc, accuracy, folds=2)[1][1] > 0.9
    @test optimize!(nnc, accuracy, folds=3)[1][1] > 0.9
    @test optimize!(nnc, accuracy, folds=5)[1][1] > 0.93
    @test optimize!(nnc, accuracy, folds=10)[1][1] > 0.93
    @show optimize!(nnc, accuracy, folds=5)
end

@testset "Clustering and centroid computation (with cosine)" begin
    using KernelMethods: fftraversal, sqrt_criterion, invindex, centroid!
    using SimilaritySearch: l2_distance, l1_distance, angle_distance, cosine_distance
    X, y = loadiris()
    dist = l2_distance
    refs = Vector{typeof(X[1])}()
    dmax = 0.0

    function callback(c, _dmax)
        push!(refs, X[c])
        dmax = _dmax
    end

    fftraversal(callback, dist, X, sqrt_criterion())
    R = fit(Sequential, refs)
    a = [centroid!(X[plist]) for plist in invindex(dist, X, R)]
    g = gaussian_kernel(dist, dmax/4)
    M = kmap(X, g, a)
    nnc = NearNeighborClassifier(cosine_distance, [normalize!(w) for w in M], y)

    @test optimize!(nnc, accuracy, folds=2)[1][1] > 0.8
    @test optimize!(nnc, accuracy, folds=3)[1][1] > 0.8
    @show optimize!(nnc, accuracy, folds=10)
end

@testset "encode with dnet" begin
    using KernelMethods
    using SimilaritySearch: l2_distance, normalize!, angle_distance
    using Statistics

    X, y = loadiris()
    dist = l2_distance
    # criterion = change_criterion(0.01)
    refs = Vector{typeof(X[1])}()
    dmax = 0.0

    function callback(c, dmaxlist)
        push!(refs, X[c])
        dmax += last(dmaxlist).dist
    end

    dnet(callback, dist, X, 14)
    _dmax = dmax / length(refs)
    g = tanh_kernel(dist, _dmax)
    M = kmap(X, g, refs)
    nnc = NearNeighborClassifier(l2_distance, M, y)
    @test optimize!(nnc, accuracy, folds=2)[1][1] > 0.9
    @test optimize!(nnc, accuracy, folds=3)[1][1] > 0.9
    @test optimize!(nnc, accuracy, folds=5)[1][1] > 0.9
    @test optimize!(nnc, accuracy, folds=10)[1][1] > 0.9
    @show optimize!(nnc, accuracy, folds=5)

    C = fftclustering(angle_distance, [normalize!(x) for x in X], 21, k=3)

    matches = 0
    for (i, res) in enumerate(C.NN)
        label = Dict{eltype(typeof(y)),Float64}()
        for (pos, p) in enumerate(res)
            l = y[p.objID]
            label[l] = get(label, l, 0) + 1 / pos
        end

        L = [(v, k) for (k, v) in label]
        sort!(L)
        if y[i] == L[end][end]
            matches += 1
        end
    end

    @info "===== accuracy by fftclustering: $(matches/length(y))"
    @test matches/length(y) > 0.9
end


text:
@testset "Interface" begin
    rng = Random.MersenneTwister(42)
    D = 5
    prior = MvNormal(0.5 * ones(D))
    likelihood = MvNormal(7.0 * ones(D))
    logintegrand = x -> logpdf(likelihood, x)
    m = BayesModel(prior, logintegrand)
    s = PriorSampling()
    bquad = BayesQuad(SEKernel() ∘ ScaleTransform(0.1))
    pZ, x = bquad(m, s; nsamples=100)
    trueZ = exp(-0.5 * (logdet(cov(prior) + cov(likelihood)) + D * log(2π)))
    @test mean(pZ) ≈ trueZ atol = 1e-3
end


text:
@testset verbose=true "Distance" begin
    @testset "Int" begin
        scalars = Int[0, 1, -1]
        quaternions = [Quaternion(a, b, c, d) for a in scalars for b in scalars for c in scalars for d in scalars]
        d(q₁, q₂) = distance2(q₁, q₂)
        for q₁ in quaternions
            for q₂ in quaternions
                d₁₁ = d(q₁, q₁)
                d₁₂ = d(q₁, q₂)
                d₂₂ = d(q₂, q₂)
                d₂₁ = d(q₂, q₁)
                # Real-valued:
                @test d₁₁ isa Int
                @test d₁₂ isa Int
                @test d₂₁ isa Int
                @test d₂₂ isa Int
                # Symmetry:
                @test d₁₂ == d₂₁
                # Invariance:
                for q₃ in [Quaternion(1), 1imx, 1imy, 1imz]
                    @test d₁₂ == d(q₃*q₁, q₃*q₂)
                    @test d₁₂ == d(q₁*q₃, q₂*q₃)
                end
                # Identity:
                @test d₁₁ == 0
                @test d₂₂ == 0
                # Positive-definiteness
                @test (q₁==q₂) ⊻ (d₁₂>0)
            end
        end
    end
    @testset "$T" for T in [BigFloat, Float64]  # FloatTypes
        scalars = [zero(T), one(T), -one(T)]
        quaternions = [a + b*imn for a in scalars for b in scalars for imn in [imx, imy, imz]]
        quaternions = [quaternions...; randn(Quaternion{T}, 23)...]
        # quaternions = [Quaternion(a, b, c, d) for a in scalars for b in scalars for c in scalars for d in scalars]
        # quaternions = [quaternions...; randn(Quaternion{T}, 19)...]
        @testset "rotor=$rotor" verbose=true for rotor in [true, false]
            @testset "sqrt=$return_sqrt" verbose=true for return_sqrt in [true, false]
                d(q₁, q₂) = return_sqrt ? distance(q₁, q₂) : distance2(q₁, q₂)
                for (i, q₁) in enumerate(quaternions)
                    if rotor
                        if iszero(q₁)
                            continue
                        end
                        q₁ = Rotor(q₁)
                    end
                    for q₂ in quaternions[i+1:end]
                        if rotor
                            if iszero(q₂)
                                continue
                            end
                            q₂ = Rotor(q₂)
                        end
                        absq₁, absq₂ = abs(q₁), abs(q₂)
                        rtol = 10eps(T)
                        atol = 10eps(max(absq₁, absq₂))
                        d₁₁ = d(q₁, q₁)
                        d₁₂ = d(q₁, q₂)
                        d₂₁ = d(q₂, q₁)
                        d₂₂ = d(q₂, q₂)
                        @test d₁₁ isa T
                        @test d₁₂ isa T
                        @test d₂₁ isa T
                        @test d₂₂ isa T
                        @test d₁₂ ≈ d₂₁ rtol=rtol atol=atol
                        for q₃ in quaternions
                            if iszero(q₃)
                                continue
                            end
                            if rotor
                                q₃ = Rotor(q₃)
                            else
                                q₃ = q₃ / abs(q₃)
                            end
                            @test d₁₂ ≈ d(q₃*q₁, q₃*q₂) rtol=rtol atol=atol
                            @test d₁₂ ≈ d(q₁*q₃, q₂*q₃) rtol=rtol atol=atol
                        end
                        @test d₁₁ ≈ zero(T) rtol=rtol atol=atol
                        @test d₂₂ ≈ zero(T) rtol=rtol atol=atol
                        if rotor
                            scalar_multiples = (≈(q₁, q₂, rtol=rtol, atol=atol) || ≈(q₁, -q₂, rtol=rtol, atol=atol))
                            @test scalar_multiples ⊻ (d₁₂ > atol)
                        else
                            @test (q₁==q₂) ⊻ (d₁₂>0)
                        end
                    end
                end
            end
        end
    end
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
wedge.jl defines methods for the wedge(x, y) function
"""

# --- Exports

export wedge, ∧

# --- Operator aliases

∧(M::AbstractMultivector, N::AbstractMultivector) = wedge(M, N)

∧(M::AbstractMultivector, x::Real) = wedge(M, x)
∧(x::Real, M::AbstractMultivector) = wedge(x, M)

∧(M::AbstractMultivector, v::Vector{<:Real}) = wedge(M, v)
∧(v::Vector{<:Real}, M::AbstractMultivector) = wedge(v, M)

∧(v::Vector{<:Real}, w::Vector{<:Real}) = wedge(v, w)

# --- Method definitions

"""
    wedge(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector
    ∧(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

Compute the outer product of `M` and `N`.
"""
function wedge end

# ------ Specializations involving an AbstractMultivector instance

# M::AbstractMultivector, B::AbstractBlade
# B::AbstractBlade, M::AbstractMultivector
wedge(M::AbstractMultivector, B::AbstractBlade) =
    Multivector(map(C -> wedge(C, B), blades(M)))

wedge(B::AbstractBlade, M::AbstractMultivector) =
    Multivector(map(C -> wedge(B, C), blades(M)))

# M::AbstractMultivector, v::Vector
# v::Vector, M::AbstractMultivector
wedge(M::AbstractMultivector, v::Vector{<:Real}) =
    Multivector(map(C -> wedge(C, v), blades(M)))

wedge(v::Vector{<:Real}, M::AbstractMultivector) =
    Multivector(map(C -> wedge(v, C), blades(M)))

# M::AbstractMultivector, B::AbstractScalar
# B::AbstractScalar, M::AbstractMultivector
wedge(B::AbstractMultivector, C::AbstractScalar) = B * C
wedge(B::AbstractScalar, C::AbstractMultivector) = B * C

# M::AbstractMultivector, B::One
# B::One, M::AbstractMultivector
wedge(M::AbstractMultivector, B::One) = M
wedge(B::One, M::AbstractMultivector) = M

# M::AbstractMultivector, B::Zero
# B::Zero, M::AbstractMultivector
wedge(B::Zero, M::AbstractMultivector) = B
wedge(M::AbstractMultivector, B::Zero) = B

# M::AbstractMultivector, x::Real
# x::Real, M::AbstractMultivector
wedge(M::AbstractMultivector, x::Real) = x * M
wedge(x::Real, M::AbstractMultivector) = x * M

# ------ Specializations involving an AbstractBlade instance

# B::AbstractBlade, C::AbstractScalar
# B::AbstractScalar, C::AbstractBlade
wedge(B::AbstractBlade, C::AbstractScalar) = B * C
wedge(B::AbstractScalar, C::AbstractBlade) = B * C

# B::AbstractBlade, C::One
# B::One, C::AbstractBlade
wedge(B::AbstractBlade, C::One) = B
wedge(B::One, C::AbstractBlade) = C

# B::AbstractBlade, C::Zero
# B::Zero, C::AbstractBlade
wedge(B::AbstractBlade, C::Zero) = zero(B)
wedge(B::Zero, C::AbstractBlade) = zero(C)

# B::AbstractBlade, x::Real
# x::Real, B::AbstractBlade
wedge(B::AbstractBlade, x::Real) = B * x
wedge(x::Real, B::AbstractBlade) = x * B

# ------ Specializations involving a Blade instance

# B::Blade, C::Blade
function wedge(B::Blade, C::Blade)
    # --- Check arguments

    assert_dim_equal(B, C)

    if grade(B) + grade(C) > dim(B)
        return zero(B)
    end

    # --- Construct new blade

    volume(B) * volume(C) * Blade(hcat(basis(B), basis(C)))
end

# B::Blade, C::Pseudoscalar
# B::Pseudoscalar, C::Blade
function wedge(B::Blade, C::Pseudoscalar)
    assert_dim_equal(B, C)
    zero(B)
end

wedge(B::Pseudoscalar, C::Blade) = wedge(C, B)

# B::Blade, v::Vector
# v::Vector, B::Blade
function wedge(B::Blade, v::Vector{<:Real})
    assert_dim_equal(B, v)

    # Note: volume(B) is incorporated into the norm of `v`
    Blade(hcat(basis(B), volume(B) * v))
end

function wedge(v::Vector{<:Real}, B::Blade)
    assert_dim_equal(v, B)

    # Note: volume(B) is incorporated into the norm of `v`
    Blade(hcat(volume(B) * v, basis(B)))
end

# v::Vector, w::Vector
wedge(v::Vector{<:Real}, w::Vector{<:Real}) = Blade(hcat(v, w))

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
function wedge(B::Pseudoscalar, C::Pseudoscalar)
    assert_dim_equal(B, C)
    zero(B)
end

# B::Pseudoscalar, v::Vector
# v::Vector, B::Pseudoscalar
function wedge(B::Pseudoscalar, v::Vector{<:Real})
    assert_dim_equal(B, v)
    zero(B)
end

wedge(v::Vector{<:Real}, B::Pseudoscalar) = B ∧ v

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::One
# B::One, C::AbstractScalar
wedge(B::AbstractScalar, C::One) = B
wedge(B::One, C::AbstractScalar) = C

# B::AbstractScalar, C::Zero
# B::Zero, C::AbstractScalar
wedge(B::AbstractScalar, C::Zero) = C
wedge(B::Zero, C::AbstractScalar) = B

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
wedge(B::AbstractScalar, x::Real) = B * x
wedge(x::Real, B::AbstractScalar) = x * B

# B::AbstractScalar, v::Vector
# v::Vector, B::AbstractScalar
wedge(B::AbstractScalar, v::Vector{<:Real}) = value(B) * Blade(v)
wedge(v::Vector{<:Real}, B::AbstractScalar) = wedge(B, v)

# ------ Specializations involving a Scalar instance

# B::Scalar, C::Scalar
wedge(B::Scalar, C::Scalar) = B * C

# ------ Specializations involving a One instance

# B::One, C::One
wedge(B::One, C::One) = B

# B::One, C::Zero
# B::Zero, C::One
wedge(B::One, C::Zero) = C
wedge(B::Zero, C::One) = B

# ------ Specializations involving a Zero instance

# B::Zero, C::Zero
wedge(B::Zero, C::Zero) = B


text:
import Thermodynamics
import StatsBase
import PrettyTables
import OrderedCollections
const TD = Thermodynamics
const TP = TD.Parameters
using JET
using Test

import UnPack
import BenchmarkTools

import CLIMAParameters
const CP = CLIMAParameters
const FT = Float64
toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
aliases = string.(fieldnames(TP.ThermodynamicsParameters))
param_pairs = CP.get_parameter_values!(toml_dict, aliases, "Thermodynamics")
const param_set = TP.ThermodynamicsParameters{FT}(; param_pairs...)

#####
##### Finding indexes in profiles satisfying certain conditions
#####

function find_freezing_index(profiles)
    i = findfirst(T -> T === TP.T_freeze(param_set), profiles.T)
    isnothing(i) && error("Freezing index not found")
    return i
end

function find_dry(profiles)
    i = findfirst(q_tot -> q_tot === 0.0, profiles.q_tot)
    isnothing(i) && error("Dry index not found")
    return i
end

function find_moist_index(profiles)
    i = findfirst(q_tot -> q_tot > 0.01, profiles.q_tot)
    isnothing(i) && error("Moist index not found")
    return i
end

function find_sat_adjust_index(profiles, f::TDC) where {TDC}
    kwargs = get_kwargs(profiles, f)
    ts = f.(param_set, values(kwargs)...) # TDC is the thermo constructor
    T_sa = TD.air_temperature.(param_set, ts)
    Z = zip(T_sa, profiles.T)
    i = findfirst(x -> abs(x[1] - x[2]) > 0.0001, collect(Z))
    isnothing(i) && error("Saturation adjustment index not found")
    # @info "T_error for i_sat_adjust = $(abs(T_sa[i] - T[i]))"
    return i
end

#=
    conditions_index

Find an index in the given `profiles` that satisfies
some condition. For example, find the index that we're
sure that saturation adjustment is performed for a given
(profiles.p[i], profiles.θ_liq_ice[i], profiles.q_tot[i]).
=#
function conditions_index(profiles, sym, constructor)
    i = if sym == :dry
        find_dry(profiles)
    elseif sym == :freezing
        find_freezing_index(profiles)
    elseif sym == :sat_adjust
        find_sat_adjust_index(profiles, constructor)
    elseif sym == :moist
        find_moist_index(profiles)
    else
        error("Bad sym given")
    end
    return i
end

function sample_args(profiles, sym, constructor)
    i = conditions_index(profiles, sym, constructor)
    kwargs = get_kwargs(profiles, constructor)
    return getindex.(values(kwargs), i)
end

#####
##### BenchmarkTools's trial utils
#####

get_summary(trial) = (;
    # Using some BenchmarkTools internals :/
    mem = BenchmarkTools.prettymemory(trial.memory),
    nalloc = trial.allocs,
    t_min = BenchmarkTools.prettytime(minimum(trial.times)),
    t_max = BenchmarkTools.prettytime(maximum(trial.times)),
    t_mean = BenchmarkTools.prettytime(StatsBase.mean(trial.times)),
    t_med = BenchmarkTools.prettytime(StatsBase.median(trial.times)),
    n_samples = length(trial),
)

function tabulate_summary(summary)
    summary_keys = collect(keys(summary))
    mem = map(k -> summary[k].mem, summary_keys)
    nalloc = map(k -> summary[k].nalloc, summary_keys)
    t_mean = map(k -> summary[k].t_mean, summary_keys)
    t_min = map(k -> summary[k].t_min, summary_keys)
    t_max = map(k -> summary[k].t_max, summary_keys)
    t_med = map(k -> summary[k].t_med, summary_keys)
    n_samples = map(k -> summary[k].n_samples, summary_keys)

    table_data = hcat(
        string.(collect(keys(summary))),
        mem,
        nalloc,
        t_min,
        t_max,
        t_mean,
        t_med,
        n_samples,
    )

    header = (
        [
            "Constructor",
            "Memory",
            "allocs",
            "Time",
            "Time",
            "Time",
            "Time",
            "N-samples",
        ],
        [
            "(+conditions)",
            "estimate",
            "estimate",
            "min",
            "max",
            "mean",
            "median",
            "",
        ],
    )

    println()
    PrettyTables.pretty_table(
        table_data;
        header,
        crop = :none,
        alignment = vcat(:l, repeat([:r], length(header[1]) - 1)),
    )
end

#####
##### Constructor-specific configurations
#####

# unpack variables from profiles into NamedTuple:
up(profiles, syms) = (; zip(syms, getproperty.(Ref(profiles), syms))...)
get_kwargs(p, ::typeof(TD.PhaseEquil_ρeq)) = up(p, :(ρ, e_int, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_ρTq)) = up(p, :(p, T, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_pTq)) = up(p, :(p, T, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_peq)) = up(p, :(p, e_int, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_phq)) = up(p, :(p, h, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_ρθq)) = up(p, :(ρ, θ_liq_ice, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_pθq)) = up(p, :(p, θ_liq_ice, q_tot).args)
get_kwargs(p, ::typeof(TD.PhaseEquil_ρpq)) = up(p, :(ρ, p, q_tot).args)

# Conditions to perform microbenchmarks and JET tests:
# note: No freezing points exist in
#       TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)!
#       so we're not testing performance of these branches.

conditions(::typeof(TD.PhaseEquil_ρeq)) = (:dry, :sat_adjust)
conditions(::typeof(TD.PhaseEquil_ρTq)) = (:dry, :moist)
conditions(::typeof(TD.PhaseEquil_pTq)) = (:dry, :moist) # no sat adjust exists!
conditions(::typeof(TD.PhaseEquil_peq)) = (:dry, :sat_adjust)
conditions(::typeof(TD.PhaseEquil_phq)) = (:dry, :sat_adjust)
conditions(::typeof(TD.PhaseEquil_ρθq)) = (:dry, :sat_adjust)
conditions(::typeof(TD.PhaseEquil_pθq)) = (:dry, :sat_adjust)
conditions(::typeof(TD.PhaseEquil_ρpq)) = (:dry, :sat_adjust)


text:
#=
These are functions which calculate the marginal likelihood of a segment of a time series
(upto a constant of proprtionality)

P(s,t) = \int \Pr(y_{s+1:t} | \theta) \pi(\theta) d\theta
where theta is the parameter that changes and \pi(\theta) is its prior
=#

# the first function assumes a constant variance and chnging mean where the mean has a normal prior
# centered on 0 and has standard deviation s
function NormalMeanSegment(data::Array{Float64})
    cd = [0,cumsum( data )]
    cd_2 = [0,cumsum( abs2(data) )]
    cost(s::Int64, t::Int64) = cd_2[t+1] - cd_2[s+1] - abs2(cd[t+1] - cd[s+1])/(t-s)
    return cost
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.20"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
dists = [
    Uniform(10,20),
    Exponential(3.5),
    Gamma(0.5,7),
    Beta(10,0.5),
    Weibull(10,0.5),
    Normal(20,3.5),
    Rayleigh(2.4),
    Cauchy(20,3.5)]

# ╔═╡ 863fa0fe-2164-11eb-221c-fbdef5124dc6
begin
	df = DataFrame()
	df.distributions = dists
	df.parameters = params.(dists)
	df.support = ((d)->(minimum(d),maximum(d))).(dists)
	df
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.20"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═863fa0fe-2164-11eb-221c-fbdef5124dc6
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
module Optimization

using Evolutionary

import BlackBoxOptim

import ..QuerySynthesisStrategies:
    check_epsilon,
    data_boundaries

include("optimization_base.jl")
include("BlackBoxOptimization.jl")
include("EvolutionaryOptimization.jl")
include("ParticleSwarmOptimization.jl")

export
    QuerySynthesisOptimizer,
    BlackBoxOptimization,
    EvolutionaryOptimization,
    ParticleSwarmOptimization,

    query_synthesis_optimize
end


text:
require("FileFind")

# The simplest scenario possible
FileFind.find()

# Specifying a directory explicitly
FileFind.find(".")

# Specifying a function explicitly
FileFind.find(fn -> println(isfile(fn)))

# Specifying both a directory and a function explicitly
function file_name_and_size(filename)
  println(strcat(filename, " ", filesize(filename)))
end
FileFind.find(".", file_name_and_size)

# Using a closure to maintain global state across function calls
total_bytes = 0

function total_size(filename)
  global total_bytes
  total_bytes += filesize(filename)
end

FileFind.find(".", total_size)

total_bytes


text:
push!(LOAD_PATH, joinpath(homedir(),"GitHub/KenyaCoV/src"))
using Plots,Parameters,Distributions,JLD2,DataFrames,StatsPlots,FileIO,MAT
using Statistics: median, quantile
gr()
treatment_rates = [(0.,1),(0.,0.5),(1/7.,1.),(1/7,0.5),(1/3.5,1.),(1/3.5,0.5)]
reducted_treatment_rates = [(0.,1),(1/3.5,0.5)]
reduced_rel_transmission_perc = [10,25]
rel_transmission_perc = [0,10,25,50,100]

"""
Plots ---
1) Early growth for both scenarios
2) Peaks
3) Size and age distribution
"""

@load joinpath("KenyaCoVOutputs/results_1.jld2") results_1
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_2.jld2") results_2
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3.jld2") results_3
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_4.jld2") results_4
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_5.jld2") results_5
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_2S.jld2") results_2S
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3S.jld2") results_3S
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_2SI.jld2") results_2SI
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3SI.jld2") results_3SI
    # @load joinpath(homedir(),"Github/KenyaCoVOutputs/results_2SI_short.jld2") results_2SI_short
    # @load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3SI_short.jld2") results_3SI_short
    # @load joinpath(homedir(),"Github/KenyaCoVOutputs/results_4SI_short.jld2") results_4SI_short
    # @load joinpath(homedir(),"Github/KenyaCoVOutputs/results_5SI_short.jld2") results_5SI_short

@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_2SI_SL_three_months.jld2") results_2SI_SL_three_months
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3SI_SL_three_months.jld2") results_3SI_SL_three_months
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_4SI_SL_three_months.jld2") results_4SI_SL_three_months
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_5SI_SL_three_months.jld2") results_5SI_SL_three_months

@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_2SI_SL_DC_three_months.jld2") results_2SI_SL_DC_three_months
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_3SI_SL_DC_three_months.jld2") results_3SI_SL_DC_three_months
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_4SI_SL_DC_three_months.jld2") results_4SI_SL_DC_three_months
@load joinpath(homedir(),"Github/KenyaCoVOutputs/results_5SI_SL_DC_three_months.jld2") results_5SI_SL_DC_three_months

# gr()
include("plotting_functions.jl");
"""
Doubling time
"""

"""
Country wide incidence
"""

scenario_group = [results_1,results_2,results_3,results_4,results_5]
plt_no_control = plot_total_incidence_group(scenario_group,treatment_rates,1,rel_transmission_perc)
title!(plt_no_control,"Kenya-wide incidence: no intervention")
plot!(plt_no_control,size = (700,500))
xlabel!("Days after established SARS-CoV-2 transmission")
ylabel!("Daily new sympomatic infections")
savefig(plt_no_control,"plotting/baseline_scenarios.pdf")
savefig(plt_no_control,"plotting/baseline_scenarios.png")

savefig(plt_no_control,"plotting/baseline_scenarios_eps_test.eps")

plt_controls = plot_total_incidence_group(scenario_group,treatment_rates,6,rel_transmission_perc)
title!(plt_controls,"Kenya-wide incidence: case isolation")
plot!(plt_controls,size = (700,400))
savefig(plt_controls,"plotting/target_scenarios.pdf")

scenario_group_SD = [results_2S,results_3S]
plt_no_control_SD = plot_total_incidence_group(scenario_group_SD,reducted_treatment_rates,1,reduced_rel_transmission_perc)
title!(plt_no_control_SD,"Social distancing only")
plot!(plt_no_control_SD,size = (700,400))
plot!(plt_no_control_SD,legend = :topleft)
# xlims!(0.,30.)
savefig(plt_no_control_SD,"plotting/baseline_scenarios_with_social_distancing.pdf")

scenario_group_SD = [results_2S,results_3S]
plt_control_SD = plot_total_incidence_group(scenario_group_SD,reducted_treatment_rates,2,reduced_rel_transmission_perc)
title!(plt_control_SD,"Social distancing and direct isolation")
plot!(plt_control_SD,size = (700,400))
plot!(plt_control_SD,legend = :topleft)
# xlims!(0.,30.)
savefig(plt_control_SD,"plotting/direct_control_scenarios_with_social_distancing.pdf")

scenario_group_SDI = [results_2SI,results_3SI]
plt_no_control_SD = plot_total_incidence_group(scenario_group_SDI,reducted_treatment_rates,1,reduced_rel_transmission_perc)
title!(plt_no_control_SD,"Social distancing only")
plot!(plt_no_control_SD,size = (700,400))
plot!(plt_no_control_SD,legend = :topright)
# xlims!(0.,30.)
savefig(plt_no_control_SD,"plotting/baseline_scenarios_with_social_distancing.pdf")

plt_control_SD = plot_total_incidence_group(scenario_group_SDI,reducted_treatment_rates,2,reduced_rel_transmission_perc)
title!(plt_control_SD,"Social distancing and case isolation")
plot!(plt_control_SD,size = (700,400))
plot!(plt_control_SD,legend = :topright)
savefig(plt_control_SD,"plotting/baseline_scenarios_with_social_distancing_and_case_isolation.pdf")

plt_1 = plot([1. for i =1:(2*365)],
        fontfamily="Helvetica",
        lw = 3,
        lab= "rel. infect. of subclinical: 0%",
        legend = :topright,
        fillalpha = 0.15,
        yticks = ([1,10,100,1000,10000],["0" "10" "100" "1000" "10000"]),
        yscale = :log10,
        xlabel = "Days after detecting established CoV transmission",
        ylabel = "Daily incidence")
plt_control_SD = plot_total_incidence_group(plt_1,scenario_group_SDI,reducted_treatment_rates,2,reduced_rel_transmission_perc)
title!(plt_control_SD,"Social distancing and case isolation")
plot!(plt_control_SD,size = (700,400))
plot!(plt_control_SD,legend = :topright)
plot!(plt_control_SD,results_2[6][1][21,:,1].+1,ls=:dash,color = :red,lab = "")
plot!(plt_control_SD,results_3[6][1][21,:,1].+1,ls=:dash,color = :green,lab = "")
savefig(plt_control_SD,"plotting/baseline_scenarios_with_social_distancing_and_case_isolation.pdf")

scenario_group_lockdown = [results_2SI_SL_three_months,
                            results_3SI_SL_three_months,
                            results_4SI_SL_three_months,
                            results_5SI_SL_three_months]

plt_2 = plot([1. for i =1:(2*365)],
    fontfamily="Helvetica",
    lw = 3,
    lab= "rel. infect. of subclinical: 0%",
    legend = :topright,
    fillalpha = 0.15,
    yticks = ([1,11,101,1001,10001],["0" "10" "100" "1000" "10000"]),
    yscale = :log10,
    xlabel = "Days after detecting established CoV transmission",
    ylabel = "Daily incidence")
plt_lockdown = plot_total_incidence_group(plt_2,scenario_group_lockdown,[(0.5,1/3.5)],1,[10,25,50,100])
plot!([90,90],[1,2.5e4],color = :black,ls = :dot,lw = 3,lab="")
xlims!(0,400)
title!("Daily incidence before and after 90 days restrictions (SD+MR+CI)")
plot!(plt_lockdown,results_2[6][1][21,:,1].+1,ls=:dash,color = :red,lab = "")
plot!(plt_lockdown,results_3[6][1][21,:,1].+1,ls=:dash,color = :green,lab = "")
plot!(plt_lockdown,results_4[6][1][21,:,1].+1,ls=:dash,color = :purple,lab = "")
plot!(plt_lockdown,results_5[6][1][21,:,1].+1,ls=:dash,color = :brown,lab = "")

plot!(size = (700,400))
savefig(plt_lockdown,"plotting/baseline_scenarios_with_lockdown.pdf")

scenario_group_lockdown_reduction = [results_2SI_SL_DC_three_months,
                            results_3SI_SL_DC_three_months,
                            results_4SI_SL_DC_three_months,
                            results_5SI_SL_DC_three_months]

plt_3 = plot([1. for i =1:(2*365)],
    fontfamily="Helvetica",
    lw = 3,
    lab= "rel. infect. of subclinical: 0%",
    legend = :topright,
    fillalpha = 0.15,
    yticks = ([1,11,101,1001,10001],["0" "10" "100" "1000" "10000"]),
    yscale = :log10,
    xlabel = "Days after detecting established CoV transmission",
    ylabel = "Daily incidence")
plt_lockdown_red = plot_total_incidence_group(plt_3,scenario_group_lockdown_reduction,[(0.5,1/3.5)],1,[10,25,50,100])
plot!([90,90],[1,2.5e4],color = :black,ls = :dot,lw = 3,lab="")
xlims!(0,400)
title!("Daily incidence before and after 90 days restrictions (SD+MR+CI)")
plot!(plt_lockdown_red,results_2[6][1][21,:,1].+1,ls=:dash,color = :red,lab = "")
plot!(plt_lockdown_red,results_3[6][1][21,:,1].+1,ls=:dash,color = :green,lab = "")
plot!(plt_lockdown_red,results_4[6][1][21,:,1].+1,ls=:dash,color = :purple,lab = "")
plot!(plt_lockdown_red,results_5[6][1][21,:,1].+1,ls=:dash,color = :brown,lab = "")

plot!(size = (700,400))
xlabel!("Days after established SARS-CoV-2 transmission")
ylabel!("Daily new sympomatic infections")
savefig(plt_lockdown_red,"plotting/baseline_scenarios_with_lockdown_and_reduction.pdf")

"""
Spatial plots
"""
ordering = sortperm(results_3[1][1][1:20,40,1],rev = true)
results_3[1][1][1:20,40,1]
plot(results_3[1][1][5,:,1].+1,yscale = :log10)
plt_no_control_1 = plot_incidence_spatial(results_1,treatment_rates,1)
title!(plt_no_control_1,"Uncontrolled epidemic: 0% rel. trans. undetecteds ")
plot!(plt_no_control_1,size = (700,500))
savefig(plt_no_control_1,"plotting/spatial_baseline_scenarios_tau_0.pdf")

plt_control_1 = plot_incidence_spatial(results_1,treatment_rates,6)
title!(plt_control_1,"Targetted intervention: 0% rel. trans. undetecteds ")
plot!(plt_control_1,size = (700,500))
savefig(plt_control_1,"plotting/spatial_targetting_scenarios_tau_0.pdf")

plt_no_control_3 = plot_incidence_spatial(results_3,treatment_rates,1)
title!(plt_no_control_3,"Spatial distribution of incidence: (25% rel infect.)")
plot!(size = (700,500))
xlabel!("Days after established SARS-CoV-2 transmission")
ylabel!("Daily new sympomatic infections")
xlims!(0.,60.)

savefig(plt_no_control_3,"plotting/spatial_baseline_scenarios_epsilon_A_25_ribbon_plot.pdf")
plt_no_control_3 = plot_incidence_spatial(results_3,treatment_rates,1,rr2county_mat)
title!(plt_no_control_3,"Spatial distribution of incidence: (25% rel infect.)")
plot!(size = (1000,1050),
        xtickfont = font(14),
        ytickfont = font(14),
        guidefont = font(14),
        titlefont = font(22)
        )
xlabel!("")
ylabel!("")
xlabel!("Days after established SARS-CoV-2 transmission")
ylabel!("Daily new sympomatic infections")
xlims!(0.,60.)

savefig(plt_no_control_3,"plotting/spatial_baseline_scenarios_epsilon_A_25_ribbon_plot_county.pdf")

plt_control_2 = plot_incidence_spatial(results_2,treatment_rates,6)
title!(plt_control_2,"Targetted intervention: 10% rel. trans. undetecteds ")
plot!(plt_control_2,size = (700,500))
savefig(plt_control_2,"plotting/spatial_targetting_scenarios_tau_01.pdf")

plt_control_SD_2 = plot_incidence_spatial(results_2S,treatment_rates,1)
xlims!(plt_control_SD_2,(0.,150))
title!(plt_control_SD_2,"Social distancing: 10% rel. trans. undetecteds ")
plot!(plt_control_SD_2,size = (700,500))
savefig(plt_control_SD_2,"plotting/spatial_targetting_scenarios_tau_01_social_distancing.pdf")

"""
PEAK TIMING
"""
collect_no_control_peak_timing = hcat(results_1[1][3][:,21],
                                      results_2[1][3][:,21],
                                      results_3[1][3][:,21],
                                      results_4[1][3][:,21],
                                      results_5[1][3][:,21])
median_estimates_peak_uncontrolled = [median(collect_no_control_peak_timing[:,i]) for i = 1:5]
lb_ub = [quantile(collect_no_control_peak_timing[:,i],[0.025,0.975]) for i = 1:5]

peak_plt_no_control = boxplot(collect_no_control_peak_timing,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:5,["0%","10%","25%","50%","100%"]))
title!(peak_plt_no_control,"Uncontrolled epidemic peak timing")
ylabel!(peak_plt_no_control,"Time to peak in days")
xlabel!(peak_plt_no_control,"Rel. infectiousness of subclinical cases")
plot!(size = (700,400))
savefig(peak_plt_no_control,"plotting/time_to_peak_uncontrolled.pdf")

collect_control_peak_timing = hcat(results_1[6][3][:,21],
                                      results_2[6][3][:,21],
                                      results_3[6][3][:,21],
                                      results_4[6][3][:,21],
                                      results_5[6][3][:,21])

median_estimates_peak_CI = [median(collect_control_peak_timing[:,i]) for i = 1:5]
lb_ub = [quantile(collect_control_peak_timing[:,i],[0.025,0.975]) for i = 1:5]

peak_plt_control = boxplot(collect_control_peak_timing,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:5,["0%","10%","25%","50%","100%"]))
title!(peak_plt_control,"Epidemic peak timing with case isolation")
ylabel!(peak_plt_control,"Time to peak in days")
xlabel!(peak_plt_control,"Rel. infectiousness of subclinical cases")
ylims!((0.,300.))
plot!(size = (700,400))
savefig(peak_plt_control,"plotting/time_to_peak_controlled.pdf")

collect_SDI_peak_timing = hcat(results_2SI[1][3][:,21],
                                      results_3SI[1][3][:,21],
                                      results_2SI[2][3][:,21],
                                      results_3SI[2][3][:,21])

peak_plt_SD = boxplot(collect_SDI_peak_timing,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:4,["10% - SD","25% - SD","10% - SD+CI","25% - SD+CI"]))
title!(peak_plt_SD,"Epidemic peak timing with social distancing and case isolation")
ylabel!(peak_plt_SD,"Time to peak in days")
xlabel!(peak_plt_SD,"Rel. infectiousness of undetected cases")
plot!(size = (700,400))
savefig(peak_plt_SD,"plotting/time_to_peak_social_distancing_and_case_isolation.pdf")

collect_lockdown_peak_timing = hcat(results_2SI_SL_three_months[1][3][:,21],
                                      results_3SI_SL_three_months[1][3][:,21],
                                      results_4SI_SL_three_months[1][3][:,21],
                                      results_5SI_SL_three_months[1][3][:,21])

peak_plt_SD = boxplot(collect_lockdown_peak_timing,lab="",
                          # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                          xticks = (1:4,["10%","25%","50%","100%"]))
title!(peak_plt_SD,"Epidemic peak timing with social distancing and case isolation")
ylabel!(peak_plt_SD,"Time to peak in days")
xlabel!(peak_plt_SD,"Rel. infectiousness of subclinical cases")
plot!(size = (700,400))
savefig(peak_plt_SD,"plotting/time_to_peak_lockdown.pdf")

"""
Final size plots
"""

sum(results_1[1][4],dims= [1,2])[:]

collect_uncontrol_final_size = hcat(sum(results_1[1][4],dims= [1,2])[:],
                                    sum(results_2[1][4],dims= [1,2])[:],
                                    sum(results_3[1][4],dims= [1,2])[:],
                                    sum(results_4[1][4],dims= [1,2])[:],
                                    sum(results_5[1][4],dims= [1,2])[:])

median_estimates_uncontrolled = [median(collect_uncontrol_final_size[:,i]) for i = 1:5]
lb_ub = [quantile(collect_uncontrol_final_size[:,i],[0.025,0.975]) for i = 1:5]

collect_uncontrol_final_size_asymptoms = hcat(sum(results_1[1][5],dims= [1,2])[:],
                                    sum(results_2[1][5],dims= [1,2])[:],
                                    sum(results_3[1][5],dims= [1,2])[:],
                                    sum(results_4[1][5],dims= [1,2])[:],
                                    sum(results_5[1][5],dims= [1,2])[:])
median_estimates_uncontrolled_total = [median(collect_uncontrol_final_size_asymptoms[:,i].+ collect_uncontrol_final_size[:,i]) for i = 1:5]
lb_ub = [quantile(collect_uncontrol_final_size_asymptoms[:,i].+ collect_uncontrol_final_size[:,i],[0.025,0.975]) for i = 1:5]

size_plt_no_control = boxplot(collect_uncontrol_final_size./1e6,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:5,["0%","10%","25%","50%","100%"]))
title!(size_plt_no_control,"Uncontrolled epidemic number of cases")
ylabel!(size_plt_no_control,"Total cases (millions)")
xlabel!(size_plt_no_control,"Rel. infectiousness of undetected cases")
plot!(size = (700,400))
savefig(size_plt_no_control,"plotting/epidemic_size_no_control.pdf")

collect_control_final_size = hcat(sum(results_1[6][4],dims= [1,2])[:],
                                    sum(results_2[6][4],dims= [1,2])[:],
                                    sum(results_3[6][4],dims= [1,2])[:],
                                    sum(results_4[6][4],dims= [1,2])[:],
                                    sum(results_5[6][4],dims= [1,2])[:])

median_estimates_case_isolation = [median(collect_control_final_size[:,i]) for i = 1:5]
lb_ub = [quantile(collect_control_final_size[:,i],[0.025,0.975]) for i = 1:5]

collect_control_final_size_asymptoms = hcat(sum(results_1[6][5],dims= [1,2])[:],
                                    sum(results_2[6][5],dims= [1,2])[:],
                                    sum(results_3[6][5],dims= [1,2])[:],
                                    sum(results_4[6][5],dims= [1,2])[:],
                                    sum(results_5[6][5],dims= [1,2])[:])
median_estimates_controlled_total = [median(collect_control_final_size_asymptoms[:,i].+ collect_control_final_size[:,i]) for i = 1:5]
lb_ub = [quantile(collect_control_final_size_asymptoms[:,i].+ collect_control_final_size[:,i],[0.025,0.975]) for i = 1:5]

size_plt_control = boxplot(collect_control_final_size./1e6,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:5,["0%","10%","25%","50%","100%"]))
title!(size_plt_control,"Epidemic number of cases with case isolation")
ylabel!(size_plt_control,"Total cases (millions)")
xlabel!(size_plt_control,"Rel. infectiousness of undetected cases")
ylims!((0.,1.5))
plot!(size = (700,400))
savefig(size_plt_control,"plotting/epidemic_size_with_control.pdf")

collect_SD_size = hcat(sum(results_2S[1][4],dims= [1,2])[:],
                            sum(results_3S[1][4],dims= [1,2])[:],
                            sum(results_2S[2][4],dims= [1,2])[:],
                            sum(results_3S[2][4],dims= [1,2])[:])

median_estimates_case_isolation_and_perfect_social_distancing = [median(collect_SD_size[:,i]) for i = 1:4]

size_plt_SD = boxplot(collect_SD_size,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:4,["10% - no targeting","25% - no targeting","10% - targetting","25% - targetting"]))

collect_SDI_size = hcat(sum(results_2SI[1][4],dims= [1,2])[:],
                            sum(results_3SI[1][4],dims= [1,2])[:],
                            sum(results_2SI[2][4],dims= [1,2])[:],
                            sum(results_3SI[2][4],dims= [1,2])[:])

median_estimates_case_isolation_and_social_distancing = [median(collect_SDI_size[:,i]) for i = 1:4]

ratio_10_rel_infect_with_SD_CI = median_estimates_case_isolation_and_social_distancing[3]/median_estimates_case_isolation[2]
ratio_10_rel_infect_with_SD_CI = median_estimates_case_isolation_and_social_distancing[4]/median_estimates_case_isolation[3]

size_plt_SD = boxplot(collect_SDI_size./1e6,lab="",
                            # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                            xticks = (1:4,["10% - SD","25% - SD","10% - SD+CI","25% - SD+CI"]))
title!(size_plt_SD,"Epidemic number of cases with social distancing and case isolation")
ylabel!(size_plt_SD,"Total cases (millions)")
xlabel!(size_plt_SD,"Rel. infectiousness of undetected cases")
ylims!((0.,1.5))
plot!(size = (700,400))
savefig(size_plt_SD,"plotting/epidemic_size_with_SD_and_case_isolation.pdf")

collect_SDI_comparison_size = hcat(sum(results_2[6][4],dims= [1,2])[:],
                                sum(results_2SI[2][4],dims= [1,2])[:],
                                sum(results_3[6][4],dims= [1,2])[:],
                            sum(results_3SI[2][4],dims= [1,2])[:])
size_plt_SD_comparison = boxplot(collect_SDI_comparison_size./1e6,lab="",
                        # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                                xticks = (1:4,["10% - CI","10% - CI+SD","25% - CI","25% - CI+SD"]))
title!(size_plt_SD_comparison,"Epidemic number of cases with long-lasting social distancing and case isolation")
ylabel!(size_plt_SD_comparison,"Total cases (millions)")
xlabel!(size_plt_SD_comparison,"Rel. infectiousness of undetected cases and controls")
ylims!((0.,1.5))
plot!(size = (700,400))
savefig(size_plt_SD_comparison,"plotting/epidemic_size_comparison_SD_and_case_isolation.pdf")

collect_lockdown_comparison_size = hcat(sum(results_2[6][4],dims= [1,2])[:],
                            sum(results_2SI_SL_three_months[1][4],dims= [1,2])[:],
                            sum(results_3[6][4],dims= [1,2])[:],
                            sum(results_3SI_SL_three_months[1][4],dims= [1,2])[:],
                            sum(results_4[6][4],dims= [1,2])[:],
                            sum(results_4SI_SL_three_months[1][4],dims= [1,2])[:],
                            sum(results_5[6][4],dims= [1,2])[:],
                            sum(results_5SI_SL_three_months[1][4],dims= [1,2])[:])
size_plt_SD_compare = boxplot(collect_lockdown_comparison_size./1e6,lab="",
                        # lab = ["rel. infect. undetecteds = 0%" "rel. infect. undetecteds = 10%" "rel. infect. undetecteds = 25%" "rel. infect. undetecteds = 50%" "rel. infect. undetecteds = 100%"],
                        xticks = (1:8,["10% - CI","10% - CI+tSD+tMR","25% - CI","25% - CI+tSD+tMR",
                                "50%","50% - CI+tSD+tMR","100%","100% - CI+tSD+tMR"]))

title!(size_plt_SD_compare,"Epidemic number of cases with 90 days controls")
ylabel!(size_plt_SD_compare,"Total cases (millions)")
xlabel!(size_plt_SD_compare,"Rel. infectiousness of subclinical cases and controls")
ylims!((0.,1.5))
plot!(size = (900,400))
savefig(size_plt_SD_compare,"plotting/epidemic_size_comparison_lockdown.pdf")
"""
Age distribution of cases
"""

scenario_group = [results_1,results_2,results_3,results_4,results_5]
res_D = results_1[1][4]
res_A = results_1[1][5]
cases_age = zeros(17)
cases_age_A = zeros(17)
for a = 1:17
    cases_age[a] = median([sum(res_D[:,a,k]) for k = 1:1000])
end
for a = 1:17
    cases_age_A[a] = median([sum(res_A[:,a,k]) for k = 1:1000])
end
sum(cases_age_A .+ cases_age)
cases_age./(cases_age+ cases_age_A)

plt_age_distrib_uncontrolled = plot_total_incidence_by_age(scenario_group,treatment_rates,rel_transmission_perc,1)
plot!(size = (800,500))
xlabel!("Age of infected individual")
ylabel!("Number of symptomatic infections (thousands)")
title!("Age distribution of cases: uncontrolled epidemic")
savefig(plt_age_distrib_uncontrolled,"plotting/age_distrib_cases_uncontrolled.pdf")

plt_age_distrib_controlled = plot_total_incidence_by_age(scenario_group,treatment_rates,rel_transmission_perc,6)
plot!(size = (800,500),legend = :topleft)

xlabel!("Age of case individual")
ylabel!("Number of cases (millions)")
title!("Age distribution of cases: targetted controls")
savefig(plt_age_distrib_controlled,"plotting/age_distrib_cases_targetted.pdf")

scenario_group_SD = [results_1,results_3,results_3SI]
plt_age_distrib_SD = plot_total_incidence_by_age(scenario_group_SD,1,
                    ["0% rel. infectiousness" "25% rel. infectiousness" "25% rel. infectiousness + SD"])
plot!(size = (800,500))
xlabel!("Age of case individual")
ylabel!("Number of cases (thousands)")
title!("Age profile of cases")
savefig(plt_age_distrib_SD,"plotting/age_distrib_cases_selected.pdf")


text:
const DEMCMC = DifferentialEvolutionMCMC

"""
    sample(model::DEModel, de::DE, n_iter::Int; progress=false, kwargs...)

Samples from the posterior distribution.

# Arguments

- `model`: a model containing likelihood function with data and priors
- `de`: differential evolution object
- `n_iter`: number of iterations or samples

# Keywords

- `progress=false`: show progress of sampler
- `kwargs...`: optional keyword arguments
"""
sample(model::DEModel, de::DE, n_iter::Int; progress=false, kwargs...) = _sample(model::DEModel, de::DE, n_iter::Int; progress, stepfun=step!, kwargs...)

function _sample(model::DEModel, de::DE, n_iter::Int; progress=false, stepfun=step!, kwargs...)
    meter = Progress(n_iter)
    # initialize particles based on prior distribution
    groups = sample_init(model, de, n_iter)
    for iter in 1:n_iter
        de.iter = iter + de.n_initial
        # explicitly pass groups so parallel works
        groups = stepfun(model, de, groups)
        progress ? next!(meter) : nothing
    end
    # convert to chain object
    chain = bundle_samples(model, de, groups, n_iter)
    return chain
end

"""
    sample(model::DEModel, de::DE, ::MCMCThreads, n_iter::Int; progress=false, kwargs...)

Samples from the posterior distribution with each group of particles on a seperarate thread for
the mutation and crossover steps.

# Arguments

- `model`: a model containing likelihood function with data and priors
- `de`: differential evolution object
- `MCMCThreads`: pass MCMCThreads() object to run on multiple threads
- `n_iter`: number of iterations or samples

# Keywords

- `progress=false`: show progress of sampler
- `kwargs...`: optional keyword arguments
"""
function sample(model::DEModel, de::DE, ::MCMCThreads, n_iter::Int; progress=false, kwargs...)
    _sample(model::DEModel, de::DE, n_iter::Int; progress, stepfun=pstep!, kwargs...)
end

"""
    step!(model::DEModel, de::DE, groups)

Perform a single step for DE-MCMC.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: DE-MCMC sampler object
- `groups`: Array of vectors of particles
"""
function step!(model::DEModel, de::DE, groups)
    rand() <= de.α ? migration!(de, groups) : nothing
    groups = update!(model, de, groups)
    store_samples!(de, groups)
    return groups
end

"""
    pstep!(model::DEModel, de::DE, groups)

Perform a single step for DE-MCMC with each particle group on a separate thread.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: DE-MCMC sampler object
- `groups`: Array of vectors of particles
"""
function pstep!(model::DEModel, de::DE, groups)
    rand() <= de.α ? migration!(de, groups) : nothing
    groups = p_update!(model, de, groups)
    store_samples!(de, groups)
    return groups
end

"""
    p_update!(model, de, groups)

Multithreaded update of particles. Particles are updated by block or simultaneously.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a vector of interacting particles (e.g. chains)
"""
function p_update!(model, de, groups)
    seeds = rand(UInt, length(groups))
    if de.blocking_on(de)
        Threads.@threads for g in 1:de.n_groups
            group = block_update!(model, de, groups[g], seeds[g])
        end
        return groups
    else
        Threads.@threads for g in 1:de.n_groups
            group = mutate_or_crossover!(model, de, groups[g], seeds[g])
        end
        return groups
    end
end

"""
    update!(model, de, groups)(model, de, groups)

Particles are updated by block or simultaneously on a single thread.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a vector of interacting particles (e.g. chains)
"""
function update!(model, de, groups)
    if de.blocking_on(de)
        return map(g -> block_update!(model, de, g), groups)
    else
        return map(g -> mutate_or_crossover!(model, de, g), groups)
    end
end

function block_update!(model, de, group, seed)
    Random.seed!(seed)
    return block_update!(model, de, group)
end

function block_update!(model, de, group)
    for block_idx in de.blocks
        mutate_or_crossover!(model, de, group, block_idx)
    end
    return group
end

"""
    mutate_or_crossover!(model, de, group, seed::Number)

Selects between mutation and crossover step with probability β.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a vector of interacting particles (e.g. chains)
- `seed::Number`: RNG seed
"""
function mutate_or_crossover!(model, de, group, seed::Number)
    Random.seed!(seed)
    mutate_or_crossover!(model, de, group)
    return group
end

function mutate_or_crossover!(model, de, group)
    rand() <= de.β ? mutation!(model, de, group) : crossover!(model, de, group)
    return group
end

function mutate_or_crossover!(model, de, group, block_idx)
    rand() <= de.β ? mutation!(model, de, group) : crossover!(model, de, group, block_idx)
    return group
end

"""
    bundle_samples(model::DEModel, de::DE, groups, n_iter)

Converts group particles to a chain object capable of generating convergence diagnostics
and posterior summaries.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `groups`: a vector of groups of particles
- `n_iter`: number of iterations
"""
function bundle_samples(model::DEModel, de::DE, groups, n_iter)
    @unpack burnin, discard_burnin = de
    particles = vcat(groups...)
    Np = length(particles)
    Ns = discard_burnin ? n_iter - burnin : n_iter
    offset = discard_burnin ? burnin : 0
    all_names = get_names(model, particles[1])
    n_parms = length(all_names)
    n_names = length(model.names)
    v = fill(0.0, Ns, n_parms, Np)
    for (c,p) in enumerate(particles)
        for s in 1:Ns
            sΔ = s + offset
            temp = Float64[]
            for ni in 1:n_names
                push!(temp, de.samples[sΔ,ni,c]...)
            end
            push!(temp, p.accept[sΔ], p.lp[sΔ])
            v[s,:,c] = temp'
        end
    end
    chains = Chains(v, all_names, (parameters = [model.names...],
        internals = ["acceptance", "lp"]))
    return chains
end

"""
    sample_init(model::DEModel, de::DE, n_iter)

Creates vectors of particles and samples initial parameter values from priors.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `n_iter`: number of iterations
"""
function sample_init(model::DEModel, de::DE, n_iter)
    de.samples = initialize_samples(de, model, n_iter)
    id = 0
    groups = [[init_particle(de, model, id+=1, n_iter) for p in 1:de.Np]
        for c in 1:de.n_groups]
    return groups
end


text:
using DashBootstrapComponents, DashHtmlComponents

text_input = html_div([
    dbc_input(id = "input", placeholder = "Type something...", type = "text"),
    html_br(),
    html_p(id = "output"),
]);

callback!(app, Output("output", "children"), Input("input", "value")) do value
    return value
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module OutlierDetectionNetworks
    using OutlierDetectionInterface
    using OutlierDetectionInterface:SCORE_UNSUPERVISED, SCORE_SUPERVISED
    const OD = OutlierDetectionInterface

    include("utils.jl")
    include("models/ae.jl")
    include("models/dsad.jl")
    include("models/esad.jl")
    include("templates/templates.jl")

    const UUID = "c7f57e37-4fcb-4a0b-a36c-c2204bc839a7"
    const MODELS = [:AEDetector,
                    :DSADDetector,
                    :ESADDetector]

    for model in MODELS
        @eval begin
            OD.@default_frontend $model
            OD.@default_metadata $model $UUID
            export $model
        end
    end
end

max_stars_repo_name: davnn/OutlierDetectionNetworks

max_issues_repo_name: davnn/OutlierDetectionNetworks

max_forks_repo_name: davnn/OutlierDetectionNetworks