text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using Plots
using Distributions

"""empirical cumulative distribution function"""
ecdf(sample, x) = count(≤(x), sample)/length(sample)

λ = 3
dist_true = Poisson(λ)
n = 2^10
sample = rand(dist_true, n)
m, u² = mean(sample), var(sample)
gamma = Gamma(m^2/u², u²/m)
a, b = 0, 10

P = plot();
histogram!(sample; norm=true, alpha=0.3,
    bin=a-0.5:b+0.5, xtick=a:b,
    label="size-$n sample of Poisson($λ)");
plot!(x -> pdf(gamma, x), a, b;
    label="Gamma distribution approximation");

Q = plot(; legend=:bottomright, xtick=a:b);
plot!(x -> ecdf(sample, x), a-0.5, b+0.5;
    label="ecdf of size-$n sample of Poisson($λ)");
plot!(x -> cdf(gamma, x), a-0.5, b+0.5;
    label="cdf of Gamma distribution approximation");

plot(P, Q; size=(500, 500), layout=(2, 1))

# %%
using Plots
using Distributions
using QuadGK

"""empirical cumulative distribution function"""
ecdf(sample, x) = count(≤(x), sample)/length(sample)

λ = 3
dist_true = Poisson(λ)
n = 2^10
sample = rand(dist_true, n)
m, u² = mean(sample), var(sample)
gamma = Gamma(m^2/u², u²/m)

W_1 = quadgk(x -> abs(ecdf(sample, x) - cdf(gamma, x)), 0, Inf)[1]

a, b = 0, 10
xs = range(a-0.5, b+0.5; length=300)
zs = pdf.(gamma, xs)
Ys = ecdf.(Ref(sample), xs)
Zs = cdf.(gamma, xs)

P = plot();
histogram!(sample; norm=true, alpha=0.3, lw=2,
    bin=a-0.5:b+0.5, xtick=a:b,
    label="size-$n sample of Poisson($λ)");
plot!(xs, zs; lw=2,
    label="Gamma distribution approximation");

Q = plot(; legend=:bottomright, xtick=a:b, ytick=0:0.1:1);
plot!(xs, Ys; linealpha = 0, fillrange = Zs, fillcolor = :yellow,
    label="area = W_1 = Wasserstein distance for p = 1");
plot!(xs, Ys; lw = 2, color = 1,
    label="ecdf of size-$n sample of Poisson($λ)");
plot!(xs, Zs; lw = 2, color = 2,
    label="cdf of Gamma distribution approximation")
title!("W_1 = $(round(W_1; digits=3))")

plot(P, Q; size=(600, 800), layout=(2, 1))

# %%
using Plots
using Distributions
using QuadGK

"""empirical cumulative distribution function"""
ecdf(sample, x) = count(≤(x), sample)/length(sample)

λ = 3
dist_true = Poisson(λ)
n = 2^10
sample = rand(dist_true, n)
normal = fit(Normal, sample)

AIC = -2loglikelihood(normal, sample) + 2
W_1 = quadgk(x -> abs(ecdf(sample, x) - cdf(normal, x)), -Inf, Inf)[1]

a, b = -1, 10
xs = range(a-0.5, b+0.5; length=300)
zs = pdf.(normal, xs)
Ys = ecdf.(Ref(sample), xs)
Zs = cdf.(normal, xs)

P = plot();
histogram!(sample; norm=true, alpha=0.3, lw=2,
    bin=a-0.5:b+0.5, xtick=a:b,
    label="size-$n sample of Poisson($λ)");
plot!(xs, zs; lw=2,
    label="normal distribution approximation");
title!("AIC = $(round(AIC; digits=3))");

Q = plot(; legend=:bottomright, xtick=a:b, ytick=0:0.1:1);
plot!(xs, Ys; linealpha = 0, fillrange = Zs, fillcolor = :yellow,
    label="area = W_1 = Wasserstein distance for p = 1");
plot!(xs, Ys; lw = 2, color = 1,
    label="ecdf of size-$n sample of Poisson($λ)");
plot!(xs, Zs; lw = 2, color = 2,
    label="cdf of normal distribution approximation")
title!("W_1 = $(round(W_1; digits=3))")

plot(P, Q; size=(600, 800), layout=(2, 1))

# %%
#\using StatsFuns
AIC = -2loglikelihood(normal, sample) + 2
#S = -2sum((x = count(==(k), sample)/length(sample); xlogy(x, x)) for k in 0:maximum(sample))*length(sample)
S = -2sum((x = pdf(dist_true, k); xlogy(x, x)) for k in 0:20) * n
@show AIC S AIC-S;

# %%
KL = -2sum((x = pdf(dist_true, k); y = pdf(normal, k); xlogy(x, y/x)) for k in 0:20) * n

# %%

max_stars_repo_path: 0017/Comparing discrete and continuous distributions.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Test
include("solution.jl")

input = [
  0 0 1 0 0;
  1 1 1 1 0;
  1 0 1 1 0;
  1 0 1 1 1;
  1 0 1 0 1;
  0 1 1 1 1;
  0 0 1 1 1;
  1 1 1 0 0;
  1 0 0 0 0;
  1 1 0 0 1;
  0 0 0 1 0;
  0 1 0 1 0;
]

@test power_consumption(input) == 198
@test life_support_rating(input) == 230

ext: jl


text:
# The plant node
# This file is part of the openBuildNet framework developed at EPFL.
# (C) 2016 LA-EPFL
# Author: Truong X. Nghiem (truong.nghiem@gmail.com, xuan.nghiem@epfl.ch)

# See: http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Examples.StandardMPC
# The plant model:
# x(t+1) = A*x(t) + B*u(t)
# y(t) = x(t)

using openbuildnet

# Model data
const A = [2. -1.;1. 0.2]
const B = [1. 0.]'
const x0 = [3., 1.];  # Initial state

x = x0  # The plant's state
xs = Vector{typeof(x)}()  # To store the state trajectory
push!(xs, x)

# Create the node
node = OBNNode("plant", "simplempc")

# Create the ports
u = create_input(node, "u", Vector{Float64})
y = create_output(node, "y", Vector{Float64})

# Create Callbacks
on_init(node) do
  println("INIT")
end

on_term(node) do
  println("TERMINATE")
end

const MAINBLOCK = 0

on_block_output(node, MAINBLOCK) do
  set(y, x)
end

on_block_state(node, MAINBLOCK) do
  global x
  x = A*x + B*get(u)
  push!(xs, x)  # save the state
end

println("Ready to start simulation...")
status = runsafe(node, 20)
println("Final status: ", status)

# Plot the results
if status == 2
  using PyPlot
  xsmat = hcat(xs...)' # convert to matrix
  plot(xsmat[:,1], xsmat[:,2])
  show()
end


text:
@with_kw struct CarterMethodBL{T} <: AbstractFirstOrderMetricParams{T}
    @deftype T
    M = 1.0
    a = 0.0
    E = 1.0
end

function carter_velocity(u, E, M, a, p)
    let r = u[2], θ = u[3], L = p.L, Q = p.Q
        Σ₀ = CarterMethodBLCoords.Σ(r, a, θ)
        (
            CarterMethodBLCoords.Σδt_δλ(E, L, M, r, a, θ) / Σ₀,
            p.r * CarterMethodBLCoords.Σδr_δλ(E, L, M, Q, r, a) / Σ₀,
            p.θ * CarterMethodBLCoords.Σδθ_δλ(E, L, Q, a, θ) / Σ₀,
            CarterMethodBLCoords.Σδϕ_δλ(E, L, M, r, a, θ) / Σ₀,
        )
    end
end

calc_lq(m::CarterMethodBL{T}, pos, vel) where {T} =
    CarterMethodBLCoords.LQ(m.M, pos[2], m.a, pos[3], vel[3], vel[4])

inner_radius(m::CarterMethodBL{T}) where {T} = m.M + √(m.M^2 - m.a^2)
constrain(::CarterMethodBL{T}, u, v; μ = 0.0) where {T} = v[1]

function integrator_problem(
    m::CarterMethodBL{T},
    pos::StaticVector{S,T},
    vel::StaticVector{S,T},
    time_domain,
) where {S,T}
    L, Q = calc_lq(m, pos, vel)
    ODEProblem{false}(pos, time_domain, make_parameters(L, Q, vel[2], T)) do u, p, λ
        SVector(carter_velocity(u, m.E, m.M, m.a, p)...)
    end
end

function integrator_problem(
    m::CarterMethodBL{T},
    pos::AbstractVector{T},
    vel::AbstractVector{T},
    time_domain,
) where {T}
    L, Q = calc_lq(m, pos, vel)
    ODEProblem{true}(pos, time_domain, make_parameters(L, Q, vel[2], T)) do du, u, p, λ
        du .= carter_velocity(u, m.E, m.M, m.a, p)
    end
end

function alpha_beta_to_vel(m::CarterMethodBL{T}, u, α, β) where {T}
    sinΦ, sinΨ = CarterMethodBLCoords.sinΦsinΨ(m.M, u[2], m.a, u[3], α, β)
    (β < 0.0 ? -1.0 : 1.0, sinΦ, sinΨ)
end

function Vr(m::CarterMethodBL{T}, u, p) where {T}
    L, Q, _, _ = p
    CarterMethodBLCoords.Vr(m.E, L, m.M, Q, u[2], m.a)
end
function Vθ(m::CarterMethodBL{T}, u, p) where {T}
    L, Q, _, _ = p
    CarterMethodBLCoords.Vθ(m.E, L, Q, m.a, u[3])
end

export CarterMethodBL


text:
module WeibullCount

export weibull_count_pdf,
    weibull_count_cdf,
    WeibullCountModel,
    pdf, cdf

using SpecialFunctions
using LRUCache

N_TERMS = 500

function loggamma(x)
    log(gamma(x))
end

# By using parametric types, we can use Float32 or BigFloat, and we will return the same type
"""
    `van_wijngaarden_transformation(a::AbstractArray{T})::T where T <: AbstractFloat`

Computes the alternating sum of the terms `a` (assuming first term to be have positive coefficient).

This method is an example ofa "series acceleration", and will in most cases increase the convergence rate, i.e. requiring fewer terms
to obtain converge for a given tolerance, as compared to just computing the partial sum.
Van-Wijngaarden transformation is a particular implementation of the Euler transformation of the alternating series given by terms `a`.
"""
function van_wijngaarden_transformation(a::AbstractArray{T})::T where T <: AbstractFloat
    n = length(a)

    first_row::Array{T} = zeros(n)
    first_row[1] = a[1]

    for i = 2:n
        first_row[i] = first_row[i - 1] + (-1)^(i - 1) * a[i]
    end

    row = first_row

    for i = 1:Int(floor(2 * n / 3))
        row = 0.5 * (first_row[1:end - 1] + first_row[2:end])
    end

    return row[end]
end

function _weibull_alpha(j::Integer, i::Integer, c::Number)::Number
    if i == 0
        return exp(loggamma(c * j + 1) - loggamma(j + 1))
    end

    return sum(exp(log(weibull_alpha(m, i - 1, c)) + loggamma(c * j - c * m + 1) - loggamma(j - m + 1))
    for m = i - 1:j - 1)
end

const WEIBULL_ALPHA_CACHE_SIZE = 10000
const _weibull_alpha_cache = LRU{Tuple{Integer, Integer, Number}, Number}(WEIBULL_ALPHA_CACHE_SIZE)

cache_clear!() = empty!(_weibull_alpha_cache)
cache_size!(n) = resize!(_weibull_alpha_cache, n)
cache_info() = length(_weibull_alpha_cache)

"""
    `weibull_alpha(j::Integer, i::Integer, c::Number)::Number`

This method used repeadly
"""
function weibull_alpha(j::Integer, i::Integer, c::Number)::Number
    @get! _weibull_alpha_cache (j, i, c) _weibull_alpha(j, i, c)
end

# HACK: only have cache for the highest
# const _weibull_alpha_cache = LRU{Tuple{BigInt, BigInt, BigFloat}, BigFloat}(32000)
# function weibull_alpha(j::BigInt, i::BigInt, c::BigFloat)::BigFloat
#     @get! _weibull_alpha_cache (j, i, c) _weibull_alpha(j, i, c)
# end

"DEPRECATED."
function weibull_count_pdf_old(x::T1, λ::T2, c::T2, t::T2; k::T1=convert(T1, N_TERMS), transform::Bool=true)::T2 where {T1 <: Integer, T2 <: AbstractFloat}
    js = collect(x:x + k)

    terms::Array{T2} = exp.(js .* (log(λ) .+ c * log(t)) .+ log.(weibull_alpha.(js, x, c)) .- loggamma.(c .* js .+ 1))

    if transform
        s = van_wijngaarden_transformation(terms)

        # TODO: start with fewer terms, and then keep adding until we have "convergence"
    else
        s = sum((-1).^(x .+ js) .* terms)
    end

    return s
end

"""
    `weibull_count_pdf_approx(x::T1, λ::T2, c::T2, t::T2; k::T1 = convert(T1, $N_TERMS), transform::Bool = true, increment::T1 = convert(T1, 10), tol::T2 = convert(T2, 0.001))::Tuple{T2, Bool} where {T1 <: Integer, T2 <: AbstractFloat}`

Computes approximation of Eq. 11 in [1] using partial sum and, if `transform = true`, Van-Wijngaarden transformation.

## References
[1] Adrian, M., Bradlow, E., Fader, P., & McShane, B., Count models based on weibull interarrival times, CoRR, (),  (2013).

"""
function weibull_count_pdf_approx(x::T1, λ::T2, c::T2, t::T2; k::T1 = convert(T1, N_TERMS), transform::Bool = true, increment::T1 = convert(T1, 10), tol::T2 = convert(T2, 0.001))::Tuple{T2, Bool} where {T1 <: Integer, T2 <: AbstractFloat}
    # endpoint
    maximum = x + k

    js = collect(x:x + increment)

    terms::Array{T2} = exp.(js .* (log(λ) .+ c * log(t)) .+ log.(weibull_alpha.(js, x, c)) .- loggamma.(c .* js .+ 1))

    if transform
        s = van_wijngaarden_transformation(terms)
    else
        s = sum((-1).^(x .+ js) .* terms)
    end

    converged = false
    valid = false

    # incremently add more terms
    for i = increment:increment:maximum
        js = collect(x + i + 1: min(maximum, x + i + increment))
        append!(terms, exp.(js .* (log(λ) .+ c * log(t)) .+ log.(weibull_alpha.(js, x, c)) .- loggamma.(c .* js .+ 1)))

        # print(i)
        # print(js)

        # compute difference
        if transform
            s_new = van_wijngaarden_transformation(terms)
        else
            s_new = sum((-1).^(x .+ js) .* terms)
        end

        Δ = abs(s_new - s)

        # println("Delta: $Δ")
        # println("s: $s")

        converged = Δ ≤ tol
        valid = (s ≥ T2(0.0)) & (s ≤ T2(1.0))

        if converged & valid
            break
        end

        # update
        s = s_new
    end

    # if !converged
    #     println("failed to converge")
    # end
    # if !valid
    #     println("invalid probability")
    # end

    return s, (valid & converged)
end

struct DivergedError{T1, T2} <: Exception where {T1 <: Integer, T2 <: AbstractFloat}
    x::T1
    λ::T2
    c::T2
    t::T2
    k::T1
    transform::Bool
    increment::T1
    tol::T2
    result::T2
end

Base.showerror(io::IO, e::DivergedError) = print(io, "series approximation failed to converge; try increasing `k`:\nx = $(e.x), λ = $(e.λ), c = $(e.c), t = $(e.t), k = $(e.k), transform = $(e.transform), increment = $(e.increment), tol = $(e.tol), result = $(e.result)")

"""
    `weibull_count_pdf(x::T1, λ::T2, c::T2, t::T2; k::T1 = convert(T1, $N_TERMS), transform::Bool = true, increment::T1 = convert(T1, 10), tol::T2 = convert(T2, 0.001))::T2 where {T1 <: Integer, T2 <: AbstractFloat}`

Computes an approximation to `P(X(t) = x)` for time `t`, for a Weibull counting model with parameters `λ` and `c`.

The main computation is performed by `WeibullCount.weibull_count_pdf_approx`

## See also
`WeibullCount.weibull_count_pdf_approx`

"""
function weibull_count_pdf(x::T1, λ::T2, c::T2, t::T2; k::T1 = convert(T1, N_TERMS), transform::Bool = true, increment::T1 = convert(T1, 10), tol::T2 = convert(T2, 0.001))::T2 where {T1 <: Integer, T2 <: AbstractFloat}
    p, converged = weibull_count_pdf_approx(x, λ, c, t; k = k, transform = transform, increment = increment, tol = tol)

    if !converged
        throw(DivergedError(x, λ, c, t, k, transform, increment, tol, p))
    end

    return p
end

function weibull_count_pdf(x::T1, λ::T2, c::T2; k::T1 = convert(T1, N_TERMS), transform::Bool = true, increment::T1 = convert(T1, 10), tol::T2 = convert(T2, 0.001))::T2 where {T1 <: Integer, T2 <: AbstractFloat}
    weibull_count_pdf(x, λ, c, T2(1.0); k = k, transform = transform, increment = increment, tol = tol)
end

function weibull_count_cdf(x::T1, λ::T2, c::T2; k::T1=convert(T1, N_TERMS))::T2 where {T1 <: Integer, T2 <: AbstractFloat}
    if x < 0
        return 0
    elseif x == 0
        return sum(weibull_count_pdf(T1(0), λ, c; k=k))
    end

    return sum(weibull_count_pdf(y, λ, c; k=k) for y = 0:x)
end

# includes
include("weibull_count_distribution.jl")

end # module


text:
using Match

struct Instruction
    action::Char
    value::Int
end

mutable struct Ship
    direction::Complex{Int}
    position::Complex{Int}
end

function move!(s::Ship, i::Instruction)
    @match i.action begin
        'N' => (s.position += i.value * im)
        'S' => (s.position -= i.value * im)
        'E' => (s.position += i.value)
        'W' => (s.position -= i.value)
        'L' => (s.direction *= (im)^(i.value ÷ 90))
        'R' => (s.direction *= (-im)^(i.value ÷ 90))
        'F' => (s.position += s.direction * i.value)
        _ => error(i)
    end
end

function part1(lines)
    ship = Ship(1, 0)
    for line in lines
        inst = Instruction(line[1], parse(Int, line[2:end]))
        move!(ship, inst)
    end
    return abs(real(ship.position)) + abs(imag(ship.position))
end

const DATA = readlines("day12.txt")

println("Part 1: $(part1(DATA))")

function move2!(s::Ship, i::Instruction, waypoint::Complex{Int})
    @match i.action begin
        'N' => (waypoint += i.value * im)
        'S' => (waypoint -= i.value * im)
        'E' => (waypoint += i.value)
        'W' => (waypoint -= i.value)
        'L' => (waypoint *= (im)^(i.value ÷ 90))
        'R' => (waypoint *= (-im)^(i.value ÷ 90))
        'F' => (s.position += waypoint * i.value)
        _ => error(i)
    end
    return waypoint
end

function part2(lines)
    ship = Ship(1, 0)
    waypoint = complex(10, 1)
    for line in lines
        inst = Instruction(line[1], parse(Int, line[2:end]))
        waypoint = move2!(ship, inst, waypoint)
    end
    return abs(real(ship.position)) + abs(imag(ship.position))
end

println("Part 2: $(part2(DATA))")


text:
const vertex_data = Ref{Any}(nothing)
const edge_data = Ref{Any}(nothing)
const mode_data = Ref{Any}(nothing)
const time_data = Ref{Any}(nothing)

E() = sort!(collect(keys(edge_data[]))) #define edge set caller
V() = sort!(collect(keys(vertex_data[]))) #define vertex set caller
M() = sort!(collect(keys(mode_data[]))) #define mode set caller
T() = time_data[] #define time set caller

"""
    read_vertex(vertex_file)
accept strings of file names in .csv format to be processed into a dictionary of `satker` type. Also sets the `V()` constant.
"""
function read_vertex(vertex_file::String,demand_file::String)
    vertex = CSV.File(vertex_file)
    demand = CSV.File(demand_file)

    time_data[] = [1:1:length(demand);] #set time period
    vertex_dict = Dict{String,satker}()
    for r in vertex
        vertex_dict[r.name] = satker(
            r.x, r.y, r.MAX, r.MIN,
            r.START, getproperty(demand,Symbol(r.name))
        )
    end
    vertex_data[] = vertex_dict

    return nothing
end

"""
    read_edge(edge_list)
accept edge list files in .csv to turn into vector of 'lin' type. Sets the `E` constant.
"""
function read_edge(edge_list::String)
    trayek = CSV.File(edge_list)

    edge_dict = Dict{Int64,lin}()
    idx = 1
    for r in trayek
        edge_dict[idx] = lin(
            r.ori, r.dst,
            r.md, r.w
        )
        idx += 1
    end
    edge_data[] = edge_dict

    return nothing
end

"""
    read_mode(mode_file)
accept modes of transportation and its parameters in .csv and turns it into a dict of `mode` type. Also sets the `M` constant.
"""
function read_mode(mode_file::String)
    mode = CSV.File(mode_file)

    mode_dict = Dict{String,moda}()
    for r in mode
        mode_dict[r.name] = moda(
            r.Q, r.con, r.dis, r.var
        )
    end
    mode_data[] = mode_dict

    return nothing
end


text:
import Base: promote_rule

# SparseMatrixCSC
promote_rule(::Type{SparseMatrixCSC{Tv,Ti}}, ::Type{Matrix{T}}) where {Tv,Ti,T} =
    Matrix{promote_type(T, Tv)}

# IMatrix
promote_rule(
    ::Type{IMatrix{N,T}},
    ::Type{PermMatrix{Tv,Ti,Vv,Vi}},
) where {N,T,Tv,Ti,Vv,Vi} = (TT = promote_type(T, Tv); PermMatrix{TT,Ti,Vector{TT},Vi})
promote_rule(::Type{IMatrix{N,T}}, ::Type{SparseMatrixCSC{Tv,Ti}}) where {N,T,Tv,Ti} =
    SparseMatrixCSC{promote_type(T, Tv),Ti}
promote_rule(::Type{IMatrix{M,TA}}, ::Type{Matrix{TB}}) where {M,TA,TB} = Array{TB,2}

# PermMatrix
promote_rule(
    ::Type{PermMatrix{TvA,TiA}},
    ::Type{SparseMatrixCSC{TvB,TiB}},
) where {TvA,TiA,TvB,TiB} = SparseMatrixCSC{promote_type(TvA, TvB),promote_type(TiA, TiB)}
promote_rule(::Type{PermMatrix{Tv,Ti}}, ::Type{Matrix{T}}) where {Tv,Ti,T} =
    Array{promote_type(Tv, T),2}


text:
using Statistics
import Tables

include("reckoner_types.jl")

alpha(dist::Beta{Float64})::Float64 = params(dist)[1]
beta(dist::Beta{Float64})::Float64 = params(dist)[2]
update(prior::Beta, new_a::Real, new_b::Real)::Beta = Beta(alpha(prior) + new_a, beta(prior) + new_b)
update(left::Beta, right::Beta)::Beta = update(left, alpha(right), beta(right))
update(priors::Vector{<:Beta})::Beta = Beta(sum(alpha.(priors)), sum(beta.(priors)))
reflect(dist::Beta)::Beta = Beta(beta(dist), alpha(dist))

geom_mean(vals::Vector{Float64})::Float64 = exp(sum(log.(vals)) / length(vals))

geom_mean(dists::Vector{Beta{Float64}})::Beta{Float64} = Beta(geom_mean(alpha.(dists)), geom_mean(beta.(dists)))

scale(dist::Beta{Float64}, scale::Real)::Beta{Float64} = Beta(alpha(dist) * scale, beta(dist) * scale)

struct DefaultMatch <: AbstractMatch
    challenge::Beta{Float64}
    timestamp::Int64
    win::Bool
    team_id::Int16
end

function DefaultMatch(inrow)::DefaultMatch
    DefaultMatch(Beta(0.5, 0.5), inrow.timestamp, inrow.win, inrow.team_id)
end

challenge(match::DefaultMatch) = match.challenge
timestamp(match::DefaultMatch) = match.timestamp
win(match::DefaultMatch) = match.win
team_id(match::DefaultMatch) = match.team_id

struct DefaultMatches <: AbstractMatches
    challenge::Vector{Beta{Float64}}
    timestamp::Vector{Int64}
    win::Vector{Bool}
    win_chance::Vector{Float64}
end

challenge(matches::DefaultMatches) = matches.challenge
timestamp(matches::DefaultMatches) = matches.timestamp
win(matches::DefaultMatches) = matches.win
win_chance(matches::DefaultMatches) = matches.win_chance

function DefaultMatches(intable)::DefaultMatches
    columns = Tables.columns(intable)

    DefaultMatches(columns.challenge, columns.timestamp, columns.win, columns.rating_weight)
end

const RowDefaultMatches = Vector{DefaultMatch}

function default_aup(curr::DefaultMatch)::DefaultMatches
    DefaultMatches([Beta(1,1), Beta(1,1)], [curr.timestamp, curr.timestamp], [true, false], [1.0, 1.0])
end

function default_weight(curr::AbstractMatch, prev)::Float64

    function time_penalty(timestamp_1::Int64, timestamp_2::Int64)::Float64
        # The time penalty is e^(rt) where r is -0.02 and t is in days
        # Thus, a game becomes worth 2% less towards a rank for every day.
        rate::Float64 = -0.02
        time::Float64 = (timestamp_1 - timestamp_2) / (24 * 60 * 60)
        penalty::Float64 = exp(rate * time)
    end

    if (curr.timestamp < prev.timestamp)
        # print(curr.timestamp, ", ", prev.timestamp, "\n")
        return 0.0 # future games do not count to your rank
    end

    weight = time_penalty(curr.timestamp, prev.timestamp)
end

function default_challenge_window(curr::AbstractMatch, prev)::Float64
    challenge_1::Beta{Float64} = challenge(curr)
    challenge_2::Beta{Float64} = challenge(prev)

    if sum(params(challenge_1)) < sum(params(challenge_2))
        penalty::Float64 = cdf(challenge_1, mean(challenge_2))
    else
        penalty = 1 - cdf(challenge_2, mean(challenge_2))
    end

    if !prev.win
        penalty = 1 - penalty
    end

    penalty
end

function default_skill(wins::Vector{Bool}, weights::Vector{<:Real}, challenge_windows::Vector{<:Real})::Beta{Float64}
    weights .*= challenge_windows

    a::Float64 = sum(weights[wins])
    b::Float64 = sum(weights[.!wins])

    if (isnan(a) || isnan(b)) print(weights, "\n") end

    if ((a <= 0) || (b <= 0)) print(a, ", ", b, "\n") end

    Beta(a, b)
end

function default_rating(wins::Vector{Bool}, weights::Vector{<:Real}, win_chances::Vector{<:Real})::Beta{Float64}
    win_chances[wins] .= 1 .- win_chances[wins]

    weights .*= win_chances

    a::Float64 = sum(weights[wins])
    b::Float64 = sum(weights[.!wins])

    if (isnan(a) || isnan(b)) print(weights, "\n") end

    if ((a <= 0) || (b <= 0)) print(a, ", ", b, "\n") end

    Beta(a, b)
end

function allocate_losses(skills::Vector{Beta{Float64}}, teams::Vector{<:Integer})::Dirichlet{Float64}
    n::Int64 = length(skills)

    totals::Vector{Float64} = zeros(n)

    for i in 1:n
        team_size::Integer = sum(teams .== teams[i])
        totals[teams .!= teams[i]] .+= (beta(skills[i]) / (n - team_size))
        totals[i] += alpha(skills[i])
    end

    Dirichlet(totals)
end

function default_eff_challenge(ratings::Vector{Beta{Float64}}, teams::Vector{<:Integer})::Vector{Beta{Float64}}
    # Effectively measures the strength of the opponents "minus" the strength of teammates
    n::Int64 = length(ratings)

    challenges::Vector{Beta{Float64}} = Vector{Beta{Float64}}(undef, n)

    av_sz::Float64 = n / length(unique(teams))

    for i in 1:n
        opp::BitArray = teams .!= teams[i]
        a_opp::Float64 = sum(alpha.(ratings[opp]))
        b_opp::Float64 = sum(beta.(ratings[opp]))

        team::BitArray = .!opp .& (1:n .!= i)
        a_ally::Float64 = sum(alpha.(ratings[team]))
        b_ally::Float64 = sum(beta.(ratings[team]))

        a_eff::Float64 = (a_opp + b_ally) * av_sz / (sum(team) + 1)
        b_eff::Float64 = (b_opp + a_ally) * av_sz / (sum(opp))

        if (!(a_eff > 0.0) || !(b_eff > 0.0))
            println(ratings)
            println(teams)
            println("opp: $a_opp, $b_opp")
            println("allies: $a_ally, $b_ally")
        end

        challenges[i] = Beta(a_eff, b_eff)
    end

    challenges
end

function default_win_chances(local_skills::Vector{Beta{Float64}}, teams::Vector{<:Integer})::Vector{Beta{Float64}}
    # We're actually constructing a Dirichlet distribution here and then taking the marginal distributions per team.
    # We've just skipped the part where we explictly call anything a Dirichlet distribution
    n::Int64 = maximum(teams)

    alphas::Vector{Float64} = Vector{Float64}(undef, n)

    for i in 1:n
        a_team::Float64 = sum(alpha.(local_skills[teams .== i]))

        b_opp::Float64 = sum(beta.(local_skills[teams .!= i])) / (n - 1)

        alphas[i] = a_team + b_opp
    end

    total::Float64 = sum(alphas)

    chances::Vector{Beta{Float64}} = Beta.(alphas, total .- alphas)
end

function default_rank_interval(skill::Beta{Float64})::Tuple{Float64, Float64}
    (quantile(skill, .25), quantile(skill, .75))
end

function default_display_rank(win_chance::Real)::Float64
    base_display::Int32 = 1500

    display_rank::Float64 = base_display / (1 - win_chance) - base_display

    threshold::Int32 = 3000

    if (display_rank > threshold)
        display_rank = log(display_rank) / log(threshold) * threshold
    end

    display_rank
end

# function default_display_rank(win_chance::Real, ratio::Real)::Float64
#     base_display::Int32 = 1000

#     display_rank::Float64 = base_display * (1 / (1 - win_chance) - ratio)

#     threshold::Int32 = 2000

#     if (display_rank > threshold)
#         display_rank = log(display_rank) / log(threshold) * threshold
#     end

#     display_rank
# end

# function bradley_terry_display(win_chance::Real)::Float64
#     1500 * (1/win_chance - 1)^-1
# end

# function elo_display(win_chance::Real)::Float64
#     1500 - 400(log10((1/win_chance) - 1))
# end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__); Pkg.instantiate()

include("prinz_neuron_setup.jl")
include("helpers.jl")

system = neuron.sys
@named simulation = ODESystem([D(system.Isyn) ~ 0]; systems = [system])
simplified = structural_simplify(simulation)

sim_exp = ModelingToolkit.build_torn_function(simplified; expression=true)
cleaned = clean_expr!(sim_exp, simplified)

clipboard(format_text(string(prettify(cleaned)))) # can paste readable output below

function prinz_neuron!(du, u, p, t)
    @inbounds let Caᵢ         = u[1],
                  Vₘ          = u[2],
                  NaV₊m       = u[3],
                  NaV₊h       = u[4],
                  CaT₊m       = u[5],
                  CaT₊h       = u[6],
                  CaS₊m       = u[7],
                  CaS₊h       = u[8],
                  KA₊m        = u[9],
                  KA₊h        = u[10],
                  KCa₊m       = u[11],
                  Kdr₊m       = u[12],
                  H₊m         = u[13],

                  cₘ          = p[1],
                  aₘ          = p[2],
                  τCa         = p[3],
                  Ca∞         = p[4],
                  f           = p[5],
                  ENa         = p[6],
                  EK          = p[7],
                  EK          = p[8],
                  EK          = p[9],
                  EH          = p[10],
                  El          = p[11],
                  NaV₊gbar    = p[12],
                  CaT₊gbar    = p[13],
                  CaS₊gbar    = p[14],
                  KA₊gbar     = p[15],
                  KCa₊gbar    = p[16],
                  Kdr₊gbar    = p[17],
                  H₊gbar      = p[18],
                  leak₊g      = p[19]

        let
            du[1] =
                (
                    Ca∞ +
                    -1Caᵢ +
                    -1000 *
                    f *
                    (
                        CaS₊gbar *
                        aₘ *
                        CaS₊h *
                        CaS₊m^3.0 *
                        (
                            Vₘ +
                            12.199306599503455 *
                            log(0.0003333333333333333 * max(Caᵢ, 0.001))
                        ) +
                        CaT₊gbar *
                        aₘ *
                        CaT₊h *
                        CaT₊m^3.0 *
                        (
                            Vₘ +
                            12.199306599503455 *
                            log(0.0003333333333333333 * max(Caᵢ, 0.001))
                        )
                    )
                ) * inv(τCa)
            du[2] =
                (
                    -1 * aₘ * leak₊g * (Vₘ + -1El) +
                    -1 *
                    H₊gbar *
                    aₘ *
                    H₊m *
                    (Vₘ + -1EH) +
                    -1 *
                    KCa₊gbar *
                    aₘ *
                    KCa₊m^4.0 *
                    (Vₘ + -1EK) +
                    -1 *
                    Kdr₊gbar *
                    aₘ *
                    Kdr₊m^4.0 *
                    (Vₘ + -1EK) +
                    -1 *
                    CaS₊gbar *
                    aₘ *
                    CaS₊h *
                    CaS₊m^3.0 *
                    (
                        Vₘ +
                        12.199306599503455 *
                        log(0.0003333333333333333 * max(Caᵢ, 0.001))
                    ) +
                    -1 *
                    CaT₊gbar *
                    aₘ *
                    CaT₊h *
                    CaT₊m^3.0 *
                    (
                        Vₘ +
                        12.199306599503455 *
                        log(0.0003333333333333333 * max(Caᵢ, 0.001))
                    ) +
                    -1 *
                    KA₊gbar *
                    aₘ *
                    KA₊h *
                    KA₊m^3.0 *
                    (Vₘ + -1EK) +
                    -1 *
                    NaV₊gbar *
                    aₘ *
                    NaV₊h *
                    NaV₊m^3.0 *
                    (Vₘ + -1ENa)
                ) *
                inv(aₘ) *
                inv(cₘ)
            du[3] =
                (
                    -1NaV₊m +
                    inv(1.0 + exp(-4.8204158790170135 + -0.1890359168241966Vₘ))
                ) * inv(2.64 + -2.52 * inv(1 + exp(-4.8 + -0.04Vₘ)))
            du[4] =
                0.7462686567164178 *
                (1.0 + exp(-6.29 + -0.1Vₘ)) *
                (
                    -1NaV₊h +
                    inv(1.0 + exp(9.44015444015444 + 0.19305019305019305Vₘ))
                ) *
                inv(1.5 + inv(1.0 + exp(9.694444444444445 + 0.2777777777777778Vₘ)))
            du[5] =
                (
                    -1CaT₊m +
                    inv(1.0 + exp(-3.7638888888888893 + -0.1388888888888889Vₘ))
                ) * inv(
                    43.4 +
                    -42.6 *
                    inv(1.0 + exp(-3.321951219512195 + -0.04878048780487805Vₘ)),
                )
            du[6] =
                (
                    -1CaT₊h +
                    inv(1.0 + exp(5.836363636363637 + 0.18181818181818182Vₘ))
                ) * inv(
                    210.0 +
                    -179.6 *
                    inv(1.0 + exp(-3.2544378698224854 + -0.0591715976331361Vₘ)),
                )
            du[7] =
                (
                    -1CaS₊m +
                    inv(1.0 + exp(-4.074074074074074 + -0.1234567901234568Vₘ))
                ) * inv(
                    2.8 +
                    14.0 * inv(
                        exp(-5.384615384615385 + -0.07692307692307693Vₘ) +
                        exp(2.7 + 0.1Vₘ),
                    ),
                )
            du[8] =
                (
                    -1CaS₊h +
                    inv(1.0 + exp(9.67741935483871 + 0.16129032258064516Vₘ))
                ) * inv(
                    120.0 +
                    300.0 * inv(
                        exp(-4.0625 + -0.0625Vₘ) +
                        exp(6.111111111111111 + 0.1111111111111111Vₘ),
                    ),
                )
            du[9] =
                (
                    -1KA₊m +
                    inv(1.0 + exp(-3.1264367816091956 + -0.1149425287356322Vₘ))
                ) * inv(
                    23.2 +
                    -20.8 *
                    inv(1.0 + exp(-2.164473684210526 + -0.06578947368421052Vₘ)),
                )
            du[10] =
                (
                    -1KA₊h +
                    inv(1.0 + exp(11.612244897959181 + 0.2040816326530612Vₘ))
                ) * inv(
                    77.2 +
                    -58.4 *
                    inv(1.0 + exp(-1.4679245283018867 + -0.03773584905660377Vₘ)),
                )
            du[11] =
                (
                    -1KCa₊m +
                    Caᵢ *
                    inv(3.0 + Caᵢ) *
                    inv(1.0 + exp(-2.246031746031746 + -0.07936507936507936Vₘ))
                ) * inv(
                    180.6 +
                    -150.2 *
                    inv(1.0 + exp(-2.026431718061674 + -0.04405286343612335Vₘ)),
                )
            du[12] =
                (
                    -1Kdr₊m +
                    inv(1.0 + exp(-1.0423728813559323 + -0.0847457627118644Vₘ))
                ) * inv(
                    14.4 +
                    -12.8 *
                    inv(1.0 + exp(-1.4739583333333335 + -0.052083333333333336Vₘ)),
                )
            du[13] =
                1 // 2 *
                (
                    exp(-1.8671328671328669 + 0.06993006993006992Vₘ) +
                    exp(-14.629310344827585 + -0.08620689655172414Vₘ)
                ) *
                (
                    -1H₊m +
                    inv(1.0 + exp(13.636363636363637 + 0.18181818181818182Vₘ))
                )
            nothing
        end
    end
end


text:
# determine if any of an ordered sequence are missing

isnothing(x::Nothing) = true
isnothing(x::T) where {T} = false

issomething(x::Nothing) = false
issomething(x::T) where {T} = true

isNaN(x::AbstractFloat) = isnan(x)
isNaN(x::T) where {T} = false

Base.findall(::Type{Missing}, x::AbstractArray{T,N}) where {N,T} = findmissings(x)
Base.findall(::Type{Nothing}, x::AbstractArray{T,N}) where {N,T} = findnothings(x)
Base.findlast(::Type{Missing}, x::AbstractArray{T,N}) where {N,T} = findlast(missing, x)
Base.findlast(::Type{Nothing}, x::AbstractArray{T,N}) where {N,T} = findlast(nothing, x)
Base.findfirst(::Type{Missing}, x::AbstractArray{T,N}) where {N,T} = findfirst(missing, x)
Base.findfirst(::Type{Nothing}, x::AbstractArray{T,N}) where {N,T} = findfirst(nothing, x)

for F in (:Float64, :Float32, :Float16)
  @eval begin
    Base.findall(::Type{$F}, x::AbstractArray{T,N}) where {N,T} = findNaNs(x)
    Base.findlast(::Type{$F}, x::AbstractArray{T,N}) where {N,T} = findNaNs(x)
    Base.findfirst(::Type{$F}, x::AbstractArray{T,N}) where {N,T} = findNaNs(x)
 end
end

findmissings(x::AbstractArray{T,N}) where {N,T} = findall(x .=== missing)
findnothings(x::AbstractArray{T,N}) where {N,T} = findall(x .=== nothing)
findNaNs(x::AbstractArray{T,N})     where {N,T} = findall(map(isNaN,x))

findpresents(x::AbstractArray{T,N})   where {N,T} = findall(x .!== missing)
findsomethings(x::AbstractArray{T,N}) where {N,T} = findall(x .!== nothing)
findnonNaNs(x::AbstractArray{T,N})    where {N,T} = findall(map(!isNaN,x))

for T in (:Int8, :Int16, :Int32, :Int64, :Int128,
          :UInt8, :UInt16, :UInt32, :UInt64, :UInt128,
          :Float16, :Float32, :Float64, :String, :Char)
  @eval begin
     removemissings(x::AbstractArray{$T,N}) where {N} = x
     removenothings(x::AbstractArray{$T,N}) where {N} = x
     function removemissings(x::AbstractArray{Union{Missing,$T},N}) where {N}
        result::$T = filter(!ismissing, x)
        return result
     end
     function removenothings(x::AbstractArray{Union{Nothing,$T},N}) where {N}
        result::$T = filter(!isnothing, x)
        return result
     end
     function removeNaNs(x::AbstractArray{Union{Nothing,$T},N}) where {N}
        result::$T = filter(!isNaN, x)
        return result
     end
  end
end


text:
export OperatorFunction

abstract OperatorFunction{BT,FF,T} <: Operator{T}

immutable ConcreteOperatorFunction{BT<:Operator,FF,T} <: OperatorFunction{BT,FF,T}
    op::BT
    f::FF
end

ConcreteOperatorFunction(op::Operator,f::Function) =
    ConcreteOperatorFunction{typeof(op),typeof(f),eltype(op)}(op,f)
OperatorFunction(op::Operator,f::Function) = ConcreteOperatorFunction(op,f)

for op in (:domainspace,:rangespace,:domain,:bandinds)
    @eval begin
        $op(OF::ConcreteOperatorFunction) = $op(OF.op)
    end
end

function getindex(OF::ConcreteOperatorFunction,k::Integer,j::Integer)
    @assert isdiag(OF.op)
    if k==j
        OF.f(OF.op[k,k])::eltype(OF)
    else
        zero(eltype(OF))
    end
end

function Base.convert{T}(::Type{Operator{T}},D::ConcreteOperatorFunction)
    if T==eltype(D)
        D
    else
        ConcreteOperatorFunction{typeof(D.op),T}(D.op,D.f)
    end
end

for OP in (:(Base.inv),:(Base.sqrt))
    @eval begin
        $OP(D::DiagonalOperator) = OperatorFunction(D,$OP)
        $OP(C::ConstantTimesOperator) = $OP(C.λ)*$OP(C.op)
        function $OP(D::ConcreteOperatorFunction)
            @assert isdiag(D)
            OperatorFunction(D.op,x->$OP(D.f(x)))
        end
        $OP(A::Operator) = isdiag(A) ? OperatorFunction(A,$OP) :
                                error("Not implemented.")
    end
end

Base.sqrt(S::SpaceOperator) = SpaceOperator(sqrt(S.op),S.domainspace,S.rangespace)
Base.inv(S::SpaceOperator) = SpaceOperator(inv(S.op),S.rangespace,S.domainspace)


text:
using MolecularGraph: sdfatom, sdfbond

@testset "sdfilereader" begin

@testset "sdfatom" begin
    full = sdfatom(
        "    1.1763    0.6815    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0"
    )
    @test full[1] == :C
    @test full[2] == 0
    @test full[3] == 1
    @test full[4] == [1.1763, 0.6815, 0.0]
    atomBr = sdfatom("   -0.1041    2.3896    0.0000 Br  0  0")
    @test atomBr[1] == :Br
    @test atomBr[4] == [-0.1041, 2.3896, 0.0]
    charged = sdfatom("    2.5488    1.2083    0.0000 O   0  3")
    @test charged[2] == 1
    radical = sdfatom("    1.6514    2.7627    0.0000 C   0  4")
    @test radical[2] == 0
    @test radical[3] == 2
end

@testset "sdfbond" begin
    @test sdfbond("  1  2  2  0  0  0  0") == (1, 2, 2, 0)
    @test sdfbond("  5  4  1  6  0  0  0") == (5, 4, 1, 6)
end

@testset "sdfmol" begin
    demomol = joinpath(dirname(@__FILE__), "..", "assets", "test", "demo.mol")
    mol = sdftomol(demomol)
    @test edgecount(mol) == 37
    @test nodecount(mol) == 37
    @test mol isa SDFile

    # attributes
    aspirinfile = joinpath(dirname(@__FILE__), "..", "assets", "test", "aspirin_3d.sdf")
    mol = sdftomol(aspirinfile)
    pcharges = mol.attributes[:PUBCHEM_MMFF94_PARTIAL_CHARGES]
    @test length(pcharges) == parse(Int, pcharges[1]) + 1
    for i = 2:length(pcharges)
        node, charge = split(pcharges[i])
        @test 1 <= parse(Int, node) <= nodecount(mol)
        @test isa(parse(Float32, charge), Float32)
    end
    @test mol.attributes[:PUBCHEM_COORDINATE_TYPE] == ["2", "5", "10"]  # ensure the last attribute is read
end

end # sdfilereader


text:
using MXNet
using Base.Test

# run test in the whole directory, latest modified files
# are run first, this makes waiting time shorter when writing
# or modifying unit-tests
function test_dir(dir)
  jl_files = sort(filter(x -> ismatch(r".*\.jl$", x), readdir(dir)), by = fn -> stat(joinpath(dir,fn)).mtime)
  map(reverse(jl_files)) do file
    include("$dir/$file")
  end
end

include(joinpath(dirname(@__FILE__), "common.jl"))
test_dir(joinpath(dirname(@__FILE__), "unittest"))

# run the basic MNIST mlp example
if haskey(ENV, "CONTINUOUS_INTEGRATION")
  include(joinpath(Pkg.dir("MXNet"), "examples", "mnist", "mlp-test.jl"))
end


text:
function ALback(; verbose=true)
    (activehistory.current == 1) && (replthrow("Nothing to go back to"); return)

    back!(activehistory)
    setreplset!(activelogicset(activehistory))
    verbose && println(activecommand(activehistory) * " - " * reportfeasible(replset))

end

function ALnext(; verbose=true)
    (activehistory.current == length(activehistory.logicsets)) && (replthrow("Nothing to go forward to"); return)

    next!(activehistory)
    setreplset!(activelogicset(activehistory))
    verbose &&  println(activecommand(activehistory) * " - " * reportfeasible(replset))

end


text:
# Internal dynamic model
function g(s_t_min::Vector{Float64}, u_t::Vector{Float64}) # Transition function
    ϕ_t_min = s_t_min[1]
    x_t_min = s_t_min[2]
    y_t_min = s_t_min[3]

    Δϕ_t = u_t[1]
    r_t = u_t[2]

    # Compute next state
    ϕ_t = ϕ_t_min + Δϕ_t
    x_t = x_t_min + r_t*cos(ϕ_t)
    y_t = y_t_min + r_t*sin(ϕ_t)

    s_t = [ϕ_t, x_t, y_t]

    return s_t
end

function g_inv_in1(s_t::Vector{Float64}, u_t::Vector{Float64}) # Inverse transition function towards previous state
    ϕ_t = s_t[1]
    x_t = s_t[2]
    y_t = s_t[3]

    Δϕ_t = u_t[1]
    r_t = u_t[2]

    # Compute previous state
    ϕ_t_min = ϕ_t - Δϕ_t
    x_t_min = x_t - r_t*cos(ϕ_t)
    y_t_min = y_t - r_t*sin(ϕ_t)

    s_t_min = [ϕ_t_min, x_t_min, y_t_min]

    return s_t_min
end

function g_inv_in2(s_t::Vector{Float64}, s_t_min::Vector{Float64}) # Inverse transition function towards control
    ϕ_t_min = s_t_min[1]
    x_t_min = s_t_min[2]
    y_t_min = s_t_min[3]

    ϕ_t = s_t[1]
    x_t = s_t[2]
    y_t = s_t[3]

    # Compute corresponding control
    Δϕ_t = ϕ_t - ϕ_t_min
    if Δϕ_t > 0.0
        Δϕ_t_tor = Δϕ_t - 2*pi
    else
        Δϕ_t_tor = Δϕ_t + 2*pi
    end
    if abs(Δϕ_t_tor) < abs(Δϕ_t)
        Δϕ_t = Δϕ_t_tor
    end
    Δx_t = x_t - x_t_min
    Δy_t = y_t - y_t_min
    r_t = sign(Δx_t*cos(ϕ_t) + Δy_t*sin(ϕ_t))*sqrt(Δx_t^2 + Δy_t^2)

    u_t = [Δϕ_t, r_t]

    return u_t
end

function initializeAgent()
    m_s_t_min = [ϕ_0, x_0, y_0] # Set initial brain state prior
    V_s_t_min = 0.01*diageye(3)

    # Initialize messages and marginals dictionary
    messages = Vector{Message}(undef, 59)
    marginals = Dict{Symbol, ProbabilityDistribution}()
    function infer(u_t::Vector{Float64}, o_t::Vector{Float64})
        x_t_c = m_s_t_min[2]
        y_t_c = m_s_t_min[3]
        θ_t = atan(y_t_c, x_t_c)
        if θ_t < 0
            m = [θ_t + pi, 0.0, 0.0] # Goal prior mean with orientation against polar angle
        else
            m = [θ_t - pi, 0.0, 0.0]
        end

        data = Dict(:o_t       => o_t,
                    :u_t       => u_t,
                    :m_o       => [m for k=2:T],
                    :V_o       => [Sigma for k=2:T],
                    :m_s_t_min => m_s_t_min,
                    :V_s_t_min => V_s_t_min)

        step!(data, marginals, messages)
    end

    marginals[:u_2] = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=zeros(2), v=tiny*diageye(2)) # Register initial action
    act() = mode(marginals[:u_2]) # Choose the mode of the current control state as action

    function slide(slide_msg_idx=7)
        (m_s_t_min, V_s_t_min) = ForneyLab.unsafeMeanCov(messages[slide_msg_idx].dist) # Reset prior state statistics
    end

    return (infer, act, slide)
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
abstract type AbstractCase end
export AbstractCase, EnumCase, SwitchCase, Leaf, Shaped

Continuation = Int
Branch = Pair{PatternInfo, Tuple{LineNumberNode, Continuation}}

"""for generating patterns with one-by-one checks
"""
struct EnumCase <: AbstractCase
    cases::Vector{AbstractCase}
end

"""for generating patterns with a optimized switch
"""
struct SwitchCase <: AbstractCase
    cases::Dict{TypeObject, AbstractCase}
end

"""specifying the code body of a case
"""
struct Leaf <: AbstractCase
    cont::Continuation
end

"""checking the shape of the datum with the defined patterns
"""
struct Shaped <: AbstractCase
    pattern::PatternInfo
    ln::LineNumberNode
    case::AbstractCase
end

@nospecialize

function case_split!(result::Vector{Branch}, branches::Vector{Branch})
    for (p, (ln, branch)) in branches
        if p.pattern isa Or
            case_split!(result, Branch[info => (ln, branch) for info in p.pattern.ps])
        else
            push!(result, p => (ln, branch))
        end
    end
end

function build_dyn(top::TypeObject, branches::Vector{Branch})::AbstractCase
    length(branches) === 1 && return begin
        br = branches[1]
        Shaped(br.first, br.second[1], Leaf(br.second[2]))
    end
    @assert !isempty(branches)
    enum_cases = AbstractCase[]
    labels = Vector{TypeObject}(undef, length(branches))
    groups = [1]
    for i in eachindex(branches)
        current_type = branches[i].first.typetag
        @assert current_type <: top
        if current_type === top
            # Detected forcely specified matching order.
            # e.g:
            #    match val with
            #    | (_: Int) -> ...
            #    | _ -> ...
            #    | (_ : String) -> ..
            # Although Int and String are orthogonal,
            # the second pattern `_` breaks their merging
            push!(groups, i)
            push!(groups, i + 1)
            labels[i] = top
            continue
        end

        last_group_start = groups[end]
        non_orthogonal_types = TypeObject[]
        non_orthogonal_indices = Int[]
        for j = last_group_start:i-1
            if typeintersect(current_type, labels[j]) !== Base.Bottom
                # Merging branches
                # e.g:
                #    match val with
                #   | (_ : Int) ->
                #   | (_ : String) ->
                #   | (_ : Float64) ->
                # We gives the 3 patterns such a label array:
                #   [Real, String, Real]
                # where, the 1st and 3rd branches merged
                push!(non_orthogonal_indices, j)
                push!(non_orthogonal_types, labels[j])
            end
        end
        merged_type = reduce(typejoin, non_orthogonal_types, init = current_type)
        for j in non_orthogonal_indices
            labels[j] = merged_type
        end
        labels[i] = merged_type
    end

    push!(groups, length(branches) + 1)
    n_groups = length(groups)
    # fix-point check
    if n_groups === 2 && all(labels .=== top)
        return EnumCase(
            AbstractCase[Shaped(pat, ln, Leaf(cont)) for (pat, (ln, cont)) in branches],
        )
    end
    for i_group = 1:n_groups-1
        start = groups[i_group]
        final = groups[i_group+1] - 1

        if start === final
            br = branches[start]
            push!(enum_cases, Shaped(br.first, br.second[1], Leaf(br.second[2])))
            continue
        elseif start > final
            continue
        end

        switch_map = Dict{TypeObject, Vector{Branch}}()
        for i = start:final
            vec = get!(switch_map, labels[i]) do
                Branch[]
            end
            push!(vec, branches[i])
        end
        switch = Pair{TypeObject, AbstractCase}[
            top′ => build_dyn(top′, branches′) for (top′, branches′) in switch_map
        ]
        push!(enum_cases, SwitchCase(Dict(switch)))
    end

    if length(enum_cases) === 1
        enum_cases[1]
    else
        EnumCase(enum_cases)
    end
end

function case_merge(branches::Vector{Branch})
    top = reduce(typejoin, TypeObject[pattern.typetag for (pattern, _) in branches])
    build_dyn(top, branches)
end

@specialize


text:
## Pattern matching needs a new implementation. The ad hoc approach can't go further.
## This is a start ...

### Types for patterns

abstract type AbstractBlanks end

abstract type ABlank  <: AbstractBlanks end
abstract type ABlankSequence <: AbstractBlanks end
abstract type ABlankNullSequence <: AbstractBlanks end

mutable struct BlankNoHead <: ABlank end
mutable struct BlankWithHead{T} <: ABlank
    head::T
end

mutable struct BlankSequenceNoHead <: ABlankSequence end
mutable struct BlankSequenceWithHead{T} <: ABlankSequence
    head::T
end

mutable struct BlankNullSequenceNoHead <: ABlankNullSequence end
mutable struct BlankNullSequenceWithHead{T} <: ABlankNullSequence
    head::T
end

const BlanksNoHead = Union{BlankNoHead,BlankSequenceNoHead,BlankNullSequenceNoHead}
const BlanksWithHead = Union{BlankWithHead,BlankSequenceWithHead,BlankNullSequenceWithHead}

topattern(mx::Mxpr) = mxpra(topattern(mhead(mx)), mapmargs(x -> topattern(x), margs(mx)))
topattern(x) = x

## replace this later with BlankT. after removing the old BlankT
topattern(b::Mxpr{:Blank}) = isempty(b) ? BlankNoHead() : BlankWithHead(_make_blank_head(b,margs(b)...))
topattern(b::Mxpr{:BlankSequence}) = isempty(b) ? BlankSequenceNoHead() : BlankSequenceWithHead(_make_blank_head(b,margs(b)...))
topattern(b::Mxpr{:BlankNullSequence}) = isempty(b) ? BlankNullSequenceNoHead() :
          BlankNullSequenceWithHead(_make_blank_head(b,margs(b)...))

_make_blank_head(b,head,args...) = error("More than one argument in Blank.")
_make_blank_head(b,head::Symbol) = @isdefined(head) ? eval(head) : head
_make_blank_head(b,head) = head

#### Matching1

pmatch(b::BlanksNoHead,ex) = true
pmatch(b::BlankWithHead{T},ex) where {T<:DataType} = isa(ex,b.head)
pmatch(b::BlankSequenceWithHead{T},ex) where {T<:DataType} = isa(ex,b.head)
pmatch(b::BlankNullSequenceWithHead{T},ex) where {T<:DataType} = isa(ex,b.head)
pmatch(b::BlanksWithHead,ex) = isa(ex,Mxpr{b.head})
pmatch(b::AbstractBlanks,ex) = symerror("match: Can't match Blank of type ", typeof(b))

#### Matching

mutable struct Capt
    c::Dict{SJSym,Any}
end

Capt() = Capt(Dict{SJSym,Any}())

abstract type  AbstractMatchData end

mutable struct SearchRange
    start::Int
    stop::Int
end

SearchRange(r::UnitRange) = SearchRange(r.start, r.stop)

mutable struct MatchDataIn
    _expr::Mxpr
    _indr::SearchRange
    _pat ## pattern in general sense
    _capt::Capt
end

## Neither flat nor orderless
mutable struct MatchData <: AbstractMatchData
    d::MatchDataIn
end

## flat
mutable struct MatchDataF <: AbstractMatchData
    d::MatchDataIn
end

## orderless
mutable struct MatchDataO <: AbstractMatchData
    d::MatchDataIn
end

## flat and orderless
mutable struct MatchDataFO <: AbstractMatchData
    d::MatchDataIn
end

getexpr(d::AbstractMatchData) = d.d._expr
getcap(d::AbstractMatchData) = d.d._capt
getstartind(d::AbstractMatchData) = d.d._indr.start
getstopind(d::AbstractMatchData) = d.d._indr.stop
setind(d::AbstractMatchData,i::Integer) = d.d._indr.start = i
incrind(d::AbstractMatchData) = d.d._indr.start += 1
getpat(d::AbstractMatchData) = d.d._pat

function setsearchrange(d::AbstractMatchData, r::UnitRange)
    d.d._indr.start = r.start
    d.d._indr.stop = r.stop
end

startsearchind(d::AbstractMatchData) = getstartind(d) + 1
stopsearchind(d::AbstractMatchData) = getstopind(d) + 1
firstel(d::AbstractMatchData) = d[startsearchind(d)]

Base.getindex(d::AbstractMatchData,i::Integer) =  getexpr(d)[i]
Base.length(d::AbstractMatchData) = length(getexpr(d))

matchdata(expr::Mxpr, pat) = matchdata(expr,pat,Capt())

function matchdata(expr::Mxpr, pat, capt::Capt)
    if isFlat(expr)
        if isOrderless(expr)
            MatchDataFO(expr,pat,capt)
        else
            MatchDataF(expr,pat,capt)
        end
    elseif isOrderless(expr)
        MatchDataO(expr,pat,capt)
    else
        MatchData(expr,pat,capt)
    end
end

_matchdatain(expr,pat,capt) = MatchDataIn(expr,SearchRange(0:0),pat,capt)

MatchData(expr,pat,capt) = MatchData(_matchdatain(expr,pat,capt))
MatchDataF(expr,pat,capt) = MatchDataF(_matchdatain(expr,pat,capt))
MatchDataO(expr,pat,capt) = MatchDataO(_matchdatain(expr,pat,capt))
MatchDataFO(expr,pat,capt) = MatchDataFO(_matchdatain(expr,pat,capt))

####

abstract type AbstractMatchIndex end

#abstract MatchIndexCompound <: AbstractMatchIndex

mutable struct MatchIndexSingle <: AbstractMatchIndex
    i::Int
end

mutable struct MatchIndexAbstractArray{T<:AbstractArray} <: AbstractMatchIndex
    a::T
end

MatchIndex(a::AbstractArray) = MatchIndexAbstractArray(a)
MatchIndex(a::Int) = MatchIndexSingle(a)

getstart(mi::MatchIndexAbstractArray) = start(mi.a)
getstop(mi::MatchIndexAbstractArray) = lastindex(mi.a)

getstart(mi::MatchIndexSingle) = mi.i
getstop(mi::MatchIndexSingle) = mi.i

####

abstract type AbstractMatchType end

mutable struct MatchOne{T<:AbstractMatchIndex} <: AbstractMatchType
    a::T
end

MatchOne(a) = MatchOne(MatchIndex(a))

mutable struct MatchNone <: AbstractMatchType
end

mutable struct MatchFail <: AbstractMatchType
end

struct MatchSequence{T<:AbstractMatchIndex} <: AbstractMatchType
    a::T
end

MatchSequence(a) = MatchSequence(MatchIndex(a))

struct MatchNullSequence{T<:AbstractMatchIndex} <: AbstractMatchType
    a::T
end

MatchNullSequence(a) = MatchNullSequence(MatchIndex(a))

getstart(mr::AbstractMatchType) = getstart(mr.a)
getstop(mr::AbstractMatchType) = getstop(mr.a)

function setsearchrange(d::AbstractMatchData, mr::AbstractMatchType)
    d.d._indr.start = getstart(mr)
    d.d._indr.stop = getstop(mr)
end

function setsearchstop(d::AbstractMatchData, mr::AbstractMatchType)
    d.d._indr.stop = getstop(mr)
end

#### Matching 2

function matchrange(data::MatchData, pat::ABlank)
#    pmatch(pat,data[startsearchind(data)]) ?
    #       MatchOne(startsearchind(data)) : MatchFail()
    (j1,j) = _matchrange_blank(data,pat)
    return j1 == 0 ? MatchFail() : MatchOne(j1:j)
end

function _matchrange_blank(data,pat)
    local j1
    matchedone::Bool = false
    for i in startsearchind(data):stopsearchind(data)
        j1 = i
        if pmatch(pat,data[i])
            matchedone = true
            break
        end
    end
    matchedone || return (0,0)
    j = j1
    #    for i in j1+1:length(data)
    maxi = min(length(data),stopsearchind(data))
#    for i in j1+1:stopsearchind(data)
    for i in j1+1:maxi
        pmatch(pat,data[i]) || break
        j = i
    end
    return (j1,j)
end

function matchrange(data::MatchData, pat::ABlankSequence)
    (j1,j) = _matchrange_blanksequence(data,pat)
    return j1 == 0 ? MatchFail() : MatchSequence(j1:j)
end

function matchrange(data::MatchData, pat::ABlankNullSequence)
    (j1,j) = _matchrange_blanksequence(data,pat)
    return j1 == 0 ? MatchNone() : MatchNullSequence(j1:j)
end

function _matchrange_blanksequence(data,pat)
    local j1
    matchedone::Bool = false
    for i in startsearchind(data):stopsearchind(data)
        j1 = i
        if pmatch(pat,data[i])
            matchedone = true
            break
        end
    end
    matchedone || return (0,0)
    j = j1
    for i in j1+1:length(data)
#    for i in j1+1:stopsearchind(data)
        pmatch(pat,data[i]) || break
        j = i
    end
    return (j1,j)
end

function advanceind(d::AbstractMatchData, mr::AbstractMatchType)
    nothing
end

function advanceind(d::AbstractMatchData, mr::MatchOne)
    setsearchrange(d,mr)
end

function advanceind(d::AbstractMatchData, mr::MatchSequence)
    setsearchrange(d,mr)
end

function advanceind(d::AbstractMatchData, mr::MatchNullSequence)
    setsearchstop(d,mr)
end

function matchranges(expr::Mxpr, pat)
    pat = topattern(pat)
    matchranges(matchdata(expr,pat))
end

function matchranges(d::AbstractMatchData)
    pat = getpat(d)
    a = Any[]
    for p in pat
        mr = matchrange(d,p)
        advanceind(d,mr)
        push!(a,mr)
    end
    a
end


text:
"""
    CompGraph
    CompGraph(input::AbstractVertex, output::AbstractVertex)
    CompGraph(input::AbstractVector{<:AbstractVertex}, output::AbstractVertex)
    CompGraph(input::AbstractVertex, output::AbstractVector{<:AbstractVertex})

Basic graph for computation. While not strictly neccessary to compute anything,
it makes it easier to keep track of things.

# Examples
```jldoctest
julia> using NaiveNASlib

julia> v1 = inputvertex("in1", 1) + inputvertex("in2", 1);

julia> v2 = invariantvertex(x -> 3x, v1);

julia> CompGraph(inputs(v1), v2)(2,3) # (2 + 3) * 3
15

julia> CompGraph(inputs(v1), [v1, v2])(2,3)
(5, 15)
```

"""
struct CompGraph{I<:AbstractVector{<:AbstractVertex}, O<:AbstractVector{<:AbstractVertex}}
    inputs::I
    outputs::O
end
CompGraph(input::AbstractVertex, output::AbstractVertex) = CompGraph([input], [output])
CompGraph(input::AbstractVector{<:AbstractVertex}, output::AbstractVertex) = CompGraph(input, [output])
CompGraph(input::AbstractVertex, output::AbstractVector{<:AbstractVertex}) = CompGraph([input], output)

@functor CompGraph

function (g::CompGraph)(x...)
    @assert length(x) == length(g.inputs) "Must supply one input for each input vertex!"
    memo = Dict{AbstractVertex, Any}(zip(inputs(g), x))
    if length(outputs(g)) == 1
        return output!(memo, first(outputs(g)))
    end
    return Tuple(map(v -> output!(memo, v), outputs(g)))
end

"""
    inputs(g::CompGraph)

Return the inputs vertices of `g`.
"""
inputs(g::CompGraph) = g.inputs

"""
    outputs(g::CompGraph)

Return the output vertices of `g`.
"""
outputs(g::CompGraph) = g.outputs

"""
    output!(memo::AbstractDict{K, V}, v::AbstractVertex) where {K,V}

Return the output from v given any input in memo by traversing the graph.
Intermediate results from all visited vertices will be stored in memo after
function exits.

# Examples
```jldoctest
julia> using NaiveNASlib, NaiveNASlib.Advanced, NaiveNASlib.Extend

julia> ivs = InputVertex.(1:2);

julia> v1 = CompVertex(*, ivs...);

julia> v2 = CompVertex(-, v1, ivs[1]);

julia> results = Dict{AbstractVertex, Any}(zip(ivs, [2,3]));

julia> output!(results, v2)
4
julia> Pair{AbstractVertex, Int}[v=>results[v] for v in ancestors(v2)]
4-element Vector{Pair{AbstractVertex, Int64}}:
                                         InputVertex(1) => 2
                                         InputVertex(2) => 3
 CompVertex(*, inputs=[InputVertex(1), InputVertex(2)]) => 6
  CompVertex(-, inputs=[CompVertex(*), InputVertex(1)]) => 4
```
"""
function output!(memo::AbstractDict{K,V}, v::AbstractVertex) where {K,V}
    # Calculate outputs which are not already calculated
    return get!(memo, v) do
        inpt = map(iv -> output!(memo, iv), inputs(v))
        v(inpt...)
    end::V
end

"""
    nvertices(g::CompGraph)

Return the number of vertices in the graph.
"""
nvertices(g::CompGraph) = length(vertices(g))

Base.getindex(g::CompGraph, args...) = getindex(vertices(g), args...)
Base.firstindex(g::CompGraph) = firstindex(vertices(g))
Base.lastindex(g) = lastindex(vertices(g))

"""
    findvertices(vname::AbstractString, g::CompGraph)

Return all vertices for which [`name(v)`](@ref) == vname`.
"""
findvertices(vname::AbstractString, g::CompGraph) = findvertices(v -> name(v) == vname, g)
"""
    findvertices(vpat::Regex, g::CompGraph)

Return all vertices for which `vpat` matches [`name(v)`](@ref).
"""
findvertices(vpat::Regex, g::CompGraph) = findvertices(v -> occursin(vpat, name(v)), g)

"""
    findvertices(predicate, g::CompGraph)

Return all vertices for which `predicate(v)` return `true`.
"""
findvertices(predicate, g::CompGraph) = filter(predicate, vertices(g))

"""
    ancestors(v::AbstractVertex)

Return an array of all ancestors of `v`, including `v` itself.

# Examples
```jldoctest
julia> using NaiveNASlib, NaiveNASlib.Advanced, NaiveNASlib.Extend

julia> ancestors(invariantvertex(+, inputvertex("in", 1)))
2-element Vector{AbstractVertex}:
 InputSizeVertex(InputVertex(in, outputs=[CompVertex(+)]), 1)
 MutationVertex(CompVertex(+, inputs=[in], outputs=[]), SizeInvariant())
```
"""
ancestors(v::AbstractVertex,args...) = collect_vertices_from(inputs, v, args...)

"""
    descendants(v::AbstractVertex)

Return an array of all descendants of `v`, including `v` itself.

# Examples
```jldoctest
julia> using NaiveNASlib, NaiveNASlib.Advanced, NaiveNASlib.Extend

julia> descendants(invariantvertex(+, inputvertex("in", 1)) |> inputs |> first)
2-element Vector{AbstractVertex}:
 MutationVertex(CompVertex(+, inputs=[in], outputs=[]), SizeInvariant())
 InputSizeVertex(InputVertex(in, outputs=[CompVertex(+)]), 1)
```
"""
descendants(v::AbstractVertex,args...) = collect_vertices_from(outputs, v, args...)

function collect_vertices_from(f, v::AbstractVertex, vertices::Vector{AbstractVertex} = Vector{AbstractVertex}(), visited::Vector{AbstractVertex} = Vector{AbstractVertex}())
    v in vertices && return vertices
    if !(v in visited)
        push!(visited, v)
        fvs = f(v)
        fvs === nothing && return vertices
        foreach(fv -> collect_vertices_from(f, fv, vertices), fvs)
    end
    push!(vertices, v)
    return vertices
end

"""
    vertices(g::CompGraph)

Return an topologically sorted array of all vertices in the graph `g`.

# Examples
```julia-repl
julia> ins = InputVertex.(1:2);

julia> v1 = CompVertex(+, ins...);

julia> v2 = CompVertex(*, v1, ins[2]);

julia> graph = CompGraph(ins, v2);

julia> vertices(graph)
4-element Array{AbstractVertex,1}:
 InputVertex(1)
 InputVertex(2)
 CompVertex(+), inputs=[InputVertex(1), InputVertex(2)]
 CompVertex(*), inputs=[CompVertex(+), InputVertex(2)]
```
"""
vertices(g::CompGraph) = unique(mapfoldl(ancestors, vcat, g.outputs))


text:
# i include this just incase
using Revise

using GeodesicRendering
using AccretionGeometry
using CarterBoyerLindquist
using AccretionFormulae
using StaticArrays

using Plots
gr()

m = CarterMethodBL(M=1.0, a=1.0)

risco = AccretionFormulae.r_isco(m.a, m.M)
d = GeometricThinDisc(risco, 50.0, deg2rad(90))

u = @SVector [0.0, 1000.0, deg2rad(30), 0.0]

cache = @time prerendergeodesics(
    m,
    u,
    2000.0,
    d,
    fov_factor = 9.0 * 4,
    abstol = 1e-9,
    reltol = 1e-9,
    image_width = 350 * 4,
    image_height = 250 * 4
)

# define a common filter
filter = ConstPointFunctions.filter_intersected

function flux(m, gp, max_time)
    g = AccretionFormulae.redshift(m, gp, max_time)
    r = gp.u[2]
    g^4 / r^3
end

flux_vf = (
    PointFunction(flux) ∘ filter
)

redshift_vf = (
    PointFunction(AccretionFormulae.redshift) ∘ filter
)

flux_img = GeodesicRendering.apply(flux_vf, cache)
redshift_img = GeodesicRendering.apply(redshift_vf, cache)

# select non-NaN
mask = @. ! isnan(redshift_img)
line_flux = flux_img[mask]
energy = (redshift_img .* 6.4)[mask]

# use `vec` to flatten the matrix into a 1d array
his = histogram(vec(energy), weights=vec(line_flux), legend=false)
yaxis!("Line Flux")
xaxis!("Energy")


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

name = "SDPA_GMP-QD-DD_Builder"
version = v"7.1.2"

# Collection of sources required to build SDPABuilder
sources = [
    "https://github.com/ericphanson/sdpa-gmp/archive/v7.1.3-patch-2.tar.gz" =>
    "198cf7449ae751911f40b501f4b9a4e87c5cefd0110442067e16d6c212d996ea",
    "https://gmplib.org/download/gmp/gmp-6.1.2.tar.bz2" =>
    "5275bb04f4863a13516b2f39392ac5e272f5e1bb8057b18aec1c9b79d73d8fb2",
    "https://github.com/ericphanson/sdpa-qd/archive/qd-v7.1.2-patch-1.tar.gz" =>
    "592ba5f75b79ab8ff36cd39b37b446800bdc36de349efb0274bbbb01d5e596b6",
    "https://github.com/ericphanson/sdpa-dd/archive/dd-7.1.2-patch-1.tar.gz" =>
    "0cbb2fd7c7e583553c771afacc36d1bcb561212f6a0a1e482e8eb3c8572aa125",
]

# Bash recipe for building across all platforms
script = raw"""
# Build GMP
cd $WORKSPACE/srcdir/gmp-*

update_configure_scripts

if [ $target = "x86_64-apple-darwin14" ]; then
  # seems static linking requires apple's ar
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

flags=(--enable-cxx --disable-shared --enable-static --with-pic)

# On x86_64 architectures, build fat binary
if [[ ${proc_family} == intel ]]; then
    flags+=(--enable-fat)
fi

./configure --prefix=$prefix --build=x86_64-linux-gnu --host=$target  ${flags[@]}
make
make install

# Build SDPA-GMP

cd $WORKSPACE/srcdir/sdpa-gmp-*/

update_configure_scripts

if [ $target = "x86_64-apple-darwin14" ]; then
  # seems static linking requires apple's ar
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

mv configure.in configure.ac
autoreconf -i

CXXFLAGS="-std=c++03"; export CXXFLAGS

./configure --prefix=$prefix --with-gmp-includedir=$prefix/include --with-gmp-libdir=$prefix/lib --host=$target lt_cv_deplibs_check_method=pass_all

make

cp sdpa_gmp $prefix/bin/sdpa_gmp
cp COPYING $prefix/bin/COPYING

# Build SDPA-QD

cd $WORKSPACE/srcdir/sdpa-qd-*/

update_configure_scripts

if [ $target = "x86_64-apple-darwin14" ]; then
  # seems static linking requires apple's ar
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

mv configure.in configure.ac
autoreconf -i

CXXFLAGS="-std=c++03"; export CXXFLAGS

./configure --prefix=$prefix --with-qd-includedir=$prefix/include --with-qd-libdir=$prefix/lib --host=$target lt_cv_deplibs_check_method=pass_all

make

cp sdpa_qd $prefix/bin/sdpa_qd

# Build SDPA-QD

cd $WORKSPACE/srcdir/sdpa-dd-*/

update_configure_scripts

if [ $target = "x86_64-apple-darwin14" ]; then
  # seems static linking requires apple's ar
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

mv configure.in configure.ac
autoreconf -i

CXXFLAGS="-std=c++03"; export CXXFLAGS

./configure --prefix=$prefix --with-qd-includedir=$prefix/include --with-qd-libdir=$prefix/lib --host=$target lt_cv_deplibs_check_method=pass_all

make

cp sdpa_dd $prefix/bin/sdpa_dd
"""

platforms = Platform[
   MacOS(:x86_64),
   Linux(:x86_64),
   # Windows(:x86_64), # doesn't work :(
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "sdpa_gmp", :sdpa_gmp),
    ExecutableProduct(prefix, "sdpa_qd", :sdpa_qd),
    ExecutableProduct(prefix, "sdpa_dd", :sdpa_dd),
]

# Dependencies that must be installed before this package can be built
dependencies = [ "https://github.com/ericphanson/QD_Builder/releases/download/v2.3.22/build_QD.v2.3.22.jl" ]

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


text:
##############################
# Turing Interop
##############################

using Turing.DynamicPPL

function Turing.DynamicPPL.Sampler(
    alg::AnISHMC,
    model::Turing.Model,
    s::Turing.DynamicPPL.Selector=Turing.DynamicPPL.Selector()
)
    info = Dict{Symbol,Any}()
    initial_state = Turing.Inference.SamplerState(Turing.VarInfo(model))

    return Turing.DynamicPPL.Sampler(alg, info, s, initial_state)
end

# From AIS code in Turing PR.
"""
    gen_logjoint(v, model, spl)
Return the log joint density function corresponding to model.
"""
function gen_logjoint(vi, model, spl)
    # TODO: Would this be thread-safe if we would copy the sampler N times?
    vis = [deepcopy(vi) for _ in 1:Threads.nthreads()]
    mdls = [deepcopy(model) for _ in 1:Threads.nthreads()]
    spls = [deepcopy(spl) for _ in 1:Threads.nthreads()]

    function logjoint(z)::Float64
        tid = Threads.threadid()
        vis[tid][spls[tid]] = z
        mdls[tid](vis[tid], spls[tid])
        lj = getlogp(vis[tid])
        return lj
    end
    return logjoint
end

function gen_logjoint_grad(vi, model, spl)
    vis = [deepcopy(vi) for _ in 1:Threads.nthreads()]
    mdls = [deepcopy(model) for _ in 1:Threads.nthreads()]
    spls = [deepcopy(spl) for _ in 1:Threads.nthreads()]

    function logjoint_grad(x)
        tid = Threads.threadid()
        return Turing.Core.gradient_logp(
            x, vis[tid], mdls[tid], spls[tid]
        )
    end
    return logjoint_grad
end

"""
    gen_logprior(v, model, spl)
Return the log prior density function corresponding to model.
"""
function gen_logprior(vi, model, spl)
    vis = [deepcopy(vi) for _ in 1:Threads.nthreads()]
    mdls = [deepcopy(model) for _ in 1:Threads.nthreads()]

    function logprior(z)::Float64
        tid = Threads.threadid()
        vis[tid][Turing.SampleFromPrior()] = z
        mdls[tid](
            vis[tid],
            Turing.SampleFromPrior(),
            Turing.DynamicPPL.PriorContext()
        )
        lj = getlogp(vis[tid])
        return lj
    end
    return logprior
end

function gen_logprior_grad(vi, model, spl)
    vis = [deepcopy(vi) for _ in 1:Threads.nthreads()]
    mdls = [deepcopy(model) for _ in 1:Threads.nthreads()]
    spls = [deepcopy(spl) for _ in 1:Threads.nthreads()]

    function logprior_grad(x)
        tid = Threads.threadid()
        return Turing.Core.gradient_logp(
            x, vis[tid], mdls[tid], spls[tid], Turing.DynamicPPL.PriorContext()
        )
    end
    return logprior_grad
end

function sample_from_prior(rng, model)
    # TODO: Need to make use of rng.
    vi = Turing.VarInfo(model)

    # Extract a NamedTuple from VarInfo which has variable names as keys and
    # the sampled values as values.
    vns = Turing.DynamicPPL._getvns(vi, Turing.SampleFromPrior())
    value_tuple = _val_tuple(vi, vns)
    return value_tuple
end

function gen_prior_sample(model)
    spl = Turing.SampleFromPrior()
    function prior_sampling()
        vi = Turing.VarInfo(model)
        z = vi[spl]
        return z
    end

    return prior_sampling
end

"""
Returns a function which is the prior density.
"""
function make_log_prior_density(model, vi=nothing)
    if isnothing(vi)
        vi = Turing.VarInfo(model)
    end

    return function logprior_density(named_tuple)
        set_namedtuple!(vi, named_tuple)
        model(vi, Turing.SampleFromPrior(), Turing.PriorContext())
        return Turing.getlogp(vi)
    end
end

function make_log_joint_density(model, vi=nothing)
    if isnothing(vi)
        vi = Turing.VarInfo(model)
    end

    return function logjoint_density(named_tuple)
        set_namedtuple!(vi, named_tuple)
        model(vi)
        return Turing.getlogp(vi)
    end
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::T,
    alg::AnIS,
    num_samples::Int;
    kwargs...
) where {T <: Union{Turing.Model,AnISModel}}
    vi = Turing.VarInfo(model)
    anis = AnnealedISSampler(model, alg, vi)
    samples, diagnostics = ais_sample(rng, anis, num_samples; kwargs...)
    samples = make_turing_samples(samples, model)
    return to_mcmcchains(samples, diagnostics)
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::T,
    alg::AnISHMC,
    num_samples::Int;
    kwargs...
) where {T <: Union{Turing.Model,AnISModel}}
    spl = Turing.DynamicPPL.Sampler(alg, model)
    # Transform variables from X -> R
    #Turing.DynamicPPL.link!(spl.state.vi, spl)

    anis = AnnealedISSamplerHMC(model, alg, spl)
    samples, diagnostics = ais_sample(rng, anis, num_samples; kwargs...)
    #Turing.DynamicPPL.link!(spl.state.vi, spl)
    samples = make_turing_samples(samples, model, spl)
    return to_mcmcchains(samples, diagnostics)
end

function make_turing_samples(samples, model::Turing.Model)
    vi = Turing.VarInfo(model)
    return map(samples) do sample
        set_namedtuple!(vi, sample.params)
        nt = DynamicPPL.tonamedtuple(vi)
        (θ = nt, log_weight = sample.log_weight)
    end
end

function make_turing_samples(samples, model::Turing.Model, spl::Turing.DynamicPPL.Sampler)
    vi = spl.state.vi
    return map(samples) do sample
        vi[spl] = sample.params
        nt = DynamicPPL.tonamedtuple(vi)
        (θ = nt, log_weight = sample.log_weight)
    end
end

function to_mcmcchains(samples, diagnostics)
    # Part of this function is adapted from the generic Turing bundle_samples at
    # https://github.com/TuringLang/Turing.jl/blob/60563a64c51f2ea465f85d819344be00d0186d1b/src/inference/Inference.jl

    # Convert transitions to array format.
    # Also retrieve the variable names.
    nms, vals = _params_to_array(samples)

    # Get the values of the extra parameters in each transition.
    extra_params, extra_values = get_transition_extras(samples)

    # Extract names & construct param array.
    nms = [nms; extra_params]
    parray = hcat(vals, extra_values)

    # Calculate logevidence.
    le = logsumexp(map(s -> s.log_weight, samples)) - log(length(samples))

    info = Dict{Symbol,Any}(
        :ess => diagnostics[:ess],
    )
    if haskey(diagnostics, :intermediate_samples)
        info[:intermediate_samples] = diagnostics[:intermediate_samples]
        info[:intermediate_log_weights] = diagnostics[:intermediate_log_weights]
    end
    if haskey(diagnostics, :intermediate_acceptance_rate)
        info[:intermediate_acceptance_rate] = diagnostics[:intermediate_acceptance_rate]
        info[:acceptance_rate] = diagnostics[:acceptance_rate]
    end
    info = (;info...)

    # Conretize the array before giving it to MCMCChains.
    parray = MCMCChains.concretize(parray)

    return MCMCChains.Chains(
        parray,
        string.(nms),
        (internals = ["log_weight"],);
        evidence=le,
        info=info
    )
end

##############################
# Functions adapted from https://github.com/TuringLang/Turing.jl/blob/60563a64c51f2ea465f85d819344be00d0186d1b/src/inference/Inference.jl
##############################

"""
    getparams(t)
Return a named tuple of parameters.
"""
getparams(t) = t.θ

function _params_to_array(ts::Vector)
    names = Vector{String}()
    # Extract the parameter names and values from each transition.
    dicts = map(ts) do t
        nms, vs = flatten_namedtuple(getparams(t))
        for nm in nms
            if !(nm in names)
                push!(names, nm)
            end
        end
        # Convert the names and values to a single dictionary.
        return Dict(nms[j] => vs[j] for j in 1:length(vs))
    end
    # names = collect(names_set)
    vals = [get(dicts[i], key, missing) for i in eachindex(dicts),
        (j, key) in enumerate(names)]

    return names, vals
end

function flatten_namedtuple(nt::NamedTuple)
    names_vals = mapreduce(vcat, keys(nt)) do k
        v = nt[k]
        if length(v) == 1
            return [(String(k), v)]
        else
            return mapreduce(vcat, zip(v[1], v[2])) do (vnval, vn)
                return collect(FlattenIterator(vn, vnval))
            end
        end
    end
    return [vn[1] for vn in names_vals], [vn[2] for vn in names_vals]
end

additional_parameters(::Type{<:NamedTuple}) = [:log_weight]

function get_transition_extras(ts::AbstractVector)
    # Get the extra field names from the sampler state type.
    # This handles things like :lp or :weight.
    extra_params = additional_parameters(eltype(ts))

    # Get the values of the extra parameters.
    local extra_names
    all_vals = []

    # Iterate through each transition.
    for t in ts
        extra_names = Symbol[]
        vals = []

        # Iterate through each of the additional field names
        # in the struct.
        for p in extra_params
            # Check whether the field contains a NamedTuple,
            # in which case we need to iterate through each
            # key/value pair.
            prop = getproperty(t, p)
            if prop isa NamedTuple
                for (k, v) in pairs(prop)
                    push!(extra_names, Symbol(k))
                    push!(vals, v)
                end
            else
                push!(extra_names, Symbol(p))
                push!(vals, prop)
            end
        end
        push!(all_vals, vals)
    end

    # Convert the vector-of-vectors to a matrix.
    valmat = [all_vals[i][j] for i in 1:length(ts), j in 1:length(all_vals[1])]

    return extra_names, valmat
end

##############################
# Functions from https://github.com/TuringLang/Turing.jl/blob/master/src/inference/mh.jl
##############################

# Code from https://github.com/TuringLang/Turing.jl/blob/master/src/inference/mh.jl
function set_namedtuple!(vi, nt::NamedTuple)
    for (n, vals) in pairs(nt)
        vns = vi.metadata[n].vns
        nvns = length(vns)

        # if there is a single variable only
        if nvns == 1
            # assign the unpacked values
            if length(vals) == 1
                vi[vns[1]] = [vals[1];]
            # otherwise just assign the values
            else
                vi[vns[1]] = [vals;]
            end
        # if there are multiple variables
        elseif vals isa AbstractArray
            nvals = length(vals)
            # if values are provided as an array with a single element
            if nvals == 1
                # iterate over variables and unpacked values
                for (vn, val) in zip(vns, vals[1])
                    vi[vn] = [val;]
                end
            # otherwise number of variables and number of values have to be equal
            elseif nvals == nvns
                # iterate over variables and values
                for (vn, val) in zip(vns, vals)
                    vi[vn] = [val;]
                end
            else
                error("Cannot assign `NamedTuple` to `VarInfo`")
            end
        else
            error("Cannot assign `NamedTuple` to `VarInfo`")
        end
    end
end

# unpack a vector if possible
unvectorize(dists::AbstractVector) = length(dists) == 1 ? first(dists) : dists

# possibly unpack and reshape samples according to the prior distribution
reconstruct(dist::Distribution, val::AbstractVector) = Turing.DynamicPPL.reconstruct(dist, val)
function reconstruct(
    dist::AbstractVector{<:UnivariateDistribution},
    val::AbstractVector
)
    return val
end
function reconstruct(
    dist::AbstractVector{<:MultivariateDistribution},
    val::AbstractVector
)
    offset = 0
    return map(dist) do d
        n = length(d)
        newoffset = offset + n
        v = val[(offset + 1):newoffset]
        offset = newoffset
        return v
    end
end

@generated function _val_tuple(
    vi,
    vns::NamedTuple{names}
) where {names}
    isempty(names) === 0 && return :(NamedTuple())
    expr = Expr(:tuple)
    expr.args = Any[
        :($name = reconstruct(unvectorize(Turing.DynamicPPL.getdist.(Ref(vi), vns.$name)),
                              Turing.DynamicPPL.getval(vi, vns.$name)))
        for name in names]
    return expr
end

##############################
# Functions from https://github.com/TuringLang/Turing.jl/blob/515c4d5316d52619388557af8f850740ae310122/src/utilities/helper.jl
##############################

struct FlattenIterator{Tname, Tvalue}
    name::Tname
    value::Tvalue
end

Base.length(iter::FlattenIterator) = _length(iter.value)
_length(a::AbstractArray) = sum(_length, a)
_length(a::AbstractArray{<:Number}) = length(a)
_length(::Number) = 1

Base.eltype(iter::FlattenIterator{String}) = Tuple{String, _eltype(typeof(iter.value))}
_eltype(::Type{TA}) where {TA <: AbstractArray} = _eltype(eltype(TA))
_eltype(::Type{T}) where {T <: Number} = T

@inline function Base.iterate(iter::FlattenIterator{String, <:Number}, i = 1)
    i === 1 && return (iter.name, iter.value), 2
    return nothing
end
@inline function Base.iterate(
    iter::FlattenIterator{String, <:AbstractArray{<:Number}},
    ind = (1,),
)
    i = ind[1]
    i > length(iter.value) && return nothing
    name = getname(iter, i)
    return (name, iter.value[i]), (i+1,)
end
@inline function Base.iterate(
    iter::FlattenIterator{String, T},
    ind = startind(T),
) where {T <: AbstractArray}
    i = ind[1]
    i > length(iter.value) && return nothing
    name = getname(iter, i)
    local out
    while i <= length(iter.value)
        v = iter.value[i]
        out = iterate(FlattenIterator(name, v), Base.tail(ind))
        out !== nothing && break
        i += 1
    end
    if out === nothing
        return nothing
    else
        return out[1], (i, out[2]...)
    end
end

@inline startind(::Type{<:AbstractArray{T}}) where {T} = (1, startind(T)...)
@inline startind(::Type{<:Number}) = ()
@inline startind(::Type{<:Any}) = throw("Type not supported.")
@inline function getname(iter::FlattenIterator, i::Int)
    name = string(ind2sub(size(iter.value), i))
    name = replace(name, "(" => "[");
    name = replace(name, ",)" => "]");
    name = replace(name, ")" => "]");
    name = iter.name * name
    return name
end
@inline ind2sub(v, i) = Tuple(CartesianIndices(v)[i])

####
#### Compiler interface, i.e. tilde operators.
####
#### Taken from https://github.com/TuringLang/Turing.jl/blob/master/src/inference/hmc.jl
# function DynamicPPL.assume(
#     rng,
#     spl::Sampler{<:AnISHMC},
#     dist::Distribution,
#     vn::VarName,
#     vi,
# )
#     updategid!(vi, vn, spl)
#     r = vi[vn]
#     # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
#     # r
#     return r, Bijectors.logpdf_with_trans(dist, r, istrans(vi, vn))
# end

# function DynamicPPL.dot_assume(
#     rng,
#     spl::Sampler{<:AnISHMC},
#     dist::MultivariateDistribution,
#     vns::AbstractArray{<:VarName},
#     var::AbstractMatrix,
#     vi,
# )
#     @assert length(dist) == size(var, 1)
#     updategid!.(Ref(vi), vns, Ref(spl))
#     r = vi[vns]
#     var .= r
#     return var, sum(Bijectors.logpdf_with_trans(dist, r, istrans(vi, vns[1])))
# end
# function DynamicPPL.dot_assume(
#     rng,
#     spl::Sampler{<:AnISHMC},
#     dists::Union{Distribution, AbstractArray{<:Distribution}},
#     vns::AbstractArray{<:VarName},
#     var::AbstractArray,
#     vi,
# )
#     updategid!.(Ref(vi), vns, Ref(spl))
#     r = reshape(vi[vec(vns)], size(var))
#     var .= r
#     return var, sum(Bijectors.logpdf_with_trans.(dists, r, istrans(vi, vns[1])))
# end

# function DynamicPPL.observe(
#     spl::Sampler{<:AnISHMC},
#     d::Distribution,
#     value,
#     vi,
# )
#     return DynamicPPL.observe(SampleFromPrior(), d, value, vi)
# end

# function DynamicPPL.dot_observe(
#     spl::Sampler{<:AnISHMC},
#     ds::Union{Distribution, AbstractArray{<:Distribution}},
#     value::AbstractArray,
#     vi,
# )
#     return DynamicPPL.dot_observe(SampleFromPrior(), ds, value, vi)
# end


text:
# Gathers network information used for generating AC OPF samples

"Loads a PowerModels network dictionary from the net file and parses out the file dir"
function load_net(net::AbstractString, net_path, print_level=0)
	(print_level > 0) && println("Loading PowerModules.jl model data...")
	net = joinpath(net_path, net)  		# Join net_path and net if net_path is given
	net_path, net_name = splitdir(net)  # Get the net path in case it was given with the network
	net = PM.parse_file(net)  			# Load PowerModels data
	return net, net_path
end

"Returns the PowerModels.jl network data information needed to generate AC OPF samples"
function gather_net_info(net)
	net_name = net["name"]

	# Get system information for relaxed AC OPF runs
	net_ref = PM.build_ref(net)
	bus_ref = net_ref[:it][:pm][:nw][NW_DEFAULT][:bus]

	load_buses = sort([load["load_bus"] for load in values(net["load"])])
	gen_buses = sort([gen["gen_bus"] for gen in values(net["gen"])
                 if bus_ref[gen["gen_bus"]]["bus_type"] != 3])
    ext_bus = sort([gen["gen_bus"] for gen in values(net["gen"])
                 if bus_ref[gen["gen_bus"]]["bus_type"] == 3])
    gen_buses_w_ext = sort([ext_bus; gen_buses])
	num_loads = length(load_buses)
    num_gens = length(gen_buses)
    num_buses = length(net["bus"])
    num_ext = length(ext_bus)
	num_branches = length(net["branch"])

	pg_max = [gen["pmax"] for gen in values(net["gen"])]

	# Get base load
	load_ref = net_ref[:it][:pm][:nw][NW_DEFAULT][:load]

    base_pd = [load_ref[bus]["pd"] for bus in 1:length(load_ref)]
	base_qd = [load_ref[bus]["qd"] for bus in 1:length(load_ref)]
	base_load = vcat(base_pd, base_qd)

	keys_ = ("net_name", "load_buses", "gen_buses", "ext_bus", "num_loads", "num_gens", "num_buses",
			 "num_ext", "num_branches", "pg_max", "base_load")

	values_ = (net_name, load_buses, gen_buses, ext_bus, num_loads, num_gens, num_buses,
			   num_ext, num_branches, pg_max, base_load)

	net_info = Dict(zip(keys_, values_))
	return net_info
end

"Creates a copy of the given network dictionary and stores additional constraint information"
function create_net_r(net, pf_min, pf_lagging, base_load)
	net_r = copy(net)  # PowerModels data dict to use for convex formulations
	num_loads = length(net["load"])
	pf_min isa Number && (pf_min = ones(num_loads) .* pf_min)
	net_r["pf_min"] = pf_min  # Set pf_min for relaxed formulation load constraints
	net_r["pf_lagging"] = pf_lagging
	net_r["base_load"] = base_load
	return net_r
end

""
function find_max_loads(pd_max, net_r, net_path, net_name, model_type, r_solver,
						save_max_load, results, save_while, print_level)
	if isnothing(pd_max)
		max_load_file = joinpath(net_path, net_name * "_found_max_load.csv")
		if isfile(max_load_file )
			pd_max = readdlm(max_load_file , ',')
		else
			pm_pd_max = PM.instantiate_model(net_r, model_type, build_opf_var_load)
			pd_max, status = find_max_loads(pm_pd_max, print_level=print_level, solver=r_solver)

			save_while && writedlm(max_load_file, pd_max, ',')
		end
	end
	net_r["pd_max"] = pd_max
	if save_max_load
		!haskey(results, "load_constraints") && (results["load_constraints"]=Dict())
		results["load_constraints"]["pd_max"] = pd_max
	end
	return pd_max
end

"""
Sets the load values in the net to the given load profile
specified in a 2*n_load x 1 array of n_load real then n_load
reactive load powers.
"""
function set_network_load(net, new_load; scale_load=true)
    if scale_load
		new_load = new_load * 1 / net["baseMVA"]
	end

    num_loads = length(net["load"])
    for load in values(net["load"])
        load_index = load["index"]
        load["pd"] = float(new_load[load_index][1])
        load["qd"] = float(new_load[load_index + num_loads][1])
	end
end


text:
#returns solution for both parts
function x()
	f = open("6.txt")
	c1 = 0
	c2 = 0
	s1 = Set()
	s2 = Set()
	b = true
	for i=1:2149
		l = readline(f)
		if l == ""
			c1 += length(s1)
			c2 += length(s2)
			s1 = Set()
			s2 = Set()
			b = true
			continue
		end
		g = Set()
		for char in l
			push!(g, char)
		end
		union!(s1, g)
		if !b
			intersect!(s2, g)
		else
			s2 = g
			b = false
		end
	end
	(c1, c2)
end


text:
function part1(input)
    risk = reduce(hcat, (parse.(Int, collect(line))
                         for line in eachline(input)))
    return shortest_path(risk)
end

function part2(input)
    risk = reduce(hcat, (parse.(Int, collect(line))
                         for line in eachline(input)))
    h, w = size(risk)
    risk = repeat(risk, 5, 5)
    risk .+= div.(0:(size(risk, 1) - 1), h)
    risk .+= permutedims(div.(0:(size(risk, 2) - 1), w))
    risk .= mod1.(risk, 9)
    return shortest_path(risk)
end

function shortest_path(risk)
    queue = [(CartesianIndex((1, 1)), 0)]
    neighbors = CartesianIndex.(((1, 0), (-1, 0), (0, 1), (0, -1)))
    total_risk = fill(typemax(Int), size(risk))
    C = CartesianIndices(risk)
    finish = last(C)
    while true
        v, i = findmin(last, queue)
        pos = first(queue[i])
        pos == finish && return v
        deleteat!(queue, i)
        for neighbor in Ref(pos) .+ neighbors
            neighbor in C || continue
            neighbor_risk = v + risk[neighbor]
            neighbor_risk < total_risk[neighbor] || continue
            total_risk[neighbor] = neighbor_risk
            pushfirst!(queue, (neighbor, neighbor_risk))
        end
    end
end


text:
using TimeSeries,  MarketData, Base.Dates
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("collapse operations") do

    context("collapse squishes correctly") do

        @fact collapse(cl, week, first).values[2]     --> 97.75
        @fact collapse(cl, week, first).timestamp[2]  --> Date(2000,1,10)
        @fact collapse(cl, week, first, last).values[2]     --> 100.44
        @fact collapse(cl, week, first, last).timestamp[2]  --> Date(2000,1,10)
        @fact collapse(cl, month, first).values[2]    --> 100.25
        @fact collapse(cl, month, first).timestamp[2] --> Date(2000,2,1)

        @fact collapse(ohlc, week, first).values[2, :]  --> [102.0, 102.25, 94.75, 97.75]
        @fact collapse(ohlc, week, first).timestamp[2]  --> Date(2000,1,10)
        @fact collapse(ohlc, week, first, last).values[2, :]  --> [100.0, 102.25, 99.38, 100.44]
        @fact collapse(ohlc, week, first, last).timestamp[2]  --> Date(2000,1,10)
        @fact collapse(ohlc, month, first).values[2, :] --> [104.0, 105.0, 100.0, 100.25]
        @fact collapse(ohlc, month, first).timestamp[2] --> Date(2000,2,1)

    end
end

facts("merge works correctly") do

    cl1  = cl[1:3]
    op1  = cl[2:4]
    aapl = tail(AAPL)
    ba   = tail(BA)

    context("takes colnames kwarg correctly") do
        @fact merge(cl, ohlc["High", "Low"], colnames=["a","b","c"]).colnames --> ["a", "b", "c"]
        @fact merge(cl, op, colnames=["a","b"]).colnames                      --> ["a", "b"]
        @fact_throws merge(cl, op, colnames=["a"])
        @fact_throws merge(cl, op, colnames=["a","b","c"])

        @fact merge(cl, ohlc["High", "Low"], :inner, colnames=["a","b","c"]).colnames --> ["a", "b", "c"]
        @fact merge(cl, op, :inner, colnames=["a","b"]).colnames         --> ["a", "b"]
        @fact_throws merge(cl, op, :inner, colnames=["a"])
        @fact_throws merge(cl, op, :inner, colnames=["a","b","c"])

        @fact merge(cl, ohlc["High", "Low"], :left, colnames=["a","b","c"]).colnames --> ["a", "b", "c"]
        @fact merge(cl, op, :left, colnames=["a","b"]).colnames          --> ["a", "b"]
        @fact_throws merge(cl, op, :left, colnames=["a"])
        @fact_throws merge(cl, op, :left, colnames=["a","b","c"])

        @fact merge(cl, ohlc["High", "Low"], :right, colnames=["a","b","c"]).colnames --> ["a", "b", "c"]
        @fact merge(cl, op, :right, colnames=["a","b"]).colnames         --> ["a", "b"]
        @fact_throws merge(cl, op, :right, colnames=["a"])
        @fact_throws merge(cl, op, :right, colnames=["a","b","c"])

        @fact merge(cl, ohlc["High", "Low"], :outer, colnames=["a","b","c"]).colnames --> ["a", "b", "c"]
        @fact merge(cl, op, :outer, colnames=["a","b"]).colnames         --> ["a", "b"]
        @fact_throws merge(cl, op, :outer, colnames=["a"])
        @fact_throws merge(cl, op, :outer, colnames=["a","b","c"])
    end

    context("returns correct alignment with Dates and values") do
        @fact merge(cl,op).values --> merge(cl,op, :inner).values
        @fact merge(cl,op).values[2,1] --> cl.values[2,1]
        @fact merge(cl,op).values[2,2] --> op.values[2,1]

    end

    context("aligns with disparate sized objects") do
        @fact merge(cl, op[2:5]).values[1,1]  --> cl.values[2,1]
        @fact merge(cl, op[2:5]).values[1,2]  --> op.values[2,1]
        @fact merge(cl, op[2:5]).timestamp[1] --> Date(2000,1,4)
        @fact length(merge(cl, op[2:5]))      --> 4

        @fact length(merge(cl1, op1, :inner))    --> 2
        @fact merge(cl1,op1, :inner).values[2,1] --> cl1.values[3,1]
        @fact merge(cl1,op1, :inner).values[2,2] --> op1.values[2,1]

        @fact length(merge(cl1, op1, :left))     --> 3
        @fact merge(cl1,op1, :left).values[1,2]  --> isnan
        @fact merge(cl1,op1, :left).values[2,1]  --> cl1.values[2,1]
        @fact merge(cl1,op1, :left).values[2,2]  --> op1.values[1,1]

        @fact length(merge(cl1, op1, :right))    --> 3
        @fact merge(cl1,op1, :right).values[2,1] --> cl1.values[3,1]
        @fact merge(cl1,op1, :right).values[2,2] --> op1.values[2,1]
        @fact merge(cl1,op1, :right).values[3,1] --> isnan

        @fact length(merge(cl1, op1, :outer))    --> 4
        @fact merge(cl1,op1, :outer).values[1,2] --> isnan
        @fact merge(cl1,op1, :outer).values[2,1] --> cl1.values[2,1]
        @fact merge(cl1,op1, :outer).values[2,2] --> op1.values[1,1]
        @fact merge(cl1,op1, :outer).values[4,1] --> isnan
    end

    context("column names match the correct values") do
        @fact merge(cl, op[2:5]).colnames               --> ["Close", "Open"]
        @fact merge(op[2:5], cl).colnames               --> ["Open", "Close"]

        @fact merge(cl, op[2:5], :inner).colnames  --> ["Close", "Open"]
        @fact merge(op[2:5], cl, :inner).colnames  --> ["Open", "Close"]

        @fact merge(cl, op[2:5], :left).colnames   --> ["Close", "Open"]
        @fact merge(op[2:5], cl, :left).colnames   --> ["Open", "Close"]

        @fact merge(cl, op[2:5], :right).colnames  --> ["Close", "Open"]
        @fact merge(op[2:5], cl, :right).colnames  --> ["Open", "Close"]

        @fact merge(cl, op[2:5], :outer).colnames  --> ["Close", "Open"]
        @fact merge(op[2:5], cl, :outer).colnames  --> ["Open", "Close"]
    end
end

facts("vcat works correctly") do
    context("concatenates time series correctly in 1D") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], ["Number"])
        b = TimeArray([Date(2015, 12, 01)], [17], ["Number"])
        c = vcat(a, b)

        @fact length(c)  --> length(a) + length(b)
        @fact c.colnames --> a.colnames
        @fact c.colnames --> b.colnames
        @fact c.values   --> [15, 16, 17]
    end

    context("concatenates time series correctly in 2D") do
        a = TimeArray([Date(2015, 09, 01), Date(2015, 10, 01), Date(2015, 11, 01)], [[15 16]; [17 18]; [19 20]], ["Number 1", "Number 2"])
        b = TimeArray([Date(2015, 12, 01)], [18 18], ["Number 1", "Number 2"])
        c = vcat(a, b)

        @fact length(c)  --> length(a) + length(b)
        @fact c.colnames --> a.colnames
        @fact c.colnames --> b.colnames
        @fact c.values   --> [[15 16]; [17 18]; [19 20]; [18 18]]
    end

    context("rejects when column names do not match") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], ["Number"])
        b = TimeArray([Date(2015, 12, 01)], [17], ["Data does not match number"])

        @fact_throws vcat(a, b)
    end

    context("rejects when metas do not match") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], ["Number"], :FirstMeta)
        b = TimeArray([Date(2015, 12, 01)], [17], ["Number"], :SecondMeta)

        @fact_throws vcat(a, b)
    end

    context("rejects when dates overlap") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], ["Number"])
        b = TimeArray([Date(2015, 11, 01)], [17], ["Number"])

        @fact_throws vcat(a, b)
    end

    context("still works when dates are mixed") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 12, 01)], [15, 17], ["Number"])
        b = TimeArray([Date(2015, 11, 01)], [16], ["Number"])
        c = vcat(a, b)

        @fact length(c)    --> length(a) + length(b)
        @fact c.colnames   --> a.colnames
        @fact c.colnames   --> b.colnames
        @fact c.values     --> [15, 16, 17]
        @fact c.timestamp  --> issorted
    end
end

facts("map works correctly") do
    context("works on both time stamps and 1D values") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 11, 01)], [15, 16], ["Number"], :Something)
        b = map((timestamp, values) -> (timestamp + Dates.Year(1), values - 1), a)

        @fact length(b)                  --> length(a)
        @fact b.colnames                 --> a.colnames
        @fact Dates.year(b.timestamp[1]) --> Dates.year(a.timestamp[1]) + 1
        @fact b.values[1]                --> a.values[1] - 1
        @fact b.meta                     --> a.meta
    end

    context("works on both time stamps and 2D values") do
        a = TimeArray([Date(2015, 09, 01), Date(2015, 10, 01), Date(2015, 11, 01)], [[15 16]; [17 18]; [19 20]], ["Number 1", "Number 2"])
        b = map((timestamp, values) -> (timestamp + Dates.Year(1), [values[1] + 2, values[2] - 1]), a)

        @fact length(b)                  --> length(a)
        @fact b.colnames                 --> a.colnames
        @fact Dates.year(b.timestamp[1]) --> Dates.year(a.timestamp[1]) + 1
        @fact b.values[1, 1]             --> a.values[1, 1] + 2
        @fact b.values[1, 2]             --> a.values[1, 2] - 1
    end

    context("works with order of elements that varies after modifications") do
        a = TimeArray([Date(2015, 10, 01), Date(2015, 12, 01)], [15, 16], ["Number"])
        b = map((timestamp, values) -> (timestamp + Dates.Year((timestamp >= Date(2015, 11, 01)) ? -1 : 1), values), a)

        @fact length(b)    --> length(a)
        @fact b.timestamp  --> issorted
    end
end


text:
# Sample densities from the c-posterior, for the skew-normal example
module Densities
using Distributions
using PyPlot
draw() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())

# Settings
n = 10000  # sample size to use
alphas = [Inf,100]  # robustification params to use
n_samples = 10^4  # number of MCMC iterations

# Model parameters and sampler code
include("setup.jl")

for (i_a,alpha) in enumerate(alphas)
    x = [(rand()<p1 ? skewrnd(1,l1,s1,a1) : skewrnd(1,l2,s2,a2))[1]::Float64 for j=1:n] # Sample data

    # Run sampler
    zeta = (1/n)/(1/n + 1/alpha)
    p,theta,k_r,v_r,art,arv,m_r,s_r = sampler(x,n_samples,m,c,sigma,zeta)

    # Plot
    figure(i_a,figsize=(8,3.2)); clf(); hold(true)
    subplots_adjust(bottom=0.2)
    xmin,xmax = -5,5
    xs = linspace(xmin,xmax,1000)
    density = zeros(length(xs))
    for i = 1:m
        if p[i]>0
            pixs = p[i]*normpdf(xs,theta[i][1],theta[i][2])
            density += pixs
            plot(xs,pixs)
        end
    end
    plot(xs,density,"k--",linewidth=3)
    xlim(xmin,xmax)
    ylim(0,0.4)
    if alpha<Inf; title("\$\\mathrm{Coarsened\\,posterior}\$",fontsize=17)
    else; title("\$\\mathrm{Standard\\,posterior}\$",fontsize=17)
    end
    xlabel("\$x\$",fontsize=16)
    ylabel("\$\\mathrm{density}\$",fontsize=15)
    draw()
    savefig("density-alpha=$alpha.png",dpi=150)
end

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
module EMLcolumn
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...stage_separators.EMLtray
	import EMLtypes.length
	@reexport using ...stage_separators.EMLreboiler
	import EMLtypes.length
	@reexport using ...stage_separators.EMLcondenser
	import EMLtypes.length
	@reexport using ...mixers_splitters.EMLsplitter
	import EMLtypes.length
	@reexport using ...pressure_changers.EMLpump
	import EMLtypes.length
	include("column/ControlSection.jl")
	include("column/InitializeSection.jl")
	include("column/InitializeStage.jl")
	include("column/SectionColumnBasic.jl")
	include("column/Section_Column.jl")
	include("column/Section_Column2.jl")
	include("column/Section_Column3.jl")
	include("column/Packed_Section_ColumnBasic.jl")
	include("column/Packed_Section_Column.jl")
	include("column/Packed_Section_Column2.jl")
end