text:
""" SparseIndex """
struct SparseIndex{B<:BoolType,I}
    issparse::B
    index::I
end

@inline function unsafe_getvalue(x, i::SparseIndex{Bool})
    if getfield(i, :issparse)
        return zero(eltype(x))
    else
        return unsafe_getvalue(x, getfield(i, :index))
    end
end

"""
    SparseLinear{T}(nonzeros::Vector{T})

Transforms to sparsely populated linear space.
"""
struct SparseLinear{T<:Integer} <: CoordinateTransform{1,1}
    nonzeros::Vector{T}

    SparseLinear(x::AbstractSparseVector{<:Any,T}) where {T} = new{T}(nonzeros(x))
end

SparseArrays.nonzeros(x::SparseLinear) = getfield(x, :nonzeros)
SparseArrays.nnz(x::SparseLinear) = length(nonzeros(x))

@inline function (t::SparseLinear{T})(i::CanonicalInt) where {T}
    nzind = nonzeros(t)
    ii = searchsortedfirst(nzind, convert(T, i))
    if ii <= length(nzind) && (nzind[ii] == i)
        return SparseIndex(false, @inbounds(nzval[ii]))
    else
        return SparseIndex(true, ii)
    end
end

"""
    CompactSparseColumn{T}(colptr::Vector{T}, rowval::Vector{T})

Transforms two dimensional coordinates to linearly stored values, representing a compact
sparse column matrix.
"""
struct CompactSparseColumn{T<:Integer} <: CoordinateTransform{2,1}
    colptr::Vector{T}      # Column i is in colptr[i]:(colptr[i+1]-1)
    rowval::Vector{T}      # Row indices of stored values

    global function _CompactSparseColumn(m::Int, n::Int, c::Vector{T}, r::Vector{T}) where {T}
        new{T}(m, n, c, r)
    end
    CompactSparseColumn{T}() where {T} = new{T}(T[1], T[])
    function CompactSparseColumn{T}(n::Integer) where {T}
        (n < 0) && throw(ArgumentError("size of SquareCompactSparseColumn cannot be negative"))
        return new{T}(fill(one(T), n + 1), T[], m)
    end
    function CompactSparseColumn{T}(c::Vector{T}, r::Vector{T}) where {T}
        if 0 ≤ n && (!isbitstype(Ti) || n ≤ typemax(Ti))
            throw(ArgumentError("number of columns (c = $n) does not fit in eltype = $(T)"))
        end
        if !((c[end] - 1) === length(r))
            throw(ArgumentError("Illegal buffers for AdjacencyMap construction $n $colptr $rowval $nzval"))
        end
        return new{T}(c, r)
    end
    function CompactSparseColumn(x::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
        n, m = size(x)
        n === m || throw(ArgumentError("cannot derive SquareCompactSparseColumn from a matrix with unequal rows and columns"))
        return new{Ti}(copy(getcolptr(x)), copy(rowvals(x)))
    end
end

SparseArrays.getcolptr(x::CompactSparseColumn) = getfield(x, :colptr)
SparseArrays.rowvals(x::CompactSparseColumn) = getfield(x, :rowval)
SparseArrays.nnz(x::CompactSparseColumn) = Int(getcolptr(x)[end]) - 1

@inline function (t::CompactSparseColumn)(index::AbstractCartesianIndex{2})
    i0, i1 = Tuple(index)
    r1 = Int(@inbounds(getcolptr(t)[i1]))
    r2 = Int(@inbounds(getcolptr(t)[i1+1])-1)
    if r1 > r2
        return SparseIndex(true, r1)
    else
        r1 = searchsortedfirst(rowvals(t), i0, r1, r2, Base.Forward)
        if ((r1 > r2) || (@inbounds(rowvals(t)[r1]) != i0))
            return SparseIndex(true, r1)
        else
            return SparseIndex(false, r1)
        end
    end
end


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

#

export Workspace

"""
    struct Workspace … end

A collection of variables.
"""
struct Workspace
    contents::OrderedDict{Symbol,Any}
    Workspace() = new(OrderedDict{Symbol,Any}())
end

@inline _c(w::Workspace) = getfield(w, :contents)

Base.propertynames(w::Workspace) = tuple(keys(_c(w))...)
Base.getproperty(w::Workspace, sym::Symbol) = getindex(_c(w), sym)
Base.setproperty!(w::Workspace, sym::Symbol, val) = setindex!(_c(w), val, sym)

Base.getindex(w::Workspace, Args...) = getindex(_c(w), Args...)
Base.setindex!(w::Workspace, Args...) = setindex!(_c(w), Args...)

@inline Base.isempty(w::Workspace) = isempty(_c(w))
@inline Base.in(name, w::Workspace) = Symbol(name) ∈ keys(_c(w))

Base.keys(w::Workspace) = keys(_c(w))
Base.values(w::Workspace) = values(_c(w))
Base.iterate(w::Workspace, args...) = iterate(_c(w), args...)

function Base.show(io::IO, ::MIME"text/plain", w::Workspace)

    if isempty(w)
        return print(io, "Empty Workspace")
    end

    nvars = length(_c(w))
    println(io, "Workspace with ", nvars, "-variables")

    limit = get(io, :limit, true)
    io = IOContext(io, :SHOWN_SET => _c(w),
        :typeinfo => eltype(_c(w)),
        :compact => get(io, :compact, true),
        :limit => limit)

    dheight, dwidth = displaysize(io)

    if limit && nvars + 5 > dheight
        # we're printing some but not all rows (no room on the screen)
        top = div(dheight - 5, 2)
        bot = nvars - dheight + 7 + top
    else
        top, bot = nvars + 1, nvars + 1
    end

    max_align = 0
    prows = Vector{String}[]
    for (i, (k, v)) ∈ enumerate(_c(w))
        top < i < bot && continue

        sk = sprint(print, k, context = io, sizehint = 0)
        sv = sprint(summary, v, context = io, sizehint = 0)
        max_align = max(max_align, length(sk))

        push!(prows, [sk, sv])
        i == top && push!(prows, ["⋮", "⋮"])
    end

    cutoff = dwidth - 5 - max_align

    for (sk, sv) ∈ prows
        lv = length(sv)
        sv = lv <= cutoff ? sv : sv[1:cutoff-1] * "…"
        println(io, "  ", lpad(sk, max_align), " ⇒ ", sv)
    end

end

ext: jl


text:
using ULID
using Test

@testset "Pseudorandom number generation" begin
    @test isa(ULID.prng(), Float32) || isa(ULID.prng(), Float64)
    @test isfinite(ULID.prng())
    @test 0 < ULID.prng() < 1
end

@testset "Encoding based on a given system time" begin
    let t = 1469918176385
        @test ULID.encodetime(t, 10) == "01ARYZ6S41"
        @test ULID.encodetime(t, 12) == "0001ARYZ6S41"
        @test ULID.encodetime(t, 8) == "ARYZ6S41"
    end
end

@testset "Encode to a given length" begin
    @test length(ULID.encoderandom(12)) == 12
end

@testset "The reason we're here" begin
    @test length(ulid()) == 26
end


text:
import IncrementalInference: selectFactorType, getDomain

export getDomain

## =============================================================================
## Needs a home
## =============================================================================

# Still experimental
# export BearingRange2
@defVariable BearingRange2 2 (:Circular, :Euclid)

## ============================================================================
# Starting integration with Manifolds.jl, via ApproxManifoldProducts.jl first
## ============================================================================

getDomain(::InstanceType{Point2Point2}) = Point2
getDomain(::InstanceType{Pose2Point2}) = Point2
getDomain(::InstanceType{Pose2Pose2}) = Pose2
# getDomain(::InstanceType{Pose3Point3}) = Point3
getDomain(::InstanceType{Pose3Pose3}) = Pose3
getDomain(::InstanceType{Pose2Point2BearingRange}) = BearingRange2

# defines on which variable manifold the measurements for a factor type are represented
getManifolds(fctType::InstanceType{Pose2Pose2}) = getManifolds(getDomain(fctType))
getManifolds(fctType::InstanceType{Pose2Point2}) = getManifolds(getDomain(fctType))
getManifolds(fctType::InstanceType{Point2Point2}) = getManifolds(getDomain(fctType))
getManifolds(fctType::InstanceType{Pose2Point2BearingRange}) = getManifolds(getDomain(fctType))

function selectFactorType(T1::Type{<:InferenceVariable}, T2::Type{<:InferenceVariable})
  # hacky version
  if T1 == Pose2 && T2 == Pose2
    return Pose2Pose2
  elseif T1 == Pose2 && T2 == Point2
    return Pose2Point2
  elseif T1 == DynPose2 && T2 == DynPose2
      return DynPose2DynPose2
  elseif T1 == Point2 && T2 == Point2
    return Point2Point2
  elseif T1 == Point3 && T2 == Point3
    return Point3Point3
  elseif T1 == Pose3 && T2 == Pose3
    return Pose3Pose3
  else
    error("dont know which Factor type to select between $T1 and $T2")
  end
end

# FIXME, much consolidation required here, see RoME #244
convert(::Type{<:AMP.Manifold}, ::InstanceType{Point2}) = AMP.Euclid2
convert(::Type{<:AMP.Manifold}, ::InstanceType{Pose2}) = AMP.SE2_Manifold

convert(::Type{<:AMP.Manifold}, ::InstanceType{Point2Point2}) = AMP.Euclid2
convert(::Type{<:AMP.Manifold}, ::InstanceType{Pose2Pose2}) = AMP.SE2_Manifold


text:
# Diffusion maps
# --------------
# Diffusion maps,
# Coifman, R. & Lafon, S., Applied and Computational Harmonic Analysis, Elsevier, 2006, 21, 5-30

#### DiffMap type
struct DiffMap{T <: Real} <: AbstractDimensionalityReduction
    t::Int
    ɛ::T
    K::AbstractMatrix{T}
    proj::Projection{T}

    DiffMap{T}(t::Int, ɛ::T, K::AbstractMatrix{T}, proj::Projection{T})  where T = new(t, ɛ, K, proj)
end

## properties
outdim(R::DiffMap) = size(R.proj, 1)

## custom
kernel(R::DiffMap) = R.K

## show
summary(io::IO, R::DiffMap) = print(io, "Diffusion Maps(outdim = $(outdim(R)), t = $(R.t), ɛ = $(R.ɛ))")
function show(io::IO, R::DiffMap)
    summary(io, R)
    io = IOContext(io, :limit=>true)
    println(io)
    println(io, "Kernel: ")
    Base.print_matrix(io, R.K, "[", ",","]")
    println(io)
    println(io, "Embedding:")
    Base.print_matrix(io, transform(R), "[", ",","]")
end

## interface functions
function fit(::Type{DiffMap}, X::AbstractMatrix{T}; maxoutdim::Int=2, t::Int=1, ɛ::Real=1.0) where {T<:Real}
    # rescale data
    Xtr = standardize(StatsBase.UnitRangeTransform, X)
    Xtr[findall(isnan, Xtr)] .= 0

    # compute kernel matrix
    sumX = sum(Xtr.^ 2, dims=1)
    K = exp.(-( transpose(sumX) .+ sumX .- 2*transpose(Xtr) * Xtr ) ./ convert(T, ɛ))

    p = transpose(sum(K, dims=1))
    K ./= (p * transpose(p)) .^ convert(T, t)
    p = transpose(sqrt.(sum(K, dims=1)))
    K ./= p * transpose(p)

    U, S, V = svd(K, full=true)
    U ./= U[:,1]
    Y = U[:,2:(maxoutdim+1)]

    return DiffMap{T}(t, convert(T, ɛ), K, transpose(Y))
end

transform(R::DiffMap) = R.proj


text:
# -*- mode: julia; -*-
module WAVPlay
import ..wavplay

import Libdl

# typedef enum pa_sample_format
const PA_SAMPLE_U8        =  0 # Unsigned 8 Bit PCM
const PA_SAMPLE_ALAW      =  1 # 8 Bit a-Law
const PA_SAMPLE_ULAW      =  2 # 8 Bit mu-Law
const PA_SAMPLE_S16LE     =  3 # Signed 16 Bit PCM, little endian (PC)
const PA_SAMPLE_S16BE     =  4 # Signed 16 Bit PCM, big endian
const PA_SAMPLE_FLOAT32LE =  5 # 32 Bit IEEE floating point, little endian (PC), range -1.0 to 1.0
const PA_SAMPLE_FLOAT32BE =  6 # 32 Bit IEEE floating point, big endian, range -1.0 to 1.0
const PA_SAMPLE_S32LE     =  7 # Signed 32 Bit PCM, little endian (PC)
const PA_SAMPLE_S32BE     =  8 # Signed 32 Bit PCM, big endian
const PA_SAMPLE_S24LE     =  9 # Signed 24 Bit PCM packed, little endian (PC). \since 0.9.15
const PA_SAMPLE_S24BE     = 10 # Signed 24 Bit PCM packed, big endian. \since 0.9.15
const PA_SAMPLE_S24_32LE  = 11 # Signed 24 Bit PCM in LSB of 32 Bit words, little endian (PC). \since 0.9.15
const PA_SAMPLE_S24_32BE  = 12 # Signed 24 Bit PCM in LSB of 32 Bit words, big endian. \since 0.9.15

struct pa_sample_spec
    format::Int32
    rate::UInt32
    channels::UInt8
end

struct pa_channel_map
    channels::UInt8

    # map data (max 32 channels)
    map0::Cint
    map1::Cint
    map2::Cint
    map3::Cint
    map4::Cint
    map5::Cint
    map6::Cint
    map7::Cint
    map8::Cint
    map9::Cint
    map10::Cint
    map11::Cint
    map12::Cint
    map13::Cint
    map14::Cint
    map15::Cint
    map16::Cint
    map17::Cint
    map18::Cint
    map19::Cint
    map20::Cint
    map21::Cint
    map22::Cint
    map23::Cint
    map24::Cint
    map25::Cint
    map26::Cint
    map27::Cint
    map28::Cint
    map29::Cint
    map30::Cint
    map31::Cint

    pa_channel_map() = new(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
end

struct pa_buffer_attr
    maxlength::UInt32
    tlength::UInt32
    prebuf::UInt32
    minreq::UInt32
    fragsize::UInt32
end

const pa_simple = Ptr{Cvoid}
const LibPulseSimple = Libdl.find_library(["libpulse-simple", "libpulse-simple.so.0"])
const PA_STREAM_PLAYBACK = 1
const PA_CHANNEL_MAP_AIFF = 0
const PA_CHANNEL_MAP_DEFAULT = PA_CHANNEL_MAP_AIFF

function wavplay(data, fs)
    nChannels = size(data,2)
    ss = pa_sample_spec(PA_SAMPLE_FLOAT32LE, fs, nChannels)

    # Manually layout the samples.
    # convert doesn't lay out the samples as pulse audio expects
    samples = Array{Float32, 1}(undef, size(data, 1) * size(data, 2))
    idx = 1
    for i = 1:size(data, 1)
        for j = 1:size(data, 2)
            samples[idx] = convert(Float32, data[i, j])
            idx += 1
        end
    end

    s = ccall((:pa_simple_new, LibPulseSimple),
              pa_simple,
              (Cstring,
               Cstring,
               Cint,
               Cstring,
               Cstring,
               Ptr{pa_sample_spec},
               Ptr{pa_channel_map},
               Ptr{pa_buffer_attr},
               Ptr{Cint}),
              C_NULL, # Use the default server
              "Julia WAV.jl",  # Application name
              PA_STREAM_PLAYBACK,
              C_NULL, # Use the default device
              "wavplay", # description of stream
              Ref(ss),
              C_NULL, # Use default channel map
              C_NULL, # Use default buffering attributes
              C_NULL) # Ignore error code
    if s == C_NULL
        error("pa_simple_new failed")
    end

    write_ret = ccall((:pa_simple_write, LibPulseSimple),
                      Cint,
                      (pa_simple, Ptr{Cvoid}, Csize_t, Ptr{Cint}),
                      s, samples, sizeof(samples), C_NULL)
    if write_ret != 0
        error("pa_simple_write failed with $write_ret")
    end

    drain_ret = ccall((:pa_simple_drain, LibPulseSimple),
                      Cint,
                      (pa_simple, Ptr{Cint}), s, C_NULL)
    if drain_ret != 0
        error("pa_simple_drain failed with $drain_ret")
    end

    ccall((:pa_simple_free, LibPulseSimple), Cvoid, (pa_simple,), s)
end
end # module


text:
@testset "Containers" begin
    dv = @data ones(3)
    push!(dv, 3.0)
    push!(dv, missing)

    @test isequal(dv, (@data [1.0, 1.0, 1.0, 3.0, missing]))

    a, b = pop!(dv), pop!(dv)
    @test ismissing(a)
    @test b == 3.0

    pushfirst!(dv, 3.0)
    pushfirst!(dv, missing)

    @test isequal(dv, (@data [missing, 3.0, 1.0, 1.0, 1.0]))

    a, b = popfirst!(dv), popfirst!(dv)
    @test ismissing(a)
    @test b == 3.0

    ## SPLICE
    function test_splice(dv, spliceout, splicein...)
        dv = copy(dv)
        v = Any[x for x in dv]
        expval = dv[spliceout]
        retval = splice!(dv, spliceout, splicein...)
        isequal(retval, expval) ||
            error("splice!$(tuple(dv, spliceout, splicein...)) gave incorrect return value $retval (expected $expval)")
        splice!(v, spliceout, splicein...)
        for i = 1:length(v)
            isequal(dv[i], v[i]) ||
                error("splice!$(tuple(dv, spliceout, splicein...)) gave incorrect output")
        end
    end

    function test_deleteat(dv, spliceout)
        dv = copy(dv)
        v = Any[x for x in dv]
        @test deleteat!(dv, spliceout) === dv
        deleteat!(v, spliceout)
        for i = 1:length(v)
            isequal(dv[i], v[i]) ||
                error("deleteat!($(tuple(dv, spliceout)) gave incorrect output")
        end
    end

    dv1 = @data [1.0, 2.0, missing, 2.0, missing, 3.0]
    for dv in (dv1, convert(DataVector{Number}, dv1), convert(PooledDataArray, dv1))
        for spliceout in (2, 3, 2:3, 5:6)
            test_splice(dv, spliceout)
            test_deleteat(dv, spliceout)
            for splicein in ([], [3], @data([3]), @pdata([3]),
                             [3, 4, 5], [3., 4., 5.], @data([3, missing, 4]),
                             @pdata([3, missing, 4]), @data([missing, 3.0, 4.0]),
                             @pdata([missing, 3.0, 4.0]))
                test_splice(dv, spliceout, splicein)
            end
        end
    end

    ## sizehint
    sizehint!(@data([1.0, 2.0]), 5)
    sizehint!(@pdata([1.0, 2.0]), 5)
end


text:
using Distributed
using Base.Filesystem
using DataFrames
using CSV
using Query
using Statistics
using UnicodePlots
using ClusterManagers
using Dates
using DelimitedFiles

## load the packages by covid19abm
using Parameters, Distributions, StatsBase, StaticArrays, Random, Match, DataFrames
using covid19abm
const cv=covid19abm

#addprocs(30, exeflags="--project=.")
#@everywhere using covid19abm

addprocs(SlurmManager(512), N=16, topology=:master_worker, exeflags="--project=.")
@everywhere using covid19abm

function run(myp::ModelParameters, nsims=500, folderprefix="./")
    println("starting $nsims simulations...\nsave folder set to $(folderprefix)")
    dump(myp)
    myp.calibration && error("can not run simulation, calibration is on.")
    # will return 6 dataframes. 1 total, 4 age-specific
    cdr = pmap(1:nsims) do x
            cv.runsim(x, myp)
    end

    println("simulations finished")
    println("total size of simulation dataframes: $(Base.summarysize(cdr))")
    ## write the infectors
    writedlm("$(folderprefix)/infectors.dat", [cdr[i].infectors for i = 1:nsims])

    ## write contact numbers
    writedlm("$(folderprefix)/ctnumbers.dat", [cdr[i].ct_numbers for i = 1:nsims])
    ## stack the sims together
    allag = vcat([cdr[i].a  for i = 1:nsims]...)
    ag1 = vcat([cdr[i].g1 for i = 1:nsims]...)
    ag2 = vcat([cdr[i].g2 for i = 1:nsims]...)
    ag3 = vcat([cdr[i].g3 for i = 1:nsims]...)
    ag4 = vcat([cdr[i].g4 for i = 1:nsims]...)
    ag5 = vcat([cdr[i].g5 for i = 1:nsims]...)
    #mydfs = Dict("all" => all, "ag1" => ag1, "ag2" => ag2, "ag3" => ag3, "ag4" => ag4, "ag5" => ag5)
    mydfs = Dict("all" => allag)

    ## save at the simulation and time level
    ## to ignore for now: miso, iiso, mild
    c1 = Symbol.((:LAT, :ASYMP, :INF, :IISO, :HOS, :ICU, :DED), :_INC)
    c2 = Symbol.((:LAT, :ASYMP, :INF, :IISO, :HOS, :ICU, :DED), :_PREV)
    for (k, df) in mydfs
        println("saving dataframe sim level: $k")
        # simulation level, save file per health status, per age group
        for c in vcat(c1..., c2...)
            udf = unstack(df, :time, :sim, c)
            fn = string("$(folderprefix)/simlevel_", lowercase(string(c)), "_", k, ".dat")
            CSV.write(fn, udf)
        end
        println("saving dataframe time level: $k")
        # time level, save file per age group
        yaf = compute_yearly_average(df)
        fn = string("$(folderprefix)/timelevel_", k, ".dat")
        CSV.write(fn, yaf)
    end
    return mydfs
end

function compute_yearly_average(df)
    ya = df |> @groupby(_.time) |> @map({time=key(_), cnt=length(_),
              sus_prev=mean(_.SUS_PREV),
              lat_prev=mean(_.LAT_PREV),
              pre_prev=mean(_.PRE_PREV),
              asymp_prev=mean(_.ASYMP_PREV),
              mild_prev=mean(_.MILD_PREV),
              miso_prev=mean(_.MISO_PREV),
              inf_prev=mean(_.INF_PREV),
              iiso_prev=mean(_.IISO_PREV),
              hos_prev=mean(_.HOS_PREV),
              icu_prev=mean(_.ICU_PREV),
              rec_prev=mean(_.REC_PREV),
              ded_prev=mean(_.DED_PREV),
              sus_inc=mean(_.SUS_INC),
              lat_inc=mean(_.LAT_INC),
              pre_inc=mean(_.PRE_INC),
              asymp_inc=mean(_.ASYMP_INC),
              mild_inc=mean(_.MILD_INC),
              miso_inc=mean(_.MISO_INC),
              inf_inc=mean(_.INF_INC),
              iiso_inc=mean(_.IISO_INC),
              hos_inc=mean(_.HOS_INC),
              icu_inc=mean(_.ICU_INC),
              rec_inc=mean(_.REC_INC),
              ded_inc=mean(_.DED_INC)
              }) |> DataFrame
    return ya
end

function savestr(p::ModelParameters, custominsert="/", customstart="")
    datestr = (Dates.format(Dates.now(), dateformat"mmdd_HHMM"))
    ## setup folder name based on model parameters
    taustr = replace(string(p.τmild), "." => "")
    fstr = replace(string(p.fmild), "." => "")
    rstr = replace(string(p.β), "." => "")
    prov = replace(string(p.prov), "." => "")
    eldr = replace(string(p.eldq), "." => "")
    eldqag = replace(string(p.eldqag), "." => "")
    #fpre = replace(string(p.fpre), "." => "")
    fasymp = replace(string(p.fasymp), "." => "")
    fpreiso = replace(string(p.fpreiso), "." => "")
    tpreiso = replace(string(p.tpreiso), "." => "")
    fsev = replace(string(p.fsevere), "." => "")
    frelasymp = replace(string(p.frelasymp), "." => "")
    pct = replace(string(p.fctcapture), "." => "")
    cct = replace(string(p.fcontactst), "." => "")
    idt = replace(string(p.cidtime), "." => "")
    tback = replace(string(p.cdaysback), "." => "")
    #fldrname = "/data/covid19abm/simresults/$(custominsert)/b$(rstr)_$(prov)_pct$(pct)_cct$(cct)_idt$(idt)_tback$(tback)_fsev$(fsev)_tau$(taustr)_fmild$(fstr)_q$(eldr)_qag$(eldqag)_relasymp$(frelasymp)_asymp$(fasymp)_tpreiso$(tpreiso)_preiso$(fpreiso)/"
    fldrname = "/data/covid19abm/simresults/$(custominsert)/$(customstart)$(prov)_pct$(pct)_cct$(cct)_idt$(idt)_tback$(tback)_fsev$(fsev)_tau$(taustr)_fmild$(fstr)_q$(eldr)_qag$(eldqag)_relasymp$(frelasymp)_asymp$(fasymp)_tpreiso$(tpreiso)_preiso$(fpreiso)/"
    mkpath(fldrname)
end

function _calibrate(nsims, myp::ModelParameters)
    myp.calibration != true && error("calibration parameter not turned on")
    vals = zeros(Int64, nsims)
    println("calibrating with beta: $(myp.β), total sims: $nsims, province: $(myp.prov)")
    println("calibration parameters:")
    dump(myp)
    cdr = pmap(1:nsims) do i
        h = main(myp) ## gets the entire model.
        val = sum(_get_column_incidence(h, covid19abm.LAT))
        return val
    end
    return mean(cdr), std(cdr)
end

function calibrate(beta, nsims, prov=:ontario)
    myp = ModelParameters()
    myp.β = beta
    myp.prov = prov
    myp.calibration = true
    myp.fmild = 0.0
    myp.fsevere = 0.0
    myp.fpreiso = 0.0
    myp.fasymp = 0.5
    myp.initialinf = 1
    m, sd = _calibrate(nsims, myp)
    println("mean R0: $(m) with std: $(sd)")
    myp.calibration = false
    return m
end

function calibrate_robustness(beta, reps, prov=:ontario)
    #[:ontario, :alberta, :bc, :manitoba, :newbruns, :newfdland, :nwterrito, :novasco, :nunavut, :pei, :quebec, :saskat, :yukon]
    # once a beta is found based on nsims simulations,
    # see how robust it is. run calibration with same beta 100 times
    # to see the variation in R0 produced.
    nsims = [500, 1000]
    means = zeros(Float64, reps, length(nsims))
    for (i, ns) in enumerate(nsims)
        cd = map(1:reps) do x
            println("iter: $x, sims: $ns")
            mval = calibrate(beta, ns, prov)
            return mval
        end
        means[:, i] = cd
    end
    # for i in 2:nworkers()
    #     ## mf defined as: @everywhere mg() = covid19abm.p.β
    #     rpr = remotecall_fetch(mf,  i+1).prov
    #     rpr != prov && error("province didn't get set in the remote workers")
    # end
    return means
end
#provs = []
#for x in [:ontario, :alberta, :bc, :manitoba, :newbruns]
#    m = calibrate_robustness(0.101, x)
#    push!(provs, m)
#end

#@everywhere mg() = covid19abm.p.β
#@everywhere mf(b) = covid19abm.p.β = b
# function hmatrix_analysis()
#     ## function takes the output of the model
#     ## and splits it into an array of 2d matrices for plotting purposes.
#     myp = cv.ModelParameters()
#     myp.β = 0.038
#     reset_params($myp)
#     cd, ags = main(1)

#     #mats = Array{Int64, 3}(undef, 500)
#     mats = Array{Int64, 3}(undef, 100,100,500)
#     for i = 1:500
#         col = cd[:, i]
#         mats[:, :, i] = reshape(col, (100, 100))
#     end

#     colorarray = IndirectArray(mats, [colorant"green", colorant"teal", colorant"blue", colorant"azure",
#     colorant"red", colorant"lightsalmon", colorant"magenta", colorant"magenta", colorant"grey60", colorant"black"])
#     save("myanimation.gif", colorarray; fps=5)

#     idxarray = rand(1:3, 100, 100, 50);
#     colorarray = IndirectArray(idxarray, [colorant"red", colorant"purple", colorant"indigo"])

#     return mats
# end

# function makeheatmaps(bufs)
#     heatmaps = map(bufs) do buf
#         heatmap(
#             buf, scale_plot = false, show_axis = false)

#     end
#     scene = hbox(map(i-> vbox(heatmaps[i, :]), 1:size(bufs, 1))...)
#     scene, last.(heatmaps)
# end
#datarows = 500; datacols = 500
#plotrows = 4; plotcols = 4
#bufs = [fill(0.0f0, datarows, datacols) for _ in 1:plotrows, _ in 1:plotcols]

# http://juliaplots.org/MakieReferenceImages/gallery//lots_of_heatmaps/index.html
# http://juliaplots.org/MakieReferenceImages/gallery//chess_game/index.html


text:
__precompile__()
module NeighbourLists

const MAX_THREADS = [1]   # temporarily hard-coded single thread!
set_maxthreads!(n) = (NeighbourLists.MAX_THREADS[1] = n)

include("types.jl")

# this contains the cell-list data structures and assembly
include("cell_list.jl")

# this contains the different iterators over sites, bonds, etc
include("iterators.jl")

# alternative assembly protocol more akin to mapreduce
include("mapreduce.jl")

end # module


text:
# models for OSM data

struct Way
    id::Int64
    nodes::Vector{Int64}
    tags::Dict{String, String}
end

struct Node
    id::Int64
    lat::Float64
    lon::Float64
    tags::Dict{String, String}
end

@enum RelationMemberType way node relation

struct RelationMember
    id::Int64
    type::RelationMemberType
    role::String
end

struct Relation
    id::Int64
    members::Vector{RelationMember}
    tags::Dict{String, String}
end

# used by load_pbf but not by scan_pbf
struct OSMFile
    ways::Dict{Int64, Way}
    nodes::Dict{Int64, Node}
    relations::Dict{Int64, Node}
end


text:
module Devectorize

export
    # fun_traits
    TFun,
    TCallSig,
    result_type,
    sqr, rcp, blend,

    # texpr

    DeError,

    TExpr, TEWise, TScalar, TGeneralVar, TFunCall,
    TEmpty, TNum, TScalarVar, TVar, TQVar, TGeneralScalar,
    TRef, TIndex, TRange, TColon, TInterval,
    TGeneralRef1, TGeneralRef2, TRef1D, TRef2D, TRefRow, TRefCol,
    TMap, TReduc, TColwiseReduc, TRowwiseReduc,
    TLValue, TRValue, TAssign, TBlock,

    texpr, tnum, tscalarvar, tvar, tqvar,
    tref, tcall, tassign, topassign, tblock,
    is_trivial_assignment, ju_expr,

    tmode, ScalarMode, EWiseMode, ReducMode,
    ColwiseReducMode, RowwiseReducMode,

    # scalar_backend

    ScalarContext, dump_devec,
    compile, @devec, @inspect_devec,

    # extensions

    @devec_transform

import Base.==, Base.!=

include("fun_traits.jl")
include("texpr.jl")
include("meta_tools.jl")

include("compile_base.jl")
include("scalar_backend.jl")

include("extensions.jl")

end