text:
"""
```
function impulse_responses(m::AbstractDSGEVARModel{S}, paras::Matrix{S},
                           data::Matrix{S}, input_type::Symbol, method::Symbol;
                           parallel::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           n_obs_shock::Int = 1, draw_shocks::Bool = false,
                           flip_shocks::Bool = false,
                           X̂::AbstractMatrix{S} = Matrix{S}(undef, 0, 0),
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false, test_meansbands::Bool = false,
                           minimize::Bool = true,
                           forecast_string::String = "",
                           verbose::Symbol = :high) where {S<:Real}
```
computes the VAR impulse responses identified by the state space system
```
sₜ = TTT × sₜ₋₁ + RRR × impact[:, i],
yₜ = ZZ × sₜ + DD + MM × impact[:, i],
```
where `impact[:, i]` is a linear combination of
(orthogonal) structural shocks `ϵₜ ∼ 𝒩 (0, I)`, and
`MM × impact[:, i]` are the correlated measurement errors.

The VAR impulse responses are computed according to
```
ŷₜ₊₁ = X̂ₜ₊₁β + uₜ₊₁,
```
where `X̂ₜ₊₁` are the lags of observables in period `t + 1`, i.e. `yₜ, yₜ₋₁, ..., yₜ₋ₚ`.

If the method is `:rotation`, the shock `uₜ₊₁` is identified via
```
Σᵤ = 𝔼[u × u'] = chol(Σᵤ) × Ω × ϵₜ,
```
where the rotation matrix `Ω` is the `Q` matrix from a QR decomposition
of the impact response matrix corresponding to the state space system, i.e.
```
Ω, _ = qr(∂yₜ / ∂ϵₜ').
```

Otherwise, we draw a β and Σᵤ from the posterior implied by the DSGE
and data, and we then compute normal VAR impulse responses given those
coefficients and innovations variance-covariance matrix.

****
NOTE: this function generally involves taking random draws from
probability distributions, so seeds need to be set
to achieve reproducibility.
****

### Inputs
* `m::Union{AbstractDSGEModel,AbstractDSGEVARModel}`: DSGE/DSGEVAR model object
* `paras::Matrix{S}` or `paras::Vector{S}`: parameters to calibrate the model
* `input_type::Symbol`: `:mode` specifies a modal impulse response, and
    `:full` specifies a full-distribution forecast if `paras` is not given.
    This argument is also used to construct the file names of computed `MeansBands`.
* `method::Symbol`: type of impulse response to compute. The options are
    `:cholesky` and `:rotation`. For the first, see `?cholesky_shock`,
    and for the latter, we use the DSGE model to identify the rotation matrix
    which maps the DSGE's structural shocks to the innovations in the VAR's observables.
* `lags::Int`: number of lags in the VAR(p) approximation, i.e. p = lags
* `observables::Vector{Symbol}`: observables to be used in the VAR. These can be
    any of the observables or pseudo-observables in `m`.
* `shocks::Vector{Symbol}`: (structural) exogenous shocks to be used in the DSGE-VAR.
    These shocks must be in `m`.
* `n_obs_shock::Int`: the index of the observable corresponding to the orthogonalized shock causing the impulse response.

### Keywords
* `parallel::Bool`: use parallel workers or not
* `frequency_band::Tuple{S,S}`: See `?maxBC_shock`.
* `n_obs_shock::Int`: Index of observable to be shocked when using a Cholesky-based impulse response
* `draw_shocks::Bool`: true if you want to draw shocks along the entire horizon
* `flip_shocks::Bool`: impulse response shocks are negative by default. Set to `true` for
    a positive signed shock.
* `X̂::AbstractMatrix{S}`: matrix stacking the intercept and lags of the data for
    rotation IRFs. Set to a vector of zeros with length `1 + n_observables * p`
    to compute the rotation IRFs in deviations from the baseline forecast.
* `deviations::Bool`: set true to compute the impulse response in deviations
    rather than as a forecast. Mechnically, we ignore `X̂` (treated as zeros)
    and the intercept term.
* `normalize_rotation::Bool`: set to true to normalize the rotation
    so that rows have the correct sign. This requires as many structural shocks
    as there are observables in the DSGE-VAR.
* `density_bands::Vector{Float64}`: bands for full-distribution IRF computations
* `create_meansbands::Bool`: set to `true` to save output as a `MeansBands` object.
* `minimize::Bool`: choose shortest interval if true, otherwise just chop off lowest and
    highest (percent/2)
* `forecast_string::String`: string tag for identifying this impulse response
* `verbose::Symbol`: quantity of output desired

"""
function impulse_responses(m::AbstractDSGEVARModel{S}, paras::Matrix{S},
                           data::Matrix{S}, input_type::Symbol, method::Symbol;
                           parallel::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           n_obs_shock::Int = 1, draw_shocks::Bool = false,
                           flip_shocks::Bool = false,
                           X̂::AbstractMatrix{S} = Matrix{S}(undef, 0, 0),
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false, test_meansbands::Bool = false,
                           minimize::Bool = true,
                           forecast_string::String = "",
                           verbose::Symbol = :high) where {S<:Real}
    @assert !isempty(data) "A non-empty data matrix with dimensions (n_observables x n_periods) must be passed to use $(string(method))"

    if !(method in [:rotation, :cholesky, :cholesky_long_run, :choleskyLR, :maxBC,
                      :maximum_business_cycle_variance])
        error("A VAR IRF identified by a DSGE's rotation using method"
              * " $(string(method)) has not been implemented.")
    end

    # Compute dimensions of needed objects
    lags = get_lags(m)
    nobs = size(data, 1)
    k = nobs * lags + 1
    h = impulse_response_horizons(m)

    # Preompute X̂ and MM for the rotation identification
    if method == :rotation
        if isempty(X̂)
            XX = lag_data(data, lags; use_intercept = true)
            X̂ = vcat(1, data[:, end], XX[end, 1+1:k - nobs])
        end

        # Get measurement error
        if hasmethod(measurement_error, (typeof(m),))
            _, MM = measurement_error(m)
        else
            MM = zeros(S, n_observables(m), n_shocks(m))
        end
    end

    dsgevarrotationirf_method = if method == :rotation
        function _dsgevar_λ_rotation_irf_(para)
            DSGE.update!(m, para)
            return impulse_responses(m, data, X̂; horizon = h, MM = MM,
                                     flip_shocks = flip_shocks, draw_shocks = draw_shocks,
                                     deviations = deviations, normalize_rotation = normalize_rotation,
                                     verbose = verbose)
        end
    elseif method in [:cholesky, :cholesky_long_run, :choleskyLR, :maxBC,
                      :maximum_business_cycle_variance]
        function _dsgevar_λ_irf_(para)
            DSGE.update!(m, para)
            return impulse_responses(m, data, method, n_obs_shock; horizon = h,
                                     flip_shocks = flip_shocks,
                                     frequency_band = frequency_band)
        end
    end

    mapfcn = parallel ? pmap : map
    paras = mapslices(x -> [vec(x)], paras, dims = 2)
    irf_output =
        mapfcn(para -> dsgevarrotationirf_method(para), paras)

    if create_meansbands
        # Set up metadata and output from IRFs computation
        metadata = Dict{Symbol,Any}()
        metadata[:para] = input_type
        metadata[:cond_type] = :none
        metadata[:product] = :dsgevarlambdairf
        metadata[:class] = :obs # We default everything to an observable
        metadata[:date_inds] = OrderedDict()

        # Set up for loop over variable names
        means = DataFrame()
        bands = Dict{Symbol,DataFrame}()
        metadata[:indices] = get_observables(m)

        # Means and Bands for each variable in a class
        if method == :rotation && !draw_shocks
            for (shock, shock_i) in get_shocks(m)
                for (name, name_i) in get_observables(m)
                    # irf_output is Vector{nperiod x nobs x nshocks} -> for each observable,
                    # we want to select its specific IRF, i.e. map(x -> x[:,obs_index, shock_index]).
                    # This creates a nperiod x ndraws matrix, which we want to transpose
                    # to get a ndraws x nperiod matrix
                    single_var = Matrix(reduce(hcat, map(x -> x[:, name_i, shock_i], irf_output))')
                    means[!, Symbol(name, :__, shock)] = vec(mean(single_var, dims = 1))
                    bands[Symbol(name, :__, shock)]    = find_density_bands(single_var, density_bands;
                                                                            minimize = minimize)
                end
            end
        else
            for (name, name_i) in get_observables(m)
                # irf_output is Vector{nperiod x nobs} -> for each observable,
                # we want to select its specific IRF, i.e. map(x -> x[:,obs_index]).
                # This creates a nperiod x ndraws matrix, which we want to transpose
                # to get a ndraws x nperiod matrix
                single_var = Matrix(reduce(hcat, map(x -> x[:, name_i], irf_output))')
                means[!, name] = vec(mean(single_var, dims = 1))
                bands[name]    = find_density_bands(single_var, density_bands;
                                                    minimize = minimize)
            end
        end
        mb = MeansBands(metadata, means, bands)

        # Save MeansBands
        if !test_meansbands
            tail = if method == :cholesky
                :cholesky
            elseif method == :maxBC || method == :maximum_business_cycle_variance
                :maxBC
            elseif method == :rotation
                :rotation
            else
                :choleskyLR
            end

            fp = get_meansbands_output_file(m, input_type, :none,
                                            Symbol(metadata[:product], :obs_, tail),
                                            forecast_string = forecast_string)
            dirpath = dirname(fp)
            isdir(dirpath) || mkpath(dirpath)
            JLD2.jldopen(fp, true, true, true, IOStream) do file
                write(file, "mb", mb)
            end
            println(verbose, :high, "  " * "wrote " * basename(fp))
        end
        return mb
    else
        if method in [:rotation]
            # Reshape irf_output to nobs x nperiod x nshock x ndraw
            return cat(irf_output..., dims = ndims(irf_output[1]) + 1)
        else
            # Reshape irf_output to nobs x nperiod x ndraw
            return cat(irf_output..., dims = 3)
        end
    end
end

function impulse_responses(m::AbstractDSGEVARModel{S}, paras::Vector{S},
                           data::Matrix{S}, input_type::Symbol, method::Symbol;
                           parallel::Bool = false,
                           frequency_band::Tuple{S,S} = (2*π/32, 2*π/6),
                           n_obs_shock::Int = 1, draw_shocks::Bool = false,
                           flip_shocks::Bool = false,
                           X̂::AbstractMatrix{S} = Matrix{S}(undef, 0, 0),
                           deviations::Bool = false, normalize_rotation::Bool = false,
                           density_bands::Vector{Float64} = [.5, .6, .7, .8, .9],
                           create_meansbands::Bool = false,
                           minimize::Bool = true,
                           forecast_string::String = "",
                           verbose::Symbol = :high) where {S<:Real}
    return impulse_responses(m, reshape(paras, 1, length(paras)),
                             data, input_type, method; parallel = parallel,
                             frequency_band = frequency_band, n_obs_shock = n_obs_shock,
                             draw_shocks = draw_shocks, flip_shocks = flip_shocks,
                             X̂ = X̂, deviations = deviations, normalize_rotation = normalize_rotation,
                             density_bands = density_bands,
                             create_meansbands = create_meansbands, minimize = minimize,
                             forecast_string = forecast_string, verbose = verbose)
end


text:
using CompScienceMeshes
using BEAST
using LinearAlgebra
using Test
using StaticArrays

o, x, y, z = euclidianbasis(3)

p1 = 2x
p2 = y
p3 = 3z
p4 = o
tet = simplex(p1, p2, p3, p4)

q = abs(CompScienceMeshes.relorientation([1,2,3],[1,2,3,4]))
tri = simplex(p1, p2, p3)
n = normal(tri)

i = 3
j = 3
edg = simplex(p1, p2)

T = Float64
x = BEAST.NDLCDRefSpace{T}()
y = BEAST.NDRefSpace{T}()

p = neighborhood(tet, [0.5, 0.5, 0.0])
r = neighborhood(tri, [0.5, 0.5])

@test carttobary(edg, cartesian(p)) ≈ [0.5]
@test carttobary(edg, cartesian(r)) ≈ [0.5]

xp = x(p)[j].value
yr = y(r)[i].value

a, b = extrema((n × xp) ./ yr)
@test a ≈ b

tgt = p2 - p1
a2 = dot(n × xp, tgt) / dot(yr, tgt)
@test a ≈ a2

z = a * yr
@test z ≈ n × xp

volume(simplex(p1,p2,p4))
volume(simplex(p1,p2))

Q = BEAST.ttrace(x, tet, q, tri)

p = neighborhood(tet, [1/3, 1/3, 1/3])
r = neighborhood(tri, [1/3, 1/3])
@test cartesian(p) ≈ cartesian(r)

xp = n × x(p)[j].value
yr = y(r)[i].value
@test xp ≈ yr * Q[i,j]

x = BEAST.NDLCCRefSpace{Float64}()
y = BEAST.RTRefSpace{Float64}()
for q in 1:4
    q = 3
    fc = BEAST.faces(tet)[q]
    Q = BEAST.ttrace(x, tet, q, fc)

    nbdi = CompScienceMeshes.center(fc)
    nbdj = neighborhood(tet, carttobary(tet, cartesian(nbdi)))

    xvals = x(nbdj)
    yvals = y(nbdi)

    for j in 1:6
        trc = sum(Q[i,j]*yvals[i].value for i in 1:3)
        @test isapprox(trc, normal(fc) × xvals[j].value, atol=1e-4)
    end
end

# test the case where intrinsic and extrinsic orientations differ
o, x, y, z = euclidianbasis(3)
fc = simplex(z,o,y)
x = BEAST.NDLCCRefSpace{Float64}()
y = BEAST.RTRefSpace{Float64}()
Q = BEAST.ttrace(x, tet, 3000, fc)

nbdi = CompScienceMeshes.center(fc)
nbdj = neighborhood(tet, carttobary(tet, cartesian(nbdi)))

@test cartesian(nbdi) ≈ cartesian(nbdj)

xvals = x(nbdj)
yvals = y(nbdi)

for j in 1:6
    trc = sum(Q[i,j]*yvals[i].value for i in 1:3)
    @test isapprox(trc, -normal(fc) × xvals[j].value, atol=1e-4)
end

o, x, y, z = euclidianbasis(3)
m = Mesh([x,y,z,o], [@SVector[1,2,3,4]])

m1 = skeleton(m,1)
X = BEAST.nedelecc3d(m, m1)
@test numfunctions(X) == 6

# m2 = skeleton(m,2)
m2 = boundary(m)
Y = BEAST.ttrace(X,m2)
@test numfunctions(Y) == 6

pa = Y.fns[1][1].cellid
pb = Y.fns[1][2].cellid

tria = chart(m2, cells(m2)[pa])
trib = chart(m2, cells(m2)[pb])

ra = Y.fns[1][1].refid
rb = Y.fns[1][2].refid

ctra = cartesian(CompScienceMeshes.center(CompScienceMeshes.edges(tria)[ra]))
ctrb = cartesian(CompScienceMeshes.center(CompScienceMeshes.edges(trib)[rb]))
@test ctra ≈ ctrb
# tri1 = chart(m, cells(m)[Y.fns[1][1].cellid])

# ctr1 = cartesian(center(CompScienceMeshes.edges(tri)[]))

ext: jl


text:
# copied from Pluto, should split this out
# at some point?
@option mutable struct CompilerOptions
    compile::Union{Nothing,String} = nothing
    sysimage::Union{Nothing,String} = nothing
    banner::Union{Nothing,String} = nothing
    optimize::Union{Nothing,Int} = nothing
    math_mode::Union{Nothing,String} = nothing

    # notebook specified options
    # the followings are different from
    # the default julia compiler options

    # we use nothing to represent "@v#.#"
    # unlike Pluto, we use parsed content of the Project.toml here
    # so that we can transfer the TOML file through sockets.
    project::Dict{String, Any} = TOML.parsefile(Base.load_path_expand("@."))
    # we don't load startup file in notebook
    startup_file::Union{Nothing,String} = "no"
    # we don't load history file in notebook
    history_file::Union{Nothing,String} = "no"

    threads::Union{Nothing,String,Int} = get_nthreads()
end

function _convert_to_exeflags(options::CompilerOptions)::Vector{String}
    option_list = String[]

    workspace_deps = Dict{String, Any}(
        "Suppressor" => "fd094767-a336-5f1f-9728-57cf17d0bbfb"
    )

    for name in fieldnames(CompilerOptions)
        name === :project && continue # skip project
        flagname = string("--", replace(String(name), "_" => "-"))
        value = getfield(options, name)
        if value !== nothing
            if !(VERSION <= v"1.5.0-" && name === :threads)
                push!(option_list, string(flagname, "=", value))
            end
        end
    end

    # handle in memory project
    # only deps and compat matters for workspace
    temp_project_dir = mktempdir(;prefix="workspace_")
    temp_project = joinpath(temp_project_dir, "Project.toml")
    open(temp_project, "w+") do io
        deps = get(options.project, "deps", Dict{String, Any}())
        compat = get(options.project, "compat", Dict{String, Any}())

        merge!(deps, workspace_deps)
        TOML.print(io, Dict("deps"=>deps, "compat"=>compat);
        sorted=true, by=key -> (Pkg.Types.project_key_order(key), key))
    end

    push!(option_list, "--project=$temp_project")
    return option_list
end

@option struct DistributedOptions
    dir::String = pwd()
    enable_threaded_blas::Bool = false
    exename::String = joinpath(Sys.BINDIR, "julia")
    topology::Symbol = :all_to_all
    lazy::Bool = true
end

function addprocs_flags(options::DistributedOptions)
    kw_list = []
    for name in fieldnames(DistributedOptions)
        value = getfield(options, name)
        isnothing(value) && continue
        push!(kw_list, name => value)
    end
    return kw_list
end

@option struct WorkspaceInstance
    name::Symbol = gensym(:workspace)
    uuid::UUID = uuid1()
    compiler::CompilerOptions = CompilerOptions()
    distributed::DistributedOptions = DistributedOptions()
    # scripts that is evaluated in this workspace instance
    linked_scripts::Vector{String} = String[]
end

StructTypes.StructType(::Type{WorkspaceInstance}) = StructTypes.Struct()

function Configurations.convert_to_option(::Type{WorkspaceInstance}, ::Type{Symbol}, x::String)
    return Symbol(x)
end

function Configurations.convert_to_option(::Type{WorkspaceInstance}, ::Type{UUID}, x::String)
    return UUID(x)
end

struct WorkspaceTable
    workspaces::Dict{UUID, WorkspaceInstance}
    process::Dict{UUID, Int}
end

WorkspaceTable() = WorkspaceTable(Dict{UUID, WorkspaceInstance}(), Dict{UUID, Int}())

function get_nthreads()
    haskey(ENV, "JULIA_NUM_THREADS") || return
    isempty(ENV["JULIA_NUM_THREADS"]) && return
    return parse(Int, ENV["JULIA_NUM_THREADS"])
end

function add_workspace_process!(table::WorkspaceTable, ins::WorkspaceInstance)
    pid = addprocs(1;addprocs_flags(ins.distributed)..., exeflags=_convert_to_exeflags(ins.compiler))[]

    # init workspace project env
    Distributed.remotecall_eval(Main, [pid], quote
        using Pkg; Pkg.instantiate()
    end)
    # create the workspace module
    Distributed.remotecall_eval(Main, [pid], Expr(:toplevel, :(module $(ins.name) end)))
    Distributed.remotecall_eval(Main, [pid], Expr(:toplevel, :(using Suppressor)))
    table.workspaces[ins.uuid] = ins
    table.process[ins.uuid] = pid
    return table
end

rm_workspace!(table::WorkspaceTable, ins::WorkspaceInstance) = rm_workspace!(table, ins.uuid)

function rm_workspace!(table::WorkspaceTable, uuid::UUID)
    haskey(table.workspaces, uuid) || error("cannot find workspace: $uuid")
    wait(rmprocs(table.process[uuid]; waitfor=5))
    delete!(table.workspaces, uuid)
    delete!(table.process, uuid)
    return table
end

function eval_in_workspace(table::WorkspaceTable, ins::WorkspaceInstance, expr)
    haskey(table.process, ins.uuid) || error("workspace process is not spawned, call add_workspace_process")
    pid = table.process[ins.uuid]
    Distributed.remotecall_fetch(Core.eval, pid, Main, quote
        Suppressor.@capture_out begin
            Core.eval($(ins.name), $(QuoteNode(expr)))
        end
    end)
end


text:
# add more julia processes
using Distributed
nprocs() < 3 ? addprocs(4-nprocs()) : nothing

# tell Julia that you want to use these modules/namespaces
using RoME, Distributions

# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*Matrix{Float64}(LinearAlgebra.I, 3,3))))

# Drive around in a hexagon a number of times
for i in 0:12
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal( [0.1;0.1;0.1].^2 ) )))
    addFactor!(fg, [psym;nsym], pp )
end

# perform inference, and remember first runs are slower owing to Julia's just-in-time compiling
batchSolve!(fg)

## Inter-operating visualization packages for Caesar/RoME/IncrementalInference exist
using RoMEPlotting
using Compose

# For Juno/Jupyter style use
pl = drawPoses(fg)
# For scripting use-cases you can export the image
Gadfly.push_theme(:default) # light background, where Juno uses dark background
Gadfly.draw(Gadfly.PDF("/tmp/test1.pdf", 20cm, 10cm),pl)  # or PNG(...)

# Add a landmark l1
addVariable!(fg, :l1, Point2, labels=["LANDMARK"])

bear1 = atan(10,20)
bear2 = atan(10,10) - pi/3

# Add landmarks with Bearing range measurements at x0, x6, x12, x18, x24, x30...
vars = ls(fg)[1] # Get variables
for xIndex in 1:6:length(vars)
    info("Creating factor between $(vars[xIndex]) and l1...")
    p2br = Pose2Point2Bearing(Normal(bear1,0.05))
    addFactor!(fg, [@show vars[xIndex]; :l1], p2br)
    p2br = Pose2Point2Bearing(Normal(bear2,0.05))
    addFactor!(fg, [@show vars[xIndex+1]; :l1], p2br)
end

# plotKDE(fg, [:x0;:x6;:x12], dims=[1;2])
# plotKDE(fg, [:x0;:x6;:x12], dims=[3])
# plotKDE(fg, [:x1;:x7;:x13], dims=[1;2])

# writeGraphPdf(fg)

# Initialize :l1 numerical values but do not rerun solver
ensureAllInitialized!(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(Gadfly.PDF("/tmp/test2.pdf", 20cm, 10cm),pl)  # or PNG(...)

# solve
batchSolve!(fg)

# redraw to see that additional information aids in refining the results
pl = drawPosesLandms(fg)
Gadfly.draw(Gadfly.PDF("/tmp/test3.pdf", 20cm, 10cm),pl)  # or PNG(...)


text:
"""
    RSparseBitSet{T}(size, trailer)

Create a Reversible Sparse BitSet with `size` elements.

This datastructure follows the implementation from Demeulenaere J. et al.
(2016) Compact-Table: Efficiently Filtering Table Constraints with Reversible Sparse Bit-Sets.
https://doi.org/10.1007/978-3-319-44953-1_14 .
"""
struct RSparseBitSet{T}
    words::Vector{StateObject{T}}
    index::Vector{Int}
    limit::StateObject{Int}
    mask::Vector{T}

    function RSparseBitSet{T}(size::Int, trailer) where T <: Unsigned
        n = sizeof(T) * 8
        p = Int(ceil(size/n))
        rest = size % n
        words = Vector{StateObject{T}}(undef, p)
        for i = 1:(p-1)
            words[i] = StateObject{T}(typemax(T), trailer)
        end
        if rest == 0
            words[p] = StateObject{T}(typemax(T), trailer)
        else
            words[p] = StateObject{T}(parse(T, "0b" * "1"^rest * "0"^(n-rest)), trailer)
        end

        index = collect(Int, 1:p)
        limit = StateObject{Int}(p, trailer)
        mask = Vector{T}(undef, p)

        return new(words, index, limit, mask)
    end
end

RSparseBitSet(size::Int, trailer) = RSparseBitSet{UInt64}(size::Int, trailer)

"""
    clearMask!(::RSparseBitSet{T})

Set all the masks to 0.
"""
function clearMask!(set::RSparseBitSet{T})::Nothing where T <: Unsigned
    offsets = set.index[1:set.limit.value]
    set.mask[offsets] .= T(0)
    return
end

"""
    reverseMask!(::RSparseBitSet)

Set the mask to its opposite: Boolean NOT.
"""
function reverseMask!(set::RSparseBitSet{T})::Nothing where T <: Unsigned
    offsets = set.index[1:set.limit.value]
    set.mask[offsets] .= .~ set.mask[offsets]
    return
end

"""
    addToMask!(::RSparseBitSet{T}, m)

Add the content of m to the mask: Boolean OR.
"""
function addToMask!(set::RSparseBitSet{T}, m::Vector{T})::Nothing where T <: Unsigned
    offsets = set.index[1:set.limit.value]
    set.mask[offsets] .= set.mask[offsets] .| m[offsets]
    return
end

"""
    intersectWithMask!(::RSparseBitSet{T})

Apply the mask to words: Bolean AND.
"""
function intersectWithMask!(set::RSparseBitSet{T})::Nothing where T <: Unsigned
    for i = set.limit.value:-1:1
        offset = set.index[i]
        word = set.words[offset].value & set.mask[offset]
        setValue!(set.words[offset], word)
        if word == 0
            set.index[i], set.index[set.limit.value] = set.index[set.limit.value], set.index[i]
            setValue!(set.limit, set.limit.value-1)
        end
    end
    return
end

"""
    intersectIndex(::RSparseBitSet{T}, m)

Return the index of the first word intersecting m, or -1 if none exist.
"""
function intersectIndex(set::RSparseBitSet{T}, m::Vector{T})::Int where T <: Unsigned
    for i = 1:set.limit.value
        offset = set.index[i]
        if set.words[offset].value & m[offset] != 0
            return offset
        end
    end
    return -1
end

"""
    bitVectorToUInt64Vector(bitset)::Vector{UInt64}

Convert a Julia BitVector to a vector preformatted for the RSparseBitSet{UInt64}.

# Arguments
- `bitset::BitVector`: the BitVector to convert.
"""
function bitVectorToUInt64Vector(bitset::BitVector)::Vector{UInt64}
    return [bitreverse(chunk) for chunk in bitset.chunks]
end

function Base.BitVector(set::RSparseBitSet{UInt64})
    bitvector = BitVector(undef, 8*sizeof(UInt64)*length(set.words))
    bitvector .= false
    for i = set.limit.value:-1:1
        offset = set.index[i]
        bitvector.chunks[offset] = bitreverse(set.words[offset].value)
    end
    return bitvector
end

function Base.isempty(set::RSparseBitSet{T})::Bool where T <: Unsigned
    return set.limit.value == 0
end

function Base.getindex(set::RSparseBitSet{T}, idx::Integer) where T <: Unsigned
    n = sizeof(T) * 8
    wordIndex = Int(ceil((idx)/n))
    offset = (n - idx%n)%n
    return (set.words[wordIndex].value & (T(1) << offset)) > 0
end

function Base.show(io::IO, ::MIME"text/plain", set::RSparseBitSet{T}) where T <: Unsigned
    n = sizeof(T) * 8
    println(io, string(typeof(set)), ": index = ", set.index, ", limit = ", set.limit.value)
    println(io, "   words: ", join([string(x.value, base=16, pad=Int(n/4)) for x in set.words], " "))
    println(io, "   mask:  ", join([string(x, base=16, pad=Int(n/4)) for x in set.mask], " "))
end

function Base.show(io::IO, set::RSparseBitSet{T}) where T <: Unsigned
    n = sizeof(T) * 8
    println(io, string(typeof(set)), ": index = ", set.index, ", limit = ", set.limit.value)
    println(io, "   words: ", join([string(x.value, base=16, pad=Int(n/4)) for x in set.words], " "))
end


text:
using Hip
using Base.Test

using Compat, BenchmarkTools
const hippath = Pkg.dir("Hip") * "/src";
const lib_hip_jl = Libdl.dlopen(hippath * "/hip_jl.so")

# write your own tests here
const N = 4096
a = randn(Float32, N, N);
b = randn(Float32, N, N);
c = zeros(a);

ah = HipArray(a);
bh = HipArray(b);
#ch = HipArray(c);
cht = HipArray(c);

tilesgemm!(cht, 1f0, ah, bh, 0f0)
A_mul_B!(c, a, b)
synchronize()
cres = Array(cht)

chh = HipArray(c);
sgemm!(chh, 1f0, ah, bh, 0f0)
synchronize()
cres2 = Array(chh)

free!(ah)
free!(bh)
free!(cht)
N2 = 32
a = randn(Float32, N2, N2);
b = randn(Float32, N2, N2);
c = zeros(a);
ah = HipArray(a);
bh = HipArray(b);
ch = HipArray(c);
tilesgemm!(ch, 1f0, ah, bh, 0f0)
A_mul_B!(c, a, b)
synchronize(); Array(ch)

A_mul_Bt!(c, a, b)
At_mul_B!(c, a, b)
At_mul_Bt!(c, a, b)
Array(ch)

function set_ab!(a, b, ia = 0.25f0, ja = 1f0, ib = 0.5f0, jb = 2f0)
    for i in 1:32
        for j in 1:32
            a[j,i] = ia*i + ja*j
            b[j,i] = ib*i + jb*j
        end
    end
end
free!(ah)
free!(bh)
ah = HipArray(a);
bh = HipArray(b);
tilesgemm!(ch, 1f0, ah, bh, 0f0)
A_mul_B!(c, a, b)
copy!(cha, ch); cha

chn = HipArray(c);
naivesgemm!(chn, 1f0, ah, bh, 0f0)

chh = HipArray(c);
sgemm!(chh, 1f0, ah, bh, 0f0)

function sync_func(f, c, alpha, a, b, beta)
    f(c, alpha, a, b, beta)
    synchronize()
end

copy!(cres, chn); cres
copy!(cres, chh); cres

@benchmark A_mul_B!($c, $a, $b)
@benchmark sync_func($naivesgemm!, $cht, 1f0, $ah, $bh, 0f0)
@benchmark sync_func($tilesgemm!, $cht, 1f0, $ah, $bh, 0f0)
@benchmark sync_func($sgemm!, $chh, 1f0, $ah, $bh, 0f0)

const THREADS_PER_BLOCK_X = 32
const THREADS_PER_BLOCK_Y = 32
const hipBlockDim_x = 32
const hipBlockDim_y = 32
const hipBlockIdx_x = 0
const hipBlockIdx_y = 0

function run_kern!(c, alpha, a, b, beta)

    M, K = size(a)
    N = size(b,2)

    c .*= beta

    num_tiles = div(K,THREADS_PER_BLOCK_X);

    tileA = Vector{Float32}(THREADS_PER_BLOCK_X * THREADS_PER_BLOCK_Y)
    tileB = Vector{Float32}(THREADS_PER_BLOCK_X * THREADS_PER_BLOCK_Y)
    for t in 0:num_tiles-1

        for i in 0:THREADS_PER_BLOCK_X-1, j in 0:THREADS_PER_BLOCK_Y-1
            kern_shared!(tileA, tileB, i, j, t, M, K, a, b)
        end
        for i in 0:THREADS_PER_BLOCK_X-1, j in 0:THREADS_PER_BLOCK_Y-1
            kern_tile!(c, a, b, i, j, alpha, tileA, tileB, M)
        end

    end
    c, tileA, tileB
end
function kern_shared!(tileA, tileB, hipThreadIdx_x, hipThreadIdx_y, t, M, K, A, B)
    idx_x = hipBlockIdx_x * hipBlockDim_x + hipThreadIdx_x;
    idx_y = hipBlockIdx_y * hipBlockDim_y + hipThreadIdx_y;

    id_local = hipThreadIdx_y + THREADS_PER_BLOCK_Y * hipThreadIdx_x;
    id_localT = hipThreadIdx_x + THREADS_PER_BLOCK_X * hipThreadIdx_y;

    t_col = t*THREADS_PER_BLOCK_X;
    tileA[id_localT+1] = A[(t_col + hipThreadIdx_x)*M + idx_y+1];
    tileB[id_local+1] = B[t_col + hipThreadIdx_y + K*idx_x+1];
end

function kern_tile!(C, A, B, hipThreadIdx_x, hipThreadIdx_y, alpha, tileA, tileB, M)

    idx_x = hipBlockIdx_x * hipBlockDim_x + hipThreadIdx_x;
    idx_y = hipBlockIdx_y * hipBlockDim_y + hipThreadIdx_y;

    myIdx = idx_y + idx_x * M;

    local_c = zero(eltype(C))

    idxA = hipThreadIdx_y * THREADS_PER_BLOCK_X ;
    idxB = hipThreadIdx_x * THREADS_PER_BLOCK_X ;

    for k in 0:THREADS_PER_BLOCK_X-1
        local_c += tileA[k + idxA+1] * tileB[k + idxB+1];
    end

    C[myIdx+1] += alpha * local_c;
end


text:
### higher level problem-solve functions

function _solve{T<:AbstractFloat}(
    pb::Problem{T},
    reg::Regularizer,
    θ::Array{T},
    solver::DescentSolver,
    options::Options,
    callback::Function)

    if has_bias(pb)
        rmodel = riskmodel(pred_with_bias(pb), loss(pb))
        f = RegRiskFun(rmodel, reg, inputs(pb), outputs(pb))
        solve!(solver, f, θ, options, callback)

    else
        rmodel = riskmodel(pred_without_bias(pb), loss(pb))
        f = RegRiskFun(rmodel, reg, inputs(pb), outputs(pb))
        solve!(solver, f, θ, options, callback)

    end::Solution{typeof(θ)}
end

function solve{T<:AbstractFloat}(
    pb::Problem{T};
    reg::Regularizer=ZeroReg(),
    init::StridedArray{T}=T[],
    solver::Solver=default_solver,
    options::Options=Options(),
    callback::Function=no_op)

    θ = isempty(init) ? initsol(pb) : copy(init)
    _solve(pb, reg, θ, solver, options, callback)
end

### Specific solvers

prep_dir(::Solver, g::Array) = similar(g)
has_accelerate(::Solver) = false

function solve!{T<:AbstractFloat}(solver::DescentSolver,
    f::Functional{T},       # the objective function
    θ::Array{T},            # the solution (which would be updated inplace)
    options::Options,       # options to control the procedure
    callback::Function)     # callback function

    ## extract arguments and options

    maxiter = options.maxiter
    vbose = verbosity_level(options.verbosity)::Int

    ## prepare storage

    θ2 = similar(θ)    # tempoarily new parameter (linear search)
    g = similar(θ)     # gradient
    g2 = similar(θ)    # another gradient
    p = prep_dir(solver, g)

    buffer = ndims(θ) == 2 ? zeros(size(θ,1), size(f.X,2)) : zeros(size(f.X,2))

    ## main loop
    t = 0
    converged = false
    v, _ = value_and_grad!(buffer, f, g, θ)

    if vbose >= VERBOSE_ITER
        print_iter_head()
        print_iter(t, v)
    end

    states = init_states(solver, θ, g)

    while !converged && t < maxiter
        t += 1
        v_pre = v

        # accelerate solution
        if has_accelerate(solver)
            accelerate!(solver, t, states, θ)
            v, _ = value_and_grad!(buffer, f, g, θ)
        end

        # descent direction: g
        descent_dir!(solver, t, states, p, θ, θ2, g, g2)

        # backtracking
        α = backtrack!(buffer, f, θ2, θ, v, g, p, one(T), options)
        θ, θ2 = θ2, θ  # swap current solution and previous solution

        # compute new gradient
        v, _ = value_and_grad!(buffer, f, g2, θ)
        g, g2 = g2, g  # swap current gradient with previous gradient

        # test convergence
        converged = test_convergence(θ, θ2, v, v_pre, g, options)

        # print iteration
        if vbose >= VERBOSE_ITER
            print_iter(t, v, v_pre, g, α)
        end

        # invoke callback (when requested)
        callback(t, θ, v, g)
    end

    if vbose >= VERBOSE_FINAL
        print_final(t, v, converged)
    end

    return Solution(θ, v, t, converged)
end

## Gradient Descent Solver

type GD <: DescentSolver
end

Base.show(io::IO, ::GD) = print(io, "GD")
init_states(::GD, θ, g) = nothing

type AGD <: DescentSolver   # Accelerated GD
end

type AGDStates{T}
    d::Int
    τ::T
    xpre::Vector{T}
    y::Vector{T}
end

Base.show(io::IO, ::AGD) = print(io, "AGD")
function init_states{T<:Real}(::AGD, θ::Array{T}, g::Array{T})
    d = length(θ)
    length(g) == d || throw(DimensionMismatch())
    AGDStates{T}(d, one(T), Array(T, d), Array(T, d))
end

prep_dir(::Union{GD,AGD}, g::Array) = g

descent_dir!{T<:Real}(::Union{GD,AGD}, t::Int, states, p::Array{T},
                      θ::Array{T}, θp::Array{T}, g::Array{T}, gp::Array{T}) =
    (is(p, g) || copy!(p, g); p)

has_accelerate(::AGD) = true
function accelerate!{T}(::AGD, t::Int, states::AGDStates{T}, x::Array{T})
    τ = states.τ
    xpre = states.xpre
    y = states.y

    # compute y
    if t > 1
        d = states.d
        τ2 = (one(T) + sqrt(one(T) + 4 * abs2(τ))) * convert(T, 0.5)
        c = (τ - one(T)) / τ2
        @inbounds for i = 1:d
            x_i = x[i]
            y[i] = x_i + c * (x_i - xpre[i])
        end
        # @show τ2
        states.τ = τ2
    end

    # store current x as xpre
    copy!(xpre, x)

    # x <- y
    if t > 1
        copy!(x, y)
    end
    x
end

## BFGS Solver

type BFGS <: DescentSolver
end

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

immutable BFGSStates{T}
    d::Int
    Λ::Matrix{T}
    s::Vector{T}
    y::Vector{T}
    z::Vector{T}
end

function init_states{T<:Real}(::BFGS, θ::Array{T}, g::Array{T})
    d = length(θ)
    length(g) == d || throw(DimensionMismatch())
    BFGSStates{T}(d, eye(T, d), Array(T, d), Array(T, d), Array(T, d))
end

function descent_dir!{T<:Real}(::BFGS, t::Int, states::BFGSStates{T}, p::Array{T},
                               θ::Array{T}, θp::Array{T}, g::Array{T}, gp::Array{T})
    if t > 1
        d = states.d
        Λ = states.Λ
        s = states.s
        y = states.y
        z = states.z

        # compute s <- θ - θ_pre
        @inbounds for i = 1:d
            s[i] = θ[i] - θp[i]
        end

        # compute y <- g - g_pre
        @inbounds for i = 1:d
            y[i] = g[i] - gp[i]
        end
        A_mul_B!(z, Λ, y)

        u = dot(s, y)
        c1 = (u + dot(y, z)) / (u * u)
        c2 = one(T) / u

        # update Λ
        @inbounds for j = 1:d
            s_j, z_j = s[j], z[j]
            for i = 1:d
                s_i, z_i = s[i], z[i]
                Λ[i,j] += (c1 * s_i * s_j - c2 *(z_i * s_j + z_j * s_i))
            end
        end

        # compute direction
        A_mul_B!(vec(p), Λ, vec(g))
    else
        copy!(p, g)
    end
    p
end

## L-BFGS Solver

immutable LBFGS <: DescentSolver
    m::Int    # the numbe of memorized steps
end

Base.show(io::IO, s::LBFGS) = print(io, "LBFGS($(s.m))")

immutable LBFGSStates{T}
    d::Int
    Δx::Matrix{T}
    Δg::Matrix{T}
    ρ::Vector{T}
    α::Vector{T}
    q::Vector{T}
end

function init_states{T<:Real}(solver::LBFGS, θ::Array{T}, g::Array{T})
    d = length(θ)
    length(g) == d || throw(DimensionMismatch())
    m = solver.m
    LBFGSStates{T}(d,
        zeros(T, d, m),     # Δx
        zeros(T, d, m),     # Δg
        zeros(T, m),        # ρ
        zeros(T, m),        # α
        zeros(T, d))        # q
end

function descent_dir!{T<:Real}(solver::LBFGS, t::Int, states::LBFGSStates{T}, p::Array{T},
                               θ::Array{T}, θp::Array{T}, g::Array{T}, gp::Array{T})

    if t > 1
        # extract fields
        m = solver.m
        d = states.d
        Δx = states.Δx
        Δg = states.Δg
        ρ = states.ρ
        α = states.α
        q = states.q
        p_ = vec(p)

        # q <- g
        copy!(q, vec(g))

        # backward pass
        h = min(t-1, m)
        for idx = t-1:-1:t-h
            i = mod1(idx, m)
            α[i] = ρ[i] * dot(view(Δx,:,i), q)
            axpy!(-α[i], view(Δg,:,i), q)
        end

        # p <- q
        copy!(p_, q)

        # forward pass
        for idx = t-h:t-1
            i = mod1(idx, m)
            β = ρ[i] * dot(view(Δg,:,i), p_)
            axpy!(α[i] - β, p_, view(Δx, :, i))
        end
    else
        copy!(p, g)
    end
    p
end


text:
using GeoStatsBase
using DelimitedFiles
using Distances
using Distributions
using DataFrames
using Plots, VisualRegressionTests
using Test, Pkg, Random

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
datadir = joinpath(@__DIR__,"data")
visualtests = !istravis || (istravis && islinux)
if !istravis
  Pkg.add("Gtk")
  using Gtk
end

# dummy variables for testing
include("dummy.jl")

# list of tests
testfiles = [
  "distances.jl",
  "distributions.jl",
  "data.jl",
  "domains.jl",
  "collections.jl",
  "partitioning.jl",
  "paths.jl",
  "neighborhoods.jl",
  "mappers.jl",
  "problems.jl",
  "solvers.jl",
  "solutions.jl",
  "utils.jl",
  "statistics.jl",
  "learning.jl",
  "plotrecipes.jl"
]

@testset "GeoStatsBase.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
# compare the time taken to do simple subtractions compared to <two different implementations of floor
# with a view to making some routiens that use subtraction a lot faster
using IndexedTables
using SurrealNumbers

# old code for reference, e.g., in case you want to rerun some timing tests to understand
# why current code is written the way it is

function floor2(s::SurrealFinite)
    if zero(s) <= s < one(s)
        return zero(s)
    elseif s < zero(s)
        return floor(s + one(s)) - one(s)
    elseif s >= one(s)
        return floor(s - one(s)) + one(s)
    end
end

function isinteger2(s::SurrealFinite)
    if s ≅ zero(s)
        return true
    elseif s <= -one(s)
        return isinteger(s + one(s))
    elseif s >= one(s)
        return isinteger(s - one(s))
    else
        return false
    end
end

# not the more general form of rounding defined in Julia -- should fix
function round2(s::SurrealFinite)
    if s ≅ zero(s)
        return s
    elseif s <= -one(s)
        return round(s + one(s)) - one(s)
    elseif s >= one(s)
        return round(s - one(s)) + one(s)
    elseif s >= convert(SurrealFinite, 1//2)
        return one(s)
    elseif s >= convert(SurrealFinite, -1//2)
        return zero(s)
    else
        return -one(s)
    end
end

convert(SurrealFinite, 2) * convert(SurrealFinite, 1)

#ns = 6:-1:0
ns = 0:4
n = length(ns)
m = 100 * ones(size(ns))
x1 = Array{SurrealFinite,1}(n)
x2 = Array{SurrealFinite,1}(n)
t1 = Array{Float64,1}(n)
t2 = Array{Float64,1}(n)
bytes2 = Array{Float64,1}(n)
bytes1 = Array{Float64,1}(n)
gctime1 = Array{Float64,1}(n)
gctime2 = Array{Float64,1}(n)
s = convert(SurrealFinite, 1)
floor(s)
floor2(s)
for i=1:n
    println("  working on ", ns[i])
    t1[i] = 0
    t2[i] = 0
    b1 = 0
    g1 = 0
    mem1 = 0
    val1 = 0
    b2 = 0
    g2 = 0
    mem2 = 0
    val2 = 0
    for j=1:m[i]
        s = convert(SurrealFinite, ns[i] + 0.5)
        # val, time, b, g, mem = @timed convert(Rational, s)
        val1, time1, b1, g1, mem1 = @timed floor( s )
        val2, time2, b2, g2, mem2 = @timed floor2( s )
        t1[i] += time1
        t2[i] += time2
    end
    x1[i] = val1
    bytes1[i] = b1
    gctime1[i] = g1
    x2[i] = val2
    bytes2[i] = b2
    gctime2[i] = g2
end
# s = size.(x)
# d = depth_av.(x)
# f = n_zeros.(x)
t1 = t1./m[1:n]
bytes1 = bytes1./m[1:n]
t2 = t2./m[1:n]
bytes2 = bytes2./m[1:n]

t1_p_ms = 1.0e6 * t1
t2_p_ms = 1.0e6 * t2

tab1 = table(@NT( i=0:n-1,
                time1=t1_p_ms,
                time2=t2_p_ms,
                bytes1=bytes1,
                bytes2=bytes2);
           pkey = [:i])


text:
# load functions and packages
try
  cd("DWPSDE model")
catch
  warn("Already in the DWPSDE model folder")
end

# load case models
cd("..")
include(pwd()*"\\select case\\selectcase.jl")
cd("DWPSDE model")

include("set_up.jl")
#include("run_pf_paralell.jl")
using PyPlot

# set parameters
nbr_iterations = 1000 # nbr of iteration of the MC algorithm
nbr_particels = 25  # nbr of particel on each core
burn_in = 500 # nbr of iterations for the burn in

# create a problem type where we estimate twp parameters and use the AM aglorithm for
# the adaptive updating
data_set = "old"
problem = set_up_problem(nbr_of_unknown_parameters=7,use_sim_data = true,data_set=data_set)
problem.alg_param.nbr_of_cores = 4

problem.alg_param.R = nbr_iterations
problem.alg_param.N = nbr_particels
problem.alg_param.burn_in = burn_in
problem.adaptive_update =  AMUpdate(eye(2), 2.4/sqrt(2), 2, 0.5, 50) # 2.38/sqrt(2)

Z = problem.data.Z
N = 5000
R = problem.alg_param.R
burn_in = problem.alg_param.burn_in
pf_alg = problem.alg_param.pf_alg
nbr_of_cores = problem.alg_param.nbr_of_cores
nbr_x0  = problem.alg_param.nbr_x0
nbr_x = problem.alg_param.nbr_x
subsample_interval = problem.alg_param.subsample_int
dt = 0.035 #problem.alg_param.dt # 0.5 if data new
dt_U = 1. #problem.alg_param.dt # 1 if data new

print_on = true
return_weigths_and_particels = false
# model parameters
theta_true = problem.model_param.theta_true
theta_known = problem.model_param.theta_known
theta_0 = problem.model_param.theta_0

# run pf at theta_true
theta = theta_true #log([1.00556;2.66301;33.8242;19.5643;1.3644;1.87421;4.10332])
#theta = [log(0.1) log(100) log(40)]

nbr_parallel = 1
nbr_itr_pf = 100*nbr_parallel
loglik_vector = zeros(100)

N = 200

(Κ, Γ, A, A_sign, B,c,d,g,f,power1,power2,sigma) = set_parameters(theta, theta_known,length(theta))

A = A*A_sign

# set value for constant b function
b_const = sqrt(2.*sigma^2 / 2.)

# use long data
problem.data.Z = convert(Array, readtable("data_new_long.csv")[:,1])
dt = 0.5
Z = problem.data.Z

# set values needed for calculations in Float64
(subsample_interval_calc, nbr_x0_calc, nbr_x_calc,N_calc) = map(Float64, (subsample_interval, nbr_x0, nbr_x,N))
loglik = @time run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, true, false, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)

nbr_of_proc = set_nbr_cores(4, pf_alg)
loglik_vec = SharedArray{Float64}(nbr_of_proc)

loglik = pf_paralell(Z, theta,theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,true,false, nbr_of_proc,loglik_vec)

println(theta)
println(theta_known)
println(N)
println(dt)
println(dt_U)
println(nbr_x0)
println(nbr_x)
println(subsample_interval)
println(nbr_of_proc)
println(loglik_vec)

tic()
for i = 1:100
  if mod(i,10) == 0
    println(i)
    println(loglik_vec)
  end
  #loglik_vector[i] = run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, false, false, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)
  loglik_vector[i] = pf_paralell(Z, theta,theta_known, N,dt,dt_U,nbr_x0, nbr_x,subsample_interval,false,false, nbr_of_proc,loglik_vec)
end
toc()

std(loglik_vector)
mean(loglik_vector)

PyPlot.figure()
h = PyPlot.plt[:hist](loglik_vector,10)

addprocs(4)
if pf_alg == "parallel_bootstrap"
  @everywhere include("run_pf_paralell.jl")
else
  @everywhere include("run_apf_paralell.jl")
end

loglik_vec = SharedArray{Float64}(4)

tic()
@parallel for i = 1:4
  loglik_vec[i] = run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, true, false, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)
end
println(loglik_vec)
logsumexp(loglik_vec) - log(4)
toc()

tic()
for j = 1:nbr_itr_pf
  if mod(j,100) == 0
    println(j)
  end

  @sync begin

  @parallel for i = 1:4
    loglik_vec[i] = run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, false, false, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)
  end

  end

  loglik_vector[j] = logsumexp(loglik_vec) - log(4)
end
toc()

PyPlot.figure()
h = PyPlot.plt[:hist](loglik_vector,10)

std(loglik_vector)
mean(loglik_vector)

N_vec = [25 50 100 500 1000]
loglik_store = zeros(length(N_vec),nbr_itr_pf)
w_stor = []
x_stor = []

# var < 2 for N > 1000

for j = 1:length(N_vec)

  N = N_vec[j]

  println(N)

  #(Κ, Γ, A,B,c,d,g,f,power1,power2,sigma) = set_parameters(theta, theta_known,length(theta))

  # set value for constant b function
  #b_const = sqrt(2.*sigma^2 / 2.)

  # set values needed for calculations in Float64
  #(subsample_interval_calc, nbr_x0_calc, nbr_x_calc,N_calc) = map(Float64, (subsample_interval, nbr_x0, nbr_x,N))

  (loglik, w, x) = @time run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, print_on, return_weigths_and_particels, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)

  for i = 1:nbr_itr_pf
    if mod(i,100) == 0
      println(i)
    end
    loglik_vector[i] = run_pf_paralell(Z,theta,theta_known, N, N_calc, dt, dt_U, nbr_x0, nbr_x0_calc, nbr_x, nbr_x_calc, subsample_interval, subsample_interval_calc, false, false, Κ, Γ, A, B, c, d, f, g, power1, power2, b_const)
  end

  loglik_store[j,:] = loglik_vector

  push!(w_stor,w)
  push!(x_stor,x)

  #=
  for i = 1:nbr_parallel:nbr_itr_pf
    loglik_vector[i] = 1/nbr_parallel*sum(loglik_vector[i:i+nbr_parallel-1])
  end

  loglik_vector = loglik_vector[1:nbr_parallel:nbr_itr_pf]

  loglik_store[i,:] = loglik_vector
  =#

end

n = 5
w = w_stor[n]
x = x_stor[n]
loglik_vector = loglik_store[n,:]
loglik_vector_print = zeros(1, length(loglik_vector))

for i = 1:length(loglik_vector)
  loglik_vector_print[i] = loglik_vector[i]
end

export_pf_diagnostics(loglik_vector_print, w, x, Z, Z)

#=
res_problem_normal_prior_est_7_AM_gen = @time MCMC(problem_normal_prior_est_7_AM_gen)
res_problem_normal_prior_est_9_AM_gen = @time MCMC(problem_normal_prior_est_9_AM_gen)

max_stars_repo_path: DWPSDE model/pf_diagnostics.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: DWPSDE model/pf_diagnostics.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: DWPSDE model/pf_diagnostics.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Structure for Parameters
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
mutable struct Str_Param
    num_act::Int
    num_obs::Int
    count_decay::Float64            # An additonal parameter for decaying counts for novelty. It is equal to 1 for all analysis in the papaer.
    p_change::Float64               # m in the paper is equal to p_change/(1-p_change)
    ϵ::Float64
    α::Float64                      # ρ_b in the paper
    Δα::Float64                     # δρ in the paper
    β_R::Float64
    β_ratio_N2R::Float64
    γ_RL_R::Float64                 # λ_R in the paper
    γ_RL_N::Float64                 # λ_N in the paper
    λ_elig_R::Float64               # μ_R in the paper
    λ_elig_N::Float64               # μ_R in the paper
    Q_R_0::Float64
    Q_N_0::Float64
    N_prior_sweep::Int              # T_PS in the paper
    w_MF_to_MB_1::Float64           # ω_11 or ω_12 in the paper
    w_MF_to_MB_2::Float64           # ω_0 in the paper
    w_MF_to_MB_Scaler::Float64      # ω_MF in the paper
    γ_RL_R_MB::Float64
    γ_RL_N_MB::Float64
    C_thresh::Int                   # For the case of "trap detector", this is the threshold
    β_ratio_N2R_exploit::Float64    # For testing the effect of novelty for exploitation phases
end
function Str_Param(; num_act=4, num_obs=11, count_decay=1,
                    p_change = 0.9, ϵ = 0.01,
                    β_R=3, β_ratio_N2R=0.5, α = 0.1, Δα = 0,
                    γ_RL_R=0.9, γ_RL_N=0.9, λ_elig_R = 0.9, λ_elig_N = 0.9,
                    Q_R_0 = 0, Q_N_0 = 0, N_prior_sweep = 20,
                    w_MF_to_MB_1=1, w_MF_to_MB_2=1, w_MF_to_MB_Scaler=1,
                    γ_RL_R_MB=-1, γ_RL_N_MB=-1, C_thresh = 50, β_ratio_N2R_exploit = 0)

    Str_Param(num_act, num_obs, count_decay,
               p_change, ϵ, α, Δα,
               β_R, β_ratio_N2R, γ_RL_R, γ_RL_R,
               λ_elig_R, λ_elig_N, Q_R_0, Q_N_0,
               N_prior_sweep, w_MF_to_MB_1, w_MF_to_MB_2, w_MF_to_MB_Scaler,
               γ_RL_R_MB, γ_RL_N_MB, C_thresh, β_ratio_N2R_exploit)
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Structure for hyper parameters: Defining what learning rule is used
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
struct Str_HyperParam
    gamma_form::String                  # Form of adaptation
    surprise_form::String               # Definition of surprise
    novelty_form::String                # Definition of novelty
    exploration_form::String            # Whether uses novelty or MB_Curiosity for exploration
    update_form::String                 # Update rule
    update_form_model_free::String      # Update rule for model free
    MF_surprise_modul::String           # Surprise modulation form of MF
    reward_func::Array{Float64,1}       # Reward indecator
    First_Epi_reward_update::Bool       # Whether update MB Q-values in 1st epi
end
function Str_HyperParam(; gamma_form="satur",
                         surprise_form="BF",
                         novelty_form="log-count",
                         exploration_form="Novelty",
                         update_form="VarSMiLe",
                         update_form_model_free="Q_learning",
                         MF_surprise_modul = "Adapt",
                         reward_func=(R = zeros(11); R[1] = 1; R),
                         First_Epi_reward_update=false)

    Str_HyperParam(gamma_form, surprise_form, novelty_form, exploration_form,
                   update_form, update_form_model_free, MF_surprise_modul,
                   reward_func, First_Epi_reward_update)
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Structure for Inputs
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
mutable struct Str_Input
    obs::Array{Int,1}
    act::Array{Int,1}
end
function Str_Input(; obs=zeros(1), act=zeros(1))

    Str_Input(obs, act)
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Structure for Agent
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
mutable struct Str_Agent
    Sub::Int
    Epi::Int
    Block::Int

    Param::Str_Param
    HyperParam::Str_HyperParam
    Input::Str_Input

    Π::Array{Float64,3}                 # Belief = 3D matrix of αs
    Counter::Array{Float64,2}           # Counts for computing novelty
    Novelty::Array{Float64,2}

    Q_N_MF::Array{Float64,2}            # Model free Q-values for novelty
    Q_R_MF::Array{Float64,2}            # Model free Q-values for reward
    E_N::Array{Float64,2}               # Eligibility trace for novelty
    E_R::Array{Float64,2}               # Eligibility trace for reward

    Q_N_MB::Array{Float64,2}            # Model based Q-values for novelty
    Q_R_MB::Array{Float64,2}            # Model based Q-values for reward
    U_N::Array{Float64,2}               # U-values for prioritized sweeping
    U_R::Array{Float64,2}               # U-values for prioritized sweeping

    Surprise::Array{Float64,1}          # Surprise time-seriese
    Novelty_Seq::Array{Float64,1}       # Novelty time-seriese
    δ_N_Seq::Array{Float64,1}           # NPE time-seriese
    δ_R_Seq::Array{Float64,1}           # RPE time-seriese
    U_N_Seq::Array{Float64,1}           # U_N time-seriese
    U_R_Seq::Array{Float64,1}           # U_R time-seriese
    Gamma::Array{Float64,1}             # γ time-seriese
    t_since_change::Array{Float64,1}    # the expected sequence of time-points since the last change-point

    State::Array{Int,1}                 # "State" time-seriese, defined as combinations of actions and states

    P_Action::Array{Float64,2}          # Output of policy

    LogL::Float64                       # Log likelihood
    LogL_Unif::Float64                  # Log likelihood for random choice model

    Acc_Rate::Float64                   # Accuracy rate

    Est_Act::Array{Int,1}               # The sequence of estimated actions
    Est_Act_P::Array{Float64,1}         # The sequence of probabilities of estimated actions

    Sort_Dist_Model::Array{Float64,2}   # Sorted P_Action
    Sort_Dist_Data::Array{Float64,2}    # Sorted P_Action corresponding to data

    First_Epi::Bool                     # Indicator of 1st episode
end
function Str_Agent(; Sub=1, Epi=1, Block=1,
                    Param=Str_Param(),
                    HyperParam=Str_HyperParam(), Input=Str_Input(),
                    Π=ones(11,44,1), Counter=zeros(11,1), Novelty=zeros(11,1),
                    Q_N_MF=zeros(44,1), Q_R_MF=zeros(44,1),
                    E_N=zeros(44,1), E_R=zeros(44,1),
                    Q_N_MB=zeros(44,1), Q_R_MB=zeros(44,1),
                    U_N=zeros(11,1), U_R=zeros(11,1),
                    Surprise=-ones(1), Novelty_Seq=-ones(1),
                    δ_N_Seq=-ones(1), δ_R_Seq=-ones(1),
                    U_N_Seq=-ones(1), U_R_Seq=-ones(1),Gamma=-ones(1),
                    t_since_change = ones(1),
                    State=zeros(1), P_Action=ones(4,1)/4,
                    LogL=0, LogL_Unif=0, Acc_Rate=0,
                    Est_Act=zeros(1), Est_Act_P=zeros(1),
                    Sort_Dist_Model=ones(4,1)/4, Sort_Dist_Data=ones(4,1)/4,
                    First_Epi=true)

   Str_Agent(Sub, Epi, Block, Param, HyperParam, Input, Π,
             Counter, Novelty, Q_N_MF, Q_R_MF, E_N, E_R,
             Q_N_MB, Q_R_MB, U_N, U_R,
             Surprise, Novelty_Seq, δ_N_Seq, δ_R_Seq, U_N_Seq, U_R_Seq,
             Gamma, t_since_change, State, P_Action, LogL, LogL_Unif, Acc_Rate,
             Est_Act, Est_Act_P, Sort_Dist_Model, Sort_Dist_Data,
             First_Epi)
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Structure for the graph of the environment
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
struct Str_EnvGraph
    TranMat::Array{Int,2}
    InitObs::Array{Int,2}
    Goal::Int
end
function Str_EnvGraph(; TranMat = zeros(2), InitObs = zeros((1,2)), Goal=0)

    Str_EnvGraph(TranMat, InitObs, Goal)
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Structure of the agent, switible for saving as .MAT for MATLAB and EEG analysis
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
mutable struct Str_Agent_MAT
    Sub::Int
    Epi::Int
    Block::Int

    HyperParam::Str_HyperParam
    Input::Str_Input

    Pi::Array{Float64,3}
    Counter::Array{Float64,2}
    Novelty::Array{Float64,2}

    Q_N_MF::Array{Float64,2}
    Q_R_MF::Array{Float64,2}
    E_N::Array{Float64,2}
    E_R::Array{Float64,2}

    Q_N_MB::Array{Float64,2}
    Q_R_MB::Array{Float64,2}
    U_N::Array{Float64,2}
    U_R::Array{Float64,2}

    Surprise::Array{Float64,1}
    Novelty_Seq::Array{Float64,1}
    delta_N_Seq::Array{Float64,1}
    delta_R_Seq::Array{Float64,1}
    U_N_Seq::Array{Float64,1}
    U_R_Seq::Array{Float64,1}
    Gamma::Array{Float64,1}

    State::Array{Int,1}

    P_Action::Array{Float64,2}

    LogL::Float64
    LogL_Unif::Float64

    Acc_Rate::Float64

    Est_Act::Array{Int,1}
    Est_Act_P::Array{Float64,1}

    Sort_Dist_Model::Array{Float64,2}
    Sort_Dist_Data::Array{Float64,2}

    First_Epi::Bool
end
function Str_Agent_MAT(; Agent=Str_Agent)

   Str_Agent_MAT(Agent.Sub, Agent.Epi, Agent.Block, Agent.HyperParam, Agent.Input,
                 Agent.Π, Agent.Counter, Agent.Novelty, Agent.Q_N_MF, Agent.Q_R_MF, Agent.E_N, Agent.E_R,
                 Agent.Q_N_MB, Agent.Q_R_MB, Agent.U_N, Agent.U_R,
                 Agent.Surprise, Agent.Novelty_Seq, Agent.δ_N_Seq, Agent.δ_R_Seq, Agent.U_N_Seq, Agent.U_R_Seq,
                 Agent.Gamma, Agent.State, Agent.P_Action, Agent.LogL, Agent.LogL_Unif, Agent.Acc_Rate,
                 Agent.Est_Act, Agent.Est_Act_P, Agent.Sort_Dist_Model, Agent.Sort_Dist_Data,
                 Agent.First_Epi)
end

mutable struct Str_EEG
    Sub::Int
    Epi::Int
    Block::Int
    obs::Array{Int,1}
    EEG::Array{Float64,2}
end
function Str_EEG(; Sub=1, Epi=1, Block=1,
                   obs=zeros(1), EEG=zeros(1,1) )
    Str_EEG(Sub, Epi, Block,obs, EEG)
end


text:
"""
    title(
        [scene=current_scene(), ], string;
        align = (:center, :bottom), textsize = 30, parent = Scene(), formatter = string, kw...
    )

Add a title with content `string` to `scene`.  Pass a `Function` to the `formatter` kwarg
if the value passed to `string` isn't actually a String.
"""
function title(scene, tstring; align = (:center, :bottom), textsize = 30, parent = Scene(), formatter = string, kw...)

    string = to_node(tstring)

    pos = lift(pixelarea(scene)) do area
        x = widths(area)[1] ./ 2
        Vec2f0(x, 10) # offset 10px, to give it some space
    end
    t = text(
        lift(formatter, string),
        position = pos,
        camera = campixel!,
        raw = true, align = align,
        textsize = textsize;
        kw...
    )
    hbox(scene, t; parent = parent)
end

function title(string; kw...)
    scene = current_scene()
    title(scene, string; kw...)
end


text:
module InvertedPendulum

using MDPs
using DelimitedFiles
using PGFPlots

export state_space, action_space, transition, reward, simulate,
       viz_policy, viz_trajectory
export PMIN, PMAX, VMIN, VMAX, T

# mdp
const MP = 50
const MV = 50
const N = 1000

# state
const PMIN = float(-pi)
const PMAX = float(pi)
const VMIN = -10.0
const VMAX = 10.0

# action
const AMIN = -1.0
const AMAX = 1.0

# transition
const A = 1.0
const B = 1.0
const DT = MV / MP * pi / VMAX  # time intervals
const NOISE = deg2rad(0.5)

# reward
const K = 1.0
const R = 0.1

# test
const T = 200

# visualization
const EPS = 1e-4

function state_space(mp::Int64=MP, mv::Int64=MV)
    p = range(PMIN, stop=PMAX, length=mp)
    v = range(VMIN, stop=VMAX, length=mv)
    return RectangleGrid(p, v)
end

function action_space(n::Int64=N)
    return range(AMIN, stop=AMAX, length=n)
end

function wrap_around(angle::Float64)
    while angle < -pi
        angle += 2 * pi
    end
    while angle > pi
        angle -= 2 * pi
    end
    return angle
end

function transition(state::Vector{Float64}, action::Float64, noise::Float64=0.0)
    snext = zeros(length(state))
    snext[1] = wrap_around(state[1] + DT * state[2])
    snext[2] = state[2] + DT * (A * sin(state[1]) - B * state[2] + action) + noise
    return snext
end

function reward(state::Vector{Float64}, action::Float64)
    return exp(K * (cos(state[1]) - 1)) - R * action^2 - 1
end

function simulate(mdp::MDP, policy::Policy, state::Vector{Float64})
    trajectory = zeros(T, dimensions(mdp.S))
    actions = zeros(T)
    for t = 1:T
        trajectory[t, :] = state
        if state[1] == 0.0
            trajectory = trajectory[1:t, :]
            actions = actions[1:t - 1]
            break
        end
        action = policy!(policy, get_belief(mdp, state))
        state = mdp.transition(state, action, NOISE * randn())
        actions[t] = action
    end
    return trajectory, actions
end

function viz_policy(mdp::MDP, policy::Policy, title::String="Policy heatmap", saveplot::Bool=false, name::String="")
    function getmap(p::Float64, v::Float64)
        return policy!(policy, get_belief(mdp, [p, v]))
    end

    p = Axis([Plots.Image(getmap, (PMIN + EPS, PMAX - EPS),
                      (VMIN + EPS, VMAX - EPS),
                      xbins = 250, ybins = 250,
                      colormap = ColorMaps.Named("Blues"))],
          width="14cm", height="14cm",
          xlabel="angle", ylabel="angular speed",
          title=title)
    p
    if saveplot
      PGFPlots.save(name, p)
    end
end

function rad2deg_vec(trajectory::Vector{Float64})
    l = length(trajectory)
    out = zeros(l)
    for i = 1:l
        out[i] = rad2deg(trajectory[i])
    end
    return out
end

function viz_trajectory(trajectory::Matrix{Float64}, actions::Vector{Float64},
                        title1::String="angular position over time", title2::String="control input over time",
                        saveplot::Bool=false, name::String="")
    g = GroupPlot(2, 1, groupStyle="horizontal sep=1.5cm")
    push!(g, Axis([Plots.Linear(rad2deg_vec(trajectory[:, 1]))],
                  width="12cm", height="12cm",
                  xlabel="time", ylabel="angle (deg)",
                  title=title1))
    push!(g, Axis([Plots.Linear(actions)],
                  width="12cm", height="12cm",
                  xlabel="time", ylabel="input",
                  title=title2))
    g
    if saveplot
      PGFPlots.save(name, g)
    end
end

end

max_stars_repo_name: tsamtsam/SV-RL

max_issues_repo_name: tsamtsam/SV-RL

max_forks_repo_name: tsamtsam/SV-RL


text:
using Test
using JuMP
using MultivariateMoments
using PolyJuMP
using DynamicPolynomials

function zero_polynomial_test(optimizer::MOI.AbstractOptimizer,
                              config::MOI.Test.TestConfig)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)
    @variable(model, β ≤ 1)
    @variable(model, γ)

    @polyvar x y
    cref = @constraint(model, (α - β) * x * y == 0)
    cγ = @constraint(model, γ * x * y in PolyJuMP.ZeroPoly())

    @objective(model, Max, α + γ)
    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 1.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 1.0 atol=atol rtol=rtol
    @test value(β) ≈ 1.0 atol=atol rtol=rtol
    @test value(UpperBoundRef(β)) ≈ 1.0 atol=atol rtol=rtol
    @test value(cref) isa MultivariatePolynomials.AbstractPolynomial{Float64}
    @test value(cref) ≈ 0.0 * x * y atol=atol rtol=rtol

    @test dual_status(model) == MOI.FEASIBLE_POINT
    @test dual(UpperBoundRef(β)) ≈ -1.0 atol=atol rtol=rtol

    for μ in [dual(cref), moments(cref), dual(cγ), moments(cγ)]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 1
        @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol
        @test monomial(moments(μ)[1]) == x*y
    end

    F = MOI.VectorAffineFunction{Float64}
    S = PolyJuMP.ZeroPolynomialSet{FullSpace, MB.MonomialBasis}
    @test_broken Set(MOI.get(model, MOI.ListOfConstraints())) == Set([
        (MOI.SingleVariable, MOI.LessThan{Float64}), (F, S), (MOI.VectorOfVariables, S)])
    @testset "Delete" begin
        test_delete_bridge(model, cref, 3, ((F, MOI.Zeros, 0),))
    end
end

linear_tests["zero_polynomial"] = zero_polynomial_test


text:
# This file contains functions for economic optimization problems.
# These functions are intended for internal use and package extensions.
function deamincost(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}; rts::Symbol = :VRS, dispos::Symbol = :Strong,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end

    if dispos != :Strong && dispos != :Weak
        throw(ArgumentError("`disposX` must be :Strong or :Weak"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx

    Xtarget = zeros(n,m)
    Ytarget = Y[:,:]
    clambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        y0 = Y[i,:]
        w0 = W[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, Xeff[1:m])
        @variable(deamodel, lambda[1:n] >= 0)

        @objective(deamodel, Min, sum(w0[j] .* Xeff[j] for j in 1:m))

        @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) <= Xeff[j])

        if dispos == :Strong
            @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) >= y0[j])
        elseif dispos == :Weak
            @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) == y0[j])
        end

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        Xtarget[i,:]  = JuMP.value.(Xeff)
        clambdaeff[i,:] = JuMP.value.(lambda)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    return Xtarget::Matrix, clambdaeff::SparseMatrixCSC{Float64, Int64}

end

function deamaxrevenue(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    P::Union{Matrix,Vector}; rts::Symbol = :VRS, dispos::Symbol = :Strong,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    if dispos != :Strong && dispos != :Weak
        throw(ArgumentError("`disposY` must be :Strong or :Weak"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx

    Xtarget = X[:,:]
    Ytarget = zeros(n,s)
    rlambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        p0 = P[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, Yeff[1:s])
        @variable(deamodel, lambda[1:n] >= 0)

        @objective(deamodel, Max, sum(p0[j] .* Yeff[j] for j in 1:s))

        if dispos == :Strong
            @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) <= x0[j])
        elseif dispos == :Weak
            @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) == x0[j])
        end

        @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) >= Yeff[j])

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        Ytarget[i,:]  = JuMP.value.(Yeff)
        rlambdaeff[i,:] = JuMP.value.(lambda)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    return Ytarget::Matrix, rlambdaeff::SparseMatrixCSC{Float64, Int64}

end

function deamaxprofit(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, P::Union{Matrix,Vector};
    optimizer::Union{DEAOptimizer,Nothing} = nothing)

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)
    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx

    Xtarget = zeros(n,m)
    Ytarget = zeros(n,s)
    plambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        w0 = W[i,:]
        p0 = P[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, Xeff[1:m])
        @variable(deamodel, Yeff[1:s])
        @variable(deamodel, lambda[1:n] >= 0)

        @objective(deamodel, Max, (sum(p0[j] .* Yeff[j] for j in 1:s)) - (sum(w0[j] .* Xeff[j] for j in 1:m)))

        @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) <= Xeff[j])
        @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) >= Yeff[j])

        @constraint(deamodel, sum(lambda) == 1)

        # Optimize and return results
        JuMP.optimize!(deamodel)

        Xtarget[i,:]  = JuMP.value.(Xeff)
        Ytarget[i,:]  = JuMP.value.(Yeff)
        plambdaeff[i,:] = JuMP.value.(lambda)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    return Xtarget::Matrix, Ytarget::Matrix, plambdaeff::SparseMatrixCSC{Float64, Int64}

end


text:
using NumericalAnalysis
using Test

@testset "NumericalAnalysis.jl" begin
    #2x + 4y
    # Write your tests here.
end

@testset "basic test" begin
    @test 1<=ξ(1, 2)<=2
    @test absolute_error(1.0, 2) == 1.0
    @test relative_error(1, 2) == 1.0
    @test NthDerivative(exp, 4, 4)  == exp(4)
    @test NthDerivative(x->x^3, 4, 4) == 0
    @test 0.9<=TaylorPolynomials(cos, 0.1, 0, 6)<=1.0

end


text:
using CUDA, Libdl, Tar, SHA

NVCC = nothing
CFLAGS = Sys.iswindows() ? ["/Ox","/LD"] : ["-O3","-Wall","-fPIC","-std=c++11"]
NVCCFLAGS = ["-O3","--use_fast_math","-Wno-deprecated-gpu-targets","--default-stream", "per-thread"]

# For compatibility with older chips, from NVIDIA samples Makefiles:
SMS = [35, 37, 50, 52, 53, 60, 61, 62, 70, 72, 75, 80]
for sm in SMS; push!(NVCCFLAGS,"-gencode=arch=compute_$sm,code=sm_$sm"); end
push!(NVCCFLAGS,"-gencode=arch=compute_$(SMS[end]),code=compute_$(SMS[end])")

const OBJEXT = Sys.iswindows() ? ".obj" : ".o"
const LIBKNET8 = "libknet8."*Libdl.dlext
const DLLEXPORT = Sys.iswindows() ? "__declspec(dllexport)" : "" # this needs to go before function declarations
inforun(cmd)=(@info(cmd);run(cmd))

# NVCC must be on the path:

try
    run(`nvcc -V`)
catch
    error("nvcc not found, libknet8 will not be built.")
end

NVCC = "nvcc"
push!(NVCCFLAGS,"--compiler-options",join(CFLAGS,' '))

# If CUDA is available add architecture optimization flags
# Uncomment this for better compiler optimization
# We keep it commented to compile for multiple gpu types

# if NVCC !== nothing && CUDA.functional()
#     dev = CuDevice(0)
#     cap = capability(dev)
#     arch = "sm_$(cap.major)$(cap.minor)"
#     push!(NVCCFLAGS,"--gpu-architecture",arch)
# end

# Build scripts

function build_nvcc()
    SRC = [("cuda1","gamma","../knetarrays/unary","ops"),
           ("cuda01","../knetarrays/binary","ops"),
           ("cuda11","../knetarrays/binary","../knetarrays/unary","ops"),
           ("cuda111","../knetarrays/binary","../knetarrays/unary","ops"),
           ("cuda12","../knetarrays/binary","ops"),
           ("cuda13","../knetarrays/binary","ops"),
           ("cuda16","../knetarrays/binary","ops"),
           ("cuda17","../knetarrays/binary","ops"),
           ("cuda20","../knetarrays/reduction","ops"),
           ("cuda21","../knetarrays/reduction","ops"),
           ("cuda22","../knetarrays/reduction","ops"),
           ("relu","ops"),
           ]

    OBJ = []

    include("ops.jl")

    for names in SRC
        for name in names
            if !isfile("$name.jl")
                error("$name.jl not found")
            end
        end
        name1 = names[1]
        if !isfile("$name1.cu") || any(d->(mtime("$d.jl") > mtime("$name1.cu")), names)
            @info("$name1.jl")
            include("$name1.jl")     # outputs name1.cu
        end
        obj = name1*OBJEXT
        if !isfile(obj) || mtime("$name1.cu") > mtime(obj)
            inforun(`$NVCC $NVCCFLAGS -c $name1.cu`)
        end
        push!(OBJ,obj)
    end

    if any(f->(mtime(f) > mtime(LIBKNET8)), OBJ)
        inforun(`$NVCC $NVCCFLAGS --shared -o $LIBKNET8 $OBJ`)
    end
end

function build()
    @assert NVCC !== nothing "no compilers found, libknet8 will not be built."
    build_nvcc()
    if Sys.iswindows()
        run(`tar cf libknet8.tar libknet8.dll libknet8.lib libknet8.exp`)
    else
        run(`tar cf libknet8.tar $LIBKNET8`)
    end
    sha1 = Tar.tree_hash("libknet8.tar")
    run(`gzip libknet8.tar`)
    sha2 = open("libknet8.tar.gz") do f; bytes2hex(sha256(f)); end
    @info "git-tree-sha1 = \"$sha1\""
    @info "sha256 = \"$sha2\""
end

build()


text:
module SchoolOutbreak
using Distributions, LinearAlgebra, StatsFuns
Base.isless(n::Nothing,f)=false;Base.isless(f, n::Nothing)=true;Base.isless(n1::Nothing, n2::Nothing)=false

# utils
updateArray!(x::A,y) where A<:AbstractArray =begin x.=y end
updateArray!(x::A,y) where A<:AbstractArray{<:AbstractArray} =begin updateArray!.(x,y) end
interval(x,lo=-Inf,up=Inf)=min(max(x,lo),up)

function getstdistance(from,to)
    Int.(to.stratum.==from.stratum) |> cumprod |> sum |> x->x+one(x) # look for the match level betweeen stX and stY
end
function getstratum(student)
    student.stratum
end
function getlogNbylevels(stratum, nclassmat, nstudentsmat, denom=ones(5))
    begin
        N=Vector{Float64}(undef,5)
        N[1]=sum(nstudentsmat)/denom[1] |> log # outside school
        N[2]=(sum(@view nstudentsmat[stratum[1],:])/sum(nclassmat[stratum[1],:]))/denom[2] |> log # different grade
        N[3]=nstudentsmat[stratum[1:2]...]/nclassmat[stratum[1:2]...]/denom[3] |> log # different class same grade
        N[4] = N[5] = nstudentsmat[stratum[1:2]...]/nclassmat[stratum[1:2]...]/denom[4] |> log # same class
        N .|>fill
    end
end
function getlogMbylevels(stratum, nclassmat, nstudentsmat, denom=ones(5))
    begin
        N=Vector{Float64}(undef,5)
        N[1]=sum(nstudentsmat)/denom[1] |> log # outside school
        N[2]=(sum(@view nclassmat[stratum[1],:]))/denom[2] |> log # different grade
        N[3]=nclassmat[stratum[1:2]...]/denom[3] |> log # different class same grade
        N[4] = N[5] = nclassmat[stratum[1:2]...]/denom[4] |> log # same class
        N .|>fill
    end
end
function linregval(covar,coef)
    exp(dot(first.(covar),first.(coef)))
end

## Estimation
# Types and related functions
Scalar{T}=Array{T,0}
AbstractScalar{T}=AbstractArray{T,0}
Base.:*(x::Scalar,y::Scalar)=x[]*y[]
Base.:*(x::Number,y::Scalar)=x*y[]
Base.:*(x::Scalar,y::Number)=x[]*y
mutable struct
    Student{I<:Integer, R<:Real, VI<:AbstractVector,VR<:AbstractVector,SR<:AbstractScalar}
    isinfected::Bool
    onset::I
    stratum::VI
    suscovar::VR
    infcovar::VR
    logN::Vector{SR}
    logM::Vector{SR}
    HHllkh::SR
    susceptibility::SR
    infectiousness::SR
    pdfcommunity::SR
    cdfcommunity::SR
    βs::Vector{SR}
    γs::Vector{SR}
    δs::Vector{SR}
    pdfgtime::Vector{SR}
    cdfgtime::Vector{SR}
    sampleweight::R
end
function studentsfromdata(parameters,isinfecteds,onsets,strata,
            suscovars,infcovars,nclasses,nstudents,HHlkhs,samplingrate)
    logNs=getindex.(Ref(parameters.logN),strata)
    logMs=getindex.(Ref(parameters.logM),strata)
    students=Student.(isinfecteds,onsets,strata,suscovars,infcovars,logNs,logMs, fill.(HHlkhs),
        fill.(linregval.(suscovars,Ref(parameters.suscoef))),fill.(linregval.(infcovars,Ref(parameters.infcoef))),
        fill.(cdf.(parameters.rcomdist,onsets).-cdf.(parameters.rcomdist,onsets.-1)), fill.(cdf.(parameters.rcomdist,onsets.-1)),
        Ref(parameters.β), Ref(parameters.γ),Ref(parameters.δ), Ref(parameters.pdfgtime), Ref(parameters.cdfgtime), 1 ./samplingrate)
end
struct StudentsData{VStudent,NT1,NT2,MR<:Matrix,VR<:Vector,V<:Vector,M<:Matrix}
    members::VStudent
    parameters::NT1
    data::NT2
    isinfected::Vector{Bool}
    βmat::MR
    γmat::MR
    δmat::MR
    logNmat::MR
    logMmat::MR
    pdfgtimemat::MR
    cdfgtimemat::MR
    pdfcommunity::VR
    cdfcommunity::VR
    infectiousness::VR
    susceptibility::VR
    HHllkh::VR
    sampleweight::V
    temps::NamedTuple{(:cdfmat,:pdfmat,:pdfvec),Tuple{M,M,V}}
    StudentsData(members, parameters = NamedTuple(), data = NamedTuple())=new{typeof(members),typeof(parameters),typeof(data),Matrix{Scalar{Float64}},Vector{Scalar{Float64}},Vector{Float64},Matrix{Float64}}(members, parameters, data, makemats(members)...)
end
function makemats(members::T where T<:Vector{<:Student})
    (isinfected=[from.isinfected for from in members],
    βmat=[to.βs[getstdistance(from,to)] for from in members, to in members],
    γmat=[to.γs[getstdistance(from,to)] for from in members, to in members],
    δmat=[to.δs[getstdistance(from,to)] for from in members, to in members],
    logNmat=[to.logN[getstdistance(from,to)] for from in members, to in members],
    logMmat=[to.logM[getstdistance(from,to)] for from in members, to in members],
    pdfgtimemat=[from.pdfgtime[interval(to.onset-from.onset+1,1,length(from.pdfgtime))] for from in members, to in members],
    cdfgtimemat=[from.cdfgtime[interval(to.onset-from.onset+1,1,length(from.cdfgtime))] for from in members, to in members],
    pdfcommunity=[to.pdfcommunity for to in members],
    cdfcommunity=[to.cdfcommunity for to in members],
    infectiousness=[from.infectiousness for from in members],
    susceptibility=[to.susceptibility for to in members],
    HHllkh=[to.HHllkh for to in members],
    sampleweight=[to.sampleweight for to in members],
    temps=(cdfmat=[0.0 for from in members[getfield.(members,:isinfected)], to in members],
            pdfmat=[0.0 for from in members[getfield.(members,:isinfected)], to in members[getfield.(members,:isinfected)]],
            pdfvec=[0.0 for from in members[getfield.(members,:isinfected)]]))
end

# main functions
function makeparameters(βs,γs,δs,rcom,suscoef,infcoef,gtimedist,gtimemax,rcomdist,strata,nclasses,nstudents,Ndenom,Mdenom,overdispersion)
    uniquestrata=unique(eachcol(strata))
    logN=getlogNbylevels.(uniquestrata,Ref(nclasses),Ref(nstudents),Ref(Ndenom))
    dlogN=Dict{eltype(uniquestrata),eltype(logN)}(uniquestrata.=>logN)
    logM=getlogMbylevels.(uniquestrata,Ref(nclasses),Ref(nstudents),Ref(Mdenom))
    dlogM=Dict{eltype(uniquestrata),eltype(logM)}(uniquestrata.=>logM)
    β=fill.(βs)
    γ=fill.(γs)
    δ=fill.(δs)
    pdfgtime=[0;cdf.(gtimedist,1:gtimemax).-cdf.(gtimedist,0:gtimemax-1)]
    cdfgtime=[0;cdf.(gtimedist,0:gtimemax-1)]
    push!(pdfgtime,0.0);push!(cdfgtime,1.0) # add 0 to the end
    rcom=fill.(rcom)
    (β=β, γ=γ,δ=δ, rcom=rcom,suscoef=suscoef,infcoef=infcoef, pdfgtime=fill.(pdfgtime),cdfgtime=fill.(cdfgtime), rcomdist=fill(rcomdist),logN=dlogN,logM=dlogM,Ndenom=Ndenom,Mdenom=Mdenom,overdispersion=fill.(overdispersion))
end
function ll!(studentsdata::StudentsData)
    updatesusinf!.(studentsdata.members, Ref(studentsdata.parameters)) # reflect covariates

    @views studentsdata.temps.cdfmat.=first.(studentsdata.βmat[studentsdata.isinfected,:]).*exp.(.-first.(studentsdata.γmat[studentsdata.isinfected,:]).*first.(studentsdata.logNmat[studentsdata.isinfected,:]).-first.(studentsdata.δmat[studentsdata.isinfected,:]).*first.(studentsdata.logMmat[studentsdata.isinfected,:])) # βN⁻ᵞ
    @views studentsdata.temps.pdfmat.=studentsdata.temps.cdfmat[:,studentsdata.isinfected] # copy for pdf computation

    @views studentsdata.temps.cdfmat.*=first.(studentsdata.cdfgtimemat[studentsdata.isinfected,:])

    @views ll= - first.(studentsdata.infectiousness[studentsdata.isinfected])'*(studentsdata.temps.cdfmat)*(first.(studentsdata.susceptibility).*first.(studentsdata.sampleweight)) # siβN⁻ᵞcdf
    ll-= (studentsdata.parameters.rcom)sum(first.(studentsdata.cdfcommunity).*first.(studentsdata.susceptibility).*first.(studentsdata.sampleweight))

    studentsdata.temps.pdfmat.*=first.(studentsdata.pdfgtimemat[studentsdata.isinfected,studentsdata.isinfected])
    studentsdata.temps.pdfvec.= (studentsdata.temps.pdfmat)'*first.(studentsdata.infectiousness[studentsdata.isinfected])

    @views studentsdata.temps.pdfvec.+=(studentsdata.parameters.rcom).*first.(studentsdata.pdfcommunity[studentsdata.isinfected])
    @views studentsdata.temps.pdfvec.*=first.(studentsdata.susceptibility[studentsdata.isinfected])
    @views studentsdata.temps.pdfvec.=log.((1 .-exp.(.-studentsdata.temps.pdfvec)).*(1 .-exp.(first.(studentsdata.HHllkh[studentsdata.isinfected]))).+exp.(first.(studentsdata.HHllkh[studentsdata.isinfected])))

    ll+=dot(studentsdata.temps.pdfvec, first.(@view studentsdata.sampleweight[studentsdata.isinfected]))
    ll
end
function llfunc!(studentsdatasets::AS) where AS<:Array{<:StudentsData}
    sum((ll!(dataset) for dataset in studentsdatasets))
end
function updatesusinf!(student::Student,parameters::NT where NT<:NamedTuple)
    student.susceptibility.=linregval(student.suscovar,parameters.suscoef)
    student.infectiousness.=linregval(student.infcovar,parameters.infcoef)
end

## Simulation

mutable struct StudentPair{I<:Integer,SR<:AbstractScalar,VR<:AbstractVector,ST<:Student}
    isinfected::Bool
    serialinterval::I
    stratadistance::I
    γ::SR
    δ::SR
    pairwiseβ::SR
    pairwiselogN::SR
    pairwiselogM::SR
    infectiousness::SR
    susceptibility::SR
    pdfgtime::VR
    cdfgtime::VR
    from::ST
    to::ST
    StudentPair(from,to)=begin
        stdistance=getstdistance(from,to)
        new{typeof(stdistance),typeof(to.susceptibility),typeof(to.pdfgtime),typeof(from)}(
            from.isinfected,
            to.onset.-from.onset,
            stdistance,
            to.γs[stdistance],
            to.δs[stdistance],
            to.βs[stdistance],
            to.logN[stdistance],
            to.logM[stdistance],
            from.infectiousness,
            to.susceptibility,
            from.pdfgtime,
            from.cdfgtime,
            from, to
        )
    end
end

mutable struct Students{NT1<:NamedTuple,D<:Dict,ST<:Student,SP<:StudentPair}
    members::Vector{ST}
    pairs::Matrix{SP}
    parameters::NT1
    data::D
    Students(members, parameters = NamedTuple(), data = Dict())=new{typeof(parameters),typeof(data),eltype(members),typeof(StudentPair(members[1],members[1]))}(members, StudentPair.(members,permutedims(members[:,:],(2,1))), parameters, data)
end
function Rmean(students)
    size(students.pairs,1)mean((pair->pair.pairwiseβ[]*exp(-pair.pairwiselogN[]*pair.γ[]-pair.pairwiselogM[]*pair.δ[])),students.pairs)
end
function llpair(pair::StudentPair)
    first(.-pair.pairwiseβ.*exp.(.-pair.γ.*pair.pairwiselogN.-pair.δ.*pair.pairwiselogM)
        .*(@view pair.cdfgtime[interval(pair.serialinterval+1,1,length(pair.cdfgtime))]).*pair.infectiousness.*pair.susceptibility
  .+ (pair.isinfected ? (
        begin
            prob=pair.pdfgtime[interval(pair.serialinterval+1,1,length(pair.cdfgtime))][]
            log((1-exp(-prob)))#*(1-exp(pair.to.HHllkh[]))+exp(pair.to.HHllkh[]))
        end) : 0) )
    ## sampling, hhfoi
end
function initialize!(students::Students, updateparameters=NamedTuple())
    for par in keys(updateparameters)
        updateArray!(students.parameters[par],updateparameters[par])
    end
    onsets=Vector{Union{Nothing,typeof(students.members[1].onset)}}(undef,length(students.members))
    onsets.=getfield.(students.members,:onset)
    onsets[.!getfield.(students.members,:isinfected)].=nothing
    setfield!.(students.members,:isinfected,false)
    setfield!.(students.members,:onset,typemax(Int16)|>Int)
    (x->x.infcovar[end]=0.0).(students.members) # remove winterbreak flag
    setfield!.(students.pairs,:isinfected,false)
    setfield!.(students.pairs,:serialinterval,typemax(Int16)|>Int)
    updatesusinf!.(students.members, Ref(students.parameters))
    onsets # returns the currently recorded onset dates
end
function initialcases!(students::Students, id, onset)
    @views setfield!.(students.members[id[1:min(length(id),length(onset))]],:onset,onset[1:min(length(id),length(onset))])
    setfield!.(students.members[@view id[1:min(length(id),length(onset))]],:isinfected,true)
end
function communitytransmission!(students::Students)
    infected=rand.(fill(Bernoulli(1-exp(-first(students.parameters.rcom))),length(students.members)))
    setfield!.(students.members[infected], :isinfected, true)
    setfield!.(students.members[infected],:onset,floor.(typeof(students.members[1].onset),rand(students.parameters.rcomdist[],sum(infected))))
end
function pairtransmit!(pair, parameters, data)
    # update winterbreak
    if haskey(data,:duringbreak) && (pair.from.infcovar[end]>0)!=(pair.from.onset in data.duringbreak)
        pair.from.infcovar[end]=(pair.from.onset in data.duringbreak)
        pair.from.infectiousness.*=exp(parameters.infcoef[end]*(2pair.from.infcovar[end]-1)) #updatesusinf!(pair.from, parameters)
    end

    if rand(Bernoulli(exp(llpair(pair)))) || pair.to.onset≤pair.from.onset
        return end
    pair.isinfected=true
    pair.to.isinfected=true
    if(haskey(parameters,:overdispersion)&&0<parameters.overdispersion[]<Inf)
        pair.to.infectiousness.*=rand(Gamma(parameters.overdispersion[],1/parameters.overdispersion[]))
    end
    SI=1:length(pair.cdfgtime)-1
    lls=[begin pair.serialinterval=si
        llpair(pair) end
            for si in SI]
    lls.=exp.(lls.-logsumexp(lls))
    pair.serialinterval=rand(Categorical(lls))
    pair.to.onset=min(pair.to.onset,pair.from.onset+pair.serialinterval) # update onset
end
function transmission!(students::Students, onset::Integer, statuscheck!)
    statuscheck!(students,onset)
    infectorids=findall(x->getfield(x,:onset) == onset ,students.members)
    @views pairtransmit!.(students.pairs[infectorids,:], Ref(students.parameters), Ref(students.data))
    return(infectorids)
end
function transmission!(students::Students, onsets, warning, statuscheck!)
    if hasfield(typeof(statuscheck!),:closuretime) updateArray!(statuscheck!.closuretime,0) end
    incidence=transmission!.(Ref(students), onsets, statuscheck!)
    if warning && ((x->ifelse(x.isinfected, x.onset, -Inf)).(students.members) |>maximum > last(onsets)) @warn "Simulation ended halfway (time step "*string(last(onsets))*", unprocessed onset at "* ((x->ifelse(x.isinfected, x.onset, -Inf)).(students.members) |>maximum |> string) *"). Restart simulation from the following step if necessary." end
    if hasfield(typeof(statuscheck!),:closuretime)
            return((incidence=incidence, closuretime=(x->min.(1,x)).(statuscheck!.closuretime))) end
    return(incidence)
end
function posterior(chain, parameter::Symbol)
    mat=chain[:,parameter,1].value[:,:,1]
    names=chain[:,parameter,1].names
    (names=names ,samples=mat)
end
function posterior(chain, parameter::Symbol,iter)
    mat=chain[(x->first(x):last(x))(chain.range[iter]),parameter,1].value[:,:,1]
    names=chain[:,parameter,1].names
    (names=names ,samples=mat)
end
function simulateoutbreaks!(students, parameters, times; initialize=true, initcases=0, warning=false, statuscheck! = (x...)->nothing)
    if initialize
        onsets=SchoolOutbreak.initialize!.(students,Ref(parameters))
        seedonsets=getindex.(sort.(onsets),range.(1,min.(initcases,length.(onsets)),step=1))#minimum.(onsets)
        initialcases=findall.((x->!isnothing(x) && x in seedonset for seedonset in seedonsets),onsets)
        SchoolOutbreak.initialcases!.(students,initialcases,seedonsets)
        SchoolOutbreak.communitytransmission!.(students)
    end
    result=SchoolOutbreak.transmission!.(students,Ref(times), warning, statuscheck!)
    #incidences=(x->length.(x)).(result)
    #totalincidence=sum(incidences)
    result
end
function incidence(outbreak)
    length(outbreak)
end
function incidence(outbreak::AbstractArray{<:AbstractArray})
    incidence.(outbreak)
end
end


text:
# ---
# title: 17. Letter Combinations of a Phone Number
# id: problem17
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Backtracking, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/letter-combinations-of-a-phone-number/description/>
# hidden: true
# ---
#
# Given a string containing digits from `2-9` inclusive, return all possible
# letter combinations that the number could represent. Return the answer in
# **any order**.
#
# A mapping of digit to letters (just like on the telephone buttons) is given
# below. Note that 1 does not map to any letters.
#
# ![](https://upload.wikimedia.org/wikipedia/commons/thumb/7/73/Telephone-
# keypad2.svg/200px-Telephone-keypad2.svg.png)
#
#
#
# **Example 1:**
#
#
#
#     Input: digits = "23"
#     Output: ["ad","ae","af","bd","be","bf","cd","ce","cf"]
#
#
# **Example 2:**
#
#
#
#     Input: digits = ""
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: digits = "2"
#     Output: ["a","b","c"]
#
#
#
#
# **Constraints:**
#
#   * `0 <= digits.length <= 4`
#   * `digits[i]` is a digit in the range `['2', '9']`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end