text:
export
ruleSPProbitOutNG,
ruleSPProbitIn1PN,
ruleEPProbitIn1BG,
ruleEPProbitIn1CG,
ruleEPProbitIn1PG

function ruleSPProbitOutNG(msg_out::Nothing,
                           msg_in1::Message{F, Univariate}) where F<:Gaussian

    d_in1 = convert(ProbabilityDistribution{Univariate, GaussianMeanVariance}, msg_in1.dist)

    p = normcdf(d_in1.params[:m] / sqrt(1 + d_in1.params[:v]))
    isnan(p) && (p = 0.5)

    Message(Univariate, Bernoulli, p=p)
end

function ruleSPProbitIn1PN(msg_out::Message{PointMass, Univariate},
                           msg_in1::Nothing)
    p = msg_out.dist.params[:m]
    if p == 1.0
        log_pdf = normlogcdf
    elseif p == 0.0
        log_pdf = (z)->normlogcdf(-z)
    else
        log_pdf = (z)->log(p*normcdf(z) + (1.0-p)*normcdf(-z))
    end

    return Message(Univariate, Function, log_pdf=log_pdf)
end

function ruleSPProbitIn1BN(msg_out::Message{Bernoulli, Univariate},
                           msg_in1::Nothing)

    p = msg_out.dist.params[:p]
    if p == 1.0
        log_pdf = normlogcdf
    elseif p == 0.0
        log_pdf = (z)->normlogcdf(-z)
    else
        log_pdf = (z)->log(p*normcdf(z) + (1.0-p)*normcdf(-z))
    end

    return Message(Univariate, Function, log_pdf=log_pdf)
end

function ruleEPProbitIn1BG(msg_out::Message{Bernoulli, Univariate},
                           msg_in1::Message{F, Univariate}) where F<:Gaussian

    # Calculate approximate (Gaussian) message towards i[:in1]
    # The approximate message is an 'expectation' under the context (cavity distribution) encoded by incoming message msg_cavity.
    # Propagating the resulting approximate msg through the factor graph results in the expectation propagation (EP) algorithm.
    # Approximation procedure:
    #  1. Calculate exact (non-Gaussian) message towards i[:in1].
    #  2. Combine exact outbound msg on i[:in1] with exact inbound msg (cavity distribution) to find exact marginal.
    #  3. Approximate the exact (non-Gaussian) marginal with a Gaussian one using moment matching, under the constraint VAR[approximate marginal] > VAR[cavity].
    #  4. Calculate back the Gaussian outbound msg on i[:in1] that yields this approximate Gaussian marginal.
    # IMPORTANT NOTES:
    #  - This calculation results in an implicit cycle in the factor graph since the outbound message depends on the inbound message (cavity dist.).

    # Shorthand notations
    p = msg_out.dist.params[:p]
    (μ, σ2) = unsafeMeanCov(msg_in1.dist) # Moments of cavity distribution
    (ξ, γ) = unsafeWeightedMeanPrecision(msg_in1.dist)

    # Calculate first and second moment (mp_1, mp_2) of the 'true' marginal p(x) on edge connected to i[:in1]
    # p(x) = f(x) / Z
    # f(x) = (1-p)*N(x|μ,σ2) + (2p-1)*Φ(x)*N(x|μ,σ2)
    #      = (1-p)*N(x|μ,σ2) + (2p-1)*Φ(z)*(Φ(x)*N(x|μ,σ2)/Φ(z))
    #      = (1-p)*N(x|μ,σ2) + (2p-1)*Φ(z)*g(x)
    # See paper for detailed derivation

    z = μ / sqrt(1 + σ2)
    N = exp(-0.5*z^2)./sqrt(2*pi) # 𝓝(z)

    # Moments of g(x)
    mg_1 = normcdf(z)*μ + σ2*N / sqrt(1+σ2)  # First moment of g
    mg_2 = 2*μ*mg_1 + normcdf(z)*(σ2 - μ^2) - σ2^2*z*N / (1+σ2)  # Second moment of g

    # Moments of f(x) (exact marginal)
    Z = 1 - p + (2*p-1)*normcdf(z)
    mp_1 = ((1-p)*μ + (2*p-1)*mg_1) / Z
    mp_2 = ((1-p)*(μ^2+σ2) + (2*p-1)*mg_2) / Z

    # Calculate Gaussian marginal with identical first and second moments (moment matching approximation)
    marginal_v = mp_2 - mp_1^2
    marginal_v = clamp(marginal_v, tiny, σ2-tiny) # ensure variance of marginal is not larger than variance of cavity distribution
    marginal_w = 1.0 / marginal_v
    marginal_xi = marginal_w * mp_1

    # Calculate the approximate message towards i[:in1]
    outbound_dist_w = marginal_w - γ
    outbound_dist_xi = marginal_xi - ξ

    return Message(Univariate, GaussianWeightedMeanPrecision, xi=outbound_dist_xi, w=outbound_dist_w)
end

function ruleEPProbitIn1PG(msg_out::Message{PointMass, Univariate}, msg_in1::Message{F, Univariate}) where F<:Gaussian
    p = msg_out.dist.params[:m]
    isnan(p) && (p = 0.5)
    (0.0 <= p <= 1.0) || error("Binary input $p must be between 0 and 1")

    return ruleEPProbitIn1BG(Message(Univariate, Bernoulli, p=p), msg_in1)
end

function ruleEPProbitIn1CG(msg_out::Message{Categorical, Univariate}, msg_in1::Message{F, Univariate}) where F<:Gaussian
    (length(msg_out.dist.params[:p]) == 2) || error("Probit node only supports categorical messages with 2 categories")
    p = msg_out.dist.params[:p][1]

    return ruleEPProbitIn1BG(Message(Univariate, Bernoulli, p=p), msg_in1)
end


text:
"""
    Poll(p::DSProblem{T})::IterationOutcome where T

Generate points and call evaluate on them.
"""
function Poll(p::DSProblem{T})::IterationOutcome where T
    points = GeneratePollPoints(p, p.mesh)
    return EvaluatePoint!(p, points)
end

"""
    GeneratePollPoints(p::DSProblem{T}, ::AbstractMesh)::Vector{Vector{T}} where T

Generate a set of directions with the configured polling algorithm, then return
the set of points these directions give from the incumbent points.
"""
function GeneratePollPoints(p::DSProblem{T}, ::AbstractMesh
                           )::Vector{Vector{T}} where T
    points = []
    dirs = GenerateDirections(p)

    if !isnothing(p.x)
        append!(points, [p.x+(p.mesh.Δᵐ*p.meshscale.*d) for d in eachcol(dirs)])
    end
    if !isnothing(p.i)
        append!(points, [p.i+(p.mesh.Δᵐ*p.meshscale.*d) for d in eachcol(dirs)])
    end

    return points
end

GenerateDirections(p::DSProblem) = GenerateDirections(p, p.poll)

ext: jl


text:
using Pkg

Pkg.test(ARGS[1], coverage=true)

try
    import Coverage
catch
    exit()
end

coverage = Coverage.process_folder()

Coverage.LCOV.writefile("lcov.info", coverage)

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode


text:
using ACE1pack
using ArtifactUtils
using Pkg.Artifacts

function safe_add_artifact!(label, url, artifacts_toml)
    data_hash = artifact_hash(label, artifacts_toml)
    if data_hash == nothing || !artifact_exists(data_hash)
        add_artifact!(artifacts_toml, label, url, lazy=true, force=true)
    end
end

artifacts_toml = joinpath(pathof(ACE1pack)[1:end-16], "Artifacts.toml")

label = "TiAl_tiny_dataset"
url = "https://github.com/ACEsuit/ACEData/blob/master/trainingsets/TiAl_tiny.tar.gz?raw=true"
safe_add_artifact!(label, url, artifacts_toml)

label = "ACE1pack_test_files"
url = "https://github.com/ACEsuit/ACEData/blob/master/tests/ACE1pack_test_files.tar.gz?raw=true"
safe_add_artifact!(label, url, artifacts_toml)


text:
"""
# Dynamic programming

- Shortest paths in dags, revisited
- Longest increasing subsequences
    - Recursion? No, thanks.
    - Programming?
- Edit distance
    - A dynamic programming solution
    - The underlying dag
    - Common subproblems
    - Of mice and men
- Knapsack
    - Knapsack with repetition
    - Knapsack without repetition
    - Memoization
- Chain matrix multiplication
- Shortest paths
    - Shortest reliable paths
    - All-pairs shortest paths
    - The traveling salesman problem
- Independent sets in trees
    - On time and memory
"""
module Chapter6

export
    lcs

"""
    lcs(A, B)

Find the longest common subsequence of two sequences `A` and `B`.

# Examples
```jldoctest
julia> lcs("This is a test", "This is not a test")
"This is a test"
```
"""
function lcs(A, B)
    T = zeros(length(A) + 1, length(B) + 1)

    for i ∈ 1:(length(A) + 1), j ∈ 1:(length(B) + 1)
        if i == 1 || j == 1
            T[i, j] = 0
        elseif A[i - 1] == B[j - 1]
            T[i, j] = T[i - 1, j - 1] + 1
        else
            T[i, j] = max(T[i - 1, j], T[i, j - 1])
        end
    end

    return lcs(A, B, T)
end

function lcs(a, b, T)
    result = []
    i, j = size(T)
    while i > 1 && j > 1
        if a[i - 1] == b[j - 1]
            pushfirst!(result, a[i - 1])
            i -= 1
            j -= 1
        elseif T[i - 1, j] > T[i, j - 1]
            i -= 1
        else
            j -= 1
        end
    end
    return join(result)
end

end # module


text:
using ForwardDiff, DiffResults
using FiniteDiff
include("preconditioners.jl")

function standard_optimization!(parameters::Array{Float64, 1}, loss_function::Function, epochs::Integer, optim)
    result = DiffResults.GradientResult(parameters)
    # cfg = ForwardDiff.GradientConfig(loss_function, parameters, ForwardDiff.Chunk{length(parameters)}())
    history = zeros(epochs)
    for i = 1:epochs
        ForwardDiff.gradient!(result, loss_function, parameters)
        update!(optim, parameters, DiffResults.gradient(result))
        history[i] = DiffResults.value(result)
    end
    return history
end

function early_stopping_optimization!(parameters::Array{Float64, 1}, loss_function::Function, epochs::Integer, optim; check_time::Integer=10, tol::Real=1e-2)
    result = DiffResults.GradientResult(parameters)
    # cfg = ForwardDiff.GradientConfig(loss_function, parameters, ForwardDiff.Chunk{length(parameters)}())
    history = []
    j = 1
    running_average = []
    while j<=epochs
        ForwardDiff.gradient!(result, loss_function, parameters)
        update!(optim, parameters, DiffResults.gradient(result))
        push!(history, DiffResults.value(result))
        # println(DiffResults.gradient(result))
        if (j%check_time) == 0
            push!(running_average, sum(history[(end-check_time+1):end])/check_time)
            if length(running_average)>=2
                difference = abs(running_average[end] - running_average[end-1])
                if difference<tol
                    break
                end
            end
        end
        j += 1
    end
    return history
end

function tracking_optimization!(parameters::Array{Float64, 1}, prec::preconditioner, params::optimization_parameters, loss_function::Function)
    result = similar(parameters)
    history = zeros(params.epochs)
    condition_number = zeros(params.epochs)
    for i = 1:params.epochs
        ForwardDiff.gradient!(result, loss_function, parameters)
        value = loss_function(parameters)
        update!(params.optim, parameters, result)
        history[i] = value
        # condition_number[i] = measure_s_condition_number(parameters, prec, params)
        condition_number[i] = measure_condition_number(parameters, prec)
    end
    return history, condition_number
end

function tracking_optimization!(parameters::Array{Float64, 1}, prec::T_preconditioner, params::optimization_parameters, loss_function::Function)
    result = similar(parameters)
    history = zeros(params.epochs)
    condition_number = zeros(params.epochs)
    for i = 1:params.epochs
        ForwardDiff.gradient!(result, loss_function, parameters)
        value = loss_function(parameters)
        update!(params.optim, parameters, result)
        history[i] = value
        # condition_number[i] = measure_s_condition_number(parameters, prec, params)
        condition_number[i] = measure_condition_number(parameters, prec)
    end
    return history, condition_number
end

function spectral_optimization!(parameters::Array{Float64, 1}, loss_function::Function, ρ::Function, epochs::Integer, optim)
    result = DiffResults.GradientResult(parameters)
    result_2 = DiffResults.GradientResult(parameters)
    grads = Array{Float64, 1}(undef, length(parameters))
    history = zeros(epochs)
    for i = 1:epochs
        ForwardDiff.gradient!(result_2, ρ, parameters)
        parameters[end] = 1/(DiffResults.value(result_2))
        ForwardDiff.gradient!(result, loss_function, parameters)
        grads .= DiffResults.gradient(result) - parameters[end]^2 * DiffResults.gradient(result)[end] * DiffResults.gradient(result_2)
        update!(optim, parameters, grads)
        history[i] = DiffResults.value(result)
    end
    return history
end

function tracking_spectral_optimization!(parameters::Array{Float64, 1}, prec::T_preconditioner, params::optimization_parameters, loss_function::Function, ρ::Function, epochs::Integer, optim)
    result = DiffResults.GradientResult(parameters)
    result_2 = DiffResults.GradientResult(parameters)
    grads = Array{Float64, 1}(undef, length(parameters))
    condition_number = zeros(epochs)
    history = zeros(epochs)
    for i = 1:epochs
        ForwardDiff.gradient!(result_2, ρ, parameters)
        parameters[end] = 1/(DiffResults.value(result_2))
        ForwardDiff.gradient!(result, loss_function, parameters)
        grads .= DiffResults.gradient(result) - parameters[end]^2 * DiffResults.gradient(result)[end] * DiffResults.gradient(result_2)
        update!(optim, parameters, grads)
        history[i] = DiffResults.value(result)
        # condition_number[i] = measure_s_condition_number(parameters, prec, params)
        condition_number[i] = measure_condition_number(parameters, prec)
    end
    return history, condition_number
end

function tracking_spectral_optimization!(parameters::Array{Float64, 1}, prec::preconditioner, params::optimization_parameters, loss_function::Function, ρ::Function, epochs::Integer, optim)
    result = DiffResults.GradientResult(parameters)
    result_2 = DiffResults.GradientResult(parameters)
    grads = Array{Float64, 1}(undef, length(parameters))
    condition_number = zeros(epochs)
    history = zeros(epochs)
    for i = 1:epochs
        ForwardDiff.gradient!(result_2, ρ, parameters)
        parameters[end] = 1/(DiffResults.value(result_2))
        ForwardDiff.gradient!(result, loss_function, parameters)
        grads .= DiffResults.gradient(result) - parameters[end]^2 * DiffResults.gradient(result)[end] * DiffResults.gradient(result_2)
        update!(optim, parameters, grads)
        history[i] = DiffResults.value(result)
        # condition_number[i] = measure_s_condition_number(parameters, prec, params)
        condition_number[i] = measure_condition_number(parameters, prec)
    end
    return history, condition_number
end

max_stars_repo_name: VLSF/neuralBPX

max_issues_repo_name: VLSF/neuralBPX

max_forks_repo_name: VLSF/neuralBPX


text:
using TestImages
using Base.Test

img = testimage("cameraman")


text:
"""
SkewNormal(ξ, ω, α)
    The *skew normal distribution* is a continuous probability distribution
    that generalises the normal distribution to allow for non-zero skewness.
External links
* [Skew normal distribution on Wikipedia](https://en.wikipedia.org/wiki/Skew_normal_distribution)
* [Discourse](https://discourse.julialang.org/t/skew-normal-distribution/21549/7)
* [SkewDist.jl](https://github.com/STOR-i/SkewDist.jl)
"""
struct SkewNormal{T<:Real} <: ContinuousUnivariateDistribution
    ξ::T
    ω::T
    α::T
    SkewNormal{T}(ξ::T, ω::T, α::T) where {T} = new{T}(ξ, ω, α)
end

function SkewNormal(ξ::T, ω::T, α::T; check_args::Bool=true) where {T <: Real}
    @check_args SkewNormal (ω, ω > zero(ω))
    return SkewNormal{T}(ξ, ω, α)
end

SkewNormal(ξ::Real, ω::Real, α::Real; check_args::Bool=true) = SkewNormal(promote(ξ, ω, α)...; check_args=check_args)
SkewNormal(ξ::Integer, ω::Integer, α::Integer; check_args::Bool=true) = SkewNormal(float(ξ), float(ω), float(α); check_args=check_args)
SkewNormal(α::Real=0.0) = SkewNormal(zero(α), one(α), α; check_args=false)

@distr_support SkewNormal -Inf Inf

#### Conversions
convert(::Type{SkewNormal{T}}, ξ::S, ω::S, α::S) where {T <: Real, S <: Real} = SkewNormal(T(ξ), T(ω), T(α))
convert(::Type{SkewNormal{T}}, d::SkewNormal{S}) where {T <: Real, S <: Real} = SkewNormal(T(d.ξ), T(d.ω), T(d.α), check_args=false)

#### Parameters
params(d::SkewNormal) = (d.ξ, d.ω, d.α)
@inline partype(d::SkewNormal{T}) where {T<:Real} = T

#### Statistics
delta(d::SkewNormal) = d.α / √(1 + d.α^2)
mean_z(d::SkewNormal) = √(2/π) * delta(d)
std_z(d::SkewNormal) = 1 - (2/π) * delta(d)^2

mean(d::SkewNormal) = d.ξ + d.ω * mean_z(d)
var(d::SkewNormal) = abs2(d.ω) * (1 - mean_z(d)^2)
std(d::SkewNormal) = √var(d)
skewness(d::SkewNormal) = ((4 - π)/2) * (mean_z(d)^3/(1 - mean_z(d)^2)^(3/2))
kurtosis(d::SkewNormal) = 2 * (π-3) * ((delta(d) * sqrt(2/π))^4/(1-2 * (delta(d)^2)/π)^2)

# no analytic expression for max m_0(d) but accurate numerical approximation
m_0(d::SkewNormal) = mean_z(d) - (skewness(d) * std_z(d))/2 - (sign(d.α)/2) * exp(-2π/abs(d.α))
mode(d::SkewNormal) = d.ξ + d.ω * m_0(d)

#### Evalution
pdf(d::SkewNormal, x::Real) = (2/d.ω) * normpdf((x-d.ξ)/d.ω) * normcdf(d.α * (x-d.ξ)/d.ω)
logpdf(d::SkewNormal, x::Real) = log(2) - log(d.ω) + normlogpdf((x-d.ξ) / d.ω) + normlogcdf(d.α * (x-d.ξ) / d.ω)
#cdf requires Owen's T function.
#cdf/quantile etc

mgf(d::SkewNormal, t::Real) = 2 * exp(d.ξ * t + (d.ω^2 * t^2)/2 ) * normcdf(d.ω * delta(d) * t)

cf(d::SkewNormal, t::Real) = exp(im * t * d.ξ - (d.ω^2 * t^2)/2) * (1 + im * erfi((d.ω * delta(d) * t)/(sqrt(2))) )

#### Sampling
function rand(rng::AbstractRNG, d::SkewNormal)
    u0 = randn(rng)
    v = randn(rng)
    δ = delta(d)
    u1 = δ * u0 + √(1-δ^2) * v
    return d.ξ + d.ω * sign(u0) * u1
end

## Fitting  # to be added see: https://github.com/STOR-i/SkewDist.jl/issues/3


text:
function test_vector_value(l::MarginLoss, t, y)
    ref = [value(l,t[i],y[i]) for i in CartesianIndices(size(y))]
    yt = t .* y
    buf = fill!(similar(ref), 0)
    @test @inferred(value!(buf, l, t, y)) == ref
    @test buf == ref
    buf1 = fill!(similar(ref), 0)
    @test @inferred(value!(buf1, l, yt)) == ref
    @test buf1 == ref
    buf2 = fill!(similar(ref), 0)
    @test @inferred(value!(buf2, l, t, y, AggMode.None())) == ref
    @test buf2 == ref
    buf3 = fill!(similar(ref), 0)
    @test @inferred(value!(buf3, l, yt, AggMode.None())) == ref
    @test buf3 == ref
    @test @inferred(value(l, t, y, AggMode.None())) == ref
    @test @inferred(value(l, t, y)) == ref
    @test @inferred(value(l, yt, AggMode.None())) == ref
    @test @inferred(value(l, yt)) == ref
    @test value.(Ref(l), t, y) == ref
    @test value.(Ref(l), yt) == ref
    # Sum
    s = sum(ref)
    @test @inferred(value(l, t, y, AggMode.Sum())) ≈ s
    @test @inferred(value(l, yt, AggMode.Sum())) ≈ s
    # Mean
    m = mean(ref)
    @test @inferred(value(l, t, y, AggMode.Mean())) ≈ m
    @test @inferred(value(l, yt, AggMode.Mean())) ≈ m
    # Obs specific
    n = size(t, 1)
    k = size(t, ndims(t))
    ## Weighted Mean
    @test_throws DimensionMismatch value(l, t, y, AggMode.WeightedMean(ones(n-1)), ObsDim.First())
    @test_throws DimensionMismatch value(l, yt, AggMode.WeightedMean(ones(n-1)), ObsDim.First())
    @test @inferred(value(l, t, y, AggMode.WeightedMean(ones(n)), ObsDim.First())) ≈ m
    @test @inferred(value(l, yt, AggMode.WeightedMean(ones(n)), ObsDim.First())) ≈ m
    @test @inferred(value(l, t, y, AggMode.WeightedMean(ones(k)), ObsDim.Last())) ≈ m
    @test @inferred(value(l, yt, AggMode.WeightedMean(ones(k)), ObsDim.Last())) ≈ m
    @test @inferred(value(l, t, y, AggMode.WeightedMean(ones(k)))) ≈ m
    @test @inferred(value(l, yt, AggMode.WeightedMean(ones(k)))) ≈ m
    ## Weighted Sum
    @test_throws DimensionMismatch value(l, t, y, AggMode.WeightedSum(ones(n-1)), ObsDim.First())
    @test_throws DimensionMismatch value(l, yt, AggMode.WeightedSum(ones(n-1)), ObsDim.First())
    @test @inferred(value(l, t, y, AggMode.WeightedSum(ones(n)), ObsDim.First())) ≈ s
    @test @inferred(value(l, yt, AggMode.WeightedSum(ones(n)), ObsDim.First())) ≈ s
    @test @inferred(value(l, t, y, AggMode.WeightedSum(ones(k)), ObsDim.Last())) ≈ s
    @test @inferred(value(l, yt, AggMode.WeightedSum(ones(k)), ObsDim.Last())) ≈ s
    @test @inferred(value(l, t, y, AggMode.WeightedSum(ones(k)))) ≈ s
    @test @inferred(value(l, yt, AggMode.WeightedSum(ones(k)))) ≈ s
    if typeof(t) <: AbstractVector
        @test_throws ArgumentError value(l, t, y, AggMode.Sum(), ObsDim.First())
        @test_throws ArgumentError value(l, t, y, AggMode.Mean(), ObsDim.First())
    else
        # Sum per obs
        sv = vec(sum(ref, dims=1:(ndims(ref)-1)))
        @test @inferred(value(l, t, y, AggMode.Sum(), ObsDim.Last())) ≈ sv
        @test @inferred(value(l, yt, AggMode.Sum(), ObsDim.Last())) ≈ sv
        buffer1 = fill!(similar(sv), 0)
        @test @inferred(value!(buffer1, l, t, y, AggMode.Sum(), ObsDim.Last())) ≈ sv
        @test buffer1 ≈ sv
        buffer2 = fill!(similar(sv), 0)
        @test @inferred(value!(buffer2, l, yt, AggMode.Sum(), ObsDim.Last())) ≈ sv
        @test buffer2 ≈ sv
        # Weighted sum compare
        @test @inferred(value(l, t, y, AggMode.WeightedSum(1:k), ObsDim.Last())) ≈ sum(sv .* (1:k))
        @test @inferred(value(l, yt, AggMode.WeightedSum(1:k), ObsDim.Last())) ≈ sum(sv .* (1:k))
        @test round(@inferred(value(l, t, y, AggMode.WeightedSum(1:k,normalize=true), ObsDim.Last())), digits=3) ≈
            round(sum(sv .* ((1:k)/(sum(1:k)))), digits=3)
        @test round(@inferred(value(l, yt, AggMode.WeightedSum(1:k,normalize=true), ObsDim.Last())), digits=3) ≈
            round(sum(sv .* ((1:k)/(sum(1:k)))), digits=3)
        # Mean per obs
        mv = vec(mean(ref, dims=1:(ndims(ref)-1)))
        @test @inferred(value(l, t, y, AggMode.Mean(), ObsDim.Last())) ≈ mv
        @test @inferred(value(l, yt, AggMode.Mean(), ObsDim.Last())) ≈ mv
        buffer3 = fill!(similar(mv), 0)
        @test @inferred(value!(buffer3, l, t, y, AggMode.Mean(), ObsDim.Last())) ≈ mv
        @test buffer3 ≈ mv
        buffer4 = fill!(similar(mv), 0)
        @test @inferred(value!(buffer4, l, yt, AggMode.Mean(), ObsDim.Last())) ≈ mv
        @test buffer4 ≈ mv
        # Weighted mean compare
        @test @inferred(value(l, t, y, AggMode.WeightedMean(1:k,normalize=false), ObsDim.Last())) ≈ sum(mv .* (1:k))
        @test @inferred(value(l, yt, AggMode.WeightedMean(1:k,normalize=false), ObsDim.Last())) ≈ sum(mv .* (1:k))
    end
end

function test_vector_value(l::DistanceLoss, t, y)
    ref = [value(l,t[i],y[i]) for i in CartesianIndices(size(y))]
    yt = y .- t
    buf = fill!(similar(ref), 0)
    @test @inferred(value!(buf, l, t, y)) == ref
    @test buf == ref
    buf1 = fill!(similar(ref), 0)
    @test @inferred(value!(buf1, l, yt)) == ref
    @test buf1 == ref
    buf2 = fill!(similar(ref), 0)
    @test @inferred(value!(buf2, l, t, y, AggMode.None())) == ref
    @test buf2 == ref
    buf3 = fill!(similar(ref), 0)
    @test @inferred(value!(buf3, l, yt, AggMode.None())) == ref
    @test buf3 == ref
    @test @inferred(value(l, t, y, AggMode.None())) == ref
    @test @inferred(value(l, t, y)) == ref
    @test @inferred(value(l, yt, AggMode.None())) == ref
    @test @inferred(value(l, yt)) == ref
    @test value.(Ref(l), t, y) == ref
    @test value.(Ref(l), yt) == ref
    # Sum
    s = sum(ref)
    @test @inferred(value(l, t, y, AggMode.Sum())) ≈ s
    @test @inferred(value(l, yt, AggMode.Sum())) ≈ s
    # Mean
    m = mean(ref)
    @test round(@inferred(value(l, t, y, AggMode.Mean())), digits=5) ≈ round(m, digits=5)
    @test round(@inferred(value(l, yt, AggMode.Mean())), digits=5) ≈ round(m, digits=5)
    # Obs specific
    n = size(t, 1)
    k = size(t, ndims(t))
    ## Weighted Mean
    @test_throws DimensionMismatch value(l, t, y, AggMode.WeightedMean(ones(n-1)), ObsDim.First())
    @test_throws DimensionMismatch value(l, yt, AggMode.WeightedMean(ones(n-1)), ObsDim.First())
    @test @inferred(value(l, t, y, AggMode.WeightedMean(ones(n)), ObsDim.First())) ≈ m
    @test @inferred(value(l, yt, AggMode.WeightedMean(ones(n)), ObsDim.First())) ≈ m
    @test @inferred(value(l, t, y, AggMode.WeightedMean(ones(k)), ObsDim.Last())) ≈ m
    @test @inferred(value(l, yt, AggMode.WeightedMean(ones(k)), ObsDim.Last())) ≈ m
    ## Weighted Sum
    @test_throws DimensionMismatch value(l, t, y, AggMode.WeightedSum(ones(n-1)), ObsDim.First())
    @test_throws DimensionMismatch value(l, yt, AggMode.WeightedSum(ones(n-1)), ObsDim.First())
    @test @inferred(value(l, t, y, AggMode.WeightedSum(ones(n),normalize=false), ObsDim.First())) ≈ s
    @test @inferred(value(l, yt, AggMode.WeightedSum(ones(n),normalize=false), ObsDim.First())) ≈ s
    @test @inferred(value(l, t, y, AggMode.WeightedSum(ones(k),normalize=false), ObsDim.Last())) ≈ s
    @test @inferred(value(l, yt, AggMode.WeightedSum(ones(k),normalize=false), ObsDim.Last())) ≈ s
    if typeof(t) <: AbstractVector
        @test_throws ArgumentError value(l, t, y, AggMode.Sum(), ObsDim.First())
        @test_throws ArgumentError value(l, t, y, AggMode.Mean(), ObsDim.First())
    else
        # Sum per obs
        sv = vec(sum(ref, dims=1:(ndims(ref)-1)))
        @test @inferred(value(l, t, y, AggMode.Sum(), ObsDim.Last())) ≈ sv
        @test @inferred(value(l, yt, AggMode.Sum(), ObsDim.Last())) ≈ sv
        buffer1 = fill!(similar(sv), 0)
        @test @inferred(value!(buffer1, l, t, y, AggMode.Sum(), ObsDim.Last())) ≈ sv
        @test buffer1 ≈ sv
        buffer2 = fill!(similar(sv), 0)
        @test @inferred(value!(buffer2, l, yt, AggMode.Sum(), ObsDim.Last())) ≈ sv
        @test buffer2 ≈ sv
        # Weighted sum compare
        @test @inferred(value(l, t, y, AggMode.WeightedSum(1:k,normalize=false), ObsDim.Last())) ≈ sum(sv .* (1:k))
        @test @inferred(value(l, yt, AggMode.WeightedSum(1:k,normalize=false), ObsDim.Last())) ≈ sum(sv .* (1:k))
        # Mean per obs
        mv = vec(mean(ref, dims=1:(ndims(ref)-1)))
        @test @inferred(value(l, t, y, AggMode.Mean(), ObsDim.Last())) ≈ mv
        @test @inferred(value(l, yt, AggMode.Mean(), ObsDim.Last())) ≈ mv
        buffer3 = fill!(copy(mv), 0)
        @test @inferred(value!(buffer3, l, t, y, AggMode.Mean(), ObsDim.Last())) ≈ mv
        @test buffer3 ≈ mv
        buffer4 = fill!(copy(mv), 0)
        @test @inferred(value!(buffer4, l, yt, AggMode.Mean(), ObsDim.Last())) ≈ mv
        @test buffer4 ≈ mv
        # Weighted mean compare
        @test @inferred(value(l, t, y, AggMode.WeightedMean(1:k,normalize=false), ObsDim.Last())) ≈ sum(mv .* (1:k))
        @test @inferred(value(l, yt, AggMode.WeightedMean(1:k,normalize=false), ObsDim.Last())) ≈ sum(mv .* (1:k))
    end
end

function test_vector_deriv(l::MarginLoss, t, y)
    ref = [LossFunctions.deriv(l,t[i],y[i]) for i in CartesianIndices(size(y))]
    yt = t .* y
    @test @inferred(LossFunctions.deriv(l, t, y, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv(l, t, y)) == ref
    @test t .* @inferred(LossFunctions.deriv(l, yt, AggMode.None())) == ref
    @test t .* @inferred(LossFunctions.deriv(l, yt)) == ref
    @test LossFunctions.deriv.(Ref(l), t, y) == ref
    @test t .* LossFunctions.deriv.(Ref(l), yt) == ref
end

function test_vector_deriv(l::DistanceLoss, t, y)
    ref = [LossFunctions.deriv(l,t[i],y[i]) for i in CartesianIndices(size(y))]
    yt = y .- t
    @test @inferred(LossFunctions.deriv(l, t, y, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv(l, t, y)) == ref
    @test @inferred(LossFunctions.deriv(l, yt, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv(l, yt)) == ref
    @test LossFunctions.deriv.(Ref(l), t, y) == ref
    @test LossFunctions.deriv.(Ref(l), yt) == ref
end

function test_vector_deriv2(l::MarginLoss, t, y)
    ref = [LossFunctions.deriv2(l,t[i],y[i]) for i in CartesianIndices(size(y))]
    yt = t .* y
    @test @inferred(LossFunctions.deriv2(l, t, y, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv2(l, t, y)) == ref
    @test @inferred(LossFunctions.deriv2(l, yt, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv2(l, yt)) == ref
    @test LossFunctions.deriv2.(Ref(l), t, y) == ref
    @test LossFunctions.deriv2.(Ref(l), yt) == ref
end

function test_vector_deriv2(l::DistanceLoss, t, y)
    ref = [LossFunctions.deriv2(l,t[i],y[i]) for i in CartesianIndices(size(y))]
    yt = y .- t
    @test @inferred(LossFunctions.deriv2(l, t, y, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv2(l, t, y)) == ref
    @test @inferred(LossFunctions.deriv2(l, yt, AggMode.None())) == ref
    @test @inferred(LossFunctions.deriv2(l, yt)) == ref
    @test LossFunctions.deriv2.(Ref(l), t, y) == ref
    @test LossFunctions.deriv2.(Ref(l), yt) == ref
end

@testset "Vectorized API" begin
    for T in (Float32, Float64)
        for O in (Float32, Float64)
            @testset "Margin-based $T -> $O" begin
                for (targets,outputs) in (
                        (rand(T[-1, 1], 4), (rand(O, 4) .- O(.5)) .* O(20)),
                        (rand(T[-1, 1], 4, 5), (rand(O, 4, 5) .- O(.5)) .* O(20)),
                    )
                    for loss in (LogitMarginLoss(),ModifiedHuberLoss(),
                                 L1HingeLoss(),SigmoidLoss())
                        @testset "$(loss): " begin
                            test_vector_value(loss, targets, outputs)
                            test_vector_deriv(loss, targets, outputs)
                            test_vector_deriv2(loss, targets, outputs)
                        end
                    end
                end
            end
            println("<HEARTBEAT>")
            @testset "Distance-based $T -> $O" begin
                for (targets,outputs) in (
                        ((rand(T, 4) .- T(.5)) .* T(20), (rand(O, 4) .- O(.5)) .* O(20)),
                        ((rand(T, 4, 5) .- T(.5)) .* T(20), (rand(O, 4, 5) .- O(.5)) .* O(20)),
                    )
                    for loss in (QuantileLoss(0.75),L2DistLoss(),
                                 EpsilonInsLoss(1))
                        @testset "$(loss): " begin
                            test_vector_value(loss, targets, outputs)
                            test_vector_deriv(loss, targets, outputs)
                            test_vector_deriv2(loss, targets, outputs)
                        end
                    end
                end
            end
            println("<HEARTBEAT>")
        end
    end
end

@testset "Broadcasting higher-order arrays" begin
    for f in (value,deriv,deriv2)
        @testset "$f" begin

            m,n,k = 100,99,98
            loss = L2DistLoss()

            targ1 = randn(m)
            targ2 = repeat(targ1,outer=(1,n))
            targ3 = repeat(targ1,outer=(1,n,k))

            out1 = randn(m,n)
            out2 = randn(m,n,k)

            for avg in (AggMode.None(),AggMode.Mean(),AggMode.Sum())
                @test f(loss,targ1,out1,avg) ≈ f(loss,targ2,out1,avg)
                @test f(loss,targ1,out2,avg) ≈ f(loss,targ2,out2,avg)
                @test f(loss,targ1,out2,avg) ≈ f(loss,targ3,out2,avg)
                @test f(loss,targ2,out2,avg) ≈ f(loss,targ3,out2,avg)
                @test f(loss,out1,targ1,avg) ≈ f(loss,out1,targ2,avg)
                @test f(loss,out2,targ1,avg) ≈ f(loss,out2,targ2,avg)
                @test f(loss,out2,targ1,avg) ≈ f(loss,out2,targ3,avg)
                @test f(loss,out2,targ2,avg) ≈ f(loss,out2,targ3,avg)
            end
        end
    end
end


text:
using ImageFiltering, ImageCore, OffsetArrays, Colors, FFTViews, ColorVectorSpace, ComputationalResources
using Base.Test

@testset "FIR/FFT" begin
    f32type(img) = f32type(eltype(img))
    f32type{C<:Colorant}(::Type{C}) = base_colorant_type(C){Float32}
    f32type{T<:Number}(::Type{T}) = Float32
    zerona!(a) = (a[isnan.(a)] = zero(eltype(a)); a)
    zerona!(a, nanflag) = (a[nanflag] = zero(eltype(a)); a)
    ## Images for which the boundary conditions will be irrelevant
    imgf = zeros(5, 7); imgf[3,4] = 1
    imgi = zeros(Int, 5, 7); imgi[3,4] = 1
    imgg = fill(Gray(0), 5, 7); imgg[3,4] = 1
    imgc = fill(RGB(0,0,0), 5, 7); imgc[3,4] = RGB(1,0,0)
    # Dense inseparable kernel
    kern = [0.1 0.2; 0.4 0.5]
    kernel = OffsetArray(kern, -1:0, 1:2)
    for img in (imgf, imgi, imgg, imgc)
        targetimg = zeros(typeof(img[1]*kern[1]), size(img))
        targetimg[3:4,2:3] = rot180(kern)*img[3,4]
        @test @inferred(imfilter(img, kernel)) ≈ targetimg
        @test @inferred(imfilter(img, (kernel,))) ≈ targetimg
        @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32.(targetimg)
        ret = imfilter(img, kernel)
        fill!(ret, zero(eltype(ret)))
        @test @inferred(imfilter!(ret, img, kernel)) ≈ targetimg
        fill!(ret, zero(eltype(ret)))
        @test @inferred(imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel)) ≈ targetimg
        for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(img))))
            @test @inferred(imfilter(img, kernel, border)) ≈ targetimg
            @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32.(targetimg)
            fill!(ret, zero(eltype(ret)))
            @test @inferred(imfilter!(ret, img, kernel, border)) ≈ targetimg
            for alg in (Algorithm.FIR(), Algorithm.FFT())
                @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg
                @test @inferred(imfilter(img, (kernel,), border, alg)) ≈ targetimg
                @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32.(targetimg)
                fill!(ret, zero(eltype(ret)))
                @test @inferred(imfilter!(CPU1(alg), ret, img, kernel, border)) ≈ targetimg
            end
            @test_throws MethodError imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel, border, Algorithm.FFT())
        end
        targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5)
        @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner
        @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32.(targetimg_inner)
        for alg in (Algorithm.FIR(), Algorithm.FFT())
            @test @inferred(imfilter(img, kernel, Inner(), alg)) ≈ targetimg_inner
            @test @inferred(imfilter(f32type(img), img, kernel, Inner(), alg)) ≈ float32.(targetimg_inner)
            @test @inferred(imfilter(CPU1(alg), img, kernel, Inner())) ≈ targetimg_inner
        end
    end
    # Factored kernel
    kernel = (OffsetArray([0.2,0.8], -1:0), OffsetArray([0.3 0.6], 0:0, 1:2))
    kern = parent(kernel[1]).*parent(kernel[2])
    for img in (imgf, imgi, imgg, imgc)
        targetimg = zeros(typeof(img[1]*kern[1]), size(img))
        targetimg[3:4,2:3] = rot180(kern)*img[3,4]
        ret = similar(targetimg)
        @test @inferred(imfilter(img, kernel)) ≈ targetimg
        @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32.(targetimg)
        for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(img))))
            @test @inferred(imfilter(img, kernel, border)) ≈ targetimg
            @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32.(targetimg)
            for alg in (Algorithm.FIR(), Algorithm.FFT())
                @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg
                @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32.(targetimg)
            end
            @test_throws MethodError imfilter!(CPU1(Algorithm.FIR()), ret, img, kernel, border, Algorithm.FFT())
        end
        targetimg_inner = OffsetArray(targetimg[2:end, 1:end-2], 2:5, 1:5)
        @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner
        @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32.(targetimg_inner)
        for alg in (Algorithm.FIR(), Algorithm.FFT())
            @test @inferred(imfilter(img, kernel, Inner(), alg)) ≈ targetimg_inner
            @test @inferred(imfilter(f32type(img), img, kernel, Inner(), alg)) ≈ float32.(targetimg_inner)
        end
    end
    # Rational filter coefficients
    kernel = (centered([1//3, 1//3, 1//3]), centered([1//3, 1//3, 1//3]'))
    kern = kernel[1].*kernel[2]
    for img in (imgf, imgi, imgg, imgc)
        targetimg = zeros(typeof(img[1]*kern[1]), size(img))
        targetimg[2:4,3:5] = img[3,4]*(1//9)
        @test @inferred(imfilter(img, kernel)) ≈ targetimg
        @test @inferred(imfilter(f32type(img), img, kernel)) ≈ float32.(targetimg)
        for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(img))))
            @test @inferred(imfilter(img, kernel, border)) ≈ targetimg
            @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32.(targetimg)
            for alg in (Algorithm.FIR(), Algorithm.FFT())
                if alg == Algorithm.FFT() && eltype(img) == Int
                    @test @inferred(imfilter(Float64, img, kernel, border, alg)) ≈ targetimg
                else
                    @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg
                end
                @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32.(targetimg)
            end
        end
        targetimg_inner = OffsetArray(targetimg[2:end-1, 2:end-1], 2:4, 2:6)
        @test @inferred(imfilter(img, kernel, Inner())) ≈ targetimg_inner
        @test @inferred(imfilter(f32type(img), img, kernel, Inner())) ≈ float32.(targetimg_inner)
        for alg in (Algorithm.FIR(), Algorithm.FFT())
            if alg == Algorithm.FFT() && eltype(img) == Int
                @test @inferred(imfilter(Float64, img, kernel, Inner(), alg)) ≈ targetimg_inner
            else
                @test @inferred(imfilter(img, kernel, Inner(), alg)) ≈ targetimg_inner
            end
            @test @inferred(imfilter(f32type(img), img, kernel, Inner(), alg)) ≈ float32.(targetimg_inner)
        end
    end
    ## Images for which the boundary conditions matter
    imgf = zeros(5, 7); imgf[1,2] = 1
    imgi = zeros(Int, 5, 7); imgi[1,2] = 1
    imgg = fill(Gray(0), 5, 7); imgg[1,2] = 1
    imgc = fill(RGB(0,0,0), 5, 7); imgc[1,2] = RGB(1,0,0)
    # Dense kernel
    kern = [0.1 0.2; 0.4 0.5]
    kernel = OffsetArray(kern, -1:0, 1:2)
    function target1{T}(img::AbstractArray{T}, kern, border)
        if border ∈ ("replicate", "symmetric")
            ret = float64.(zero(img))
            x = img[1,2]
            ret[1,1] = (kern[1,1]+kern[2,1])*x
            ret[2,1] = kern[1,1]*x
            ret
        elseif border == "circular"
            a = FFTView(float64.(zero(img)))
            x = img[1,2]
            a[0:1, -1:0] = rot180(kern).*x
            parent(a)
        elseif border == "reflect" || isa(border, Fill)
            ret = float64.(zero(img))
            x = img[1,2]
            ret[1,1] = kern[2,1]*x
            ret[2,1] = kern[1,1]*x
            ret
        end
    end
    function target1{T}(img::AbstractArray{T}, kern, ::NA)
        ret = float64.(zero(img))
        x = img[1,2]
        ret[1,1] = kern[2,1]*x/(kern[2,1]+kern[2,2])
        ret[2,1] = kern[1,1]*x/sum(kern)
        ret[:,end] = nan(eltype(ret))  # for the last column, this kernel is entirely in the padding region
        ret
    end
    for img in (imgf, imgi, imgg, imgc)
        for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(img))))
            targetimg = target1(img, kern, border)
            @test @inferred(imfilter(img, kernel, border)) ≈ targetimg
            @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32.(targetimg)
            for alg in (Algorithm.FIR(), Algorithm.FFT())
                @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg
                @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32.(targetimg)
            end
        end
        border = NA()
        targetimg0 = target1(img, kern, border)
        nanflag = isnan.(targetimg0)
        targetimg = zerona!(copy(targetimg0))
        @test @inferred(zerona!(imfilter(img, kernel, border))) ≈ targetimg
        @test @inferred(zerona!(imfilter(f32type(img), img, kernel, border))) ≈ float32.(targetimg)
        for alg in (Algorithm.FIR(), Algorithm.FFT())
            @test @inferred(zerona!(imfilter(img, kernel, border, alg), nanflag)) ≈ targetimg
            @test @inferred(zerona!(imfilter(f32type(img), img, kernel, border, alg), nanflag)) ≈ float32.(targetimg)
        end
    end
    # Factored kernel
    kernel = (OffsetArray([0.2,0.8], -1:0), OffsetArray([0.3 0.6], 0:0, 1:2))
    kern = parent(kernel[1]).*parent(kernel[2])
    for img in (imgf, imgi, imgg, imgc)
        for border in ("replicate", "circular", "symmetric", "reflect", Fill(zero(eltype(img))))
            targetimg = target1(img, kern, border)
            @test @inferred(imfilter(img, kernel, border)) ≈ targetimg
            @test @inferred(imfilter(f32type(img), img, kernel, border)) ≈ float32.(targetimg)
            for alg in (Algorithm.FIR(), Algorithm.FFT())
                @test @inferred(imfilter(img, kernel, border, alg)) ≈ targetimg
                @test @inferred(imfilter(f32type(img), img, kernel, border, alg)) ≈ float32.(targetimg)
            end
        end
        border = NA()
        targetimg0 = target1(img, kern, border)
        nanflag = isnan.(targetimg0)
        targetimg = zerona!(copy(targetimg0))
        @test @inferred(zerona!(imfilter(img, kernel, border))) ≈ targetimg
        @test @inferred(zerona!(imfilter(f32type(img), img, kernel, border))) ≈ float32.(targetimg)
        for alg in (Algorithm.FIR(), Algorithm.FFT())
            @test @inferred(zerona!(imfilter(img, kernel, border, alg), nanflag)) ≈ targetimg
            @test @inferred(zerona!(imfilter(f32type(img), img, kernel, border, alg), nanflag)) ≈ float32.(targetimg)
        end
    end
    # filtering with a 0d kernel
    a = rand(3,3)
    # @test_approx_eq imfilter(a, reshape([2])) 2a
    # imfilter! as a copy
    ret = zeros(3, 3)
    b = @inferred(imfilter!(CPU1(Algorithm.FIR()), ret, a, (), NoPad()))
    @test b == a
    @test !(b === a)
end

nothing


text:
function show(io::IO, dp::DPresults)
    println("######################################")
    println("Number of clones: $(dp.nclones)")
    println("Frequency of clones: $(dp.clonefrequencies)")
    println("Weight of clones: $(dp.cloneweights)")
    println("#######################################")
end


text:
using SimulatedNeuralMoments
using Flux, Turing, MCMCChains, AdvancedMH
using StatsPlots, DelimitedFiles, LinearAlgebra
using BSON:@save
using BSON:@load

# the model-specific code
include("MNlib.jl")

function main()

# fill in the structure that defines the model
lb, ub = PriorSupport() # bounds of support
model = SNMmodel("Mixture of normals example", lb, ub, InSupport, PriorDraw, auxstat)

# train the net, and save it and the transformation info
transf = bijector(@Prior) # transforms draws from prior to draws from  ℛⁿ
transformed_prior = transformed(@Prior, transf) # the transformed prior
#nnmodel, nninfo = MakeNeuralMoments(model, transf)
#@save "neuralmodel.bson" nnmodel nninfo  # use this line to save the trained neural net
@load "neuralmodel.bson" nnmodel nninfo # use this to load a trained net

# draw a sample at the design parameters, or use an existing data set
y = MNmodel(TrueParameters()) # draw a sample of 500 obsns. at design parameters
#y = readdlm("svdata.txt") # load a data set
n = size(y,1)
p1 = plot(y)
p2 = density(y)
plot(p1, p2, layout=(2,1))
#savefig("data.png")

# define the neural moments using the real data
m = NeuralMoments(auxstat(y), nnmodel, nninfo)
# the raw NN parameter estimate
θhat = invlink(@Prior, m)

# setting for sampling
names = ["μ1", "μ2 ", "σ1", "σ2", "prob"]
S = 100
covreps = 1000
length = 1500
nchains = 4
burnin = 0
tuning = 1.0
# the covariance of the proposal (scaled by tuning)
junk, Σp = mΣ(θhat, covreps, model, nnmodel, nninfo)

@model function MSM(m, S, model)
    θt ~ transformed_prior
    if !InSupport(invlink(@Prior, θt))
        Turing.@addlogprob! -Inf
        return
    end
    # sample from the model, at the trial parameter value, and compute statistics
    mbar, Σ = mΣ(invlink(@Prior,θt), S, model, nnmodel, nninfo)
    m ~ MvNormal(mbar, Symmetric(Σ))
end

chain = sample(MSM(m, S, model),
    MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(size(m,1)), tuning*Σp))),
    MCMCThreads(), length, nchains;  init_params=Iterators.repeated(m), discard_initial=burnin)

# single thread
#=
chain = sample(MSM(m, S, model),
    MH(:θt => AdvancedMH.RandomWalkProposal(MvNormal(zeros(size(m,1)), tuning*Σp))),
    length*nchains; init_params=m, discard_initial=burnin)
=#
# transform back to original domain
chain = Array(chain)
acceptance = size(unique(chain[:,1]),1)[1] / size(chain,1)
println("acceptance rate: $acceptance")
for i = 1:size(chain,1)
    chain[i,:] = invlink(@Prior, chain[i,:])
end
chain = Chains(chain, names)
chain
end
chain = main()

display(chain)
display(plot(chain))
#savefig("chain.png")

max_stars_repo_name: mcreel/SimulatedNeuralMoments

max_issues_repo_name: mcreel/SimulatedNeuralMoments

max_forks_repo_name: mcreel/SimulatedNeuralMoments


text:
# -*- coding: utf-8 -*-

# cd ~/.julia/v0.4/DirectX
# julia test/dx9adl_run_demos.jl

include("dx9adl_demo_0002.jl")

println("test demo_0002")
d9 = DirectX.connect(80 * 16, 80 * 9)
Demo_0002.demo_0002(d9)
DirectX.close(d9)
println("ok")


text:
#  Copyright 2016-2019  Eric S. Tellez <eric.tellez@infotec.mx>
#
#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.

export Performance, PerformanceResult, probe
using Statistics: mean

mutable struct PerformanceResult
    precision::Float64
    recall::Float64
    macrof1::Float64
    seconds::Float64
    exhaustive_search_seconds::Float64
    evaluations::Float64
    distances_sum::Float64
    PerformanceResult() = new(0, 0, 0, 0, 0, 0)
end

mutable struct Performance{T}
    db::AbstractVector{T}
    queries::AbstractVector{T}
    results::Vector{Set{Int}}
    distances_sum::Vector{Float64}
    expected_k::Int
    queries_from_db::Bool
    exhaustive_search_seconds::Float64

    function Performance(dist::Function, db::AbstractVector{T}, queries::AbstractVector{T}; queries_from_db=false, expected_k=10, create_index=Sequential) where T
        results = Vector{Set{Int}}(undef, length(queries))
        distances_sum = Vector{Float64}(undef, length(queries))

        s = create_index(db)
        start = time()
        for i in 1:length(queries)
            res = search(s, dist, queries[i], KnnResult(expected_k))
            if queries_from_db
                popfirst!(res)
            end
            results[i] = Set(item.objID for item in res)
            distances_sum[i] = sum(item.dist for item in res)
        end

        elapsed = time() - start
        new{T}(db, queries, results, distances_sum, expected_k, queries_from_db, elapsed / length(queries))
    end

    function Performance(db::AbstractVector{T}, dist::Function; num_queries::Int=128, expected_k::Int=10) where T
        queries = rand(db, num_queries)
        Performance(dist, db, queries, queries_from_db=true, expected_k=expected_k)
    end
end

function probe(perf::Performance, index::Index, dist::Function; repeat::Int=1, aggregation=:mean, field=:seconds)
    if repeat == 1
        return _probe(perf, index, dist)
    end

    if aggregation == :mean
        p = _probe(perf, index, dist)
        for i in 2:repeat
            q = _probe(perf, index, dist)
            p.recall += q.recall
            p.seconds += q.seconds
            p.evaluations += q.evaluations
        end
        p.recall /= repeat
        p.seconds /= repeat
        p.evaluations /= repeat
        return p
    end

    M = [_probe(perf, index, dist) for i in 1:repeat]
    sort!(M, by=(x) -> getfield(x, field))

    if aggregation == :median
        return M[ceil(Int, length(M) / 2)]
    elseif aggregation == :min
        return M[1]
    elseif aggregation == :max
        return M[end]
    else
        error("Unknown aggregation strategy: $aggregation")
    end
end

function _probe(perf::Performance, index::Index, dist::Function)
    eval_counter = 0
    function dist_(a_, b_)
        eval_counter += 1
        dist(a_, b_)
    end

    p = PerformanceResult()
    m = length(perf.queries)
    p.evaluations = eval_counter
    p.seconds = 0.0
    p.distances_sum = 0.0

    for i in 1:m
        start = time()
        res = search(index, dist_, perf.queries[i], KnnResult(perf.expected_k))
        p.seconds += time() - start
        if perf.queries_from_db
            popfirst!(res)
        end
        base = perf.results[i]
        curr = Set(item.objID for item in res)
        for item in res
            p.distances_sum += item.dist
        end
        tp = intersect(base, curr)  # tn
        fp = length(setdiff(curr, tp)) # |fp| == |ft| when |curr| == |base|
        fn = length(setdiff(base, tp))
        _tp = length(tp)
        precision = _tp / (_tp + fp)
        recall = _tp / (_tp + fn)
        p.precision += precision
        p.recall += recall
        p.macrof1 += 2 * precision * recall / (precision + recall)
    end

    p.evaluations = (eval_counter - p.evaluations) / m
    p.seconds = p.seconds / m
    p.precision /= m
    p.recall /= m
    p.macrof1 /= m
    p.exhaustive_search_seconds = perf.exhaustive_search_seconds

    p
end


text:
"""
    far_field_drag(system)

Computes induced drag using the Trefftz plane (far field method).

Note that this function assumes that the circulation distribution has already
been computed and is present in `system`

# Arguments
 - `system`: Pre-allocated system properties
"""
function far_field_drag(system)

    # unpack system
    surfaces = system.surfaces
    trefftz = system.trefftz
    ref = system.reference[]
    fs = system.freestream[]
    symmetric = system.symmetric
    Γ = system.Γ

    # construct trefftz panels
    trefftz_panels!(trefftz, surfaces, fs, Γ)

    # perform far field analysis
    nsurf = length(surfaces)
    CD = zero(eltype(system))
    for i = 1:nsurf, j = 1:nsurf
        CD += far_field_drag(trefftz[i], trefftz[j], ref, symmetric[j])
    end

    return CD
end

# one surface on another surface
"""
    far_field_drag(receiving, sending, reference, symmetric)

Computes the induced drag on `receiving` from `sending` using the Trefftz
plane analysis.

# Arguments
 - `receiving`: Vector of receiving Trefftz panels (see [`TrefftzPanel`](@ref))
 - `sending`: Vector of sending Trefftz panels (see [`TrefftzPanel`](@ref))
 - `reference`: Reference parameters (see [`Reference`](@ref))
 - `symmetric`: Flag indicating whether a mirror image of the panels in `surface`,
    should be used when calculating induced velocities
"""
@inline function far_field_drag(receiving, sending, ref, symmetric)

    # get float type
    TF = promote_type(eltype(eltype(receiving)), eltype(eltype(sending)), eltype(ref))

    # get number of receiving and sending panels
    Nr = length(receiving)
    Ns = length(sending)

    # add up drag
    Di = zero(TF)
    for j = 1:Ns, i = 1:Nr
        Di += trefftz_panel_induced_drag(receiving[i], sending[j]; symmetric)
    end

    # apply symmetry
    if symmetric
        Di *= 2
    end

    # reference dynamic pressure
    q = 1/2*RHO*ref.V^2

    # normalize
    CDi = Di / (q*ref.S)

    return CDi
end


text:
################################################################################
#
#             NfOrd/NfOrdGen.jl : Generic orders in Number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export EquationOrder, Order

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::NfOrdGenSet)
  print(io, "Set of orders of the number field ")
  print(io, a.nf)
end

function show(io::IO, a::NfOrdGen)
  print(io, "Order of ")
  println(io, a.nf)
  print(io, "with Z-basis ")
  print(io, basis(a))
end

################################################################################
#
#  Construction
#
################################################################################

doc"""
    Order(B::Array{nf_elem, 1}, check::Bool = true) -> NfOrd

> Returns the order with $\mathbf Z$-basis $B$. If `check` is set, it is checked
> whether $B$ defines an order.
"""
function Order(a::Array{nf_elem, 1}, check::Bool = true)
  # We should check if it really is a basis and the elements are integral
  return NfOrdGen(nf_elem[ deepcopy(x) for x in a])
end

doc"""
    Order(K::AnticNumberField, A::FakeFmpqMat, check::Bool = true) -> NfOrd

> Returns the order which has basis matrix $A$ with respect to the power basis of $K$.
> If `check` is set, it is checked whether $A$ defines an order.
"""
function Order(K::AnticNumberField, a::FakeFmpqMat, check::Bool = true)
  # We should check if a has full rank and the elements are integral?
  return NfOrdGen(K, deepcopy(a))
end

doc"""
    EquationOrder(K::AnticNumberField) -> NfOrd

> Returns the equation of the number field $K$.
"""
function EquationOrder(K::AnticNumberField)
  z = NfOrdGen(K)
  z.isequationorder = true
  return z
end

################################################################################
#
#  Creation from fractional ideals
#
################################################################################

function NfOrdGen(a::NfOrdGenFracIdl)
  z = NfOrdGen(nf(order(a)), a.basis_mat*order(a).basis_mat)
  return z
end


text:
module complete

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
using DynamicGrids, Test
using DynamicGrids: SimData, WritableGridData, Extent,
    _getreadgrids, _getwritegrids, _readcell, _writecell!, grids, source, dest

rule = Cell{Tuple{:c,:a,:b},Tuple{:a,:c}}(identity)
init = (a=fill(1, 4, 4), b=fill(2, 4, 4), c=fill(3, 4, 4))
simdata = SimData(Extent(;init=init, tspan=1:1), Ruleset(rule))

@testset "_getreadgrids gets read grids for a Rule" begin
    rkeys, rgrids = _getreadgrids(rule, simdata)
    @test rkeys === (Val(:c), Val(:a), Val(:b))
    @test rgrids === (simdata[:c], simdata[:a], simdata[:b])
end

@testset "_getwritegrids gets write grids for a Rule" begin
    wkeys, wgrids = _getwritegrids(rule, simdata)
    @test wkeys === (Val(:a), Val(:c))
    @test wgrids === map(WritableGridData, (simdata[:a], simdata[:c]))
end

@testset "_readcell read from specified grids" begin
    @test _readcell(simdata, Val(:a), 1, 1) == 1
    @test _readcell(simdata, Val(:b), 1, 1) == 2
    @test _readcell(simdata, Val(:c), 1, 1) == 3
    @test _readcell(simdata, (Val(:c), Val(:a), Val(:b)), 1, 1) == (c=3, a=1, b=2)
    @test _readcell(simdata, (Val(:a), Val(:c)), 1, 1) == (a=1, c=3)
end

@testset "_writecell writes to source for CellRule" begin
    simdata = SimData(Extent(;init=init, tspan=1:1), Ruleset(rule))
    _writecell!(simdata, Val(CellRule), (Val(:c), Val(:a), Val(:b)), (8, 6, 7), 1, 2)
    @test map(g -> source(g)[1, 2], grids(simdata)) == (a=6, b=7, c=8)
    @test map(g -> dest(g)[1, 2], grids(simdata)) == (a=1, b=2, c=3)
    simdata = SimData(Extent(;init=init, tspan=1:1), Ruleset(rule))
    _writecell!(simdata, Val(CellRule), Val(:c), 99, 4, 3)
    @test source(simdata[:c])[4, 3] == 99
    @test dest(simdata[:c])[4, 3] == 3
end

@testset "_writecell writes to dest for other Rules" begin
    simdata = SimData(Extent(;init=init, tspan=1:1), Ruleset(rule))
    _writecell!(simdata, Val(Rule), (Val(:b), Val(:a)), (11, 10), 4, 4)
    @test map(g -> source(g)[4, 4], grids(simdata)) == (a=1, b=2, c=3)
    @test map(g -> dest(g)[4, 4], grids(simdata)) == (a=10, b=11, c=3)
    simdata = SimData(Extent(;init=init, tspan=1:1), Ruleset(rule))
    _writecell!(simdata, Val(Rule), Val(:c), 99, 4, 3)
    @test source(simdata[:c])[4, 3] == 3
    @test dest(simdata[:c])[4, 3] == 99
end

nothing


text:
using Test, SimpleHypergraphs
import LightGraphs

h1 = Hypergraph{Float64}(5,4)
h1[1:3,1] .= 1.5
h1[3,4] = 2.5
h1[2,3] = 3.5
h1[4,3:4] .= 4.5
h1[5,4] = 5.5
h1[5,2] = 6.5

@testset "SimpleHypergraphs Hypergraph    " begin

    h = hg_load("data/test1.hgf", Int)
    @test size(h) == (4, 4)
    m = Matrix(h)
    @test m == h
    @test h == [1       nothing 4       nothing
                2       3       nothing nothing
                nothing nothing 5       nothing
                nothing nothing 6       nothing]
    mktemp("data") do path, _
        println(path)
        hg_save(path, h)

        loaded_hg = replace(read(path, String), r"\n*$" => "")

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test1.hgf", String)) #no comments

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test_singlelinecomment.hgf", String)) #single line comment

        @test loaded_hg ==
            reduce(replace,
                ["\r\n"=>"\n",
                r"^\"\"\"(?s).*\"\"\"\n"=>"", #remove initial comments
                r"\n*$"=>""], #remove final \n*
                init=read("data/test_multiplelinescomment.hgf", String)) #multiple lines comment
    end

    @test_throws ArgumentError hg_load("data/test_malformedcomment.hgf", Int)

    h2 = Hypergraph{Float64}(0,0)
    for i in 1:4 add_vertex!(h2) end
    add_hyperedge!(h2;vertices=Dict(1:3 .=> 1.5))
    add_hyperedge!(h2)
    add_vertex!(h2;hyperedges=Dict(2=>6.5))
    add_hyperedge!(h2;vertices=Dict(2 => 3.5, 4 => 4.5))
    add_hyperedge!(h2;vertices=Dict(3:5 .=> (2.5,4.5,5.5)))
    @test h1 == h2
    m = Matrix(h1)
    @test  m == Matrix(h2)
    @test h1 == Hypergraph(m)
    @test h1 == Hypergraph{Nothing}(m)
    @test h1 == Hypergraph{Nothing, Nothing}(m)
    @test getindex(h1,3,1) == 1.5

    h3 = Hypergraph{Float64,String,Nothing}(1,1)
    @test add_vertex!(h3,vertex_meta="test") == 2
    @test set_vertex_meta!(h3,"t",1) == ["t","test"]
    @test get_vertex_meta(h3,2) == "test"
    @test get_hyperedge_meta(h3,1) == nothing
    @test_throws BoundsError get_hyperedge_meta(h3,2)

    h4 = Hypergraph{Float64,Nothing,String}(1,1)
    @test add_hyperedge!(h4,hyperedge_meta="test") == 2
    @test set_hyperedge_meta!(h4,"t",1) == ["t","test"]
    @test get_hyperedge_meta(h4,2) == "test"
    @test get_vertex_meta(h4,1) == nothing
    @test_throws BoundsError get_vertex_meta(h4,2)
end;

@testset "SimpleHypergraphs BipartiteView " begin
    b = BipartiteView(h1)

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(b))) == 18

    @test sort(collect(LightGraphs.outneighbors(b,5))) == [7,9]
    @test sort(collect(LightGraphs.outneighbors(b,1))) == [6]
    @test sort(collect(LightGraphs.inneighbors(b,9))) == [3,4,5]

    @test Set(LightGraphs.vertices(b)) == Set(1:LightGraphs.nv(b))

    @test shortest_path(b,1,5) == [1,3,5]
    @test LightGraphs.is_weakly_connected(b) == true

    @test SimpleHypergraphs.add_vertex!(h1) == 6
    @test add_hyperedge!(h1) == 5
    h1[5,5] = 1
    h1[6,5] = 1

    @test shortest_path(b,1,6) == [1,3,5,6]

end;

@testset "SimpleHypergraphs TwoSectionView" begin

    t = TwoSectionView(h1)

    @test LightGraphs.nv(t) == 6
    @test LightGraphs.ne(t) == 8

    @test sort(LightGraphs.all_neighbors(t, 1)) == [2,3]
    @test sort(LightGraphs.outneighbors(t, 5)) == [3,4,6]
    @test sort(LightGraphs.inneighbors(t, 4)) == [2,3,5]
    @inferred LightGraphs.all_neighbors(t, 1)

    @test LightGraphs.has_edge(t, 1, 2) == true
    @test LightGraphs.has_edge(t, 1, 5) == false

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(t))) == 16
    @test shortest_path(t,1,5) == [1,3,5]
    @test LightGraphs.is_weakly_connected(t) == true

    @test SimpleHypergraphs.add_vertex!(h1) == 7
    h1[7,5] = 1

    @test shortest_path(t,1,6) == [1,3,5,6]

    @test LightGraphs.nv(t) == 7
    @test LightGraphs.ne(t) == 10
    @test sort(LightGraphs.outneighbors(t, 5)) == [3,4,6,7]

    @test sum(LightGraphs.adjacency_matrix(LightGraphs.SimpleGraph(t))) == 20

end;


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 849e15e4-f786-11ea-3e89-315057d10f40
using Pkg, DrWatson

# ╔═╡ 849e4d66-f786-11ea-1fff-71edc16a1bdb
begin
	@quickactivate "SR2TuringPluto"
	using StatisticalRethinking
end

# ╔═╡ 2b094a8a-f786-11ea-3688-3769dfa0c5f8
md"## Fig3.2t.jl"

# ╔═╡ 849eb8c8-f786-11ea-1b05-e7a8c5c1d987
begin
	p_grid = range(0, stop=1, length=10000)
	prior = ones(length(p_grid))
	likelihood = pdf.(Binomial.(9, p_grid), 6)
	posterior = likelihood .* prior
	posterior = posterior / sum(posterior)
	samples = sample(p_grid, Weights(posterior), length(p_grid));
end;

# ╔═╡ 84abc770-f786-11ea-3a46-b5f847b0f5b8
begin
	N = 10000
	samples2 = sample(p_grid, Weights(posterior), N)
end;

# ╔═╡ 84ac5e4c-f786-11ea-37bc-3318feecfff9
begin
	b1 = mapreduce(p -> p < 0.5 ? 1 : 0, +, samples2) / N
	b2 = mapreduce(p -> (p > 0.5 && p < 0.75) ? 1 : 0, +, samples2) / N
	b3 = quantile(samples2, 0.8)
	b4 = quantile(samples2, [0.1, 0.9])

	p1 = plot_density_interval(samples2, [0.0, 0.5],
	  xlab="Proportion water (p)");
	p2 = plot_density_interval(samples2, [0.5, 0.75],
	  xlab="Proportion water (p)");
	p3 = plot_density_interval(samples2, [0.0, b3],
	  xlab="Proportion water (p)");
	p4 = plot_density_interval(samples2, b4,
	  xlab="Proportion water (p)")
end;

# ╔═╡ 84ba27de-f786-11ea-00e5-ff77781b2b88
plot(p1, p2, p3, p4, layout=(2, 2))

# ╔═╡ 84c23eb0-f786-11ea-262c-2f906a343b13
md"## End of Fig3.2t.jl"

# ╔═╡ Cell order:
# ╟─2b094a8a-f786-11ea-3688-3769dfa0c5f8
# ╠═849e15e4-f786-11ea-3e89-315057d10f40
# ╠═849e4d66-f786-11ea-1fff-71edc16a1bdb
# ╠═849eb8c8-f786-11ea-1b05-e7a8c5c1d987
# ╠═84abc770-f786-11ea-3a46-b5f847b0f5b8
# ╠═84ac5e4c-f786-11ea-37bc-3318feecfff9
# ╠═84ba27de-f786-11ea-00e5-ff77781b2b88
# ╟─84c23eb0-f786-11ea-262c-2f906a343b13


text:
mutable struct XWindowHandler <: AbstractWindowHandler
    conn::Connection
    windows::Vector{XCBWindow}
    keymap::Keymap
    callbacks::Dict{XCBWindow, WindowCallbacks}
end

XWindowHandler(conn::Connection, windows::Vector{XCBWindow}) = XWindowHandler(conn, windows, Keymap(conn), Dict())

function XWindowHandler(conn::Connection, windows::Vector{XCBWindow}, callbacks::Dict{XCBWindow, WindowCallbacks})
    wh = XWindowHandler(conn, windows, Keymap(conn))
    for (win, cb) ∈ callbacks
        set_callbacks!(wh, win, cb)
    end
end

function set_callbacks!(wh::XWindowHandler, win::XCBWindow, callbacks::WindowCallbacks)
    wh.callbacks[win] = callbacks
    set_event_mask(win, callbacks)
end

function poll_for_event(wh::XWindowHandler)
    while true
        event = xcb_poll_for_event(wh.conn)
        event ≠ C_NULL && return event
        yield()
    end
end

function wait_for_event(wh::XWindowHandler)
    event = xcb_wait_for_event(wh.conn)
    event == C_NULL ? nothing : event
end

function terminate_window!(wh::XWindowHandler, win::XCBWindow)
    deleteat!(wh.windows, window_index(wh, win))
    finalize(win)
end

window_index(wh::XWindowHandler, win::XCBWindow) = findfirst(x -> x.id == win.id, wh.windows)
window_index(wh::XWindowHandler, id::Integer) = findfirst(x -> x.id == id, wh.windows)

get_window(wh::XWindowHandler, id::Integer) = wh.windows[window_index(wh, id)]
get_window(wh::XWindowHandler, event::xcb_xkb_state_notify_event_t) = nothing
get_window(wh::XWindowHandler, event::xcb_keymap_notify_event_t) = nothing
get_window(wh::XWindowHandler, event) = get_window(wh, window_id(event))

callbacks(wh::XWindowHandler, win::XCBWindow) = get(wh.callbacks, win, WindowCallbacks())

"""
Extract a keycode from an key event.
"""
keycode(wh::XWindowHandler, details::EventDetails{<:KeyEvent}) = keycode(wh.keymap, details.data.key_name)

"""
Return condensed information regarding a keystroke as a `String`. Includes key name, keycode, character input and key symbol.
"""
keystroke_info(wh::XWindowHandler, event_details::EventDetails) = keystroke_info(wh.keymap, event_details)


text:
import Base: ==, !=

const RANGE_DELIMITER = '/'
const DATE_STRING_LENGTHS = (4, 7, 10)
const DATETIME_STRING_LENGTHS = (4, 7, 10, 13, 16, 19)

function thrudt(s::AbstractString, t::Vector{Date})
    n = length(s)
    @assert n in DATE_STRING_LENGTHS "Invalid indexing string: Unable to parse $s"
    if n == 4  # yyyy given
        y = parse(Int, s)
        ymd = Date(y, 12, 31)
    elseif n == 7  # yyyy-mm given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        ymd = lastdayofmonth(Date(y, m, 1))
    elseif n == 10  # yyyy-mm-dd given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        d = parse(Int, a[3])
        ymd = Date(y, m, d)
    end
    return t .<= ymd
end

function fromdt(s::AbstractString, t::Vector{Date})
    n = length(s)
    @assert n in DATE_STRING_LENGTHS "Invalid indexing string: Unable to parse $s"
    if n == 4  # yyyy given
        y = parse(Int, s)
        ymd = Date(y, 1, 1)
    elseif n == 7  # yyyy-mm given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        ymd = Date(y, m, 1)
    elseif n == 10  # yyyy-mm-dd given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        d = parse(Int, a[3])
        ymd = Date(y, m, d)
    end
    return t .>= ymd
end

function thisdt(s::AbstractString, t::Vector{Date})
    n = length(s)
    @assert n in DATE_STRING_LENGTHS "Invalid indexing string: Unable to parse $s"
    if n == 4  # yyyy given
        y = parse(Int, s)
        return year.(t) .== y
    elseif n == 7  # yyyy-mm given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        return (year.(t) .== y) .& (month.(t) .== m)
    elseif n == 10  # yyyy-mm-dd given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        d = parse(Int, a[3])
        idx = (year.(t) .== y) .& (month.(t) .== m) .& (day.(t) .== d)
    end
end

function fromthru(from::AbstractString, thru::AbstractString, t::Vector{Date})
    fromidx = fromdt(from, t)
    thruidx = thrudt(thru, t)
    return fromidx .& thruidx
end

function thrudt(s::AbstractString, t::Vector{DateTime})
    n = length(s)
    @assert n in DATETIME_STRING_LENGTHS || n >= 20 "Invalid indexing string: Unable to parse $s"
    if n == 4  # yyyy given
        y = parse(Int, s)
        ymdhms = DateTime(y, 12, 31)
    elseif n == 7  # yyyy-mm given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        ymdhms = lastdayofmonth(DateTime(y, m, 1))
    elseif n == 10  # yyyy-mm-dd given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        d = parse(Int, a[3])
        ymdhms = DateTime(y, m, d)
    elseif n == 13  # yyyy-mm-ddTHH given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = 59
        sec = 59
        ymdhms = DateTime(y, m, d, hr, min, sec)
    elseif n == 16  # yyyy-mm-ddTHH:MM given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = 59
        ymdhms = DateTime(y, m, d, hr, min, sec)
    elseif n == 19  # yyyy-mm-ddTHH:MM:SS given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = parse(Int, c[3])
        ymdhms = DateTime(y, m, d, hr, min, sec)
    elseif n >= 20  # milliseconds given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        f = split(c[3], '.')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = parse(Int, f[1])
        milli = parse(Int, f[2])
        ymdhms = DateTime(y, m, d, hr, min, sec, milli)
    end
    return t .<= ymdhms
end

function fromdt(s::AbstractString, t::Vector{DateTime})
    n = length(s)
    @assert n in DATETIME_STRING_LENGTHS || n >= 20 "Invalid indexing string: Unable to parse $s"
    if n == 4  # yyyy given
        y = parse(Int, s)
        ymdhms = DateTime(y)
    elseif n == 7  # yyyy-mm given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        ymdhms = DateTime(y, m)
    elseif n == 10  # yyyy-mm-dd given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        d = parse(Int, a[3])
        ymdhms = DateTime(y, m, d)
    elseif n == 13  # yyyy-mm-ddTHH given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        ymdhms = DateTime(y, m, d, hr)
    elseif n == 16  # yyyy-mm-ddTHH:MM given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        ymdhms = DateTime(y, m, d, hr, min)
    elseif n == 19  # yyyy-mm-ddTHH:MM:SS given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = parse(Int, c[3])
        ymdhms = DateTime(y, m, d, hr, min, sec)
    elseif n >= 20  # milliseconds given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        f = split(c[3], '.')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = parse(Int, f[1])
        milli = parse(Int, f[2])
        ymdhms = DateTime(y, m, d, hr, min, sec, milli)
    end
    return t .>= ymdhms
end

function thisdt(s::AbstractString, t::Vector{DateTime})
    n = length(s)
    if n == 4  # yyyy given
        y = parse(Int, s)
        return year.(t) .== y
    elseif n == 7  # yyyy-mm given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        return (year.(t) .== y) .& (month.(t) .== m)
    elseif n == 10  # yyyy-mm-dd given
        a = split(s, '-')
        y = parse(Int, a[1])
        m = parse(Int, a[2])
        d = parse(Int, a[3])
        return (year.(t) .== y) .& (month.(t) .== m) .& (day.(t) .== d)
    elseif n == 13  # yyyy-mm-ddTHH given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        return (year.(t) .== y) .& (month.(t) .== m) .& (day.(t) .== d) .& (hour.(t) .== hr)
    elseif n == 16  # yyyy-mm-ddTHH:MM given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        return (year.(t) .== y) .& (month.(t) .== m) .& (day.(t) .== d) .& (hour.(t) .== hr) .& (minute.(t) .== min)
    elseif n == 19  # yyyy-mm-ddTHH:MM:SS given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = parse(Int, c[3])
        return (year.(t) .== y) .& (month.(t) .== m) .& (day.(t) .== d) .& (hour.(t) .== hr) .& (minute.(t) .== min) .& (second.(t) .== sec)
    elseif n >= 20  # milliseconds given
        a = split(s, 'T')
        b = split(a[1], '-')
        c = split(a[2], ':')
        f = split(c[3], '.')
        y = parse(Int, b[1])
        m = parse(Int, b[2])
        d = parse(Int, b[3])
        hr = parse(Int, c[1])
        min = parse(Int, c[2])
        sec = parse(Int, f[1])
        milli = parse(Int, f[2])
        return (year.(t) .== y) .& (month.(t) .== m) .& (day.(t) .== d) .& (hour.(t) .== hr) .& (minute.(t) .== min) .& (second.(t) .== sec) .& (millisecond.(t) .== milli)
    end
end

function fromthru(from::AbstractString, thru::AbstractString, t::Vector{DateTime})
    fromidx = fromdt(from, t)
    thruidx = thrudt(thru, t)
    return fromidx .& thruidx
end

function dtidx(s::AbstractString, t::Vector{Date})
    if s == "/" || s == ""
        return trues(length(t))
    end
    @assert !in('T', s) "Cannot index Date type with sub-daily frequency."
    bds = split(s, RANGE_DELIMITER)
    @assert length(bds) in [1,2]
    if length(bds) == 1  # single date
        return thisdt(s, t)
    elseif length(bds) == 2  # date range
        n1 = length(bds[1])
        n2 = length(bds[2])
        if n1 == 0 && n2 != 0  # thru date given
            return thrudt(bds[2], t)
        elseif n1 != 0 && n2 == 0  # from date given
            return fromdt(bds[1], t)
        elseif n1 != 0 && n2 != 0  # from and thru date given
            return fromthru(bds[1], bds[2], t)
        end
    end
end

function dtidx(s::AbstractString, t::Vector{DateTime})
    if s == "/" || s == ""
        return trues(length(t))
    end
    bds = split(s, RANGE_DELIMITER)
    @assert length(bds) in [1,2]
    if length(bds) == 1  # single date
        return thisdt(s, t)
    elseif length(bds) == 2  # date range
        n1 = length(bds[1])
        n2 = length(bds[2])
        if n1 == 0 && n2 != 0  # thru date given
            return thrudt(bds[2], t)
        elseif n1 != 0 && n2 == 0  # from date given
            return fromdt(bds[1], t)
        elseif n1 != 0 && n2 != 0  # from and thru date given
            return fromthru(bds[1], bds[2], t)
        end
    end
end


text:
@testset "Test linear_index" begin
    a = rand(10,10,10);
    @test a[3,4,5]  == a[433]
    @test linear_index((3,4,5), (10,10,10)) == 433
end
@testset "Test cast_iter" begin
    i1 = NDTools.cast_iter([1;2;3;4])
    i2 = NDTools.cast_iter([1,2,3,4])
    i3 = NDTools.cast_iter([(1,),(2,),(3,),(4,)])
    for (i,j,k) in zip(i1,i2,i3)
        @test i == j
        @test j == k[1]
    end
end

@testset "Test cast_number_iter" begin
    i1 = NDTools.cast_number_iter([1;2;3;4])
    i2 = NDTools.cast_number_iter([1,2,3,4])
    i3 = NDTools.cast_number_iter((1,2,3,4))
    for (i,j,k) in zip(i1,i2,i3)
        @test i == j
        @test j == k
    end
end

@testset "Test optional_mat_to_iter" begin
    @test all(NDTools.optional_mat_to_iter([1;2;3;4]) .== NDTools.cast_iter([1;2;3;4]))
end

@testset "Test mat_to_tvec" begin
    @test NDTools.mat_to_tvec([1 3 5; 2 4 6]) == [(1,2),(3,4),(5,6)]
end

@testset "Test apply_dims" begin
    @test NDTools.apply_dims((1,2,3),(1,2), 4)  == (1, 2, 0, 0)
    @test NDTools.apply_dims([(1,2,3),(1,2,3),(1,2,3)],(1,2), 4)  == ((1, 2, 0, 0),(1, 2, 0, 0),(1, 2, 0, 0))
end


text:
export
    HMM,
    LinearGaussian,
    NonLinearGaussian,
    GaussianHMM,
    generate

abstract type AbstractHMM end
abstract type GaussianHMM end

struct HMM <: AbstractHMM
    transmean::Function
    transloglik::Function # log transition function:  f(xk|xkm1)
    obsmean::Function
    obsloglik::Function   # log observation function: g(y|xk)
    dimx::Int
    dimy::Int
end

struct LinearGaussian <: GaussianHMM
    #=
        x <- Ax + chol(Q)' * randn
        y <- Bx + chol(R)' * randn
    =#
    A::Matrix{Float}
    B::Matrix{Float}
    Q::Union{Float,Matrix{Float}}
    R::Union{Float,Matrix{Float}}
    # implicit
    dimx::Int
    dimy::Int
    cholQ::Union{Float,UpperTriangular{Float, Matrix{Float}}}
    cholR::Union{Float,UpperTriangular{Float, Matrix{Float}}}
    transmean::Function
    obsmean::Function
    function LinearGaussian(A,B,Q,R)
        dimx=size(A,1)
        dimy=size(B,1)
        @assert dimx==size(A,2)==size(Q,1)==size(Q,2) "dimensions don't match"
        @assert dimy==size(R,1)==size(R,2) "dimensions don't match"
        @assert issymmetric(Q) && issymmetric(R) "cov mat must be symmetric"
        @assert isposdef(Q) && isposdef(R) "cov mat must be pos def"
        new(A,B,Q,R,dimx,dimy,chol(Q),chol(R),
            (k,xkm1) -> A*xkm1,
            (k,xk)   -> B*xk )
    end
end

struct NonLinearGaussian <: GaussianHMM
    #=
        x <- f(x) + chol(Q)' * randn
        y <- g(x) + chol(R)' * randn
    =#
    transmean::Function
    obsmean::Function
    Q::Union{Float,Matrix{Float}}
    R::Union{Float,Matrix{Float}}
    # implicit
    dimx::Int
    dimy::Int
    cholQ::Union{Float,UpperTriangular{Float, Matrix{Float}}}
    cholR::Union{Float,UpperTriangular{Float, Matrix{Float}}}
    function NonLinearGaussian(f, g, Q, R)
        dimx = size(Q, 1)
        dimy = size(R, 1)
        @assert issymmetric(Q) && issymmetric(R) "cov mat must be symmetric"
        @assert isposdef(Q) && isposdef(R) "cov mat must be pos def"
        new((k,xkm1)-> f(k,xkm1),
            (k,xk)  -> g(k,xk),
            Q,R,dimx,dimy,chol(Q),chol(R)  )
    end
end

function HMM(g::GaussianHMM)
    transloglik = (k,xkm1,xk) -> -norm(g.cholQ'\(xk - g.transmean(k,xkm1)))^2/2
    obsloglik   = (k,yk,  xk) -> -norm(g.cholR'\(yk - g.obsmean(k,xk)))^2/2
    HMM(g.transmean, transloglik, g.obsmean, obsloglik, g.dimx, g.dimy)
end

### Generation of observations

"""
    generate(lg, x0, K)

Generate observations following a given (non)linear Gaussian dynamic for `K`
time steps.
"""
function generate(g::GaussianHMM, x0::Vector{Float}, K::Int
                    )::Tuple{Matrix{Float},Matrix{Float}}
    @assert length(x0)==g.dimx "dimensions don't match"
    # allocate states/observations
    states, observations = zeros(g.dimx, K), zeros(g.dimy, K)
    # assign first state
    states[:,1] = x0
    # pre-generate noise
    noisex = g.cholQ' * randn(g.dimx,K)
    noisey = g.cholR' * randn(g.dimy,K)
    # use noise in iterative linear system
    for k = 1:(K-1)
        observations[:,k] = g.obsmean(k, states[:,k]) + noisey[:,k]
        states[:,k+1]     = g.transmean(k+1, states[:,k]) + noisex[:,k+1]
    end
    # last observation
    observations[:,K] = g.obsmean(K, states[:,K]) + noisey[:,K]
    # package and return
    return (states, observations)
end


text:
# script to be run in Julia REPL at start of the session making it
# easier to test bits of code
using Revise # reloading modules now check for changes
using JSON
using JuMP
using Gurobi
push!(LOAD_PATH, "./") # modules below are defined in the current folder
using Chen1995
using Kurpel2018
using Check3DPackings
using CornerPoints

instances_folder = "../../../instances/hbd_basic_tests/"
instances_fnames = filter(x -> endswith(x, ".json"), readdir(instances_folder))
json_instances = []
for f in instances_fnames
  s = open(x -> read(x, String), instances_folder * f)
  push!(json_instances, JSON.parse(s))
end

function json2model(json;
  no_rotation = false,
  max_packed_volume = false,
  formulation = Chen1995
)
  model = Model(solver = GurobiSolver())
  pqrLWH = map(l -> json[l], split("pqrLWH", ""))

  formulation.build(model, pqrLWH...;
    no_rotation = no_rotation, max_packed_volume = max_packed_volume
  )
end


text:
#=
Define a discrete state space for the PedCarMDP problem
=#

function n_states(mdp::PedCarMDP)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    n_ped = n_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    routes = get_car_routes(mdp.env)
    n_cars = 0
    for route in routes
        n_cars += n_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
    end
    return n_ego*(n_cars + 1)*(n_ped + 1) # do not forget absent state
end

function POMDPs.states(mdp::PedCarMDP)
    state_space = PedCarMDPState[]
    ego_states = get_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    ped_states = get_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    push!(ped_states, get_off_the_grid(mdp))
    routes = get_car_routes(mdp.env)
    for route in routes
        for ego in ego_states
            for ped in ped_states
                car_states = get_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
                for car in car_states
                    collision =  is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(car, mdp.car_type, CAR_ID)) || is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(ped, mdp.ped_type, PED_ID))
                    # collision = crash(mdp, ego, car, ped)
                    # enumerate all possible routes
                    lane = get_lane(mdp.env.roadway, car)
                    push!(state_space, PedCarMDPState(collision, ego, ped, car, SVector{2, LaneTag}(route[1], route[end])))
                end
            end
        end
    end
    for ego in ego_states
        for ped in ped_states
            # add absent states
            collision = is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(ped, mdp.ped_type, PED_ID))
            # collision = collision_checker(ego, ped, mdp.ego_type, mdp.ped_type)
            push!(state_space, PedCarMDPState(collision, ego, ped, get_off_the_grid(mdp), SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0, 0))))
        end
    end
    return state_space
end

function POMDPs.stateindex(mdp::PedCarMDP, s::PedCarMDPState)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    n_ped = n_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    routes = get_car_routes(mdp.env)
    # step 1: find ego_index
    ego_i = ego_stateindex(mdp.env, s.ego, mdp.pos_res, mdp.vel_res)
    # step 2: find ped index
    if s.ped.posG == mdp.off_grid
        ped_i = n_ped + 1
    else
        ped_i = ped_stateindex(mdp.env, s.ped, mdp.pos_res, mdp.vel_ped_res)
    end

    # step 3 find route index
    route_i = 0
    for (i, route) in enumerate(routes)
        if s.route[end] == route[end] && s.route[1] == route[1]
            route_i = i
        end
    end

    # handle off the grid case
    if s.car.posG == mdp.off_grid || s.route == [LaneTag(0, 0), LaneTag(0, 0)]
        si = 0
        for route in routes
            si += n_ego * (n_ped + 1) * n_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
        end
        si += LinearIndices((n_ped + 1, n_ego))[ped_i, ego_i]
    else
        n_car = n_car_states(mdp.env, find_route(mdp.env, s.route), mdp.pos_res, mdp.vel_res)
        # step 2: find route_index
        route_i = 0
        for (i, route) in enumerate(routes)
            if  s.route[end] == route[end] && s.route[1] == route[1]
                route_i = i
            end
        end
        # step 3: find car_index in car states
        car_i = car_stateindex(mdp.env, s.car, find_route(mdp.env, s.route), mdp.pos_res, mdp.vel_res)
        # use linear/Cartesian magic
        si = LinearIndices((n_car, n_ped + 1, n_ego))[car_i, ped_i, ego_i]

        for i=2:route_i
            size_r = n_ego * (n_ped + 1) * n_car_states(mdp.env, routes[i-1], mdp.pos_res, mdp.vel_res)
            si += size_r
        end
    end
    return si
end

function ind2state(mdp::PedCarMDP, si::Int64)
    n_ego = n_ego_states(mdp.env, mdp.pos_res, mdp.vel_res)
    n_ped = n_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    routes = get_car_routes(mdp.env)
    n_routes = length(routes)
    car, ped, ego = nothing, nothing, nothing
    # find route first
    ns = 0
    route_ind = nothing
    route_shift = 0
    for (i, route) in enumerate(routes)
        n_cars = n_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
        route_shift = ns
        ns += n_cars*n_ego*(n_ped + 1)
        if ns >= si
            route_ind = i
            break
        end
    end
    # find car, ped, ego
    if route_ind == nothing # route was not found, car is off the grid
        si_ = si - ns # shift by all the states that were added before
        car = get_off_the_grid(mdp)
        # retrieve ped and ego
        ped_i, ego_i = Tuple(CartesianIndices((n_ped + 1, n_ego))[si_])
        ego = ind2ego(mdp.env, ego_i, mdp.pos_res, mdp.vel_res)
        if ped_i > n_ped
            ped = get_off_the_grid(mdp)
        else
            ped = ind2ped(mdp.env, ped_i, mdp.pos_res, mdp.vel_ped_res)
        end
        collision = crash(mdp, ego, car, ped)
        return PedCarMDPState(collision, ego, ped, car, SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0, 0)))
    else
        si_ = si - route_shift
        route = routes[route_ind]
        sroute = SVector{2, LaneTag}(route[1], route[end])
        n_cars = n_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
        # cartesian/linear magic
        car_i, ped_i, ego_i = Tuple(CartesianIndices((n_cars, n_ped + 1, n_ego))[si_])
        car = ind2car(mdp.env, car_i, route, mdp.pos_res, mdp.vel_res)
        ego = ind2ego(mdp.env, ego_i, mdp.pos_res, mdp.vel_res)
        if ped_i > n_ped
            ped = get_off_the_grid(mdp)
        else
            ped = ind2ped(mdp.env, ped_i, mdp.pos_res, mdp.vel_ped_res)
        end
        collision = crash(mdp, ego, car, ped)
        return PedCarMDPState(collision, ego, ped, car, sroute)
    end
end

#### INITIAL STATES

function POMDPs.initialstate(mdp::PedCarMDP)
    ego = initial_ego_state(mdp)
    init_car_states, init_car_routes = initial_car_state_distribution(mdp)
    init_ped_dist = initial_ped_state_distribution(mdp)
    init_states = Vector{PedCarMDPState}()
    for i=1:length(init_car_states)
        for j=1:length(init_ped_dist.vals)
            ped = init_ped_dist.vals[j]
            car = init_car_states[i]
            route = init_car_routes[i]
            crash = is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(car, mdp.car_type, CAR_ID)) || is_colliding(Vehicle(ego, mdp.ego_type, EGO_ID), Vehicle(ped, mdp.ped_type, PED_ID))
            push!(init_states, PedCarMDPState(crash, ego, ped, car, SVector{2, LaneTag}(route[1], route[end])))
        end
    end
    # uniform
    probs = ones(length(init_states))
    normalize!(probs, 1)
    return SparseCat(init_states, probs)
end

function initial_ego_state(mdp::PedCarMDP)
    lanes = get_ego_route(mdp.env)
    posF = Frenet(mdp.env.roadway[lanes[1]], mdp.ego_start)
    v0 = 0.
    return VehicleState(posF, mdp.env.roadway, v0)
end

function initial_car_state_distribution(mdp::PedCarMDP)
    routes = get_car_routes(mdp.env)
    init_car_routes = []
    init_car_states = VehicleState[]
    for route in routes
        car_states = get_car_states(mdp.env, route, mdp.pos_res, mdp.vel_res)
        for cs in car_states
            push!(init_car_states, cs)
            push!(init_car_routes, route)
        end
    end
    push!(init_car_states, get_off_the_grid(mdp))
    push!(init_car_routes, SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0,0)))
    return init_car_states, init_car_routes
end

function initial_car_state(mdp::PedCarMDP, rng::AbstractRNG)
    init_car_states, init_car_routes = initial_car_state_distribution(mdp)
    return rand(rng, init_car_states), rand(rng, init_car_routes)
end

function car_starting_states(mdp::PedCarMDP, min_speed::Float64 = 6.0)
    # list of car starting states
    routes = get_car_routes(mdp.env)
    v_space = min_speed:mdp.vel_res:mdp.env.params.speed_limit
    N_states = 0
    for route in routes
        N_states += length(v_space)
    end
    car_states = Vector{VehicleState}(N_states)
    start_routes = Vector{SVector{2, LaneTag}}(N_states)
    i = 1
    for route in routes
        lane = route[1]
        for v in v_space
            car_states[i] = VehicleState(Frenet(mdp.env.roadway[lane], 0.), mdp.env.roadway, v)
            start_routes[i] = SVector{2, LaneTag}(route[1], route[end])
            i += 1
        end
    end
    return car_states, start_routes
end

function initial_ped_state_distribution(mdp::PedCarMDP)
    init_ped_states = get_ped_states(mdp.env, mdp.pos_res, mdp.vel_ped_res)
    push!(init_ped_states, get_off_the_grid(mdp))
    # uniform (maybe add more weights to the states when pedestrians are not there?)
    probs = ones(length(init_ped_states))
    normalize!(probs, 1)
    return SparseCat(init_ped_states, probs)
end

function initial_ped_state(mdp::PedCarMDP, rng::AbstractRNG)
    init_dist = initial_ped_state_distribution(mdp)
    return rand(rng, init_dist)
end

function pedestrian_starting_states(mdp::PedCarMDP)
    # list of pedestrian starting states
    n_headings = 2
    lanes = get_ped_lanes(mdp.env)
    v_space = get_ped_vspace(mdp.env, mdp.vel_ped_res)
    ped_states = Vector{VehicleState}(length(lanes)*length(v_space)*n_headings)
    i = 1
    for lane in lanes
        dlane = get_discretized_lane(lane, mdp.env.roadway, mdp.pos_res)
        for v in v_space
            ped_states[i] = VehicleState(Frenet(mdp.env.roadway[lane], 0., 0., 0.), mdp.env.roadway, v)
            i += 1
            ped_states[i] = VehicleState(Frenet(mdp.env.roadway[lane], dlane[end], 0., float(pi)), mdp.env.roadway, v)
            i += 1
        end
    end
    return ped_states
end


text:
# This file was generated, do not modify it. # hide
using MLJ, PyPlot, PrettyPrinting, Random,
      DataFrames

X, y = @load_boston
sch = schema(X)
p = length(sch.names)
n = sch.nrows
@show (n, p)
describe(y)


text:
# mri/z-test.jl

using Test: @testset

list = [
"exp_mult.jl"
"exp_xform.jl"
"coil_compress.jl"
"image_geom_mri.jl"
"kspace.jl"
"kspace-spiral.jl"
"mri_objects.jl"
"mri_trajectory.jl"
"sensemap-sim.jl"
]

for file in list
	@testset "$file" begin
		include(file)
	end
end


text:
import TPT

using Base.Test

using Plots; pyplot()

println("--- Unit LWCA ---")

#
# Lennard-Jones fluid
# Ref: A. Meyer et al. Chemical Physics 49 (1980) 147-152
#
kB = TPT.kB
ϵ = 119.8kB
r₀ = 3.405 / 0.5292
n = 0.844 / r₀^3
T = 0.723 * ϵ / kB

lj = TPT.LennardJones(ϵ, r₀)

u = TPT.pairpotential(lj)[1,1]
rmin = TPT.pairpotential_minimizer(lj)[1,1]
umin = u(rmin)
u′ = TPT.pairpotential_derivative(lj)[1,1]
plot([u, u′], 5, 20, ylims=(umin, -umin))

ahs = TPT.AHS(σ = r₀, ρ = n, approx=:RFA)
lwca = TPT.LWCA(ahs, T, struct=:linear)
wca = TPT.WCA(ahs, T)

sys = TPT.TPTSystem(lwca, lj)
sys2 = TPT.TPTSystem(wca, lj)

σ = sys.ref.trial.σ[1]
σ₀ = 3.474 / 0.5292

σ2 = sys2.ref.trial.σ[1]

B = TPT.blipfunction(sys.ref)[1,1]
B2 = TPT.blipfunction(sys2.ref)[1,1]
plot([B, B2], 5, 8)

g0 = TPT.paircorrelation(ahs)[1,1]
g1 = TPT.paircorrelation(sys)[1,1]
g2 = TPT.paircorrelation(sys2)[1,1]
plot([g0, g1, g2], 5, 15, labels = ["AHS" "LWCA" "WCA"])

S0 = TPT.structurefactor(ahs)[1,1]
S1 = TPT.structurefactor(sys)[1,1]
S2 = TPT.structurefactor(sys2)[1,1]
plot([S0, S1, S2], 0.1, 5, labels=["AHS" "LWCA" "WCA"])

@testset "Unit LWCA" begin
  @test isapprox(σ, 6.55, atol=1e-2)
end


text:
"""
    compute_laplace_eig(mesh, matrices, pde, eiglim, neig_max)

Compute the Laplace eigenvalues, eigenfunctions and first order moments of products of pairs of eigenfunctions.
"""
function compute_laplace_eig(model, matrices, eiglim = Inf, neig_max = Inf)

    # Measure function evaluation time
    starttime = Base.time()

    # Extract parameters
    @unpack mesh, D, T₂ = model
    @unpack M, S, R, Mx, Q = matrices
    ncompartment = length(mesh.points)

    # Compute at most all eigenvalues in the given domain
    neig = Int(min(neig_max, size(M, 1)))

    println("Solving Laplace eigenvalue problem, computing $neig eigenvalues.")
    println("Problem size: $(size(M, 1)) points.")

    # Solve generalized eigenproblem, computing the smallest eigenvalues only.
    # If 2 * neig_max_domain >= nnode, a full decomposition is performed,
    # calling the eig function inside eigs
    λ, ϕ = eigs(S + Q, M, nev = neig, which = :SR)
    # λ, ϕ = eigs(Hermitian(S + Q), Hermitian(M), nev = neig, which = :SR);
    # λ, ϕ = eigen(Hermitian(Matrix(S + Q)), Hermitian(Matrix(M)))

    # All Laplace eigenvalues are nonnegative
    all(0 .≤ λ) ||
        @warn "Obtained negative eigenvalues for Laplace operator." findall(λ .< 0) λ[λ.<0]

    # Only keep modes with length scales larger than minimum
    inds = λ .≤ eiglim
    λ = λ[inds]
    ϕ = ϕ[:, inds]
    isinf(eiglim) ||
        length(λ) < neig ||
        @warn "No eigenvalues were outside the interval. Consider increasing `neig_max`." eiglim neig_max

    # Normalize eigenfunctions with mass weighting
    ϕ ./= .√sum(ϕ .* (M * ϕ), dims = 1)

    # Compute first order moments of product of pairs of eigenfunctions
    moments = cat([ϕ' * Mx[dim] * ϕ for dim = 1:3]..., dims = 3)

    # Compute Laplace relaxation matrix
    massrelax = ϕ' * R * ϕ

    values = λ
    funcs = ϕ
    totaltime = Base.time() - starttime

    (; values, funcs, moments, massrelax, totaltime)
end


text:
module ALPS

# package code goes here

include("write_inputfiles.jl")

end # module


text:
## Efficient implementation of node and edge attributes ##

length(x::AbstractAttribute) = x.getlen()
lastindex(x::AbstractAttribute) = x.getlen()
size(x::AbstractAttribute) = (x.getlen(),)

"""
    ConstantAttribute{T}

Attribute always returning the same (default) value.
"""
mutable struct ConstantAttribute{T,N,F} <: AbstractAttribute{T,N,F}
    default :: T
    getlen  :: F
end

ConstantAttribute(d::T, f::F) where {T,F<:Function} = ConstantAttribute{T,1,F}(d, f)

getindex(x::ConstantAttribute, ::Integer) = x.default
setindex!(::ConstantAttribute, ::Any, ::Integer) =
    error("constant attributes don't support element assignments")
setdefault!(x::ConstantAttribute{T}, d::T) where {T} = x.default = d
deleteat!(::ConstantAttribute, ::Integer) = nothing

"""
    SparseAttribute{T}

Attribute having most of its elements equal to the default value.
Under the hood uses `Dict` to store non-default elements.
"""
mutable struct SparseAttribute{T,N,F} <: AbstractAttribute{T,N,F}
    default :: T
    data    :: Dict{Int,T}
    getlen  :: F
end

SparseAttribute(d::T, f::F) where {T,F<:Function} = SparseAttribute{T,1,F}(d, Dict{Int,T}(), f)
SparseAttribute(x::ConstantAttribute{T,N,F}) where {T,N,F} =
    SparseAttribute{T,N,F}(x.default, Dict{Int,T}(), x.getlen)
density(x::SparseAttribute) = length(x.data) / x.getlen()

getindex(x::SparseAttribute, i::Integer) = get(x.data, i, x.default)
function setindex!(x::SparseAttribute{T}, v, i::Integer) where {T}
    @boundscheck checkbounds(1:x.getlen(), i)
    if v != x.default
        x.data[i] = v
    else
        delete!(x.data, i)
    end
end
setdefault!(x::SparseAttribute{T}, d::T) where {T} = x.default = d
deleteat!(x::SparseAttribute, i::Integer) = delete!(x.data, i)

"""
    DenseAttribute{T}

Attribute with most of its values being unique. Under the hood, uses `Vector`
to store all elements.
"""
mutable struct DenseAttribute{T,N,F} <: AbstractAttribute{T,N,F}
    default :: T
    data    :: Vector{T}
    getlen  :: F
end

DenseAttribute(d::T, f::F) where {T,F<:Function} = DenseAttribute{T,1,F}(d, T[], f)
DenseAttribute(d::T, x::Vector{T}, f::F) where {T,F<:Function} =
    DenseAttribute{T,1,F}(d, x, f)
DenseAttribute(x::AbstractVector{T}, f::F) where {T,F<:Function} =
    DenseAttribute(defval(T), collect(x), f)
DenseAttribute(x::AbstractAttribute{T,N,F}) where {T,N,F<:Function} =
    DenseAttribute{T,N,F}(x.default, collect(x), x.getlen)

getindex(x::DenseAttribute, i::Integer) = checkbounds(Bool, x.data, i) ? x.data[i] : x.default
function setindex!(x::DenseAttribute{T}, v::T, i::Integer) where {T}
    @boundscheck checkbounds(1:x.getlen(), i)
    while length(x.data) < i
        push!(x.data, x.default)
    end
    x.data[i] = v
end
function setdefault!(x::DenseAttribute{T}, d::T) where {T}
    for (i, v) = enumerate(x.data)
        if v == x.default
            x.data[i] = d
        end
    end
    x.default = d
end
deleteat!(x::DenseAttribute, i) = checkbounds(Bool, x, i) && deleteat!(x.data[i], i)

"""
    attribute(x, f)

Create an appropriate attribute instance depending on the type of `x`.
Anonymous function `f()` should return the length.
"""
attribute(x::AbstractAttribute, f::Function) = x  # TODO: use this to reassign attributes to different graphs
attribute(x::AbstractArray, f::Function) = DenseAttribute(x, f)
attribute(x, f::Function) = ConstantAttribute(x, f)

nodeattr(g::Graph, v) = attribute(v, ()->nodecount(g))
edgeattr(g::Graph, v) = attribute(v, ()->edgecount(g))

# TODO: 1. get rid of setindex!(d::AttributeDict, v, s::Symbol, i) in graph_operations.jl
# TODO: 2. introduce an `AttributeBuilder` class
# TODO: 3. find most common value on at least 100 fist ones and do an informed conversion
# TODO: 4. create a method that substantiates attributes from all `AttributeBuilder`s
# TODO: 5. replace currrent attribute assignment in all IO files