text:
function timeseriesindices(t, identifiers...)
    N = length(t)
    @assert all(x->length(x)==N,identifiers)
    combinedIds = collect(zip(identifiers...))
    uniqueIds = unique(combinedIds)
    repInd = indexin(combinedIds, uniqueIds) # assign the same index to all samples in the same replicate

    ind = Vector{Int}[findall(repInd.==i) for i=1:length(uniqueIds)] # each entry is a vector with all indices for that replicate
    ind = map(x->x[sortperm(t[x])], ind) # sort each time series by increasing time
    ind, uniqueIds
end

_indexdim(N,dim,ind) = [ i==dim ? ind : Colon() for i=1:N ] # utility to put : as index for every dimension except the chosen one

function timeseries(X, timeSeriesIndices, dim=1)
    @assert dim<=ndims(X)
    map( ind->X[_indexdim(ndims(X),dim,ind)...], timeSeriesIndices )
end

function _interpolate(tStart, tEnd, A, B, t)
    t==tStart && return A
    t==tEnd && return B
    α = (t-tStart)/(tEnd-tStart)
    (1-α)*A + α*B
end

function timeseriesat(keyTimes::AbstractVector{T},V::AbstractMatrix{S},t::Number; clamp=false) where {T,S}
    @assert all(diff(keyTimes).>0)
    if clamp
        t = Base.clamp(t,keyTimes[1],keyTimes[end])
    else
        @assert keyTimes[1] <= t <= keyTimes[end]
    end

    ind = searchsortedlast(keyTimes,t)
    keyTimes[ind]   == t && return V[:,ind]   # avoids unnecessary interpolation and handles edge case
    keyTimes[ind+1] == t && return V[:,ind+1] # avoids unnecessary interpolation
    _interpolate(keyTimes[ind],keyTimes[ind+1], V[:,ind], V[:,ind+1], t)
end

timeseriesat(keyTimes::AbstractVector,V::AbstractMatrix,t::AbstractVector{T}; clamp=false) where {T<:Number} =
    hcat( [timeseriesat(keyTimes,V,s;clamp=clamp) for s in t]... )

timeseriesat(keyTimes::AbstractVector{T},V::AbstractVector{S},i,t; kwargs...) where {T<:AbstractVector,S<:AbstractMatrix} =
    timeseriesat(keyTimes[i],V[i],t; kwargs...)

function timeseriesmean(keyTimes::AbstractVector,V::AbstractMatrix)
    w = diff(keyTimes)/2
    w = vcat(w,0)+vcat(0,w)
    w = w/sum(w)
    sum(V.*w'; dims=2)
end
timeseriesmean(keyTimes::AbstractVector{T},V::AbstractVector{S}) where {T<:AbstractVector,S<:AbstractMatrix} = map(timeseriesmean,keyTimes,V)

# outputs d(V1,V2)(t), the distance of the multivariate time series V1 and V2 (DxT1, DxT2 matrices) as a function of time
# output given as a piecewise 2nd degree polynomial, 1st row is constant factor, 2nd row is linear factor, 3rd row is square factor (the polynomial describes squared distance)
# each piece is paramterized for the interval [0, tEnd-tStart]
function distancecurve(t1::AbstractVector{T},t2::AbstractVector{T},V1::AbstractMatrix{S},V2::AbstractMatrix{S}) where {T,S}
    #@assert issorted(t1) && issorted(t2)
    @assert all(diff(t1).>0) && all(diff(t2).>0) # check that the values are strictyly increasing

    # do not deal with degenerate cases
    @assert maximum(t1)>minimum(t2) && maximum(t2)>minimum(t1)
    @assert length(t1)>1 && length(t2)>1

    t = Vector{T}()
    params = Vector{S}()

    i=j=1 # index in t1 and t2
    while true
        tStart = max(t1[i], t2[j])
        tEnd = min(t1[i+1], t2[j+1])

        # handle interval
        if tEnd>tStart
            # starting vertices
            A1 = _interpolate(t1[i],t1[i+1],V1[:,i],V1[:,i+1],tStart)
            A2 = _interpolate(t2[j],t2[j+1],V2[:,j],V2[:,j+1],tStart)

            # ending vertices
            B1 = _interpolate(t1[i],t1[i+1],V1[:,i],V1[:,i+1],tEnd)
            B2 = _interpolate(t2[j],t2[j+1],V2[:,j],V2[:,j+1],tEnd)

            # l: P+sQ   # s∈[0,tEnd-tStart]
            P = A2-A1   # difference in starting pos
            Q = B2-B1-P # difference in direction
            Q = Q/(tEnd-tStart) # rescale s.t. parameter ranges from 0 to tEnd-tStart

            # ||l(s)||² = ||P+sQ||² = t²||Q||² + 2t<P,Q> + ||P||² = at² + bt + c
            a = dot(Q,Q)
            b = 2*dot(P,Q)
            c = dot(P,P)

            push!(t, tStart)
            push!(params, c, b, a)
        end

        # move one or both of i and j forward
        t1[i+1]==tEnd && (i+=1)
        t2[j+1]==tEnd && (j+=1)

        if i>=length(t1) || j>=length(t2)
            push!(t,tEnd)
            break
        end
    end

    params = reshape(params, 3, length(t)-1)
    t, params
end
distancecurve(t1::AbstractVector,t2::AbstractVector,V1::AbstractMatrix,V2::AbstractMatrix) = distancecurve(promote(t1,t2)...,promote(V1,V2)...)

# find all values of t for which f(t)==y, where f(t) is the square root of a positive piecewise 2nd-degree polynomial
function distancecurveequals(t::AbstractVector{T}, P::AbstractVecOrMat, y) where {T}
    @assert y>=0
    y = y^2

    tOut = Vector{promote_type(Float64,T)}()

    P[1,1]==y && push!(tOut, t[1]) # edge case

    for i=1:length(t)-1
        # solve in this segment
        c,b,a = (P[:,i]...,)
        c -= y

        if abs(a)<1e-16
            # linear case;  bx+c=0

            abs(b)<1e-16 && continue # constant case, ignore solutions since would consist of the entire interval

            x = -c/b

            # check that the solution is within the interval
            0<x<=t[i+1]-t[i] && push!(tOut, t[i]+x)
        else
            # solve 2nd degree polynomial
            w = b^2-4a*c
            w<0 && continue # no solution in this segment

            xmid = -b/(2a)
            xpm  = sqrt(w)/(2a)

            x1 = xmid - xpm
            x2 = xmid + xpm

            # println(i," ",x1," ", x2)

            # check that the solutions are within the interval
            0<x1<=t[i+1]-t[i] &&        push!(tOut, t[i]+x1)
            0<x2<=t[i+1]-t[i] && w>0 && push!(tOut, t[i]+x2) # avoid dubble root for w==0
        end
    end
    tOut
end

# first time the distance reaches a given threshold
function bifurcationtime(t::AbstractVector{T}, P::AbstractVecOrMat, y) where {T}
    RT = promote_type(Float64,T)
    y^2<=P[1,1] && return convert(RT,t[1])
    tEq = distancecurveequals(t,P,y)
    isempty(tEq) ? convert(RT,Inf) : tEq[1] # first time
end

function evaluatedistancecurve(t::AbstractVector, P::AbstractVecOrMat, x::AbstractVector{T}) where {T}
    @assert issorted(t)
    @assert issorted(x)
    @assert size(P,2)==length(t)-1

    y = convert(T,NaN)*ones(x)

    xStartInd = searchsortedfirst(x,t[1])

    # for each segment
    for i=1:length(t)-1
        tStart = t[i]
        tEnd = t[i+1]
        xEndInd = searchsortedlast(x,tEnd)

        # evaluate
        if xEndInd>=xStartInd
            xl = x[xStartInd:xEndInd] - tStart #
            y[xStartInd:xEndInd] = sqrt( [ones(xl) xl xl.^2]*P[:,i] )
        end

        xStartInd = xEndInd+1
    end
    y
end

function distancecurvemean(t::AbstractVector, P::AbstractVecOrMat{T}) where {T}
    # integrate distance

    integral = zero(T)

    # for each segment
    for i=1:length(t)-1
        # find primitive function

        # the 2nd degree polynomial for the segment is by construction nonnegative for all t (not only in the current segment)
        # use formula ∫√(x²+u²)dx = 1/2 x√(x²+u²) + 1/2u²ln(x + √(x²+u²))

        c,b,a = (P[:,i]...,)

        if abs(a)<1e-16
            # linear case
            error("Linear case not yet implemented, but should be very rare.")
        end

        b,c = b/a, c/a # x²+bx+c

        # complete the squares (x+α)²+β
        α = b/2
        β = c-α^2

        β<0 && error("Unexpected negative distance encountered (numerical problems?)")

        # translate function: now √(x²+β)
        # integrate from x0 to x1
        x0 = α
        x1 = α + t[i+1]-t[i]

        sq1 = sqrt(x1^2+β)
        Fx1 = x1*sq1 + β*log(x1+sq1)

        sq0 = sqrt(x0^2+β)
        Fx0 = x0*sq0 + β*log(x0+sq0)

        integral += sqrt(a) * 1/2*(Fx1-Fx0)
    end

    integral/(t[end]-t[1]) # mean
end


text:
# Matrix Product Operators
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Matrix%20Product%20Operators][Matrix Product Operators:1]]
"""
    MPO{L, T<:Number}

Matrix product operator on L sites. The `i`th tensor in the operator
has indices `[aⁱ⁻¹, aⁱ, σⁱ, σ′ⁱ]` where `(σⁱ, σ′ⁱ)` are the physical
indices and `(aⁱ⁻¹, aⁱ)` are bond indices.

A four site MPS would be diagrammatically represented

    σ¹          σ²          σ³          σ⁴
    |           |           |           |
    •--(a¹ a¹)--•--(a² a²)--•--(a³ a³)--•
    |           |           |           |
    σ′¹         σ′²         σ′³         σ′⁴

Note that `a⁰` and `aᴸ` must be of dimension 1.
"""
struct MPO{L, T<:Number}
    tensors::Vector{Array{T,4}}
end

"""
    MPO(W::Array{T,4}, L)
Create an `MPO` for `L` sites with all interior sites containing the tensor
`W`. The tensor is assumed to have the usual matrix-of-operators structure,
with the first two indices being the bond (matrix) dimension and the last two
indices being the physical (operator) dimension. The first and last sites only
use the last row and first column of `W`, respectively.

For example, the MPO form of the Hamiltonian for the TFIM is
constructed as with coupling `g` and length `L` is constructed as
follows:

    id = [1 0
          0 1]

    σᶻ = [1  0
          0 -1]

    σˣ = [0 1
          1 0]

    σʸ = [0  -im
          im   0]

    W = zeros(3, 3, 2, 2)
    W[1, 1, :, :] = id
    W[2, 1, :, :] = σᶻ
    W[3, 1, :, :] = -g*σˣ
    W[3, 2, :, :] = -σᶻ
    W[3, 3, :, :] = id

returning

    Ĥ::MPO = Ŵ¹ Ŵ² Ŵ³ ⋅⋅⋅ Ŵᴸ⁻¹ Wᴸ
"""
function MPO(W::Array{T,4}, L) where {T}
    L >= 2 || throw(DomainError(L, "At least 2 sites."))

    tensors = Vector{Array{T,4}}(undef, L)

    tensors[1] = W[end:end, :, :, :] # Row vector.
    for i in 2:(L-1)
        tensors[i] = W # Matrix
    end
    tensors[L] = W[:, 1:1, :, :] # Column vector.

    MPO{L,T}(tensors)
end

Base.:(==)(O::MPO, U::MPO) = O.tensors == U.tensors
Base.:(≈)(O::MPO, U::MPO)  = O.tensors ≈ U.tensors
Base.getindex(O::MPO, args...) = getindex(O.tensors, args...)

function Base.show(io::IO, ::MIME"text/plain", O::MPO{L, T}) where {L, T}
    d = length(O[2][1, 1, 1, :])
    bonddims = [size(O[i][:, :, 1, 1]) for i in 1:L]
    println(io, "Matrix product Operator on $L sites")
    _show_mpo_dims(io, L, d, bonddims)
end

function _show_mpo_dims(io::IO, L, d, bonddims)
    println(io, "  Physical dimension: $d")
    print(io, "  Bond dimensions:   ")
    if L > 8
        for i in 1:8
            print(io, bonddims[i], " × ")
        end
        print(io, " ... × ", bonddims[L])
    else
        for i in 1:(L-1)
            print(io, bonddims[i], " × ")
        end
        print(io, bonddims[L])
    end
end

function Base.show(io::IO, O::MPO{L, T}) where {L, T}
    print(io, "MPO on $L sites")
end
# Matrix Product Operators:1 ends here

ext: jl


text:
struct PeriodicHMM{F,T} <: AbstractHMM{F}
    a::Vector{T}
    A::Array{T,3}
    B::Matrix{<: Distribution{F}}
    PeriodicHMM{F,T}(a, A, B) where {F,T} = assert_periodichmm(a, A, B) && new(a, A, B)
end

PeriodicHMM(
    a::AbstractVector{T},
    A::AbstractArray{T,3},
    B::AbstractMatrix{<:Distribution{F}},
) where {F,T} = PeriodicHMM{F,T}(a, A, B)

PeriodicHMM(A::AbstractArray{T,3}, B::AbstractMatrix{<:Distribution{F}}) where {F,T} =
    PeriodicHMM{F,T}(ones(size(A, 1)) ./ size(A, 1), A, B)

function assert_periodichmm(a::AbstractVector, A::AbstractArray{T,3} where T, B::AbstractMatrix{<:Distribution})
    @argcheck isprobvec(a)
    @argcheck istransmats(A)
    @argcheck all(length.(B) .== length(B[1])) ArgumentError("All distributions must have the same dimensions")
    @argcheck size(A,3) == size(B,2) ArgumentError("Period length must be the same for transition matrix and distribution")
    @argcheck length(a) == size(A,1) == size(B,1) ArgumentError("Number of transition rates must match length of chain")
    return true
end

function assert_periodichmm(hmm::PeriodicHMM)
    assert_periodichmm(hmm.a, hmm.A, hmm.B)
end

istransmats(A::AbstractArray{T,3} where {T}) = all(i->istransmat(@view A[:,:,i]), 1:size(A,3))

# ==(h1::PeriodicHMM, h2::PeriodicHMM) = (h1.a == h2.a) && (h1.A == h2.A) && (h1.B == h2.B)

function rand(
    rng::AbstractRNG,
    hmm::PeriodicHMM,
    N::Integer;
    init = rand(rng, Categorical(hmm.a)),
    seq = false,
)
    T = size(hmm.B, 2)
    z = Vector{Int}(undef, N)
    (T >= 1) && (z[1] = init)
    n2t = CyclicArray(1:T, "1D")
    for n = 2:N
        tm1 = n2t[n-1] # periodic t-1
        z[n] = rand(rng, Categorical(hmm.A[z[n-1],:,tm1]))
    end
    y = rand(rng, hmm, z)
    seq ? (z, y) : y
end

function rand(rng::AbstractRNG, hmm::PeriodicHMM{Univariate}, z::AbstractVector{<:Integer})
    y = Vector{Float64}(undef, length(z))
    T = size(hmm.B, 2)
    n2t = CyclicArray(1:T, "1D")
    for n in eachindex(z)
        t = n2t[n] # periodic t
        y[n] = rand(rng, hmm.B[z[n],t])
    end
    y
end

function rand(
    rng::AbstractRNG,
    hmm::PeriodicHMM{Multivariate},
    z::AbstractVector{<:Integer},
)
    y = Matrix{Float64}(undef, length(z), size(hmm,2))
    T = size(hmm.B, 2)
    n2t = CyclicArray(1:T, "1D")
    for n in eachindex(z)
        t = n2t[n] # periodic t
        y[n, :] = rand(rng, hmm.B[z[n],t])
    end
    y
end

# should not be necessary with AbstractHMM type
# rand(hmm::PeriodicHMM, N::Integer; kwargs...) = rand(GLOBAL_RNG, hmm, N; kwargs...)
#
# rand(hmm::PeriodicHMM, z::AbstractVector{<:Integer}) = rand(GLOBAL_RNG, hmm, z)

"""
    size(hmm, [dim]) -> Int | Tuple
Return the number of states in `hmm`, the dimension of the observations and the length of the chain.
"""
size(hmm::PeriodicHMM, dim = :) = (size(hmm.B, 1), length(hmm.B[1]), size(hmm.B, 2))[dim]

copy(hmm::PeriodicHMM) = PeriodicHMM(copy(hmm.a), copy(hmm.A), copy(hmm.B))

function nparams(hmm::PeriodicHMM)
    (length(hmm.a) - 1) + (size(hmm.A,1) * size(hmm.A,2) - size(hmm.A,1)) * size(hmm.A,3)
end

###
function likelihoods!(L::AbstractMatrix, hmm::PeriodicHMM{Univariate}, observations)
    N, K, T = size(observations, 1), size(hmm, 1), size(hmm, 3)
    @argcheck size(LL) == (N, K)
    n2t = CyclicArray(1:T, "1D")
    @inbounds for i in OneTo(K), n in OneTo(T)
        t = n2t[n] # periodic t
        L[n, i] = pdf(hmm.B[i,t], observations[n])
    end
end

function likelihoods!(L::AbstractMatrix, hmm::PeriodicHMM{Multivariate}, observations)
    N, K, T = size(observations, 1), size(hmm, 1), size(hmm, 3)
    @argcheck size(LL) == (N, K)
    n2t = CyclicArray(1:T, "1D")
    @inbounds for i in OneTo(K), n in OneTo(T)
        t = n2t[n] # periodic t
        L[n, i] = pdf(hmm.B[i,t], view(observations, n, :))
    end
end

function loglikelihoods!(LL::AbstractMatrix, hmm::PeriodicHMM{Univariate}, observations)
    N, K, T = size(observations, 1), size(hmm, 1), size(hmm, 3)
    @argcheck size(LL) == (N, K)
    n2t = CyclicArray(1:T, "1D")
    @inbounds for i in OneTo(K), n in OneTo(N)
        t = n2t[n] # periodic t
        LL[n,i] = logpdf(hmm.B[i,t], observations[n])
    end
end

function loglikelihoods!(LL::AbstractMatrix, hmm::PeriodicHMM{Multivariate}, observations)
    N, K, T = size(observations, 1), size(hmm, 1), size(hmm, 3)
    @argcheck size(LL) == (N, K)
    n2t = CyclicArray(1:T, "1D")
    @inbounds for i in OneTo(K), n in OneTo(N)
        t = n2t[n] # periodic t
        LL[n,i] = logpdf(hmm.B[i,t], view(observations, n, :))
    end
end

####

function update_A!(
    A::AbstractArray{T,3} where {T},
    ξ::AbstractArray,
    α::AbstractMatrix,
    β::AbstractMatrix,
    LL::AbstractMatrix,
)
    @argcheck size(α, 1) == size(β, 1) == size(LL, 1) == size(ξ, 1)
    @argcheck size(α, 2) ==
              size(β, 2) ==
              size(LL, 2) ==
              size(A, 1) ==
              size(A, 2) ==
              size(ξ, 2) ==
              size(ξ, 3)

    N, K = size(LL)
    T = size(A,3)
    n2t = CyclicArray(1:T, "1D")
    @inbounds for n in OneTo(N - 1)
        t = n2t[n] # periodic t
        m = vec_maximum(view(LL, n + 1, :))
        c = 0.0

        for i in OneTo(K), j in OneTo(K)
            ξ[n, i, j] = α[n, i] * A[i, j, t] * exp(LL[n+1, j] - m) * β[n+1, j]
            c += ξ[n, i, j]
        end

        for i in OneTo(K), j in OneTo(K)
            ξ[n, i, j] /= c
        end
    end

    fill!(A, 0.0)
    ## For periodicHMM only the n observation corresponding to A(t) are used to update A(t)
    t_n = n2t[1:N]
    n_in_t = [findall(t_n .== t) for t in 1:T] # could probably be speeded up

    @inbounds for t in OneTo(T)
            for i in OneTo(K)
                c = 0.0

                for j in OneTo(K)
                    for n in setdiff(n_in_t[t],N)
                        A[i, j, t] += ξ[n, i, j]
                    end
                    c += A[i, j, t]
                end

                for j in OneTo(K)
                    A[i, j, t] /= c
                end
            end
    end
end

# # In-place update of the observations distributions.
function update_B!(B::AbstractMatrix, γ::AbstractMatrix, observations, estimator)
    @argcheck size(γ, 1) == size(observations, 1)
    @argcheck size(γ, 2) == size(B, 1)
    N = size(γ, 1)
    K = size(B, 1)
    T = size(B, 2)
    ## For periodicHMM only the n observation corresponding to A(t) are used to update A(t)
    n2t = CyclicArray(1:T, "1D")
    t_n = n2t[1:N]
    n_in_t = [findall(t_n .== t) for t in 1:T] # could probably be speeded up. For exemple computed outside the function only once

    @inbounds for t in OneTo(T)
        n_t = n_in_t[t]
        for i in OneTo(K)
            if sum(γ[n_t, i]) > 0
                B[i,t] = estimator(typeof(B[i,t]), permutedims(observations[n_t,:]), γ[n_t, i])
            end
        end
    end
end

function fit_mle!(
    hmm::PeriodicHMM,
    observations;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    estimator = fit_mle,
    fit_dists = false
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T = size(observations, 1), size(hmm, 1), size(hmm, 3)
    @argcheck T == size(hmm.B, 2)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates
    c = zeros(N)
    α = zeros(N, K)
    β = zeros(N, K)
    γ = zeros(N, K)
    ξ = zeros(N, K, K)
    LL = zeros(N, K)

    loglikelihoods!(LL, hmm, observations)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL)
    backwardlog!(β, c, hmm.a, hmm.A, LL)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, ξ, α, β, LL)
        fit_dists && update_B!(hmm.B, γ, observations, estimator)

        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmats(hmm.A)

        if fit_dists
            loglikelihoods!(LL, hmm, observations)
            robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
        end

        forwardlog!(α, c, hmm.a, hmm.A, LL)
        backwardlog!(β, c, hmm.a, hmm.A, LL)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history
end

# In-place forward pass, where α and c are allocated beforehand.
function forwardlog!(
    α::AbstractMatrix,
    c::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    LL::AbstractMatrix,
)
    @argcheck size(α, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(α, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    N, K = size(LL)
    T = size(A,3)

    fill!(α, 0.0)
    fill!(c, 0.0)

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        α[1, j] = a[j] * exp(LL[1, j] - m)
        c[1] += α[1, j]
    end

    for j in OneTo(K)
        α[1, j] /= c[1]
    end

    c[1] = log(c[1]) + m

    n2t = CyclicArray(1:T, "1D")
    @inbounds for n = 2:N
        tm1 = n2t[n-1] # periodic t-1
        m = vec_maximum(view(LL, n, :))

        for j in OneTo(K)
            for i in OneTo(K)
                α[n,j] += α[n-1, i] * A[i,j,tm1]
            end
            α[n, j] *= exp(LL[n, j] - m)
            c[n] += α[n, j]
        end

        for j in OneTo(K)
            α[n, j] /= c[n]
        end

        c[n] = log(c[n]) + m
    end
end

# In-place backward pass, where β and c are allocated beforehand.
function backwardlog!(
    β::AbstractMatrix,
    c::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    LL::AbstractMatrix,
)
    @argcheck size(β, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(β, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    N, K = size(LL)
    T = size(A,3)
    L = zeros(K)
    (T == 0) && return

    fill!(β, 0.0)
    fill!(c, 0.0)

    for j in OneTo(K)
        β[end, j] = 1.0
    end
    n2t = CyclicArray(1:T, "1D")

    @inbounds for n = N-1:-1:1
        t = n2t[n] # periodic t
        m = vec_maximum(view(LL, n + 1, :))

        for i in OneTo(K)
            L[i] = exp(LL[n+1, i] - m)
        end

        for j in OneTo(K)
            for i in OneTo(K)
                β[n,j] += β[n+1, i] * A[j,i,t] * L[i]
            end
            c[n+1] += β[n, j]
        end

        for j in OneTo(K)
            β[n, j] /= c[n+1]
        end

        c[n+1] = log(c[n+1]) + m
    end

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        c[1] += a[j] * exp(LL[1, j] - m) * β[1, j]
    end

    c[1] = log(c[1]) + m
end

function forward(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix)
    m = Matrix{Float64}(undef, size(LL))
    c = Vector{Float64}(undef, size(LL, 1))
    forwardlog!(m, c, a, A, LL)
    m, sum(c)
end

function backward(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix)
    m = Matrix{Float64}(undef, size(LL))
    c = Vector{Float64}(undef, size(LL, 1))
    backwardlog!(m, c, a, A, LL)
    m, sum(c)
end

function posteriors(a::AbstractVector, A::AbstractArray{T,3} where {T}, LL::AbstractMatrix; kwargs...)
    α, _ = forward(a, A, LL; kwargs...)
    β, _ = backward(a, A, LL; kwargs...)
    posteriors(α, β)
end

function viterbi!(
    T1::AbstractMatrix,
    T2::AbstractMatrix,
    z::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    L::AbstractMatrix,
)
    N, K = size(L)
    T = size(A,3)
    (T == 0) && return

    fill!(T1, 0.0)
    fill!(T2, 0)

    c = 0.0

    for i in OneTo(K)
        T1[1, i] = a[i] * L[1, i]
        c += T1[1, i]
    end

    for i in OneTo(K)
        T1[1, i] /= c
    end

    n2t = CyclicArray(1:T, "1D")
    @inbounds for n = 2:N
        tm1 = n2t[n-1] # t-1
        c = 0.0

        for j in OneTo(K)
            # TODO: If there is NaNs in T1 this may
            # stay to 0 (NaN > -Inf == false).
            # Hence it will crash when computing z[t].
            # Maybe we should check for NaNs beforehand ?
            amax = 0
            vmax = -Inf

            for i in OneTo(K)
                v = T1[n-1, i] * A[i, j, tm1]
                if v > vmax
                    amax = i
                    vmax = v
                end
            end

            T1[n, j] = vmax * L[n, j]
            T2[n, j] = amax
            c += T1[n, j]
        end

        for i in OneTo(K)
            T1[n, i] /= c
        end
    end

    z[N] = argmax(T1[N, :])
    for n = N-1:-1:1
        z[n] = T2[n+1, z[n+1]]
    end
end

function viterbilog!(
    T1::AbstractMatrix,
    T2::AbstractMatrix,
    z::AbstractVector,
    a::AbstractVector,
    A::AbstractArray{T,3} where {T},
    LL::AbstractMatrix,
)
    N, K = size(LL)
    T = size(A,3)
    (T == 0) && return

    fill!(T1, 0.0)
    fill!(T2, 0)

    al = log.(a)
    Al = log.(A)

    for i in OneTo(K)
        T1[1, i] = al[i] + LL[1, i]
    end
    n2t = CyclicArray(1:T, "1D")
    @inbounds for n = 2:N
        tm1 = n2t[n-1] # t-1
        for j in OneTo(K)
            amax = 0
            vmax = -Inf

            for i in OneTo(K)
                v = T1[n-1, i] + Al[i, j, tm1]
                if v > vmax
                    amax = i
                    vmax = v
                end
            end

            T1[n, j] = vmax + LL[n, j]
            T2[n, j] = amax
        end
    end

    z[N] = argmax(T1[N, :])
    for n = N-1:-1:1
        z[n] = T2[n+1, z[n+1]]
    end
end


text:
# Autogenerated wrapper script for libgraphicsmagic_jll for x86_64-apple-darwin
export gm, libGraphicsMagic, libGraphicsMagicPlusPlus, libGraphicsMagicWand

using CompilerSupportLibraries_jll
using libwebp_jll
using libpng_jll
using FreeType2_jll
using Ghostscript_jll
using Xorg_libXext_jll
using Xorg_libSM_jll
using Libtiff_jll
using LittleCMS_jll
using Bzip2_jll
using XZ_jll
using Libtool_jll
JLLWrappers.@generate_wrapper_header("libgraphicsmagic")
JLLWrappers.@declare_executable_product(gm)
JLLWrappers.@declare_library_product(libGraphicsMagic, "@rpath/libGraphicsMagick.3.dylib")
JLLWrappers.@declare_library_product(libGraphicsMagicPlusPlus, "@rpath/libGraphicsMagick++.12.dylib")
JLLWrappers.@declare_library_product(libGraphicsMagicWand, "@rpath/libGraphicsMagickWand.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, libwebp_jll, libpng_jll, FreeType2_jll, Ghostscript_jll, Xorg_libXext_jll, Xorg_libSM_jll, Libtiff_jll, LittleCMS_jll, Bzip2_jll, XZ_jll, Libtool_jll)
    JLLWrappers.@init_executable_product(
        gm,
        "bin/gm",
    )

    JLLWrappers.@init_library_product(
        libGraphicsMagic,
        "lib/libGraphicsMagick.3.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGraphicsMagicPlusPlus,
        "lib/libGraphicsMagick++.12.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libGraphicsMagicWand,
        "lib/libGraphicsMagickWand.2.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# a 2D triangulation example
# this example uses a factor graph to represent the data
# Factor graphs can be used to construct either fixed lag smoothing or full blown SLAM solutions.

# Load required libraries
using Caesar, RoMEPlotting
# using Distributions
# using ApproxManifoldProducts

## Environment setup -- in Juno you can run cells with Alt+Enter

# true pose position (assume hidden)  - [X Y θ]
p0 = [0.0;0.0;0.0]

# noise level for Gaussian assumption -- see https://github.com/JuliaStats/Distributions.jl/issues/584
pcov = Matrix(Diagonal([0.01; 0.01; 0.05].^2))

# true landmark positions
l0 = [2.0;0.15]
l1 = [2.0;-0.15]

d0 = l0-p0[1:2]
m0 = atan(d0[2],d0[1])

d1 = l1-p0[1:2]
m1 = atan(d1[2],d1[1])

mstd = 0.1

## build the factor graph

fg = initfg()

# add two beacons
addVariable!(fg, :l0, Point2, tags=[:BEACON;])
# addFactor!(fg, [:l0], Prior(MvNormal(l0, lcov)))

addVariable!(fg, :l1, Point2, tags=[:BEACON;])
# addFactor!(fg, [:l1], Prior(MvNormal(l1, lcov)))

addFactor!(fg, [:l0; :l1], Point2Point2(MvNormal(l1-l0, [0.01 0; 0 0.01].^2)))

# add unknown pose location
addVariable!(fg, :x0, Pose2, tags=[:POSE;])
addFactor!(fg, [:x0], PriorPose2(MvNormal(p0, pcov)))

# add two bearing only measurements
addFactor!(fg, [:x0; :l0], Pose2Point2BearingRange(Normal(m0, mstd), Rayleigh(2.0)))
addFactor!(fg, [:x0; :l1], Pose2Point2BearingRange(Normal(m1, mstd), Rayleigh(2.0)))

# addFactor!(fg, [:l1], PartialPrior(Normal(0.0,0.1), (1,)))

# # add unknown pose location
# addVariable!(fg, :x1, Pose2, tags=[:POSE;])
# addFactor!(fg, [:x1], PriorPose2(MvNormal(p1, pcov)))
#
# # add two bearing only measurements
# addFactor!(fg, [:x1; :l0], Pose2Point2Bearing(Normal(m2, mstd)))
# addFactor!(fg, [:x1; :l1], Pose2Point2Bearing(Normal(m3, mstd)))
#
#
# addFactor!(fg, [:x0; :x1], Pose2Pose2(MvNormal([0.0;1.0;0.0], [0.01 0 0; 0 0.01 0; 0 0 0.01].^2)))

drawGraph(fg, show=true)

# ensureAllInitialized!(fg)
# setValKDE!(fg, :l0, AMP.manikde!(randn(2,100), Point2().manifolds))
# setValKDE!(fg, :l1, AMP.manikde!(randn(2,100), Point2().manifolds))

# solve the system -- remember first runs are slow (just in time compile code)
solveTree!(fg)

## Plot values to see result

Gadfly.set_default_plot_size(30cm,20cm)
plotPose(Pose2(), [getKDE(fg, :x0)]) #, axis=[-0.5 2.0; -0.5 1.5])

# plotPose(Pose2(), [getKDE(fg, :x1)]) #, axis=[-0.5 2.0; -0.5 1.5])

plotKDE(fg, [:l0; :l1],dims=[1;2],levels=3) # ,axis=[0 2.0; 0.0 1.5])

## Building a test

fg = initfg()

addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3),Matrix(Diagonal(0.001.*ones(3))))))

addVariable!(fg, :l0, Point2)
addFactor!(fg, [:l0], PriorPoint2(MvNormal([1.0;0.0],Matrix(Diagonal(0.001.*ones(2))))))
addFactor!(fg, [:x0;:l0], Pose2Point2Bearing(Normal(0, 0.01)))

addVariable!(fg, :l1, Point2)
addFactor!(fg, [:l1], PriorPoint2(MvNormal([0.0;1.0],Matrix(Diagonal(0.001.*ones(2))))))
addFactor!(fg, [:x0;:l1], Pose2Point2Bearing(Normal(pi/2, 0.01)))

fg = initfg()

addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3),Matrix(Diagonal(0.001.*ones(3))))))

addVariable!(fg, :l0, Point2)
# addFactor!(fg, [:l0], PriorPoint2(MvNormal([1.0;0.0],Matrix(Diagonal(0.001.*ones(2))))))
addFactor!(fg, [:x0;:l0], Pose2Point2Bearing(Normal(0, 0.01)))

addVariable!(fg, :l1, Point2)
# addFactor!(fg, [:l1], PriorPoint(MvNormal([0.0;1.0],Matrix(Diagonal(0.001.*ones(2))))))
addFactor!(fg, [:x0;:l1], Pose2Point2Bearing(Normal(pi/2, 0.01)))

addFactor!(fg, [:l0;:l1], Point2Point2(MvNormal([-2.0; 2.0], Matrix(Diagonal(0.001.*ones(2))))))

## util development

## more debugging

stuff = localProduct(fg, :x0)

plotKDE(stuff[2],dims=[1;2],levels=3, c=["red","green","blue"])

stuff = localProduct(fg, :l0)

plotKDE([stuff[1];stuff[2]],dims=[1;2],levels=3, c=["black","red","green","blue"])

pcm = [marginal(stuff[2][1],[3]);marginal(stuff[2][2],[3])]
AMP.plotKDECircular(pcm)

np = manifoldProduct(stuff[2], Pose2().manifolds)

AMP.plotKDECircular(marginal(np,[3]))

## more direct calculation

pts0 = approxConv(fg, :x0l0f1, :x0)
X0 = manikde!(pts0, Pose2().manifolds)

plotPose(Pose2(), [X0])

pts1 = approxConv(fg, :x0l1f1, :x0)
X1 = manikde!(pts1, Pose2().manifolds)

plotPose(Pose2(), [X1])

## Major error here

ptsl0 = approxConv(fg, :x0l0f1, :l0)
L0 = manikde!(ptsl0, Point2().manifolds)
plotKDE(L0, levels=3)

stuff = IIF.localProduct(fg, :l0)

plotKDE([stuff[2];stuff[1]])

## more debugging below

pts0_0 = approxConv(fg, :x0l0f1, :l0)
L00 = manikde!(pts0_0, Point2().manifolds)

plotKDE(L00, levels=3, axis=[-0.5 2; -0.5 1.5])

pts0_1 = approxConv(fg, :x0l1f1, :l1)
L01 = manikde!(pts0_1, Point2().manifolds)

plotKDE(L01, levels=3, axis=[-0.5 2; -0.5 1.5])

#

getSolverParams(fg).drawtree=true
getSolverParams(fg).showtree=true
tree = solveTree!(fg)

spyCliqMat(tree,:l0)

writeGraphPdf(fg, show=true)

stuff = IIF.treeProductUp(fg, tree, :l0, :l0)

Gadfly.plot(x=stuff[1][1,:], y=stuff[1][2,:], Geom.histogram2d)

## Util development

##


text:
const score_table = (1 => 10, 5 => 5, 10 => 1)

"""
    score(x, y)

Return the score earned by a dart landing at coordinates (x, y).
"""
function score(x, y)
    z = hypot(x, y)
    for (radius, points) in score_table
        z <= radius && return points
    end
    return 0
end


text:
observation_variance = 1
transition_variance = 10

function lnprob(theta_i)
	lnprob_theta_i(self, theta_i, theta_t_minus_1, time_series,t,iter_):
        #ln poisson observations
            lnprob_theta_i = 1.0/(np.sqrt(2*np.pi*observation_variance))*np.exp(-.5*(1.0/observation_variance)*((time_series[t] - theta_i )**2))

            transition_sum = 0
	    for theta_t_minus_1_i in theta_t_minus_1:

                transition_sum += 1.0/(np.sqrt(2*np.pi*transition_variance))*np.exp(-.5*(1.0/transition_variance)*((theta_i - theta_t_minus_1_i )**2))

            return (np.log(lnprob_theta_i)+np.log(transition_sum))

time_series = [1,2,3,4]

max_stars_repo_name: gcgibson/ssvgd

max_issues_repo_name: gcgibson/ssvgd

max_forks_repo_name: gcgibson/ssvgd


text:
An example of using a context menu with cairo graphics.

Cairo graphics only work with `Tk` (not `Qt`) so we specify that, then load our packages:

```
ENV["toolkit"] = "Tk"
using JGUI, Winston
```

Our GUI is basic, consisting of a top-level window and a graphic window.

```
w = window()
c = cairographic(w)
push!(w, c)
```

We can plot to the `cairographic` window via `Winston.display`:

```
f(x) = sin(x)
p = plot(f, 0, pi)
Winston.display(c[:widget], p)
```

Adding a context menu is easy. We just pass a widget to the `menu`
constructor instead of a top-level window (which is how we add a menu
bar).

```
m = menu(c)
```

The menu actions are added as with any menu.
The context value is the pixel coordinates of the mouse click (this needs work!):

```
a1 = action(m, "Position", () -> println(c[:context]))
addAction(m, a1)
```

Radio groups can be easily added to a menu bar. Just pass the context menu to the constructor as the parent. Checkboxes can also be added this way.

```
rb = radiogroup(m, ["one", "two"])
connect(rb, "valueChanged") do value
    value == "one" && println("one")
    value == "two" && println("two")
end

addAction(m, separator(m))
addAction(m, rb)
```


text:
import Base: match
using Combinatorics: combinations, permutations

export match

mutable struct Match <: AbstractSet{AbstractDict{Variable,Any}}
    matches::Set{Dict{Variable,Any}}
end
Match(xs::Union{Pair,Dict}...) = Match(Set(Dict.(xs)))
Base.zero(::Type{Match}) = Match()
Base.one(::Type{Match}) = Match(Dict())
Base.length(Θ::Match) = length(Θ.matches)
Base.iterate(Θ::Match) = iterate(Θ.matches)
Base.iterate(Θ::Match, state) = iterate(Θ.matches, state)
Base.push!(Θ::Match, items...) = (push!(Θ.matches, items...); Θ)
Base.copy(Θ::Match) = Match(copy(Θ.matches))
Base.union(Θ₁::Match, Θ₂::Match) = Match(union(Θ₁.matches, Θ₂.matches))
(Θ::Match)(t::Term) = Set{Term}([replace(t, σ) for σ ∈ Θ])

function Base.merge!(Θ::Match, Θs::Match...)
    for Θ′ ∈ Θs
        result = Match()
        for σ′ ∈ Θ′
            foreach(Θ) do σ
                res = copy(σ)
                for (k, v) in σ′
                    if haskey(σ, k)
                        σ[k] == v || return
                    end
                    res[k] = v
                end
                push!(result, res)
            end
        end
        Θ.matches = result
    end
    Θ
end
Base.merge(σ::Match, σs::Match...) = merge!(one(Match), σ, σs...)

"""
    match(pattern::Term, subject::Term) -> Match

Match term `t` to `pattern`, producing a `Match` if the process succeeds.

# Examples
```jldoctest
julia> @vars x;
       @syms a b;

julia> match(@term(x), @term(f(a)))
Match(Set(Dict{Variable,Any}[Dict(x=>:((f)(a)))]))

julia> match(@term(f(x)), @term(a))
Match(Set(Dict{Variable,Any}[]))

julia> match(@term(f(x, x)), @term(f(a, b)))
Match(Set(Dict{Variable,Any}[]))

julia> match(@term(f(x, x)), @term(f(a, a)))
Match(Set(Dict{Variable,Any}[Dict(x=>a)]))
```
"""
match(pattern::Term, subject::Term) = match(Term, get(pattern), get(subject), one(Match))

match(::Type{Term}, x::Variable, t, Θ) = merge(Θ, Match(x => t))
function match(::Type{Term}, p::Expr, s::Expr, Θ)
    p.head === s.head   || return zero(Match)

    if p.head === :call
        f = s.args[1]
        isvalid(Commutative(f)) && return match(Commutative, p, s, Θ)
        isvalid(Associative(f)) && return match(Associative, p, s, Θ)
    end

    _match(p, s, Θ)
end
match(::Type{Term}, p, s, Θ) = (typeof(s) <: typeof(p) && p == s) ? Θ : zero(Match)

function _match(f::Expr, g::Expr, Θ)
    length(f.args) == length(g.args) || return zero(Match)

    for (x, y) ∈ zip(f.args, g.args)
        Θ = match(Term, x, y, Θ)
    end

    Θ
end

"""
    match(::Type{Associative}, p::Expr, s::Expr, Θ::Match) -> Match

Match an associative function call to another associative function call, based on the
algorithm by [Krebber](https://arxiv.org/abs/1705.00907). Requires `p` and `s` to have
head `:call`.
"""
function match(::Type{Associative}, p::Expr, s::Expr, Θ)
    @assert p.head === s.head === :call

    pname, pargs = p.args[1], p.args[2:end]
    sname, sargs = s.args[1], s.args[2:end]
    Θ = match(Term, pname, sname, Θ)

    m, n = length(pargs), length(sargs)
    m > n && return zero(Match)
    n_free = n - m
    n_vars = count(x -> x isa Variable, pargs)
    Θᵣ = zero(Match)

    for k ∈ Iterators.product((0:n_free for i ∈ 1:n_vars)...)
        (isempty(k) ? 0 : sum(k)) == n_free || continue  # FIXME: efficiency
        i, j = 1, 1
        Θ′ = Θ
        for pₗ ∈ pargs
            l_sub = 0
            if pₗ isa Variable
                l_sub += k[j]
                j += 1
            end
            s′ = l_sub > 0 ? Expr(:call, sname, sargs[i:i+l_sub]...) : sargs[i]
            Θ′ = match(Term, pₗ, s′, Θ′)
            isempty(Θ′) && break
            i += l_sub + 1
        end
        Θᵣ = Θᵣ ∪ Θ′
    end
    Θᵣ
end
"""
    match(::Type{Commutative}, p::Expr, s::Expr, Θ::Match) -> Match

Match a commutative function call to another commutative function call.
Requires `p` and `s` to have head `:call`.
"""
function match(::Type{Commutative}, p::Expr, s::Expr, Θ)
    @assert p.head === s.head === :call

    matches = map(perms(s)) do s′  # FIXME: efficiency
        isvalid(Associative(s.args[1])) && return match(Associative, p, s′, Θ)
        _match(p, s′, Θ)
    end
    reduce(union, matches)
end
function perms(ex::Expr)
    @assert ex.head === :call

    f, args = ex.args[1], ex.args[2:end]

    map(permutations(args)) do perm
        Expr(:call, f, perm...)
    end |> unique
end


text:
# synthetic_ECG_include.jl

include("../tcm_code/tcm_estimateForest_L1.jl")
include("../tcm_code/tcm_estimateForest_L1_MIO.jl")
include("../tcm_code/tcm_randomCompleteTree.jl")
include("../tcm_code/tcm_estimateForest_L1_RTS.jl")
include("../tcm_code/tcm_predictForest.jl")

include("synthetic_ECG_helperfn.jl")

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
module HDDriverModule

using Profile
using HierarchicalStateMachineModule: HierarchicalStateMachine
using HDModule: RiverParameters,RiverPrognosticFields,RiverPrognosticFieldsOnly,RunHD,handle_event
using HDModule: SetDrainage,SetRunoff,PrintResults,PrognosticFields,WriteRiverInitialValues
using HDModule: WriteRiverFlow,AccumulateRiverFlow,ResetCumulativeRiverFlow,WriteMeanRiverFlow
using HDModule: PrintGlobalValues, water_to_lakes,water_from_lakes,get_river_parameters
using HDModule: SetLakeEvaporation
using FieldModule: Field
using LakeModule: LakeParameters,LakePrognostics,LakeFields,RiverAndLakePrognosticFields,RunLakes
using LakeModule: PrintSection,WriteLakeNumbers,WriteLakeVolumes,SetupLakes,DistributeSpillover
using LakeModule: WriteDiagnosticLakeVolumes,CheckWaterBudget
using LakeModule: water_to_lakes,handle_event,water_from_lakes
using GridModule: get_number_of_cells

FloatFieldOrNothing = Union{Field{Float64},Nothing}

function drive_hd_model_with_or_without_lakes(prognostic_fields::PrognosticFields,
                                              drainages::Array{Field{Float64},1},
                                              runoffs::Array{Field{Float64},1},
                                              lake_evaporations::Array{Field{Float64},1},
                                              timesteps::Int64,run_lakes_flag::Bool,
                                              process_initial_lake_water::Bool=false,
                                              initial_water_to_lake_centers::FloatFieldOrNothing=
                                              nothing,
                                              initial_spillover_to_rivers::FloatFieldOrNothing=
                                              nothing;
                                              print_timestep_results::Bool=true,
                                              output_timestep::Int64=100,
                                              write_output::Bool=true)
  ncells::Int64 = get_number_of_cells(get_river_parameters(prognostic_fields).grid)
  hsm::HierarchicalStateMachine = HierarchicalStateMachine(prognostic_fields)
  runHD::RunHD = RunHD()
  if run_lakes_flag
    run_lakes::RunLakes = RunLakes()
    if process_initial_lake_water
      setup_lakes::SetupLakes = SetupLakes(initial_water_to_lake_centers)
      distribute_spillover::DistributeSpillover =
        DistributeSpillover(initial_spillover_to_rivers)
      handle_event(hsm,setup_lakes)
      # Run lakes once to cascade water from any overflowing lakes
      # and process through flow if not using instant through flow
      handle_event(hsm,run_lakes)
      # This must go after the initial lake run as the lake run
      # will reset the water to hd.
      handle_event(hsm,distribute_spillover)
      handle_event(hsm,CheckWaterBudget(sum(initial_water_to_lake_centers)+
                                        sum(initial_spillover_to_rivers)))
    end
  end
  for i in 1:timesteps
    if size(drainages) == (1,)
      set_drainage = SetDrainage(deepcopy(drainages[1]))
    else
      set_drainage = SetDrainage(deepcopy(drainages[convert(Int64,ceil(convert(Float64,i)/30.0))]))
    end
    handle_event(hsm,set_drainage)
    if size(runoffs) == (1,)
      set_runoff = SetRunoff(deepcopy(runoffs[1]))
    else
      set_runoff = SetRunoff(deepcopy(runoffs[convert(Int64,ceil(convert(Float64,i)/30.0))]))
    end
    handle_event(hsm,set_runoff)
    if size(lake_evaporations) == (1,)
      set_lake_evaporation =
        SetLakeEvaporation(deepcopy(lake_evaporations[1]))
    else
      set_lake_evaporation =
        SetLakeEvaporation(deepcopy(lake_evaporations[convert(Int64,ceil(convert(Float64,i)/30.0))]))
    end
    handle_event(hsm,set_lake_evaporation)
    handle_event(hsm,runHD)
    if run_lakes_flag
      handle_event(hsm,run_lakes)
      handle_event(hsm,CheckWaterBudget())
    end
    if print_timestep_results
      if i%output_timestep == 0 || i == 1
        if ncells < 1000
          print_results::PrintResults = PrintResults(i)
          handle_event(hsm,print_results)
        else
          print_section::PrintSection = PrintSection()
          handle_event(hsm,print_section)
        end
      end
    end
    if write_output
      accumulate_river_flow::AccumulateRiverFlow = AccumulateRiverFlow()
      handle_event(hsm,accumulate_river_flow)
      if i%output_timestep == 0 || i == 1
        print_global_values::PrintGlobalValues = PrintGlobalValues()
        handle_event(hsm,print_global_values::PrintGlobalValues)
        #write_mean_river_flow::WriteMeanRiverFlow = WriteMeanRiverFlow(i,output_timestep)
        #handle_event(hsm,write_mean_river_flow)
        write_river_flow::WriteRiverFlow = WriteRiverFlow(i)
        handle_event(hsm,write_river_flow)
        reset_cumulative_river_flow::ResetCumulativeRiverFlow = ResetCumulativeRiverFlow()
        handle_event(hsm,reset_cumulative_river_flow)
        if run_lakes_flag
           #write_lake_numbers::WriteLakeNumbers = WriteLakeNumbers(i)
           #handle_event(hsm,write_lake_numbers)
           write_diagnostic_lake_volumes::WriteDiagnosticLakeVolumes =
              WriteDiagnosticLakeVolumes(i)
           handle_event(hsm,write_diagnostic_lake_volumes)
         end
      end
    end
  end
  if write_output
    write_river_initial_values::WriteRiverInitialValues = WriteRiverInitialValues()
    handle_event(hsm,write_river_initial_values)
    if run_lakes_flag
      write_lake_volumes::WriteLakeVolumes = WriteLakeVolumes()
      handle_event(hsm,write_lake_volumes)
    end
  end
end

function drive_hd_model(river_parameters::RiverParameters,river_fields::RiverPrognosticFields,
                        drainages::Array{Field{Float64},1},runoffs::Array{Field{Float64},1},
                        lake_evaporations::Array{Field{Float64},1},
                        timesteps::Int64;print_timestep_results::Bool=false,
                        output_timestep::Int64=100,
                        write_output::Bool=true,
                        return_output::Bool=false)
  prognostic_fields::PrognosticFields = RiverPrognosticFieldsOnly(river_parameters,river_fields)
  drive_hd_model_with_or_without_lakes(prognostic_fields,drainages,runoffs,
                                       lake_evaporations,timesteps,false,
                                       print_timestep_results=print_timestep_results,
                                       output_timestep=output_timestep,
                                       write_output=write_output)
  if return_output
    return river_fields.water_to_ocean,river_fields.river_inflow
  end
end

function drive_hd_model(river_parameters::RiverParameters,drainages::Array{Field{Float64},1},
                        runoffs::Array{Field{Float64},1},
                        lake_evaporations::Array{Field{Float64},1},timesteps::Int64;
                        print_timestep_results::Bool=false,
                        output_timestep::Int64=100,
                        write_output::Bool=true,
                        return_output::Bool=false)
  river_fields::RiverPrognosticFields = RiverPrognosticFields(river_parameters)
  drive_hd_model(river_parameters,river_fields,drainages,runoffs,lake_evaporations,
                 timesteps,print_timestep_results=print_timestep_results,
                 output_timestep=output_timestep,
                 write_output=write_output,
                 return_output=false)
  if return_output
    return river_fields.water_to_ocean,river_fields.river_inflow
  end
end

function drive_hd_and_lake_model(river_parameters::RiverParameters,river_fields::RiverPrognosticFields,
                                 lake_parameters::LakeParameters,lake_prognostics::LakePrognostics,
                                 lake_fields::LakeFields,drainages::Array{Field{Float64},1},
                                 runoffs::Array{Field{Float64},1},
                                 lake_evaporations::Array{Field{Float64},1},timesteps::Int64,
                                 process_initial_lake_water::Bool=false,
                                 initial_water_to_lake_centers::FloatFieldOrNothing=nothing,
                                 initial_spillover_to_rivers::FloatFieldOrNothing=nothing;
                                 print_timestep_results::Bool=false,
                                 output_timestep::Int64=100,
                                 write_output::Bool=true,
                                 return_output::Bool=true)
  prognostic_fields::PrognosticFields = RiverAndLakePrognosticFields(river_parameters,river_fields,
                                                                     lake_parameters,lake_prognostics,
                                                                     lake_fields)
  drive_hd_model_with_or_without_lakes(prognostic_fields,drainages,runoffs,
                                       lake_evaporations,timesteps,
                                       true,process_initial_lake_water,
                                       initial_water_to_lake_centers,
                                       initial_spillover_to_rivers,
                                       print_timestep_results=print_timestep_results,
                                       output_timestep=output_timestep,
                                       write_output=write_output)
  if return_output
    return river_fields,lake_prognostics,lake_fields
  end
end

function drive_hd_and_lake_model(river_parameters::RiverParameters,
                                 river_fields::RiverPrognosticFields,
                                 lake_parameters::LakeParameters,
                                 drainages::Array{Field{Float64},1},
                                 runoffs::Array{Field{Float64},1},
                                 lake_evaporations::Array{Field{Float64},1},
                                 timesteps::Int64,
                                 process_initial_lake_water::Bool=false,
                                 initial_water_to_lake_centers::FloatFieldOrNothing=nothing,
                                 initial_spillover_to_rivers::FloatFieldOrNothing=nothing;
                                 print_timestep_results::Bool=false,
                                 output_timestep::Int64=100,
                                 write_output::Bool=true,
                                 return_output::Bool=false)
  lake_fields::LakeFields = LakeFields(lake_parameters)
  lake_prognostics::LakePrognostics = LakePrognostics(lake_parameters,lake_fields)
  drive_hd_and_lake_model(river_parameters,river_fields,lake_parameters,lake_prognostics,
                          lake_fields,drainages,runoffs,lake_evaporations,timesteps,
                          process_initial_lake_water,initial_water_to_lake_centers,
                          initial_spillover_to_rivers,
                          print_timestep_results=print_timestep_results,
                          output_timestep=output_timestep,
                          write_output=write_output,
                          return_output=false)
  if return_output
    return river_fields,lake_prognostics,lake_fields
  end
end

function drive_hd_and_lake_model(river_parameters::RiverParameters,lake_parameters::LakeParameters,
                                 lake_prognostics::LakePrognostics,lake_fields::LakeFields,
                                 drainages::Array{Field{Float64},1},
                                 runoffs::Array{Field{Float64},1},
                                 lake_evaporations::Array{Field{Float64},1},timesteps::Int64,
                                 process_initial_lake_water::Bool=false,
                                 initial_water_to_lake_centers::FloatFieldOrNothing=nothing,
                                 initial_spillover_to_rivers::FloatFieldOrNothing=nothing;
                                 print_timestep_results::Bool=false,
                                 output_timestep::Int64=100,
                                 write_output::Bool=true,
                                 return_output::Bool=false)
  river_fields::RiverPrognosticFields = RiverPrognosticFields(river_parameters)
  drive_hd_and_lake_model(river_parameters,river_fields,lake_parameters,lake_prognostics,
                          lake_fields,drainages,runoffs,lake_evaporations,timesteps,
                          process_initial_lake_water,initial_water_to_lake_centers,
                          initial_spillover_to_rivers,
                          print_timestep_results=print_timestep_results,
                          output_timestep=output_timestep,
                          write_output=write_output,
                          return_output=false)
  if return_output
    return river_fields,lake_prognostics,lake_fields
  end
end

function drive_hd_and_lake_model(river_parameters::RiverParameters,lake_parameters::LakeParameters,
                                 drainages::Array{Field{Float64},1},runoffs::Array{Field{Float64},1},
                                 lake_evaporations::Array{Field{Float64},1},timesteps::Int64,
                                 process_initial_lake_water::Bool=false,
                                 initial_water_to_lake_centers::FloatFieldOrNothing=nothing,
                                 initial_spillover_to_rivers::FloatFieldOrNothing=nothing;
                                 print_timestep_results::Bool=false,
                                 output_timestep::Int64=100,
                                 write_output::Bool=true,
                                 return_output::Bool=false)
  river_fields::RiverPrognosticFields = RiverPrognosticFields(river_parameters)
  lake_fields::LakeFields = LakeFields(lake_parameters)
  lake_prognostics::LakePrognostics = LakePrognostics(lake_parameters,lake_fields)
  drive_hd_and_lake_model(river_parameters,river_fields,lake_parameters,lake_prognostics,
                          lake_fields,drainages,runoffs,lake_evaporations,timesteps,
                          process_initial_lake_water,initial_water_to_lake_centers,
                          initial_spillover_to_rivers,
                          print_timestep_results=print_timestep_results,
                          output_timestep=output_timestep,
                          write_output=write_output,
                          return_output=false)
  if return_output
    return river_fields,lake_prognostics,lake_fields
  end
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
"""
    SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer

The `SingleBridgeOptimizer` bridges any constraint supported by the bridge `BT`.
This is in contrast with the [`LazyBridgeOptimizer`](@ref) which only bridges the constraints that are unsupported by the internal model, even if they are supported by one of its bridges.
"""
mutable struct SingleBridgeOptimizer{BT<:AbstractBridge, OT<:MOI.ModelLike} <: AbstractBridgeOptimizer
    model::OT
    con_to_name::Dict{CI, String}
    name_to_con::Union{Dict{String, MOI.ConstraintIndex}, Nothing}
    # Constraint Index of bridged constraint -> Bridge.
    # It is set to `nothing` when the constraint is deleted.
    bridges::Vector{Union{Nothing, AbstractBridge}}
    # Constraint Index of bridged constraint -> Constraint type.
    constraint_types::Vector{Tuple{DataType, DataType}}
    # For `SingleVariable` constraints: (variable, set type) -> bridge
    single_variable_constraints::Dict{Tuple{Int64, DataType}, AbstractBridge}
end
function SingleBridgeOptimizer{BT}(model::OT) where {BT, OT <: MOI.ModelLike}
    SingleBridgeOptimizer{BT, OT}(
        model, Dict{CI, String}(), nothing,
        Union{Nothing, AbstractBridge}[], Tuple{DataType, DataType}[],
        Dict{Tuple{Int64, DataType}, AbstractBridge}())
end

function supports_bridging_constraint(b::SingleBridgeOptimizer{BT},
                                      F::Type{<:MOI.AbstractFunction},
                                      S::Type{<:MOI.AbstractSet}) where BT
    return MOI.supports_constraint(BT, F, S)
end
function is_bridged(b::SingleBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
                    S::Type{<:MOI.AbstractSet})
    return supports_bridging_constraint(b, F, S)
end
bridge_type(b::SingleBridgeOptimizer{BT}, ::Type{<:MOI.AbstractFunction}, ::Type{<:MOI.AbstractSet}) where BT = BT


text:
# Copyright (c) 2014: Daniel C. Jones
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

# the following tests are copied from Switch.jl.

# basic switch constructs
function switch_with_fallthrough(x::Int)
    xs = Int[]
    @cswitch x begin
        @case 1
            push!(xs, 1)

        @case 2
            push!(xs, 2)

        @case 3
            push!(xs, 3)

        @default
            push!(xs, 4)
    end
    return xs
end

# fallthrough with weird default case placement
function switch_with_leading_default(x::Int)
    xs = Int[]
    @cswitch x begin
        @default
            push!(xs, 4)

        @case 1
            push!(xs, 1)

        @case 2
            push!(xs, 2)

        @case 3
            push!(xs, 3)
    end
    return xs
end

# switching with expressions
function switch_with_expressions(x::Int)
    xs = Int[]
    @cswitch x + 1 begin
        @case 1 + 1
            push!(xs, 1)

        @case 2 + 1
            push!(xs, 2)

        @case 3 + 1
            push!(xs, 3)

        @default
            push!(xs, 4)
    end
    return xs
end

# switching with breaks
function switch_with_break(x::Int)
    xs = Int[]
    @cswitch x begin
        @case 1
            push!(xs, 1)
            break

        @case 2
            push!(xs, 2)
            break

        @case 3
            push!(xs, 3)
            break

        @default
            push!(xs, 4)
    end
    return xs
end


text:
mutable struct MyStruct a end
x = MyStruct(2)
x.a = 3 + true

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct Part
    base::Vector{Point}
    color::Color
    entrances::Vector{Point}
    exits::Vector{Point}
    junctions::Vector{Point}
    shape::Int
end

part(color::Color, shape::Int) = part(color, Val(shape))

base(p) = p.base
color(p) = p.color
entrances(p) = p.entrances
exits(p) = p.exits
junctions(p) = p.junctions
shape(p) = p.shape

is_base(p) = !isempty(base(p))
is_junction(p) = !isempty(junctions(p))
is_root(p) = isempty(entrances(p)) && !isempty(exits(p))
is_leaf(p) = !isempty(entrances(p)) && isempty(exits(p))

# Part 1
function part(color, ::Val{1})
    entrances = [Point(0, 0, 1)]
    exits = [Point(0, 0, 0)]
    return Part(empty_points, color, entrances, exits, empty_points, 1)
end

# Part 2
function part(color, ::Val{2})
    base = [Point(0, 0, 0)]
    junctions = [Point(0, 0, 1)]
    return Part(base, color, empty_points, empty_points, junctions, 2)
end


text:
include("BPX1D.jl")

function Schwarz_matrix_2D(A::SparseMatrixCSC{Float64, Int64}, L::Integer, BCs::AbstractArray{String, 1})
    B = spdiagm(0 => diag(A).^-1)
    for l=1:(L-1)
        ψ_x, χ_x = BPX_weights(L, l, BCs[1:2])
        P_x = BPX_layer_half_matrix(L, ψ_x, χ_x, BCs[1:2])
        ψ_y, χ_y = BPX_weights(L, l, BCs[3:4])
        P_y = BPX_layer_half_matrix(L, ψ_y, χ_y, BCs[3:4])
        P = kron(P_y, P_x)
        if l == 1
            A_P = sparse(inv(Array(transpose(P)*A*P + 1e-6*I)))
        else
            A_P = spdiagm(0 => diag(transpose(P)*A*P).^-1)
        end
        B += P*A_P*transpose(P)
    end
    return B
end

function BPX_matrix_2D(L::Integer, BCs::AbstractArray{String, 1})
    N = 2^L+1
    h = 1/(N-1)
    DOF_x, DOF_y = N - 2 + sum(BCs[1:2] .== "N"), N - 2 + sum(BCs[3:4] .== "N")
    B = spdiagm(0 => ones(DOF_x*DOF_y))
    for l=1:(L-1)
        ψ_x, χ_x = BPX_weights(L, l, BCs[1:2])
        P_x = BPX_layer_half_matrix(L, ψ_x, χ_x, BCs[1:2])
        ψ_y, χ_y = BPX_weights(L, l, BCs[3:4])
        P_y = BPX_layer_half_matrix(L, ψ_y, χ_y, BCs[3:4])
        P = kron(P_y, P_x)
        # B += (P*transpose(P))*h
        N_ = 2^l+1
        h_ = 1/(N_-1)
        B += (P*transpose(P))*(h/h_)
    end
    return B
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
##==============================================================================
## Abstract Types
##==============================================================================

abstract type InferenceVariable end

##==============================================================================
## VariableNodeData
##==============================================================================

"""
$(TYPEDEF)
Data container for solver-specific data.

  ---
Fields:
$(TYPEDFIELDS)
"""
mutable struct VariableNodeData{T<:InferenceVariable}
    val::Array{Float64,2}
    bw::Array{Float64,2}
    BayesNetOutVertIDs::Array{Symbol,1}
    dimIDs::Array{Int,1} # Likely deprecate
    dims::Int
    eliminated::Bool
    BayesNetVertID::Symbol #  Union{Nothing, }
    separator::Array{Symbol,1}
    softtype::T
    initialized::Bool
    inferdim::Float64
    ismargin::Bool
    dontmargin::Bool
    solveInProgress::Int
    solvedCount::Int
    events::Dict{Symbol,Threads.Condition}
    VariableNodeData{T}() where {T <:InferenceVariable} =
    new{T}(zeros(1,1), zeros(1,1), Symbol[], Int[], 0, false, :NOTHING, Symbol[], T(), false, 0.0, false, false, 0, 0, Dict{Symbol,Threads.Condition}())
    VariableNodeData{T}(val::Array{Float64,2},
                        bw::Array{Float64,2},
                        BayesNetOutVertIDs::Array{Symbol,1},
                        dimIDs::Array{Int,1},
                        dims::Int,eliminated::Bool,
                        BayesNetVertID::Symbol,
                        separator::Array{Symbol,1},
                        softtype::T,
                        initialized::Bool,
                        inferdim::Float64,
                        ismargin::Bool,
                        dontmargin::Bool,
                        solveInProgress::Int=0,
                        solvedCount::Int=0,
                        events::Dict{Symbol,Threads.Condition}=Dict{Symbol,Threads.Condition}()) where T <: InferenceVariable =
                            new{T}(val,bw,BayesNetOutVertIDs,dimIDs,dims,
                                   eliminated,BayesNetVertID,separator,
                                   softtype::T,initialized,inferdim,ismargin,
                                   dontmargin, solveInProgress, solvedCount, events)
end

##------------------------------------------------------------------------------
## Constructors

VariableNodeData(val::Array{Float64,2},
                 bw::Array{Float64,2},
                 BayesNetOutVertIDs::Array{Symbol,1},
                 dimIDs::Array{Int,1},
                 dims::Int,eliminated::Bool,
                 BayesNetVertID::Symbol,
                 separator::Array{Symbol,1},
                 softtype::T,
                 initialized::Bool,
                 inferdim::Float64,
                 ismargin::Bool,
                 dontmargin::Bool,
                 solveInProgress::Int=0,
                 solvedCount::Int=0) where T <: InferenceVariable =
                   VariableNodeData{T}(val,bw,BayesNetOutVertIDs,dimIDs,dims,
                                       eliminated,BayesNetVertID,separator,
                                       softtype::T,initialized,inferdim,ismargin,
                                       dontmargin, solveInProgress, solvedCount)

VariableNodeData(softtype::T) where T <: InferenceVariable =
    VariableNodeData{T}(zeros(1,1), zeros(1,1), Symbol[], Int[], 0, false, :NOTHING, Symbol[], softtype, false, 0.0, false, false, 0, 0)

##==============================================================================
## PackedVariableNodeData.jl
##==============================================================================

"""
$(TYPEDEF)
Packed VariabeNodeData structure for serializing DFGVariables.

  ---
Fields:
$(TYPEDFIELDS)
"""
mutable struct PackedVariableNodeData
    vecval::Array{Float64,1}
    dimval::Int
    vecbw::Array{Float64,1}
    dimbw::Int
    BayesNetOutVertIDs::Array{Symbol,1} # Int
    dimIDs::Array{Int,1}
    dims::Int
    eliminated::Bool
    BayesNetVertID::Symbol # Int
    separator::Array{Symbol,1} # Int
    softtype::String
    initialized::Bool
    inferdim::Float64
    ismargin::Bool
    dontmargin::Bool
    solveInProgress::Int
    solvedCount::Int
    PackedVariableNodeData() = new()
    PackedVariableNodeData(x1::Vector{Float64},
                         x2::Int,
                         x3::Vector{Float64},
                         x4::Int,
                         x5::Vector{Symbol}, # Int
                         x6::Vector{Int},
                         x7::Int,
                         x8::Bool,
                         x9::Symbol, # Int
                         x10::Vector{Symbol}, # Int
                         x11::String,
                         x12::Bool,
                         x13::Float64,
                         x14::Bool,
                         x15::Bool,
                         x16::Int,
                         solvedCount::Int) = new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16, solvedCount)
end

##==============================================================================
## PointParametricEst
##==============================================================================

##------------------------------------------------------------------------------
## AbstractPointParametricEst interface
##------------------------------------------------------------------------------

abstract type AbstractPointParametricEst end

##------------------------------------------------------------------------------
## MeanMaxPPE
##------------------------------------------------------------------------------
"""
    $TYPEDEF

Data container to store Parameteric Point Estimate (PPE) for mean and max.
"""
struct MeanMaxPPE <: AbstractPointParametricEst
    solverKey::Symbol #repeated because of Sam's request
    suggested::Vector{Float64}
    max::Vector{Float64}
    mean::Vector{Float64}
    lastUpdatedTimestamp::DateTime
end
##------------------------------------------------------------------------------
## Constructors

MeanMaxPPE(solverKey::Symbol, suggested::Vector{Float64}, max::Vector{Float64},mean::Vector{Float64}) = MeanMaxPPE(solverKey, suggested, max, mean, now())

##==============================================================================
## DFG Variables
##==============================================================================

##------------------------------------------------------------------------------
## DFGVariable lv2
##------------------------------------------------------------------------------
"""
$(TYPEDEF)
Complete variable structure for a DistributedFactorGraph variable.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGVariable{T<:InferenceVariable} <: AbstractDFGVariable
    """Variable label, e.g. :x1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref), [`setTimestamp`](@ref)"""
    timestamp::DateTime
    """Nano second time, for more resolution on timestamp (only subsecond information)"""
    nstime::Nanosecond
    """Variable tags, e.g [:POSE, :VARIABLE, and :LANDMARK].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Dictionary of parametric point estimates keyed by solverDataDict keys
    Accessors: [`addPPE!`](@ref), [`updatePPE!`](@ref), and [`deletePPE!`](@ref)"""
    ppeDict::Dict{Symbol, <: AbstractPointParametricEst}
    """Dictionary of solver data. May be a subset of all solutions if a solver key was specified in the get call.
    Accessors: [`addVariableSolverData!`](@ref), [`updateVariableSolverData!`](@ref), and [`deleteVariableSolverData!`](@ref)"""
    solverDataDict::Dict{Symbol, VariableNodeData{T}}
    """Dictionary of small data associated with this variable.
    Accessors: [`getSmallData`](@ref), [`setSmallData!`](@ref)"""
    smallData::Dict{String, String}#Ref{Dict{String, String}} #why was Ref here?
    """Dictionary of large data associated with this variable.
    Accessors: [`addBigDataEntry!`](@ref), [`getBigDataEntry`](@ref), [`updateBigDataEntry!`](@ref), and [`deleteBigDataEntry!`](@ref)"""
    bigData::Dict{Symbol, AbstractBigDataEntry}
    """Solvable flag for the variable.
    Accessors: [`getSolvable`](@ref), [`setSolvable!`](@ref)"""
    solvable::Base.RefValue{Int}
end

##------------------------------------------------------------------------------
## Constructors

"""
    $SIGNATURES
The default DFGVariable constructor.
"""
DFGVariable(label::Symbol, softtype::T;
            timestamp::DateTime=now(),
            nstime::Nanosecond = Nanosecond(0),
            tags::Set{Symbol}=Set{Symbol}(),
            estimateDict::Dict{Symbol, <: AbstractPointParametricEst}=Dict{Symbol, MeanMaxPPE}(),
            solverDataDict::Dict{Symbol, VariableNodeData{T}}=Dict{Symbol, VariableNodeData{T}}(),
            smallData::Dict{String, String}=Dict{String, String}(),
            bigData::Dict{Symbol, AbstractBigDataEntry}=Dict{Symbol,AbstractBigDataEntry}(),
            solvable::Int=1) where {T <: InferenceVariable} =
    DFGVariable{T}(label, timestamp, nstime, tags, estimateDict, solverDataDict, smallData, bigData, Ref(solvable))

DFGVariable(label::Symbol,
            solverData::VariableNodeData{T};
            timestamp::DateTime=now(),
            nstime::Nanosecond = Nanosecond(0),
            tags::Set{Symbol}=Set{Symbol}(),
            estimateDict::Dict{Symbol, <: AbstractPointParametricEst}=Dict{Symbol, MeanMaxPPE}(),
            smallData::Dict{String, String}=Dict{String, String}(),
            bigData::Dict{Symbol, AbstractBigDataEntry}=Dict{Symbol,AbstractBigDataEntry}(),
            solvable::Int=1) where {T <: InferenceVariable} =
    DFGVariable{T}(label, timestamp, nstime, tags, estimateDict, Dict{Symbol, VariableNodeData{T}}(:default=>solverData), smallData, bigData, Ref(solvable))

Base.getproperty(x::DFGVariable,f::Symbol) = begin
    if f == :solvable
        getfield(x,f)[]
    else
        getfield(x,f)
    end
end

Base.setproperty!(x::DFGVariable,f::Symbol, val) = begin
    if f == :solvable
        getfield(x,f)[] = val
    else
        setfield!(x,f,val)
    end
end

##------------------------------------------------------------------------------
# TODO: can't see the reason to overwrite copy, leaving it here for now
# function Base.copy(o::DFGVariable)::DFGVariable
#     return DFGVariable(o.label, getSofttype(o)(), tags=copy(o.tags), estimateDict=copy(o.estimateDict),
#                         solverDataDict=copy(o.solverDataDict), smallData=copy(o.smallData),
#                         bigData=copy(o.bigData), solvable=getSolvable(o))
# end

##------------------------------------------------------------------------------
## DFGVariableSummary lv1
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Summary variable structure for a DistributedFactorGraph variable.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct DFGVariableSummary <: AbstractDFGVariable
    """Variable label, e.g. :x1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable timestamp.
    Accessors: [`getTimestamp`](@ref), [`setTimestamp`](@ref)"""
    timestamp::DateTime
    """Variable tags, e.g [:POSE, :VARIABLE, and :LANDMARK].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
    """Dictionary of parametric point estimates keyed by solverDataDict keys
    Accessors: [`addPPE!`](@ref), [`updatePPE!`](@ref), and [`deletePPE!`](@ref)"""
    ppeDict::Dict{Symbol, <:AbstractPointParametricEst}
    """Symbol for the softtype for the underlying variable.
    Accessor: [`getSofttype`](@ref)"""
    softtypename::Symbol
    """Dictionary of large data associated with this variable.
    Accessors: [`addBigDataEntry!`](@ref), [`getBigDataEntry`](@ref), [`updateBigDataEntry!`](@ref), and [`deleteBigDataEntry!`](@ref)"""
    bigData::Dict{Symbol, AbstractBigDataEntry}
end

##------------------------------------------------------------------------------
## SkeletonDFGVariable.jl
##------------------------------------------------------------------------------

"""
$(TYPEDEF)
Skeleton variable structure for a DistributedFactorGraph variable.

  ---
Fields:
$(TYPEDFIELDS)
"""
struct SkeletonDFGVariable <: AbstractDFGVariable
    """Variable label, e.g. :x1.
    Accessor: [`getLabel`](@ref)"""
    label::Symbol
    """Variable tags, e.g [:POSE, :VARIABLE, and :LANDMARK].
    Accessors: [`getTags`](@ref), [`mergeTags!`](@ref), and [`removeTags!`](@ref)"""
    tags::Set{Symbol}
end

SkeletonDFGVariable(label::Symbol) = SkeletonDFGVariable(label, Set{Symbol}())

##==============================================================================
# Define variable levels
##==============================================================================
const VariableDataLevel0 = Union{DFGVariable, DFGVariableSummary, SkeletonDFGVariable}
const VariableDataLevel1 = Union{DFGVariable, DFGVariableSummary}
const VariableDataLevel2 = Union{DFGVariable}

##==============================================================================
## Convertion constructors
##==============================================================================

DFGVariableSummary(v::DFGVariable) =
        DFGVariableSummary(v.label, v.timestamp, deepcopy(v.tags), deepcopy(v.ppeDict), Symbol(typeof(getSofttype(v))), v.bigData)

SkeletonDFGVariable(v::VariableDataLevel1) =
            SkeletonDFGVariable(v.label, deepcopy(v.tags))


text:
function modifier(ability)
    return floor((ability - 10) / 2)
end

function ability()
    return rand(3:18)
end

mutable struct DNDCharacter
    strength::Int
    dexterity::Int
    constitution::Int
    intelligence::Int
    wisdom::Int
    charisma::Int
    hitpoints::Int
end

function DNDCharacter()
    strength = ability()
    dexterity = ability()
    constitution = ability()
    intelligence = ability()
    wisdom = ability()
    charisma = ability()
    hitpoints = 10 + modifier(constitution)

    return DNDCharacter(
        strength,
        dexterity,
        constitution,
        intelligence,
        wisdom,
        charisma,
        hitpoints,
    )
end


text:
### Hastings & Powell 1991
## Examining behavior when varying b2 value
using LaTeXStrings
### Figure S1 (Bifurcation diagram)

## Wrap in function
# b1val : value of b1
# b2low : lower limit of b2 values to examine
# b2up : upper limit of b2 values to examine
# out : desired output plot, either "maxloc" for local maxima only (default) or "all" for all values
function b2_behavior(;b1val::Float64, b2low::Float64=1.5, b2up::Float64=3.2,out::String="maxloc")
  ## Define initial parameters
  p = a1, a2, b1, b2, d1, d2 = [5.0, 0.1, b1val, b2low, 0.4, 0.01]
  u0 = [1.0, 1.0, 1.0]
  tspan = (0.0, 10000.0)
  # b1 values interval
  bs = collect(b2low:0.01:b2up)
  # number of values to collect
  N = 1000

  ## Loop to select all values of z for all b1 values
  Z = zeros(Float64, (N, length(bs)))
  Zmax = zeros(Float64, length(bs))
  Zmin = zeros(Float64, length(bs))
  errors = zeros(Float64, length(bs))
  for j in 1:length(bs)
    # Change b1 value
    p[4] = b2 = bs[j]
    # Solve for given b1 value
    prob = ODEProblem(parameterized_model!, u0, tspan, p)
    sol = solve(prob)
    # Save values
    if length(sol[3,:]) > N
      Z[:,j] = sol[3, end-(N-1):end]
      Zmax[j] = maximum(Z[:,j])
      Zmin[j] = minimum(Z[:,j])
    else
      # Safeguard for when true solution is unstable
      errors[j] = j
    end
  end
  Z
  Zmax
  Zmin
  errors[errors .> 0]

  ## Select only local maxima
  Zmaxloc = zeros(Float64, size(Z))
  for i in 2:(N-1), j in 1:length(bs)
    # Select maxima values(greater than neighbors)
    if Z[i, j] > Z[i-1, j] && Z[i, j] > Z[i+1, j]
      # Select first local maxima only with arbitrary interval
      Zrange = Zmax[j] - Zmin[j]
      if Z[i, j] > (Zmin[j] + 0.66*(Zrange))
         Zmaxloc[i, j] = Z[i, j]
      end
    end
  end
  Zmaxloc

  ## Visualize results
  # Vector of b1 values repeated N times
  B = repeat(bs, inner = N)
  # Match local maxima with respective b1 value
  pts_all = unique(hcat(B, vec(Z)); dims = 1)
  pts_maxloc = unique(hcat(B, vec(Zmaxloc)); dims = 1)
  # Remove null values
  pts_all = pts_all[pts_all[:, 2] .> 0.00, :]
  pts_maxloc = pts_maxloc[pts_maxloc[:,2] .> 0.00, :]
  # Select desired output (out parameter value)
  if out == "maxloc" # default
    results = pts_maxloc
  end
  if out == "all" # default
    results = pts_all
  end
  # Plot results (bifurcation diagram)
  scatter(results[:,1], results[:,2],
          markersize = 1, c = :black,
          xlabel = L"b_2", xlim = (b2low,b2up), xticks = b2low:0.25:b2up,
          ylabel = L"z_{max}", yticks = (0.0:2.0:maximum(results[:, 2]), string.(0.0:2.0:maximum(results[:, 2]))),
          legend = false, grid = :none, framestyle = :box,
          guidefontsize = 10)
end

# Examine behavior when changing b2 at b1 = 3.0
figS1_maxloc = b2_behavior(b1val = 3.0, b2low = 1.5, b2up = 3.2,
                           out = "maxloc") # default parameters
figS1_all = b2_behavior(b1val = 3.0, out = "all")

## Export figures
savefig(figS1_maxloc, joinpath("..", "article", "figures", "figS1.pdf"))

# Additionnal examinations
b2_behavior(b1val = 2.0)
b2_behavior(b1val = 4.0)
b2_behavior(b1val = 5.0)
b2_behavior(b1val = 6.0)
b2_behavior(b1val = 1.5)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Coalesce(value)

Replaces all missing values from the table with `value`.

# Examples

```julia
Coalesce(0)
```

## Notes

* The transform can alter the element type of columns from `Union{Missing,T}` to `T`.
"""
struct Coalesce{T} <: Colwise
  value::T
end

isrevertible(::Type{<:Coalesce}) = true

colcache(::Coalesce, x) = findall(ismissing, x)

colapply(tramsform::Coalesce, x, c) = coalesce.(x, tramsform.value)

colrevert(::Coalesce, x, c) = [i ∈ c ? missing : x[i] for i in 1:length(x)]


text:
using DelimitedFiles

"Load numbers from file."
function loaddata(file::String)
    return readdlm(file, '\t', Int, '\n')
end

"Find the two entries that sum to 2020."
function findtwo(data::Array{Int64,2})
    return [x * y for x in data, y in data if x + y == 2020][1]
end

"Find the three entries that sum to 2020."
function findthree(data::Array{Int64,2})
    return [x * y * z for x in data, y in data, z in data if x + y + z == 2020][1]
end

function main()
    data = loaddata("../assets/data/01.txt")
    println(typeof(data))

    # part one
    println("The product of the two numbers is ", findtwo(data), ".")

    # part two
    println("The product of the three numbers is ", findthree(data), ".")
end

main()


text:
module Metatheory

using RuntimeGeneratedFunctions
using Base.Meta
using Reexport

include("docstrings.jl")

RuntimeGeneratedFunctions.init(@__MODULE__)

include("options.jl")

macro log(args...)
    quote options.verbose && @info($(args...)) end |> esc
end

export options

include("Util/Util.jl")
using .Util
export Util

include("rgf.jl")
export @metatheory_init

include("Patterns/Patterns.jl")
@reexport using .Patterns

include("ematch_compiler.jl")
@reexport using .EMatchCompiler

include("Rules/Rules.jl")
@reexport using .Rules

include("Classic/Classic.jl")
@reexport using .Classic: gettheory

include("EGraphs/EGraphs.jl")
@reexport using .EGraphs

include("Library/Library.jl")
export Library

end # module


text:
function calc_cstar!(c,cm1,Δt,Δt2)
    @mthreads for j in TRANGE
        r = (Δt/Δt2)
        @inbounds @simd for i in REAL_RANGES[j]
            cstar = muladd(r,(c[i]-cm1[i]),c[i])
            cm1[i] = c[i]
            #c[i] = max(0.0,cstar)
            c[i] = cstar
        end
    end
end

function average_c!(c,Δ2)
    fullfourier!(c)
    ind = CartesianIndices((Base.OneTo(NY),Base.OneTo(NZ)))
    @mthreads for jk in ind
        j,k = Tuple(jk)
            @inbounds @msimd for i in Base.OneTo(NX)
                G = Gaussfilter(Δ2,i,j,k)
                c[i,j,k] *= G
            end
    end

    real!(c)
    return nothing
end


text:
using Pkg
Pkg.develop(PackageSpec(path=splitdir(@__DIR__)[1]))
pkg"instantiate"
using Documenter, GridArrays

const format = Documenter.HTML(
        prettyurls = ("deploy" in ARGS),
    )

makedocs(sitename="GridArrays.jl",
    modules = [GridArrays],
    format = format,
    pages = [
        "Home" => "index.md",
        "Manual" => "man/GridArrays.md"
        ],
    doctest=true
)

if "deploy" in ARGS && Sys.ARCH === :x86_64 && Sys.KERNEL === :Linux
    deploydocs(
        repo = "github.com/JuliaApproximation/GridArrays.jl.git",
    )
end


text:
using Test
using ManifoldsBase
@testset "ManifoldsBase" begin

    if VERSION >= v"1.1"
        num_ambiguities = length(Test.detect_ambiguities(ManifoldsBase))
        #num_ambiguities > 0 && @warn "The number of ambiguities in ManifoldsBase is $(num_ambiguities)."
        if VERSION >= v"1.6"
            @test num_ambiguities <= 1
        else
            @test num_ambiguities == 0
        end
    end
    include("allocation.jl")
    include("numbers.jl")
    include("bases.jl")
    include("decorator_manifold.jl")
    include("empty_manifold.jl")
    include("errors.jl")
    include("default_manifold.jl")
    include("complex_manifold.jl")
    include("validation_manifold.jl")
    include("embedded_manifold.jl")
    include("power.jl")
    include("domain_errors.jl")
    include("vector_transport_along.jl")
end


text:
dbinom(k, n, p) = Distributions.pdf(Distributions.Binomial(n, p), k)

qnorm(p) = Distributions.quantile(Distributions.Normal(0, 1), p)

function findgrid(prior, l, u, n)

  dp = 1 / (n + 1)
  modprior_(p) = .9 * prior(p) + .1 * Distributions.pdf(Distributions.Beta(1, 1), p)
  zmodprior = QuadGK.quadgk(modprior_, l, u)[1]
  modprior(p) = modprior_(p) / zmodprior * (p < l ? 0 : 1)
  zprior = QuadGK.quadgk(prior, l, u)[1]
  modcdf(p) = QuadGK.quadgk(modprior, l, p)[1]
  cprior(p) = prior(p) / zprior * (p < l ? 0 : 1)
  cdf(p) = QuadGK.quadgk(cprior, l, p)[1]
  pivots  = [l + 10^-6.0]
  modcdfs = [0.0]
  cdfs    = [0.0]
  for i in 2:(n + 1)
    push!(pivots, Roots.fzero(p -> modcdf(p) - dp - modcdfs[i - 1], pivots[i - 1], 1))
    push!(modcdfs, modcdf(pivots[i]))
    push!(cdfs, cdf(pivots[i]))
  end
  pivots = (pivots[2:(n + 1)] .+ pivots[1:n]) / 2 # midpoint pivots!
  dcdf = cdfs[2:(n + 1)] - cdfs[1:n] # vector of first order differences
  dcdf = dcdf / sum(dcdf)
  return pivots, dcdf

end

macro checkprob(x...)

  tmp = Expr(:block)
  for xx in x
    name = string(xx)
    push!(tmp.args, :($xx < 0 ? error(@sprintf("%s must be >= 0, is %.5f", $name, $xx)) : true))
    push!(tmp.args, :($xx > 1 ? error(@sprintf("%s must be <= 1, is %.5f", $name, $xx)) : true))
  end
  return esc(tmp)

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
abstract type AbstractReservoirMemoryMachine <: AbstractEchoStateNetwork end

mutable struct RMM{T, S<:AbstractArray{T}, I, B, F, N} <: AbstractReservoirMemoryMachine
    res_size::I
    in_size::I
    out_size::I
    train_data::S
    alpha::T
    nla_type::N
    activation::F
    W::S
    W_in::S
    states::S
    extended_states::B
    W_out::S
    memory_size::I
    write_matrix::S
    read_matrix::S
end

"""
    RMM(W::AbstractArray{T}, in_data::AbstractArray{T}, out_data::AbstractArray{T}, W_in::AbstractArray{T},
    memory_size::Int, beta::Float64 [, activation::Any, alpha::T, nla_type::NonLinearAlgorithm, extended_states::Bool])

Create a Reservoir Memory Machine struct, as described in [1].

[1] Paaßen, Benjamin, and Alexander Schulz. "Reservoir memory machines." arXiv preprint arXiv:2003.04793 (2020).
"""
function RMM(W::AbstractArray{T},
        in_data::AbstractArray{T},
        out_data::AbstractArray{T},
        W_in::AbstractArray{T},
        memory_size::Int,
        beta::Float64;
        activation::Any = tanh,
        alpha::T = 1.0,
        nla_type::NonLinearAlgorithm = NLADefault(),
        extended_states::Bool = false) where T<:AbstractFloat

    in_size = size(in_data, 2)
    out_size = size(out_data, 2)
    res_size = size(W, 1)

    if size(W_in, 1) != res_size
        throw(DimensionMismatch("size(W_in, 1) must be equal to size(W, 1)"))
    elseif size(W_in, 2) != in_size
        throw(DimensionMismatch("size(W_in, 2) must be equal to in_size"))
    end

    states = states_matrix_inplace(W, W_in, collect(in_data'), alpha, activation, extended_states)
    write_matrix, read_matrix = identity_init(in_data, out_data, memory_size, states, beta)#forward_init(states, beta) #identity_init(input, output, memory_size, states, beta)
    memory_states = apply_write(in_data, states, write_matrix, memory_size)
    reads = apply_read(in_data, states, write_matrix, read_matrix, memory_size)
    wout = RMMtrain(in_data, out_data, states, write_matrix, read_matrix, memory_states, memory_size, beta)

    return RMM{T, typeof(in_data),
        typeof(res_size),
        typeof(extended_states),
        typeof(activation),
        typeof(nla_type)}(res_size, in_size, out_size, in_data,
    alpha, nla_type, activation, W, W_in, states, extended_states, wout, memory_size, write_matrix, read_matrix)
end

function RMMtrain(input, output, states, write_matrix, read_matrix, memory_states, memory_size, beta)

    states_pinv = (states*inv(states'*states+beta*Matrix(1.0I, size(states, 2), size(states, 2))))
    last_loss = Inf
    last_vs = nothing

    for epoch=1:30
        reads = apply_read(input, states, write_matrix, read_matrix, memory_size)
        wout = (output'*reads)*inv(add_reg(reads'*reads, beta))
        loss = sqrt(mean(sum((output - reads*wout').^2, dims=2)))

        if loss < last_loss - 1*10^(-3)
            last_loss = loss
            last_vs = (write_matrix, read_matrix, wout)
        else
            write_matrix, read_matrix, wout = last_vs

            break
        end

        write_actions = align_write(input, output, wout, memory_size)
        write_matrix = write_actions'*states_pinv

        a, read_actions = align_read(memory_states, output, memory_size, wout)
        read_matrix = collect(reduce(hcat, read_actions))*states_pinv
    end

    read = apply_read(input, states, write_matrix, read_matrix, memory_size)
    hr = hcat(states, read)
    wout = (output'*hr)*inv(add_reg(hr'*hr, beta))
    return wout

end

"""
    RMMdirect_predict(rmm::AbstractReservoirMemoryMachine, input)

Given the input data return the corresponding predicted output, as described in [1].

[1] Paaßen, Benjamin, and Alexander Schulz. "Reservoir memory machines." arXiv preprint arXiv:2003.04793 (2020).
"""
function RMMdirect_predict(rmm::AbstractReservoirMemoryMachine, input)

    predict_len = size(input, 1)
    predict_states = zeros(Float64, predict_len, rmm.res_size)
    #compute first state
    predict_states[1, :] = rmm.activation.(rmm.W_in*input[1, :])

    if rmm.extended_states == false
        for i=2:predict_len
            predict_states[i, :] = leaky_fixed_rnn(rmm.activation, rmm.alpha, rmm.W, rmm.W_in, predict_states[i-1, :], input[i,:])
        end
    end

    reads = apply_read(input, predict_states, rmm.write_matrix, rmm.read_matrix, rmm.memory_size)
    predict_states = hcat(predict_states, reads)
    output = predict_states*rmm.W_out'
    return output
end

function states_matrix_inplace(W::AbstractArray{Float64},
        W_in::AbstractArray{Float64},
        train_data::AbstractArray{Float64},
        alpha::Float64,
        activation::Function,
        extended_states::Bool)

    train_len = size(train_data, 2)
    res_size = size(W, 1)
    in_size = size(train_data, 1)

    states = zeros(Float64, train_len, res_size)
    states[1, :] = activation.(W_in*train_data[:, 1])

    for i=2:train_len
        states[i, :] = leaky_fixed_rnn(activation, alpha, W, W_in, states[i-1, :], train_data[:, i])
    end

    if extended_states == true
        ext_states = vcat(states, hcat(zeros(Float64, in_size), train_data[:,1:end-1]))
        return ext_states
    else
        return states
    end
end

function apply_write(input, states, write_matrix, memory_size)

    in_size = size(input, 2)
    train_len = size(input, 1)
    write_actions = states*write_matrix'
    memory_states = zeros(Float64, memory_size, in_size, train_len)
    write_loc = 1

    for i=1:train_len
        if i > 1
            memory_states[:, :, i] = memory_states[:, :, i-1]
        end

        if write_actions[i] > 0
            memory_states[write_loc, :, i] = input[i,:]
            write_loc += 1
            if write_loc >= memory_size
                write_loc = 1
            end
        end
    end
    return memory_states
end

function apply_read(input, states, write_matrix, read_matrix, memory_size)

    train_len = size(input, 1)

    write_actions = states*write_matrix'
    read_actions =states*read_matrix'
    read_actions = mapslices(argmax, read_actions, dims=2)
    memory = zeros(memory_size, size(input, 2))

    write_loc = 1
    read_loc = 1
    reads = Array{Float64}(undef, train_len, size(input, 2))
    for i=1:train_len
        if write_actions[i] > 0
            memory[write_loc, :] = input[i,:]
            write_loc += 1
            if write_loc >= memory_size
                write_loc = 1
            end
        end

        if read_actions[i] == 2
            read_loc += 1
            if read_loc >= memory_size
                read_loc = 1
            end
        elseif read_actions[i] == 3
            read_loc = 1
        end

        reads[i, :] =  memory[read_loc, :]
    end

    return reads
end

function align_write(input, output, wout, memory_size)#  permit_duplicates always false. Input is Transpose than my implementation

    train_len = size(input, 1)
    delta = pairwise(Euclidean(), (input*transpose(wout)), output, dims = 1)
    delta_min = minimum(delta, dims = 1)

    writes = Array{Int}(undef, train_len)

    for i=1:train_len
        small_min =  findall(delta[:,i] .< minimum(delta[:,i])+1*10^(-3))[1]
        writes[i] = small_min#push!(writes, small_min)
    end

    write_loc = 1
    write_actions = zeros(Int, train_len)

    for i=1:train_len
        if i in writes
            write_loc += 1
            if write_loc >= memory_size
                write_loc = 1
            end
            write_actions[i] = 1
        else
            write_actions[i] = -1
        end
    end

    return write_actions
end

function align_read(memory_states, output, memory_size, wout)

    train_len = size(output, 1)
    delta = zeros(memory_size, train_len)
    for i=1:memory_size
        for j=1:train_len
            delta[i, j] = norm(wout*memory_states[i, :, j] - output[j, :])
        end
    end

    dyn_matrix = zeros(memory_size, train_len)
    for i=1:memory_size
        dyn_matrix[i, train_len-1] = delta[i, train_len-1]
    end

    for i in train_len-2: -1:1
        for j=1:memory_size-1
            dyn_matrix[j, i] = delta[j, i] +min(dyn_matrix[j, i+1], min(dyn_matrix[j+1, i+1], dyn_matrix[1, i+1]))
        end
        dyn_matrix[memory_size-1, i] = delta[memory_size-1, i] + min(dyn_matrix[memory_size-1, i+1], dyn_matrix[1, i+1])
    end

    actions = Array{Float64}[]
    if dyn_matrix[1, 1] < dyn_matrix[2, 1]+ 1*10^(-3)
        k = 1
        d = dyn_matrix[1, 1]
        push!(actions, [1, 0, 0])
    else
        k = 2
        d = dyn_matrix[2, 1]
        push!(actions, [0, 1, 0])
    end

    for i=1:train_len-1
        if delta[k, i] + dyn_matrix[k, i+1] < dyn_matrix[k, i] + 1*10^(-3)
            push!(actions, [1, 0, 0])
        elseif k < memory_size-1 && delta[k, i] + dyn_matrix[k+1, i+1] < dyn_matrix[k, i] + 1*10^(-3)
            push!(actions, [0, 1, 0])
            k+=1
        elseif delta[k, i] + dyn_matrix[1, i+1] < dyn_matrix[k, i] + 1*10^(-3)
            if k == memory_size-1
                push!(actions, [0, 1, 0])
            else
                push!(actions, [0, 0, 1])
            end
            k = 1
        else
            println("error")
        end
    end
    return min(dyn_matrix[1, 1], dyn_matrix[2, 1]), actions
end

function identity_init(input, output, memory_size, states, beta)

    states_pinv = (states*inv(states'*states+beta*Matrix(1.0I, size(states, 2), size(states, 2))))

    in_size = size(input, 2)
    train_len = size(input, 1)
    wout = Matrix{Float64}(I, in_size, in_size)
    #write_actions = []

    waj = align_write(input, output, wout, memory_size)
    #push!(write_actions, waj)
    write_matrix = collect(waj')*states_pinv

    memory_states = apply_write(input, states, write_matrix, memory_size)
    a, raj = align_read(memory_states, output, memory_size, wout)
    read_matrix = collect(reduce(hcat, raj))*states_pinv

    return write_matrix, read_matrix
end

function forward_init(states, beta)

    write_actions = ones(Float64, size(states, 1), 1)
    read_actions = zeros(Float64, size(states, 1), 3)
    read_actions[:,1] .= 1
    states_pinv = (states*inv(states'*states+beta*Matrix(1.0I, size(states, 2), size(states, 2))))
    write_matrix = write_actions'*states_pinv
    read_matrix  = read_actions'*states_pinv

    return write_matrix, read_matrix
end


text:
module Simm

import ..Clarus

function history(;params...)
  return Clarus.api_request("SIMM","History";params...)
end

function impact(;params...)
  return Clarus.api_request("SIMM","Impact";params...)
end

function margin(;params...)
  return Clarus.api_request("SIMM","Margin";params...)
end

function sensitivity(;params...)
  return Clarus.api_request("SIMM","Sensitivity";params...)
end

end #MODULE-END


text:
f = float16(2.)
g = float16(1.)

@test f >= g
@test f > g
@test g < f
@test g <= g
@test all([g g] .< [f f])
@test all([g g] .<= [f f])
@test all([f f] .> [g g])
@test all([f f] .>= [g g])
@test isless(g, f)
@test !isless(f, g)

@test -f === float16(-2.)

@test f+g === float16(3f0)
@test f-g === float16(1f0)
@test f*g === float16(2f0)
@test f/g === float16(2f0)
@test f^g === float16(2f0)
@test f^-g === float16(0.5f0)

@test f + 2 === float32(4f0)
@test f - 2 === float32(0f0)
@test f*2 === float32(4f0)
@test f/2 === float32(1f0)
@test f + 2. === 4.
@test f - 2. === 0.
@test f*2. === 4.
@test f/2. === 1.

@test_approx_eq sin(f) sin(2f0)

@test isnan(NaN16)
@test isnan(-NaN16)
@test !isnan(Inf16)
@test !isnan(-Inf16)
@test !isnan(float16(2.6))
@test NaN16 != NaN16
@test isequal(NaN16, NaN16)
@test repr(NaN16) == "NaN16"
@test sprint(showcompact, NaN16) == "NaN"

@test isinf(Inf16)
@test isinf(-Inf16)
@test !isinf(NaN16)
@test !isinf(-NaN16)
@test !isinf(float16(2.6))
@test Inf16 == Inf16
@test Inf16 != -Inf16
@test -Inf16 < Inf16
@test isequal(Inf16, Inf16)
@test repr(Inf16) == "Inf16"
@test sprint(showcompact, Inf16) == "Inf"

for z1 in (float16(0.0), float16(-0.0)), z2 in (float16(0.0), float16(-0.0))
    @test z1 == z2
    @test isequal(z1, z1)
    @test z1 === z1
    for elty in (Float32, Float64)
        z3 = convert(elty, z2)
        @test z1==z3
    end
end

@test float16(2.5) == float16(2.5)
@test float16(2.5) != float16(2.6)
@test isequal(float16(0.0), float16(0.0))
@test !isequal(float16(-0.0), float16(0.0))
@test !isequal(float16(0.0), float16(-0.0))

@test isnan(reinterpret(Float16,0x7c01))
@test !isinf(reinterpret(Float16,0x7c01))

@test nextfloat(Inf16) === Inf16
@test prevfloat(-Inf16) === -Inf16

# rounding in conversions
let
    for f in [.3325f0, -.3325f0]
        f16 = float16(f)
        # need to round away from 0. make sure we picked closest number.
        @test abs(f-f16) < abs(f-nextfloat(f16))
        @test abs(f-f16) < abs(f-prevfloat(f16))
    end
    # halfway between and last bit is 1
    f = reinterpret(Float32,                           0b00111110101010100011000000000000)
    @test float32(float16(f)) === reinterpret(Float32, 0b00111110101010100100000000000000)
    # halfway between and last bit is 0
    f = reinterpret(Float32,                           0b00111110101010100001000000000000)
    @test float32(float16(f)) === reinterpret(Float32, 0b00111110101010100000000000000000)
end

# issue #5948
@test string(reinterpret(Float16, 0x7bff)) == "65504.0"

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia