text:
## Common:

### backward compat

function DecisionRule(exo_grid, endo_grid, values::Vector{Matrix{Float64}})
    N,n_x = size(values[1])
    vals = [to_LOP(v) for v in values]
    DecisionRule(exo_grid, endo_grid, vals)
end

function DecisionRule(exo_grid, endo_grid, values::Vector{Matrix{Float64}}, tt)
    N,n_x = size(values[1])
    vals = [reshape(reinterpret(Value{n_x}, vec(v')), (N,)) for v in values]
    DecisionRule(exo_grid, endo_grid, vals, tt)
end

# EmptyGrid x UCGrid

function evaluate(dr::AbstractDecisionRule{EmptyGrid{d1}, UCGrid{d}, n_x}, z::AbstractMatrix{Float64}) where n_x where d where d1
    N = size(z,1)
    @assert size(z,2)==d
    points = to_LOP(Point{d}, z)
    out = evaluate(dr,copy(points))
    return reshape(reinterpret(Float64, vec(out)), (n_x,N))'
end

evaluate(dr::AbstractDecisionRule{EmptyGrid{d1}, <:UCGrid}, z::Vector{Float64})  where d1 = vec(evaluate(dr,z'))
evaluate(dr::AbstractDecisionRule{EmptyGrid{d1}, <:UCGrid}, i::Int, y::Union{Vector,AbstractMatrix})  where d1 = evaluate(dr,y) # kind of nonsensical
evaluate(dr::AbstractDecisionRule{EmptyGrid{d1}, <:UCGrid}, x::Union{Vector,AbstractMatrix}, y::Union{Vector,AbstractMatrix})  where d1 = evaluate(dr,y)
UCGrid
function set_values!(dr::AbstractDecisionRule{EmptyGrid{d1}, <:UCGrid, n_x}, vals::Vector{Matrix{Float64}}) where n_x where d1
    N = size(vals[1],1)
    dims = tuple(dr.grid_endo.n...)
    V = [reshape(reinterpret(Value{n_x}, vec(vals[1]')), dims)]
    set_values!(dr, V)
end

# UCGrid x UCGrid

function evaluate(dr::AbstractDecisionRule{UCGrid{d1}, UCGrid{d2}, n_x}, x::AbstractMatrix, y::AbstractMatrix) where d1 where d2 where n_x
    N = size(x,1)
    @assert size(y,1)==N
    xx = to_LOP(Point{d1}, x)
    yy = to_LOP(Point{d2}, y)
    res = evaluate(dr,copy(xx),copy(yy))
    return reshape(reinterpret(Float64,vec(res)), (n_x, N))'
end
#
evaluate(dr::AbstractDecisionRule{<:UCGrid, <:UCGrid}, x::Vector{Float64}, y::Matrix{Float64}) = evaluate(dr, repeat(x', size(y,1)), y)
evaluate(dr::AbstractDecisionRule{<:UCGrid, <:UCGrid}, x::Vector{Float64}, y::Vector{Float64}) = evaluate(dr, vector_to_matrix(x), vector_to_matrix(y))
evaluate(dr::AbstractDecisionRule{UCGrid{d1}, UCGrid{d2}}, z::AbstractMatrix) where d1 where d2 = evaluate(dr, [z[:,1:d1] z[:,1:(d1+1:end)]])
evaluate(dr::AbstractDecisionRule{<:UCGrid, <:UCGrid}, z::Vector{Float64}) = vec(evaluate(dr,vector_to_matrix(z)))
evaluate(dr::AbstractDecisionRule{UCGrid{d1}, UCGrid{d2}}, i::Int, y::Union{Vector{Float64},AbstractMatrix{Float64}}) where d1 where d2  = evaluate(dr,node(dr.grid_exo,i), y)

function set_values!(dr::AbstractDecisionRule{<:UCGrid, <:UCGrid, n_x}, vals::Vector{Matrix{Float64}}) where n_x
    N = size(vals[1],1)
    dims = tuple(dr.grid_endo.n...)
    V = [reshape(reinterpret(Value{n_x}, vec(v')), dims) for v in vals]
    set_values!(dr, V)
end

# UnstructuredGrid x UCGrid

function evaluate(dr::AbstractDecisionRule{UnstructuredGrid{d1}, UCGrid{d2},n_x}, i::Int, z::AbstractMatrix{Float64}) where n_x where d1 where d2
    N = size(z,1)
    @assert size(z,2)==d2
    points = to_LOP(Point{d2}, z)
    # TODO: remove copy()
    out = evaluate(dr,i,copy(points))
    reshape(reinterpret(Float64, vec(out)), (n_x, N))'
end

evaluate(dr::AbstractDecisionRule{UnstructuredGrid{d1}, UCGrid{d2}, n_x}, i::Int, x::Vector{Float64}) where n_x where d1 where d2 = vec(evaluate(dr,i,vector_to_matrix(x)))
function evaluate(dr::AbstractDecisionRule{UnstructuredGrid{d1}, UCGrid{d2}, n_x}, inds::Vector{Int}, x::AbstractMatrix{Float64}) where n_x where d1 where d2
    l = [evaluate(dr,inds[n],x[n,:])' for n=1:size(x,1)]
    vcat(l...)
end

function evaluate(dr::AbstractDecisionRule{UnstructuredGrid{d1}, UCGrid{d2}, n_x}, inds::Vector{Int}, x::Vector{Float64}) where n_x where d1 where d2
    vcat([evaluate(dr,i,x)' for i in inds]...)
end

function set_values!(dr::AbstractDecisionRule{<:UnstructuredGrid, <:UCGrid, n_x}, vals::Vector{Matrix{Float64}}) where n_x
    N = size(vals[1],1)
    dims = tuple(dr.grid_endo.n...)
    V = [to_LOP(Value{n_x}, v) for v in vals]
    set_values!(dr, V)
end


text:
using Test
using ForwardDiff: Dual, gradient, partials

using DiffEqBase: ODE_DEFAULT_NORM
const internalnorm = ODE_DEFAULT_NORM

val = rand(10)
par = rand(10)
u = Dual.(val, par)
reference(val, par) = sqrt((sum(abs2,val) + sum(abs2,par)) / (length(val) + length(par)))
dual_real = internalnorm(u, 1)
dual_dual = internalnorm(u, u[1])
@test reference(val, par) ≈ dual_real
@test reference(val, par) ≈ dual_dual
@test partials(dual_dual, 1) ≈ gradient(x->internalnorm(x, x[1]), val)'par

ext: jl


text:
##################################################################
"""
$(TYPEDEF)

Extendable sparse matrix. A nonzero  entry of this matrix is contained
either in cscmatrix, or in lnkmatrix, never in both.

$(TYPEDFIELDS)
"""
mutable struct ExtendableSparseMatrix{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    """
    Final matrix data
    """
    cscmatrix::SparseMatrixCSC{Tv,Ti}

    """
    Linked list structure holding data of extension
    """
    lnkmatrix::Union{SparseMatrixLNK{Tv,Ti},Nothing}

    """
    Pattern hash
    """
    phash::UInt64
end

"""
```
ExtendableSparseMatrix(Tv,Ti,m,n)
ExtendableSparseMatrix(Tv,m,n)
ExtendableSparseMatrix(m,n)
```
Create empty ExtendableSparseMatrix. This is equivalent to `spzeros(m,n)` for
`SparseMartrixCSC`.

"""
function ExtendableSparseMatrix{Tv,Ti}(m, n) where{Tv,Ti<:Integer}
    ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,0)
end

function ExtendableSparseMatrix(valuetype::Type{Tv},indextype::Type{Ti},m, n) where{Tv,Ti<:Integer}
    ExtendableSparseMatrix{Tv,Ti}(m,n)
end

ExtendableSparseMatrix(valuetype::Type{Tv},m, n) where{Tv}=ExtendableSparseMatrix{Tv,Int}(m,n)
ExtendableSparseMatrix(m, n)=ExtendableSparseMatrix{Float64,Int}(m,n)

"""
$(SIGNATURES)

 Create ExtendableSparseMatrix from SparseMatrixCSC
"""
function ExtendableSparseMatrix(csc::SparseMatrixCSC{Tv,Ti}) where{Tv,Ti<:Integer}
    return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, phash(csc))
end

"""
$(SIGNATURES)

Create similar but emtpy extendableSparseMatrix
"""
Base.similar(m::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}=ExtendableSparseMatrix{Tv,Ti}(size(m)...)

"""
$(SIGNATURES)

Update element of the matrix  with operation `op`.
This can replace the following code and save one index
search during acces:

```@example
using ExtendableSparse # hide
A=ExtendableSparseMatrix(3,3)
A[1,2]+=0.1
A
```

```@example
using ExtendableSparse # hide

A=ExtendableSparseMatrix(3,3)
updateindex!(A,+,0.1,1,2)
A
```

If `v` is zero, no new entry is created.
"""
function updateindex!(ext::ExtendableSparseMatrix{Tv,Ti}, op,v, i,j) where {Tv,Ti<:Integer}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        ext.cscmatrix.nzval[k]=op(ext.cscmatrix.nzval[k],v)
    else
        if ext.lnkmatrix==nothing
            ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
        end
        updateindex!(ext.lnkmatrix,op,v,i,j)
    end
    ext
end

"""
$(SIGNATURES)
Like [`updateindex!`](@ref) but without
checking if v is zero.
"""
function rawupdateindex!(ext::ExtendableSparseMatrix{Tv,Ti}, op,v, i,j) where {Tv,Ti<:Integer}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        ext.cscmatrix.nzval[k]=op(ext.cscmatrix.nzval[k],v)
    else
        if ext.lnkmatrix==nothing
            ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
        end
        rawupdateindex!(ext.lnkmatrix,op,v,i,j)
    end
    ext
end

"""
$(SIGNATURES)

Find index in CSC matrix and set value if it exists. Otherwise,
set index in extension if `v` is nonzero.
"""
function Base.setindex!(ext::ExtendableSparseMatrix{Tv,Ti}, v, i,j) where {Tv,Ti}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        ext.cscmatrix.nzval[k]=v
    else
        if ext.lnkmatrix==nothing
            ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
        end
        ext.lnkmatrix[i,j]=v
    end
end

"""
$(SIGNATURES)

Find index in CSC matrix and return value, if it exists.
Otherwise, return value from extension.
"""
function Base.getindex(ext::ExtendableSparseMatrix{Tv,Ti},i, j)  where{Tv,Ti<:Integer}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        return ext.cscmatrix.nzval[k]
    elseif ext.lnkmatrix==nothing
        return zero(Tv)
    else
        return ext.lnkmatrix[i,j]
    end
end

"""
$(SIGNATURES)

Size of ExtendableSparseMatrix.
"""
Base.size(ext::ExtendableSparseMatrix) = (ext.cscmatrix.m, ext.cscmatrix.n)

"""
$(SIGNATURES)

Show matrix
"""
function Base.show(io::IO,::MIME"text/plain",ext::ExtendableSparseMatrix)
    flush!(ext)
    xnnz = nnz(ext)
    m, n = size(ext)
    print(io, m, "×", n, " ", typeof(ext), " with ", xnnz, " stored ",
          xnnz == 1 ? "entry" : "entries")
    if !(m == 0 || n == 0 || xnnz==0)
        print(io, ":")
        show(IOContext(io), ext.cscmatrix)
    end
end

"""
$(SIGNATURES)

If there are new entries in extension, create new CSC matrix by adding the
cscmatrix and linked list matrix and reset the linked list based extension.
"""
function flush!(ext::ExtendableSparseMatrix)
    if ext.lnkmatrix!=nothing && nnz(ext.lnkmatrix)>0
        ext.cscmatrix=ext.lnkmatrix+ext.cscmatrix
        ext.lnkmatrix=nothing
        ext.phash=phash(ext.cscmatrix)
    end
    return ext
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return number of nonzeros in ext.cscmatrix.
"""
function SparseArrays.nnz(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return nnz(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return nonzeros in ext.cscmatrix.
"""
function SparseArrays.nonzeros(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return nonzeros(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return rowvals in ext.cscmatrix.
"""
function SparseArrays.rowvals(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return rowvals(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return colptr of  in ext.cscmatrix.
"""
function colptrs(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return ext.cscmatrix.colptr
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return findnz(ext.cscmatrix).
"""
function SparseArrays.findnz(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return findnz(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`\\`](@ref) for extmatrix
"""
function LinearAlgebra.:\(ext::ExtendableSparseMatrix,B::AbstractVecOrMat{T} where T)
    flush!(ext)
    ext.cscmatrix\B
end

"""
$(SIGNATURES)

[`\\`](@ref) for Symmetric{ExtendableSparse}
"""
function LinearAlgebra.:\(symm_ext::Symmetric{Tm, ExtendableSparseMatrix{Tm, Ti}}, B::AbstractVecOrMat{T} where T) where{Tm,Ti}
    flush!(symm_ext.data)
    symm_csc=Symmetric(symm_ext.data.cscmatrix,Symbol(symm_ext.uplo))
    symm_csc\B
end

"""
$(SIGNATURES)

[`\\`](@ref) for Hermitian{ExtendableSparse}
"""
function LinearAlgebra.:\(symm_ext::Hermitian{Tm, ExtendableSparseMatrix{Tm, Ti}}, B::AbstractVecOrMat{T} where T) where{Tm,Ti}
    flush!(symm_ext.data)
    symm_csc=Hermitian(symm_ext.data.cscmatrix,Symbol(symm_ext.uplo))
    symm_csc\B
end

"""
$(SIGNATURES)

[`flush!`](@ref) and ldiv with ext.cscmatrix
"""
function  LinearAlgebra.ldiv!(r, ext::ExtendableSparse.ExtendableSparseMatrix, x)
    @inbounds flush!(ext)
    return LinearAlgebra.ldiv!(r,ext.cscmatrix,x)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and multiply with ext.cscmatrix
"""
function  LinearAlgebra.mul!(r,ext::ExtendableSparse.ExtendableSparseMatrix, x)
    @inbounds flush!(ext)
    return LinearAlgebra.mul!(r,ext.cscmatrix,x)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and calculate norm from cscmatrix
"""
function LinearAlgebra.norm(A::ExtendableSparseMatrix, p::Real=2)
    @inbounds flush!(A)
    return LinearAlgebra.norm(A.cscmatrix,p)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and calculate opnorm from cscmatrix
"""
function LinearAlgebra.opnorm(A::ExtendableSparseMatrix, p::Real=2)
    @inbounds flush!(A)
    return LinearAlgebra.opnorm(A.cscmatrix,p)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and calculate cond from cscmatrix
"""
function LinearAlgebra.cond(A::ExtendableSparseMatrix, p::Real=2)
    @inbounds flush!(A)
    return LinearAlgebra.cond(A.cscmatrix,p)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and check for symmetry of cscmatrix
"""
function LinearAlgebra.issymmetric(A::ExtendableSparseMatrix)
    @inbounds flush!(A)
    return LinearAlgebra.issymmetric(A.cscmatrix)
end

"""
$(SIGNATURES)

Add SparseMatrixCSC matrix and [`ExtendableSparseMatrix`](@ref)  ext.
"""
function Base.:+(ext::ExtendableSparseMatrix,csc::SparseMatrixCSC)
    @inbounds flush!(ext)
    return ext.cscmatrix+csc
end

"""
$(SIGNATURES)

Subtract  SparseMatrixCSC matrix from  [`ExtendableSparseMatrix`](@ref)  ext.
"""
function Base.:-(ext::ExtendableSparseMatrix,csc::SparseMatrixCSC)
    @inbounds flush!(ext)
    return ext.cscmatrix-csc
end

"""
$(SIGNATURES)

Subtract  [`ExtendableSparseMatrix`](@ref)  ext from  SparseMatrixCSC.
"""
function Base.:-(csc::SparseMatrixCSC,ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return csc - ext.cscmatrix
end

"""
$(SIGNATURES)
"""
function SparseArrays.dropzeros!(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    dropzeros!(ext.cscmatrix)
end


text:
include("basic_omnivory_module.jl")
include("press.jl")
using PyPlot
pygui(true)

# PULSE SIMULATIONS for ranges of K, aCP, eCP and mP.

# range of K
println("Simulations for a range of K values")
res_K = []
rg_K = 2:0.05:4
for k in rg_K
    par = ModelPar()
    par.K_base = par.K = k
    res_K = push!(res_K, press(par, 0.1, 0.5))
end

# range of a_CP
println("Simulations for a range of aCP values")
res_aCP = []
rg_aCP = 0.3:.0025:.55
for a in rg_aCP
    par = ModelPar()
    par.a_CP = a
    res_aCP = push!(res_aCP, press(par, 0.1, 0.5))
end

# range of e_CP
println("Simulations for a range of eCP values")
res_eCP = []
rg_eCP = 0.35:.0025:0.7
for e in rg_eCP
    par = ModelPar()
    par.e_CP = e
    res_eCP = push!(res_eCP, press(par, 0.1, 0.5))
end

# range of m_P
println("Simulations for a range of mP values")
res_mP = []
rg_mP = 0.15:.005:.3
for m in rg_mP
    par = ModelPar()
    par.m_P = m
    res_mP = push!(res_mP, press(par, 0.1, 0.5))
end

# As this takes some time to run you may want to save this (it takes ~500Mo),
# if so uncomment the lines below.
# using FileIO, JLD2
# # mkdir("res")
# save(
#     "res/res_S2.jld2",
#     "rg_K", rg_K, "res_K", res_K,
#     "rg_aCP", rg_aCP, "res_aCP", res_aCP,
#     "rg_eCP", rg_eCP, "res_eCP", res_eCP,
#     "rg_mP", rg_mP, "res_mP", res_mP
#     )

# FIGURE
println("Drawing Figure S2")

fig = figure(figsize = (12, 8))
ind = [2 5 8]

##----------- K
tly = ["-1/Re(λ1)" "Degree of Overshoot" "Max-Min"]
tlx = ["" "" "K"]
for i in 1:3
    subplot(3, 4, (i - 1) * 4 + 1)
    plot_sa_unit(rg_K, res_K, ind[i], false, tlx[i], tly[i])
end

##----------- aCP
tlx = ["" "" L"a_{CP}"]
for i in 1:3
    subplot(3, 4, (i-1) * 4 + 2)
    plot_sa_unit(rg_aCP, res_aCP, ind[i], false, tlx[i], "")
end

##----------- eCP
tlx = ["" "" L"e_{CP}"]
for i in 1:3
    subplot(3, 4, (i-1) * 4 + 3)
    plot_sa_unit(rg_eCP, res_eCP, ind[i], false, tlx[i], "")
end

##----------- mP
tlx = ["" "" L"m_P"]
for i in 1:3
    subplot(3, 4, (i-1) * 4 + 4)
    plot_sa_unit_rev(rg_mP, res_mP, ind[i], i == 1, tlx[i], "")
    val = [0.15, 0.20, 0.25, 0.30]
    xticks(val, reverse(val))
end

tight_layout()
savefig("fig/fig_S2.svg")

max_stars_repo_name: McCannLab/Labmnivory

max_issues_repo_name: McCannLab/Labmnivory

max_forks_repo_name: McCannLab/Labmnivory


text:
function DistributedSeqModel(
    # constants
    num_partitions::Int64,
    max_time::Float64,
    max_sequence_length::Float64,

    # warp parameters
    num_warp_values::Int64,
    max_warp::Float64,
    warp_variance::Float64,

    # priors
    seq_event_rate::Float64,
    seq_type_proportions::SymmetricDirichlet,
    seq_event_amplitude::RateGamma,
    neuron_response_proportions::SymmetricDirichlet,
    neuron_response_profile::NormalInvChisq,
    bkgd_amplitude::RateGamma,
    bkgd_proportions::SymmetricDirichlet
)

    primary_model = SeqModel(
        max_time,
        max_sequence_length,
        num_warp_values,
        max_warp,
        warp_variance,
        seq_event_rate,
        seq_type_proportions,
        seq_event_amplitude,
        neuron_response_proportions,
        neuron_response_profile,
        bkgd_amplitude,
        bkgd_proportions
    )

    return DistributedSeqModel(
        primary_model,
        num_partitions
    )
end

function DistributedSeqModel(primary_model::SeqModel, num_partitions::Int64)

    # Number of neurons (N), sequence types (R), and warps (W).
    N = num_neurons(primary_model)
    R = num_sequence_types(primary_model)
    W = num_warp_values(primary_model)

    # Even though the submodels are each allocated a time interval of
    # (max_time / num_partitions), we set "max_time", "bkgd_log_prob",
    # and "new_cluster_log_prob" to be the same as the primary model.
    # Thus, all the probability calculations within the submodels match
    # the primary model.

    submodels = SeqModel[]
    for part = 1:num_partitions
        push!(
            submodels,
            SeqModel(
                primary_model.max_time,
                primary_model.max_sequence_length,

                primary_model.priors,
                primary_model.globals,
                SeqEventList(R, W),
                primary_model.new_cluster_log_prob,
                primary_model.bkgd_log_prob,
                zeros(R),    # _R_buffer
                zeros(R, W), # _RW_buffer
                Float64[],   # _K_buffer
            )
        )
    end

    return DistributedSeqModel(
        primary_model,
        num_partitions,
        submodels
    )
end

num_neurons(model::DistributedSeqModel) =
    model.primary_model.priors.bkgd_proportions.dim
num_sequence_types(model::DistributedSeqModel) =
    model.primary_model.priors.seq_type_proportions.dim
num_sequence_events(model::DistributedSeqModel) =
    length(model.primary_model.sequence_events)

log_prior(model::DistributedSeqModel) =
    log_prior(model.primary_model)

log_p_latents(model::DistributedSeqModel) =
    log_p_latents(model.primary_model)

log_like(model::DistributedSeqModel, spikes::Vector{Spike}) =
    log_like(model.primary_model, spikes)

log_joint(model::DistributedSeqModel, spikes::Vector{Spike}) =
    log_joint(model.primary_model, spikes)

sample(model::DistributedSeqModel; kwargs...) = sample(model.primary_model; kwargs...)

priors(model::DistributedSeqModel) =
    return model.primary_model.priors

function set_priors!(model::DistributedSeqModel, priors::SeqPriors)
    model.primary_model.priors = priors
    for submodel in model.submodels
        submodel.priors = priors
    end
end

function set_new_cluster_log_prob!(model::DistributedSeqModel, prob::Float64)
    # Since we constructed each submodel.max_time == primary_model.max_time, the
    # probability of forming new clusters is the same.
    model.primary_model.new_cluster_log_prob = prob
    for submodel in model.submodels
        submodel.new_cluster_log_prob = prob
    end
end

get_max_time(model::DistributedSeqModel) = model.primary_model.max_time

get_bkgd_amplitude(model::DistributedSeqModel) = model.primary_model.globals.bkgd_amplitude

function set_bkgd_log_prob!(model::DistributedSeqModel, prob::Float64)
    # Since we constructed each submodel.max_time == primary_model.max_time, the
    # probability of attributing spikes to the background is the same.
    model.primary_model.bkgd_log_prob = prob
    for submodel in model.submodels
        submodel.bkgd_log_prob = prob
    end
end


text:
#=
# Stress on Truss
In this example we'll plot an animation of the stress on some truss structure
using `GaphMakie.jl` and [`NetworkDynamics.jl`](https://github.com/PIK-ICoN/NetworkDynamics.jl)

![truss animation](truss.mp4)
=#
using NetworkDynamics
using OrdinaryDiffEq
using Graphs
using GraphMakie
using LinearAlgebra: norm, ⋅
using Printf
using CairoMakie
CairoMakie.activate!()

#=
## Definition of the dynamical system
Define dynamic model for NetworkDynamics
=#
function edge_f!(e, v_s, v_d, (K, L), _)
    Δd = view(v_d, 1:2) - view(v_s, 1:2)
    nd = norm(Δd)
    F = K * ( L - nd )
    e[1:2] = F .* view(Δd, 1:2) ./ nd
end

function vertex_free!(du, u, edges, (M, γ, g), _)
    F = sum(edges) - γ * view(u, 3:4) + [0, -M*g]
    du[3:4] .= F ./ M
    du[1:2] = view(u, 3:4)
end

function vertex_fixed!(du, u, edges, _, _)
    du[1:2] .= 0.0
end

edge = StaticEdge(f = edge_f!, dim=2, coupling=:antisymmetric)
vertex_free = ODEVertex(f = vertex_free!, dim=4, sym=[:x, :y, :v, :w])
vertex_fix  = ODEVertex(f = vertex_fixed!, dim=2, sym=[:x, :y, :v, :w])
nothing #hide

#=
Set up graph topology and initial positions.
=#
N = 5
dx = 1.0
shift = 0.2
g = SimpleGraph(2*N + 1)
for i in 1:N
    add_edge!(g, i, i+N); add_edge!(g, i, i+N)
    if i < N
        add_edge!(g, i+1, i+N); add_edge!(g, i, i+1); add_edge!(g, i+N, i+N+1)
    end
end
add_edge!(g, 2N, 2N+1)
pos0 = zeros(Point2f, 2N + 1)
pos0[1:N] = [Point((i-1)dx,0) for i in 1:N]
pos0[N+1:2*N] = [Point(i*dx + shift, 1) for i in 1:N]
pos0[2N+1] = Point(N*dx + 1, -1)

fixed = [1,4] # set fixed vertices
nothing #hide

# create networkdynamics object
vertices = ODEVertex[vertex_free for i in 1:nv(g)]
for i in fixed
    vertices[i] = vertex_fix # use the fixed vertex for the fixed
end
nd = network_dynamics(vertices, edge, g);
nothing #hide

#=
write some auxiliary functions to map between state vector of DGL `u` and graph data
=#
x_idxs = [findfirst(isequal(Symbol("x_$i")), nd.syms) for i in 1:nv(g)]

"Set positions `pos` inside dgl state `u`"
function set_pos!(u, pos)
    for (i, idx) in enumerate(x_idxs)
        u[idx] = pos[i][1]
        u[idx+1] = pos[i][2]
    end
end

"Extract vector of Points from dgl state `u`"
function to_pos(u)
    pos = Vector{Point2f}(undef, nv(g))
    for (i, idx) in enumerate(x_idxs)
        pos[i] = Point(u[idx], u[idx+1])
    end
    return pos
end

"Calculate load on edge for given state."
function get_load(u, p, t=0.0)
    gd_nd = nd(u, p, t, GetGD) # exposes underlying graph data struct
    force = Vector{Float64}(undef, ne(g))
    pos = to_pos(u)
    for (i,e) in enumerate(edges(g))
        edgeval = get_edge(gd_nd, i)
        fvec = Point(edgeval[1], edgeval[2])
        dir = pos[e.dst] .- pos[e.src]
        force[i] = sign(fvec ⋅ dir) * norm(fvec)
    end
    return force
end
nothing #hide

#=
Set parameters.
=#
M = [10 for i in 1:nv(g)] # mass of the nodes
M[end] = 200 # heavy mass
gc = [9.81 for i in 1:nv(g)] # gravitational constant
γ = [200.0 for i in 1:nv(g)] # damping parameter
γ[end] = 100.0

L = [norm(pos0[e.src] - pos0[e.dst]) for e in edges(g)] # length of edges
K = [0.5e6 for i in 1:ne(g)] # spring constant of edges

## bundle parameters for NetworkDynamics
para = collect.((zip(M, γ, gc),
                 zip(K, L)))
nothing #hide

#=
Set initial state and solve the system
=#
u0 = zeros(length(nd.syms))
set_pos!(u0, pos0)

tspan = (0.0, 12.0)
prob = ODEProblem(nd, u0, tspan, para)
sol  = solve(prob, Tsit5());
nothing #hide

#=
## Plot the solution
=#

fig = Figure(resolution=(1000,550))
fig[1,1] = title = Label(fig, "Stress on truss", textsize=30)
title.tellwidth = false

fig[2,1] = ax = Axis(fig)

## calculate some values for colorscaling
(fmin, fmax) = 0.3 .* extrema([(map(u->get_load(u, para), sol)...)...])
nlabels = [" " for i in 1:nv(g)]
nlabels[fixed] .= "Δ"
elabels = ["edge $i" for i in 1:ne(g)]
p = graphplot!(ax, g;
               edge_width=4.0,
               node_size=sqrt.(M)*3,
               nlabels=nlabels,
               nlabels_align=(:center,:top),
               nlabels_textsize=30,
               elabels=elabels,
               elabels_opposite=[ne(g)],
               edge_color=[0.0 for i in 1:ne(g)],
               edge_attr=(colorrange=(fmin,fmax),
                          colormap=:diverging_bkr_55_10_c35_n256))
hidespines!(ax); hidedecorations!(ax); p[:node_pos][]=to_pos(u0); ax.aspect = DataAspect()
limits!(ax, -0.1, pos0[end][1]+0.3, pos0[end][2]-0.5, 1.15)

## draw colorbar
fig[3,1] = cb = Colorbar(fig, p.plots[1], label = "Axial force", vertical=false)

T = tspan[2]
fps = 30
trange = range(0.0, sol.t[end], length=Int(T * fps))
record(fig, "truss.mp4", trange; framerate=fps) do t
    title.text = @sprintf "Stress on truss (t = %.2f )" t
    p[:node_pos][] = to_pos(sol(t))
    load = get_load(sol(t), para)
    p.elabels = [@sprintf("%.0f", l) for l in load]
    p.edge_color[] = load
end
nothing #hide


text:
struct NPlayerUnicycleCost{TIC,TICR,TSC,TSCR,TPC,TGC}<:NPlayerNavigationCost
    "the index of the player this cost applies to"
    player_id::Int
    inputcost::TIC
    inputconstr::TICR
    statecost::TSC
    stateconstr::TSCR
    proximitycost::TPC
    goalcost::TGC
end

function NPlayerUnicycleCost(player_id, xg, t_final, np;
         inputcost::TIC=QuadCost(SMatrix{2,2}([1. 0.; 0. 1.]) * 10),
         inputconstr::TICR=(SoftConstr(1, deg2rad(-10), deg2rad(10), 50),
                            SoftConstr(2, -9.81, 9.81, 50)),
         statecost::TSC=QuadCost(SMatrix{4,4}(diagm([0., 0., 0., 1.])) * 30),
         stateconstr::TSCR=(SoftConstr(4, -0.05, 2., 50),),
         proximitycost::TPC=ProximityCost(2.0, 50.0, np),
         goalcost::TGC=GoalCost(t_final, xg, SMatrix{4,4}(diagm([1.,1.,1.,0.])) *
                                300)) where {TXI,TUI,TIC,TICR,TSC,TSCR,TPC,TGC}

    return NPlayerUnicycleCost(player_id, inputcost, inputconstr, statecost,
                               stateconstr, proximitycost, goalcost)
end

"-----------------Implementing the NPlayerNavigationCost interface-----------------"

player_id(c::NPlayerUnicycleCost) = c.player_id
inputcost(c::NPlayerUnicycleCost) = c.inputcost
inputconstr(c::NPlayerUnicycleCost) = c.inputconstr
statecost(c::NPlayerUnicycleCost) = c.statecost
stateconstr(c::NPlayerUnicycleCost) = c.stateconstr
proximitycost(c::NPlayerUnicycleCost) = c.proximitycost
goalcost(c::NPlayerUnicycleCost) = c.goalcost
ravoid(g) = mapreduce(c->minimum(proximitycost(c).rs), min, player_costs(g))

"------------------ Implementing Feedback Linearization Interface -----------------"

"Here, we don't transform the cost explicitly but as suggested in
https://arxiv.org/abs/1910.00681 we formulate a new cost directly in ξ coordinates."
function transformed_cost(cs::Unicycle4D, c::NPlayerUnicycleCost, np::Int)
    @unpack t_active, xg = goalcost(c)
    ξg = ξ_from(cs, xg)
    if λ_issingular(cs, ξg)
        @warn "State conversion map is singular at provided goal state."
    end
    return NPlayer2DDoubleIntegratorCost(player_id(c), ξg, t_active, np)
end


text:
using Celeste: Model, SensitiveFloats, DeterministicVI
using DeterministicVI: KLDivergence
using Distributions, DiffBase
using Base.Test

"""
Use Monte Carlo to check whether KL(q_dist || p_dist) matches exact_kl

Args:
    q_dist, p_dist: Distribution objects
    exact_kl: The expected exact KL
"""
function test_kl_value(q_dist, p_dist, exact_kl::Float64)
    sample_size = 2_000_000
    q_samples = rand(q_dist, sample_size)
    empirical_kl_samples = logpdf(q_dist, q_samples) - logpdf(p_dist, q_samples)
    empirical_kl = mean(empirical_kl_samples)
    tol = 4 * std(empirical_kl_samples) / sqrt(sample_size)
    min_diff = 1e-2 * std(empirical_kl_samples) / sqrt(sample_size)
    @test isapprox(empirical_kl, exact_kl, atol=tol)
end

function test_beta_kl_value()
    alpha1, beta1 = 4.1, 3.9
    alpha2, beta2 = 3.5, 4.3
    p1_dist = Beta(alpha1, beta1)
    p2_dist = Beta(alpha2, beta2)
    kl = KLDivergence.beta_kl(alpha1, beta1, alpha2, beta2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_categorical_kl_value()
    p1 = Float64[1, 2, 3, 4]
    p2 = Float64[5, 6, 2, 1]
    p1 = p1 ./ sum(p1)
    p2 = p2 ./ sum(p2)
    p1_dist = Categorical(p1)
    p2_dist = Categorical(p2)
    kl = KLDivergence.categorical_kl(p1, p2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_diagmvn_mvn_kl_value()
    K = 4
    mean1 = rand(K)
    var1 = rand(K)
    var1 = var1 .* var1
    mean2 = rand(K)
    cov2 = rand(K, K)
    cov2 = 0.2 * cov2 * cov2' + eye(K)
    cov2 = 0.5 * (cov2 + cov2')
    p1_dist = MvNormal(mean1, diagm(var1))
    p2_dist = MvNormal(mean2, cov2)
    kl = KLDivergence.diagmvn_mvn_kl(mean1, var1, mean2, cov2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_gaussian_kl_value()
    mean1, var1 = 0.5, 2.0
    mean2, var2 = 0.8, 1.8
    p1_dist = Normal(mean1, sqrt(var1))
    p2_dist = Normal(mean2, sqrt(var2))
    kl = KLDivergence.gaussian_kl(mean1, var1, mean2, var2)
    test_kl_value(p1_dist, p2_dist, kl)
end

function test_subtract_kl()
    sf = SensitiveFloat{Float64}(KLDivergence.PARAM_LENGTH, 1, true, true)
    vs = rand(MersenneTwister(1), KLDivergence.PARAM_LENGTH)
    kl_result = DiffBase.DiffResult(0.0, sf.d)
    kl_helper = KLDivergence.get_kl_helper(Float64)
    KLDivergence.subtract_kl_source!(sf, kl_result, vs, kl_helper)
end

println("Running KL divergence tests.")
test_beta_kl_value()
test_categorical_kl_value()
test_diagmvn_mvn_kl_value()
test_gaussian_kl_value()
test_subtract_kl()


text:
using ConstrainedDynamics

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1
box = Box(width, depth, length1, length1)

p2 = [0.0;0.0;length1 / 2] # joint connection point

# Initial orientation
ϕ1 = π / 2
q1 = UnitQuaternion(RotX(ϕ1))

# Links
origin = Origin{Float64}()
link1 = Body(box)

# Constraints
joint_between_origin_and_link1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=p2))

links = [link1]
constraints = [joint_between_origin_and_link1]
shapes = [box]

mech = Mechanism(origin, links, constraints, shapes = shapes)
setPosition!(origin,link1,p2 = p2,Δq = q1)


text:
function fullneighborlist(x, a, b, c, pbc, rcutmax)

latticemethod = 0
if (pbc[1] == 1) & (rcutmax >= 1.0*a[1])
    latticemethod = 1
end
if (pbc[2] == 1) & (rcutmax >= 1.0*b[2])
    latticemethod = 1
end
if (length(c)>2) & (length(pbc) > 2)
    if (pbc[3] == 1) & (rcutmax >= 1.0*c[3])
        latticemethod = 1
    end
end

if (latticemethod == 1)
    y, alist, neighlist, neighnum = latticecoords(x, rcutmax, pbc, a, b, c)
    return y, alist, neighlist, neighnum
end

# offset simulation box for periodic boundary conditions
boxoffset = [rcutmax; rcutmax; rcutmax];

# bounding box
B2R, R2B, v, w, vr, wr, pimages = prebox(pbc[:], boxoffset, a[:], b[:], c[:]);

y, alist = atomlist(x, pimages, wr, B2R);

neighi, neighj, neighnum = neighborlist(y, rcutmax);

n = size(x,2);
neighi = neighi[neighi .<= n];
ne = length(neighi);
neighlist = neighj[1:ne];
neighnum = neighnum[1:n];
neighnum = [0; cumsum(neighnum)];

return y, alist, neighlist, neighnum

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
# example, following here

inline_png2 = b"\x89PNG\r\n\x1a\n\x00\x00\x00\rIHDR\x00\x00\x00d\x00\x00\x00d\x08\x02\x00\x00\x00\xff\x80\x02\x03\x00\x00\x00\x06bKGD\x00\xff\x00\xff\x00\xff\xa0\xbd\xa7\x93\x00\x00\x01=IDATx\x9c\xed\xdc\xc1i\x85@\x14@\xd1IH!\x96bi\x96f)\x96\x92E@!\x01\x7f\xeeF\xfdp\xce\xca\x9d\xc3\xe5\xb9{\xce\xc7\xb2,\x83\xff\xf9\xbc\xfb\x00\xefD\xac@\xac@\xac@\xac@\xac@\xac@\xac@\xac\xe0\xeb\xee\x03\xfc\xb1\xaec]\xc7\x18c\x9e\xc7<\xdf{\x96_LV V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V V\xf0j[\xf9gq\xf8J\xdbv<\\\xff\xf6\xd3\xfd\xe8\xe7\xc5\xdam\xdb\x11\xee2\xa7\xb1|\x86\xc1\xab\xc9\xba~m\x7f\x1f\xa8i\x1a\xd3t\xf5\xdbO=/\xd6\xba\x1e\xb1\xfca\xf1\xbe\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\xc4\n\x9ew\x1f\xfc\xf3\xae\x81\xdf\x99\xac@\xac@\xac@\xac@\xac@\xac@\xac@\xac@\xac\xe0\x1b\xaf\xef\x1c\xd4\x17\xfe\xcaX\x00\x00\x00\x00IEND\xaeB`\x82"
png_stream = IOBuffer();
write(png_stream,inline_png2);
seekstart(png_stream);

imagesurface = read_from_png(png_stream);
set_source_surface(cr, imagesurface, 0, 0);
paint(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_imagestream.png");


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
using SparseArrays
Random.seed!(233)

function sparse_to_dense(A)
    A = constant(A)
    ij = A.o.indices
    vv = values(A)
    m, n = size(A)
    sparse_to_dense_ = load_op_and_grad("./build/libSparseToDense","sparse_to_dense_ad")
    m_, n_ = convert_to_tensor(Any[m,n], [Int64,Int64])
    out = sparse_to_dense_(ij, vv, m_,n_)
    set_shape(out, (m, n))
end

A = sprand(10,10,0.3)
# TODO: specify your input parameters
u = sparse_to_dense(A)
sess = Session(); init(sess)
@show run(sess, u) - Array(A)

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    B = SparseTensor(ij[:,1]+1, ij[:,2]+1, x, m, n)
    return sum(sparse_to_dense(B)^2)
end

# TODO: change `m_` and `v_` to appropriate values
A = constant(A)
m, n = size(A)
K = length(values(A))
ij = A.o.indices
m_ = constant(rand(K))
v_ = rand(K)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
"""
*AbstractTask*

An abstract type for a task. A task requires the following fields:

- `visible`: shows GUI if true
- `realtime`: executes model in realtime if true
- `speed`: speed of realtime model execution
- `scheduler`: a reference to the event scheduler

A task may optionally contain the following fields:

- `screen`: a vector of visual objects on the screen
- `canvas`: a GUI object component
- `window`: a window for the GUI

*Example*

The following is an example of the task object for the PVT.

```julia
mutable struct PVT{T,W,C,F1,F2} <: AbstractTask
    n_trials::Int
    trial::Int
    lb::Float64
    ub::Float64
    width::Float64
    hight::Float64
    scheduler::T
    screen::Vector{VisualObject}
    canvas::C
    window::W
    visible::Bool
    realtime::Bool
    speed::Float64
end
```
The constructor for `PVT` is

```julia
function PVT(;
    n_trials=10,
    trial=1,
    lb=2.0,
    ub=10.0,
    width=600.0,
    height=600.0,
    scheduler=nothing,
    screen=Vector{VisualObject}(),
    window=nothing,
    canvas=nothing,
    visible=false,
    realtime=false,
    speed=1.0,
    )
    visible ? ((canvas,window) = setup_window(width)) : nothing
    visible ? Gtk.showall(window) : nothing
    return PVT(n_trials, trial, lb, ub, width, height, scheduler, screen, canvas, window, visible,
        realtime, speed)
end

function setup_window(width)
	canvas = @GtkCanvas()
    window = GtkWindow(canvas, "PVT", width, width)
    Gtk.visible(window, true)
    @guarded draw(canvas) do widget
        ctx = getgc(canvas)
        rectangle(ctx, 0.0, 0.0, width, width)
        set_source_rgb(ctx, .8, .8, .8)
        fill(ctx)
    end
	return canvas,window
end
```
"""
abstract type AbstractTask end

"""
    ACTRScheduler <: AbstractScheduler

An ACT-R event scheduler object.

# Fields

- `events`: a priority queue of events
- `time`: current time of the system
- `running`: simulation can run if true
- `model_trace`: will print out model events if true
- `task_trace`: will print out task events if true
- `store`: will store a vector of completed events if true
- `complete_events`: an optional vector of completed events

"""
mutable struct ACTRScheduler{PQ<:PriorityQueue,E} <: AbstractScheduler
    events::PQ
    time::Float64
    running::Bool
    model_trace::Bool
    task_trace::Bool
    store::Bool
    complete_events::E
end

"""
    ACTRScheduler(;event=Event, time=0.0, running=true, model_trace=false, task_trace=false, store=false)

Constructor for Scheduler with default keyword values:

# Keywords

- `event`: a function that is executed at the specified time
- `time`: the time at which the event will be execute
- `running`: whether the model is running
- `model_trace`: prints model trace if true
- `task_trace`: prints task trace if true
- `store`: stores the executed events for replay if set to true
"""
function ACTRScheduler(;event=Event, time=0.0, running=true, model_trace=false, task_trace=false, store=false)
    events = PriorityQueue{event,Float64}()
    return ACTRScheduler(events, time, running, model_trace, task_trace, store, Vector{event}())
end

run!(actr::AbstractACTR, task::AbstractTask, until=Inf) = run!([actr], task, until)

"""
    run!(actr, task::AbstractTask, until=Inf)

Simulate an ACT-R model

# Arguments

- `models`: a dictionary of ACT-R model objects
- `task`: a task that is a subtype of `AbstractTask`
- `until=Inf`: a specified termination time unless terminated sooner manually
"""
function run!(models, task::AbstractTask, until=Inf)
    s = task.scheduler
    last_event!(s, until)
    start!(task, models)
    start!.(models)
    fire!.(models)
    while is_running(s, until)
        next_event!(s, models, task)
    end
    s.task_trace && !s.running ? print_event(s.time, "", "Stopped") : nothing
    return nothing
end

next_event!(actr, task) = next_event!(actr.scheduler, actr, task)

"""
    next_event!(s, actr, task)

Dequeue and process a single event.

# Arguments

- `s`: an event scheduler
- `actr`: an ACT-R model object
- `task`: a task that is a subtype of `AbstractTask`
"""
function next_event!(s, models, task)
    event = dequeue!(s.events)
    pause(task, event)
    s.time = event.time
    event.fun()
    fire!.(models)
    s.store ? push!(s.complete_events, event) : nothing
    s.model_trace && event.type == "model" ? print_event(event) : nothing
    s.task_trace && event.type !== "model" ? print_event(event) : nothing
    return nothing
end

"""
    pause(task, event)

Pauses simulation for specified `speed` if `realtime` in `task` is true.

# Arguments

- `task`: a task that is a subtype of `AbstractTask`
- `event`: a task event or an internal event of the model
"""
function pause(task, event)
    !task.realtime ? (return nothing) : nothing
    t = (event.time - task.scheduler.time) / task.speed
    sleep(t)
    return nothing
end

"""
    start!(actr)

A function that initializes the simulation for the model.

# Arguments

- `task`: a task that is a subtype of `AbstractTask`
- `actr`: an ACT-R model object
"""
function start!(actr::AbstractACTR)
    register!(actr, ()->(), now; description="Starting", id=get_name(actr))
    return nothing
end


text:
tomlpath = joinpath(@__DIR__, "flextopo_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_flextopo_model(config)
@unpack network = model

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    row = csv_first_row(model.writer.csv_path)

    @test row.time == DateTime("2010-01-01T00:00:00")
    @test row.Q_16 ≈ 0.0000370817920883859f0
    @test row.percentageH_16 ≈ 0.340213f0
    @test row.QfP_1011 ≈ 0.0f0
    @test row.kfW_10 ≈ 0.136334478855133f0
    @test row.kfH_10 ≈ 0.0454448275268077f0
    @test row.Qs_503 ≈ 0.1575031550601124f0
    @test row.Si_16 ≈ 0.0f0
    @test row.Ss_16 ≈ 29.89794354323524f0
    @test row.Ea_16 ≈ 0.0110627892408438f0
end

@testset "first timestep" begin
    flextopo = model.vertical
    @test flextopo.tt[3500] ≈ 1.3f0
    @test model.clock.iteration == 2
    @test flextopo.rootzonestorage[3500] ≈ [147.11238663084805f0, 79.08369375691255f0, 79.23637697443984f0]
    @test flextopo.runoff[3500] ≈ 0.19008129369467497f0
    @test flextopo.rootevap[3500] ≈ [0.009225917980074883f0, 0.009225917980074883f0, 0.009225917980074883f0]
    @test flextopo.snow[3500] == 0.0
end

# run the second timestep
Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "second timestep" begin
    flextopo = model.vertical
    @test flextopo.rootzonestorage[3500] ≈ [146.73222385533154f0, 78.55190222246516f0, 78.85739340140256f0]
    @test flextopo.runoff[3500] ≈ 0.18887692333975817f0
    @test flextopo.rootevap[3500] ≈ [0.38016277551651f0, 0.38016277551651f0, 0.38016277551651f0]
    @test flextopo.snow[3500] == 0.0f0
end

@testset "overland domain" begin
    q = model.lateral.land.q_av
    land = model.lateral.land
    @test sum(q) ≈ 202.1162497378859f0
    @test q[10354] ≈ 0.001624108000162103f0
    @test land.volume[10354] ≈ 66.41398729880544f0
    @test land.inwater[10354] ≈ 0.001658182220814705f0
    @test q[network.land.order[end]] ≈ 0.0033337667005815565f0
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 108.45732176546284f0
    @test q[651] ≈ 0.014755293483684526f0
    @test q[1056] ≈ 0.0030933658056867823f0
    @test q[network.river.order[end]] ≈ 0.0005541137960384742f0
end

#change the config to use other functions for the storages for several classes:
config = Wflow.Config(tomlpath)
config["model"]["select_snow"] = ["common_snow_no_storage"]
config["model"]["select_interception"] = ["interception_no_storage", "interception_overflow", "interception_overflow"]
config["model"]["select_hortonponding"] = ["hortonponding_no_storage", "hortonponding", "hortonponding_no_storage"]
config["model"]["select_hortonrunoff"] = ["hortonrunoff_no_storage", "hortonrunoff", "hortonrunoff_no_storage"]
config["model"]["select_rootzone"] = ["rootzone_storage", "rootzone_storage", "rootzone_no_storage"]
config["model"]["select_fast"] = ["fast_no_storage", "fast_storage", "fast_storage"]
config["model"]["select_slow"] = ["common_slow_storage"]

model = Wflow.initialize_flextopo_model(config)
@unpack network = model

Wflow.load_dynamic_input!(model)
model = Wflow.update(model)

@testset "first timestep" begin
    flextopo = model.vertical
    @test flextopo.rootzonestorage[3500] ≈ [147.11238663084805f0, 79.08369375691255f0, 0.0f0]
    @test flextopo.runoff[3500] ≈ 0.19008129369467497f0
    @test flextopo.rootevap[3500] ≈ [0.009225917980074883f0, 0.009225917980074883f0, 0.0f0]
    @test flextopo.snow[3500] == 0.0f0
end

@testset "overland domain" begin
    q = model.lateral.land.q_av
    land = model.lateral.land
    @test sum(q) ≈ 69.51434665542678f0
    @test q[10354] ≈ 0.0009376627801121855f0
    @test land.volume[10354] ≈ 63.46997463442048f0
    @test land.inwater[10354] ≈ 0.0016722689680105701f0
    @test q[network.land.order[end]] ≈ 0.0006175366185364747f0
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 22.236556345676156f0
    @test q[651] ≈ 0.00035906990499781237f0
    @test q[1056] ≈ 0.0002647264689318855f0
    @test q[network.river.order[end]] ≈ 0.000015647838277869928f0
end

Wflow.close_files(model, delete_output = false)


text:
module Hooks

export @hook_str, handle, unhandle

struct HookName{S} end

mutable struct Hook
    handlers::Set{Any}
    hasrun::Bool
    args::Tuple
end

@inline Hook(args...) = Hook(Set{Any}(args), false, ())

"""
    hook"some-event"

String macro to represent the name of a hook.
"""
macro hook_str(name)
    :(HookName{Symbol($name)}())
end

hooks = Dict{Symbol, Hook}()

function Base.notify(::HookName{S}, args...) where S
    hook = get(hooks, S, nothing)
    if hook === nothing
        hook = Hook()
        hooks[S] = hook
    end

    for h in hook.handlers
        h(args...)
    end

    hook.hasrun = true
    hook.args = args

    nothing
end

function handle(handler, ::HookName{S}) where S
    if S in keys(hooks)
        hook = hooks[S]
        push!(hook.handlers, handler)
        if hook.hasrun
            handler(hook.args...)
        end
    else
        hooks[S] = Hook(handler)
    end

    handler
end

function unhandle(handler, ::HookName{S}) where S
    hook = get(hooks, S, nothing)
    if hook === nothing
        throw(ErrorException("Hook $S not found"))
    end
    if handler in hook.handlers
        delete!(hook.handlers, handler)
    else
        throw(ErrorException("Hook $S has no handler $handler"))
    end

    nothing
end

function Base.reset(::HookName{S}) where S
    hook = get(hooks, S, nothing)
    if hook === nothing
        throw(ErrorException("Hook $S not found"))
    end
    hook.handlers = Set()
    hook.hasrun = false
    hook.args = ()

    nothing
end

end # module


text:
#####
##### Outer functions for filling halo regions for open boundary conditions.
#####

# TODO: support true open boundary conditions.
# For this we need to have separate functions for each of the six boundaries,
# and need to unroll a loop over the boundary normal direction.
# The syntax for `getbc` is also different for OpenBoundaryCondition than for others,
# because the boundary-normal index can vary (and array boundary conditions need to be
# 3D in general).

#=
@kernel function set_west_or_east_u!(u, i_boundary, bc, loc, grid, args...)
    j, k = @index(Global, NTuple)
    @inbounds u[i_boundary, j, k] = getbc(bc, j, k, grid, args...)
end

@kernel function set_south_or_north_v!(v, j_boundary, bc, loc, grid, args...)
    i, k = @index(Global, NTuple)
    @inbounds v[i, j_boundary, k] = getbc(bc, i, k, grid, args...)
end

@kernel function set_bottom_or_top_w!(w, k_boundary, bc, loc, grid, args...)
    i, j = @index(Global, NTuple)
    @inbounds w[i, j, k_boundary] = getbc(bc, i, j, grid, args...)
end

  @inline fill_west_halo!(u, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :yz, set_west_or_east_u!,   u,           1, bc, grid, args...; dependencies=dep, kwargs...)
  @inline fill_east_halo!(u, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :yz, set_west_or_east_u!,   u, grid.Nx + 1, bc, grid, args...; dependencies=dep, kwargs...)
 @inline fill_south_halo!(v, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xz, set_south_or_north_v!, v,           1, bc, grid, args...; dependencies=dep, kwargs...)
 @inline fill_north_halo!(v, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xz, set_south_or_north_v!, v, grid.Ny + 1, bc, grid, args...; dependencies=dep, kwargs...)
@inline fill_bottom_halo!(w, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xy, set_bottom_or_top_w!,  w,           1, bc, grid, args...; dependencies=dep, kwargs...)
   @inline fill_top_halo!(w, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xy, set_bottom_or_top_w!,  w, grid.Nz + 1, bc, grid, args...; dependencies=dep, kwargs...)
=#

@inline   _fill_west_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[1, j, k]           = getbc(bc, j, k, grid, args...)
@inline   _fill_east_halo!(j, k, grid, c, bc::OBC, loc, args...) = @inbounds c[grid.Nx + 1, j, k] = getbc(bc, j, k, grid, args...)
@inline  _fill_south_halo!(i, k, grid, c, bc::OBC, loc, args...) = @inbounds c[i, 1, k]           = getbc(bc, i, k, grid, args...)
@inline  _fill_north_halo!(i, k, grid, c, bc::OBC, loc, args...) = @inbounds c[i, grid.Ny + 1, k] = getbc(bc, i, k, grid, args...)
@inline _fill_bottom_halo!(i, j, grid, c, bc::OBC, loc, args...) = @inbounds c[i, j, 1]           = getbc(bc, i, j, grid, args...)
@inline    _fill_top_halo!(i, j, grid, c, bc::OBC, loc, args...) = @inbounds c[i, j, grid.Nz + 1] = getbc(bc, i, j, grid, args...)


text:
import LibGEOS
import GeoInterface

function squares()
    A_T = LibGEOS.readgeom("POLYGON((1 1,3 1,3 3,1 3,1 1))")
    A_C = LibGEOS.readgeom("POLYGON((3 1,5 1,5 3,3 3,3 1))")
    border = GeoRDD.get_border(A_T, A_C, 0.0)
    @test GeoInterface.coordinates(border) == [[3.0, 1.0], [3.0, 3.0]]
    border = GeoRDD.get_border(A_T, A_C, 1.0)
    @test GeoInterface.coordinates(border) == [[2.0, 1.0], [3.0, 1.0], [3.0, 3.0], [2.0, 3.0]]
end

function twosquaresT()
    A_T1 = LibGEOS.readgeom("POLYGON((3 1,5 1,5 3,3 3,3 1))")
    A_T2 = LibGEOS.readgeom("POLYGON((1 3,3 3,3 5,1 5,1 3))")
    A_T = LibGEOS.union(A_T1, A_T2)
    A_C = LibGEOS.readgeom("POLYGON((1 1,3 1,3 3,1 3,1 1))")
    border = GeoRDD.get_border(A_T, A_C, 0.0)
    @test GeoInterface.coordinates(border) == [[3.0, 1.0], [3.0, 3.0], [1.0, 3.0]]
    border = GeoRDD.get_border(A_T, A_C, 1.0)
    # this isn't ideal:
    @test GeoInterface.coordinates(border) == [[[4.0, 1.0], [3.0, 1.0], [3.0, 3.0], [4.0, 3.0]],
                                               [[3.0, 4.0], [3.0, 3.0], [1.0, 3.0], [1.0, 4.0]]]
end

function twosquaresC()
    A_C1 = LibGEOS.readgeom("POLYGON((3 1,5 1,5 3,3 3,3 1))")
    A_C2 = LibGEOS.readgeom("POLYGON((1 3,3 3,3 5,1 5,1 3))")
    A_C = LibGEOS.union(A_C1, A_C2)
    A_T = LibGEOS.readgeom("POLYGON((1 1,3 1,3 3,1 3,1 1))")
    border = GeoRDD.get_border(A_T, A_C, 0.0)
    @test GeoInterface.coordinates(border) == [[3.0, 1.0], [3.0, 3.0], [1.0, 3.0]]
    border = GeoRDD.get_border(A_T, A_C, 1.0)
    # this one is better...
    @test GeoInterface.coordinates(border) == [[2.0, 1.0], [3.0, 1.0], [3.0, 3.0], [1.0, 3.0], [1.0, 2.0]]
end

function twosquaresC_clock()
    A_C1 = LibGEOS.readgeom("POLYGON((3 1,5 1,5 3,3 3,3 1))")
    A_C2 = LibGEOS.readgeom("POLYGON((1 3,1 5,3 5,3 3,1 3))")
    A_C = LibGEOS.union(A_C1, A_C2)
    A_T = LibGEOS.readgeom("POLYGON((1 1,3 1,3 3,1 3,1 1))")
    border = GeoRDD.get_border(A_T, A_C, 0.0)
    @test GeoInterface.coordinates(border) == [[3.0, 1.0], [3.0, 3.0], [1.0, 3.0]]
    border = GeoRDD.get_border(A_T, A_C, 1.0)
    # this one is better...
    @test GeoInterface.coordinates(border) == [[2.0, 1.0], [3.0, 1.0], [3.0, 3.0], [1.0, 3.0], [1.0, 2.0]]
end

function wrapsquares()
    A_T = LibGEOS.readgeom("POLYGON((3 1,5 1,5 5,1 5,1 3,3 3,3 1))")
    A_C = LibGEOS.readgeom("POLYGON((1 1,3 1,3 3,1 3,1 1))")
    border = GeoRDD.get_border(A_T, A_C, 0.0)
    @test GeoInterface.coordinates(border) == [[1.0, 3.0], [3.0, 3.0], [3.0, 1.0]]
end

@testset "test geometry" begin
    squares()
    twosquaresT()
    twosquaresC()
    twosquaresC_clock()
    wrapsquares()
end


text:
export E
export b00, b01, b10, b11, plus, minus, iplus, iminus

"""
    E(n::Integer, m::Integer=n, a::Integer, b::Integer)

Computational-basis vectors of `n`x`m` matrices.
"""
function E(n::Integer, m::Integer, a::Integer, b::Integer)
    (1 <= a <= n && 1 <= b <= m) || error("E: Invalid indices, violate 1 <= $a <= $n && 1 <= $b <= $m")
    Em = zeros(Int8,  n, m)
    Em[a, b] = 1
    return Em
end
E(n::Integer, a::Integer, b::Integer) = E(n, n, a, b)

####
#### Bell States
####

"""
    b00, b01, b10, b11

Bell states.
"""
const b00 = (ket(0, 0) + ket(1, 1)) / sqrt(2)
const b10 = (ket(0, 0) - ket(1, 1)) / sqrt(2)
const b01 = (ket(0, 1) + ket(1, 0)) / sqrt(2)
const b11 = (ket(0, 1) - ket(1, 0)) / sqrt(2)
@doc (@doc b00) b01
@doc (@doc b00) b10
@doc (@doc b00) b11

####
#### Plus and minus states
####

"""
    plus, minux, iplus, iminus

The states |+⟩, |-⟩, |+i⟩, |-i⟩. These are eigenstates of
the Pauli X and Y operators.
"""
const plus = (ket(0) + ket(1)) / sqrt(2)
const minus = (ket(0) - ket(1)) / sqrt(2)
const iplus = (ket(0) + im * ket(1)) / sqrt(2)
const iminus = (ket(0) - im * ket(1)) / sqrt(2)
@doc (@doc plus) minus
@doc (@doc plus) iplus
@doc (@doc plus) iminus


text:
mutable struct DropoutDesc
    ptr::Cptr
    states

    function DropoutDesc(droprate::Float64)
        seed = 0
        ref = Ref{Cptr}()
        @cudnn :cudnnCreateDropoutDescriptor (Ptr{Cptr},) ref
        desc = new(ref[])

        h = gethandle()
        ref = Ref{Csize_t}()
        @cudnn :cudnnDropoutGetStatesSize (Cptr,Ptr{Csize_t}) h ref
        states = CuArray{UInt8}(Int(ref[]))
        @cudnn(:cudnnSetDropoutDescriptor,
            (Cptr,Cptr,Cfloat,Cptr,Csize_t,Culonglong),
            desc, h, droprate, states, length(states), seed)
        desc.states = states

        finalizer(desc) do x
            @cudnn :cudnnDestroyDropoutDescriptor (Cptr,) x.ptr
        end
        desc
    end
end

Base.cconvert(::Type{Cptr}, desc::DropoutDesc) = desc.ptr

const DICT_DropoutDesc = Dict()

function dropout(x::CuArray{T,N}, droprate::Float64) where {T,N}
    h = gethandle()
    dropdesc = get!(DICT_DropoutDesc, (h,droprate)) do
        DropoutDesc(droprate)
    end
    xdesc = TensorDesc(x, 4)

    ref = Ref{Csize_t}()
    @cudnn :cudnnDropoutGetReserveSpaceSize (Cptr,Ptr{Csize_t}) xdesc ref
    reserve_space = CuArray{UInt8}(Int(ref[]))

    y = similar(x)
    @cudnn(:cudnnDropoutForward,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Csize_t),
        h, dropdesc, xdesc, x, xdesc, y, reserve_space, length(reserve_space))
    y, (dropdesc,xdesc,reserve_space)
end

function ∇dropout!(dy::CuArray, dx::CuArray, work)
    dropdesc, xdesc, reserve_space = work
    h = gethandle()
    @cudnn(:cudnnDropoutBackward,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Csize_t),
        h, dropdesc, xdesc, dy, xdesc, dx, reserve_space, length(reserve_space))
end


text:
module CUDAapi2
using Compat
include("logging.jl")
include("compatibility.jl")
include("discovery.jl")
export find_library, libnvml
end


text:
naolaplacian00(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = ud2/(2*sqrt(pi)*r)
naolaplacian1m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(3/pi)*sin(θ)*sin(ϕ)*(-2*u + r^2*ud2))/(2*r^3)
naolaplacian10(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(3/pi)*cos(θ)*(-2*u + r^2*ud2))/(2*r^3)
naolaplacian11(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(3/pi)*cos(ϕ)*sin(θ)*(-2*u + r^2*ud2))/ (2*r^3)
naolaplacian2m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((sqrt(15/pi)*sin(θ)^2*sin(2*ϕ)*(-6*u + r^2*ud2))/(4*r^3))
naolaplacian2m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(15/pi)*sin(2*θ)*sin(ϕ)*(-6*u + r^2*ud2))/(4*r^3)
naolaplacian20(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(5/pi)*(1 + 3*cos(2*θ))*(-6*u + r^2*ud2))/(8*r^3)
naolaplacian21(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(15/pi)*cos(ϕ)*sin(2*θ)*(-6*u + r^2*ud2))/(4*r^3)
naolaplacian22(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(15/pi)*cos(2*ϕ)*sin(θ)^2*(-6*u + r^2*ud2))/(4*r^3)
naolaplacian3m3(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(35/(2*pi))*sin(θ)^3*sin(3*ϕ)*(-12*u + r^2*ud2))/(4*r^3)
naolaplacian3m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((sqrt(105/pi)*cos(θ)*sin(θ)^2*sin(2*ϕ)*(-12*u + r^2*ud2))/(4*r^3))
naolaplacian3m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(21/(2*pi))*(3 + 5*cos(2*θ))*sin(θ)*sin(ϕ)*(-12*u + r^2*ud2))/(8*r^3)
naolaplacian30(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(7/pi)*(3*cos(θ) + 5*cos(3*θ))*(-12*u + r^2*ud2))/(16*r^3)
naolaplacian31(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(21/(2*pi))*(3 + 5*cos(2*θ))*cos(ϕ)*sin(θ)*(-12*u + r^2*ud2))/(8*r^3)
naolaplacian32(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(105/pi)*cos(θ)*cos(2*ϕ)*sin(θ)^2*(-12*u + r^2*ud2))/(4*r^3)
naolaplacian33(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(35/(2*pi))*cos(3*ϕ)*sin(θ)^3*(-12*u + r^2*ud2))/(4*r^3)
naolaplacian4m4(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(35/pi)*sin(θ)^4*sin(4*ϕ)*(-20*u + r^2*ud2))/(16*r^3))
naolaplacian4m3(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(35/(2*pi))*cos(θ)*sin(θ)^3*sin(3*ϕ)*(-20*u + r^2*ud2))/(4*r^3)
naolaplacian4m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(5/pi)*(5 + 7*cos(2*θ))*sin(θ)^2*sin(2*ϕ)*(-20*u + r^2*ud2))/(16*r^3))
naolaplacian4m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(5/(2*pi))*(2*sin(2*θ) + 7*sin(4*θ))*sin(ϕ)*(-20*u + r^2*ud2))/(32*r^3)
naolaplacian40(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*(9 + 20*cos(2*θ) + 35*cos(4*θ))*(-20*u + r^2*ud2))/(128*sqrt(pi)*r^3)
naolaplacian41(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(5/(2*pi))*cos(ϕ)*(2*sin(2*θ) + 7*sin(4*θ))*(-20*u + r^2*ud2))/(32*r^3)
naolaplacian42(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(5/pi)*(5 + 7*cos(2*θ))*cos(2*ϕ)*sin(θ)^2*(-20*u + r^2*ud2))/(16*r^3)
naolaplacian43(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(35/(2*pi))*cos(θ)*cos(3*ϕ)*sin(θ)^3*(-20*u + r^2*ud2))/(4*r^3)
naolaplacian44(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(35/pi)*cos(4*ϕ)*sin(θ)^4*(-20*u + r^2*ud2))/(16*r^3)
naolaplacian5m5(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(77/(2*pi))*sin(θ)^5*sin(5*ϕ)*(-30*u + r^2*ud2))/(16*r^3)
naolaplacian5m4(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(385/pi)*cos(θ)*sin(θ)^4*sin(4*ϕ)*(-30*u + r^2*ud2))/(16*r^3))
naolaplacian5m3(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(385/(2*pi))*(7 + 9*cos(2*θ))*sin(θ)^3*sin(3*ϕ)*(-30*u + r^2*ud2))/(32*r^3)
naolaplacian5m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((sqrt(1155/pi)*(5*cos(θ) + 3*cos(3*θ))*sin(θ)^2*sin(2*ϕ)*(-30*u + r^2*ud2))/(32*r^3))
naolaplacian5m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(165/pi)*(15 + 28*cos(2*θ) + 21*cos(4*θ))*sin(θ)*sin(ϕ)*(-30*u + r^2*ud2))/(128*r^3)
naolaplacian50(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(11/pi)*(30*cos(θ) + 35*cos(3*θ) + 63*cos(5*θ))*(-30*u + r^2*ud2))/(256*r^3)
naolaplacian51(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(165/pi)*(15 + 28*cos(2*θ) + 21*cos(4*θ))*cos(ϕ)*sin(θ)*(-30*u + r^2*ud2))/(128*r^3)
naolaplacian52(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(1155/pi)*(5*cos(θ) + 3*cos(3*θ))*cos(2*ϕ)*sin(θ)^2*(-30*u + r^2*ud2))/(32*r^3)
naolaplacian53(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(385/(2*pi))*(7 + 9*cos(2*θ))*cos(3*ϕ)*sin(θ)^3*(-30*u + r^2*ud2))/(32*r^3)
naolaplacian54(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(385/pi)*cos(θ)*cos(4*ϕ)*sin(θ)^4*(-30*u + r^2*ud2))/(16*r^3)
naolaplacian55(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(77/(2*pi))*cos(5*ϕ)*sin(θ)^5*(-30*u + r^2*ud2))/(16*r^3)
naolaplacian6m6(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((sqrt(3003/(2*pi))*sin(θ)^6*sin(6*ϕ)*(-42*u + r^2*ud2))/(32*r^3))
naolaplacian6m5(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(1001/(2*pi))*cos(θ)*sin(θ)^5*sin(5*ϕ)*(-42*u + r^2*ud2))/(16*r^3)
naolaplacian6m4(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(91/pi)*(9 + 11*cos(2*θ))*sin(θ)^4*sin(4*ϕ)*(-42*u + r^2*ud2))/(64*r^3))
naolaplacian6m3(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(1365/(2*pi))*(21*cos(θ) + 11*cos(3*θ))*sin(θ)^3*sin(3*ϕ)*(-42*u + r^2*ud2))/(64*r^3)
naolaplacian6m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((sqrt(1365/(2*pi))*(35 + 60*cos(2*θ) + 33*cos(4*θ))*sin(θ)^2*sin(2*ϕ)*(-42*u + r^2*ud2))/(256*r^3))
naolaplacian6m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(273/pi)*(50*cos(θ) + 45*cos(3*θ) + 33*cos(5*θ))*sin(θ)*sin(ϕ)*(-42*u + r^2*ud2))/(256*r^3)
naolaplacian60(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(13/pi)*(50 + 105*cos(2*θ) + 126*cos(4*θ) + 231*cos(6*θ))*(-42*u + r^2*ud2))/(1024*r^3)
naolaplacian61(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(273/pi)*(50*cos(θ) + 45*cos(3*θ) + 33*cos(5*θ))*cos(ϕ)*sin(θ)*(-42*u + r^2*ud2))/(256*r^3)
naolaplacian62(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(1365/(2*pi))*(35 + 60*cos(2*θ) + 33*cos(4*θ))*cos(2*ϕ)*sin(θ)^2*(-42*u + r^2*ud2))/(256*r^3)
naolaplacian63(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(1365/(2*pi))*(21*cos(θ) + 11*cos(3*θ))*cos(3*ϕ)*sin(θ)^3*(-42*u + r^2*ud2))/(64*r^3)
naolaplacian64(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(91/pi)*(9 + 11*cos(2*θ))*cos(4*ϕ)*sin(θ)^4*(-42*u + r^2*ud2))/(64*r^3)
naolaplacian65(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(1001/(2*pi))*cos(θ)*cos(5*ϕ)*sin(θ)^5*(-42*u + r^2*ud2))/(16*r^3)
naolaplacian66(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(3003/(2*pi))*cos(6*ϕ)*sin(θ)^6*(-42*u + r^2*ud2))/(32*r^3)
naolaplacian7m7(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(715/pi)*sin(θ)^7*sin(7*ϕ)*(-56*u + r^2*ud2))/(64*r^3)
naolaplacian7m6(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(5005/(2*pi))*cos(θ)*sin(θ)^6*sin(6*ϕ)*(-56*u + r^2*ud2))/(32*r^3))
naolaplacian7m5(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(385/pi)*(11 + 13*cos(2*θ))*sin(θ)^5*sin(5*ϕ)*(-56*u + r^2*ud2))/(128*r^3)
naolaplacian7m4(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(385/pi)*(27*cos(θ) + 13*cos(3*θ))*sin(θ)^4*sin(4*ϕ)*(-56*u + r^2*ud2))/(128*r^3))
naolaplacian7m3(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(35/pi)*(189 + 308*cos(2*θ) + 143*cos(4*θ))*sin(θ)^3*sin(3*ϕ)*(-56*u + r^2*ud2))/(512*r^3)
naolaplacian7m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(35/(2*pi))*(350*cos(θ) + 275*cos(3*θ) + 143*cos(5*θ))*sin(θ)^2*sin(2*ϕ)*(-56*u + r^2*ud2))/(512*r^3))
naolaplacian7m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(105/pi)*(350 + 675*cos(2*θ) + 594*cos(4*θ) + 429*cos(6*θ))*sin(θ)*sin(ϕ)*(-56*u + r^2*ud2))/(2048*r^3)
naolaplacian70(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(15/pi)*(175*cos(θ) + 189*cos(3*θ) + 231*cos(5*θ) + 429*cos(7*θ))*(-56*u + r^2*ud2))/(2048*r^3)
naolaplacian71(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(105/pi)*(350 + 675*cos(2*θ) + 594*cos(4*θ) + 429*cos(6*θ))*cos(ϕ)*sin(θ)*(-56*u + r^2*ud2))/(2048*r^3)
naolaplacian72(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(35/(2*pi))*(350*cos(θ) + 275*cos(3*θ) + 143*cos(5*θ))*cos(2*ϕ)*sin(θ)^2*(-56*u + r^2*ud2))/(512*r^3)
naolaplacian73(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(35/pi)*(189 + 308*cos(2*θ) + 143*cos(4*θ))*cos(3*ϕ)*sin(θ)^3*(-56*u + r^2*ud2))/(512*r^3)
naolaplacian74(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(385/pi)*(27*cos(θ) + 13*cos(3*θ))*cos(4*ϕ)*sin(θ)^4*(-56*u + r^2*ud2))/(128*r^3)
naolaplacian75(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(385/pi)*(11 + 13*cos(2*θ))*cos(5*ϕ)*sin(θ)^5*(-56*u + r^2*ud2))/(128*r^3)
naolaplacian76(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(5005/(2*pi))*cos(θ)*cos(6*ϕ)*sin(θ)^6*(-56*u + r^2*ud2))/(32*r^3)
naolaplacian77(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(715/pi)*cos(7*ϕ)*sin(θ)^7*(-56*u + r^2*ud2))/(64*r^3)
naolaplacian8m8(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(12155/pi)*sin(θ)^8*sin(8*ϕ)*(-72*u + r^2*ud2))/(256*r^3))
naolaplacian8m7(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(12155/pi)*cos(θ)*sin(θ)^7*sin(7*ϕ)*(-72*u + r^2*ud2))/(64*r^3)
naolaplacian8m6(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((sqrt(7293/(2*pi))*(13 + 15*cos(2*θ))*sin(θ)^6*sin(6*ϕ)*(-72*u + r^2*ud2))/(128*r^3))
naolaplacian8m5(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(17017/pi)*(11*cos(θ) + 5*cos(3*θ))*sin(θ)^5*sin(5*ϕ)*(-72*u + r^2*ud2))/(256*r^3)
naolaplacian8m4(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(1309/pi)*(99 + 156*cos(2*θ) + 65*cos(4*θ))*sin(θ)^4*sin(4*ϕ)*(-72*u + r^2*ud2))/(1024*r^3))
naolaplacian8m3(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(19635/pi)*(126*cos(θ) + 91*cos(3*θ) + 39*cos(5*θ))*sin(θ)^3*sin(3*ϕ)*(-72*u + r^2*ud2))/(1024*r^3)
naolaplacian8m2(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = -((3*sqrt(595/(2*pi))*(210 + 385*cos(2*θ) + 286*cos(4*θ) + 143*cos(6*θ))*sin(θ)^2*sin(2*ϕ)*(-72*u + r^2*ud2))/(2048*r^3))
naolaplacian8m1(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(17/pi)*(70*sin(2*θ) + 154*sin(4*θ) + 286*sin(6*θ) + 715*sin(8*θ))*sin(ϕ)*(-72*u + r^2*ud2))/(8192*r^3)
naolaplacian80(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(17/pi)*(1225 + 2520*cos(2*θ) + 2772*cos(4*θ) + 3432*cos(6*θ) + 6435*cos(8*θ))*(-72*u + r^2*ud2))/(32768*r^3)
naolaplacian81(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(17/pi)*cos(ϕ)*(70*sin(2*θ) + 154*sin(4*θ) + 286*sin(6*θ) + 715*sin(8*θ))*(-72*u + r^2*ud2))/(8192*r^3)
naolaplacian82(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(595/(2*pi))*(210 + 385*cos(2*θ) + 286*cos(4*θ) + 143*cos(6*θ))*cos(2*ϕ)*sin(θ)^2*(-72*u + r^2*ud2))/(2048*r^3)
naolaplacian83(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(19635/pi)*(126*cos(θ) + 91*cos(3*θ) + 39*cos(5*θ))*cos(3*ϕ)*sin(θ)^3*(-72*u + r^2*ud2))/(1024*r^3)
naolaplacian84(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(1309/pi)*(99 + 156*cos(2*θ) + 65*cos(4*θ))*cos(4*ϕ)*sin(θ)^4*(-72*u + r^2*ud2))/(1024*r^3)
naolaplacian85(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(17017/pi)*(11*cos(θ) + 5*cos(3*θ))*cos(5*ϕ)*sin(θ)^5*(-72*u + r^2*ud2))/(256*r^3)
naolaplacian86(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (sqrt(7293/(2*pi))*(13 + 15*cos(2*θ))*cos(6*ϕ)*sin(θ)^6*(-72*u + r^2*ud2))/(128*r^3)
naolaplacian87(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(12155/pi)*cos(θ)*cos(7*ϕ)*sin(θ)^7*(-72*u + r^2*ud2))/(64*r^3)
naolaplacian88(u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = (3*sqrt(12155/pi)*cos(8*ϕ)*sin(θ)^8*(-72*u + r^2*ud2))/(256*r^3)

naolaplacian(l::Int, m::Int, u::Float64, ud2::Float64, r::Float64, θ::Float64, ϕ::Float64) = naolaplaciantab[l,m](u,ud2,r,θ,ϕ)

naolaplaciantab = Dict(
    (0,0) => naolaplacian00,
    (1,-1) => naolaplacian1m1,
    (1,0) => naolaplacian10,
    (1,1) => naolaplacian11,
    (2,-2) => naolaplacian2m2,
    (2,-1) => naolaplacian2m1,
    (2,0) => naolaplacian20,
    (2,1) => naolaplacian21,
    (2,2) => naolaplacian22,
    (3,-3) => naolaplacian3m3,
    (3,-2) => naolaplacian3m2,
    (3,-1) => naolaplacian3m1,
    (3,0) => naolaplacian30,
    (3,1) => naolaplacian31,
    (3,2) => naolaplacian32,
    (3,3) => naolaplacian33,
    (4,-4) => naolaplacian4m4,
    (4,-3) => naolaplacian4m3,
    (4,-2) => naolaplacian4m2,
    (4,-1) => naolaplacian4m1,
    (4,0) => naolaplacian40,
    (4,1) => naolaplacian41,
    (4,2) => naolaplacian42,
    (4,3) => naolaplacian43,
    (4,4) => naolaplacian44,
    (5,-5) => naolaplacian5m5,
    (5,-4) => naolaplacian5m4,
    (5,-3) => naolaplacian5m3,
    (5,-2) => naolaplacian5m2,
    (5,-1) => naolaplacian5m1,
    (5,0) => naolaplacian50,
    (5,1) => naolaplacian51,
    (5,2) => naolaplacian52,
    (5,3) => naolaplacian53,
    (5,4) => naolaplacian54,
    (5,5) => naolaplacian55,
    (6,-6) => naolaplacian6m6,
    (6,-5) => naolaplacian6m5,
    (6,-4) => naolaplacian6m4,
    (6,-3) => naolaplacian6m3,
    (6,-2) => naolaplacian6m2,
    (6,-1) => naolaplacian6m1,
    (6,0) => naolaplacian60,
    (6,1) => naolaplacian61,
    (6,2) => naolaplacian62,
    (6,3) => naolaplacian63,
    (6,4) => naolaplacian64,
    (6,5) => naolaplacian65,
    (6,6) => naolaplacian66,
    (7,-7) => naolaplacian7m7,
    (7,-6) => naolaplacian7m6,
    (7,-5) => naolaplacian7m5,
    (7,-4) => naolaplacian7m4,
    (7,-3) => naolaplacian7m3,
    (7,-2) => naolaplacian7m2,
    (7,-1) => naolaplacian7m1,
    (7,0) => naolaplacian70,
    (7,1) => naolaplacian71,
    (7,2) => naolaplacian72,
    (7,3) => naolaplacian73,
    (7,4) => naolaplacian74,
    (7,5) => naolaplacian75,
    (7,6) => naolaplacian76,
    (7,7) => naolaplacian77,
    (8,-8) => naolaplacian8m8,
    (8,-7) => naolaplacian8m7,
    (8,-6) => naolaplacian8m6,
    (8,-5) => naolaplacian8m5,
    (8,-4) => naolaplacian8m4,
    (8,-3) => naolaplacian8m3,
    (8,-2) => naolaplacian8m2,
    (8,-1) => naolaplacian8m1,
    (8,0) => naolaplacian80,
    (8,1) => naolaplacian81,
    (8,2) => naolaplacian82,
    (8,3) => naolaplacian83,
    (8,4) => naolaplacian84,
    (8,5) => naolaplacian85,
    (8,6) => naolaplacian86,
    (8,7) => naolaplacian87,
    (8,8) => naolaplacian88
)

max_stars_repo_name: johncsnyder/Element


text:
include("Include.jl")

# load the data -
path_data_file_TF = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF.csv")
full_data_table_TF = load(path_data_file_TF)
path_data_file_TF_TM = joinpath(_PATH_TO_DATA, "Training-Thrombin-TF-TM.csv")
full_data_table_TF_TM = load(path_data_file_TF_TM)
full_training_data_frame = vcat(full_data_table_TF, full_data_table_TF_TM)

# filter -
experimental_data_table = filter(:visitid => x -> x == 1, full_training_data_frame)

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
module TestShow

include("preamble.jl")

function nullshows(x)
    io = devnull
    # io = stdout
    show(io, x)
    show(IOContext(io, :limit => true), x)
    show(io, "text/plain", x)
    show(IOContext(io, :limit => true), "text/plain", x)
    return true
end

@testset begin
    n = 3
    C = zeros(n)
    A = rand(-10:10, n, n)
    B = rand(-10:10, n)
    CAB = ndreducible(
        instantiate(broadcasted(Ref, Ref(C), eachindex(C))) => (:i,),
        A => (:i, :j),
        B => (:j,),
    )
    @test nullshows(CAB.plan.indices[1])
    @test nullshows(CAB.plan)
    @test nullshows(CAB)
    @test nullshows(CAB.patterns)
    @test nullshows(CAB.patterns[1])
end

end  # module


text:
module TestFactorizationMachines

using FactorizationMachines
using Base.Test

# Matrix constructed as the outer product of two vectors: [.1; .2] and [.1 .2]
num_factors = 4
X = [
1. 1. 0. 0.;
0. 0. 1. 1.;
1. 0. 1. 0.;
0. 1. 0. 1.
]
y = [1., 2., 2., 4.]

info("Decomposing Rank 1 Matrix")
fm = train(sparse(X), y;
    method = Methods.sgd(alpha = 0.15),
    model_params = Models.gauss(k0 = false, k1 = false))
@test fm.model.w0 == 0
@test fm.model.w == zeros(num_factors)
@test_approx_eq_eps 0.0 evaluate(Evaluators.rmse(), fm, sparse(X), y) 1e-7

# Matrix constructed as the outer product of two vectors: [-1.; 1.] and [1. 1.]
num_factors = 4
X = [
1. 1. 0. 0.;
0. 0. 1. 1.;
1. 0. 1. 0.;
0. 1. 0. 1.
]
y = [0., 1., 0., 1.]

info("Decomposing Rank 1 Matrix")
fm = train(sparse(X), y;
    method = Methods.sgd(alpha = 0.15),
    model_params = Models.gauss(k0 = false, k1 = false),
    task_params = Tasks.classification())
@test fm.model.w0 == 0
@test fm.model.w == zeros(num_factors)
@test evaluate(Evaluators.heaviside(), fm, sparse(X), y) < 1.0

end


text:
function Progress(cur_run_count)
  cur_progress_meter = PmapProgressMeter.Progress(cur_run_count)
  PmapProgressMeter.update!(cur_progress_meter, 0)
  sleep(0.1)

  cur_progress_meter
end


text:
"""
    mblocks(X, listbl)
Make blocks from a matrix.
* `X` : X-data.
* `listbl` : A vector whose each component defines the colum numbers
    defining a block in `X` . The length of `listbl` is the number
    of blocks.

The function returns a list of blocks.

## Examples
```julia
n = 5 ; p = 10
X = rand(n, p)
listbl = [3:4, 1, [6; 8:10]]

X_bl = mblocks(X, listbl)
```
"""
function mblocks(X, listbl)
    nbl = length(listbl)
    zX = list(nbl, Matrix{Float64})
    @inbounds for i = 1:nbl
        zX[i] = ensure_mat(X[:, listbl[i]])
    end
    zX
end

"""
    blockscal(X, weights = ones(size(X, 1)); listbl, scal = nothing)
Autoscale blocks of a matrix.
* `X` : List (vector) of blocks (matrices) of X-data.
    Each component of the list is a block.
* `weights` : Weights of the observations (rows). Only used if `scal = nothing`.
* `scal` : If `nothing`, each block is autoscaled (i.e. divided) by the
    the square root of the sum of the variances of each column of the block.
    Else, `scal` must be a vector defining the scaling values dividing the blocks.

The function returns the scaled blocks, and the scaling values.

Vector `weights` (row-weighting) is internally normalized to sum to 1.

## Examples
```julia
n = 5 ; p = 10
X = rand(n, p)
listbl = [3:4, 1, [6; 8:10]]
X_bl = mblocks(X, listbl)

scal = nothing
#scal = [3.1 2 .7]
res = blockscal(X_bl; scal = scal) ;
res.X
res.scal
sum(colvars(res.X[1]))

X_concat = reduce(hcat, res.X)
```
"""
function blockscal(X, weights = ones(size(X[1], 1)); scal = nothing)
    weights = mweights(weights)
    nbl = length(X)
    if isnothing(scal)
        scal = list(nbl, Float64)
        @inbounds for i = 1:nbl
            scal[i] = sqrt(sum(colvars(X[i], weights)))
        end
    end
    Xs = list(nbl, Matrix{Float64})
    @inbounds for i = 1:nbl
        Xs[i] = X[i] / scal[i]
    end
    (X = Xs, scal = scal)
end


text:
function shuffle_and_insert(p::AbstractParticles, ind, val)
    p = deepcopy(p)
    p.particles[ind] = p.particles[1]
    p.particles[1] = val
    p
end

function shuffle_and_insert(p::StaticParticles, ind, val)
    part = p.particles
    part = setindex(part, part[1], ind)
    part = setindex(part, val, 1)
    StaticParticles(part)
end

"""
    pn = with_nominal(p, val)

Endow particles `p` with a nominal value `val`. The particle closest to `val` will be replaced with val, and moved to index 1. This operation introduces a slight bias in the statistics of `pn`, but the operation is asymptotically unbiased for large sample sizes. To obtain the nominal value of `pn`, call `nominal(pn)`.
"""
function with_nominal(p::AbstractParticles, val)
    minind = argmin(abs.(p.particles .- val))
    shuffle_and_insert(p, minind, val)
end

function with_nominal(p::MvParticles, val::AbstractVector)
    M = Matrix(p)
    minind = argmin(vec(sum(abs2, M .- val', dims=1)))
    shuffle_and_insert.(p, minind, val)
end

"""
    nominal(p)

Return the nominal value of `p` (assumes that `p` has been endowed with a nominal value using `with_nominal`).
"""
nominal(p::AbstractParticles) = p.particles[1]
nominal(p::MvParticles) = nominal.(p)
nominal(P) = replace_particles(P, replacer=nominal)


text:
module TestDataFrames
	using Base.Test
	using DataArrays
	using DataFrames

	##########
	## paste
	##########

	@assert paste(["a", "b"], "X", [1:2]) == ["aX1", "bX2"]
	@assert paste(["a", "b"], "X", [1:4]) == ["aX1", "bX2", "aX3", "bX4"]

	##########
	## cut
	##########

	@assert isequal(cut([2, 3, 5], [1, 3, 6]), PooledDataArray(["(1,3]", "(1,3]", "(3,6]"]))
	@assert isequal(cut([2, 3, 5], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]"]))
	@assert isequal(cut([2, 3, 5, 6], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]", "(3,6]"]))
	@assert isequal(cut([1, 2, 4], [1, 3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
	@assert isequal(cut([1, 2, 4], [3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
	@assert isequal(cut([1, 2, 4], [3]), PooledDataArray(["[1,3]", "[1,3]", "(3,4]"]))
	@assert isequal(cut([1, 5, 7], [3, 6]), PooledDataArray(["[1,3]", "(3,6]", "(6,7]"]))

	ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
	bins = [18, 25, 35, 60, 100]
	cats = cut(ages, bins)
	pdv = PooledDataArray(["(18,25]", "(18,25]", "(18,25]",
	                       "(25,35]", "(18,25]", "(18,25]",
	                       "(35,60]", "(25,35]", "(60,100]",
	                       "(35,60]", "(35,60]", "(25,35]"])
	@assert isequal(cats, pdv)

	##########
	## rep
	##########

	@assert rep(3, 2) == [3,3]
	@assert rep([3,4], 2) == [3,4,3,4]
	@assert rep([3,4], 1, 2) == [3,3,4,4]
	@assert rep([3,4], each = 2) == [3,3,4,4]
	@assert rep([3,4], times = 2) == [3,4,3,4]
	@assert rep([3,4], times = [2,3]) == [3,3,4,4,4]
	@assert rep([3,4], [2,3]) == [3,3,4,4,4]
	@assert isequal(rep(@data([NA, 3, 4]), 2), @data([NA, 3, 4, NA, 3, 4]))
	@assert isequal(rep(@data([NA, 3, 4]), [2, 1, 2]), @data([NA, NA, 3, 4, 4]))
	@assert isequal(rep(@data([NA, 3, 4]), [2, 1, 0]), @data([NA, NA, 3]))
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: accessanalyzer

using AWS.Compat
using AWS.UUIDs
"""
    CreateAnalyzer()

Creates an analyzer for your account.

# Required Parameters
- `analyzerName`: The name of the analyzer to create.
- `type`: The type of analyzer to create. Only ACCOUNT analyzers are supported. You can create only one analyzer per account per Region.

# Optional Parameters
- `archiveRules`: Specifies the archive rules to add for the analyzer. Archive rules automatically archive findings that meet the criteria you define for the rule.
- `clientToken`: A client token.
- `tags`: The tags to apply to the analyzer.
"""

create_analyzer(analyzerName, type; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer", Dict{String, Any}("analyzerName"=>analyzerName, "type"=>type, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_analyzer(analyzerName, type, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerName"=>analyzerName, "type"=>type, "clientToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    CreateArchiveRule()

Creates an archive rule for the specified analyzer. Archive rules automatically archive findings that meet the criteria you define when you create the rule.

# Required Parameters
- `analyzerName`: The name of the created analyzer.
- `filter`: The criteria for the rule.
- `ruleName`: The name of the rule to create.

# Optional Parameters
- `clientToken`: A client token.
"""

create_archive_rule(analyzerName, filter, ruleName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule", Dict{String, Any}("filter"=>filter, "ruleName"=>ruleName, "clientToken"=>string(uuid4())); aws_config=aws_config)
create_archive_rule(analyzerName, filter, ruleName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filter"=>filter, "ruleName"=>ruleName, "clientToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    DeleteAnalyzer()

Deletes the specified analyzer. When you delete an analyzer, Access Analyzer is disabled for the account in the current or specific Region. All findings that were generated by the analyzer are deleted. You cannot undo this action.

# Required Parameters
- `analyzerName`: The name of the analyzer to delete.

# Optional Parameters
- `clientToken`: A client token.
"""

delete_analyzer(analyzerName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
delete_analyzer(analyzerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    DeleteArchiveRule()

Deletes the specified archive rule.

# Required Parameters
- `analyzerName`: The name of the analyzer that associated with the archive rule to delete.
- `ruleName`: The name of the rule to delete.

# Optional Parameters
- `clientToken`: A client token.
"""

delete_archive_rule(analyzerName, ruleName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
delete_archive_rule(analyzerName, ruleName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    GetAnalyzedResource()

Retrieves information about a resource that was analyzed.

# Required Parameters
- `analyzerArn`: The ARN of the analyzer to retrieve information from.
- `resourceArn`: The ARN of the resource to retrieve information about.

"""

get_analyzed_resource(analyzerArn, resourceArn; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzed-resource", Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn); aws_config=aws_config)
get_analyzed_resource(analyzerArn, resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzed-resource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn), args)); aws_config=aws_config)

"""
    GetAnalyzer()

Retrieves information about the specified analyzer.

# Required Parameters
- `analyzerName`: The name of the analyzer retrieved.

"""

get_analyzer(analyzerName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)"; aws_config=aws_config)
get_analyzer(analyzerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)", args; aws_config=aws_config)

"""
    GetArchiveRule()

Retrieves information about an archive rule.

# Required Parameters
- `analyzerName`: The name of the analyzer to retrieve rules from.
- `ruleName`: The name of the rule to retrieve.

"""

get_archive_rule(analyzerName, ruleName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)"; aws_config=aws_config)
get_archive_rule(analyzerName, ruleName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", args; aws_config=aws_config)

"""
    GetFinding()

Retrieves information about the specified finding.

# Required Parameters
- `analyzerArn`: The ARN of the analyzer that generated the finding.
- `id`: The ID of the finding to retrieve.

"""

get_finding(analyzerArn, id; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/finding/$(id)", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
get_finding(analyzerArn, id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/finding/$(id)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), args)); aws_config=aws_config)

"""
    ListAnalyzedResources()

Retrieves a list of resources of the specified type that have been analyzed by the specified analyzer..

# Required Parameters
- `analyzerArn`: The ARN of the analyzer to retrieve a list of analyzed resources from.

# Optional Parameters
- `maxResults`: The maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned.
- `resourceType`: The type of resource.
"""

list_analyzed_resources(analyzerArn; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/analyzed-resource", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
list_analyzed_resources(analyzerArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/analyzed-resource", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), args)); aws_config=aws_config)

"""
    ListAnalyzers()

Retrieves a list of analyzers.

# Optional Parameters
- `maxResults`: The maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned.
- `type`: The type of analyzer.
"""

list_analyzers(; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer"; aws_config=aws_config)
list_analyzers(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer", args; aws_config=aws_config)

"""
    ListArchiveRules()

Retrieves a list of archive rules created for the specified analyzer.

# Required Parameters
- `analyzerName`: The name of the analyzer to retrieve rules from.

# Optional Parameters
- `maxResults`: The maximum number of results to return in the request.
- `nextToken`: A token used for pagination of results returned.
"""

list_archive_rules(analyzerName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule"; aws_config=aws_config)
list_archive_rules(analyzerName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/analyzer/$(analyzerName)/archive-rule", args; aws_config=aws_config)

"""
    ListFindings()

Retrieves a list of findings generated by the specified analyzer.

# Required Parameters
- `analyzerArn`: The ARN of the analyzer to retrieve findings from.

# Optional Parameters
- `filter`: A filter to match for the findings to return.
- `maxResults`: The maximum number of results to return in the response.
- `nextToken`: A token used for pagination of results returned.
- `sort`: The sort order for the findings returned.
"""

list_findings(analyzerArn; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/finding", Dict{String, Any}("analyzerArn"=>analyzerArn); aws_config=aws_config)
list_findings(analyzerArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/finding", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn), args)); aws_config=aws_config)

"""
    ListTagsForResource()

Retrieves a list of tags applied to the specified resource.

# Required Parameters
- `resourceArn`: The ARN of the resource to retrieve tags from.

"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    StartResourceScan()

Immediately starts a scan of the policies applied to the specified resource.

# Required Parameters
- `analyzerArn`: The ARN of the analyzer to use to scan the policies applied to the specified resource.
- `resourceArn`: The ARN of the resource to scan.

"""

start_resource_scan(analyzerArn, resourceArn; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/resource/scan", Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn); aws_config=aws_config)
start_resource_scan(analyzerArn, resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/resource/scan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "resourceArn"=>resourceArn), args)); aws_config=aws_config)

"""
    TagResource()

Adds a tag to the specified resource.

# Required Parameters
- `resourceArn`: The ARN of the resource to add the tag to.
- `tags`: The tags to add to the resource.

"""

tag_resource(resourceArn, tags; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes a tag from the specified resource.

# Required Parameters
- `resourceArn`: The ARN of the resource to remove the tag from.
- `tagKeys`: The key for the tag to add.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateArchiveRule()

Updates the criteria and values for the specified archive rule.

# Required Parameters
- `analyzerName`: The name of the analyzer to update the archive rules for.
- `filter`: A filter to match for the rules to update. Only rules that match the filter are updated.
- `ruleName`: The name of the rule to update.

# Optional Parameters
- `clientToken`: A client token.
"""

update_archive_rule(analyzerName, filter, ruleName; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}("filter"=>filter, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_archive_rule(analyzerName, filter, ruleName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/analyzer/$(analyzerName)/archive-rule/$(ruleName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("filter"=>filter, "clientToken"=>string(uuid4())), args)); aws_config=aws_config)

"""
    UpdateFindings()

Updates the status for the specified findings.

# Required Parameters
- `analyzerArn`: The ARN of the analyzer that generated the findings to update.
- `status`: The state represents the action to take to update the finding Status. Use ARCHIVE to change an Active finding to an Archived finding. Use ACTIVE to change an Archived finding to an Active finding.

# Optional Parameters
- `clientToken`: A client token.
- `ids`: The IDs of the findings to update.
- `resourceArn`: The ARN of the resource identified in the finding.
"""

update_findings(analyzerArn, status; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/finding", Dict{String, Any}("analyzerArn"=>analyzerArn, "status"=>status, "clientToken"=>string(uuid4())); aws_config=aws_config)
update_findings(analyzerArn, status, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = accessanalyzer("PUT", "/finding", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("analyzerArn"=>analyzerArn, "status"=>status, "clientToken"=>string(uuid4())), args)); aws_config=aws_config)


text:
"""
    volt_to_pressure(15e6; hydrophone_id = :Onda_HGL0200_2322, preamp_id = :Onda_AH2020_1238_20dB)
    volt_to_pressure(15e6; hydrophone_id = :PA_3197)

# Arguments
- `f0`: Frequency in Hz.
- `hydrophone_id`: Currently accepts two: :Onda_HGL0200_2322, :PA_3197.
- `preamp_id`: Only for :Onda_HGL0200_2322. Choices are: :Onda_AH2020_1238_0dB, :Onda_AH2020_1238_20dB.
# Returns
- `factor`: conversion factor in form of pascal/volt.
"""
function volt_to_pressure(f0, hydrophone_id::Symbol)
    f0 = enforce_hertz_unit(f0)
    calibration_data, _ = get_calibration(hydrophone_id)
    frequencies = calibration_data[:, 1] * u"MHz"

    if occursin("PA_", String(hydrophone_id))
        volt_per_pascal = calibration_data[:, 2] * u"mV/MPa"
    else
        volt_per_pascal = calibration_data[:, 2] * u"dBV/μPa"
    end

    volt_per_pascal = interp(frequencies, volt_per_pascal, f0)
    return calculate_pascal_per_volt(volt_per_pascal)
end

function calculate_pascal_per_volt(volt_per_pascal)
    volt_per_pascal = uconvertrp(u"V/Pa", volt_per_pascal)
    pascal_per_volt = 1/volt_per_pascal
    return pascal_per_volt
end

function volt_to_pressure(f0, hydrophone_id::Symbol, preamp_id::Symbol)
    factor, _ = volt_to_pressure_and_phase(f0, hydrophone_id, preamp_id)
    return factor
end

function volt_to_pressure_and_phase(f0, hydrophone_id::Symbol, preamp_id::Symbol)
    f0 = enforce_hertz_unit(f0)
    calibration_data, preamp_data = get_calibration(hydrophone_id, preamp_id)

    hydrophone_freqs = calibration_data[:, 1] * u"MHz"
    hydrophone_volt_per_pascal = calibration_data[:, 2] * u"dBV/μPa"
    hydrophone_capacitances = calibration_data[:, 5] * u"pF"
    hydrophone_volt_per_pascal_at_f0 = interp(hydrophone_freqs, hydrophone_volt_per_pascal, f0)
    hydrophone_capacitance_at_f0 = interp(hydrophone_freqs, hydrophone_capacitances, f0)

    preamp_freqs = preamp_data[:,1] * u"MHz"
    preamp_gains = uconvertrp.(NoUnits, preamp_data[:, 2] * u"dB")
    preamp_phases = preamp_data[:, 3]
    preamp_capacitances = preamp_data[:, 4] * u"pF"
    preamp_gain_at_f0 = interp(preamp_freqs, preamp_gains, f0)
    preamp_phase_at_f0 = interp(preamp_freqs, preamp_phases, f0)
    preamp_capacitance_at_f0 = interp(preamp_freqs, preamp_capacitances, f0)

    voltage_factor = hydrophone_capacitance_at_f0 ./ (hydrophone_capacitance_at_f0 + preamp_capacitance_at_f0)

    # Combine result
    combined_volt_per_pascal = preamp_gain_at_f0 * hydrophone_volt_per_pascal_at_f0 * voltage_factor

    return calculate_pascal_per_volt(combined_volt_per_pascal), preamp_phase_at_f0
end

function interp(frequency, value, f0)
    f_lower, f_upper = frequency[1], frequency[end]
    if f0 < f_lower || f0 > f_upper
        throw("f0 must be between $f_lower-$f_upper")
    end

    interpolant = LinearInterpolation(frequency, value)
    interpolated_value = interpolant(f0)
    return interpolated_value
end

function enforce_hertz_unit(center_frequency)
    if isa(center_frequency, Unitful.Frequency)
        center_frequency = uconvert(u"Hz", center_frequency)
    elseif isa(center_frequency, Real)
        center_frequency = center_frequency*u"Hz"
    else
        error("Center frequency must have a frequency unit or be a Real number.")
    end

    return center_frequency
end


text:
using HTTP
using JSON
## !!! needs Julia 1.3:w
import Base.Threads.@spawn

const TEST_URL = "http://localhost:8000/stb"
const LIVE_URL = "https://oustb.virtual-worlds.scot/oustb/stb/"

function getdata( rate1:: Real, rate2 :: Real ) :: Dict
    url = "$(LIVE_URL)/?it_allow=12500&it_rate_1=$rate1&it_rate_$rate2&it_band=50000"
    println( "fetching from URL $url" )
    resp = HTTP.request( "GET", url )
    json = JSON.parse(join((map(Char,resp.body))))
    return json;
end

function doRunBatch( max :: Integer ) :: String
    n = 0
    threadid = Threads.threadid()
    println("doRunBatch; running on thread $threadid")
    rc = @timed for r1 in 2968:30
        for r2 in 40:50
            println( "getting data r1=$r1 r2=$r2 " )
            json = getdata( r1, r2 )
            if n % 1 == 0
                print( json )
            end
            n += 1
            if n > max
                break
            end
        end # rand2
    end # rand1
    secs = rc[2]
    "total runs on thread $threadid = $n in $secs secs"
end # func doRunBatch

# FIXME something unsafe about connections
function doall_threads( threadno :: Integer )
    println(" at start" )
    out = []
    for i in 1:threadno
        @time response = @spawn doRunBatch(1_000)
        s = fetch( response )
        push!( out, s )
    end
    out
end

function doall( threadno :: Integer )
    println(" at start" )
    doRunBatch(1_000)
end

out = doall(1)
println( out )


text:
include("task_func.jl")

function draw(r)
    goToSide!(r, West)
    goToSide!(r, Nord)
    count = 0
    flag = 0
    while(!isborder(r, Sud))
        while(!isborder(r, Ost))
            move!(r, Ost)
            if(flag == 0 && isborder(r, Sud))
                flag = 1
                count += 1
            elseif(flag == 1 && !isborder(r, Sud))
                flag = 0
            end
        end
        move!(r, Sud)
        goToSide!(r, West)
    end
    print(count)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
@testset "getters" begin
b = BGEN.Bgen(example_8bits)
@test fsize(b) == 128746
@test all(samples(b) .== b.samples)
@test n_samples(b) == 500
@test n_variants(b) == 199
@test compression(b) == "Zlib"

v = first(iterator(b; from_bgen_start=true))
@test n_samples(v::Variant) == 500
@test varid(v::Variant) == "SNPID_2"
@test rsid(v::Variant) == "RSID_2"
@test chrom(v::Variant) == "01"
@test pos(v::Variant) == 2000
@test n_alleles(v::Variant) == 2
@test length(alleles(v::Variant)) == 2
@test all(alleles(v) .== ["A", "G"])

minor_allele_dosage!(b, v)
@test phased(v) == 0
@test min_ploidy(v) == 2
@test max_ploidy(v) == 2
@test all(ploidy(v) .== 2)
@test bit_depth(v) == 8
@test length(missings(v)) == 1
@test missings(v)[1] == 1
@test minor_allele(v) == "A"
@test major_allele(v) == "G"
end


text:
const updatevalue = js"function () {_webIOScope.setObservableValue('value', this.value);}"
const countchanges = js"""
function () {
    var val = _webIOScope.getObservableValue('changes') + 1;
    _webIOScope.setObservableValue('changes', val);
}
"""

function input(;
               imports = String[],
               addevents = Dict(:input => updatevalue, :change => countchanges),
               value = Observable(""),
               kwargs...)

    scp = Scope(imports = imports)
    setobservable!(scp, "changes", Observable(0))
    addreactivity!(scp, :input; addevents = addevents, value = value, kwargs...)

    Widget{:input}(
        scope = scp,
        output = value,
        layout = x -> node(:div, scope(x), className = "interact-widget")
    )
end

struct TextBox{T} <: AbstractVanillaWidget{T}
    value::AbstractObservable{T}
    props::Dict{Symbol, Any}
    function TextBox(value::AbstractObservable{T}; className = "input", kwargs...) where {T}
        props = Dict{Symbol, Any}(observify(pair) for pair in kwargs)
        new{T}(value, props)
    end
end
TextBox(value; kwargs...) = TextBox(to_observable(value); kwargs...)

value(a::TextBox) = getfield(a, 1)
WebIO.props(a::TextBox) = getfield(a, 2)

function WebIO.render(a::TextBox)
    input()