text:
#
# utils.jl --
#
# General purpose methods.
#

"""
    to_type(T, x)

yields `x` converted to type `T`, the result is asserted to be of type `T`.

"""
to_type(::Type{T}, x::T) where {T} = x
to_type(::Type{T}, x::Any) where {T} = convert(T, x)::T

"""
    to_int(x)

converts `x` to a similar object whose values are all of type `Int`.  Argument
`x` can be a scalar, an array, a tuple or a range.  Argument `x` is returned if
already of the correct type.

"""
to_int(x::Int) = x
to_int(x::Integer) = Int(x)
to_int(I::AbstractUnitRange{Int}) = I
to_int(I::AbstractUnitRange{<:Integer}) = to_int(first(I)):to_int(last(I))
to_int(I::StepRange{Int,Int}) = I
to_int(I::StepRange{<:Integer,<:Integer}) =
    to_int(first(I)):to_int(step(I)):to_int(last(I))
to_int(A::AbstractArray{Int,N}) where {N} = A
to_int(A::AbstractArray{<:Integer,N}) where {N} = map(to_int, A)
to_int(T::NTuple{N,Int}) where {N} = T
to_int(T::NTuple{N,Integer}) where {N} = map(to_int, T)

"""
    promote_eltype(A, B, ...) -> T

yields an element type `T` resulting from calling `promote_type` onto the
element types of all the arguments which must be arrays or array types.

"""
promote_eltype(args::Union{AbstractArray,Type{<:AbstractArray}}...) =
    promote_type(map(eltype, args)...)

"""
    all_match(val, f, args...) -> bool

yields as soon as possible (short-circuit) whether `f(arg) == val` for each
argument `arg` in `args...`.  The returned value is `true` if there are no
arguments after `f`.

"""
all_match(val, f::Function) = true
all_match(val, f::Function, A) = f(A) == val
@inline all_match(val, f::Function, A, B...) =
    all_match(val, f, A) && all_match(val, f::Function, B...)

"""
    allof(f, args...) -> Bool

checks whether predicate function `f` returns `true` for all arguments in
`args...`, returning `false` as soon as possible (short-circuiting).

    allof(args...) -> Bool

checks whether all arguments `args...` are `true`, returning `false` as soon as
possible (short-circuiting).  Arguments can be booleans or arrays of booleans.
The latter are considered as `true` if all their elements are `true` and are
considered as `false` otherwise (if any of their elements are `false`).
Arguments can also be iterables to check whether all their values are `true`.
An empty iterable is considered as `true`.

This method can be much faster than `all(f, args)` or `all(args)` because its
result may be determined at compile time.  However, `missing` values are not
considered as special.

See also `all`, [`anyof`](@ref), [`noneof`](@ref).

"""
allof(f::Function, a) = f(a)::Bool
@inline allof(f::Function, a, b...) = allof(f, a) && allof(f, b...)
@inline allof(a, b...) = allof(a) && allof(b...)
allof(a::Bool) = a
function allof(a::AbstractArray{Bool})
    @inbounds for i in eachindex(a)
        a[i] || return false
    end
    return true
end
function allof(itr)
    for val in itr
        allof(val) || return false
    end
    return true
end

"""
    anyof(f, args...) -> bool

checks whether predicate function `f` returns `true` for any argument
`args...`, returning `true` as soon as possible (short-circuiting).

    anyof(args...) -> bool

checks whether all arguments `args...` are `true`, returning `false` as soon as
possible (short-circuiting).  Arguments can be booleans or arrays of booleans.
The latter are considered as `true` if any of their elements are `true` and are
considered as `false` otherwise (if all their elements are `false`).  Arguments
can also be iterables to check whether any of their values are `true`.  An
empty iterable is considered as `false`.

This method can be much faster than `any(f, args)` or `any(args)` because its
result may be determined at compile time.  However, `missing` values are not
considered as special.

See also `any`, [`allof`](@ref), [`noneof`](@ref).

"""
anyof(f::Function, a) = f(a)::Bool
@inline anyof(f::Function, a, b...) = anyof(f, a) || anyof(f, b...)
@inline anyof(a, b...) = anyof(a) || anyof(b...)
anyof(a::Bool) = a
function anyof(a::AbstractArray{Bool})
    @inbounds for i in eachindex(a)
        a[i] && return true
    end
    return false
end
function anyof(itr)
    for val in itr
        anyof(val) && return true
    end
    return false
end

"""
    noneof(f, args...) -> bool

checks whether predicate `f` returns `false` for all argument `args...`, while

    noneof(args...) -> bool

checks whether all argument `args...` are false.

See also `any`, [`allof`](@ref), [`noneof`](@ref).

"""
@inline noneof(args...) = ! anyof(args...)

"""
    reversemap(f, args)

applies the function `f` to arguments `args` in reverse order and return the
result.  For now, the arguments `args` must be in the form of a simple tuple
and the result is the tuple: `(f(args[end]),f(args[end-1]),...,f(args[1])`.

Also see: `map`, `ntuple`.

"""
reversemap(f::Function, args::NTuple{N,Any}) where {N} =
    ntuple(i -> f(args[(N + 1) - i]), Val{N}())

"""
    strictmap!(dst, f, src) -> dst

does `dst[i] = f(src[i])` for all indices `i` and returns `dst`.  Arguments
`dst` and `src` must have the same axes.

Except for the strict condition on the axes, this method is similar to
`map!(f,dst,src)`.

"""
function strictmap!(dst::AbstractArray{<:Any,N}, f,
                    src::AbstractArray{<:Any,N}) where {N}
    @inbounds @simd for i in all_indices(dst, src)
        dst[i] = f(src[i])
    end
    return dst
end


text:
# This file was generated, do not modify it. # hide
X3 = hcat(X.LStat, X.LStat.^2)
machine(mdl, X3, y)
fit!(mach)
ŷ = predict(mach, X3)
round(rms(ŷ, y), sigdigits=4)

ext: jl

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
module BTGModel

using Reexport
using SpecialFunctions
using FastGaussQuadrature
using Distributions
using LinearAlgebra

include("./BTGFuncs.jl")

@reexport using .BTGFuncs

export
    Model,
    btgdensity,
    btgdistribution,
    computeweights

"""
    BTGModel

Settings for the Bayesian Transformed Gaussian model.
"""
struct Model{K<:Correlation,G<:Transform,V<:Covariate}
    k::K
    g::G
    f::V
    X::Array{Float64, 2}
    Y::Array{Float64, 1}
    knumgauss::Int
    gnumgauss::Int
end

function gaussweights(d::Uniform, n)
    x, w = gausslegendre(n)
    return ((d.b - d.a) .* (x .+ 1) ./ 2 .+ d.a), w ./ 2
end

function gaussweights(d::Normal, n)
    x, w = gausshermite(n)
    return (sqrt(2) * d.σ) .* x .+ d.μ, w ./ sqrt(π)
end

function gaussweights(d::Gamma, n)
    x, w = gausslaguerre(n, d.α - 1)
    return d.θ .* x, w ./ gamma(d.α)
end

function gaussweights(d::Exponential, n)
    x, w = gausslaguerre(n)
    return d.θ .* x, w
end

function weighttensor(f::Parametric, n)
    tups = (gaussweights(getfield(f, s), n) for s in fieldnames(typeof(f)))
    x = collect(Iterators.product((t[1] for t in tups)...))
    w = map(x -> *(x...), Iterators.product((t[2] for t in tups)...))
    return reshape(x, length(x)), reshape(w, length(w))
end

function computeweights(mdl::Model)
    kgauss = weighttensor(mdl.k, mdl.knumgauss)
    ggauss = weighttensor(mdl.g, mdl.gnumgauss)
    return (kgauss, ggauss)
end

struct HSSSolver end # TODO Yuanxi Solver

# function Base.size(s::HSSSolver) end # TODO

# function Base.:\(s::HSSSolver, y) end # TODO

# function LinearAlgebra.det(s::HSSSolver) end # TODO

function makesolver!(A)
    return cholesky!(A) # TODO Placeholder
end

# TODO Clean this function up, make it faster
function computedist(mdl::Model, x::Array{Float64, 1}, y::Float64, ((kx, kw), (gx, gw)), f)
    n, p = size(mdl.X)

    c = 0
    acc = 0

    for i = 1:length(kx)
        θ = kx[i]
        α = kw[i]

        cholΣ = makesolver!(Symmetric(kernelmatrix(mdl.k, θ, mdl.X, mdl.X)))
        cholXΣX = makesolver!(Symmetric(mdl.X' * (cholΣ \ mdl.X)))
        B = kernelmatrix(mdl.k, θ, reshape(x, 1, length(x)), mdl.X)
        D = 1.0 - (B * (cholΣ \ B'))[1]
        H = x .- B * (cholΣ \ mdl.X)
        C = D + (H * (cholXΣX \ H'))[1]

        for j = 1:length(gx)
            λ = gx[j]
            α′ = gw[j]

            Z = mdl.g.(λ, mdl.Y)
            z = mdl.g(λ, y)

            β = cholXΣX \ (mdl.X' * (cholΣ \ Z))
            q = ((Z - mdl.X * β)' * (cholΣ \ (Z - mdl.X * β)))[1]
            J = abs(prod(prime.(Ref(mdl.g), λ, mdl.Y)))
            m = (B * (cholΣ \ Z))[1] + (H * β)[1]

            h = q^(-(n - p) / 2) * J^(1 - p / n) / sqrt(det(cholΣ) * det(cholXΣX))
            T = LocationScale(m, sqrt(q * C), TDist(n - p))
            P = f(T, θ, λ, z)

            c += α * α′ * h
            acc += α * α′ * h * P
        end
    end

    return acc / c
end

"""
    btgdensity(mdl, x, y)

The probability density of the value `y` at the location `x` given a BTG model.

```math
\\mathcal{P}(y_x\\lvert m)
```
"""
function btgdensity(mdl::Model, x::Array{Float64, 1}, y::Float64, weights)
    return computedist(mdl, x, y, weights, function (T, θ, λ, z)
                       return pdf(T, z) * abs(prime(mdl.g, λ, y))
                       end)
end

"""
    btgdistribution(mdl, x, y)

The cumulative probability of the value `y` at the location `x` given a BTG model.

```math
\\Phi(y_x\\lvert m)
```
"""
function btgdistribution(mdl::Model, x::Array{Float64, 1}, y::Float64, weights)
    return computedist(mdl, x, y, weights, (T, _, _, z) -> cdf(T, z))
end

end # module

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Array Containers                  size
v = Vector{Int64}(undef, 5)         # 5
m = Matrix{Float64}(undef, 3, 4)    # 3, 4
a = Array{Real}(undef, 2, 2, 3)     # 2, 2, 3

# Tuple Containers -- immutable
t = ("hello", 3, 4.5)               # tuple
nt = (a="world", b=5, c='a')        # named tuple
println("keys ", keys(nt))
println("values ", values(t))

# Dictionary Containers -- mutable
d = Dict('a'=>3.4, "hello"=>t[2], 1=>nt.a)
println("keys ", keys(d))
println("values ", values(d))

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
## callbacks for integration
function condition(out, u, t, int)
    # end of step 1 when com == touchdown
    out[1] = pocmy(u, t, int.sol.prob.p) - pocmy(int.sol.prob.u0, 0.0, int.sol.prob.p)
    # end of step 2 when q2 == 0
    out[2] = u[2]
    # end of step 3 when pop2y == 0
    out[3] = pop2y(u, t, int.p)
end

function affect!(int, idx) # postive crossing (-ve to +ve)
    if idx == 1
        # do nothing
        nothing

    elseif idx == 2
        # do nothing
        nothing

    elseif idx == 3
        # do nothing
        nothing

    end
end

function affect_neg!(int, idx) # negative crossings (+ve to -ve)
    if idx == 1
        if int.sol.prob.p.virtual_force.flag == false
            # end of step 1
            # fit spline to force
            sol = int.sol
            T = int.t
            t = sol.t[1]:0.001:T
            splX1 = CubicSplineInterpolation(t, rx1(sol))
            splY1 = CubicSplineInterpolation(t, ry1(sol))
            splX2 = CubicSplineInterpolation(t, rx2(sol))
            splY2 = CubicSplineInterpolation(t, ry2(sol))

            println(T)
            println(t)
            println(sol.t)

            vrx1(t) = splX1(t - T)
            vry1(t) = splY1(t - T)
            vrx2(t) = splX2(t - T)
            vry2(t) = splY2(t - T)

            # update parameters
            sol.prob.p.virtual_force.flag = true
            sol.prob.p.virtual_force.vrx1 = vrx1
            sol.prob.p.virtual_force.vry1 = vry1
            sol.prob.p.virtual_force.vrx2 = vrx2
            sol.prob.p.virtual_force.vry2 = vry2

        end

    elseif idx == 2
        # end of step 2
        # int.t > 0.05 ? terminate!(int) : nothing # q2 starts above the ground

    elseif idx == 3
        # end of step 2
        # terminate!(int)

    end
end

vcb = VectorContinuousCallback(condition, affect!, affect_neg!, 3, save_positions=(false, false))

# pocmy and pop2y with u
function pocmy(u, t, p)
    @unpack l7, mc, l8, md, me, mf, mg, mh, mi, ma, mb, l10, l9, l6, l4, l2, l1, l11, l3, footang = p
    @unpack ea, fa, ha, ia, gs = p
    @inbounds q1, q2, q3, q4, q5, q6, q7 = u

    ea = ea(t)
    fa = fa(t)
    ha = ha(t)
    ia = ia(t)
    gs = gs(t)

    return ((q2 + ((l7 * mc + l8 * md + l8 * me + l8 * mf + l8 * mg + l8 * mh + l8 * mi) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7))) + ((l10 * me + l10 * mf + l10 * mg + l10 * mh + l10 * mi + l9 * md) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (cos(q6) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7))) + sin(q6) * (cos(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)) + sin(q3) * (cos(q6) * sin(q7) + sin(q6) * cos(q7)))) + ((l10 * mf + l10 * mh + me * (l10 - l9)) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) + ((l8 * mh + mf * (l8 - l7)) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (-(cos(fa)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - sin(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))) + ((mh * (0.5l6 + 0.5 * (l6 - l4))) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (sin(ha) * (sin(fa) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - cos(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))) - cos(ha) * (-(cos(fa)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - sin(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3)))) + ((mi * (l2 - l1)) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (sin(ia) * (-(cos(ha)) * (sin(fa) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - cos(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))) - sin(ha) * (-(cos(fa)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - sin(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3)))) - cos(ia) * (sin(ha) * (sin(fa) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - cos(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))) - cos(ha) * (-(cos(fa)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - sin(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))))) + ((l11 * mi + mg * gs) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * sin(q3) + 0.5 * ((l4 * mb + 2 * l6 * mc + 2 * l6 * md + 2 * l6 * me + 2 * l6 * mf + 2 * l6 * mg + 2 * l6 * mh + 2 * l6 * mi) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (cos(q4) * ((cos(q4) * cos(q5) - sin(q4) * sin(q5)) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7))) + (-(cos(q4)) * sin(q5) - sin(q4) * cos(q5)) * (cos(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)) + sin(q3) * (cos(q6) * sin(q7) + sin(q6) * cos(q7)))) + sin(q4) * ((cos(q4) * cos(q5) - sin(q4) * sin(q5)) * (cos(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)) + sin(q3) * (cos(q6) * sin(q7) + sin(q6) * cos(q7))) + (cos(q4) * sin(q5) + sin(q4) * cos(q5)) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7))))) + 0.5 * ((l3 * mb) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (((cos(q4) * cos(q5) - sin(q4) * sin(q5)) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7))) + (-(cos(q4)) * sin(q5) - sin(q4) * cos(q5)) * (cos(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)) + sin(q3) * (cos(q6) * sin(q7) + sin(q6) * cos(q7)))) * (cos(footang) * cos(q4) - sin(footang) * sin(q4)) + ((cos(q4) * cos(q5) - sin(q4) * sin(q5)) * (cos(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)) + sin(q3) * (cos(q6) * sin(q7) + sin(q6) * cos(q7))) + (cos(q4) * sin(q5) + sin(q4) * cos(q5)) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)))) * (cos(footang) * sin(q4) + sin(footang) * cos(q4)))) - ((l1 * ma + l2 * mb + l2 * mc + l2 * md + l2 * me + l2 * mf + l2 * mg + l2 * mh + l2 * mi) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * ((cos(q4) * cos(q5) - sin(q4) * sin(q5)) * (cos(q3) * (-(cos(q6)) * sin(q7) - sin(q6) * cos(q7)) + sin(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7))) + (-(cos(q4)) * sin(q5) - sin(q4) * cos(q5)) * (cos(q3) * (cos(q6) * cos(q7) - sin(q6) * sin(q7)) + sin(q3) * (cos(q6) * sin(q7) + sin(q6) * cos(q7))))) - 0.5 * ((l3 * mh) / (ma + mb + mc + md + me + mf + mg + mh + mi)) * (cos(footang) * (sin(ha) * (sin(fa) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - cos(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))) - cos(ha) * (-(cos(fa)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - sin(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3)))) - sin(footang) * (-(cos(ha)) * (sin(fa) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - cos(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3))) - sin(ha) * (-(cos(fa)) * (cos(ea) * sin(q3) - sin(ea) * cos(q3)) - sin(fa) * (cos(ea) * cos(q3) + sin(ea) * sin(q3)))))
end

function pop2y(u, t, p)
    @unpack l2 = p
    @inbounds q1, q2, q3, q4, q5, q6, q7 = u

    return q2 - l2 * sin(q3 - q4 - q5 - q6 - q7)

end

## multi-simulation callbacks
condition_step1(u, t, int) = pocmy(u, t, int.sol.prob.p) - pocmy(int.sol.prob.u0, 0.0, int.sol.prob.p)
affect_step1!(int) = nothing
affect_neg_step1!(int) = terminate!(int)
cb1 = ContinuousCallback(condition_step1, affect_step1!, affect_neg_step1!, save_positions=(false, true))

condition_step2(u, t, int) = u[2]
affect_step2!(int) = nothing
affect_neg_step2!(int) = int.t > 0.005 ? terminate!(int) : nothing # q2 starts above the ground
cb2 = ContinuousCallback(condition_step2, affect_step2!, affect_neg_step2!; save_positions=(false, true))


text:
@testset "DataSplits" begin
    @testset "calc_mask" begin
        pools = fill(:U, 5)
        pools[1] = :Lin
        pools[2] = :Lout
        init_mask = trues(5)
        init_mask[5] = false
        @test calc_mask(FullSplitStrat(), init_mask, pools) == init_mask
        @test calc_mask(UnlabeledSplitStrat(), init_mask, pools) == [false, false, true, true, false]
        @test calc_mask(UnlabeledAndLabeledInlierSplitStrat(), init_mask, pools) == [true, false, true, true, false]
        @test calc_mask(LabeledSplitStrat(), init_mask, pools) == [true, true, false, false, false]
        @test calc_mask(LabeledInlierSplitStrat(), init_mask, pools) == [true, false, false, false, false]
        @test calc_mask(LabeledOutlierSplitStrat(), init_mask, pools) == [false, true, false, false, false]
    end

    @testset "no holdout" begin
        data = reshape(collect(1:10), (2, 5))
        pools = fill(:U, 5)
        pools[1] = :Lin
        pools[2] = :Lout
        train = trues(5)
        ds = DataSplits(train, UnlabeledSplitStrat())
        train_data, train_pools, train_idx = get_train(ds, data, pools)
        @test size(train_data, 2) == length(train_pools) == length(train_idx)
        @test train_data == data[:, pools .== :U]
        @test all(train_pools .== :U)
        @test train_idx == [3, 4, 5]
    end

    @testset "with holdout" begin
        data = zeros(2, 5)
        pools = fill(:U, 5)
        train = BitArray([true, true, true, false, false])
        test = BitArray([false, false, false, true, true])
        @assert all(train .⊻ test)
        ds = DataSplits(train, test, FullSplitStrat())
        _, _, train_idx = get_train(ds, data, pools)
        _, _, test_idx = get_test(ds, data, pools)
        _, _, query_idx = get_query(ds, data, pools)
        @test train_idx == [1, 2, 3]
        @test test_idx == [4, 5]
        @test query_idx == [1, 2, 3]
    end
end

@testset "splits and init pools" begin
    data = zeros(2, 100)
    labels = fill(:inlier, 100)
    labels[1:10] .= :outlier
    @test_throws ArgumentError OneClassActiveLearning.get_splits_and_init_pools(data, labels, "foo", "Pu")
    @test_throws ArgumentError OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "foo")
    @testset "full split strategy" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pu")
        @test ds.train == ds.test == trues(size(data, 2))
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pp")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == size(data, 2) * 0.1
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pn")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pa", x=10)
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == size(data, 1) + 10
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sf", "Pa", x=5)
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == size(data, 1) + 5
    end
    @testset "80 20 holdout" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pu")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pp")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip[ds.train])) == sum(ds.train) * 0.1
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pn")
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip[ds.train])) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Sh", "Pa", x=10)
        @test (sum(ds.train) == 80) && (sum(ds.test) == 20) && all(ds.train .⊻ ds.test)
        @test length(findall(x -> x == :Lin, ip[ds.train])) == size(data, 1) + 10
    end
    @testset "inff" begin
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pu")
        @test ds.train == ds.test == trues(size(data, 2))
        @test ip == fill(:U, size(data, 2))
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pp")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == sum(ds.train) * 0.1
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pn")
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> (x == :Lin) || (x == :Lout), ip)) == 20
        ds, ip = OneClassActiveLearning.get_splits_and_init_pools(data, labels, "Si", "Pa", x=10)
        @test ds.train == ds.test == trues(size(data, 2))
        @test length(findall(x -> x == :Lin, ip)) == size(data, 1) + 10
    end
end


text:
const MOIU = MOI.Utilities
using LinearAlgebra

# Taken from https://github.com/JuliaOpt/MathOptInterface.jl/pull/473
"""
    set_dot(x::Vector, y::Vector, set::AbstractVectorSet)

Return the scalar product between a vector `x` of the set `set` and a vector
`y` of the dual of the set `s`.
"""
set_dot(x::Vector, y::Vector, set::MOI.AbstractVectorSet) = dot(x, y)

"""
    set_dot(x, y, set::AbstractScalarSet)

Return the scalar product between a number `x` of the set `set` and a number
`y` of the dual of the set `s`.
"""
set_dot(x, y, set::MOI.AbstractScalarSet) = dot(x, y)

_getconstant(set::MOI.AbstractScalarSet, T) = MOIU.getconstant(set)
_getconstant(set::MOI.Integer, T) = zero(T)

scalar_constant(T::Type, ::MOI.SingleVariable) = zero(T)
scalar_constant(::Type, f::MOI.AbstractScalarFunction) = MOI._constant(f)

function constraint_constant(model::MOI.ModelLike,
                             ci::MOI.ConstraintIndex{<:MOI.AbstractVectorFunction,
                                                     <:MOI.AbstractVectorSet},
                             T::Type)
    return MOIU._constant(MOI.get(model, MOI.ConstraintFunction(), ci))
end
function constraint_constant(model::MOI.ModelLike,
                             ci::MOI.ConstraintIndex{<:MOI.AbstractScalarFunction,
                                                     <:MOI.AbstractScalarSet},
                             T::Type)
    return scalar_constant(T, MOI.get(model, MOI.ConstraintFunction(), ci)) -
           _getconstant(MOI.get(model, MOI.ConstraintSet(), ci), T)
end

"""
    objective_bound(model::MOI.ModelLike,
                    F::Type{<:MOI.AbstractFunction},
                    S::Type{<:MOI.AbstractSet},
                    T::Type)

Return the part of `ObjectiveBound` due to the constraint of index `ci` using
scalar type `T`.
"""
function objective_bound(model::MOI.ModelLike,
                         ci::MOI.ConstraintIndex,
                         T::Type)
    return set_dot(constraint_constant(model, ci, T),
                   MOI.get(model, MOI.ConstraintDual(), ci),
                   MOI.get(model, MOI.ConstraintSet(), ci))
end

function objective_bound(model::MOI.ModelLike,
                         ci::MOI.ConstraintIndex{<:MOI.AbstractScalarFunction,
                                                 <:MOI.Interval},
                         T::Type)
    constant = scalar_constant(T, MOI.get(model, MOI.ConstraintFunction(), ci))
    set = MOI.get(model, MOI.ConstraintSet(), ci)
    dual = MOI.get(model, MOI.ConstraintDual(), ci)
    if dual < zero(dual)
        # The dual is negative so it is in the dual of the MOI.LessThan cone
        # hence the upper bound of the Interval set is tight
        constant -= set.upper
    else
        # the lower bound is tight
        constant -= set.lower
    end
    return set_dot(constant, dual, set)
end

"""
    objective_bound(model::MOI.ModelLike,
                    F::Type{<:MOI.AbstractFunction},
                    S::Type{<:MOI.AbstractSet},
                    T::Type)

Return the part of `ObjectiveBound` due to `F`-in-`S` constraints using scalar
type `T`.
"""
function objective_bound(model::MOI.ModelLike,
                         F::Type{<:MOI.AbstractFunction},
                         S::Type{<:MOI.AbstractSet},
                         T::Type)
    bound = zero(T) # sum won't work if there are now constraints
    for ci in MOI.get(model, MOI.ListOfConstraintIndices{F, S}())
        bound += objective_bound(model, ci, T)
    end
    return bound
end

function objective_bound(model::MOI.ModelLike,
                         F::Type{MOI.VectorOfVariables},
                         S::Type{<:MOI.AbstractVectorSet},
                         T::Type)
    # No constant in the function nor set so no contribution to the objective
    # bound
    return zero(T)
end

"""
    dual_objective_value(model::MOI.ModelLike, ::MOI.ObjectiveBound, T::Type)::T

Compute the objective value of the dual problem of type `T` using the
`ConstraintDual` results and the `ConstraintFunction` and `ConstraintSet` values.
"""
function dual_objective_value(model::MOI.ModelLike, T::Type)
    bound = zero(T) # sum will not work if there are zero constraints
    for (F, S) in MOI.get(model, MOI.ListOfConstraints())
        bound += objective_bound(model, F, S, T)::T
    end
    if MOI.get(model, MOI.ObjectiveSense()) != MOI.MAX_SENSE
        bound = -bound
    end
    dual_status = MOI.get(model, MOI.DualStatus())
    if dual_status == MOI.INFEASIBILITY_CERTIFICATE ||
        status == MOI.NEARLY_INFEASIBILITY_CERTIFICATE
        # The objective constant should not be present in rays
        F = MOI.get(model, MOI.ObjectiveFunctionType())
        f = MOI.get(model, MOI.ObjectiveFunction{F}())
        bound += scalar_constant(T, f)
    end
    return bound::T
end


text:
using CompileBot

# using runtests:
snoop_bot(BotConfig("TestPackage3", version = [v"1.0.5", v"1.2",v"1.5.0"], else_version = v"1.2", tmin =0.0))


text:
using Test

using Dice
using Dice: CuddMgr, DistBool

@testset "Sugar Tests" begin

    @test infer((@dice flip(0.5)), :bdd) ≈ 0.5
    @test infer((@dice DistBool(dicecontext(), 0.5)), :bdd) ≈ 0.5

    @test infer((@dice prob_equals(DistInt(42), 42)), :bdd) ≈ 1
    @test infer((@dice prob_equals(DistInt(dicecontext(), 42), 42)), :bdd) ≈ 1

end


text:
# module Reports

using Match

include("reports/balancesheet.jl")

include("export/json.jl")

export balancesheet, total

# end  # module Reports


text:
# This file was generated, do not modify it. # hide
@from_network machine(Deterministic(), X, y; predict=yhat) begin
    mutable struct MyTwoModelStack
        regressor1=model1
        regressor2=model2
        judge=judge
    end
end

my_two_model_stack = MyTwoModelStack()


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

miss20_dat = Feather.read("test/comparisons/dat_miss20_dat.feather")
miss30_dat = Feather.read("test/comparisons/dat_miss30_dat.feather")
miss50_dat = Feather.read("test/comparisons/dat_miss50_dat.feather")

miss20_par = Feather.read("test/comparisons/dat_miss20_par.feather")
miss30_par = Feather.read("test/comparisons/dat_miss30_par.feather")
miss50_par = Feather.read("test/comparisons/dat_miss50_par.feather")

miss20_par_mean = Feather.read("test/comparisons/dat_miss20_par_mean.feather")
miss30_par_mean = Feather.read("test/comparisons/dat_miss30_par_mean.feather")
miss50_par_mean = Feather.read("test/comparisons/dat_miss50_par_mean.feather")

three_path_dat = Feather.read("test/comparisons/three_path_dat.feather")

miss20_mat = Matrix(miss20_dat)
miss30_mat = Matrix(miss30_dat)
miss50_mat = Matrix(miss50_dat)

diff_fin = SemFiniteDiff(BFGS(), Optim.Options(;g_tol = 0.001))

diff_fin = SemFiniteDiff(LBFGS(
        ;alphaguess = LineSearches.InitialStatic(;scaled = true),
        linesearch = LineSearches.Static()), Optim.Options(;g_tol = 0.001))

## Model definition
@ModelingToolkit.variables x[1:31], mₓ, my[1:8], m[1:11]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

M = [mₓ, mₓ, mₓ, my[1:5]..., 3.0, my[7:8]..., 0.0, 0.0, 0.0]

M_free = [m[1:11]..., 0.0, 0.0, 0.0]

#S

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

start_val = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

start_val_mean = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3),
    [mean(Matrix(three_path_dat), dims = 1)...][[1, 4:8..., 10:11...]]
    )

start_val_mean_free = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3),
    [mean(Matrix(three_path_dat), dims = 1)...]
)

semobserved = SemObsMissing(miss20_mat)

imply = ImplySymbolic(A, S, F, [x..., m...], start_val_mean_free; M = M_free)
imply_mean = ImplySymbolic(A, S, F, [x..., mₓ, my[1:5]..., my[7:8]...], start_val_mean; M = M)

loss = Loss([SemFIML(semobserved, imply, [0.0], start_val_mean_free)])
loss_mean = Loss([SemFIML(semobserved, imply, [0.0], start_val_mean)])

model_fin = Sem(semobserved, imply, loss, diff_fin)
model_fin_mean = Sem(semobserved, imply_mean, loss_mean, diff_fin)

model_fin(start_val)

solution_fin = sem_fit(model_fin)
solution_fin_mean = sem_fit(model_fin_mean)

par_order = [collect(21:34); collect(15:20); 2;3; 5;6;7; collect(9:14);
    collect(35:45)]
all(
    abs.(solution_fin.minimizer .- miss20_par.est[par_order]
        ) .< 0.05*abs.(miss20_par.est[par_order]))

par_order_mean = [collect(25:38); collect(15:20); 2;3; 5;6;7; collect(9:14); 21;
    collect(39:45)]

all(
    abs.(solution_fin_mean.minimizer .- miss20_par_mean.est[par_order_mean]
        ) .< 0.05*abs.(miss20_par.est[par_order_mean]))

### 30% missingness

semobserved = SemObsMissing(miss30_mat)

loss = Loss([SemFIML(semobserved, [0.0], start_val_mean_free)])
loss_mean = Loss([SemFIML(semobserved, [0.0], start_val_mean)])

model_fin = Sem(semobserved, imply, loss, diff_fin)
model_fin_mean = Sem(semobserved, imply_mean, loss_mean, diff_fin)

solution_fin = sem_fit(model_fin)
solution_fin_mean = sem_fit(model_fin_mean)

all(
    abs.(solution_fin.minimizer .- miss30_par.est[par_order]
        ) .< 0.05*abs.(miss30_par_mean.est[par_order]))

all(
    abs.(solution_fin_mean.minimizer .- miss30_par_mean.est[par_order_mean]
        ) .< 0.05*abs.(miss30_par_mean.est[par_order_mean]))

### 50% missingness

semobserved = SemObsMissing(miss50_mat)

loss = Loss([SemFIML(semobserved, [0.0], start_val_mean_free)])
loss_mean = Loss([SemFIML(semobserved, [0.0], start_val_mean)])

model_fin = Sem(semobserved, imply, loss, diff_fin)
model_fin_mean = Sem(semobserved, imply_mean, loss_mean, diff_fin)

model_fin(start_val_mean_free)

solution_fin = sem_fit(model_fin)
solution_fin_mean = sem_fit(model_fin_mean)

all(
    abs.(solution_fin.minimizer .- miss50_par.est[par_order]
        ) .< 0.05*abs.(miss50_par.est[par_order]))

all(
    abs.(solution_fin_mean.minimizer .- miss50_par_mean.est[par_order_mean]
        ) .< 0.05*abs.(miss50_par_mean.est[par_order_mean]))


text:
export Jacobi

"""
    Jacobi()

Preconditioner object for the Jacobi (diagonal) left-side preconditioner.

Effective for diagonally-dominant systems.
"""
struct Jacobi{PCleft} <: AbstractPreconditionerType
    pc_side::PCleft
    function Jacobi()
        return new{PCleft}()
    end
end

mutable struct JacobiCache{AT} <: AbstractPreconditionerCache
    P_inv::AT
end

function cache(::Jacobi, linearoperator!, Q)
    x = similar(Q)
    AT = typeof(Q)
    return JacobiCache{AT}(x)
end

# Create P inverse preconditioner operator
function PCinitialize!(
    ::Jacobi,
    pc::Preconditioner,
    Q,
    Qrhs,
    args...,
)
    # extract matrix form of operator A
    n = length(Q)
    A = similar(Q, n, n)
    pc.linop!(A, I)
    # store inverse of diagonal preconditioner
    pc.cache.P_inv .= 1 ./ diag(A)
    return nothing
end

# Apply P inverse
function PCapply!(
    ::Jacobi,
    pc::Preconditioner,
    W,
    Q,
    args...,
)
    P_inv = pc.cache.P_inv
    W .= P_inv .* Q
end

max_stars_repo_name: CliMA/SystemSolvers

max_issues_repo_name: CliMA/SystemSolvers


text:
# Consider breaking this into two functions - one to create the problem and one to solve it.
function simulateMesocosm(;p = [1e-3, 1e3, 1.6, 0.0, 3.344, 1.0, 1.0], hist = 0.0, C₀ = [4.0], t₀ = 0.0, ⏰ = 1e3, reactionKinetics = "firstOrder", shape = "power law")
    # Define differential equation
    @sam_model(reactionKinetics, shape)
    # Make a tuple containing the start and end times of the simulation
    tspan = (t₀, ⏰)
    # Specify condition before simulation start
    h(p,t) = hist
    # Build the DDE problem
    prob = DDEProblem(sam_model!, C₀, h, tspan, p)
    # Solve problem
    sol = solve(prob, MethodOfSteps(RK4()))
    return(sol)
  end

max_stars_repo_name: elizabethjmohr/hyporheicBGC

max_issues_repo_name: elizabethjmohr/hyporheicBGC

max_forks_repo_name: elizabethjmohr/hyporheicBGC


text:
module ADCPtools

export hello, getzaxmap

hello(a::AbstractString) = "Hello, $a"

# function getzaxmap(rbeam::Float32, theta::Float32)
function getzaxmap(rbeam, theta)
# USAGE
# -----
# mapped_zaxis = getzaxmap(rbeam, theta)
#
# Returns the projection of the along-beam coordinate 'rbeam'
# from the local vertical axis. The along-beam direction is tilted
# 'theta' degrees from the vertical.
  Zmap = rbeam.*cos(theta.*pi/180)
  dZu = Zmap[end] - Zmap[end-1]
  rmax = rbeam[end]
  nz = length(Zmap)
  while Zmap[end]<(rmax-dZu)
    Zmap = [Zmap Zmap[end]+dZu] # Complete vertical axis to map to.
  end
end

end

max_stars_repo_name: apaloczy/ADCPtools

max_issues_repo_name: apaloczy/ADCPtools

max_forks_repo_name: apaloczy/ADCPtools


text:
@testset "aggregation" begin
    v = rand(5)
    @test aggregate(v, mav) ≈ mean(v)
    @test aggregate(v, TruePositive()) ≈ sum(v)
    @test aggregate(v, rms) ≈ sqrt(mean(v.^2))
    λ = rand()
    @test aggregate(λ, rms) === λ
    @test aggregate(aggregate(v, l2), l2) == aggregate(v, l2)
end

@testset "metadata" begin
    measures()
    measures(m -> m.target_scitype <: AbstractVector{<:Finite} &&
                  m.supports_weights)
end

@testset "coverage" begin
    # just checking that the traits work not that they're correct
    @test orientation(BrierScore()) == :score
    @test orientation(auc) == :score
    @test orientation(rms) == :loss

    @test reports_each_observation(auc) == false
    @test is_feature_dependent(auc) == false

    @test MLJBase.distribution_type(BrierScore{UnivariateFinite}) ==
        UnivariateFinite
end

true


text:
# Several devices together

export DeviceSet, MeasDataSet

mutable struct DeviceSet{DevList} <: AbstractDAQ
    "Device name associated to this `DeviceSet`"
    devname::String
    "Index of most relevant measurement device"
    iref::Int
    "List of devices"
    devices::DevList
    "Starting time of data acquisition"
    time::DateTime
    "Map from device name to device index in `devices`."
    devdict::OrderedDict{String,Int}
end

"""
`DeviceSet(dname, devices::DevList, iref=1)`

Create a meta device that handles data acquisition from several independent interfaces.

The `devices` argument specifies the individual DAQ devices that are used. It could
be a tuple (recommended) or a vector of AbstractDAQ.

The argument `iref` corresponds to the reference device (if it exists).
This reference device is simply the device that is used when checking if data acquisition
is ongoing or how many samples have been read.
"""
function DeviceSet(dname, devices::DevList, iref=1) where {DevList}

    devdict = OrderedDict{String,Int}()
    ndev = length(devices)
    for (i, dev) in enumerate(devices)
        devdict[devname(dev)] = i
    end

    return DeviceSet(dname, iref, devices, now(), devdict)
end

import Base.getindex

"""
`dev[i]`

Return the `i`-th device of a device set
"""
getindex(dev::DeviceSet, i) = dev.devices[i]
function getindex(devset::DeviceSet, dname::AbstractString)
    for dev in devset.devices
        if devname(dev) == dname
            return dev
        end
    end

    # If we got here, this the actuator doesn't exist
    # throw an exception:
    throw(KeyError(dname))
end

"""
`MeasDataSet(devname, devtype, time, data)`

Stores the data acquired by a `DeviceSet`.
"""
struct MeasDataSet <: AbstractMeasData
    "Device name"
    devname::String
    "Device type (`DeviceSet`)"
    devtype::String
    "Data acquisition time"
    time::DateTime
    "Data acquired by each device in the `DeviceSet`"
    data::OrderedDict{String,AbstractMeasData}
end

"""
`devname(d::MeasDataSet)`

Return the device name that acquired the data.
"""
devname(d::MeasDataSet) = d.devname

"""
`devtype(d::MeasDataSet)`

Return the device type that acquired the data ([`DeviceSet`](@ref) in this case).
"""
devtype(d::MeasDataSet) = d.devtype

"""
`meastime(d::MeasDataSet)`

Return the [`DateTime`](@ref) when the device started to acquire the data from
a [`DeviceSet`](@ref).
"""
meastime(d::MeasDataSet) = d.time

"""
`d["some/path/to/measurements"]`

Retrieve data acquired stored in a [`MeasDataSet`](@ref).

The data,  in this case, was acquired by several devices that
make up a [`DeviceSet](@ref). Thus sub-device might be another [`MeasDataSet`](@ref)
or, more commonly, a [`MeasData`](@ref) structure.

As an example, imagine that the [`DeviceSet`](@ref) is madeup of `dev1` and `dev2`
devices.

```
d["dev1"]
```

returns the data acquired by `dev1`. To get a specific channel, you can use

```
d["dev1/chanx"]
```

and this will return the value stored by channel `chanx` of `dev1`.

The channel can be specified independently as in the following example:

```
d["dev1", "chanx"]
```

or it can be specified by index:

```
d["dev1", 3]
```

(assuming `chanx` corresponds do channel 3)

In both of last cases, the indexing is forwarded to the `getindex` method for
[`MeasData`](@ref) for data retrieval.

"""
function getindex(d::MeasDataSet, path::String)
    println(path)
    p = split(path, '/')
    dev = p[1]
    if length(p) == 1
        return d.data[dev]
    else
        return d.data[dev][join(p[2:end], '/')]
    end

end

function getindex(d::MeasDataSet, path::String, idx...)
    p = split(path, '/')
    dev = p[1]
    dev = p[1]
    if length(p) == 1
        return getindex(d.data[dev], idx...)
    else
        return getindex(d.data[dev], join(p[2:end], '/'), idx...)
    end
end

"""
`daqchannels(d::MeasDataSet)`

Return channel names associated with each device that is acquiring data.
The device name is prepended to the channel name separated by a '/'.
"""
function daqchannels(d::MeasDataSet)
    chans = String[]
    for (k,v) in d.data
        devchans = daqchannels(v)
        for c in devchans
            push!(chans, k * "/" * c)
        end
    end
    return chans
end

"""
`daqstart(devs::DeviceSet)`

Start asynchrohous data acquisition on every device.
"""
function daqstart(devs::DeviceSet)
    devs.time = now()
    for dev in devs.devices
        daqstart(dev)
    end
    return
end

"""
`daqread(devs::DeviceSet)`

Read the data from every device in `DeviceSet`. It stores this data in a dictionary
where the key is the device name and the value is the data.
"""
function daqread(devs::DeviceSet)
    data = OrderedDict{String,MeasData}()

    for dev in devs.devices
        d = daqread(dev)
        data[devname(d)] = d
    end

    return MeasDataSet(devname(devs), "DeviceSet", devs.time, data)
end

"""
`daqacquire(devs::DeviceSet)`

Execute a synchronous data acquisition of every device.
"""
function daqacquire(devs::DeviceSet)
    daqstart(devs)
    return daqread(devs)
end

#import Base.getindex
#Base.getindex
samplesread(devs::DeviceSet) = samplesread(devs.devices[devs.iref])
isreading(devs::DeviceSet) = isreading(dev.devices[devs.iref])
isdaqfinished(devs::DeviceSet) = isdaqfinished(dev.devices[devs.iref])
issamplesavailable(devs::DeviceSet)=issamplesavailable(devs.devices[devs.iref])

"""
`savedaqdata(h5, devs::DeviceSet, data)`

Save the acquired data to a path inside a HDF5 file. It will save the data of each of the
devices.
"""
function savedaqdata(h5, devs::DeviceSet{T}, data; kw...) where {T}
    g = create_group(h5, devname(devs))
    attributes(g)["devtype"] = "DeviceSet"
    attributes(g)["devices"] = collect(keys(data))
    attributes(g)["time"] = time2ms(devs.time)
    for (k,v) in kw
        attributes(g)[string(k)] = v
    end

    for (k,v) in data
        savedaqdata(g, v)
    end
end

"""
`savedaqdata(h5, devs::DeviceSet, data)`

Save the acquired data to a path inside a HDF5 file. It will save the data of each of the
devices.
"""
function savedaqdata(h5, data::MeasDataSet; kw...) where {T}
    g = create_group(h5, data.devname)
    attributes(g)["devname"] = dev.devname
    attributes(g)["devtype"] = "DeviceSet"
    attributes(g)["devices"] = collect(keys(data.data))
    attributes(g)["time"] = time2ms(data.time)
    for (k,v) in kw
        attributes(g)[string(k)] = v
    end

    for (k,v) in data.data
        savedaqdata(g, v)
    end
end

"""
`savedaqconfig(h5, devs::DeviceSet)`

Saves the configuration of a `DeviceSet`. This configuration corresponds to the
configuration of each devince in `devs::DeviceSet`.
"""
function savedaqconfig(h5, devs::DeviceSet)

    ndevs = length(devs.devices)
    g = create_group(h5, devname(devs))
    for dev in devs.devices
        savedaqconfig(g, dev)
    end
end


text:
const default_attrs = (disabled = false, style = "", class = "", changes = 0)
const default_input_attrs = merge(default_attrs, (value = "", type = "text"))

function input(;
               oninput = js"function () {data.value = this.value;}",
               onchange = js"function () {data.changes = data.changes + 1;}",
               selector = "input.input",
               kwargs...)

    datavals = merge(default_input_attrs, values(kwargs))
    data = map(to_observable, datavals)

    attrs = Dict(key => js"data[$key]" for key in keys(data))
    attrs[:oninput] = oninput
    attrs[:onchange] = onchange
    template = js"{view: () => m($selector, $attrs)}"
    m = MithrilComponent(template, data)
    return MithrilWidget{:input}(m, data.value)
end

_parse(::Type{S}, x) where{S} = parse(S, x)
function _parse(::Type{Dates.Time}, x)
    segments = split(x, ':')
    length(segments) >= 2 && all(!isempty, segments) || return nothing
    h, m = parse.(Int, segments)
    Dates.Time(h, m)
end

function _string(x::Dates.Time)
    h = Dates.hour(x)
    m = Dates.minute(x)
    string(lpad(h, 2, "0"), ":", lpad(m, 2, "0"))
end
_string(x) = string(x)
_string(x::Colorant) = "#"*hex(x)
_string(::Nothing) = ""

for (func, typ, type) in [(:timepicker, :(Dates.Time), "time"), (:datepicker, :(Dates.Date), "date"), (:colorpicker, Colorant, "color")]
    @eval function $func(; value = nothing, type = "", kwargs...)
        (value isa AbstractObservable) || (value = Observable{Union{$typ, Nothing}}(value))
        f = _string
        g = t -> _parse($typ, t)
        string_value = Observable{String}(f(value[]))
        ObservablePair(value, string_value, f = f, g = g)
        m = input(; value = string_value, type = $type, kwargs...)
        MithrilWidget{$(Expr(:quote, func))}(component(m), value)
    end
end

function spinbox(; value = Observable{Union{Int, Nothing}}(nothing), type = "", kwargs...)
    if !isa(value, AbstractObservable)
        T = value === nothing ? Union{Int, Nothing} : Union{typeof(value), Nothing}
        value = Observable{T}(value)
    end
    jsvalue = Observable{Any}(value[])
    m = input(; value = jsvalue, type = "number", kwargs...)

    excluded = Function[]

    push!(excluded, on(value) do val
        str = jsvalue[]
        if !isa(str, AbstractString) || tryparse(Float64, str) != val
            jsvalue[notify = !in(excluded)] = val
        end
    end)

    push!(excluded, on(jsvalue) do val
        if val isa AbstractString
            res = tryparse(Float64, val)
            res !== nothing && (value[notify = !in(excluded)] = res)
        end
    end)

    MithrilWidget{:spinbox}(component(m), value)
end

const default_button_attrs = merge(default_attrs, (changes = 0,))

function button(children...; kwargs...)
    datavals = merge(default_button_attrs, values(kwargs))
    data = map(to_observable, datavals)

    attrs = Dict(key => js"data[$key]" for key in keys(data))
    attrs[:onclick] = js"function () {data.changes = data.changes + 1;}"

    template = node("button.button", children..., attributes = attrs)
    MithrilWidget{:button}(MithrilComponent(template, data), data.changes)
end


text:
"
 Test for projected newton QP solver

 author: Atsushi Sakai
"

using ArgMin
using Test

function test_solve_qp_with_projected_newton1()
    println("test_solve_qp_with_projected_newton1")
    n 		= 5
    g 		= randn(n)
    H 		= randn(n,n)
    H 		= H*H'
    lower 	= -ones(n)
    upper 	=  ones(n)

    xstar, status = solve_qp_with_projected_newton(
                                 H, g, lower, upper)

    @test status["optimal"] == true
end

function test_solve_qp_with_projected_newton2()
    println("test_solve_qp_with_projected_newton2")
    n 		= 3
    q = [3.0;4.0;-5.0]
    A = [1.0 0.0 -2.0;
         1.0 2.0 0.0;
         0.0 0.0 3.0]
    A = A*A'
    xmin = [-3.0;-5.0;-2.0]
    xmax = [1.0;1.0;1.0]

    xstar, status = solve_qp_with_projected_newton(
                                 A, q, xmin, xmax)
    @test xstar[1] ≈ 0.791667 atol=0.01
    @test xstar[2] ≈ -0.958333 atol=0.01
    @test xstar[3] ≈ 1.0 atol=0.01

end

test_solve_qp_with_projected_newton1()
test_solve_qp_with_projected_newton2()


text:
# B-spline basis function

@doc raw"""
B-spline basis function.
Right-sided limit version.
```math
\begin{aligned}
{B}_{(i,p,k)}(t)
&=
\frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t)
+\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\
{B}_{(i,0,k)}(t)
&=
\begin{cases}
    &1\quad (k_{i}\le t< k_{i+1})\\
    &0\quad (\text{otherwise})
\end{cases}
\end{aligned}
```
"""
function bsplinebasis₊₀(P::BSplineSpace, t::Real)::Vector{Float64}
    p = P.degree
    k = P.knots

    n = dim(P)
    if p == 0
        return [k[i] ≤ t < k[i+1] for i in 1:n]
    end
    K = [ifelse(k[i+p] == k[i], 0, (t - k[i]) / (k[i+p] - k[i])) for i in 1:n+1]
    B = bsplinebasis₊₀(BSplineSpace(p-1, k), t)
    return [K[i] * B[i] + (1 - K[i+1]) * B[i+1] for i in 1:n]
end

@doc raw"""
B-spline basis function.
Left-sided limit version.
```math
\begin{aligned}
{B}_{(i,p,k)}(t)
&=
\frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t)
+\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\
{B}_{(i,0,k)}(t)
&=
\begin{cases}
    &1\quad (k_{i}< t\le k_{i+1})\\
    &0\quad (\text{otherwise})
\end{cases}
\end{aligned}
```
"""
function bsplinebasis₋₀(P::BSplineSpace, t::Real)::Vector{Float64}
    p = P.degree
    k = P.knots

    n = dim(P)
    if p == 0
        return [k[i] < t ≤ k[i+1] for i in 1:n]
    end
    K = [ifelse(k[i+p] == k[i], 0, (t - k[i]) / (k[i+p] - k[i])) for i in 1:n+1]
    B = bsplinebasis₋₀(BSplineSpace(p-1, k), t)
    return [K[i] * B[i] + (1 - K[i+1]) * B[i+1] for i in 1:n]
end

@doc raw"""
B-spline basis function.
Modified version.
```math
\begin{aligned}
{B}_{(i,p,k)}(t)
&=
\frac{t-k_{i}}{k_{i+p}-k_{i}}{B}_{(i,p-1,k)}(t)
+\frac{k_{i+p+1}-t}{k_{i+p+1}-k_{i+1}}{B}_{(i+1,p-1,k)}(t) \\
{B}_{(i,0,k)}(t)
&=
\begin{cases}
    &1\quad (k_{i}\le t<k_{i+1}<k_{l})\\
    &1\quad (k_{i}\le t\le k_{i+1}=k_{l})\\
    &0\quad (\text{otherwise})
\end{cases}
\end{aligned}
```
"""
function bsplinebasis(P::BSplineSpace, t::Real)::Vector{Float64}
    p = P.degree
    k = P.knots

    n = dim(P)
    if p == 0
        return [k[i] ≤ t < k[i+1] || (k[i] ≠ k[i+1] == k[end] == t) for i in 1:n]
    end
    K = [ifelse(k[i+p] == k[i], 0, (t - k[i]) / (k[i+p] - k[i])) for i in 1:n+1]
    B = bsplinebasis(BSplineSpace(p-1, k), t)
    return [K[i] * B[i] + (1 - K[i+1]) * B[i+1] for i in 1:n]
end

"""
``i``-th B-spline basis function.
Right-sided limit version.
"""
function bsplinebasis₊₀(P::BSplineSpace, i::Integer, t::Real)::Float64
    p = P.degree
    k = P.knots

    if p == 0
        return k[i] ≤ t < k[i+1]
    else
        return (
            ((k[i+p] - k[i] ≠ 0) ? bsplinebasis₊₀(BSplineSpace(p-1, k), i, t) * (t - k[i]) / (k[i+p] - k[i]) : 0) +
            ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis₊₀(BSplineSpace(p-1, k), i+1, t) * (k[i+p+1] - t) / (k[i+p+1] - k[i+1]) : 0)
        )
    end
end

"""
``i``-th B-spline basis function.
Left-sided limit version.
"""
function bsplinebasis₋₀(P::BSplineSpace, i::Integer, t::Real)::Float64
    p = P.degree
    k = P.knots

    if p == 0
        return k[i] < t ≤ k[i+1]
    else
        return (
            ((k[i+p] - k[i] ≠ 0) ? bsplinebasis₋₀(BSplineSpace(p-1, k), i, t) * (t - k[i]) / (k[i+p] - k[i]) : 0) +
            ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis₋₀(BSplineSpace(p-1, k), i+1, t) * (k[i+p+1] - t) / (k[i+p+1] - k[i+1]) : 0)
        )
    end
end

"""
``i``-th B-spline basis function.
Modified version.
"""
function bsplinebasis(P::BSplineSpace, i::Integer, t::Real)::Float64
    p = P.degree
    k = P.knots

    if p == 0
        return k[i] ≤ t < k[i+1] || (k[i] ≠ k[i+1] == k[end] == t)
    else
        return (
            ((k[i+p] - k[i] ≠ 0) ? bsplinebasis(BSplineSpace(p-1, k), i, t) * (t - k[i]) / (k[i+p] - k[i]) : 0) +
            ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis(BSplineSpace(p-1, k), i+1, t) * (k[i+p+1] - t) / (k[i+p+1] - k[i+1]) : 0)
        )
    end
end

@doc raw"""
1st derivative of B-spline basis function.
Right-sided limit version.
```math
\dot{B}_{(i,p,k)}(t)
=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
```
"""
function bsplinebasis′₊₀(P::BSplineSpace, t)::Vector{Float64}
    p = P.degree
    k = P.knots

    n = dim(P)
    if p == 0
        return zeros(n)
    end
    K = [ifelse(k[i+p] == k[i], 0, p / (k[i+p] - k[i])) for i in 1:n+1]
    B = bsplinebasis₊₀(BSplineSpace(p-1, k), t)
    return [K[i] * B[i] - K[i+1] * B[i+1] for i in 1:n]
end

@doc raw"""
1st derivative of B-spline basis function.
Left-sided limit version.
```math
\dot{B}_{(i,p,k)}(t)
=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
```
"""
function bsplinebasis′₋₀(P::BSplineSpace, t)::Vector{Float64}
    p = P.degree
    k = P.knots

    n = dim(P)
    if p == 0
        return zeros(n)
    end
    K = [ifelse(k[i+p] == k[i], 0, p / (k[i+p] - k[i])) for i in 1:n+1]
    B = bsplinebasis₋₀(BSplineSpace(p-1, k), t)
    return [K[i] * B[i] - K[i+1] * B[i+1] for i in 1:n]
end

@doc raw"""
1st derivative of B-spline basis function.
Modified version.
```math
\dot{B}_{(i,p,k)}(t)
=p\left(\frac{1}{k_{i+p}-k_{i}}B_{(i,p-1,k)}(t)-\frac{1}{k_{i+p+1}-k_{i+1}}B_{(i+1,p-1,k)}(t)\right)
```
"""
bsplinebasis′

function bsplinebasis′(P::BSplineSpace, t)::Vector{Float64}
    p = P.degree
    k = P.knots

    n = dim(P)
    if p == 0
        return zeros(n)
    end
    K = [ifelse(k[i+p] == k[i], 0, p / (k[i+p] - k[i])) for i in 1:n+1]
    B = bsplinebasis(BSplineSpace(p-1, k), t)
    return [K[i] * B[i] - K[i+1] * B[i+1] for i in 1:n]
end

function bsplinebasis′₊₀(P::BSplineSpace, i::Integer, t::Real)::Float64
    p = P.degree
    k = P.knots

    return p * (
        ((k[i+p] - k[i] ≠ 0) ? bsplinebasis₊₀(BSplineSpace(p-1, k), i, t) / (k[i+p] - k[i]) : 0) -
        ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis₊₀(BSplineSpace(p-1, k), i+1, t) / (k[i+p+1] - k[i+1]) : 0)
    )
end

function bsplinebasis′₋₀(P::BSplineSpace, i::Integer, t::Real)::Float64
    p = P.degree
    k = P.knots

    return p * (
        ((k[i+p] - k[i] ≠ 0) ? bsplinebasis₋₀(BSplineSpace(p-1, k), i, t) / (k[i+p] - k[i]) : 0) -
        ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis₋₀(BSplineSpace(p-1, k), i+1, t) / (k[i+p+1] - k[i+1]) : 0)
    )
end

function bsplinebasis′(P::BSplineSpace, i::Integer, t::Real)::Float64
    p = P.degree
    k = P.knots

    return p * (
        ((k[i+p] - k[i] ≠ 0) ? bsplinebasis(BSplineSpace(p-1, k), i, t) / (k[i+p] - k[i]) : 0) -
        ((k[i+p+1] - k[i+1] ≠ 0) ? bsplinebasis(BSplineSpace(p-1, k), i+1, t) / (k[i+p+1] - k[i+1]) : 0)
    )
end

@doc raw"""
Return the support of ``i``-th B-spline basis function.
```math
\operatorname{supp}(B_{(i,p,k)})=[k_{i},k_{i+p+1}]
```
"""
function bsplinesupport(P::AbstractBSplineSpace, i::Integer)
    p = degree(P)
    k = knots(P)
    return k[i]..k[i+p+1]
end

function bsplinesupport(P::AbstractBSplineSpace)
    p = degree(P)
    k = knots(P)
    return [k[i]..k[i+p+1] for i in 1:dim(P)]
end

# For compatibility
bsplinebasis₊₀(i::Integer,P::AbstractBSplineSpace,t::Real) = bsplinebasis(P,i,t)
bsplinebasis₋₀(i::Integer,P::AbstractBSplineSpace,t::Real) = bsplinebasis(P,i,t)
bsplinebasis(i::Integer,P::AbstractBSplineSpace,t::Real) = bsplinebasis(P,i,t)
bsplinebasis′₊₀(i::Integer,P::AbstractBSplineSpace,t::Real) = bsplinebasis(P,i,t)
bsplinebasis′₋₀(i::Integer,P::AbstractBSplineSpace,t::Real) = bsplinebasis(P,i,t)
bsplinebasis′(i::Integer,P::AbstractBSplineSpace,t::Real) = bsplinebasis(P,i,t)
bsplinesupport(i::Integer,P::AbstractBSplineSpace) = bsplinesupport(P,i)


text:
#==============================================================================
    run a test on the REPL and get a "PASS" if its OK
    there is an existing @test macro, mine should be sweeter

    package:    DesertIslandDisk
    tested:     PASS
    depend:     NO
    author:     mahiki@
 =============================================================================#

"""
    @test

Return the variable, function, Symbol, object, et al name as a string.

@tostring(arg::Any)

# Examples
```jldoctest
julia> a = 55;
julia> bingo = "I'm the new sherrif";
```
"""
macro assert(ex)
           return :( $ex ? nothing : throw(AssertionError($(string(ex)))) )
       end

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
using Sundials
using Grid ## for interpolating grids

##
## Solving the 1D cable equation with IDA.
##
## The DAE system solved is a spatial discretization of the PDE
##
##  C_m dV/dt = d/R d^2V/x^2 - G_m V - R_m J
##
## The PDE is treated with central differences on a uniform M x 1
## grid. The values of u at the interior points satisfy ODEs, and
## equations u = 0 at the boundaries are appended, to form a DAE
## system of size M.
##

## Length of cable in um
const L = 4000

## Number of space steps to simulate
const dx = 10.0
const xsteps = round(Int, L/dx)

## Number of timesteps to simulate
const tf = 10.0
const dt = 0.1
const timesteps = round(Int, tf/dt)
const t = collect(0.0:dt:(dt*(timesteps-1)))

const d = 4.0 * 1e-4 ## Cable diameter in cm
const R_m = 2.5e11 ## Membrane resistance in Ohms/cm^2
const G_m = 1.0/R_m ## Membrane conductance
const R_i = 30.0  ## Longitudinal resistivity in Ohms*cm
const C_m = 1e-6  ## Membrane capacitance in F/cm^2

## Define the injected current J
const J = zeros(timesteps, xsteps)
const I_in = 0.1
## When to start injecting current
const I_delay = 1.0
## Duration of injected current
const I_dur = 5.0
J[max(1,round(Int,(I_delay)/dt)):round(Int,(I_delay+I_dur)/dt),round(Int,1*xsteps/2)] = I_in

const G_J = map(i -> CoordInterpGrid(t, vec(J[:,i]), 0.0, InterpQuadratic),1:xsteps)

##
## cableres: cable equation system residual function. This uses
## central differencing on the interior points, and includes algebraic
## equations for the boundary values.
##

function cableres(t, u, up, r)
    r[:] = u ## Initialize r to u, to take care of boundary equations.

    ## Loop over segments; set res = up - (central difference).
    for i = 2:(xsteps-2)
        loc = i
        r[loc] = C_m * up[loc] - (d/(4*R_i)) * (u[loc-1] + u[loc+1] -
                                                2.0 * u[loc]) +
                 G_m * u[loc] - R_m * (G_J[loc])[t]

    end

    return Sundials.CV_SUCCESS
end

function initial()
    u = zeros(xsteps)

    u[2:xsteps-2] = -60.0 ## initial value -60 mV

    id = ones(xsteps)

    up = zeros(xsteps)
    r  = zeros(xsteps)

    cableres(0.0, u, up, r)

    ## Copy -res into up to get correct interior initial up values.
    up[:] = -1.0 * r

    return (u,up,id)
end

function idabandsol(f::Function, y0::Vector{Float64}, yp0::Vector{Float64},
                    id::Vector{Float64}, t::Vector{Float64};
                    reltol::Float64=1e-4, abstol::Float64=1e-6)

    neq = length(y0)
    mem = Sundials.IDACreate()
    Sundials.@checkflag Sundials.IDAInit(mem, cfunction(Sundials.idasolfun, Cint,
                                         (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Ref{Function})),
                                         t[1], y0, yp0)
    Sundials.@checkflag Sundials.IDASetId(mem, id)
    Sundials.@checkflag Sundials.IDASetUserData(mem, f)
    Sundials.@checkflag Sundials.IDASStolerances(mem, reltol, abstol)
    Sundials.@checkflag Sundials.IDABand(mem, neq, xsteps-1, xsteps-1)
    ##Sundials.@checkflag Sundials.IDADense(mem, neq)
    rtest = zeros(neq)
    Sundials.@checkflag Sundials.IDACalcIC(mem, Sundials.IDA_YA_YDP_INIT, t[2])
    yres = zeros(Float64, length(y0), length(t))
    ypres = zeros(Float64, length(y0), length(t))
    yres[:, 1] = y0
    ypres[:, 1] = yp0
    y = copy(y0)
    yp = copy(yp0)
    tout = [0.0]
    for k in 2:length(t)
        Sundials.@checkflag Sundials.IDASolve(mem, t[k], tout, y, yp, Sundials.IDA_NORMAL)
        yres[:, k] = y
        ypres[:, k] = yp
    end
    return yres, ypres
end

u0, up0, id = initial()

yout, ypout = @time idabandsol(cableres, u0, up0, id, map(x -> x, t),
                               reltol = 1e-3, abstol = 1e-4)


text:
lap = PowerLAPACK.LAPACK

println("Test stev!")

function test_stev(N::Integer, nt::Integer)
    dv = -2ones(N)
    ev = ones(N-1)
    dv_w = copy(dv)
    ev_w = copy(ev)

    println("$(N)×$(N) symmetric tridiagonal matrix, $nt iterations")

    println("eigfact:")
    T = SymTridiagonal(copy(dv), copy(ev))
    ee = 1
    function highlevel(nn)
        for i = 1:nn
            TT = copy(T) # Ugly, I know
            ee = Base.Eigen(Base.LAPACK.stev!('V', TT.dv, TT.ev)...)
        end
    end
    highlevel(1)
    @time highlevel(nt)

    println("stev! w/ prealloc'd work, passed as splat array:")
    s_work = lap.stev!('V', dv, ev, true)
    function lowlevel_compound(nn)
        for i = 1:nn
            dv_w[:] = dv
            ev_w[:] = ev
            lap.stev!('V', dv_w, ev_w, s_work...)
        end
    end
    lowlevel_compound(1)
    @time lowlevel_compound(nt)
    @test sumabs2(ee[:values]-dv_w) == 0
    @test sumabs2(ee[:vectors]-s_work[1]) == 0

    println("stev! w/ prealloc'd work, passed separately:")
    Zmat,work,info = lap.stev!('V', dv, ev, true)
    function lowlevel_separate(nn)
        for i = 1:nn
            dv_w[:] = dv
            ev_w[:] = ev
            lap.stev!('V', dv_w, ev_w, Zmat, work, info)
        end
    end
    lowlevel_separate(1)
    @time lowlevel_separate(nt)
    @test sumabs2(ee[:values]-dv_w) == 0
    @test sumabs2(ee[:vectors]-Zmat) == 0
end

function test_stev_subspace(N,M)
    dv = -2ones(N)
    ev = ones(N-1)

    println("$(N)×$(N) symmetric tridiagonal matrix, $(M)×$(M) subspace")

    T = SymTridiagonal(copy(dv[1:M]), copy(ev[1:M-1]))
    ee = Base.Eigen(Base.LAPACK.stev!('V', T.dv, T.ev)...)

    Zmat,work,info = lap.stev!('V', dv, ev, true)
    @time lap.stev!('V', sub(dv, 1:M), sub(ev, 1:M), sub(Zmat, 1:M, 1:M), work, info)
    dv = -2ones(N)
    ev = ones(N-1)
    @time lap.stev!('V', sub(dv, 1:M), sub(ev, 1:M), sub(Zmat, 1:M, 1:M), work, info)

    @test sumabs2(ee[:values]-dv[1:M]) == 0
    @test sumabs2(ee[:vectors]-Zmat[1:M,1:M]) == 0
end

test_stev(500,10)
test_stev_subspace(20,10)

println()


text:
using Turing, MCMCChains, Test

@testset "DataFrame constructor tests" begin
    @model gdemo(x) = begin
        m ~ Normal(1, 0.01)
        s ~ Normal(5, 0.01)
    end

    model = gdemo([1.5, 2.0])
    sampler = HMC(1000, 0.01, 5)

    chns = [sample(model, sampler, save_state=true) for i in 1:4]
    chn = chainscat(chns...)

    df = DataFrame(chn, showall=true)
    @test size(df) == (4000, 8)

    df1 = DataFrame(chn, [:parameters])
    @test size(df1) == (4000, 2)

    df2 = DataFrame(chn, [:internals, :parameters])
    @test size(df2) == (4000, 8)

    df3 = DataFrame(chn[:s])
    @test size(df3) == (4000, 1)

    df4 = DataFrame(chn[:lp])
    @test size(df4) == (4000, 1)

    df5 = DataFrame(chn, [:parameters], append_chains=false)
    @test size(df5) == (4, )
    @test size(df5[1]) == (1000, 2)

    df6 = DataFrame(chn, [:parameters, :internals], append_chains=false)
    @test size(df6) == (4, )
    @test size(df6[1]) == (1000, 8)

    df7 = DataFrame(chn, [:parameters, :internals],
        remove_missing_union=false)
    @test size(df7) == (4000, 8)

    df8 = DataFrame(chn, [:parameters, :internals],
        remove_missing_union=false, append_chains=false)
    @test size(df8) == (4, )
    @test size(df8[1]) == (1000, 8)
end


text:
# already in TumorGrowth: remove
# function sampletumor_mfreqs(sampletumor)
#     mutations = unique( vcat(sampletumor.mutations...) )
#     m_ind = Dict(mutations .=> 1:length(mutations))
#     freqs = zeros(length(mutations))
#     for row in eachrow(sampletumor)
#         for (m,f) in zip(row.mutations, row.frequencies)
#             freqs[m_ind[m]] += f
#         end
#     end
#     freqs./=nrow(sampletumor)
#     return DataFrame(mutation = mutations, frequency=freqs)
# end

# function bisection(f, interval, iter; precision=1e-9)
#     xlow, xup = interval
#     sign(f(xlow)) == sign(f(xup)) && error("no zero in interval $interval")
#     any((isnan.(f.(interval)))) && error("f not defined on interval $interval")
#     for i=1:iter
#         xm = (xlow+xup)/2
#         abs(f(xm)) < precision && return (xm, i)
#         if sign(f(xm)) != sign(f(xlow))
#             xup = xm
#         else
#             xlow = xm
#         end
#     end
#     println(xup)
#     error("no convergence within $iter steps")
# end

function get_turnover(tumorinfo; useknown_N = false, useknown_T = false, Nthresh_orph, Nthresh_estr, tumor_sample_func = df -> df, mut_freqs_func = TumorGrowth.mutation_freqs, subsample_func = df -> df)
    clade_turnover = Float64[]
    clone_turnover = Float64[]

    @showprogress for sim in eachrow(tumorinfo)

        tumor, mutations, b, d= tumor_sample_func( sim.tumor ), sim.mutations, sim.b, sim.d

        htypes = unique(tumor.mutations)
        survivors = unique!(vcat(htypes...))
        freqs = if useknown_N
            Dict( survivors .=> 1 ./ mutations.N_birth[survivors] )
        elseif useknown_T
            Dict( survivors .=> @. 1 / exp( (b-d)*mutations.t_birth[survivors] ) )
        else
            mut_freqs_func(tumor) |> seq -> Dict(seq.mutation .=> seq.frequency)
        end

        orphaned_tumor = DataFrame( mutations = unique( filter.(m-> freqs[m] > 1/Nthresh_orph, htypes) ) ) |> subsample_func
        estranged_tumor = DataFrame(mutations = filter( muts -> all(  freqs[m] > 1/Nthresh_estr for m in muts), htypes) ) |> subsample_func

        W_a = Turnover.orphaned_red_treeless(orphaned_tumor) |> df -> sum(df.isorphaned)/sum(df.isgreen)

        W_o = Turnover.estranged_treeless(estranged_tumor) |> df -> sum(df.isestranged)/sum(df.isgreen)

        push!(clade_turnover, W_a)
        push!(clone_turnover, W_o)

        sleep(0.01)
    end
    return (ds = tumorinfo.d, Wa = clade_turnover, Wo = clone_turnover)
end

function infer_params( tumorinfo; Nthresh_orph, Nthresh_estr, Wa, Wo)
    dfits = []
    mufits = []

    @showprogress for i in 1:nrow(tumorinfo)
        b, mu = tumorinfo.b[i], tumorinfo.μ[i]
        W_a, W_o = Wa[i], Wo[i]

        d_solve = min(1., 2*log(Nthresh_orph)*W_a)*b

        mu_solve = missing
        try
#             mu_solve, n = bisection(x -> min(1., Turnover.W_estranged(d_solve; b=b, μ=x, T=log(Nthresh_estr)/(b-d_solve)))- W_o, [0.01,0.99], 100)
            mu_solve = fzero(x -> min(1., Turnover.W_estranged(d_solve; b=b, μ=x, T=log(Nthresh_estr)/(b-d_solve)))- W_o, 0.01, 0.99)
            catch e
        end

        push!(dfits, d_solve)
        push!(mufits, mu_solve)

        sleep(0.01)
    end
    return (ds = tumorinfo.d, dfits = dfits, mufits = mufits)
end

# function plot_infresult(ds, dfits, mufits; mu, size=(800,400))

#     p = plot(layout=(1,2), size=size, legend=:none, margin=3Plots.mm, yguidefontrotation=-90,
#              aspect_ratio=1, xaxis=(L"d", (0,1), 0.:0.2:1.), yaxis=((0,1),0:0.2:1), guidefontsize=15, tickfont=15, grid=false)

#     scatter!(p[1], ds, dfits, ylab=L"d_\mathrm{fit}", marker = (5, 0.2, :darkblue), alpha=0.4)
#     plot!(p[1], 0:1,0:1, lw=2, c=:red)
#     scatter!(p[1], unique(ds), [mean(filter(!ismissing, dfits[ds .== d])) for d in unique(ds)], yerror=[std(filter(!ismissing, dfits[ds .== d])) for d in unique(ds)], marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)

#     scatter!(p[2], ds, mufits, ylab=L"\mu_\mathrm{fit}", marker = (5, 0.2, :darkblue), alpha=0.4)
#     hline!(p[2], [mu], lw=2, c=:red)

#     mus = [filter(!ismissing, mufits[ds .== d]) for d in unique(ds)]
#     mask = .!isempty.(mus)

#     scatter!(p[2], unique(ds)[mask], mean.(mus[mask]),
#         yerror= std.(mus[mask]), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)
# end

# function plot_turnover(ds, Wa, Wo; Nthresh_orph, Nthresh_estr, b=1., mu, plotargs...)

#     p = plot(layout=(1,2), size=(600,250), legend=:none, margin=3Plots.mm, xlab=L"d", xlim=(0,1))

#     scatter!(p[1], ds, Wa, ylab=L"W_{a}", alpha=0.2, ms = 3, c=:black, ylim=(0,0.3))
#     scatter!(p[1], unique(ds), [mean(filter(!ismissing, Wa[ds .== d])) for d in unique(ds)],
#         yerror=[std(filter(!ismissing, Wa[ds .== d])) for d in unique(ds)], ms=10., marker=:hline, c=:black)
#     plot!(p[1], 0.:0.01:b, d-> Turnover.W_orphaned(d/b; N=Nthresh_orph))

#     scatter!(p[2], ds, Wo, ylab=L"W_{o}", alpha=0.5, ms = 3, c=:black, ylim=(0,1))

#     Wo_bins = [filter(!ismissing, Wo[ds .== d]) for d in unique(ds)]
#     mask = .!isempty.(Wo_bins)

#     scatter!(p[2], unique(ds)[mask], mean.(Wo_bins[mask]),
#         yerror= std.(Wo_bins[mask]), ms=10., marker=:hline, c=:black)
#     plot!(p[2], 0.:0.01:b, d-> min( 1, Turnover.W_estranged(d; b=b, μ=mu, T=log(Nthresh_estr)/(b-d))) )
# end

function plot_turnover_violin(ds, Wa, Wo; Nthresh_orph, Nthresh_estr, b=1., mu, plotargs...)
    scalex=5
    d_uni = unique(ds)
    bins = [findall(isequal(d), ds) for d in d_uni]

    p = plot(; layout=(1,2), size=(600,250), legend=:none, margin=3Plots.mm, xlab=L"d", xticks=(0:scalex, 0.:1/scalex:1.), plotargs...)

    scatter!(p[1], ds*scalex, Wa, ylab=L"W_{a}", alpha=0.2, ms = 3, c=:darkblue, ylim=(0,0.3))
    violin!(p[1], ds*scalex, Wa, marker = (5, 0.2, :darkblue), alpha=0.4, c=:lightblue)
    scatter!(p[1], d_uni*scalex, [median(Wa[bin]) for bin in bins], marker=(:hline, 3*scalex, :darkblue) )
    plot!(p[1], (0.:0.01:b)*scalex, Turnover.W_orphaned.((0.:0.01:b)/b; N=Nthresh_orph))

    scatter!(p[2], ds*scalex, Wo, ylab=L"W_{o}", alpha=0.5, ms = 3, c=:black, ylim=(0,1))
    violin!(p[2], ds*scalex, Wo, marker = (5, 0.2, :darkblue), alpha=0.4, c=:lightblue)
    scatter!(p[2], d_uni*scalex, [median(Wo[bin]) for bin in bins], marker=(:hline, 3*scalex, :darkblue))
    plot!(p[2], (0.:0.01:b)*scalex, d -> min(1, Turnover.W_estranged(d/scalex/b; b=b, μ=mu, T=log(Nthresh_estr)/(b-d/scalex))) )
end

function plot_infresult_violin(ds, dfits, mufits; mu, size=(600,300), plotargs...)
    scalex=5
    d_uni = unique(ds)
    bins = [findall(isequal(d), ds) for d in d_uni]

    p = plot(; layout=(1,2), size=size, legend=:none, margin=3Plots.mm, yguidefontrotation=-90,
             aspect = scalex, xlab=L"d", yaxis=((0,1),0:0.2:1), xticks=(0:scalex, 0.:1/scalex:1.), plotargs...)

    scatter!(p[1], ds*scalex, dfits, ylab=L"d_\mathrm{fit}", marker = (5, 0.2, :darkblue), alpha=0.4)
    violin!(p[1], ds*scalex, dfits, marker = (5, 0.2, :darkblue), alpha=0.4, c=:lightblue)
    scatter!(p[1], d_uni*scalex, [median(dfits[bin]) for bin in bins], marker=(:hline, 3*scalex, :darkblue) )
    plot!(p[1], [0, scalex], [0,1], lw=2, c=:red)

    scatter!(p[2], ds*scalex, mufits, ylab=L"\mu_\mathrm{fit}", marker = (5, 0.2, :darkblue), alpha=0.4)
    violin!(p[2], ds*scalex, mufits, marker = (5, 0.2, :darkblue), alpha=0.4, c=:lightblue)
    scatter!(p[2], d_uni*scalex, [median(mufits[bin]) for bin in bins], marker=(:hline, 3*scalex, :darkblue))
    hline!(p[2], [mu], lw=2, c=:red)
end

function plot_series_violin(ds, dfits, mufits; mu, size=(800,300))
    bins = sort!(unique(ds))
    dfits_binned = [ [dfits[findall(isequal(d), ds)]...] for d in bins]
    mufits_binned = [ [mufits[findall(isequal(d), ds)]...] for d in bins]
    length.(dfits_binned) |> println

    scalex = 5
    p = plot(layout=(1,3), legend=:none, yguidefontrotation=-90, size=size,
             aspect_ratio=scalex, xaxis=(L"b", (0,scalex)), yaxis=((0,1),0:0.2:1), guidefontsize=15, tickfont=15, xticks=(0:scalex, range(0,1,length=scalex+1)))

    plot!(p[1], 0:scalex,range(0,1,length=scalex+1), c=:red, lw=2.)
    violin!(p[1], ds*scalex, dfits, ylab=L"b_{inf}", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[1], ds*scalex, dfits, ylab=L"b_{inf}", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[1], bins*scalex, mean.(dfits_binned), yerror = std.(dfits_binned), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)

    hline!(p[2], [mu], c=:red, lw=2.)
    violin!(p[2], ds*scalex, mufits, ylab=L"\mu_{inf}", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[2], ds*scalex, mufits, ylab=L"\mu_{inf}", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[2], bins*scalex, mean.(mufits_binned), yerror = std.(mufits_binned), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)

    mubeta = [mufits_binned[i]./(1 .- dfits_binned[i]) for i =1:length(bins)]
    plot!(p[3], 0:scalex,range(0,1,length=scalex+1), c=:red, lw=2., xlab=L"\mu/(1-b/a)", yguidefontrotation=-0)
    violin!(p[3], mu./(1 .-ds)*scalex, mufits./(1 .-dfits), ylab=L"\mu_{inf}/(1-b_{inf}/a)", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[3], mu./(1 .-ds)*scalex, mufits./(1 .-dfits), ylab=L"\mu_{inf}/(1-b_{inf}/a)", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[3], mu./(1 .-bins)*scalex, mean.(mubeta), yerror = std.(mubeta), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)
end

function plot_series_scatter(ds, dfits, mufits; mu, size=(800,300))
    bins = sort!(unique(ds))
    dfits_binned = [ [dfits[findall(isequal(d), ds)]...] for d in bins]
    mufits_binned = [ [mufits[findall(isequal(d), ds)]...] for d in bins]
    length.(dfits_binned) |> println

    p = plot(layout=(1,3), legend=:none, yguidefontrotation=-90, size=size,
             aspect_ratio=1, xaxis=(L"b", (0,1), 0.:0.2:1.), yaxis=((0,1),0:0.2:1), guidefontsize=15, tickfont=15)

    plot!(p[1], 0:1,0:1, c=:red, lw=2.)
    scatter!(p[1], ds, dfits, ylab=L"b_{inf}", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[1], bins, mean.(dfits_binned), yerror = std.(dfits_binned), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)

    hline!(p[2], [mu], c=:red, lw=2.)
    scatter!(p[2], ds, mufits, ylab=L"\mu_{inf}", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[2], bins, mean.(mufits_binned), yerror = std.(mufits_binned), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)

    mubeta = [mufits_binned[i]./(1 .- dfits_binned[i]) for i =1:length(bins)]
    plot!(p[3], 0:1, 0:1, c=:red, lw=2., xlab=L"\mu/(1-b/a)", yguidefontrotation=-0)
    scatter!(p[3], mu./(1 .-ds), mufits./(1 .-dfits), ylab=L"\mu_{inf}/(1-b_{inf}/a)", marker = (5, 0.2, :darkblue), alpha=0.4)
    scatter!(p[3], mu./(1 .-bins), mean.(mubeta), yerror = std.(mubeta), marker = (:hline, 14, :darkblue), markerstrokecolor=:darkblue)
end

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
using Documenter, GlobalSensitivity

makedocs(
    sitename="GlobalSensitivity.jl",
    authors="Vaibhav Kumar Dixit",
    clean=true,
    doctest=false,
    modules=[GlobalSensitivity],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://globalsensitivity.sciml.ai/stable/"),

    pages=[
        "GlobalSensitivity.jl: Global Sensitivity Analysis (GSA)" => "index.md",
        "Methods" => Any[
          "methods/morris.md",
          "methods/sobol.md",
          "methods/regression.md",
          "methods/efast.md"
        ],
        "Examples" => Any[
          "examples/l_k_global.md",
          "examples/design_matrices.md",
          "examples/parallelized_gsa.md"
        ],
    ]
)

deploydocs(
    repo="github.com/SciML/GlobalSensitivity.jl";
    push_preview=true
)


text:
### ============== ### ============ ###
##    Extended Inertial Spin Model   ##
##    Martin Zumaya Hernandez        ##
##    EXAMPLE SIMULATION SCRIPT      ##
### ============== ### ============ ###

### ============ INCLUDE PACKAGES ============ ###

using CollectiveDynamics.InertialSpin

### =============== ### =============== ###
###   DEFINITION OF INITIAL PARAMETERS  ###
### =============== ### =============== ###

# N   = 128
# η   = 1.0
# τ   = 6
# rep = 1
# T   = 8*exp10(-5)
# δ   = 0.35

N    = parse(Int, ARGS[1]) # Number of particles
η    = parse(Float64, ARGS[2]) # Dissipation term
T    = parse(Float64, ARGS[3]) # temperature, noise
n_nl = parse(Float64, ARGS[4]) # non local interactions per particle
τ_i  = parse(Int, ARGS[5]) # number of iterations
τ_f  = parse(Int, ARGS[6]) # number of iterations
rep  = parse(Int, ARGS[7]) # repetition number

χ   = 1.25
J   = 0.8
ρ   = 0.3
v0  = 0.1
n_t = 6

pars = InertialParameters(χ, J, η, v0*sqrt(J/χ), ρ, v0, N, T, n_t, n_nl)

# σ = sqrt((2pars.d) * η * T) # noise std deviation ( square root of variance )
σ = sqrt((2pars.d) * η * T * pars.dt) # noise std deviation ( square root of variance )
L = cbrt(N / pars.ρ) # 3D

### =============== ### =============== ###
### SET UP SYSTEM AND OUTPUT STRUCTURE  ###
### =============== ### =============== ###

flock = InertialExtFlock(N, L, pars.v0)

output_path = set_output_data_structure_lr("EXTENDED_INERTIAL_SPIN_EXT", N, η, T, n_nl)

pos_file  = open(output_path * "/pos_$(rep).dat", "w+")
vel_file  = open(output_path * "/vel_$(rep).dat", "w+")
spin_file = open(output_path * "/spin_$(rep).dat", "w+")

### ============== ### ============== ### ============== ###
### INITIALIZATION FROM FILE
### ============== ### ============== ### ============== ###

# pos_data = reshape(raw_data, 3N, div(length(raw_data), 3N))
# vel_data = reshape(raw_data, 3N, div(length(raw_data), 3N))

data_path = joinpath(homedir(),"art_DATA","EXTENDED_INERTIAL_SPIN", "DATA","data_N_$(N)", "eta_$(η)", "eta_$(η)_T_$(T)_nl_$(n_nl)")

raw_data = reinterpret(Float64, read(joinpath(data_path,"pos_$(rep).dat")))
pos_data = raw_data[(end-3N+1):end]

raw_data = reinterpret(Float64, read(joinpath(data_path,"vel_$(rep).dat")))
vel_data = raw_data[(end-3N+1):end]

raw_data = reinterpret(Float64, read(joinpath(data_path,"spin_$(rep).dat")))
spin_data = raw_data[(end-3N+1):end]

for i in 1:N
    flock.pos[i]  = [pos_data[3(i-1)+1], pos_data[3(i-1)+2], pos_data[3(i-1)+3]]
    flock.vel[i]  = [vel_data[3(i-1)+1], vel_data[3(i-1)+2], vel_data[3(i-1)+3]]
    flock.spin[i] = [vel_data[3(i-1)+1], vel_data[3(i-1)+2], vel_data[3(i-1)+3]]
end

# write initial conditions
println("//////// ", 1)
write(pos_file, vcat(flock.pos...))
write(vel_file, vcat(flock.vel...))
write(spin_file, vcat(flock.spin...))

### ============== ### ============== ###
###          SYSTEM EVOLUTION         ###
### ============== ### ============== ###

times = [convert(Int, exp10(i)) for i in τ_i:τ_f]

for i in 1:(length(times) - 1)

    for t in (times[i]+1):times[i+1]

        evolve_extended_system(flock.pos, flock.vel, flock.v_t, flock.v_nl, flock.spin, flock.Rij, pars, σ)

        if t % times[i] == 0 || t % div(times[i], exp10(1)) == 0
            println("//////// ", t)
            write(pos_file, vcat(flock.pos...))
            write(vel_file, vcat(flock.vel...))
            write(spin_file, vcat(flock.spin...))
        end
    end

end

close(pos_file)
close(vel_file)
close(spin_file)

println("Done all")


text:
#=
# Infernal.jl usage example
=#

import Infernal
import Rfam

# Retrieve the path to Infernal programs

isfile(Infernal.path("cmalign"))

# Execute a command

success(Infernal.exe("cmfetch", Rfam.cm(), "RF00162"))


text:
@testset "UInt256" begin

    x = reinterpret(UInt256, UInt8[
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
        0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
    ])[1];

    @test typeof(x) == UInt256
    @test x == zero(x)
    @test x == zero(UInt256)
    @test bswap(x) == x

    @test string(x, base = 16) == "0000000000000000000000000000000000000000000000000000000000000000"

    y = UInt256("0000000000000000000000000000000000000000000000000000000000000001")
    @test y % UInt8 == 0x01
end


text:
function Position(arg0::AbsoluteDate, arg1::Vector3D, arg2::Vector{Vector{jdouble}}, arg3::jdouble, arg4::ObservableSatellite)
    return Position((AbsoluteDate, Vector3D, Vector{Vector{jdouble}}, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4)
end

function Position(arg0::AbsoluteDate, arg1::Vector3D, arg2::Vector{jdouble}, arg3::jdouble, arg4::ObservableSatellite)
    return Position((AbsoluteDate, Vector3D, Vector{jdouble}, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4)
end

function Position(arg0::AbsoluteDate, arg1::Vector3D, arg2::jdouble, arg3::jdouble, arg4::ObservableSatellite)
    return Position((AbsoluteDate, Vector3D, jdouble, jdouble, ObservableSatellite), arg0, arg1, arg2, arg3, arg4)
end

function get_correlation_coefficients_matrix(obj::Position)
    return jcall(obj, "getCorrelationCoefficientsMatrix", Vector{Vector{jdouble}}, ())
end

function get_covariance_matrix(obj::Position)
    return jcall(obj, "getCovarianceMatrix", Vector{Vector{jdouble}}, ())
end

function get_position(obj::Position)
    return jcall(obj, "getPosition", Vector3D, ())
end


text:
using Formulas
using Base.Test

## Associative property
@test Term(:(a+(b+c))) == Term(:(a+b+c))
@test Term(:((a+b)+c)) == Term(:(a+b+c))
@test Term(:(a&(b&c))) == Term(:(a&b&c))
@test Term(:((a&b)&c)) == Term(:(a&b&c))

## Distributive property
@test Term(:(a & (b+c))) == Term(:(a&b + a&c))
@test Term(:((a+b) & c)) == Term(:(a&c + b&c))
@test Term(:((a+b) & (c+d))) == Term(:(a&c + a&d + b&c + b&d))
@test Term(:(a & (b+c) & d)) == Term(:(a&b&d + a&c&d))

## Expand * to main effects + interactions
@test Term{:+}(Term(:(a*b))) == Term(:(a+b+a&b))
@test sort!(Term{:+}(Term(:(a*b*c)))) == Term(:(a+b+c+a&b+a&c+b&c+a&b&c))

## printing Terms:
@test string(Term(:a)) == "a"
@test string(Term(:(a+b))) == "+(a, b)"
@test string(Term(:(a + a&b))) == "+(a, &(a, b))"
@test string(Term(:(a+b | c))) == "(+(a, b) | c)"


text:
using JLD2, FileIO
using ClimateMachine
using ClimateMachine.SingleStackUtils
using ClimateMachine.Checkpoint
using ClimateMachine.BalanceLaws: vars_state
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
import CLIMAParameters

include(joinpath(clima_dir, "experiments", "AtmosLES", "stable_bl_model.jl"))
include(joinpath("helper_funcs", "diagnostics_configuration.jl"))
include("edmf_model.jl")
include("edmf_kernels.jl")

CLIMAParameters.Planet.T_surf_ref(::EarthParameterSet) = 265
CLIMAParameters.Atmos.EDMF.a_surf(::EarthParameterSet) = 0.0

"""
    init_state_prognostic!(
            turbconv::EDMF{FT},
            m::AtmosModel{FT},
            state::Vars,
            aux::Vars,
            localgeo,
            t::Real,
        ) where {FT}

Initialize EDMF state variables.
This method is only called at `t=0`.
"""
function init_state_prognostic!(
    turbconv::EDMF{FT},
    m::AtmosModel{FT},
    state::Vars,
    aux::Vars,
    localgeo,
    t::Real,
) where {FT}
    # Aliases:
    gm = state
    en = state.turbconv.environment
    up = state.turbconv.updraft
    N_up = n_updrafts(turbconv)
    # GCM setting - Initialize the grid mean profiles of prognostic variables (ρ,e_int,q_tot,u,v,w)
    z = altitude(m, aux)

    # SCM setting - need to have separate cases coded and called from a folder - see what LES does
    # a thermo state is used here to convert the input θ to e_int profile
    e_int = internal_energy(m, state, aux)
    param_set = parameter_set(m)
    ts = PhaseDry(param_set, e_int, state.ρ)
    T = air_temperature(ts)
    p = air_pressure(ts)
    q = PhasePartition(ts)
    θ_liq = liquid_ice_pottemp(ts)

    a_min = turbconv.subdomains.a_min
    @unroll_map(N_up) do i
        up[i].ρa = gm.ρ * a_min
        up[i].ρaw = gm.ρu[3] * a_min
        up[i].ρaθ_liq = gm.ρ * a_min * θ_liq
        up[i].ρaq_tot = FT(0)
    end

    # initialize environment covariance with zero for now
    if z <= FT(250)
        en.ρatke =
            gm.ρ *
            FT(0.4) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0)
        en.ρaθ_liq_cv =
            gm.ρ *
            FT(0.4) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0) *
            FT(1 - z / 250.0)
    else
        en.ρatke = FT(0)
        en.ρaθ_liq_cv = FT(0)
    end
    en.ρaq_tot_cv = FT(0)
    en.ρaθ_liq_q_tot_cv = FT(0)
    return nothing
end;

function main(::Type{FT}, cl_args) where {FT}

    surface_flux = cl_args["surface_flux"]

    # DG polynomial order
    N = 4
    nelem_vert = 15

    # Prescribe domain parameters
    zmax = FT(400)

    t0 = FT(0)

    # Simulation time
    timeend = FT(60)
    CFLmax = FT(0.50)

    config_type = SingleStackConfigType

    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    N_updrafts = 1
    N_quad = 3 # Using N_quad = 1 leads to norm(Q) = NaN at init.
    turbconv = EDMF(
        FT,
        N_updrafts,
        N_quad,
        param_set,
        surface = NeutralDrySurfaceModel{FT}(param_set),
    )

    model = stable_bl_model(
        FT,
        config_type,
        zmax,
        surface_flux;
        turbconv = turbconv,
    )

    # Assemble configuration
    driver_config = ClimateMachine.SingleStackConfiguration(
        "SBL_EDMF",
        N,
        nelem_vert,
        zmax,
        param_set,
        model;
        hmax = FT(40),
        solver_type = ode_solver_type,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFLmax,
    )

    # --- Zero-out horizontal variations:
    vsp = vars_state(model, Prognostic(), FT)
    horizontally_average!(
        driver_config.grid,
        solver_config.Q,
        varsindex(vsp, :turbconv),
    )
    horizontally_average!(
        driver_config.grid,
        solver_config.Q,
        varsindex(vsp, :energy, :ρe),
    )
    vsa = vars_state(model, Auxiliary(), FT)
    horizontally_average!(
        driver_config.grid,
        solver_config.dg.state_auxiliary,
        varsindex(vsa, :turbconv),
    )
    # ---

    dgn_config =
        config_diagnostics(driver_config, timeend; interval = "10ssecs")

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            (turbconv_filters(turbconv)...,),
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    diag_arr = [single_stack_diagnostics(solver_config)]
    time_data = FT[0]

    # Define the number of outputs from `t0` to `timeend`
    n_outputs = 5
    # This equates to exports every ceil(Int, timeend/n_outputs) time-step:
    every_x_simulation_time = ceil(Int, timeend / n_outputs)

    cb_data_vs_time =
        GenericCallbacks.EveryXSimulationTime(every_x_simulation_time) do
            diag_vs_z = single_stack_diagnostics(solver_config)

            nstep = getsteps(solver_config.solver)
            # Save to disc (for debugging):
            # @save "bomex_edmf_nstep=$nstep.jld2" diag_vs_z

            push!(diag_arr, diag_vs_z)
            push!(time_data, gettime(solver_config.solver))
            nothing
        end

    check_cons = (
        ClimateMachine.ConservationCheck("ρ", "3000steps", FT(0.001)),
        ClimateMachine.ConservationCheck("energy.ρe", "3000steps", FT(0.1)),
    )

    cb_print_step = GenericCallbacks.EveryXSimulationSteps(100) do
        @show getsteps(solver_config.solver)
        nothing
    end

    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        check_cons = check_cons,
        user_callbacks = (cbtmarfilter, cb_data_vs_time, cb_print_step),
        check_euclidean_distance = true,
    )

    diag_vs_z = single_stack_diagnostics(solver_config)
    push!(diag_arr, diag_vs_z)
    push!(time_data, gettime(solver_config.solver))

    return solver_config, diag_arr, time_data
end

# add a command line argument to specify the kind of surface flux
# TODO: this will move to the future namelist functionality
sbl_args = ArgParseSettings(autofix_names = true)
add_arg_group!(sbl_args, "StableBoundaryLayer")
@add_arg_table! sbl_args begin
    "--surface-flux"
    help = "specify surface flux for energy and moisture"
    metavar = "prescribed|bulk|custom_sbl"
    arg_type = String
    default = "custom_sbl"
end

cl_args = ClimateMachine.init(
    parse_clargs = true,
    custom_clargs = sbl_args,
    output_dir = get(ENV, "CLIMATEMACHINE_SETTINGS_OUTPUT_DIR", "output"),
    fix_rng_seed = true,
)

solver_config, diag_arr, time_data = main(Float64, cl_args)

## Uncomment lines to save output using JLD2
# output_dir = @__DIR__;
# mkpath(output_dir);
# function dons(diag_vs_z)
#     return Dict(map(keys(first(diag_vs_z))) do k
#         string(k) => [getproperty(ca, k) for ca in diag_vs_z]
#     end)
# end
# get_dons_arr(diag_arr) = [dons(diag_vs_z) for diag_vs_z in diag_arr]
# dons_arr = get_dons_arr(diag_arr)
# println(dons_arr[1].keys)
# z = get_z(solver_config.dg.grid; rm_dupes = true);
# save(
#     string(output_dir, "/sbl_edmf.jld2"),
#     "dons_arr",
#     dons_arr,
#     "time_data",
#     time_data,
#     "z",
#     z,
# )

include(joinpath(@__DIR__, "report_mse_sbl_edmf.jl"))

nothing


text:
__precompile__(false) # see https://github.com/samo-lin/ParallelStencil.jl/issues/39

module WhereTheWaterAlsoFlows

export flow_routing2D

const USE_GPU = false  # Use GPU? If this is set false, then the CUDA packages do not need to be installed! :)
const GPU_ID  = 0
using ParallelStencil
using ParallelStencil.FiniteDifferences2D

@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    CUDA.device!(GPU_ID) # select GPU
else
    @init_parallel_stencil(Threads, Float64, 2)
end

using Plots, Printf, Statistics, LinearAlgebra
################################
function save_infos(Lx, Ly, nx, ny, err, isave, outd)
    fid=open("$(outd)/0_infos.inf", "w"); @printf(fid,"%f %f %d %d %d %1.3e", Lx, Ly, nx, ny, isave, err); close(fid)
end
@static if USE_GPU
    function save_array(A, A_name, isave, outd)
        A_tmp = Array(A); A_tmp = convert(Array{Float32}, A_tmp)
        fid=open("$(outd)/$(isave)_$(A_name).res", "w"); write(fid, A_tmp); close(fid)
    end
else
    function save_array(A, A_name, isave, outd)
        A_tmp = convert(Array{Float32}, A)
        fid=open("$(outd)/$(isave)_$(A_name).res", "w"); write(fid, A_tmp); close(fid)
    end
end
################################
@parallel function def_err!(errD::Data.Array, D::Data.Array)

    @all(errD) = @all(D)

    return
end

@parallel function chk_err!(errD::Data.Array, D::Data.Array)

    @all(errD) = @all(D) - @all(errD)

    return
end

@parallel function mask_D!(D::Data.Array, Mask::Data.Array)

    @all(D) = @all(D) * @all(Mask)

    return
end

@parallel function compute_ϕ!(ϕ::Data.Array, rdh::Data.Number, rr::Data.Number, Zb::Data.Array, D::Data.Array, H::Data.Array)

    @all(ϕ) = @all(Zb) + rdh*@all(D) + rr*@all(H)

    return
end

@parallel_indices (ix,iy) function compute_flux!(Ux::Data.Array, Uy::Data.Array, qDx::Data.Array, qDy::Data.Array, k::Data.Number, _dx::Data.Number, _dy::Data.Number, ϕ::Data.Array, D::Data.Array)

    if (ix<=size(Ux,1)  && iy<=size(Ux,2))    Ux[ix,iy] = -k*_dx*(ϕ[ix+1,iy]-ϕ[ix,iy]) end
    if (ix<=size(Uy,1)  && iy<=size(Uy,2))    Uy[ix,iy] = -k*_dy*(ϕ[ix,iy+1]-ϕ[ix,iy]) end
    if (ix<=size(qDx,1) && iy<=size(qDx,2))  qDx[ix,iy] = Ux[ix,iy]*(Ux[ix,iy]>0.0)*D[ix,iy] + Ux[ix,iy]*(Ux[ix,iy]<0.0)*D[ix+1,iy  ] end
    if (ix<=size(qDy,1) && iy<=size(qDy,2))  qDy[ix,iy] = Uy[ix,iy]*(Uy[ix,iy]>0.0)*D[ix,iy] + Uy[ix,iy]*(Uy[ix,iy]<0.0)*D[ix  ,iy+1] end

    return
end

@parallel function compute_∂D!(∂D::Data.Array, _dx::Data.Number, _dy::Data.Number, M::Data.Number, qDx::Data.Array, qDy::Data.Array)

    @inn(∂D) = -(_dx*@d_xi(qDx) + _dy*@d_yi(qDy)) + 1.0/M

    return
end

@parallel function update_D!(∂Ddτ::Data.Array, D::Data.Array, damp::Data.Number, dτ::Data.Number, ∂D::Data.Array)

    @all(∂Ddτ) =  damp*@all(∂Ddτ) + @all(∂D)
    @all(D)    =  max(0.0, @all(D) + dτ*@all(∂Ddτ))

    return
end

@parallel_indices (ix,iy) function set_BCx!(D::Data.Array, val_BC::Data.Number)

    if (ix==1         && iy<=size(D,2)) D[ix,iy] = val_BC end
    if (ix==size(D,1) && iy<=size(D,2)) D[ix,iy] = val_BC end

    return
end

@parallel_indices (ix,iy) function set_BCy!(D::Data.Array, val_BC::Data.Number)
    # fixed height
    if (ix<=size(D,1) && iy==1        ) D[ix,iy] = val_BC end
    if (ix<=size(D,1) && iy==size(D,2)) D[ix,iy] = val_BC end
    # no flux
    # if (ix<=size(D,1) && iy==1        ) D[ix,iy] = D[ix,iy+1] end
    # if (ix<=size(D,1) && iy==size(D,2)) D[ix,iy] = D[ix,iy-1] end
    return
end

@parallel function postprocess!(Flux::Data.Array, qDx::Data.Array, qDy::Data.Array)

    @inn(Flux) = sqrt(@av_xi(qDx)*@av_xi(qDx) .+ @av_yi(qDy)*@av_yi(qDy))

    return
end

# physics scales
const ρ̂i  = 910
const ρ̂w  = 1000
const ĝ   = 9.81
const ϕ̂   = 1e7
const x̂   = 1e6
const M̂   = 1e-6
const D̂   = 1
# physics dependent
const t̂   = D̂/M̂
const Ĥ   = ϕ̂/ρ̂w/ĝ
const Ψ̂   = ϕ̂/x̂
const û   = x̂/t̂
const k̂   = û/Ψ̂
const rr  = ρ̂i/ρ̂w
const rdh = D̂/Ĥ    # fixed ice thikness routing
const s2d = 3600*24 # seconds in a day

"""

Inputs:
- xc, yc: center coordinates (as a range)
- Zb, H: bed elevation, ice thickness
- Mask: mask for input data
- D0: IC for water layer thickness (default ==1)
"""
@views function flow_routing2D(xc::AbstractRange, yc::AbstractRange, Zb, H, Mask=zero(Zb).+1, D0=zero(Zb).+1;
                               plotyes=false, saveyes=true, outdir="")

    nx     = length(xc)
    ny     = length(yc)

    @assert (nx, ny) == size(Zb) == size(H) == size(D0) "Sizes don't match"
    # fastest if multiple of 16 (as no overlength here)
    if USE_GPU && (rem(nx,32)!=0 || rem(ny,8)!=0)
        @warn "Gridpoints not divisible by 32, 8 [x, y], this leads to performance drop."
    end

    Lx, Ly = xc[end]-xc[1], yc[end]-yc[1]

    # nondim
    Lx, Ly = Lx/x̂, Ly/x̂
    xc, yc = xc/x̂, yc/x̂
    Zb     = Data.Array(Zb/Ĥ)
    H      = Data.Array(H/Ĥ)
    D      = Data.Array(D0/D̂)
    Mask   = Data.Array(Mask)
    # parameters
    k      = 100.0
    M      = M̂*s2d
    # numerics
    nsave  = 5e4
    itrMax = 3e4 #1e8
    nout   = 1000
    nmax   = 100
    ε      = 1e-7     # aboslute tolerance
    damp   = 0.9  #0.9 for nx=100
    # dτ_sc  = nx/128*30.0 #60 for nx,ny=256 #30 for nx,ny=128 #24 for nx,ny=96 #12.0 for nx,ny=64
    # dτ_sc  = nx/128*3.0 #60 for nx,ny=256 #30 for nx,ny=128 #24 for nx,ny=96 #12.0 for nx,ny=64
    dτ_sc  = 2.1 #60 for nx,ny=256 #30 for nx,ny=128 #24 for nx,ny=96 #12.0 for nx,ny=64
    dx, dy = Lx/nx, Ly/ny
    _dx, _dy = 1.0/dx, 1.0/dy
    # D      = max(Zb)-Zb + 10;
    ∂D     = @zeros(nx  ,ny  )
    ϕ      = @zeros(nx  ,ny  )
    Ux     = @zeros(nx-1,ny  )
    Uy     = @zeros(nx  ,ny-1)
    qDx    = @zeros(nx-1,ny  )
    qDy    = @zeros(nx  ,ny-1)
    ∂Ddτ   = @zeros(nx  ,ny  )
    errD   = @zeros(nx  ,ny  )
    Flux   = @zeros(nx  ,ny  )
    # preprocess
    xcp = Array(xc*x̂/1e3); ycp = Array(yc*x̂/1e3)
    save_array(xcp, "xcp", 0, outdir)
    save_array(ycp, "ycp", 0, outdir)
    save_array(Mask, "M" , 0, outdir)
    @parallel mask_D!(D, Mask)
    # action
    iter = 1; err = 2*ε;  max_U = 1.0; err1 = []; isave = 1
    while err > ε && iter < itrMax
        @parallel def_err!(errD, D)
    	# flow routing physics
    	@parallel compute_ϕ!(ϕ, rdh, rr, Zb, D, H)
    	@parallel compute_flux!(Ux, Uy, qDx, qDy, k, _dx, _dy, ϕ, D)
    	@parallel compute_∂D!(∂D, _dx, _dy, M, qDx, qDy)
    	if mod(iter,nmax)==0 || iter==1
    		max_U = max(maximum(abs.(Ux)), maximum(abs.(Uy)))
    	end
    	dτ = min(dx,dy)/max_U/dτ_sc
    	@parallel update_D!(∂Ddτ, D, damp, dτ, ∂D)
    	@parallel set_BCx!(D, 0.0)
    	@parallel set_BCy!(D, 0.0)
        @parallel mask_D!(D, Mask)
    	# Check errs
    	if mod(iter, nout)==0
    		@parallel chk_err!(errD, D)
    		push!(err1, mean(abs.(errD[:]))); err = err1[end]
            @printf("iter=%d  errD=%1.3e \n", iter, err1[end])
    	end
        if mod(iter, nsave)==0 || iter==1
            @parallel postprocess!(Flux, qDx, qDy)
            save_array(Flux, "F", isave, outdir)
            save_infos(Lx, Ly, nx, ny, err, isave, outdir)
            isave+=1
        end
    	iter+=1
    end
    # postprocess
    @parallel postprocess!(Flux, qDx, qDy)
    save_array(Flux, "F", isave, outdir)
    save_infos(Lx, Ly, nx, ny, err, isave, outdir)

    if plotyes
        Zb[D.==0.0] .= NaN
        D[D.==0.0]  .= NaN
        # ploting
        p1 = heatmap(xcp, ycp, Zb'*Ĥ, aspect_ratio=1, xlims=(xcp[1], xcp[end]), ylims=(ycp[1], ycp[end]), c=:inferno, title="Zb")
        p2 = heatmap(xcp, ycp, D'*D̂, aspect_ratio=1, xlims=(xcp[1], xcp[end]), ylims=(ycp[1], ycp[end]), c=:inferno, title="D")
        display(plot(p1, p2))
        # savefig(plot(p1, p2, p3, layout = l), joinpath(@__DIR__, "../output/o$nx.png"))
    end

    return Array(D*D̂), Array(Flux), xcp, ycp
end

end