text:
# Utilities

# Returns a dict ("Label1" => 1, "Label2" => 2, "Label3" => 3, ...)
label_index(labels) = Dict([Pair(v => k) for (k, v) in enumerate(labels)])

## Helper function. Counts the votes.
## Returns a vector of probabilities (eg. [0.2, 0.6, 0.2]) which is in the same
## order as get_labels(classifier) (eg. ["versicolor", "setosa", "virginica"])
function compute_probabilities(labels::Vector, votes::Vector, weights=1.0)
    label2ind = label_index(labels)
    counts = zeros(Float64, length(label2ind))
    for (i, label) in enumerate(votes)
        if isa(weights, Number)
            counts[label2ind[label]] += weights
        else
            counts[label2ind[label]] += weights[i]
        end
    end
    return counts / sum(counts) # normalize to get probabilities
end

# Applies `row_fun(X_row)::Vector` to each row in X
# and returns a Matrix containing the resulting vectors, stacked vertically
function stack_function_results(row_fun::Function, X::Matrix)
    N = size(X, 1)
    N_cols = length(row_fun(_squeeze(X[1,:],1))) # gets the number of columns
    out = Array{Float64}(N, N_cols)
    for i in 1:N
        out[i, :] = row_fun(_squeeze(X[i,:],1))
    end
    return out
end

################################################################################

function _split(labels::Vector, features::Matrix, nsubfeatures::Int, weights::Vector, rng::AbstractRNG)
    if weights == [0]
        _split_info_gain(labels, features, nsubfeatures, rng)
    else
        _split_neg_z1_loss(labels, features, weights)
    end
end

function _split_info_gain(labels::Vector, features::Matrix, nsubfeatures::Int,
                          rng::AbstractRNG)
    nf = size(features, 2)
    N = length(labels)

    best = NO_BEST
    best_val = -Inf

    if nsubfeatures > 0
        r = randperm(rng, nf)
        inds = r[1:nsubfeatures]
    else
        inds = 1:nf
    end

    for i in inds
        ord = sortperm(features[:,i])
        features_i = features[ord,i]
        labels_i = labels[ord]

        hist1 = _hist(labels_i, 1:0)
        hist2 = _hist(labels_i)
        N1 = 0
        N2 = N

        for (d, range) in UniqueRanges(features_i)
            value = _info_gain(N1, hist1, N2, hist2)
            if value > best_val
                best_val = value
                best = (i, d)
            end

            deltaN = length(range)

            _hist_shift!(hist2, hist1, labels_i, range)
            N1 += deltaN
            N2 -= deltaN
        end
    end
    return best
end

function _split_neg_z1_loss(labels::Vector, features::Matrix, weights::Vector)
    best = NO_BEST
    best_val = -Inf
    for i in 1:size(features,2)
        domain_i = sort(unique(features[:,i]))
        for thresh in domain_i[2:end]
            cur_split = features[:,i] .< thresh
            value = _neg_z1_loss(labels[cur_split], weights[cur_split]) + _neg_z1_loss(labels[neg(cur_split)], weights[neg(cur_split)])
            if value > best_val
                best_val = value
                best = (i, thresh)
            end
        end
    end
    return best
end

function build_stump(labels::Vector, features::Matrix, weights=[0];
                     rng=Base.GLOBAL_RNG)
    S = _split(labels, features, 0, weights, rng)
    if S == NO_BEST
        return Leaf(majority_vote(labels), labels)
    end
    id, thresh = S
    split = features[:,id] .< thresh
    return Node(id, thresh,
                Leaf(majority_vote(labels[split]), labels[split]),
                Leaf(majority_vote(labels[neg(split)]), labels[neg(split)]))
end

function build_tree(labels::Vector, features::Matrix, nsubfeatures=0, maxdepth=-1; rng=Base.GLOBAL_RNG)
    rng = mk_rng(rng)::AbstractRNG
    if maxdepth < -1
        error("Unexpected value for maxdepth: $(maxdepth) (expected: maxdepth >= 0, or maxdepth = -1 for infinite depth)")
    elseif maxdepth==0
        return Leaf(majority_vote(labels), labels)
    end
    S = _split(labels, features, nsubfeatures, [0], rng)
    if S == NO_BEST
        return Leaf(majority_vote(labels), labels)
    end
    id, thresh = S
    split = features[:,id] .< thresh
    labels_left = labels[split]
    labels_right = labels[neg(split)]
    pure_left = all(labels_left .== labels_left[1])
    pure_right = all(labels_right .== labels_right[1])
    if pure_right && pure_left
        return Node(id, thresh,
                    Leaf(labels_left[1], labels_left),
                    Leaf(labels_right[1], labels_right))
    elseif pure_left
        return Node(id, thresh,
                    Leaf(labels_left[1], labels_left),
                    build_tree(labels_right,features[neg(split),:], nsubfeatures,
                               max(maxdepth-1, -1); rng=rng))
    elseif pure_right
        return Node(id, thresh,
                    build_tree(labels_left,features[split,:], nsubfeatures,
                               max(maxdepth-1, -1); rng=rng),
                    Leaf(labels_right[1], labels_right))
    else
        return Node(id, thresh,
                    build_tree(labels_left,features[split,:], nsubfeatures,
                               max(maxdepth-1, -1); rng=rng),
                    build_tree(labels_right,features[neg(split),:], nsubfeatures,
                               max(maxdepth-1, -1); rng=rng))
    end
end

function prune_tree(tree::LeafOrNode, purity_thresh=1.0)
    function _prune_run(tree::LeafOrNode, purity_thresh::Real)
        N = length(tree)
        if N == 1        ## a Leaf
            return tree
        elseif N == 2    ## a stump
            all_labels = [tree.left.values; tree.right.values]
            majority = majority_vote(all_labels)
            matches = find(all_labels .== majority)
            purity = length(matches) / length(all_labels)
            if purity >= purity_thresh
                return Leaf(majority, all_labels)
            else
                return tree
            end
        else
            return Node(tree.featid, tree.featval,
                        _prune_run(tree.left, purity_thresh),
                        _prune_run(tree.right, purity_thresh))
        end
    end
    pruned = _prune_run(tree, purity_thresh)
    while length(pruned) < length(tree)
        tree = pruned
        pruned = _prune_run(tree, purity_thresh)
    end
    return pruned
end

apply_tree{S}(leaf::Leaf{S}, feature::Vector) = leaf.majority

function apply_tree{S,T}(tree::Node{S,T}, features::Vector)
    if features[tree.featid] < tree.featval
        return apply_tree(tree.left, features)
    else
        return apply_tree(tree.right, features)
    end
end

function apply_tree(tree::LeafOrNode, features::Matrix)
    N = size(features,1)
    predictions = Array{Any}(N)
    for i in 1:N
        predictions[i] = apply_tree(tree, _squeeze(features[i,:],1))
    end
    if typeof(predictions[1]) <: Float64
        return float(predictions)
    else
        return predictions
    end
end

"""    apply_tree_proba(::Node, features, col_labels::Vector)

computes P(L=label|X) for each row in `features`. It returns a `N_row x
n_labels` matrix of probabilities, each row summing up to 1.

`col_labels` is a vector containing the distinct labels
(eg. ["versicolor", "virginica", "setosa"]). It specifies the column ordering
of the output matrix. """
apply_tree_proba(leaf::Leaf, features::Vector, labels) =
    compute_probabilities(labels, leaf.values)

function apply_tree_proba(tree::Node, features::Vector, labels)
    if tree.featval === nothing
        return apply_tree_proba(tree.left, features, labels)
    elseif features[tree.featid] < tree.featval
        return apply_tree_proba(tree.left, features, labels)
    else
        return apply_tree_proba(tree.right, features, labels)
    end
end

apply_tree_proba(tree::LeafOrNode, features::Matrix, labels) =
    stack_function_results(row->apply_tree_proba(tree, row, labels), features)

function build_forest(labels::Vector, features::Matrix, nsubfeatures::Integer, ntrees::Integer, partialsampling=0.7, maxdepth=-1; rng=Base.GLOBAL_RNG)
    rng = mk_rng(rng)::AbstractRNG
    partialsampling = partialsampling > 1.0 ? 1.0 : partialsampling
    Nlabels = length(labels)
    Nsamples = _int(partialsampling * Nlabels)
    forest = @parallel (vcat) for i in 1:ntrees
        inds = rand(rng, 1:Nlabels, Nsamples)
        build_tree(labels[inds], features[inds,:], nsubfeatures, maxdepth;
                   rng=rng)
    end
    return Ensemble([forest;])
end

function apply_forest{S, T}(forest::Ensemble{S, T}, features::Vector)
    ntrees = length(forest)
    votes = Array{S}(ntrees)
    for i in 1:ntrees
        votes[i] = apply_tree(forest.trees[i],features)
    end
    summarize_votes(votes)
end

summarize_votes{S<:Float64}(votes::Vector{S}) = mean(votes)
summarize_votes(votes::Vector) = majority_vote(votes)

function apply_forest(forest::Ensemble, features::Matrix)
    N = size(features,1)
    predictions = Array{Any}(N)
    for i in 1:N
        predictions[i] = apply_forest(forest, _squeeze(features[i,:],1))
    end
    if typeof(predictions[1]) <: Float64
        return float(predictions)
    else
        return predictions
    end
end

"""    apply_forest_proba(forest::Ensemble, features, col_labels::Vector)

computes P(L=label|X) for each row in `features`. It returns a `N_row x
n_labels` matrix of probabilities, each row summing up to 1.

`col_labels` is a vector containing the distinct labels
(eg. ["versicolor", "virginica", "setosa"]). It specifies the column ordering
of the output matrix. """
function apply_forest_proba(forest::Ensemble, features::Vector, labels)
    votes = [apply_tree(tree, features) for tree in forest.trees]
    return compute_probabilities(labels, votes)
end

apply_forest_proba(forest::Ensemble, features::Matrix, labels) =
    stack_function_results(row->apply_forest_proba(forest, row, labels),
                           features)

function build_adaboost_stumps{S,T}(labels::Vector{S}, features::Matrix{T}, niterations::Integer; rng=Base.GLOBAL_RNG)
    N = length(labels)
    weights = ones(N) / N
    stumps = Node{S,T}[]
    coeffs = Float64[]
    for i in 1:niterations
        new_stump = build_stump(labels, features, weights; rng=rng)
        predictions = apply_tree(new_stump, features)
        err = _weighted_error(labels, predictions, weights)
        new_coeff = 0.5 * log((1.0 + err) / (1.0 - err))
        matches = labels .== predictions
        weights[neg(matches)] *= exp(new_coeff)
        weights[matches] *= exp(-new_coeff)
        weights /= sum(weights)
        push!(coeffs, new_coeff)
        push!(stumps, new_stump)
        if err < 1e-6
            break
        end
    end
    return (Ensemble(stumps), coeffs)
end

function apply_adaboost_stumps(stumps::Ensemble, coeffs::Vector{Float64}, features::Vector)
    nstumps = length(stumps)
    counts = Dict()
    for i in 1:nstumps
        prediction = apply_tree(stumps.trees[i], features)
        counts[prediction] = get(counts, prediction, 0.0) + coeffs[i]
    end
    top_prediction = stumps.trees[1].left.majority
    top_count = -Inf
    for (k,v) in counts
        if v > top_count
            top_prediction = k
            top_count = v
        end
    end
    return top_prediction
end

function apply_adaboost_stumps(stumps::Ensemble, coeffs::Vector{Float64}, features::Matrix)
    N = size(features,1)
    predictions = Array{Any}(N)
    for i in 1:N
        predictions[i] = apply_adaboost_stumps(stumps, coeffs, _squeeze(features[i,:],1))
    end
    return predictions
end

"""    apply_adaboost_stumps_proba(stumps::Ensemble, coeffs, features, labels::Vector)

computes P(L=label|X) for each row in `features`. It returns a `N_row x
n_labels` matrix of probabilities, each row summing up to 1.

`col_labels` is a vector containing the distinct labels
(eg. ["versicolor", "virginica", "setosa"]). It specifies the column ordering
of the output matrix. """
function apply_adaboost_stumps_proba(stumps::Ensemble, coeffs::Vector{Float64},
                                     features::Vector, labels::Vector)
    votes = [apply_tree(stump, features) for stump in stumps.trees]
    compute_probabilities(labels, votes, coeffs)
end

function apply_adaboost_stumps_proba(stumps::Ensemble, coeffs::Vector{Float64},
                                    features::Matrix, labels::Vector)
    stack_function_results(row->apply_adaboost_stumps_proba(stumps, coeffs, row,
                                                           labels),
                           features)
end


text:
module LanguageServer
import URIParser
using JSON, REPL, CSTParser, DocumentFormat, SymbolServer, StaticLint
using CSTParser: EXPR, Tokenize.Tokens, Tokenize.Tokens.kind, headof, parentof, valof
using StaticLint: refof, scopeof, bindingof
using UUIDs
using Base.Docs, Markdown
import JSONRPC
using JSONRPC: Outbound, @dict_readable

export LanguageServerInstance, runserver

include("exception_types.jl")
include("uri2.jl")
include("protocol/protocol.jl")
include("extensions/extensions.jl")
include("document.jl")
include("languageserverinstance.jl")
include("runserver.jl")
include("staticlint.jl")

include("requests/misc.jl")
include("requests/textdocument.jl")
include("requests/features.jl")
include("requests/hover.jl")
include("requests/completions.jl")
include("requests/workspace.jl")
include("requests/actions.jl")
include("requests/init.jl")
include("requests/signatures.jl")
include("utilities.jl")

end

ext: jl


text:
import TOML

results = TOML.parsefile("results-final.toml")

dss = ["TwoMoons", "Sphere", "Circles", "OIP300", "OIP15"]
cols = ["vietorisrips", "witness", "mahalonobis_lmc", "mahalonobis_km"]
parts = Dict(1=>"\\% success", 2=>"median relative dominance", 4=>"median number of cells")

f = open("results-table.tex", "w")

for ds in dss
    data = results[ds]
    res = """
% dataset $ds
    \\multicolumn{1}{l}{\\textbf{$ds}} & & & & \\multicolumn{1}{l}{} \\\\ \\cline{1-5}
"""

    for part in sort(collect(keys(parts)))
        vals = fill("", size(cols)...)
        for (m,r) in data
            i = findfirst(isequal(m), cols)
            if i !== nothing
                vals[i] = string(part == 4 ? round(Int, r[part]) : round(r[part], digits=2))
            end
        end

        res *= """
% result 1
    \\multicolumn{1}{r|}{$(parts[part])} &
    \\multicolumn{1}{c}{$(vals[1])} &
    \\multicolumn{1}{c}{$(vals[2])} &
    \\multicolumn{1}{c}{$(vals[3])} &
    \\multicolumn{1}{c|}{$(vals[4])} \\\\ \\cline{1-5}
"""
    end
    println(f, res)
end

close(f)


text:
module GoogleCodeSearch

using BinaryProvider
using Sockets
using JSON
using HTTP

import Base: show
export Ctx, index, search, show, indices, clear_indices, paths_indexed, run_http

const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("GoogleCodeSearch not installed properly, run Pkg.build(\"GoogleCodeSearch\"), restart Julia and try again")
end
include(depsjl_path)
include("codesearch.jl")
include("server.jl")

end # module


text:
module CellNumbersTests

using Test
using Gridap

using ..CellValuesMocks

l = 10
a = 3.5

cn = TestIterCellValue(a,l)
@test isa(cn,CellNumber)
@test isa(cn,IterCellNumber)
r = [a for i in 1:l]
test_iter_cell_number(cn,r)
@test !isa(cn,CellArray)
@test !isa(cn,IterCellArray)
@test !isa(cn,CellMap)
@test !isa(cn,IterCellMap)

cn = TestIndexCellValue(a,l)
@test isa(cn,CellNumber)
@test isa(cn,IndexCellNumber)
test_index_cell_number(cn,r)
@test !isa(cn,CellArray)
@test !isa(cn,IndexCellArray)
@test !isa(cn,CellMap)
@test !isa(cn,IndexCellMap)

end # module


text:
# if `&` not in parenthesis: `ERROR: syntax: invalid assignment location`
(&)(hits₁::Set{String}, hits₂::Set{String}) = hits₁ ∩ hits₂
|(hits₁::Set{String}, hits₂::Set{String}) = hits₁ ∪ hits₂
-(hits₁::Set{String}, hits₂::Set{String}) = setdiff(hits₁, hits₂)
!(hits::Set{String}) = _all() - hits;


text:
"Draw a grpahic representaion of a tree using [graphivz](http://www.graphviz.org/) and its [Python
Interface](https://github.com/xflr6/graphviz)"
function drawtree end

"""
    drawtree(tree_digraph::FixedGraph, show_label=false)

Draw a graphic representation of `tree_digraph`
"""
function drawtree(tree_digraph::FixedGraph, show_label=false)
    graphviz = pyimport("graphviz")
    digraph = graphviz.Graph()

    [digraph.node("$i") for i in nodes(tree_digraph)]

    for next = edges(tree_digraph)
        parent, child = next
        digraph.edge("$parent", "$child")
    end

    digraph.graph_attr = PyDict(Dict("size" => "12,12"))

    node_attr = Dict("shape" => "circle",
                                    "color"=>"red",
                                    "fillcolor"=>"orange",
                                    "style" => "filled",
                                    "fixedsize" => "true",
                                    "nodesep" => "0.3",
                                    "ranksep" => "0.3",
                                    "width" => "0.3")
    if !show_label
        node_attr["label"] = ""
    end
    digraph.node_attr = PyDict(node_attr)

    digraph
end

"""
    drawtree(deg_seq::Vector{Int}, show_label=false)

Draw a graphic representation of tree with DFS degree sequence `deg_seq`
"""
function drawtree(deg_seq::Vector{Int}, show_label=false)
    tree_digraph = treegraph(deg_seq)
    drawtree(tree_digraph , show_label)
end

"""
    drawtree(tree::FiniteTree, show_label=false)

Draw a graphic representation of the random/fixed tree `tree`.
"""
function drawtree(tree::FiniteTree, show_label=false)
    tree_digraph = treegraph(tree)
    drawtree(tree_digraph , show_label)
end


text:
push!(LOAD_PATH, pwd())

using Plotting, Filter

function main()

	fs = 8000

	low_freqs = [697, 770, 852, 941]
	high_freqs = [1209, 1336, 1477, 1633]

	for f in low_freqs
		y = get_filter(fs, 200, f)
		plot_zplane(y, f)
	end

	for f in high_freqs
		y = get_filter(fs, 200, f)
		plot_zplane(y, f)
	end
end

if PROGRAM_FILE == "ZeroPoleAnalysis.jl"

	set_plot()

	main()
end

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
# -*- coding: utf-8 -*-
using LightGraphs
using SemanticModels
using SemanticModels.ModelTools.PetriModels
using SemanticModels.ModelTools.CategoryTheory
using ModelingToolkit
using Petri

# explicit pushout definition for testing and verification of span implementation
function pushout(a::AbstractGraph, b::AbstractGraph, f::AbstractVector{Int}, g::AbstractVector{Int})
    l = f ⊔ g
    G = a ⊔ b
    map(edges(G)) do e
        s,t = e.src, e.dst
        return Edge(l[s], l[t])
    end |> Graph
end

# # Example

n = 3
a = smallgraph(:house)
b = StarGraph(n)

f = FinSetMorph([5,3,4])
F = GraphMorph(b, a, f)
@assert verify(F) == true
@assert FinSetMorph([5,4,3,2,1])(a) |>
    adjacency_matrix |>
    collect == [0 1 1 0 0;
                1 0 1 1 0;
                1 1 0 0 1;
                0 1 0 0 1;
                0 0 1 1 0]

# +
H  = pushout(a, b, 1:5, [5, 6, 4])

# Create a decorated morphism for a
g = FinSetMorph(1:5, [5,4])
# Add graph a as decoration of morphism f
dec_g = Decorated(g, a)
# Create a decorated morphism for b
f = FinSetMorph(1:3, [1,2])
# Add graph b as decoration of morphism g
dec_f = Decorated(f, b)
# Create the span of decorated morphisms
s = Span(dec_f,dec_g)
# Solve for the decorated cospan that solves the pushout defined by the span
H′ = CategoryTheory.pushout(s)

@assert H == decoration(H′)

@show collect(edges(decoration(H′)))

# +
@variables S, I, R, I′

# Create S -> I -> R model
sir = Petri.Model([S, I, R], [(S+I, 2I), (I,R)])

# Create Susceptible to 2 diseases model
sii = Petri.Model([S, I, I′], [(S+I,  2I ), (S+I′, 2I′)])

# Decorate a finite set with SIR
f = FinSetMorph(1:3, [1, 2])
dec_f = Decorated(f, sir)

# Decorate a finite set with SII
g = FinSetMorph(1:3, [1, 2])
dec_g = Decorated(g, sii)

# Create a span of decorated morphisms
s = Span(dec_f, dec_g)

# Solve the pushout that combines the two Petri decorations
H = CategoryTheory.pushout(s)

@show decoration(H).S
@show decoration(H).Δ

#@show typeof(H)

#@show collect(edges(decoration(H)))
# -


text:
# Copyright (C) 2019, by Bill Burdick, ZLIB licensed, https://github.com/zot/domdom

using Domdom

module App

using Main.Domdom

dir = dirname(dirname(dirname(@__FILE__)))

#useverbose(true)

verbose("dir: ", dir)

abstract type Backend end

struct SimpleBackend <: Backend
    accounts
    SimpleBackend() = new(Dict())
end

mutable struct AccountMgr
    editing
    currentid
    dialogs
    mode
    backend
    editingids
    login
    function AccountMgr()
        new(false, 0, [], :none, SimpleBackend(), Set(), nothing)
    end
end

struct Dialog
    dom
    okfunc
    cancelfunc
end

mutable struct Account
    id
    username::String
    password::String
end

# App API

function accounts(::Backend) end
function getaccount(::Backend, id::Symbol) end
function addaccount(::Backend, acct::Account) end
function changeaccount(::Backend, acct::Account) end
function deleteaccount(::Backend, id) end

# Simple Backend

accounts(backend::SimpleBackend) = backend.accounts
getaccount(backend::SimpleBackend, id::Symbol) = get(backend.accounts, id, nothing)
addaccount(backend::SimpleBackend, acct::Account) = backend.accounts[acct.id] = acct
changeccount(backend::SimpleBackend, acct::Account) = ()
deleteaccount(backend::SimpleBackend, id) = delete!(backend.accounts, id)

# GUI Code

const props = domproperties

app(dom::Dom) = app(connection(dom))
app(con::Connection) = get!(con.properties, :aqua, AccountMgr())
app(mgr::AccountMgr) = mgr

accounts(dom::Union{Dom, Connection, AccountMgr}) = accounts(app(dom).backend)
getaccount(dom::Union{Dom, Connection, AccountMgr}, id) = getaccount(dom, Symbol(id))
getaccount(dom::Union{Dom, Connection, AccountMgr}, id::Symbol) = getaccount(app(dom).backend, id)
addaccount(dom::Union{Dom, Connection, AccountMgr}, acct::Account) = addaccount(app(dom).backend, acct)
changeaccount(dom::Union{Dom, Connection, AccountMgr}, acct::Account) = changeaccount(app(dom).backend, acct)
deleteaccount(dom::Union{Dom, Connection, AccountMgr}, id) = deleteaccount(app(dom).backend, id)

headerdom(item) = DomObject(item, :header, heading = "", login = app(item).login == nothing ? nothing : app(item).login.username)

logindom(item) = DomObject(item, :login, username="", password="")

refdom(item) = refdom(connection(item), webpath(item))
refdom(con, ref) = DomObject(con, :ref, path = ref)
deref(ref::DomObject) = getpath(connection(ref), web2julia(ref.path))

function accountsdom(dom)
    items = map(accountdomf(dom), sort(collect(values(accounts(dom))), by=x->x.username))
    accts = DomArray(connection(dom), [], items)
    DomObject(dom, :accounts, accounts = accts)
end

accountdomf(dom) = acct-> accountdom(dom, acct)
accountdom(dom::Dom) = accountdom(dom, dom.acctId)
accountdom(dom::Dom, acctid::Symbol) = accountdom(dom, accounts(dom)[acctid])
function accountdom(dom, acct::Account)
    obj = DomObject(dom, :account)
    copy!(obj, acct)
end

function Base.copy!(dst::Union{DomObject, Account}, src::Union{DomObject, Account})
    dst.username = src.username
    dst.password = src.password
    fromId = isa(src, DomObject) ? src.acctId : src.id
    if isa(dst, DomObject)
        dst.acctId = fromId
    else
        dst.id = fromId
    end
    dst
end

function fixaccountdoms(acctdoms)
    accts = sort(collect(values(accounts(acctdoms))), by = a-> lowercase(a.username))
    verbose("ACCOUNTS: $(accts)")
    editingids = app(acctdoms).editingids
    i = 1
    while i <= max(length(acctdoms), length(accts))
        if length(accts) < i
            pop!(acctdoms)
        else
            if length(acctdoms) < i || acctdoms[i].acctId != accts[i].id
                acct = accts[i]
                acctdoms[i] = acct.id in editingids ? editaccount(acct) : accountdom(acctdoms, acct)
            end
            i += 1
        end
    end
    acctdoms
end

editaccount(dom::Dom) = editaccount(dom, accounts(dom)[dom.acctId])
editaccount(parent::Dom, acct::Account) = DomObject(parent, :view, namespace="edit", contents = accountdom(parent, acct))

function pushDialog(dom::Dom, ok, cancel)
    push!(app(dom).dialogs, Dialog(dom, ok, cancel))
end

function message(dom::Dom)
end

function displayview(main, mode::Symbol, views...)
    con = connection(main)
    app(main).mode = mode
    main[1] = headerdom(con)
    idx = 1
    for view in views
        idx += 1
        main[idx] = view
    end
    while length(main) > length(views) + 1
        verbose("DELETE LAST FROM: ", main)
        pop!(main)
    end
    cleanall!(main)
end

function arrayitem(dom)
    while !isa(path(dom)[end], Number)
        dom = parent(dom)
    end
    (parent(dom), path(dom)[end])
end

function closeview(dom)
    length(root(dom)) > 1 && pop!(root(dom).main)
    root(dom).main[1].heading = ""
end

function home(dom)
    top = root(dom).main
    while length(top) > 1
        pop!(top)
    end
    top[1].heading = ""
end

function exampleStartFunc(backend::Backend)
    start(dir * "/html", config = (dir, host, port)-> println("STARTING HTTP ON $host:$port, DIR $dir")) do con, events
        events.onset(:login, :name) do dom, key, arg, obj, event
            verbose("SETTING USERNAME TO $arg")
            dom.currentusername = dom.username
        end
        events.onclick(:header, :login) do dom, key, arg, obj, event
            top = root(dom).main
            displayview(top, :login, logindom(dom))
            top[1].heading = "LOGIN"
        end
        events.onclick(:header, :logout) do dom, key, arg, obj, event
            app(dom).login = nothing
            root(dom).main[1].login = nothing
        end
        events.onclick(:header, :edit) do dom, key, arg, obj, event
            app(dom).editing = true
            top = root(dom).main
            editor = editaccount(dom, app(dom).login)
            push!(app(dom).editingids, dom.acctId)
            acct = editor.contents
            displayview(top, :edit, editor)
            top[1].heading = "EDITING"
            cleanall!(top)
            verbose("MODE: $(app(top).mode)")
            props(acct).save = function(dom)
                (accts, index) = arrayitem(dom)
                verbose("Save: $(accts)\nACCOUNT: $(dom)")
                backendacct = accounts(dom)[dom.acctId]
                copy!(backendacct, dom)
                changeaccount(app(dom), backendacct)
                delete!(app(dom).editingids, dom.acctId)
                if acct.username != top[1].login
                    top[1].login = acct.username
                end
                home(dom)
            end
            props(acct).cancel = function(dom)
                (accts, index) = arrayitem(dom)
                verbose("Cancel: $(accounts)")
                accts[index] = accountdom(root(dom), dom.acctId)
                delete!(app(dom).editingids, dom.acctId)
                home(dom)
            end
        end
        events.onclick(:header, :home) do dom, key, arg, obj, event
            home(dom)
        end
        events.onclick(:header, :accounts) do dom, key, arg, obj, event
            verbose("CLICKED ACCOUNTS")
            top = root(dom).main
            # show a list of refs
            # not using a ref to the list here so that each ref can be replaced with an editor
            displayview(top, :accounts, DomObject(dom, :accounts, accounts = fixaccountdoms(DomArray(dom))))
            top[1].heading = "ACCOUNTS"
        end
        events.onclick(:login, :ok) do dom, key, arg, obj, event
            if (acct = getaccount(dom, dom.username)) != nothing
                app(dom).login = acct
                root(dom).main[1].login = acct.username
                closeview(dom)
            else
                root(dom).main[end] = DomObject(dom, :message, content = "Bad user name or password")
            end
        end
        events.onclick(:login, :cancel) do dom, key, arg, obj, event
            closeview(dom)
        end
        events.onclick(:message, :ok) do dom, key, arg, obj, event
            closeview(dom)
        end
        events.onclick(:accounts, :newaccount) do dom, key, arg, obj, event
            verbose("CLICKED ACCOUNTS")
            top = root(dom).main
            acct = DomObject(dom, :account, username="", password="")
            props(acct).mode = :new
            push!(top, DomObject(dom, :newaccount, account = acct))
            props(acct).save = function(dom)
                verbose("Save NEW ACCOUNT")
                mgr = app(dom)
                (id, acct) = newaccount(mgr, acct.username, acct.username, acct.password)
                addaccount(mgr, acct)
                fixaccountdoms(top[2].accounts)
                closeview(dom)
            end
            props(acct).cancel = function(dom)
                verbose("Cancel NEW ACCOUNT")
                closeview(dom)
            end
        end
        events.onclick(:account, :edit) do dom, key, arg, obj, event
            index = path(dom)[end]
            verbose("EDIT[$(index)]: $(path(dom)), $(dom), $(key)")
            editor = root(dom).main[2].accounts[index] = editaccount(dom)
            acct = editor.contents
            push!(app(dom).editingids, dom.acctId)
            props(acct).save = function(dom)
                (accts, index) = arrayitem(dom)
                verbose("Save: $(accts)\nACCOUNT: $(dom)")
                backendacct = accounts(dom)[dom.acctId]
                copy!(backendacct, dom)
                changeaccount(app(dom), backendacct)
                delete!(app(dom).editingids, dom.acctId)
                fixaccountdoms(accts)
           end
            props(acct).cancel = function(dom)
                (accts, index) = arrayitem(dom)
                verbose("Cancel: $(accounts)")
                accts[index] = accountdom(root(dom), dom.acctId)
                delete!(app(dom).editingids, dom.acctId)
            end
        end
        events.onclick(:account, :save) do dom, key, arg, obj, event
            props(dom).save(dom)
        end
        events.onclick(:account, :cancel) do dom, key, arg, obj, event
            props(dom).cancel(dom)
        end
        events.onclick(:account, :delete) do dom, key, arg, obj, event
            verbose("Delete: $(path(dom)), key = $(key), arg = $(arg), $(dom),\nArrayItem: $(arrayitem(dom))\nID: $(dom.acctId)")
            deleteaccount(dom, dom.acctId)
            (array, index) = arrayitem(dom)
            deleteat!(array, index)
        end
        initaccounts(con, backend)
        #global mainDom = DomObject(con, :document, contents = DomObject(con, :top, main = DomArray(con, [], DomValue[headerdom(con)])))
        global mainDom = DomObject(con, :document, contents = DomObject(con, :top, main = DomArray(con, [], [headerdom(con)])))
    end
end

function newaccount(mgr, username, password)
    acct = newaccount(mgr, mgr.currentid, username, password)
    mgr.currentid += 1
    acct.id => acct
end

newaccount(mgr, id, username, password) = newaccount(mgr, Symbol(id), username, password)
function newaccount(mgr, id::Symbol, username, password)
    id => Account(id, username, password)
end

function initaccounts(con, backend)
    mgr = app(con)
    mgr.backend = backend
    empty!(accounts(mgr))
    copy!(accounts(mgr), Dict([
        newaccount(mgr, "fred", "Fred Flintstone", "fred")
        newaccount(mgr, "herman", "Herman Munster", "herman")
        newaccount(mgr, "lilly", "Lilly Munster", "lilly")
    ]))
end

if isinteractive()
    @async exampleStartFunc(SimpleBackend())
else
    exampleStartFunc(SimpleBackend())
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Giflib"
version = v"5.2.1"

# Collection of sources required to build Giflib
sources = [
    ArchiveSource("https://downloads.sourceforge.net/project/giflib/giflib-$(version).tar.gz",
                  "31da5562f44c5f15d63340a09a4fd62b48c45620cd302f77a6d9acf0077879bd"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/giflib-*/
# Apply patch to make it possible to build for non Linux-like platforms.
# Adapted from also https://sourceforge.net/p/giflib/bugs/133/
atomic_patch -p1 ../patches/makefile.patch
make -j${nproc}
make install
rm "${libdir}/libgif.a"
install_license COPYING
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libgif", :libgif),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using StatsBase
using Optim
using Distributions
using StatsPlots
#------------------------------------------------------------------------------
# time-varying IAR
using Turing
using StatsFuns
include("/home/and/code/Covid19/src/visualization.jl")

num_tot = 400
num_obs = 365
start_idx = 100
iar_step = 7

function get_iar_idx(num_tot, num_obs, start_idx)
    iar_idx = zeros(Int64, num_tot)
    iar_idx[1:start_idx-1] .= 1
    idx = 2
    for (j,i) in enumerate(start_idx:num_obs)
         iar_idx[i] = idx
         j % iar_step == 0 && (idx += 1)
    end
    if num_obs < num_tot
        iar_idx[num_obs+1:num_tot] .= iar_idx[num_obs]
    end
    iar_idx
end

iar_idx = get_iar_idx(num_tot, num_obs, start_idx)
num_iar_steps = length(unique(iar_idx))

@model function iarmodel(num_iar_steps, iar_idx, last_time_step, ::Type{TV}=Vector{Float64}) where TV
    logit_iar = TV(undef, num_iar_steps)
    σ_iar = 0.1
    iar0 ~ Beta(10,100)
    logit_iar[1] = logit(iar0)
    for i in 2:num_iar_steps
        logit_iar[i] ~ Normal(logit_iar[i-1], σ_iar)
    end
    logistic.( logit_iar[iar_idx[1:last_time_step]] )
end

m = iarmodel(num_iar_steps, iar_idx, num_tot)
d = hcat([m() for _ in 1:100]...)
plot_confidence_timeseries(d)
#---------------------------------------------------------------------------
# seroconversion
using CovidSurvey

pyplot()
nbins = 300
function seropositive(nbins)
    N = 10_000_000
    # For details see Imperial Report34
    rbins = 1.5:1:(nbins + .5)
    inf2symp = rand(GammaMeanCv(5.1, .86), N)
    inf2conv = inf2symp .+ rand(GammaMeanCv(12.6,5.8/12.6), N)
    inf2rev = inf2conv .+ rand(Weibull(3.67, 132.88), N)
    seroconverted = map( ecdf(inf2conv), rbins)
    notseroreverted = 1 .- map( ecdf(inf2rev), rbins)
    serconverted .* notseroreverted
end
inf2conv = rand(GammaMeanCv(12.6,5.8/12.6), N)
seroconverted = map( ecdf(inf2conv), rbins)
inf2symp = rand(GammaMeanCv(5.1, .86), N)
inf2conv = inf2symp .+ rand(Exponential(13.3), N)
plot(1.5:1:(nbins + .5), 1 .- seroconverted)
seroconverted[20]
seropos = seropositive(nbins)
seropos[20]
plot(1.5:1:(nbins + .5), seropositive(nbins))

using QuadGK
bins = 200
λ = 13.3
μ = 143.7
κ = 3.67
f(x, ts, κ, μ, λ) = 1/λ*exp(-x/λ-((ts-x)/μ)^κ)
π = zeros(bins)
for t in 1:200
    integral, err = quadgk(x->f(x, t, κ, μ, λ), 0, t, rtol=1e-8)
    π[t] = integral
end
plot(π)
#---------------------------------------------------------------------------
# estimate parameters
function gamma(μ, cv)
    α = 1/cv^2
    θ = μ/α
    Gamma(α, θ)
end

function loss(θ)
    cv = θ
    d = gamma(2.6, cv)
    x = rand(d, 10_000)
    lq, uq = percentile(x, [2.5, 97.5])
    m = mean(x)
    sum(abs2, [lq-0.4, uq-9] )
end

op = optimize(θ->loss(θ), 0.1, 3)
cv = op.minimizer
percentile(rand(gamma(2.6, .92), 10_000), [2.5, 97.5])
##
function loss(θ, M, LQ, UP)
    m, cv = exp.(θ)
    d = gamma(m, cv)
    # lq, uq = percentile(d, [2.5, 97.5])
    lq, uq = percentile(d, [2.5, 97.5])
    sum(abs2, [m-M, lq-LQ, uq-UQ] )
end

data.seroprev[!,:m] .= 0.
data.seroprev[!,:cv] .= 0.
for v in eachrow(data.seroprev)
    M = v.country
    LQ, UQ = v.CI
    op = optimize(θ->loss(θ, M, LQ, UQ), log.([M, (UQ-M)/1.96]), NelderMead())
    v.m, v.cv = exp.(op.minimizer)
end

let
    t = data.seroprev.date
    CIs = data.seroprev.CI
    Ms  = data.seroprev.country
    le = [m-ci[1] for (m,ci) in zip(Ms, CIs)]
    ue = [ci[2]-m for (m,ci) in zip(Ms, CIs)]
    scatter(t, Ms, yerror=(le,ue))
    cvs = data.seroprev.cv
    ms  = data.seroprev.m
    le = [m-percentile(gamma(m,cv),2.5) for (m,cv) in zip(ms, cvs)]
    ue = [percentile(gamma(m,cv),97.5)-m for (m,cv) in zip(ms, cvs)]
    scatter!(t.+Day(1), ms, yerror=(le,ue))
end
##
function loss(θ, M, LQ)
    m, σ = exp.(θ)
    d = Normal(m, cv)
    # lq, uq = percentile(d, [2.5, 97.5])
    uq = percentile(d, 2.5)
    sum(abs2, [m-M, lq-LQ] )
end

data.seroprev[!,:m] .= 0.
data.seroprev[!,:std] .= 0.
for v in eachrow(data.seroprev)
    M = v.country
    LQ, UQ = v.CI
    op = optimize(θ->loss(θ, M, LQ), log.([M, (M-LQ)/1.96]), NelderMead())
    v.m, v.std = exp.(op.minimizer)
end

let
    t = data.seroprev.date
    CIs = data.seroprev.CI
    Ms  = data.seroprev.country
    le = [m-ci[1] for (m,ci) in zip(Ms, CIs)]
    ue = [ci[2]-m for (m,ci) in zip(Ms, CIs)]
    scatter(t, Ms, yerror=(le,ue))
    stds = data.seroprev.std
    ms  = data.seroprev.m
    le = [m-percentile(Normal(m,std),2.5) for (m,std) in zip(ms, stds)]
    ue = [percentile(Normal(m,std),97.5)-m for (m,std) in zip(ms, stds)]
    scatter!(t.+Day(1), ms, yerror=(le,ue))
    ms  = data.seroprev.m
    le = [m-percentile(Normal(m,(m-ci[1])/1.96),2.5) for (m,ci) in zip(ms, CIs)]
    ue = [percentile(Normal(m,(m-ci[1])/1.96),97.5)-m for (m,ci) in zip(ms, CIs)]
    scatter!(t.+Day(2), ms, yerror=(le,ue))
end

#----------------------------------------------------------------------------
# calculate kernel function
N = 10_000_000
x = rand(gamma(2.6, cv), N)
h = fit(Histogram, x, 1.5:1:30.5, closed=:right)
w = h.weights ./ N
bar(1.5:1:30.5, w)

#----------------------------------------------------------------------------
# calculate kernel function
# infection to death
N = 10_000_000
x = rand(gamma(5.1, .86), N) .+ rand(gamma(17.8, .45), N)
h = fit(Histogram, x, vcat([0], 1.5:1:100.5))

data_nature = ImperialReport13.load_data(datadir("imperial-report13", "processed.rds"))
turing_data = data_nature.turing_data
π = turing_data.π[1][1:100]
w = h.weights / N * 0.0102
# pushfirst!(w, 0.)
plotlyjs()
bar([w, π], fα=.2, lc=:match)
##
data_usa    = ImperialUSAcases.load_data(datadir("imperial-usa-cases", "processed.rds"))
π2 = data_usa.turing_data.π2[51][1:100]
h = fit(Histogram, rand(gamma(5.1, .86), N), vcat([0], 1.5:1:100.5))
w = h.weights ./ N
bar([w, π2], fα=.2, lc=:match)
#-----------------------------------------------------------------------------
# hospitalization rate
using CSV
using DataFrames

pop = begin
    pop = CSV.File("/home/and/data/population_DK.csv", delim=";", header=false) |> DataFrame
    pop = pop[:,end-1:end]
    rename!(pop, Dict(:Column4 => :age, :Column5 => :size))
    parseyear(s) = begin
        y, _ = split(s," ")
        parse(Int64, y)
    end
    pop.age = parseyear.(pop.age)
    pop
end

# Verity et al. Lancet (2020)
# 10 year bands: 0-9, 10-19, ..., ≥ 80
hosp = begin
    percent = [0.00
            0.0408
            1.04
            3.43
            4.25
            8.16
            11.8
            16.6
            18.4
    ]
    age = ["$(10*i)-$(10*(i+1)-1) years" for i in 0:8]
    age[end] = "≥ 80 year"
    DataFrame(Dict(:percent => percent, :age => age))
end

function addpop!(hosp, pop)
    hosp[!,:pop] .= 0
    for i in 0:7
        for j in 1:10
            hosp[i+1, :pop] += pop[i*10+j, :size]
        end
    end
    hosp[end, :pop] = sum(pop[81:end,:size])
    nothing
end
addpop!(hosp, pop)

ihr = hosp.percent'hosp.pop / sum(hosp.pop)
#ihr = 6.09
#----------------------------------------------------------------------------
# infection fatality ration

ifr = 100 / 100_000 # 95%CI (72–211) / 100_000
ifr*(rand(Normal(1,.3),1000) .+ rand(Normal(.2,.1),1000)) |> histogram
ihr = 0.3

rand(truncated(Normal(.25/100,.1/100),.01/100,1/100), 1000) .* 100 |> histogram
rand(Beta(8,20), 1000) |> histogram

rand(truncated(Normal(5, 5), 1e-1, Inf), 1000) |>histogram
rand(Gamma(3,1), 1000) |>histogram

rand(Beta(8,20), 1000).*.6|>histogram
rand(Beta(4,15), 1000).*.6|>histogram
rand(Beta(5,1), 1000) |>histogram

rand(truncated(Normal(0, 5), 1e-1, Inf),1000) |> histogram

function NegativeBinomial3(μ, σ2)
    p = μ/σ2
    r = μ^2 / (σ2 - μ)
    NegativeBinomial(r, p)
end
x = rand(NegativeBinomial3(30, 40),10000)
mean(x)
var(x)
std(x) / mean(x)
x|>histogram

x = rand(NegativeBinomial2(30, ((40-30)/30^3)^(-1)),10000)
mean(x)
var(x)
std(x) / mean(x)

rand(NegativeBinomial3(30, 30+30^2/rand(truncated(Normal(0, 5), 1e-1, Inf)))) for i in 1:1000]|>histogram
rand(NegativeBinomial2(30, ),1000)|>histogram
##
ts = data.dates
ds = data.cases
m  = length(ds)
n  = 1000
ys = zeros(m,n)
for j in 1:n
    for i in 1:m
        d = ds[i] .+ 1e-2
        ϕ = rand(truncated(Normal(0, 2), 0, Inf))
        ys[i,j] = rand(NegativeBinomial3(d, d+d^2/ϕ))
        # ys[i,j] = rand(NegativeBinomial2(d, ϕ))
    end
end
p = plot(ts,ds)
plot_confidence_timeseries!(p, ts, ys)
##

p0 = histogram(rand(Normal(turing_data.seroprev_mean[end], turing_data.seroprev_std[end]), 10000))
σs = rand(InverseGamma2(turing_data.seroprev_std[end], .1), 10000)
p1 = histogram( [rand(Normal(turing_data.seroprev_mean[end], σs[i])) for i in 1:10000] )
μs = rand(Normal(turing_data.seroprev_mean[end], 10000), 10000)
p2 = histogram( [rand(Normal(μs[i], σs[i])) for i in 1:10000] )
plot(p0,p1,p2,layout=(3,1), sharex=true)
# μs|>histogram

##
ts = epidemic_start:num_obs
# μs  = expected_daily_deaths[ts]
# σ2s = @. μs + μs^2 / ϕ
# deaths[ts] ~ arraydist(NegativeBinomial3.(μs, σ2s))
#
#---------------------------------------------------------------------------
n = 365
m = 100
rw = randn(n,m)*.4
rw[1,:] .= 0
rw = cumsum(rw, dims=1)
# plot_confidence_timeseries(rw)
plot(rw[:,1:10])
#--------------------------------------------------------------------------
pyplot()
0.6 * (1 .+ randn(1000)*.2) |> histogram
rand(Beta(15,35),1000)|>histogram
[rand(Exponential( rand(Exponential(1 / 0.03)) )) for i in 1:1000] |> histogram
rand(LogNormal(5, .5), 10000)|>histogram
##
# si = [ 1.832618e-02 6.659231e-02 1.019139e-01 1.177168e-01 1.183856e-01 1.096347e-01 9.612322e-02 8.105485e-02 6.638313e-02 5.314896e-02 4.178950e-02 3.237506e-02 2.477423e-02 1.876128e-02 1.408140e-02 1.048743e-02 7.758059e-03 5.704855e-03 4.172845e-03 3.037799e-03 2.202072e-03 1.590108e-03 1.144187e-03 8.206829e-04 5.869190e-04 4.186079e-04 2.978198e-04 2.113961e-04 1.497302e-04 1.058413e-04 7.467784e-05 5.259832e-05 3.698645e-05 2.596852e-05 1.820643e-05 1.274709e-05 8.913312e-06 6.224997e-06 4.342490e-06 3.025966e-06 2.106388e-06 1.464818e-06 1.017703e-06 7.064292e-07 4.899421e-07 3.395203e-07 2.350969e-07 1.626684e-07 1.124728e-07 7.771288e-08 5.366011e-08 3.702835e-08 2.553596e-08 1.760009e-08 1.212361e-08 8.346653e-09 5.743340e-09 3.949995e-09 2.715288e-09 1.865648e-09 1.281283e-09 8.795663e-10 6.035400e-10 4.139650e-10 2.838226e-10 1.945186e-10 1.332637e-10 9.126477e-11 6.248002e-11 4.275902e-11 2.925304e-11 2.000644e-11 1.367839e-11 9.348966e-12 6.388001e-12 4.363510e-12 2.979728e-12 2.034262e-12 1.388445e-12 9.474643e-13 6.462608e-13 4.407585e-13 3.004264e-13 2.049472e-13 1.395550e-13 9.514611e-14 6.472600e-14 4.418688e-14 3.008704e-14 2.042810e-14 1.387779e-14 9.547918e-15 6.439294e-15 4.329870e-15 2.997602e-15 1.998401e-15 1.443290e-15 8.881784e-16 6.661338e-16 4.440892e-16 ] |> vec
# plot(si)
# plot!(turing_data.serial_intervals)
# xlims!(0,100)

#10.92
##
μs = rand(Normal(10.92, 0.94), 100)
ϕs = rand(Normal(5.4, 0.27), 100)
i2c = NegativeBinomial2.(μs, ϕs)
let
    x = 1:50
    p = plot()
    for i in 1:length(i2c)
        plot!(x, i2c[i], c=:black, α=0.1, label=nothing, ls=:solid)
    end
    plot!(x, NegativeBinomial2.(10.92, 5.4), c=:red, lw=2, ls=:solid)
    plot!(x, turing_data.π2[x], c=:green, lw=2)
end
##
@model function i2c_sampler(::Type{V}  = Int64) where V
    μ_i2o ~ truncated(Normal(1.53, 0.051), 0, Inf)
    σ_i2o ~ truncated(Normal(0.48, 0.0759), 0, Inf)
    i2o   ~ LogNormal(μ_i2o, σ_i2o)

    μ_o2c ~ truncated(Normal(5.82, 0.68), 0, Inf)
    ϕ_o2c ~ truncated(Normal(1.57, 0.054), 0, Inf)
    o2c   ~ NegativeBinomial2(μ_o2c, ϕ_o2c)

    return V(round(i2o+o2c))
end

@model function i2h_sampler(::Type{V}  = Int64) where V
    μ_i2o ~ truncated(Normal(1.53, 0.051), 0, Inf)
    σ_i2o ~ truncated(Normal(0.48, 0.0759), 0, Inf)
    i2o   ~ LogNormal(μ_i2o, σ_i2o)

    μ_o2c ~ truncated(Normal(5.82, 0.68), 0, Inf)
    ϕ_o2c ~ truncated(Normal(1.57, 0.054), 0, Inf)
    o2c   ~ NegativeBinomial2(μ_o2c, ϕ_o2c)

    return V(round(i2o+o2c))
end
## infection 2 symptom onset
plot(0:.1:15, GammaMeanCv(5.1, .86))
plot!(0:.1:15, LogNormal(1.53, 0.48))
## serial interval
# plot(National.serialinterval(15))
plot!(0:.1:15, GammaMeanCv(5.06, 2.11/5.06))
plot!(0:.1:15, GammaMeanCv(6.5, .62))
## infection 2 death
plot(National.inf2death(90))
plot!(1:90, NegativeBinomial2(21.82, 14.26))
##
m = i2c_sampler()
i2c = [m() for _ in 1:10_000]

@model function fitNB(y)
    μ ~ LogNormal(2., .7)
    ϕ ~ LogNormal(0.5, 1.2)
    y .~ NegativeBinomial2(μ, ϕ)
end

c = sample(fitNB(i2c), NUTS(1000, 0.6), 2000, progress=true)
c = c[1001:end,:,:]

plot(c)
c["μ"]|>vec|>mean_and_std
c["ϕ"]|>vec|>mean_and_std
histogram( i2c, label=nothing, normalize=true)
plot!(NegativeBinomial2(c["μ"]|>vec|>mean, c["ϕ"]|>vec|>mean))
## infected 2 death
# from Flaxman et al.
i2d = [rand(GammaMeanCv(5.1,0.86))+rand(GammaMeanCv(17.8,0.45)) for _ in 1:10_000]
#6.49 5 .94-7 .10
@model function fitNB(y)
    μ ~ LogNormal(4.,.5)
    # ϕ ~ LogNormal(1.,1.2)
    cv ~ truncated(Normal(.5, .3), 0, Inf)
    y .~ GammaMeanCv(μ, cv)
    return μ, cv
end
m = fitNB(i2d)
let
    p = histogram( i2d, label=nothing, normalize=true)
    for i in 1:10
        plot!(1:90, GammaMeanCv(m()...))
    end
    p
end
c = sample(m, NUTS(500, 0.6), 2000, progress=true)
c = c[501:end,:,:]

plot(c)
c["μ"]|>vec|>mean_and_std # = 22.8
c["cv"]|>vec|>mean_and_std # = 0.4
histogram( i2d, label=nothing, normalize=true)
plot!(GammaMeanCv(c["μ"]|>vec|>mean, c["cv"]|>vec|>mean))
##
function inf2death(μ, nbins=90)
	r = zeros(nbins+1)
	r[2:end] = 1.5:1:(nbins + .5)
    return diff( cdf.(GammaMeanCv(μ, 0.4), r) )
end

p = let
    p = plot()
    for μ in 20.:2.:30.
        plot!(1:90, GammaMeanCv(μ, 0.4))
        plot!(1:90, inf2death(μ, 90), ls=:dash)
    end
    p
end

#----------------------------------------------------------------------------
# test capacity
function process_ssi(df)
    rename!(df, 1=>:region, 2=>:date, 3=>:counts)
    df.region = replace(df.region, dan2eng...)
    dropmissing!(df, :region)
    df = unstack(df, :date, :region, :counts)
    insert_missing_dates!(df)
    replace_missing!(df)
    DataFrames.transform(df, AsTable(Not(:date)) => sum => :country)
end

df = CSV.File(normpath( homedir(), "data/Statens-Serum-Institut/dashboard/Regionalt_DB/16_pcr_og_antigen_test_pr_region.csv" ) ) |>
    DataFrame #|> National.process_ssi

rename!(df, 1=>:date, 2=>:region, 3=>:tested, 4=>:method)
df.region = replace(df.region, National.dan2eng...)
df.region = replace(df.region, missing=>"unkown")

pcr = df[df.method .== "PCR",:]
pcr = unstack(pcr, :date, :region, :tested)
# National.insert_missing_dates!(pcr)
National.replace_missing!(pcr)
pcr = DataFrames.transform(pcr, AsTable(Not(:date)) => sum => :country)

df = unstack(df, :week, :region, :tested)
DataFrames.transform(df, AsTable(Not(:date)) => sum => :country)
insert_missing_dates!(df)
replace_missing!(df)
DataFrames.transform(df, AsTable(Not(:date)) => sum => :country)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
abstract type AbstractLocalTuple end

struct LocalTuple{T<:Tuple} <: AbstractLocalTuple
    tup::T
end

function Base.getindex(locals::LocalTuple, I::AbstractArray{<:Integer})
    arrtype = typeof(I)
    @show arrtype
    outarr = [Base.getindex(locals, i) for i in I]
    reshape(outarr, size(I))
end
function Base.getindex(locals::LocalTuple, i::Integer)
    if i > 2
        D = Differential(locals.tup[1])
        Num.(ModelingToolkit.diff2term.(ModelingToolkit.value.((D^(i - 2)).(locals.tup[2]))))
    else
        locals.tup[i]
    end
end

time(locals::LocalTuple) = locals[1]
position(locals::LocalTuple) = locals[2]
velocity(locals::LocalTuple) = locals[3]
acceleration(locals::LocalTuple) = locals[4]

function generate_generic_localtuple(qsize)
    @parameters t
    qrange = Tuple(1:qsize[i] for i in 1:length(qsize))
    @variables q[qrange...](t)
    LocalTuple((t, q))
end

function find_multiplies_to_derivative(expression)
	#@show expression
	#@show typeof(expression.val)
	#@show ModelingToolkit.isdifferential(expression.val)
	diffs = []
	premults = []
	if typeof(expression) <: SymbolicUtils.Mul
		#@show "mul!!"
		coeff = expression.coeff
		push!(premults, coeff)
		for key in keys(expression.dict)
			if ModelingToolkit.isdifferential(key)
				#@show "isa diff"
				push!(diffs, key)
			else
				#@show "isnota diff"
				push!(premults, key)
			end
			#@show key
		end
		#recursed_mul_premult, recursed_mul_deriv = recursively_find_multiplies_to_derivative(
	elseif ModelingToolkit.isdifferential(expression)
		push!(diffs, expression)
	end
	diffs, premults
end

function Lagrange_equations(Lagrangian, locals::AbstractLocalTuple)
    symbolic_L = Lagrangian(locals)
    t = time(locals)
    D = Differential(t)
    qvars = position(locals)
    qdotvars = velocity(locals)
    velocity_equations = Equation[]
    position_equations = Equation[]
	for (qvar, qdotvar) in zip(qvars, qdotvars)
		lhs_premult = MTderivative(MTderivative(symbolic_L, qdotvar), t)
		lhs_diffs, lhs_mults = find_multiplies_to_derivative(lhs_premult.val)
		rhs = MTderivative(symbolic_L, qvar)
		lhs = lhs_premult
		if length(lhs_mults) > 0
			lhs_multiplied = prod(Num.(lhs_mults))
			rhs = rhs / lhs_multiplied
			lhs = Num(lhs_diffs[1])
		end
		push!(velocity_equations, lhs ~ rhs)
		push!(position_equations, D(qvar) ~ qdotvar)
	end
	ArrayPartition(velocity_equations, position_equations)
end

rdot(v1, v2) = sum(v1 .* v2)

MTderivative(O, v) = ModelingToolkit.derivative(O, v; simplify=true)

rhs(eqs) = map(eq->eq.rhs, eqs)

function ode_problem_generator(Lagrangian_pre_param, locals::AbstractLocalTuple, params)
    Lagrangian = Lagrangian_pre_param(params...)
    qdim = length(SICMModelingToolkit.position(locals))
    current_Lagrange_equations = Lagrange_equations(Lagrangian, locals)

    f_1_rhs = Num.(rhs(current_Lagrange_equations.x[1]))
    f_2_rhs = Num.(rhs(current_Lagrange_equations.x[2]))

    qs = SICMModelingToolkit.position(locals)
    vs = velocity(locals)
    #dvs = [vs; qs]
    ps = params
    iv = SICMModelingToolkit.time(locals)

    f_1_built = build_function(f_1_rhs, vs, qs, ps, iv; expression=Val{false})[2]
    f_2_built = build_function(f_2_rhs, vs, qs, ps, iv; expression=Val{false})[2]

    function generator(qv0, q0, tspan, p0; saveat=nothing)
        if isnothing(saveat)
            prob = DynamicalODEProblem{true}(f_1_built, f_2_built, qv0, q0, tspan, p0)
        else
            prob = DynamicalODEProblem{true}(f_1_built, f_2_built, qv0, q0, tspan, p0; saveat=saveat)
        end
    end
    generator
end

max_stars_repo_name: zoemcc/SICMModelingToolkit

max_issues_repo_name: zoemcc/SICMModelingToolkit

max_forks_repo_name: zoemcc/SICMModelingToolkit


text:
"""
    FunctionCurveSpace(m, grid)

Space of curves in manifold `m`. When needed for calculations, uses grid `grid`
but does not discretize the curve when not necessary.
They have no embedding in a finite-dimensional euclidean space
(use `DiscretizedCurves` if you need this feature).
"""
struct FunctionCurveSpace{F,TM<:AbstractManifold{F},R<:AbstractRange} <: AbstractManifold{F}
    manifold::TM
    approx_grid::R
end

function FunctionCurveSpace(M::TM) where {F,TM<:AbstractManifold{F}}
    approx_grid = 0:0.01:1.0
    return FunctionCurveSpace{F,TM,typeof(approx_grid)}(M, approx_grid)
end

"""
    VectorizedFunction(f)

A Function with vector space operations defined.
"""
struct VectorizedFunction{F}
    f::F
end

*(a::Number, f::VectorizedFunction) = VectorizedFunction(t -> a * f.f(t))
function +(f1::VectorizedFunction, f2::VectorizedFunction)
    return VectorizedFunction(t -> f1.f(t) + f2.f(t))
end
function -(f1::VectorizedFunction, f2::VectorizedFunction)
    return VectorizedFunction(t -> f1.f(t) - f2.f(t))
end
-(f::VectorizedFunction) = VectorizedFunction(t -> -f.f(t))

(f::VectorizedFunction)(t) = f.f(t)

function distance(M::FunctionCurveSpace, p1, p2)
    # Calculates simple L2-like distance
    # There are also other reasonable Riemannian structures for this space
    reltol, abstol =
        concretize_tols(M, p1, p2, reltol=PARAMS.quad_rel_tol, abstol=PARAMS.quad_abs_tol)

    return sqrt(
        quadgk(
            t -> distance(M.manifold, p1(t), p2(t))^2,
            0.0,
            1.0,
            rtol=reltol,
            atol=abstol,
        )[1],
    )
end

function exp(M::FunctionCurveSpace, p, X)
    return t -> exp(M.manifold, p(t), X(t))
end

function geodesic(M::FunctionCurveSpace, p, X)
    return t -> (s -> geodesic(M.manifold, p(t), X(t), s))
end

function injectivity_radius(M::FunctionCurveSpace)
    # TODO: check this
    return injectivity_radius(M.manifold)
end

function inner(M::FunctionCurveSpace, p, X1, X2)
    reltol, abstol =
        concretize_tols(M, X1, X2, reltol=PARAMS.quad_rel_tol, abstol=PARAMS.quad_abs_tol)
    I, err = QuadGK.quadgk(
        s -> inner(M.manifold, p(s), X1(s), X2(s)),
        0.0,
        1.0,
        rtol=reltol,
        atol=abstol,
    )
    return I
end

function inverse_retract(
    M::FunctionCurveSpace,
    p,
    q,
    method::AbstractInverseRetractionMethod,
)
    return VectorizedFunction(t -> inverse_retract(M.manifold, p(t), q(t), method))
end

function inverse_retract(M::FunctionCurveSpace, p, q)
    return VectorizedFunction(t -> inverse_retract(M.manifold, p(t), q(t)))
end

function isapprox(
    M::FunctionCurveSpace,
    p1,
    p2;
    atol=atoldefault(M, p1, p2),
    rtol=rtoldefault(M, p1, p2),
)
    #=for i in M.approx_grid
        if !(isapprox(M.manifold, p1(i), p2(i), atol = atol, rtol = rtol))
            println("*** $i, $(p1(i)), $(p2(i)).")
        end
    end=#
    return all(
        isapprox(M.manifold, p1(i), p2(i), atol=atol, rtol=rtol) for i in M.approx_grid
    )
end

function isapprox(
    M::FunctionCurveSpace,
    p,
    X1,
    X2;
    atol=atoldefault(M, X1, X2),
    rtol=rtoldefault(M, X1, X2),
)
    #TODO add tolerance parameters to isapprox for this quadrature?
    I, err = quadgk(t -> norm(M.manifold, p(t), X1(t) - X2(t)), 0.0, 1.0)
    return I <= 2 * err + atol && I >= -2 * err - atol
end

function log(M::FunctionCurveSpace, p, q)
    return VectorizedFunction(t -> log(M.manifold, p(t), q(t)))
end

function manifold_dimension(M::FunctionCurveSpace)
    return Inf
end

function mid_point(M::FunctionCurveSpace, p, q)
    return t -> mid_point(M.manifold, p(t), q(t))
end

function retract(M::FunctionCurveSpace, p, X)
    return VectorizedFunction(t -> retract(M.manifold, p(t), X(t)))
end

function retract(M::FunctionCurveSpace, p, X, method::AbstractRetractionMethod)
    return VectorizedFunction(t -> retract(M.manifold, p(t), X(t), method))
end

function transport_srvf(M::FunctionCurveSpace, c_p, c_X, p)
    return t -> vector_transport_to(M.manifold, c_p(t), c_X(t), p)
end

function vector_transport_direction(
    M::FunctionCurveSpace,
    p,
    X,
    d,
    method::AbstractVectorTransportMethod,
)
    return VectorizedFunction(
        t -> vector_transport_direction(M.manifold, p(t), X(t), d(t), method),
    )
end

function zero_vector(M::FunctionCurveSpace, p)
    return VectorizedFunction(t -> zero_vector(M.manifold, p(t)))
end

function zero_vector(M::Manifolds.TangentBundleFibers{<:FunctionCurveSpace}, p)
    return VectorizedFunction(t -> zero_vector(M.manifold.manifold, p(t)))
end


text:
module Interpolations

export
    interpolate,
    interpolate!,
    extrapolate,
    scale,
    bounds,

    AbstractInterpolation,
    AbstractExtrapolation,

    OnCell,
    OnGrid,

    Flat,
    Line,
    Free,
    Periodic,
    Reflect,
    Natural,
    InPlace,
    InPlaceQ,
    Throw,

    LinearInterpolation,
    CubicSplineInterpolation

    # see the following files for further exports:
    # b-splines/b-splines.jl
    # extrapolation/extrapolation.jl
    # monotonic/monotonic.jl
    # scaling/scaling.jl

using LinearAlgebra, SparseArrays
using StaticArrays, WoodburyMatrices, Ratios, AxisAlgorithms, OffsetArrays

using Base: @propagate_inbounds
import Base: convert, size, axes, promote_rule, ndims, eltype, checkbounds, axes1

abstract type Flag end
abstract type InterpolationType <: Flag end
struct NoInterp <: InterpolationType end
abstract type GridType <: Flag end
struct OnGrid <: GridType end
struct OnCell <: GridType end

const DimSpec{T} = Union{T,Tuple{Vararg{Union{T,NoInterp}}},NoInterp}

abstract type AbstractInterpolation{T,N,IT<:DimSpec{InterpolationType}} <: AbstractArray{T,N} end
abstract type AbstractInterpolationWrapper{T,N,ITPT,IT} <: AbstractInterpolation{T,N,IT} end
abstract type AbstractExtrapolation{T,N,ITPT,IT} <: AbstractInterpolationWrapper{T,N,ITPT,IT} end

function Base.:(==)(itp1::AbstractInterpolation, itp2::AbstractInterpolation)
    propertynames(itp1) == propertynames(itp2) || return false
    for pname in propertynames(itp1)
        getproperty(itp1, pname) == getproperty(itp2, pname) || return false
    end
    return true
end

"""
    BoundaryCondition

An abstract type with one of the following values (see the help for each for details):

- `Throw(gt)`
- `Flat(gt)`
- `Line(gt)`
- `Free(gt)`
- `Periodic(gt)`
- `Reflect(gt)`
- `InPlace(gt)`
- `InPlaceQ(gt)`

where `gt` is the grid type, e.g., `OnGrid()` or `OnCell()`. `OnGrid` means that the boundary
condition "activates" at the first and/or last integer location within the interpolation region,
`OnCell` means the interpolation extends a half-integer beyond the edge before
activating the boundary condition.
"""
abstract type BoundaryCondition <: Flag end
# Put the gridtype into the boundary condition, since that's all it affects (see issue #228)
# Nothing is used for extrapolation
struct Throw{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
struct Flat{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
struct Line{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
struct Free{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
struct Periodic{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
struct Reflect{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
struct InPlace{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
# InPlaceQ is exact for an underlying quadratic. This is nice for ground-truth testing
# of in-place (unpadded) interpolation.
struct InPlaceQ{GT<:Union{GridType,Nothing}} <: BoundaryCondition gt::GT end
const Natural = Line

(::Type{BC})() where BC<:BoundaryCondition = BC(nothing)
function Base.show(io::IO, bc::BoundaryCondition)
    print(io, nameof(typeof(bc)), '(')
    bc.gt === nothing || show(io, bc.gt)
    print(io, ')')
end

Base.IndexStyle(::Type{<:AbstractInterpolation}) = IndexCartesian()

size(exp::AbstractExtrapolation) = size(exp.itp)
axes(exp::AbstractExtrapolation) = axes(exp.itp)

twotuple(r::AbstractUnitRange) = (first(r), last(r))
twotuple(x, y) = (x, y)

"""
    bounds(itp::AbstractInterpolation)

Return the `bounds` of the domain of `itp` as a tuple of `(min, max)` pairs for each coordinate. This is best explained by example:

```jldoctest
julia> itp = interpolate([1 2 3; 4 5 6], BSpline(Linear()));

julia> bounds(itp)
((1, 2), (1, 3))

julia> data = 1:3;

julia> knots = ([10, 11, 13.5],);

julia> itp = interpolate(knots, data, Gridded(Linear()));

julia> bounds(itp)
((10.0, 13.5),)
```
"""
bounds(itp::AbstractInterpolation) = map(twotuple, lbounds(itp), ubounds(itp))
bounds(itp::AbstractInterpolation, d) = bounds(itp)[d]

itptype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = IT
itptype(::Type{ITP}) where {ITP<:AbstractInterpolation} = itptype(supertype(ITP))
itptype(itp::AbstractInterpolation) = itptype(typeof(itp))
gridtype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = GT
gridtype(::Type{ITP}) where {ITP<:AbstractInterpolation} = gridtype(supertype(ITP))
gridtype(itp::AbstractInterpolation) = gridtype(typeof(itp))
ndims(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = N
ndims(::Type{ITP}) where {ITP<:AbstractInterpolation} = ndims(supertype(ITP))
ndims(itp::AbstractInterpolation) = ndims(typeof(itp))
eltype(::Type{AbstractInterpolation{T,N,IT}}) where {T,N,IT<:DimSpec{InterpolationType}} = T
eltype(::Type{ITP}) where {ITP<:AbstractInterpolation} = eltype(supertype(ITP))
eltype(itp::AbstractInterpolation) = eltype(typeof(itp))

"""
    n = count_interp_dims(ITP)

Count the number of dimensions along which type `ITP` is interpolating.
`NoInterp` dimensions do not contribute to the sum.
"""
count_interp_dims(::Type{ITP}) where {ITP<:AbstractInterpolation} = count_interp_dims(itptype(ITP), ndims(ITP))
count_interp_dims(::Type{IT}, n) where {IT<:InterpolationType} = n * count_interp_dims(IT)
count_interp_dims(it::Type{IT}, n) where IT<:Tuple{Vararg{InterpolationType,N}} where N =
    _count_interp_dims(0, it...)
@inline _count_interp_dims(c, ::IT1, args...) where IT1 =
    _count_interp_dims(c + count_interp_dims(IT1), args...)
_count_interp_dims(c) = c

"""
    BoundsCheckStyle(itp)

A trait to determine dispatch of bounds-checking for `itp`.
Can return `NeedsCheck()`, in which case bounds-checking is performed, or `CheckWillPass()`
in which case the check will return `true`.
"""
abstract type BoundsCheckStyle end
struct NeedsCheck <: BoundsCheckStyle end
struct CheckWillPass <: BoundsCheckStyle end

BoundsCheckStyle(itp) = NeedsCheck()

"""
    wi = WeightedIndex(indexes, weights)

Construct a weighted index `wi`, which can be thought of as a generalization of an
ordinary array index to the context of interpolation.
For an ordinary vector `a`, `a[i]` extracts the element at index `i`.
When interpolating, one is typically interested in a range of indexes and the output is
some weighted combination of array values at these indexes.
For example, for linear interpolation at a location `x` between integers `i` and `i+1`, we have

    ret = (1-f)*a[i] + f*a[i+1]

where `f = x-i` lies between 0 and 1. This can be represented as `a[wi]`, where

    wi = WeightedIndex(i:i+1, (1-f, f))

i.e.,

    ret = sum(a[indexes] .* weights)

Linear interpolation thus constructs weighted indices using a 2-tuple for `weights` and
a length-2 `indexes` range.
Higher-order interpolation would involve more positions and weights (e.g., 3-tuples for
quadratic interpolation, 4-tuples for cubic).

In multiple dimensions, separable interpolation schemes are implemented in terms
of multiple weighted indices, accessing `A[wi1, wi2, ...]` where each `wi` is the
`WeightedIndex` along the corresponding dimension.

For value interpolation, `weights` will typically sum to 1.
However, for gradient and Hessian computation this will not necessarily be true.
For example, the gradient of one-dimensional linear interpolation can be represented as

    gwi = WeightedIndex(i:i+1, (-1, 1))
    g1 = a[gwi]

For a three-dimensional array `A`, one might compute `∂A/∂x₂` (the second component
of the gradient) as `A[wi1, gwi2, wi3]`, where `wi1` and `wi3` are "value" weights
and `gwi2` "gradient" weights.

`indexes` may be supplied as a range or as a tuple of the same length as `weights`.
The latter is applicable, e.g., for periodic boundary conditions.
"""
abstract type WeightedIndex{L,W} end

# Type to use when array locations are adjacent. This may offer more opportunities
# for compiler optimizations (e.g., SIMD).
struct WeightedAdjIndex{L,W} <: WeightedIndex{L,W}
    istart::Int
    weights::NTuple{L,W}
end
# Type to use with non-adjacent locations. E.g., periodic boundary conditions.
struct WeightedArbIndex{L,W} <: WeightedIndex{L,W}
    indexes::NTuple{L,Int}
    weights::NTuple{L,W}
end

function WeightedIndex(indexes::AbstractUnitRange{<:Integer}, weights::NTuple{L,Any}) where L
    @noinline mismatch(indexes, weights) = throw(ArgumentError("the length of indexes must match weights, got $indexes vs $weights"))
    length(indexes) == L || mismatch(indexes, weights)
    WeightedAdjIndex(first(indexes), promote(weights...))
end
WeightedIndex(istart::Integer, weights::NTuple{L,Any}) where L =
    WeightedAdjIndex(istart, promote(weights...))
WeightedIndex(indexes::NTuple{L,Integer}, weights::NTuple{L,Any}) where L =
    WeightedArbIndex(indexes, promote(weights...))

weights(wi::WeightedIndex) = wi.weights
indexes(wi::WeightedAdjIndex) = wi.istart
indexes(wi::WeightedArbIndex) = wi.indexes
indextuple(wi::WeightedAdjIndex{L}) where L = ntuple(i->wi.istart+i-1, Val(L))
indextuple(wi::WeightedArbIndex{L}) where L = indexes(wi)

# Make them iterable just like numbers are
Base.iterate(x::WeightedIndex) = (x, nothing)
Base.iterate(x::WeightedIndex, ::Any) = nothing
Base.isempty(x::WeightedIndex) = false
Base.length(x::WeightedIndex) = 1

# Supporting arithmetic on the weights allows one to perform pre-scaling of gradient coefficients
Base.:(*)(wi::WeightedAdjIndex, x::Number) = WeightedAdjIndex(wi.istart, wi.weights .* x)
Base.:(/)(wi::WeightedAdjIndex, x::Number) = WeightedAdjIndex(wi.istart, wi.weights ./ x)
Base.:(*)(wi::WeightedArbIndex, x::Number) = WeightedArbIndex(wi.indexes, wi.weights .* x)
Base.:(/)(wi::WeightedArbIndex, x::Number) = WeightedArbIndex(wi.indexes, wi.weights ./ x)

### Indexing with WeightedIndex

# We inject indexing with `WeightedIndex` at a non-exported point in the dispatch heirarchy.
# This is to avoid ambiguities with methods that specialize on the array type rather than
# the index type.
Base.to_indices(A, I::Tuple{Vararg{Union{Int,WeightedIndex}}}) = I
@propagate_inbounds Base._getindex(::IndexLinear, A::AbstractVector, i::Int) = getindex(A, i)  # ambiguity resolution
@inline function Base._getindex(::IndexStyle, A::AbstractArray{T,N}, I::Vararg{Union{Int,WeightedIndex},N}) where {T,N}
    interp_getindex(A, I, ntuple(d->0, Val(N))...)
end

# The non-generated version is currently disabled due to https://github.com/JuliaLang/julia/issues/29117
# # This follows a "move processed indexes to the back" strategy, so J contains the yet-to-be-processed
# # indexes and I all the processed indexes.
# interp_getindex(A::AbstractArray{T,N}, J::Tuple{Int,Vararg{Any,L}}, I::Vararg{Int,M}) where {T,N,L,M} =
#     interp_getindex(A, Base.tail(J), I..., J[1])
# function interp_getindex(A::AbstractArray{T,N}, J::Tuple{WeightedIndex,Vararg{Any,L}}, I::Vararg{Int,M}) where {T,N,L,M}
#     wi = J[1]
#     interp_getindex1(A, indexes(wi), weights(wi), Base.tail(J), I...)
# end
# interp_getindex(A::AbstractArray{T,N}, ::Tuple{}, I::Vararg{Int,N}) where {T,N} =   # termination
#     @inbounds A[I...]  # all bounds-checks have already happened
#
# ## Handle expansion of a single dimension
# # version for WeightedAdjIndex
# @inline interp_getindex1(A, i::Int, weights::NTuple{K,Any}, rest, I::Vararg{Int,M}) where {M,K} =
#     weights[1] * interp_getindex(A, rest, I..., i) + interp_getindex1(A, i+1, Base.tail(weights), rest, I...)
# @inline interp_getindex1(A, i::Int, weights::Tuple{Any}, rest, I::Vararg{Int,M}) where M =
#     weights[1] * interp_getindex(A, rest, I..., i)
# interp_getindex1(A, i::Int, weights::Tuple{}, rest, I::Vararg{Int,M}) where M =
#     error("exhausted the weights, this should never happen")
#
# # version for WeightedArbIndex
# @inline interp_getindex1(A, indexes::NTuple{K,Int}, weights::NTuple{K,Any}, rest, I::Vararg{Int,M}) where {M,K} =
#     weights[1] * interp_getindex(A, rest, I..., indexes[1]) + interp_getindex1(A, Base.tail(indexes), Base.tail(weights), rest, I...)
# @inline interp_getindex1(A, indexes::Tuple{Int}, weights::Tuple{Any}, rest, I::Vararg{Int,M}) where M =
#     weights[1] * interp_getindex(A, rest, I..., indexes[1])
# interp_getindex1(A, indexes::Tuple{}, weights::Tuple{}, rest, I::Vararg{Int,M}) where M =
#     error("exhausted the weights and indexes, this should never happen")

@inline interp_getindex(A::AbstractArray{T,N}, J::Tuple{Int,Vararg{Any,K}}, I::Vararg{Int,N}) where {T,N,K} =
    interp_getindex(A, Base.tail(J), Base.tail(I)..., J[1])
@generated function interp_getindex(A::AbstractArray{T,N}, J::Tuple{WeightedAdjIndex{L,W},Vararg{Any,K}}, I::Vararg{Int,N}) where {T,N,K,L,W}
    ex = :(w[1]*interp_getindex(A, Jtail, Itail..., j))
    for l = 2:L
        ex = :(w[$l]*interp_getindex(A, Jtail, Itail..., j+$(l-1)) + $ex)
    end
    quote
        $(Expr(:meta, :inline))
        Jtail = Base.tail(J)
        Itail = Base.tail(I)
        j, w = J[1].istart, J[1].weights
        $ex
    end
end
@generated function interp_getindex(A::AbstractArray{T,N}, J::Tuple{WeightedArbIndex{L,W},Vararg{Any,K}}, I::Vararg{Int,N}) where {T,N,K,L,W}
    ex = :(w[1]*interp_getindex(A, Jtail, Itail..., ij[1]))
    for l = 2:L
        ex = :(w[$l]*interp_getindex(A, Jtail, Itail..., ij[$l]) + $ex)
    end
    quote
        $(Expr(:meta, :inline))
        Jtail = Base.tail(J)
        Itail = Base.tail(I)
        ij, w = J[1].indexes, J[1].weights
        $ex
    end
end
@inline interp_getindex(A::AbstractArray{T,N}, ::Tuple{}, I::Vararg{Int,N}) where {T,N} =   # termination
    @inbounds A[I...]  # all bounds-checks have already happened

"""
    w = value_weights(degree, δx)

Compute the weights for interpolation of the value at an offset `δx` from the "base" position.
`degree` describes the interpolation scheme.

# Example

```jldoctest
julia> Interpolations.value_weights(Linear(), 0.2)
(0.8, 0.2)
```

This corresponds to the fact that linear interpolation at `x + 0.2` is `0.8*y[x] + 0.2*y[x+1]`.
"""
function value_weights end

"""
    w = gradient_weights(degree, δx)

Compute the weights for interpolation of the gradient at an offset `δx` from the "base" position.
`degree` describes the interpolation scheme.

# Example

```jldoctest
julia> Interpolations.gradient_weights(Linear(), 0.2)
(-1.0, 1.0)
```

This defines the gradient of a linear interpolation at 3.2 as `y[4] - y[3]`.
"""
function gradient_weights end

"""
    w = hessian_weights(degree, δx)

Compute the weights for interpolation of the hessian at an offset `δx` from the "base" position.
`degree` describes the interpolation scheme.

# Example

```jldoctest
julia> Interpolations.hessian_weights(Linear(), 0.2)
(0.0, 0.0)
```

Linear interpolation uses straight line segments, so the second derivative is zero.
"""
function hessian_weights end

gradient1(itp::AbstractInterpolation{T,1}, x) where {T} = gradient(itp, x)[1]
hessian1(itp::AbstractInterpolation{T,1}, x) where {T} = hessian(itp, x)[1]

### Supporting expansion of CartesianIndex

const UnexpandedIndexTypes = Union{Number, AbstractVector, CartesianIndex}
const ExpandedIndexTypes = Union{Number, AbstractVector}

Base.to_index(::AbstractInterpolation, x::Number) = x

# Commented out because you can't add methods to an abstract type.
# @inline function (itp::AbstractInterpolation)(x::Vararg{UnexpandedIndexTypes})
#     itp(to_indices(itp, x)...)
# end
function gradient(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
    xi = to_indices(itp, x)
    xi == x && error("gradient of $itp not supported for position $x")
    gradient(itp, xi...)
end
@propagate_inbounds function gradient!(dest, itp::AbstractInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
    dest .= gradient(itp, x...)
end
function gradient!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
    gradient!(dest, itp, to_indices(itp, x)...)
end
function hessian(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
    hessian(itp, to_indices(itp, x)...)
end
@propagate_inbounds function hessian!(dest, itp::AbstractInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
    dest .= hessian(itp, x...)
end
function hessian!(dest, itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
    hessian!(dest, itp, to_indices(itp, x)...)
end

# @inline function (itp::AbstractInterpolation)(x::Vararg{ExpandedIndexTypes})
#     itp.(Iterators.product(x...))
# end
# function gradient(itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes})
#     map(y->tgradient(itp, y), Iterators.product(x...))
# end
# function hessian(itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes})
#     map(y->thessian(itp, y), Iterators.product(x...))
# end
#
# tgradient(itp, y) = gradient(itp, y...)
# thessian(itp, y) = hessian(itp, y...)

# getindex is supported only for Integer indices (deliberately)
import Base: getindex
@propagate_inbounds getindex(itp::AbstractInterpolation{T,N}, i::Vararg{Integer,N}) where {T,N} = itp(i...)
@propagate_inbounds function getindex(itp::AbstractInterpolation{T,1}, i::Integer, j::Integer) where T
    @boundscheck (j == 1 || Base.throw_boundserror(itp, (i, j)))
    itp(i)
end

@inline checkbounds(::Type{Bool}, itp::AbstractInterpolation, x::Vararg{ExpandedIndexTypes,N}) where N =
    _checkbounds(BoundsCheckStyle(itp), itp, x)

_checkbounds(::CheckWillPass, itp, x) = true
_checkbounds(::NeedsCheck, itp, x) = checklubounds(lbounds(itp), ubounds(itp), x)

checklubounds(ls, us, xs) = _checklubounds(true, ls, us, xs)
_checklubounds(tf::Bool, ls, us, xs::Tuple{Number, Vararg{Any}}) =
    _checklubounds(tf & (ls[1] <= xs[1] <= us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
_checklubounds(tf::Bool, ls, us, xs::Tuple{AbstractVector, Vararg{Any}}) =
    _checklubounds(tf & allbetween(ls[1], xs[1], us[1]), Base.tail(ls), Base.tail(us), Base.tail(xs))
_checklubounds(tf::Bool, ::Tuple{}, ::Tuple{}, ::Tuple{}) = tf

maybe_clamp(itp, xs) = maybe_clamp(BoundsCheckStyle(itp), itp, xs)
maybe_clamp(::NeedsCheck, itp, xs) = map(clamp, xs, lbounds(itp), ubounds(itp))
maybe_clamp(::CheckWillPass, itp, xs) = xs

include("nointerp/nointerp.jl")
include("b-splines/b-splines.jl")
include("gridded/gridded.jl")
include("monotonic/monotonic.jl")
include("extrapolation/extrapolation.jl")
include("scaling/scaling.jl")
include("utils.jl")
include("io.jl")
include("convenience-constructors.jl")
include("deprecations.jl")

end # module


text:
# Gumbel distribution

export Gumbel

@parameterized Gumbel()

basemeasure(::Gumbel{()}) = Lebesgue(ℝ)

@kwstruct Gumbel()

@kwstruct Gumbel(β)

@kwstruct Gumbel(μ,β)

# map affine names to those more common for Gumbel
for N in [(:μ,), (:σ,), (:μ,:σ)]
    G = tuple(replace(collect(N), :σ => :β)...)
    @eval begin
        proxy(d::Gumbel{$G}) = affine(NamedTuple{$N}(values(params(d))), Gumbel())
        logdensity(d::Gumbel{$G}, x) = logdensity(proxy(d), x)
        basemeasure(d::Gumbel{$G}) = basemeasure(proxy(d))
    end
end

function logdensity(d::Gumbel{()} , x)
    return -exp(-x) - x
end

import Base

function Base.rand(rng::AbstractRNG, d::Gumbel{()})
    u = rand(rng)
    -log(-log(u))
end

TV.as(::Gumbel) = asℝ

≪(::Gumbel, ::Lebesgue{X}) where X <: Real = true

distproxy(::Gumbel{()}) = Dists.Gumbel()


text:
module custom_module
  using LinearAlgebra

  function dummy()
    println(">> Julia says... Hello, world ! Function dummy() from module custom_module has been executed !")
  end

  function modString(inputStr)
    println(">> Julia says... Input string: ", inputStr)
    return inputStr * " This is an amazing string"
  end

  function printDict(dict, key1, val1, key2, val2)
    println(">> Julia says...", dict)
    if dict[key1] == val1 && dict[key2] == val2
      return true
    else
      return false
    end
  end

  function dictInvert!(dict, key1, val1, key2, val2)
    dict[key2] = val1
    dict[key1] = val2
    println(">> Julia says...", keys(dict))
    return dict
  end

  function addValue(val1, val2)
    return val1 + val2
  end

  function squareMeBaby(A)
    ## Square array and return the result
    println(">> Julia says...", typeof(A))
    B = A.*A
    return B
  end

  function mutateMeByTen!(A)
    ## Multiple array in place by ten
    lmul!(10, A)
  end

  function tupleTest(tt)
    ## test tuple args with specific values
    if tt.a != 123
      return false
    end
    if tt.b - (-11.11e-3) > 1e-12
      return false
    end
    if tt.c - 67.32147 > 1e-12
      return false
    end
    if tt.d != "azerty"
      return false
    end
    if tt.e != nothing
      return false
    end
    return true
  end

  function makeMyTuple()
    return (A=1, B=2, C=3,)
  end

  function nestedTuples(o)
    if (o.a.dict["A"] == 1.0 && o.a.dict["B"] == 2.0)
      if(o.a.dat == collect(1:1:10))
        if o.b.x == 36 && o.b.y == 48 && o.b.z == 60
          return true
        end
      end
    end
    return false
  end

  mutable struct MyStruct
    a::Int
    b::Float64
    c::Float64
    d::String
    e::Array{Int}
    MyStruct() = new()
  end

  function objectTest(tt::MyStruct)
    ## test tuple args with specific values
    if tt.a != 123
      return false
    end
    if tt.b - (-11.11e-3) > 1e-12
      return false
    end
    if tt.c - 67.32147 > 1e-12
      return false
    end
    if tt.d != "azerty"
      return false
    end
    if tt.e != [111, 2, 3, 4, 5, 6]
      return false
    end
    return true
  end

  export dummy
  export tupleTest
  export modString
  export printDict
  export dictInvert!
  export makeMyTuple
  export squareMeBaby
  export mutateMeByTen!
  export nestedTuples
  export addValue

  export MyStruct
  export objectTest
end

max_stars_repo_name: Clonkk/nimjl

max_issues_repo_name: Clonkk/nimjl

max_forks_repo_name: Clonkk/nimjl


text:
module Spring

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Spring (Up)" => Ahorn.EntityPlacement(
        Maple.Spring
    ),
    "Spring (Left)" => Ahorn.EntityPlacement(
        Maple.SpringRight
    ),
    "Spring (Right)" => Ahorn.EntityPlacement(
        Maple.SpringLeft
    ),
)

function Ahorn.selection(entity::Maple.Spring)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 6, y - 3, 12, 5)
end

function Ahorn.selection(entity::Maple.SpringLeft)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 1, y - 6, 5, 12)
end

function Ahorn.selection(entity::Maple.SpringRight)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 4, y - 6, 5, 12)
end

sprite = "objects/spring/00.png"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Spring, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, -8)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SpringLeft, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 9, -11, rot=pi / 2)
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SpringRight, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 3, 1, rot=-pi / 2)

end


text:
using PyCall
using ImageTransformations: imresize
using MIRTjim: jim
using Plots: plot
# you must install these packages for the python version in your environment path, i.e. PYTHON
torch = pyimport("torch") # import torch
np = pyimport("numpy") # import numpy as np
include("utils.jl") # import useful functions
include("model.jl") # import neural network
include("train.jl")

#------------- create training data and ground truth -------------#
N = 1
Cin = 1
Cout = 1
H = 64
W = 64
x = torch.randn(N, Cin, H, W) # initialize CNN input

tmp = [
   zeros(3,20);
   0 1 0 0 0 0 1 0 0 0 1 1 1 1 0 1 1 1 1 0;
   0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0;
   0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 0;
   0 0 1 1 1 1 0 0 0 0 1 1 0 0 0 0 0 1 1 0;
   zeros(3,20)
   0 0 1 0 0 1 1 1 1 0 0 1 0 0 1 1 1 1 0 0;
   0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0;
   0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0;
   0 0 1 0 0 1 1 1 1 0 0 1 0 0 0 0 0 1 0 0;
   0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0;
   0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0;
   0 0 1 0 0 1 1 1 1 0 0 1 0 0 0 0 0 1 0 0;
   zeros(3,20)
]';

# ground truth
target = kron(tmp, ones(3,3))
target = imresize(target, 64, 64) / maximum(target)
target[target .> 0] .= 1
target = Float32.(target)
# jim(target)
target_t = reshape(target, N, Cout, H, W)
target_t = py"numpy2tensor($target_t)" # convert to tensor

#------------------ training the neural network ------------------#
model = py"My_model($Cin, $Cout)" # initialize CNN, input channel is 3, output channel is 1
print("CNN training param =", py"cal_param($model)")

trainer = py"Trainer($model, num_epochs = 500)"
trainer.train(x, target_t)

#------------------ visualize training loss ------------------#
loss = torch.load("./checkpoint_res/train_loss.pt")
plot(loss, label = "training loss", xlabel = "# of epochs", ylabel = "MSE") # converged!

#------------------ test the trained neural network ------------------#
test_model = py"My_model($Cin, $Cout)" # initialize test model
ckpt_path = pwd() * "/checkpoint_res/last_checkpoint.pt"
ckpt = py"load_ckpt($ckpt_path)"
test_model.load_state_dict(ckpt)
cnn_output = py"test_model($test_model, $x)"
cnn_output = reshape(py"tensor2numpy($cnn_output)", H, W) # convert to julia arrays
jim(jim(cnn_output, "CNN"), jim(target, "GT")) # should be really good!


text:
# On Zach's machine, this outputs
#=
zach@Theresa:~/Desktop$ julia vi_allocation_test.jl
 51.437988 seconds (85.34 M allocations: 2.685 GB, 0.51% gc time)
 53.665579 seconds (74.02 M allocations: 2.062 GB, 0.39% gc time)
=#

using POMDPModels
using POMDPs
using DiscreteValueIteration
using POMDPToolbox

function solve_no_allocation{S,A}(solver::ValueIterationSolver, mdp::Union{MDP{S,A},POMDP{S,A}}, policy=create_policy(solver, mdp); verbose::Bool=false)

    # solver parameters
    max_iterations = solver.max_iterations
    belres = solver.belres
    discount_factor = discount(mdp)

    # intialize the utility and Q-matrix
    util = policy.util
    qmat = policy.qmat
    include_Q = policy.include_Q
    pol = policy.policy

    # pre-allocate the transtion distirbution and the interpolants
    dist = create_transition_distribution(mdp)

    total_time = 0.0
    iter_time = 0.0

    # create an ordered list of states for fast iteration
    states = ordered_states(mdp)

    # main loop
    for i = 1:max_iterations
        residual = 0.0
        tic()
        # state loop
        for (istate,s) in enumerate(states)
            old_util = util[istate] # for residual
            max_util = -Inf
            # action loop
            # util(s) = max_a( R(s,a) + discount_factor * sum(T(s'|s,a)util(s') )
            for (iaction, a) in enumerate(policy.action_map)
                dist = transition(mdp, s, a) # fills distribution over neighbors
                u = 0.0
                for sp in iterator(dist)
                    p = pdf(dist, sp)
                    p == 0.0 ? continue : nothing # skip if zero prob
                    r = reward(mdp, s, a, sp)
                    sidx = state_index(mdp, sp)
                    u += p * (r + discount_factor * util[sidx])
                end
                new_util = u
                if new_util > max_util
                    max_util = new_util
                    pol[istate] = iaction
                end
                include_Q ? (qmat[istate, iaction] = new_util) : nothing
            end # action
            # update the value array
            util[istate] = max_util
            diff = abs(max_util - old_util)
            diff > residual ? (residual = diff) : nothing
        end # state
        iter_time = toq()
        total_time += iter_time
        verbose ? println("Iteration : $i, residual: $residual, iteration run-time: $iter_time, total run-time: $total_time") : nothing
        residual < belres ? break : nothing
    end # main
    policy
end

function solve_pre_allocation{S,A}(solver::ValueIterationSolver, mdp::Union{MDP{S,A},POMDP{S,A}}, policy=create_policy(solver, mdp); verbose::Bool=false)

    # solver parameters
    max_iterations = solver.max_iterations
    belres = solver.belres
    discount_factor = discount(mdp)

    # intialize the utility and Q-matrix
    util = policy.util
    qmat = policy.qmat
    include_Q = policy.include_Q
    pol = policy.policy

    # pre-allocate the transtion distirbution and the interpolants
    dist = create_transition_distribution(mdp)

    total_time = 0.0
    iter_time = 0.0

    # create an ordered list of states for fast iteration
    states = ordered_states(mdp)

    # main loop
    for i = 1:max_iterations
        residual = 0.0
        tic()
        # state loop
        for (istate,s) in enumerate(states)
            old_util = util[istate] # for residual
            max_util = -Inf
            # action loop
            # util(s) = max_a( R(s,a) + discount_factor * sum(T(s'|s,a)util(s') )
            for (iaction, a) in enumerate(policy.action_map)
                dist = transition(mdp, s, a, dist) # fills distribution over neighbors
                u = 0.0
                for sp in iterator(dist)
                    p = pdf(dist, sp)
                    p == 0.0 ? continue : nothing # skip if zero prob
                    r = reward(mdp, s, a, sp)
                    sidx = state_index(mdp, sp)
                    u += p * (r + discount_factor * util[sidx])
                end
                new_util = u
                if new_util > max_util
                    max_util = new_util
                    pol[istate] = iaction
                end
                include_Q ? (qmat[istate, iaction] = new_util) : nothing
            end # action
            # update the value array
            util[istate] = max_util
            diff = abs(max_util - old_util)
            diff > residual ? (residual = diff) : nothing
        end # state
        iter_time = toq()
        total_time += iter_time
        verbose ? println("Iteration : $i, residual: $residual, iteration run-time: $iter_time, total run-time: $total_time") : nothing
        residual < belres ? break : nothing
    end # main
    policy
end

mdp = GridWorld()
solver = ValueIterationSolver()
solve_no_allocation(solver, mdp, ValueIterationPolicy(mdp, include_Q=false))
@time for i in 1:100
    solve_no_allocation(solver, mdp, ValueIterationPolicy(mdp, include_Q=false))
end

mdp = GridWorld()
solver = ValueIterationSolver()
solve_pre_allocation(solver, mdp, ValueIterationPolicy(mdp, include_Q=false))
@time for i in 1:100
    solve_pre_allocation(solver, mdp, ValueIterationPolicy(mdp, include_Q=false))
end


text:
"""
generic type for discrete function Dom -> Rng
"""
struct DiscreteFunction{D <: AbstractFiniteSet, R <: AbstractFiniteSet, A, B} <: AbstractDiscreteFunction
    m :: Dict{A, B}
    dom :: D
    rng :: R
end
const DFun = DiscreteFunction

function DiscreteFunction(d :: Dict{A, B}) where {A, B}
    min_x = minimum(keys(d))
    max_x = maximum(keys(d))
    min_y = minimum(values(d))
    max_y = maximum(values(d))
    D = Segment(min_x, max_x)
    R = Segment(min_y, max_y)
    return DiscreteFunction(d, D, R)
end

domain(df :: DiscreteFunction) = df.dom
codomain(df :: DiscreteFunction) = df.rng

(df :: DiscreteFunction)(x) = f.m[x]
Base.setindex!(f :: DiscreteFunction, value, key) = setindex!(f.m, value, key)

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
module SPCSpectra

using Dates

export SPC

mutable struct SPC
    data::Vector{Tuple{Vector{<:Number}, Vector{Float32}}} # xdata can be Float64 if calculated from range specs
    zdata::Vector{Float32}
    xlabel::String
    ylabel::String
    zlabel::String
    experimenttype::String
    timestamp::DateTime
    param_dict::Dict{String, String}
    params::Vector{String}
end

# Adaptation of the pythonpackage spc by rohanisaac

# byte positon of various parts of the file
# head_siz = 512
# old_head_siz = 256
subhead_siz = 32
log_siz = 64

"Units for x,z,w axes."
const fxtype_op = ["Arbitrary",
                    "Wavenumber (cm-1)",
                    "Micrometers (um)",
                    "Nanometers (nm)",
                    "Seconds ",
                    "Minutes", "Hertz (Hz)",
                    "Kilohertz (KHz)",
                    "Megahertz (MHz) ",
                    "Mass (M/z)",
                    "Parts per million (PPM)",
                    "Days",
                    "Years",
                    "Raman Shift (cm-1)",
                    "eV",
                    "XYZ text labels in fcatxt (old 0x4D version only)",
                    "Diode Number",
                    "Channel",
                    "Degrees",
                    "Temperature (F)",
                    "Temperature (C)",
                    "Temperature (K)",
                    "Data Points",
                    "Milliseconds (mSec)",
                    "Microseconds (uSec) ",
                    "Nanoseconds (nSec)",
                    "Gigahertz (GHz)",
                    "Centimeters (cm)",
                    "Meters (m)",
                    "Millimeters (mm)",
                    "Hours"]

"Units y-axis."
const fytype_op = ["Arbitrary Intensity",
                    "Interferogram",
                    "Absorbance",
                    "Kubelka-Munk",
                    "Counts",
                    "Volts",
                    "Degrees",
                    "Milliamps",
                    "Millimeters",
                    "Millivolts",
                    "Log(1/R)",
                    "Percent",
                    "Intensity",
                    "Relative Intensity",
                    "Energy",
                    "",
                    "Decibel",
                    "",
                    "",
                    "Temperature (F)",
                    "Temperature (C)",
                    "Temperature (K)",
                    "Index of Refraction [N]",
                    "Extinction Coeff. [K]",
                    "Real",
                    "Imaginary",
                    "Complex"]

const fytype_op2 = ["Transmission",
                    "Reflectance",
                    "Arbitrary or Single Beam with Valley Peaks",
                    "Emission"]

const fexper_op = ["General SPC",
                    "Gas Chromatogram",
                    "General Chromatogram",
                    "HPLC Chromatogram",
                    "FT-IR, FT-NIR, FT-Raman Spectrum or Igram",
                    "NIR Spectrum",
                    "UV-VIS Spectrum",
                    "X-ray Diffraction Spectrum",
                    "Mass Spectrum ",
                    "NMR Spectrum or FID",
                    "Raman Spectrum",
                    "Fluorescence Spectrum",
                    "Atomic Spectrum",
                    "Chromatography Diode Array Spectra"]

flag_bits(n) = BitVector(digits(n, base = 2, pad = 8 * sizeof(n)))

read_data(io::IO, T::DataType) = ltoh(read(io, T))
read_data(io::IO, ::Type{String}, n::Integer) = strip(String(read(io, n)), '\0')
read_data(io::IO, T::DataType, n::Integer...) = ltoh.(reshape(reinterpret(T, read(io, prod(n) * sizeof(T))), Int64.(n)...))
read_data(io::IO, TT::Union{NTuple{N, DataType} where N, Vector{DataType}}) = ltoh.(read.(Ref(io), TT))

"""
    SPC(filename::AbstractString)

Construct SPC objects.
"""
function SPC(filename::AbstractString)
    content = read(filename)
    io = IOBuffer(content)

    ftflg, fversn = content[1:2]
    # --------------------------------------------
    # NEW FORMAT (LSB)
    # --------------------------------------------
    fversn == 0x4b || return "Reading of file version $(repr(fversn)) not implemented"

    ftflg = read_data(io, UInt8)
    fversn = read_data(io, UInt8)
    fexper = read_data(io, UInt8)
    fexp = read_data(io, UInt8)
    fnpts = read_data(io, Int32)
    ffirst = read_data(io, Float64)
    flast = read_data(io, Float64)
    fnsub = read_data(io, Int32)
    fxtype = read_data(io, UInt8)
    fytype = read_data(io, UInt8)
    fztype = read_data(io, UInt8)
    fpost = read_data(io, UInt8)
    fdate = read_data(io, Int32)

    fres = read_data(io, String, 9)
    fsource = read_data(io, String, 9)
    fpeakpt = read_data(io, Int16)
    fspare  = read_data(io, String, 32)
    fcmnt = read_data(io, String, 130)
    fcatxt = read_data(io, String, 30)
    flogoff = read_data(io, Int32)
    fmods = read_data(io, Int32)
    fprocs = read_data(io, UInt8)
    flevel = read_data(io, UInt8)
    fsampin = read_data(io, Int16)
    ffactor = read_data(io, Float32)
    fmethod = read_data(io, String, 48)
    fzinc = read_data(io, Float32)
    fwplanes = read_data(io, Int32)
    fwinc = read_data(io, Float32)
    fwtype = read_data(io, UInt8)
    freser = read_data(io, String, 187)

    # Flag bits
    tsprec, tcgram, tmulti, trandm, tordrd, talabs, txyxys, txvals = flag_bits(ftflg)

    # Convert date time to appropriate format
    year = fdate >> 20
    month = (fdate >> 16) % (2^4)
    day = (fdate >> 11) % (2^5)
    hour = (fdate >> 6) % (2^5)
    minute = fdate % (2^6)
    timestamp = DateTime(year, month, day, hour, minute)

    # remove multiple spaces
    cmnt = replace(fcmnt, r"\s+" => " ")

    # figure out type of file
    dat_multi = fnsub > 1

    dat_fmt = if txyxys
        # x values are given
        "-xy"
    elseif txvals
        # only one subfile, which contains the x data
        dat_fmt = "x-y"
    else
        # no x values are given, but they can be generated
        dat_fmt = "gx-y"
    end

    println("$dat_fmt($fnsub)")

    x = if ! txyxys
        # txyxys don't have global x data
        if txvals
            # if global x data is given
            read_data(io, Float32, fnpts)
        else
            # otherwise generate them
            range(ffirst, flast; length=fnpts) |> collect
        end
    end
    # make a list of subfiles
    xydata = []
    zdata = Float32[]
    z0 = dz = 0f0

    # if subfile directory is given
    if dat_fmt == "-xy" && fnpts > 0
        directory = true
        # loop over entries in directory
        for i in 1:fnsub
            ssfposn, ssfsize, ssftime = read_data(io, (Int32, Int32, Float32))
            # add sufile, load defaults for npts and exp
            pos = position(io)
            seek(io, ssfposn) # io buffer position is zero-based!
            xloc, y, z, zinc = subFile(io, 0, 0, true, tsprec, tmulti)
            if i == 1
                z0 = z
                dz = zinc
            end
            seek(io, pos)
            push!(xydata, (isnothing(xloc) ? x : xloc, y))
            push!(zdata, z)
        end
    else
        # don't have directory, for each subfile
        for i in 1:fnsub
            xloc, y, z, zinc = subFile(io, fnpts, fexp, txyxys, tsprec, tmulti)
            if i == 1
                z0 = z
                dz = zinc
            end
            push!(xydata, (isnothing(xloc) ? x : xloc, y))
            push!(zdata, z)
        end
    end
    # if log data exists
    # flog offset to log data offset not zero (bytes)
    param_dict = Dict{String, String}()
    params = String[]  # put the rest into a list
    if flogoff > 0
        log_head_end = flogoff + log_siz
        io_log = IOBuffer(content[flogoff+1:log_head_end])
        # logstc_str = "<iiiii44s"
        logsizd, logsizm, logtxto, logbins, logdsks = read_data(io_log, Int32, 5)
        logspar = read_data(io_log, String, 44)
        log_pos = flogoff + logtxto

        log_end_pos = flogoff + logsizd

        # line endings: get rid of any '\r' and then split on '\n'
        log_content = split(strip(String(content[log_pos + 1:log_end_pos]), ['\0', '\r', '\n']), r"\r?\n")

        # split log data into dictionary based on =
        for x in log_content
            if occursin("=", x)
                # stop it from breaking if there is more than 1 =
                key, value = split(x, "=")[1:2]
                push!(param_dict, key => strip(value, '\0'))
            else
                push!(params, x)
            end
        end
    end

    labels = [
        get(fxtype_op, fxtype + 1, "Unknown"),
        get(fytype_op, fytype + 1, get(fytype_op2, fytype - 127, "Unknown")),
        get(fxtype_op, fztype + 1, "Unknown")
    ]

    # --------------------------
    # check if labels are included as text
    # --------------------------

    # split it based on '\0' character
    # format x, y, z
    if talabs
        for (i, s) in enumerate(split(fcatxt, '\0', keepempty = false))
            isempty(s) || (labels[i] = s)
        end
    end

    if (0x10 & ftflg == 0x0)
        fzinc > 0 && (dz == fzinc)
        zdata = z0 .+ collect(0:dz:((fnpts - 1) * dz))
    end
    SPC(xydata, zdata, labels..., get(fexper_op, fexper + 1, "Unknown"), timestamp, param_dict, params)
end

"""
    subFile(io::IO, fnpts, fexp, txyxy, tsprec, tmulti)

Process each subfile passed to it, extracts header information and data
information and places them in data members
Data
- x: x-data (optional)
- y: y-data
- y_int: integer y-data if y-data is not floating
"""
function subFile(io::IO, fnpts, fexp, txyxy, tsprec, tmulti)
    # extract subheader info
    subflgs, subexp, subindx, subtime, subnext, subnois, subnpts, subscan, subwlevel, subresv = read_subheader(io)

    pts = txyxy ? subnpts : fnpts

    # Choosing exponent
    # -----------------
    # choose local vs global exponent depending on tmulti
    exp = tmulti ? subexp : fexp

    # Make sure it is reasonable, if it out of range zero it
    (-128 < exp <= 128) || (exp = 0)

    # --------------------------
    # if x_data present
    # --------------------------
    x = if txyxy
        # x_str = '<' + 'i' * pts
        x_raw = read_data(io, Int32, pts)
        (2.0f0^(exp - 32)) .* x_raw
    else
        nothing
    end
    # --------------------------
    # extract y_data
    # --------------------------
    y = if exp == 128
        # Floating y-values
        read_data(io, Float32, pts)
    else
        # integer format
        if tsprec
            # 16 bit
            y_raw16 = read_data(io, Int16, pts)
            (2.0f0^(exp - 16)) .* y_raw16
        else
            # 32 bit, using size of subheader to figure out data type
            y_raw = read_data(io, Int32, pts)
            (2.0f0^(exp - 32)) .* y_raw
        end
    end

    z = subtime
    zinc = subnext - subtime
    x, y, z, zinc
end

"""
    read_subheader(io::IO)

Return the subheader as a list:
-------
10 item list with the following data members:
    [1] subflgs
    [2] subexp
    [3] subindx
    [4] subtime
    [5] subnext
    [6] subnois
    [7] subnpts
    [8] subscan
    [9] subwlevel
    [10] subresv
"""
function read_subheader(io::IO)
    subflgs = read_data(io, UInt8)
    subexp = read_data(io, UInt8)
    subindx = read_data(io, Int16)
    subtime = read_data(io, Float32)
    subnext = read_data(io, Float32)
    subnois = read_data(io, Float32)
    subnpts = read_data(io, Int32)
    subscan = read_data(io, Int32)
    subwlevel = read_data(io, Float32)
    subresv = read_data(io, String, 4)

    subflgs, subexp, subindx, subtime, subnext, subnois, subnpts, subscan, subwlevel, subresv
end

end # module