text:
plotstuff_pendcart(args...) = println("Install package Plots.jl (and call using Plots) to plot results in the end of demo_pendcart")

function care(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch
        error("R must be non-singular.")
    end
    Z = [A  -G;
    -Q  -A']

    S = schur(Z)
    S = ordschur(S, real(S.values).<0)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end
function lqr(A, B, Q, R)
    S = care(A, B, Q, R)
    K = R\B'*S
    return K
end

"""
    demo_pendcart(;kwargs...)

Run the iLQG Function to find an optimal trajectory For the "pendulum on a cart system". Requires package ControlSystems.jl

# Arguments
`x0     = [π-0.6,0,0,0]`
`goal   = [π,0,0,0]`
`Q      = Diagonal([10,1,2,1])` : State weight matrix
`R      = 1`                    : Control weight matrix
`lims   = 5.0*[-1 1]`           : control limits,
`T      = 600`                  : Number of time steps
`doplot = true`                 : Plot results
"""
function demo_pendcart(;x0 = [π-0.6,0,0,0], goal = [π,0,0,0],
    Q      = Diagonal([10.,1,2,1]), # State weight matrix
    R      = 1.,                    # Control weight matrix
    lims   = 5.0*[-1 1],            # control limits,
    T      = 600,                   # Number of time steps
    doplot = true                   # Plot results
    )

    N = T+1
    g = 9.82
    l = 0.35 # Length of pendulum
    h = 0.01 # Sample time
    d = 0.99
    A = [0 1 0 0; # Linearlized system dynamics matrix, continuous time
    g/l -d 0 0;
    0 0 0 1;
    0 0 0 0]
    B = [0, -1/l, 0, 1]
    C = eye(4) # Assume all states are measurable
    D = 4
    L = lqr(A,B,Q,R) # Calculate the optimal state feedback
    I = T

    function fsys_closedloop(t,x,L,xd)
        dx = copy(x)
        dx[1] -= pi
        u = -(L*dx)[1]
        xd[1] = x[2]
        xd[2] = -g/l * sin(x[1]) + u/l * cos(x[1]) - d*x[2]
        xd[3] = x[4]
        xd[4] = u
    end

    function fsys(t,x,u,xd)
        xd[1] = x[2]
        xd[2] = -g/l * sin(x[1]) + u/l * cos(x[1]) - d*x[2]
        xd[3] = x[4]
        xd[4] = u
    end

    dfvec = zeros(4)
    function dfsys(x,u)
        dfvec[1] = x[1]+h*x[2]
        dfvec[2] = x[2]+h*(-g/l*sin(x[1])+u[1]/l*cos(x[1])- d*x[2])
        dfvec[3] = x[3]+h*x[4]
        dfvec[4] = x[4]+h*u[1]
        dfvec
    end

    function cost_quadratic(x,u)
        local d = (x.-goal)
        0.5(d'*Q*d + u'R*u)[1]
    end

    function cost_quadratic(x::Matrix,u)
        local d = (x.-goal)
        T = size(u,2)
        c = Vector{Float64}(undef,T+1)
        for t = 1:T
            c[t] = 0.5(d[:,t]'*Q*d[:,t] + u[:,t]'R*u[:,t])[1]
        end
        c[end] = cost_quadratic(x[:,end][:],[0.0])
        return c
    end

    cx = zeros(4,T)
    cu = zeros(1,T)
    cxu = zeros(D,1)

    function dcost_quadratic(x,u)
        cx  .= Q*(x.-goal)
        cu  .= R.*u
        return cx,cu,cxu
    end

    function lin_dyn_f(x,u,i)
        u[isnan.(u)] .= 0
        f = dfsys(x,u)
    end

    fxc           = Array{Float64}(undef,D,D,I)
    fuc           = Array{Float64}(undef,D,1,I)
    fxd           = Array{Float64}(undef,D,D,I)
    fud           = Array{Float64}(undef,D,1,I)
    for ii = 1:I
        fxc[:,:,ii] = [0 1 0 0;
        0 0 0 0;
        0 0 0 1;
        0 0 0 0]
        fuc[:,:,ii] = [0, 0, 0, 1]
    end

    function lin_dyn_df(x,u)
        u[isnan.(u)] .= 0
        D            = size(x,1)
        nu,I         = size(u)
        cx,cu,cxu    = dcost_quadratic(x,u)
        cxx          = Q
        cuu          = [R]
        for ii = 1:I
            fxc[2,1,ii] = -g/l*cos(x[1,ii])-u[ii]/l*sin(x[1,ii])
            fxc[2,2,ii] = -d
            fuc[2,1,ii] = cos(x[1,ii])/l
            ABd = exp([fxc[:,:,ii]*h  fuc[:,:,ii]*h; zeros(nu, D + nu)])# ZoH sampling
            fxd[:,:,ii] = ABd[1:D,1:D]
            fud[:,:,ii] = ABd[1:D,D+1:D+nu]
        end
        fxx=fxu=fuu = []
        return fxd,fud,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu
    end

    x = zeros(4,N)
    u = zeros(1,T)

    """
    Simulate a pendulum on a cart using the non-linear equations
    """
    function simulate_pendcart(x0,L, dfsys, cost)
        x[:,1] = x0
        u[1] = 0
        for t = 2:T
            dx     = copy(x[:,t-1])
            dx[1] -= pi
            u[t]   = -(L*dx)[1]
            if !isempty(lims)
                u[t]   = clamp(u[t],lims[1],lims[2])
            end
            x[:,t] = dfsys(x[:,t-1],u[t])
        end
        dx      = copy(x[:,T])
        dx[1]  -= pi
        uT      = -(L*dx)[1]
        if !isempty(lims)
            uT   = clamp(uT,lims[1],lims[2])
        end
        x[:,T+1] = dfsys(x[:,T],uT)
        c = cost(x,u)

        return x, u, c
    end

    # Simulate the closed loop system with regular LQG control and watch it fail due to control limits
    x00, u00, cost00 = simulate_pendcart(x0, L, dfsys, cost_quadratic)

    f(x,u,i) = lin_dyn_f(x,u,i)
    df(x,u)  = lin_dyn_df(x,u)
    # plotFn(x)  = plot(squeeze(x,2)')

    println("Entering iLQG function")
    # subplot(n=4,nc=2)
    x, u, L, Vx, Vxx, cost, trace = iLQG(f,cost_quadratic, df, x0, 0*u00,
    lims      = lims,
    # plotFn  = x -> Plots.subplot!(x'),
    regType   = 2,
    α         = exp10.(range(0.2, stop=-3, length=6)),
    λmax      = 1e15,
    verbosity = 3,
    tol_fun   = 1e-8,
    tol_grad   = 1e-8,
    max_iter  = 1000);

    doplot && plotstuff_pendcart(x00, u00, x,u,cost00,cost,trace)
    println("Done")

    return x, u, L, Vx, Vxx, cost, trace
end


text:
module BetaTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, prod!, unsafeMean, unsafeLogMean, unsafeMirroredLogMean, unsafeVar, vague, dims, logPdf
using ForneyLab: SPBetaOutNMM, SPBetaMNM, SPBetaMMN, SPBetaOutMCNMM, VBBetaOut, VBBetaA, VBBetaB
using SpecialFunctions: digamma

@testset "Beta ProbabilityDistribution and Message construction" begin
    @test ProbabilityDistribution(Univariate, Beta, a=2.0, b=3.0) == ProbabilityDistribution{Univariate, Beta}(Dict(:a=>2.0, :b=>3.0))
    @test_throws Exception ProbabilityDistribution(Multivariate, Beta)
    @test ProbabilityDistribution(Beta, a=2.0, b=3.0) == ProbabilityDistribution{Univariate, Beta}(Dict(:a=>2.0, :b=>3.0))
    @test ProbabilityDistribution(Beta) == ProbabilityDistribution{Univariate, Beta}(Dict(:a=>1.0, :b=>1.0))
    @test Message(Beta) == Message{Beta, Univariate}(ProbabilityDistribution{Univariate, Beta}(Dict(:a=>1.0, :b=>1.0)))
    @test Message(Univariate, Beta) == Message{Beta, Univariate}(ProbabilityDistribution{Univariate, Beta}(Dict(:a=>1.0, :b=>1.0)))
    @test_throws Exception Message(Multivariate, Beta)
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Beta, a=2.0, b=2.0)) == 1
end

@testset "vague" begin
    @test vague(Beta) == ProbabilityDistribution(Beta, a=1.0, b=1.0)
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Beta, a=2.0, b=2.0)) == 0.5
    @test unsafeLogMean(ProbabilityDistribution(Beta, a=2.0, b=3.0)) == digamma(2.0) - digamma(5.0)
    @test unsafeMirroredLogMean(ProbabilityDistribution(Beta, a=2.0, b=3.0)) == digamma(3.0) - digamma(5.0)
    @test unsafeVar(ProbabilityDistribution(Beta, a=2.0, b=2.0)) == 0.05
end

@testset "log pdf" begin
    @test isapprox(logPdf(ProbabilityDistribution(Beta, a=2.0, b=2.0),0.3), 0.2311117209633866)
end

@testset "prod!" begin
    @test ProbabilityDistribution(Beta, a=2.0, b=2.0) * ProbabilityDistribution(Beta, a=2.0, b=3.0) == ProbabilityDistribution(Beta, a=3.0, b=4.0)
    @test ProbabilityDistribution(Univariate, Beta, a=1.0, b=2.0) * ProbabilityDistribution(Univariate, PointMass, m=0.2) == ProbabilityDistribution(Univariate, PointMass, m=0.2)
    @test ProbabilityDistribution(Univariate, PointMass, m=0.2) * ProbabilityDistribution(Univariate, Beta, a=1.0, b=2.0) == ProbabilityDistribution(Univariate, PointMass, m=0.2)
    @test_throws Exception ProbabilityDistribution(Univariate, PointMass, m=-1.0) * ProbabilityDistribution(Univariate, Beta, a=1.0, b=2.0)
end

#-------------
# Update rules
#-------------

@testset "SPBetaOutNMM" begin
    @test SPBetaOutNMM <: SumProductRule{Beta}
    @test outboundType(SPBetaOutNMM) == Message{Beta}
    @test isApplicable(SPBetaOutNMM, [Nothing, Message{PointMass}, Message{PointMass}])
    @test !isApplicable(SPBetaOutNMM, [Nothing, Message{FactorNode}, Message{PointMass}])
    @test !isApplicable(SPBetaOutNMM, [Nothing, Message{PointMass}, Message{FactorNode}])

    @test ruleSPBetaOutNMM(nothing, Message(Univariate, PointMass, m=2.0), Message(Univariate, PointMass, m=3.0)) == Message(Univariate, Beta, a=2.0, b=3.0)
end

@testset "SPBetaOutMCNMM" begin
    @test SPBetaOutMCNMM <: SumProductRule{Beta}
    @test outboundType(SPBetaOutMCNMM) == Message{SampleList}
    @test !isApplicable(SPBetaOutMCNMM, [Nothing, Message{PointMass}, Message{PointMass}])
    @test isApplicable(SPBetaOutMCNMM, [Nothing, Message{FactorNode}, Message{PointMass}])
    @test isApplicable(SPBetaOutMCNMM, [Nothing, Message{PointMass}, Message{FactorNode}])

    d = ruleSPBetaOutMCNMM(nothing, Message(Univariate, PointMass, m=2.0), Message(Univariate, Gamma, a=300.0, b=100.0)).dist
    @test 0.3<mean(d)<0.5
end

@testset "SPBetaMNM" begin
    @test SPBetaMNM <: SumProductRule{Beta}
    @test outboundType(SPBetaMNM) == Message{Function}
    @test isApplicable(SPBetaMNM, [Message,Nothing,Message])
end

@testset "SPBetaMMN" begin
    @test SPBetaMMN <: SumProductRule{Beta}
    @test outboundType(SPBetaMMN) == Message{Function}
    @test !isApplicable(SPBetaMMN, [Message,Nothing,Message])
    @test isApplicable(SPBetaMMN, [Message,Message,Nothing])
end

@testset "VBBetaOut" begin
    @test VBBetaOut <: NaiveVariationalRule{Beta}
    @test outboundType(VBBetaOut) == Message{Beta}
    @test isApplicable(VBBetaOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBBetaOut(nothing, ProbabilityDistribution(Univariate, PointMass, m=2.0), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(Univariate, Beta, a=2.0, b=3.0)
end

@testset "VBBetaA" begin
    @test VBBetaA <: NaiveVariationalRule{Beta}
    @test outboundType(VBBetaA) == Message{Function}
    @test isApplicable(VBBetaA, [ProbabilityDistribution, Nothing, ProbabilityDistribution])
end

@testset "VBBetaB" begin
    @test VBBetaB <: NaiveVariationalRule{Beta}
    @test outboundType(VBBetaB) == Message{Function}
    @test isApplicable(VBBetaB, [ProbabilityDistribution, ProbabilityDistribution, Nothing])
end

@testset "averageEnergy and differentialEntropy" begin
    @test isapprox(differentialEntropy(ProbabilityDistribution(Univariate, Beta, a=2.0, b=3.0)), averageEnergy(Beta, ProbabilityDistribution(Univariate, Beta, a=2.0, b=3.0), ProbabilityDistribution(Univariate, PointMass, m=2.0), ProbabilityDistribution(Univariate, PointMass, m=3.0)))
end

end # module

ext: jl


text:
"""
A `FujiRequest` represents an HTTP request sent to the Fuji webserver. It
currently has five fields:

* `headers`: A dictionary of the request headers (e.g. Host, User-Agent)
* `method`: The HTTP method used (e.g. GET, POST)
* `params`: A dictionary of the named parameters used (e.g. if the endpoint was
    /hello/:name and the request was /hello/jack, this would be
    ["name": "jack"])
* `query_params`: A dictionary of the query parameters used (e.g. if the
    endpoint was /hello and the request was /hello?name=jack, this would be
    ["name": "jack"])
* `resource`: The resource of the request. Will typically be the part after the
    domain name but before the query parameters (e.g. /hello/jack)
* `splat`: An array of any splat parameters used (e.g. if the endpoint was
    /hello/* and the request was /hello/jack, this would be ["jack"])
"""
type FujiRequest
    headers::Dict{String,String}
    method::String
    params::Dict{String,String}
    query_params::Dict{String,String}
    resource::String
    splat::Array{String,1}

    FujiRequest(req::Request) = new(
        Dict{String,String}(req.headers), # headers
        req.method, # method
        Dict{String,String}(), # params
        Dict{String,String}(), # query params
        req.resource, # resource
        Array{String,1}() # splat
    )

    FujiRequest(req::Request, route::Route) = begin
        request = FujiRequest(req)

        named_param_regex = r":([A-z0-9\-_]+)"
        splat_param_regex = r"\*"
        component_string = "([A-z0-9\\-_]+)"

        # if there are named parameters in the endpoint
        if ismatch(named_param_regex, route.endpoint)
            params = Dict{String,String}()

            named_names = Array{String,1}()
            named_values = Array{String,1}()

            endpoint = route.endpoint

            while ismatch(named_param_regex, endpoint)
                m = match(named_param_regex, endpoint)
                push!(named_names, m.captures[1])

                endpoint = replace(route.endpoint, named_param_regex, m.captures[1], 1)
            end

            regex_endpoint = Regex(string("^", replace(route.endpoint, named_param_regex, component_string), "\$"))
            resource_match = match(regex_endpoint, split(request.resource, "?")[1])

            for capture in resource_match.captures
                push!(named_values, capture)
            end

            for i in 1:length(named_names)
                params[named_names[i]] = named_values[i]
            end

            request.params = params
        end

        # if there are splat parameters in the endpoint
        if contains(route.endpoint, "*")
            splat = Array{String,1}()

            regex_endpoint = Regex(string("^", replace(route.endpoint, splat_param_regex, component_string), "\$"))
            resource_match = match(regex_endpoint, split(request.resource, "?")[1])

            for capture in resource_match.captures
                push!(splat, capture)
            end

            request.splat = splat
        end

        # if there are query parameters in the endpoint
        if contains(request.resource, "?")
            query_params_string = split(request.resource, "?")[2]
            query_params_array = split(query_params_string, "&")
            query_params = Dict{String,String}()

            for query_param in query_params_array
                if contains(query_param, "=")
                    data = split(query_param, "=")
                    query_params[data[1]] = data[2]
                end
            end

            request.query_params = query_params
        end

        request
    end
end


text:
using Rimu, Rimu.StatsTools
using Statistics, DataFrames, Random, Distributions
using Test
using KrylovKit, LinearAlgebra
using MonteCarloMeasurements

import Measurements

@testset "smoothen" begin
    noisy = randn(1000)
    b = 20
    smooth = smoothen(noisy, b)
    @test length(smooth) == length(noisy)
    @test mean(smooth) ≈ mean(noisy)
    @test var(noisy) > var(smooth)
    @test length(smoothen(noisy, b)) == length(noisy)
end

@testset "growth_witness" begin
    @test_throws AssertionError growth_witness(rand(10),rand(11),0.01)
    @test_throws AssertionError growth_witness(rand(5),rand(5),0.01; skip=10)
    nor = rand(200)
    shift = rand(200)
    dτ = 0.01*ones(200)
    df = DataFrame(; norm=nor, shift, dτ)
    m = mean(growth_witness(shift, nor, dτ[1]; skip=10))
    @test mean(growth_witness(df, 10; skip=10)) ≈ m
    @test mean(growth_witness(df; skip=10)) ≈ m
end

using Rimu.StatsTools: blocker
@testset "blocking" begin
    Random.seed!(17) # make sure we don't trip over rare fluctuations
    # real
    v = randn(2^10)
    br = blocking_analysis(v)
    @test MonteCarloMeasurements.:±(br) ≈ Measurements.:±(br)
    @test 0.01 < br.err < 0.04
    @test br.k ≤ 4 # should be 0+1 (independent random variables)
    brs = blocking_analysis(smoothen(v, 2^5)) # introduce correlation
    @test mean(v) ≈ br.mean ≈ brs.mean
    @test 0.01 < brs.err < 0.04
    @test 5 < brs.k ≤ 7 # should be 5+1
    bor = blocking_analysis(ones(2000)) # blocking fails
    @test bor.k == -1
    @test isnan(bor.err)

    # complex
    w = randn(ComplexF64, 2^10)
    bc = blocking_analysis(w)
    @test real(MonteCarloMeasurements.:±(bc)) ≈ real(Measurements.:±(bc))
    @test imag(MonteCarloMeasurements.:±(bc)) ≈ imag(Measurements.:±(bc))
    @test bc.k ≤ 4  # should be 0+1 (independent random variables)
    bcs = blocking_analysis(smoothen(w, 2^5))
    @test mean(w) ≈ bc.mean ≈ bcs.mean
    @test 5 < bcs.k ≤ 7 # should be 5+1

    @test blocking_analysis([1]).k == -1 == blocking_analysis(Int[]).k
end

using Rimu.StatsTools: x_by_y_linear, ratio_estimators, particles
@testset "ratio_of_means" begin
    Random.seed!(17) # make sure the tests don't trip over rare fluctuations
    n_samples = 2000
    br= ratio_of_means(rand(n_samples),rand(n_samples); α=0.05)
    @test br.blocks * 2^(br.k-1) ≤ n_samples
    @test begin show(br); true; end # does not throw error

    # Make sure all other exported functions work.
    for f in (
        pmedian, pmiddle, piterate, pextrema, pminimum, pmaximum, pmean, pcov
    )
        @test f(br) == f(br.ratio)
    end

    # badly behaved example of Cauchy distribution for ratio
    r = ratio_of_means(randn(2000),randn(2000))
    q95 = pquantile(r.ratio, [0.025,0.975])
    @test q95[1] < 0 < q95[2]

    # complex time series
    r = ratio_of_means(randn(ComplexF64,2000),randn(ComplexF64,2000))
    qr95 = pquantile(real(r.ratio), [0.025,0.975])
    @test qr95[1] < 0 < qr95[2]
    qi95 = pquantile(imag(r.ratio), [0.025,0.975])
    @test qi95[1] < 0 < qi95[2]
    @test_throws ErrorException pquantile(r, [0.025,0.975])
    @test begin show(r); true; end # does not throw error

    # mixed real and complex
    r = @inferred ratio_of_means(randn(1000),randn(ComplexF64,1000))
    @test r.k ≤ 4
    r = @inferred ratio_of_means(randn(ComplexF64,1000),randn(1000))
    @test r.k ≤ 4

    # zero variance example
    r = ratio_of_means([0,0,0,0,0,0], [1,1,1,1,1,2])
    @test Tuple(val_and_errs(r)) == (0,0,0)
    r = ratio_of_means(complex.(ones(2000)),complex.(ones(2000)))
    @test pvar(real(r.ratio)) == 0 == pvar(imag(r.ratio))

    # well behaved real example
    Random.seed!(17) # make sure the tests don't trip over rare fluctuations
    n_samples = 2000
    μ_a, μ_b = 2.0, 3.0
    σ_a, σ_b = 0.5, 0.1 # std of sample means
    f, σ_f = x_by_y_linear(μ_a, μ_b, σ_a, σ_b, 0) # expected ratio and std, uncorrelated
    @test f ≈ μ_a/μ_b
    # consistency with error propagation from Measurement.jl
    am = Measurements.:±(μ_a, σ_a)
    bm = Measurements.:±(μ_b, σ_b)
    @test σ_f ≈ (am / bm).err
    # consistency with error propagation from MonteCarloMeasurements.jl
    amc = MonteCarloMeasurements.:±(μ_a, σ_a)
    bmc = MonteCarloMeasurements.:±(μ_b, σ_b)
    @test isapprox(
        [f-σ_f, f, f+σ_f], pquantile(amc / bmc, [0.16,0.5,0.84]),
        rtol = 2/√2000
    )

    # add correlation for testing `ratio_of_means`
    Random.seed!(17) # make sure the tests don't trip over rare fluctuations
    ρ = 0.02
    f, σ_f = x_by_y_linear(μ_a, μ_b, σ_a, σ_b, ρ) # expected ratio and std
    @test f ≈ μ_a/μ_b
    ab = rand(MvNormal([μ_a, μ_b], [σ_a^2 ρ; ρ σ_b^2]*n_samples), n_samples)
    a = ab[1,:]
    b = ab[2,:]
    r = ratio_of_means(a,b) #; mc_samples = Val(10_000))
    @test 1 < r.k ≤ 4 # weakly correlated samples
    @test isapprox(μ_a/μ_b, r.f; atol = 2σ_f)
    @test isapprox(μ_a/μ_b, pmedian(r.ratio); atol = 2σ_f)
    q = pquantile(r.ratio, [0.16,0.5,0.84])
    @test isapprox(q[3]-q[2], σ_f; rtol = 4/√n_samples)
    @test isapprox(q[2]-q[1], σ_f; rtol = 4/√n_samples)
    @test isapprox(r.σ_f, σ_f; rtol = 4/√n_samples)

    # correlated time series
    rs = ratio_of_means(smoothen(a, 2^3), smoothen(b, 2^6))
    @test 6 ≤ rs.k ≤ 8 # ideally k==7 for 6 blocking steps to decorrelate
    # quantiles after blocking stay the same
    @test isapprox(q, pquantile(rs.ratio, [0.16,0.5,0.84]); atol = 2σ_f)
    @test isapprox(μ_a/μ_b, rs.f; atol = 2σ_f)
    @test isapprox(rs.σ_f, σ_f; atol = 2σ_f)

    Random.seed!(1234) # make sure the tests don't trip over rare fluctuations
    # well behaved complex example
    n_samples = 2000
    μ_a, μ_b = 2.0 + 1.0im, 3.0 - 2.0im
    σ_a, σ_b = 0.5, 0.1 # std of sample means
    f, σ_f = x_by_y_linear(μ_a, μ_b, σ_a, σ_b, 0) # expected ratio and std
    @test f ≈ μ_a/μ_b

    a = μ_a .+ √n_samples * σ_a * randn(ComplexF64, n_samples)
    b = μ_b .+ √n_samples * σ_b * randn(ComplexF64, n_samples)
    r = ratio_of_means(a,b)

    @test r.k < 3 # uncorrelated samples
    @test isapprox(μ_a/μ_b, r.f; atol = abs(2σ_f))
    @test isapprox(μ_a/μ_b, pmedian(r); atol = abs(2σ_f))
    qr = pquantile(real(r.ratio), [0.16,0.5,0.84])
    qi = pquantile(imag(r.ratio), [0.16,0.5,0.84])
    # Is this the correct way to test against linear error propagation?
    @test qr[3]-qr[2] < 2abs(σ_f)
    @test qi[3]-qi[2] < 2abs(σ_f)
    @test qr[2]-qr[1] < 2abs(σ_f)
    @test qi[2]-qi[1] < 2abs(σ_f)
    @test isapprox(r.σ_f, σ_f; rtol = 4/√n_samples)

    # type stability of Particles
    d = MvNormal([1.0,1.0],[0.1 0.01;0.01 0.1])
    @test typeof(particles(100, d)) == typeof(particles(Val(100), d))
    @test typeof(particles(100, d)) == typeof(particles(100, [1.0,1.0],[0.1 0.01;0.01 0.1]))
    @inferred particles(nothing, d)
    @inferred particles(Val(100), d)
    @inferred ratio_estimators(rand(100),rand(100),2; mc_samples = Val(100))
    @inferred ratio_estimators(
        rand(ComplexF64,100),rand(ComplexF64,100), 2;
        mc_samples = Val(100),
    )
    @inferred ratio_of_means(rand(1000), 100 .+ rand(1000))
    p = @inferred particles(nothing,[1.0,1.0],[0.1 0.01;0.01 0.0])
    @test pvar.(p)[2] == 0
end

@testset "Reweighting" begin
    ham = HubbardReal1D(BoseFS((1,1,1,1)), u=6.0, t=1.0)
    # using KrylovKit
    # fv = DVec(starting_address(ham)=>1.0; capacity=dimension(ham))
    # kkresults = eigsolve(ham, fv, 1, :SR; issymmetric = true)
    # exact_energy = kkresults[1][1]
    exact_energy = -2.869739978337469
    # run integer walker FCIQMC to get significant bias
    v = DVec(starting_address(ham)=>2; capacity=dimension(ham))
    steps_equi = 200
    steps_meas = 2^10
    p = RunTillLastStep(laststep = steps_equi+steps_meas)
    post_step = ProjectedEnergy(ham, v)
    s_strat = DoubleLogUpdate(targetwalkers=10)
    seedCRNG!(174)
    @time df = lomc!(ham, v; params=p, s_strat, post_step).df
    bs = shift_estimator(df; skip=steps_equi)
    @test bs == blocking_analysis(df.shift[steps_equi+1:end])
    pcb = bs.mean - exact_energy
    # TODO
    @test pcb > 0.0 # the shift has a large population control bias
    # test growth_estimator
    h = 2^(bs.k-1) # approximate number of steps to decorrelate the shift
    E_r = bs.mean # set up reference energy
    ge = growth_estimator(df, h; E_r, skip=steps_equi)
    pcb_est = E_r - ge.ratio # estimated PCB in the shift from reweighting
    @test 0.2 < pmedian(pcb_est) < pcb
    @inferred growth_estimator(rand(1000), 100 .+ rand(1000), 20, 0.01; change_type = to_measurement)
    @inferred growth_estimator(rand(1000), 100 .+ rand(1000), 20, 0.01)
    # fails due to type instability in MonteCarloMeasurements
    # test w_lin()
    @test ge.ratio ≈ growth_estimator(df, h; E_r, skip=steps_equi, weights = w_lin).ratio
    # projected energy
    bp = ratio_of_means(df.hproj[steps_equi+1:end],df.vproj[steps_equi+1:end])
    bp_intervals = ratio_with_errs(bp)
    @test pmedian(bp) ≈ bp_intervals.ratio
    me = @inferred mixed_estimator(
        df.hproj, df.vproj, df.shift, h, df.dτ[end];
        skip=steps_equi, E_r
    )
    @test me.ratio≈bp.ratio # reweighting has not significantly improved the projected energy
    val(projected_energy(df; skip=steps_equi)) ≈ val(mixed_estimator(df,0; skip=steps_equi))
end

using Rimu.StatsTools: replica_fidelity
@testset "Fidelity" begin
    ham = HubbardReal1D(BoseFS((1,1,1,1)), u=6.0, t=1.0)

    # get exact eigenvectors with KrylovKit
    fv = DVec(starting_address(ham)=>1.0; capacity=dimension(ham))
    kkresults = eigsolve(ham, fv, 2, :SR; issymmetric = true)
    gs = kkresults[2][1] # ground state; norm(gs) ≈ 1
    es = kkresults[2][2] # 1st excited state; norm(es) ≈ 1
    @test norm(gs) ≈ 1 ≈ norm(es) # they are normalised
    @test isapprox(gs⋅es, 0; atol = √eps(Float64)) # and orthogonal

    # Set up oblique vector at angle to ground state
    α = π/3 # 60° angle
    os = add!(cos(α)*gs, sin(α)*es)
    @test norm(os) ≈ 1

    # set up replica MC
    v = DVec(starting_address(ham)=>2; capacity=dimension(ham))
    steps_equi = 200
    steps_meas = 2^10
    p = RunTillLastStep(laststep = steps_equi+steps_meas)
    post_step = (Projector(vproj=gs), Projector(hproj=os))
    s_strat = DoubleLogUpdate(targetwalkers=10)

    # run replica fciqmc
    seedCRNG!(17)
    @time rr = lomc!(ham, v; params=p, s_strat, post_step, replica=AllOverlaps()).df

    # check fidelity with ground state
    fid_gs = replica_fidelity(rr; p_field=:vproj, skip=steps_equi)
    @test fid_gs.ratio ≈ 1
    re = ratio_with_errs(fid_gs) # extract errors from quantiles
    @test re.err1_l < 0.03 && re.err1_u < 0.03 # errors are small
    # TODO
    #=
    # check fidelity with oblique state
    fid_os = StatsTools.replica_fidelity(rr; p_field=:hproj, skip=steps_equi)
    @test fid_os.ratio ≈ cos(α)^2
    γ = acos(√fid_os.ratio) # quantum angle or Fubini-Study metric
    @test γ ≈ α
    =#
end

comp_tuples(a,b; atol=0) = mapreduce((x,y)->isapprox(x,y; atol), &, Tuple(a), Tuple(b))
@testset "convenience" begin
    v, σ = 2.0, 0.2
    m = Measurements.measurement(v, σ)
    @test comp_tuples(med_and_errs(m), (v, σ, σ, 2σ, 2σ))
    # @test map((x,y)->isapprox(x,y),Tuple(med_and_errs(m)),(v, σ, σ, 2σ, 2σ))|>all
    # @test med_and_errs(m) == (med = v, err1_l = σ, err1_u = σ, err2_l = σ, err2_u = σ)
    mp = MonteCarloMeasurements.Particles(m)
    @test comp_tuples(med_and_errs(mp), (v, σ, σ, 2σ, 2σ); atol=0.01)
    # @test map((x,y)->isapprox(x,y; atol=0.01),Tuple(med_and_errs(mp)),(v, σ, σ, 2σ, 2σ))|>all
    @test m ≈ to_measurement(mp)

    @test val(2.0) == 2.0
    @test val(m) == v
    @test val(mp) == pmedian(mp)
    @test errs(2.0) == (err_l=0, err_u=0)
    @test val_and_errs(2; name="x") == (x=2, x_l=0, x_u=0)
    @test comp_tuples(val_and_errs(m, p=0.954499736104), (v, 2σ, 2σ))
    @test comp_tuples(val_and_errs(m, n=2), (v, 2σ, 2σ))
    @test comp_tuples(val_and_errs(mp, n=2), (v, 2σ, 2σ); atol=0.01)

    r = ratio_of_means(randn(2000),randn(2000))
    @test comp_tuples(val_and_errs(r), (v = val(r), errs(r)...))
    @test NamedTuple(r).val_k == r.k
    @test NamedTuple(r; name=:frodo).frodo == pmedian(r)
    br = blocking_analysis(rand(2000))
    @test comp_tuples(val_and_errs(br), (v = val(br), errs(br)...))
    @test NamedTuple(br).val_blocks == br.blocks

    # complex time series
    r = ratio_of_means(randn(ComplexF64, 2000), randn(ComplexF64, 2000))
    nt = val_and_errs(r)
    @test imag(nt.val_l) ≠ 0 ≠ real(nt.val_l)
    r = blocking_analysis(rand(ComplexF64,2000))
    nt = val_and_errs(r)
    @test imag(nt.val_l) ≠ 0 ≠ real(nt.val_l)
end


text:
using JuliaWebAPI
using Logging
using Base.Test
using Compat

const opts = Base.JLOptions()
const inline_flag = opts.can_inline == 1 ? `` : `--inline=no`
const cov_flag = (opts.code_coverage == 1) ? `--code-coverage=user` :
                 (opts.code_coverage == 2) ? `--code-coverage=all` :
                 ``

function run_test(script, flags)
    srvrscript = joinpath(dirname(@__FILE__), script)
    srvrcmd = `$(joinpath(JULIA_HOME, "julia")) $cov_flag $inline_flag $script $flags`
    println("Running tests from ", script, "\n", "="^60)
    ret = run(srvrcmd)
    println("Finished ", script, "\n", "="^60)
    nothing
end

run_test("test_amqp_transport.jl", "--ser")
run_test("test_amqp_transport.jl", "--json")
run_test("test_apiclntsrvr.jl", "--run")

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins


text:
let
    t1 = 0.0
    t2 = 1.0
    t = 0.5
    @test_throws AssertionError verify(TimeInterval(t2,t1))
    interval = TimeInterval(t1,t2)
    @test get_start_time(interval) == t1
    @test get_end_time(interval) == t2
    @test get_Δt(interval,t) == (t - t1) / (t2 - t1)
end
let
    TrajectoryPoint()
    @test_throws AssertionError verify(TrajectoryPoint(t=NaN))

    x = [0.0,1.0]
    y = [1.0,0.0]
    θ = [0.0,π]
    vx = [0.0,1.0]
    vy = [1.0,0.0]
    ω = [1.0,2.0]
    t = [0.0,1.0]
    pt1 = TrajectoryPoint(
        pos=VecE2(x[1],y[1]),
        heading=VecE2(cos(θ[1]),sin(θ[1])),
        vel = VecE2(vx[1],vy[1]),
        yaw_rate = ω[1],
        t = t[1]
        )
    pt2 = TrajectoryPoint(
        pos=VecE2(x[2],y[2]),
        heading=VecE2(cos(θ[2]),sin(θ[2])),
        vel = VecE2(vx[2],vy[2]),
        yaw_rate = ω[2],
        t = t[2]
        )
    t_idx = 0.5
    pt = interpolate(pt1,pt2,t_idx)
    @test array_isapprox(pt.pos, (1-t_idx)*pt1.pos + t_idx*pt2.pos)
    @test isapprox(atan(pt.heading), (1-t_idx)*θ[1] + t_idx*θ[2])
    @test array_isapprox(pt.vel, (1-t_idx)*pt1.vel + t_idx*pt2.vel)
    @test isapprox(pt.yaw_rate, (1-t_idx)*pt1.yaw_rate + t_idx*pt2.yaw_rate)
    @test isapprox(pt.t, (1-t_idx)*pt1.t + t_idx*pt2.t)
end
# StraightTrajectory
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    traj = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    position = p1+dt*(p2-p1)
    velocity = (p2-p1)/(t2-t1)
    yaw_rate = 0.0

    @test get_start_time(traj) == get_start_time(interval)
    @test get_end_time(traj) == get_end_time(interval)
    @test get_Δt(traj,t) == dt

    @test get_time_from_arc_length(traj,arclength*(t-t1)/(t2-t1)) == t

    @test get_length(traj) == arclength
    @test get_dist(traj, (t1+t2)/2) == arclength/2

    @test array_isapprox(get_start_pt(traj),p1)
    @test array_isapprox(get_end_pt(traj),p2)
    @test array_isapprox(get_position(traj,t), position)
    @test array_isapprox(get_heading(traj,t), heading)
    @test array_isapprox(get_vel(traj,t), velocity)
    @test get_yaw_rate(traj,t) == yaw_rate
    @test get_time_from_pt(traj,get_position(traj,t)) == t

    traj_pt1 = get_trajectory_point_by_time(traj,t)
    @test array_isapprox(traj_pt1.pos, position)
    @test array_isapprox(traj_pt1.heading, heading)
    @test array_isapprox(traj_pt1.vel, velocity)
    @test traj_pt1.yaw_rate == yaw_rate

    traj_pt2 = get_trajectory_point_by_pt(traj, get_position(traj,t))
    @test array_isapprox(traj_pt2.pos, position)
    @test array_isapprox(traj_pt2.heading, heading)
    @test array_isapprox(traj_pt2.vel, velocity)
    @test traj_pt2.yaw_rate == yaw_rate
end
# ArcTrajectory
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)

    center = VecE2(0.0,0.0)
    radius = 1.0
    θ1 = 0.0
    Δθ = π/2
    θ2 = θ1 + Δθ
    heading1 = VecE2(-sin(θ1),cos(θ1))
    if sign(Δθ) != 0
        heading1 = VecE2(-sin(θ1),cos(θ1)) * sign(Δθ)
    end
    θ = θ1 + dt*(θ2-θ1)
    p1 = center + radius*VecE2(cos(θ1),sin(θ1))
    p2 = center + radius*VecE2(cos(θ2),sin(θ2))
    arclength = radius * abs(θ2-θ1)

    traj = ArcTrajectory(center,p1,heading1,Δθ,interval)
    # traj1 = ArcTrajectory(center,radius,θ1,Δθ,interval)
    # @test traj1.θ1 == traj.θ1

    position = center+radius*VecE2(cos(θ),sin(θ))
    heading = VecE2(-sin(θ),cos(θ))
    if sign(Δθ) != 0
        heading = VecE2(-sin(θ),cos(θ)) * sign(Δθ)
    end
    velocity = heading * arclength / (t2-t1)
    yaw_rate = (θ2-θ1) / (t2-t1)

    @test get_start_time(traj) == get_start_time(interval)
    @test get_end_time(traj) == get_end_time(interval)
    @test get_Δt(traj,t) == dt

    @test get_length(traj) == arclength
    @test get_dist(traj, (t1+t2)/2) == arclength/2

    @test array_isapprox(get_start_pt(traj),p1)
    @test array_isapprox(get_end_pt(traj),p2)
    @test array_isapprox(get_position(traj,t), position)
    @test array_isapprox(get_heading(traj,t), heading)
    @test array_isapprox(get_vel(traj,t), velocity)
    @test get_yaw_rate(traj,t) == yaw_rate
    # @test get_time_from_pt(traj,get_position(traj,t)) == t

    traj_pt1 = get_trajectory_point_by_time(traj,t)
    @test array_isapprox(traj_pt1.pos, position)
    @test array_isapprox(traj_pt1.heading, heading)
    @test array_isapprox(traj_pt1.vel, velocity)
    @test traj_pt1.yaw_rate == yaw_rate

    traj_pt2 = get_trajectory_point_by_pt(traj, get_position(traj,t))
    @test array_isapprox(traj_pt2.pos, position)
    @test array_isapprox(traj_pt2.heading, heading)
    @test array_isapprox(traj_pt2.vel, velocity)
    @test traj_pt2.yaw_rate == yaw_rate
end
# WaitTrajectory
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)
    p1 = VecE2(0.0,0.0)
    p2 = p1
    heading = VecE2(0.0,1.0)
    arclength = norm(p2-p1)
    traj = WaitTrajectory(p1,heading,TimeInterval(t1,t2))

    position = p1+dt*(p2-p1)
    velocity = (p2-p1)/(t2-t1)
    yaw_rate = 0.0

    @test get_start_time(traj) == get_start_time(interval)
    @test get_end_time(traj) == get_end_time(interval)
    @test get_Δt(traj,t) == dt

    @test get_length(traj) == arclength
    @test get_dist(traj, (t1+t2)/2) == arclength/2

    @test array_isapprox(get_start_pt(traj),p1)
    @test array_isapprox(get_end_pt(traj),p2)
    @test array_isapprox(get_position(traj,t), position)
    @test array_isapprox(get_heading(traj,t), heading)
    @test array_isapprox(get_vel(traj,t), velocity)
    @test get_yaw_rate(traj,t) == yaw_rate
    @test get_time_from_pt(traj,get_position(traj,t)) == get_start_time(traj)

    traj_pt1 = get_trajectory_point_by_time(traj,t)
    @test array_isapprox(traj_pt1.pos, position)
    @test array_isapprox(traj_pt1.heading, heading)
    @test array_isapprox(traj_pt1.vel, velocity)
    @test traj_pt1.yaw_rate == yaw_rate

    traj_pt2 = get_trajectory_point_by_pt(traj, get_position(traj,t))
    @test array_isapprox(traj_pt2.pos, position)
    @test array_isapprox(traj_pt2.heading, heading)
    @test array_isapprox(traj_pt2.vel, velocity)
    @test traj_pt2.yaw_rate == yaw_rate
end
# PivotTrajectory
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)

    center = VecE2(0.0,0.0)
    radius = 0.0
    θ1 = 0.0
    Δθ = 0.0
    θ2 = θ1 + Δθ
    heading1 = VecE2(-sin(θ1),cos(θ1))
    if sign(Δθ) != 0
        heading1 = VecE2(-sin(θ1),cos(θ1)) * sign(Δθ)
    end
    θ = θ1 + dt*(θ2-θ1)
    p1 = center + radius*VecE2(cos(θ1),sin(θ1))
    p2 = center + radius*VecE2(cos(θ2),sin(θ2))
    arclength = radius * abs(θ2-θ1)

    traj = PivotTrajectory(p1,heading1,Δθ,interval)

    position = center+radius*VecE2(cos(θ),sin(θ))
    heading = VecE2(-sin(θ),cos(θ))
    if sign(Δθ) != 0
        heading = VecE2(-sin(θ),cos(θ)) * sign(Δθ)
    end
    velocity = heading * arclength / (t2-t1)
    yaw_rate = (θ2-θ1) / (t2-t1)

    @test get_start_time(traj) == get_start_time(interval)
    @test get_end_time(traj) == get_end_time(interval)
    @test get_Δt(traj,t) == dt

    @test get_length(traj) == arclength
    @test get_dist(traj, (t1+t2)/2) == arclength/2

    @test array_isapprox(get_start_pt(traj),p1)
    @test array_isapprox(get_end_pt(traj),p2)
    @test array_isapprox(get_position(traj,t), position)
    @test array_isapprox(get_heading(traj,t), heading)
    @test array_isapprox(get_vel(traj,t), velocity)
    @test get_yaw_rate(traj,t) == yaw_rate
    # @test get_time_from_pt(traj,get_position(traj,t)) == t

    traj_pt1 = get_trajectory_point_by_time(traj,t)
    @test array_isapprox(traj_pt1.pos, position)
    @test array_isapprox(traj_pt1.heading, heading)
    @test array_isapprox(traj_pt1.vel, velocity)
    @test traj_pt1.yaw_rate == yaw_rate

    traj_pt2 = get_trajectory_point_by_pt(traj, get_position(traj,t))
    @test array_isapprox(traj_pt2.pos, position)
    @test array_isapprox(traj_pt2.heading, heading)
    @test array_isapprox(traj_pt2.vel, velocity)
    @test traj_pt2.yaw_rate == yaw_rate
end
# Reverse Trajectory
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    traj = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    reversed_traj = ReverseTrajectory(traj)
    # @test get_heading(reversed_traj,t)  == -1.0 * get_heading(traj,t)
    # @test get_vel(reversed_traj,t)      == -1.0 * get_vel(traj,t)
    @test get_position(reversed_traj,t) == get_position(reversed_traj,t)
    @test get_yaw_rate(reversed_traj,t) == get_yaw_rate(reversed_traj,t)

    pt = get_trajectory_point_by_time(traj,t)
    reversed_pt = get_trajectory_point_by_time(reversed_traj,t)
    # @test reversed_pt.heading   == -1.0 * pt.heading
    # @test reversed_pt.vel       == -1.0 * pt.vel
    @test reversed_pt.pos       == pt.pos
    @test reversed_pt.yaw_rate  == pt.yaw_rate

    dense_traj = DenseTrajectory(
        reversed_traj,
        [t1,t2],
        [0.0],
        (arclength / (t2-t1))*[1.0, 1.0],
        [0.0, arclength]
    )
    @test get_heading(reversed_traj,t)  == get_heading(dense_traj,t)
    @test get_vel(reversed_traj,t)      == get_vel(dense_traj,t)
    @test get_position(reversed_traj,t) == get_position(dense_traj,t)
    @test get_yaw_rate(reversed_traj,t) == get_yaw_rate(dense_traj,t)

end
# Trajectory
let
    t1 = 0.0
    t2 = 1.0
    t = 0.5
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    seg1 = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    t3 = 2.0
    center = VecE2(-1.0,1.0)
    radius = 1.0
    θ1 = 0.0
    Δθ = π/2
    θ2 = θ1 + Δθ
    seg2 = ArcTrajectory(center,p2,get_heading(seg1,get_end_time(seg1)),Δθ,TimeInterval(t2,t3))

    traj = Trajectory()
    push!(traj,seg1)
    push!(traj,seg2)
    verify(traj)

    get_dist(traj,t1)
    get_dist(traj,t3)
    get_position(traj,t1)
    get_position(traj,t3)
    get_heading(traj,t1)
    get_heading(traj,t3)
    get_vel(traj,t1)
    get_vel(traj,t3)
    get_yaw_rate(traj,t1)
    get_yaw_rate(traj,t3)
end
# cap_trajectory
let
    t1 = 0.0
    t2 = 1.0
    t = 0.5
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    seg1 = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    t3 = 2.0
    center = VecE2(-1.0,1.0)
    radius = 1.0
    θ1 = 0.0
    Δθ = π/2
    θ2 = θ1 + Δθ
    seg2 = ArcTrajectory(center,p2,get_heading(seg1,get_end_time(seg1)),Δθ,TimeInterval(t2,t3))

    traj = Trajectory()
    push!(traj,seg1)
    push!(traj,seg2)
    verify(traj)
    verify(cap_trajectory(traj))
end
# hierarchical Trajectory
let
end
# construct_trajectory
let
    start_pt = VecE2(0.0,0.0)
    start_time = 0.0
    cell_width = 1.0
    transition_time = 1.0

    action_sequence = [WEST,NORTH]
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)
    traj = construct_trajectory(grid_path)
    verify(traj)

    action_sequence = [WEST,WEST,NORTH,NORTH,WAIT,EAST,NORTH,EAST]
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)
    traj = construct_trajectory(grid_path)
    verify(traj)

    action_sequence = [WAIT,WEST,WEST,NORTH,SOUTH,NORTH,WAIT,EAST,NORTH,EAST]
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)
    traj = construct_trajectory(grid_path)
    verify(traj)

    action_sequence = [WEST,NORTH,SOUTH,NORTH,WAIT,EAST,NORTH,EAST,WAIT]
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)
    traj = construct_trajectory(grid_path)
    verify(traj)
end
# DenseTrajectory
let
    t1 = 0.0
    t2 = 1.0
    interval = TimeInterval(t1,t2)
    t = 0.5
    dt = get_Δt(interval,t)
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    seg1 = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    pos_vec = p1+dt*(p2-p1)
    velocity = (p2-p1)/(t2-t1)
    yaw_rate = 0.0

    t3 = 2.0
    t2_mid = (t3+t2)/2
    center = VecE2(0.0,1.0)
    radius = 0.0
    θ1 = 0.0
    Δθ = π/2
    θ2 = θ1 + Δθ
    seg2 = ArcTrajectory(center,p2,get_heading(seg1,get_end_time(seg1)),Δθ,TimeInterval(t2,t3))

    traj_ = Trajectory()
    push!(traj_,seg1)
    push!(traj_,seg2)
    verify(traj_)

    # Construct a DenseTrajectory that encodes the exact same information as the
    # piecewise constant Trajectory
    t_vec = Vector{Float64}()
    vel = Vector{Float64}()
    pos = Vector{Float64}()
    s = 0.0
    for seg in traj_.segments
        push!(t_vec, get_start_time(seg))
        push!(t_vec,   get_end_time(seg))
        push!(vel, norm(get_vel(seg,get_start_time(seg))))
        push!(vel, norm(get_vel(seg,get_end_time(seg))))
        push!(pos, s)
        s += get_length(seg)
        push!(pos, s)
    end
    accel = zeros(length(t_vec)-1)
    traj = DenseTrajectory(traj_,t_vec,accel,vel,pos)

    @test get_time_from_arc_length(traj_,arclength*(t-t1)/(t2-t1)) == t
    # @test get_time_from_arc_length(traj,arclength/2) == t
    @test get_index_time(traj,t) == t

    @test get_start_time(traj) == get_start_time(interval)
    @test get_dist(traj, (t1+t2)/2) == arclength/2

    @test array_isapprox(get_start_pt(traj),p1)
    @test array_isapprox(get_position(traj,t), pos_vec)
    @test array_isapprox(get_heading(traj,t), heading)
    @test array_isapprox(get_vel(traj,t), velocity)
    @test isapprox(get_yaw_rate(traj,t), yaw_rate)

    traj_pt1 = get_trajectory_point_by_time(traj,t)
    @test array_isapprox(traj_pt1.pos, pos_vec)
    @test array_isapprox(traj_pt1.heading, heading)
    @test array_isapprox(traj_pt1.vel, velocity)
    @test isapprox(traj_pt1.yaw_rate, yaw_rate)

    # Now verify that the yaw_rate matches in the zero-velocity segment
    @test isapprox(get_yaw_rate(traj,t2_mid), get_yaw_rate(traj_,t2_mid))
end
# Speed Profile Optimization
let
    start_pt = VecE2(0.0,0.0)
    start_time = 0.0
    action_sequence = [SOUTH,EAST,WAIT,NORTH,EAST,WEST,WEST,SOUTH,SOUTH,WEST,NORTH,WEST,SOUTH,WEST,WEST]
    # action_sequence = [EAST,NORTH,EAST]
    cell_width = 1.0
    transition_time = 2.0
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)

    traj = construct_trajectory(grid_path)
    verify(traj)

    # t_vec, accel, vel, pos = optimize_velocity_profile(traj)
    # dense_traj = DenseTrajectory(traj,t_vec,accel,vel,pos)
    dense_traj, t_vec, accel, vel, pos = optimize_velocity_profile(traj)

    verify(dense_traj)
    t = 0.5
    get_length(dense_traj)
    get_dist(dense_traj,t)
    get_position(dense_traj,t)
    get_heading(dense_traj,t)
    get_vel(dense_traj,t)
    get_yaw_rate(dense_traj,t)

    get_trajectory_point_by_time(dense_traj,t)
end
let
    start_pt = VecE2(0.0,0.0)
    start_time = 0.0
    action_sequence = [WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,NORTH,
        NORTH,NORTH,NORTH,NORTH,NORTH,NORTH,EAST,EAST,EAST,EAST,EAST,EAST,EAST,
        EAST,EAST,EAST,EAST,EAST,EAST,EAST,EAST,EAST,EAST,EAST,NORTH,NORTH,
        NORTH,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,WEST,
        WEST,WEST,WEST,WEST,WEST,SOUTH,SOUTH,SOUTH,SOUTH,SOUTH,SOUTH,SOUTH,
        SOUTH,SOUTH,WEST,WAIT,WAIT,WAIT,SOUTH,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,
        WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,WAIT,
        WAIT,]
    # action_sequence = [EAST,NORTH,EAST]
    cell_width = 0.5
    transition_time = 2.0
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)

    traj = construct_trajectory(grid_path)
    verify(traj)

    # t_vec, accel, vel, pos = optimize_velocity_profile(traj)
    # dense_traj = DenseTrajectory(traj,t_vec,accel,vel,pos)
    dense_traj, t_vec, accel, vel, pos = optimize_velocity_profile(traj;verbose=false)
    @show minimum(accel), maximum(accel)

    verify(dense_traj)
    t = 0.5
    get_length(dense_traj)
    get_dist(dense_traj,t)
    get_position(dense_traj,t)
    get_heading(dense_traj,t)
    get_vel(dense_traj,t)
    get_yaw_rate(dense_traj,t)

    get_trajectory_point_by_time(dense_traj,t)
end
let
    t1 = 0.0
    t2 = 1.0
    t = 0.5
    p1 = VecE2(0.0,0.0)
    p2 = VecE2(0.0,1.0)
    heading = (p2-p1)/norm(p2-p1)
    arclength = norm(p2-p1)
    seg1 = StraightTrajectory(p1,heading,arclength,TimeInterval(t1,t2))

    t3 = 2.0
    center = VecE2(-1.0,1.0)
    radius = 1.0
    θ1 = 0.0
    Δθ = π/2
    θ2 = θ1 + Δθ
    seg2 = ArcTrajectory(center,p2,get_heading(seg1,get_end_time(seg1)),Δθ,TimeInterval(t2,t3))

    traj_ = Trajectory()
    push!(traj_,seg1)
    push!(traj_,seg2)
    verify(traj_)

    Δt = (get_end_time(traj_)-get_start_time(traj_))/100
    SimpleTrajectory(traj_,100)
    traj = SimpleTrajectory(traj_,Δt)
    verify(traj)
    get_trajectory_point_by_time(traj,t)
    get_dist(traj,t)
    get_position(traj,t)
    get_vel(traj,t)
    get_heading(traj,t)
    get_yaw_rate(traj,t)
end
# reverse with a WAIT in the middle
let
    start_pt = VecE2(0.0,0.0)
    start_time = 0.0
    action_sequence = [EAST,WAIT,WEST]
    cell_width = 1.0
    transition_time = 2.0
    grid_path = construct_grid_world_path(start_pt,start_time,
        action_sequence,cell_width,transition_time)

    base_traj = construct_trajectory(grid_path)
    verify(base_traj)
end


text:
"""
    @htl_str -> Result

Create a `Result` object with string interpolation (`\$`) that uses
context-sensitive hypertext escaping. Unlike the `@htl` macro, this
string literal does not include escaping feature [1]. To include `\$`
within user content one must write `$`. Observe that `"` and
any other HTML ampersand escape sequence can be used as appropriate.

In this syntax, interpolation is extended beyond regular Julia strings
to handle three additional cases: tuples, named tuples (for attributes),
and generators. See Julia #38734 for the feature request so that this
could also work within the `@htl` macro syntax.

[1] There are also a few edge cases, see `@raw_str` documentation and
Julia #22926 for more detail. See Julia #38948 for a feature request
that would provide a paired unicode syntax.
"""
macro htl_str(expr::String)
    # Essentially this is an ad-hoc scanner of the string, splitting
    # it by `$` to find interpolated parts and delegating the hard work
    # to `Meta.parse`, treating everything else as a literal string.
    this = Expr(:macrocall, Symbol("@htl_str"), nothing, expr)
    (expr, notations) = extract_notation(expr)
    args = Any[]
    start = idx = 1
    strlen = lastindex(expr)
    while true
        idx = findnext(isequal('$'), expr, start)
        if idx == nothing
           chunk = expr[start:strlen]
           if occursin("htl_sentinel_htl_hack", chunk)
               throw("`htl⟪⟫` notation discovered outside interpolation")
           end
           push!(args, expr[start:strlen])
           break
        end
        push!(args, expr[start:prevind(expr, idx)])
        start = nextind(expr, idx)
        (nest, tail) = Meta.parse(expr, start; greedy=false)
        if nest == nothing
            throw("missing interpolation expression")
        end
        if !(expr[start] == '(' || nest isa Symbol)
            throw(DomainError(nest,
             "interpolations must be symbols or parenthesized"))
        end
        if Meta.isexpr(nest, :(=))
            throw(DomainError(nest,
             "assignments are not permitted in an interpolation"))
        end
        if nest isa String
            # this is an interpolated string literal
            nest = Expr(:string, nest)
        end
        if nest isa Expr
            inject_notation!(nest, notations)
        end
        push!(args, nest)
        start = tail
    end
    return interpolate(args, this)
end

#
# These functions implement Julia #38948 by discovering the notation,
# doing a temporary substitution so that Julia.parse can work, and then
# reinject the subordinate macro application.
#

function extract_notation(expr)::Tuple{String, Vector{String}}
    start = 1
    notations = String[]
    while true
        open_idx = findnext(isequal('⟪'), expr, start)
        done_idx = findnext(isequal('⟫'), expr, start)
        if done_idx != nothing
            if open_idx == nothing || done_idx < open_idx
                throw("unmatched ⟫ delimiter")
            end
        end
        if open_idx == nothing
            return (expr, notations)
        end
        done_idx = notation_extent(expr, open_idx)
        if open_idx < 3 || "htl" !=
           expr[prevind(expr,open_idx,3):prevind(expr,open_idx,1)]
            # this is content, so just move along
            start = nextind(expr, done_idx)
            continue
        end
        push!(notations, expr[nextind(expr,open_idx):prevind(expr,done_idx)])
        swap = "_sentinel_htl_hack($(length(notations)))"
        head = expr[1:prevind(expr,open_idx)]
        tail = expr[nextind(expr, done_idx):end]
        expr = head * swap * tail
        start = open_idx + length(swap)
    end
end

function notation_extent(chunk, idx)
    depth = 0
    strlen = lastindex(chunk)
    while idx <= strlen
        ch = chunk[idx]
        if ch == '⟪'
            depth += 1
        elseif ch == '⟫'
            depth -= 1
        end
        if depth == 0
            return idx
        end
        idx = nextind(chunk, idx)
    end
    throw("unmatched ⟪ delimiter")
end

function inject_notation!(expr::Expr, notations::Vector{String})
    if Meta.isexpr(expr, :call, 2) && expr.args[1] == :htl_sentinel_htl_hack
        content = notations[expr.args[2]]
        expr.head = :macrocall
        expr.args = :(htl"").args
        expr.args[3] = content
        return
    end
    for arg in expr.args
        if arg isa Expr
            inject_notation!(arg, notations)
        end
    end
end


text:
# ---
# title: 72. Edit Distance
# id: problem72
# author: Qling
# date: 2021-03-17
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/edit-distance/description/>
# hidden: true
# ---
#
# Given two strings `word1` and `word2`, return _the minimum number of
# operations required to convert`word1` to `word2`_.
#
# You have the following three operations permitted on a word:
#
#   * Insert a character
#   * Delete a character
#   * Replace a character
#
#
#
# **Example 1:**
#
#
#
#     Input: word1 = "horse", word2 = "ros"
#     Output: 3
#     Explanation:
#     horse -> rorse (replace 'h' with 'r')
#     rorse -> rose (remove 'r')
#     rose -> ros (remove 'e')
#
#
# **Example 2:**
#
#
#
#     Input: word1 = "intention", word2 = "execution"
#     Output: 5
#     Explanation:
#     intention -> inention (remove 't')
#     inention -> enention (replace 'i' with 'e')
#     enention -> exention (replace 'n' with 'x')
#     exention -> exection (replace 'n' with 'c')
#     exection -> execution (insert 'u')
#
#
#
#
# **Constraints:**
#
#   * `0 <= word1.length, word2.length <= 500`
#   * `word1` and `word2` consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

function min_distance(word1::String, word2::String)::Int32
    m, n = length(word1) + 1, length(word2) + 1
    dp = zeros(Int32, m, n)
    dp[:, 1] .= (1:m) .- 1
    dp[1, :] .= (1:n) .- 1

    for i in 2:m, j in 2:n
        dp[i, j] =
            1 + min(
                dp[i - 1, j - 1] - Int(word1[i - 1] == word2[j - 1]),
                dp[i, j - 1],
                dp[i - 1, j],
            )
    end
    return dp[end, end]
end

## @lc code=end


text:
module Clang
    export cindex
    export wrap_c
    export wrap_cpp

    include(joinpath("..", "deps", "ext.jl"))
    include("cindex.jl")
    include("wrap_c.jl")
    include("wrap_cpp.jl")
end


text:
using SafeTestsets

@safetestset "Day03" begin
    using AdventOfCode.TestUtils
    using AdventOfCode.Day03

    @testset "intersection" begin
        using AdventOfCode.Day03: closest_intersection
        @test closest_intersection("R8,U5,L5,D3", "U7,R6,D4,L4") == 6
        @test closest_intersection("R75,D30,R83,U83,L12,D49,R71,U7,L72",
                                   "U62,R66,U55,R34,D71,R55,D58,R83") == 159
        @test closest_intersection("R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51",
                                   "U98,R91,D20,R16,D67,R40,U7,R15,U6,R7") == 135
    end

    @testset "intersection" begin
        using AdventOfCode.Day03: fastest_intersection
        @test fastest_intersection("R8,U5,L5,D3", "U7,R6,D4,L4") == 30
        @test fastest_intersection("R75,D30,R83,U83,L12,D49,R71,U7,L72",
                                   "U62,R66,U55,R34,D71,R55,D58,R83") == 610
        @test fastest_intersection("R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51",
                                   "U98,R91,D20,R16,D67,R40,U7,R15,U6,R7") == 410
    end

    btest("3.1") do
        part1() == 855
    end

    btest("3.2") do
        part2() == 11238
    end
end


text:
# Error codes used in the device API.

"""
Convert a error code to a string for printing purposes.
If the error code is unrecognized, errToStr returns "UNKNOWN".
param errorCode a negative integer return code
return a pointer to a string representing the error
"""
function SoapySDR_errToStr(errorCode)
    ccall((:SoapySDR_errToStr, lib), Cstring, (Cint,), errorCode)
end

struct SoapySDRAPIError
    code::Cint
end
Base.show(io::IO, se::SoapySDRAPIError) = print(io, unsafe_string(SoapySDR_errToStr(se.code)))

"""Returned when read has a timeout."""
const SOAPY_SDR_TIMEOUT = -1

"""Returned for non-specific stream errors."""
const SOAPY_SDR_STREAM_ERROR = -2

"""
Returned when read has data corruption.
For example, the driver saw a malformed packet."""
const SOAPY_SDR_CORRUPTION = -3

"""Returned when read has an overflow condition.
For example, and internal buffer has filled."""
const SOAPY_SDR_OVERFLOW = -4

"""Returned when a requested operation or flag setting
is not supported by the underlying implementation."""
const SOAPY_SDR_NOT_SUPPORTED = -5

"""Returned when a the device encountered a stream time
which was expired (late) or too early to process."""
const SOAPY_SDR_TIME_ERROR = -6

"""Returned when write caused an underflow condition.
For example, a continuous stream was interrupted."""
const SOAPY_SDR_UNDERFLOW = -7


text:
module PlotRegion

using QuickHull
using Combinatorics
using FixedSizeArrays
using Plots

function basic_feasible_point(A::Matrix,b::Vector,set::Vector)
    m,n = size(A)
    @assert length(set) == m "need more indices to define a BFP"
    binds = set # basic variable indices
    ninds = setdiff(1:size(A,2),binds) # non-basic
    B = A[:,binds]
    N = A[:,ninds]
    #cb = c[binds]
    #cn = c[ninds]

    if rank(B) != m
        return (:Infeasible, 0)
    end

    xb = B\b
    x = zeros(eltype(xb),n)
    x[binds] = xb
    x[ninds] = zero(eltype(xb))

    if any(xb .< 0)
        return (:Infeasible, x)
    else
        #lam = B'\cb
        #sn = cn - N'*lam
        return (:Feasible, x)
    end
end

"""
Plot the feasible polytope
Ax = b, x >= 0
for the first two components of x.
"""
function plotregion(A::Matrix,b::Vector)
    m,n = size(A)
    verts = Vector{Vec{2,Float64}}()
    for inds in combinations(1:n,m)
        bfp=basic_feasible_point(A,b,inds)
        if bfp[1] == :Feasible
            push!(verts, Vec(bfp[2][1], bfp[2][2]) )
        end
    end
    hull = qhull(verts)
    plot(Shape(hull),fillalpha=0.5, fillcolor="grey", label="")
end

export plotregion
end


text:
module Data

using Random

include("dataset.jl")
export Dataset,ImageDataset, ArrayDataset, TestDataset, JLDDataset,
        JuliaDBDataset, DFDataset

include("transformer.jl")
export ImageCrop, NoisingTransfomer, onehot, OneHotEncoder,
        MiniBatch, Normalizer, Subset, Split

end


text:
@testset "projection scales" begin
    # inspired by Statutory Valuation of Life Insurance Liabilities, vol 5 chapter 21
    @testset "basic values" begin
        imp = MortalityTables.table("Projection Scale G2 – Male, ANB")
        @test imp[0] == 0.01
        @test imp[105] == 0.000
        @test imp[65] == 0.015
    end

    @testset "improved mortality" begin
        imp = MortalityTables.table("Projection Scale G2 – Male, ANB")
        mort = MortalityTables.table("2012 IAM Period Table – Male, ANB").ultimate

        @test mort[66] * prod(1 .- imp[65:65]) ≈ 8.41978 / 1000
        @test isapprox(mort[85] * prod(1 .- imp[65:84]), 44.60138 / 1000, atol=1e-5)

    end
end


text:
## fixme ** should work
# or add_dirs in project
VW_DIR="/home/graham_s/VirtualWorlds/projects/"

STB_DIR=pwd() # joinpath(VW_DIR,"action_for_children/england/")
for dr in ["src/","test"]
    push!(LOAD_PATH, joinpath(STB_DIR,dr))
end
using Revise


text:
__precompile__(true)

module ArrayViews

using Compat

import Base: eltype, ndims, size, length
import Base: to_index, getindex, setindex!, parent, similar
import Base: Ptr, pointer
import Base: iscontiguous, convert, unsafe_convert
import Base: _sub2ind, _ind2sub

export
    StridedArrayView,
    ArrayView,
    ContiguousView,
    StridedView,
    UnsafeArrayView,
    UnsafeContiguousView,
    UnsafeStridedView,

    ContiguousArray,
    ContiguousVector,
    ContiguousMatrix,

    contiguous_view,
    strided_view,
    aview,
    astride,
    astrides,
    unsafe_aview,
    ellipview,
    diagview,
    rowvec_view,
    flatten_view,
    reshape_view,
    iscontiguous,
    contiguousrank

## source files

include("common.jl")
include("arrviews.jl")
include("layout.jl")
include("indexing.jl")
include("contrank.jl")
include("subviews.jl")
include("convenience.jl")

include("deprecates.jl")

end  # module ArrayViews


text:
function PolynomialFunctionLagrangeForm(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return PolynomialFunctionLagrangeForm((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function degree(obj::PolynomialFunctionLagrangeForm)
    return jcall(obj, "degree", jint, ())
end

function evaluate(::Type{PolynomialFunctionLagrangeForm}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jdouble)
    return jcall(PolynomialFunctionLagrangeForm, "evaluate", jdouble, (Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2)
end

function get_coefficients(obj::PolynomialFunctionLagrangeForm)
    return jcall(obj, "getCoefficients", Vector{jdouble}, ())
end

function get_interpolating_points(obj::PolynomialFunctionLagrangeForm)
    return jcall(obj, "getInterpolatingPoints", Vector{jdouble}, ())
end

function get_interpolating_values(obj::PolynomialFunctionLagrangeForm)
    return jcall(obj, "getInterpolatingValues", Vector{jdouble}, ())
end

function value(obj::PolynomialFunctionLagrangeForm, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end

function verify_interpolation_array(::Type{PolynomialFunctionLagrangeForm}, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::jboolean)
    return jcall(PolynomialFunctionLagrangeForm, "verifyInterpolationArray", jboolean, (Vector{jdouble}, Vector{jdouble}, jboolean), arg0, arg1, arg2)
end


text:
using FunctionalCollections
using Test

import FunctionalCollections: SparseBitmappedTrie, SparseNode, SparseLeaf,
    bitpos, index, hasindex, arrayof, update

@testset "Sparse Bitmapped Vector Tries" begin

    @testset "bitpos and index" begin
        l = SparseLeaf{Int}([6, 11, 16, 21, 26],
                            2^5 | 2^10 | 2^15 | 2^20 | 2^25)
        @test index(l, 6) == 1
        @test index(l, 11) == 2
        @test index(l, 16) == 3
        @test index(l, 21) == 4
        @test index(l, 26) == 5

        @test !hasindex(l, 5)
        @test hasindex(l, 6)
        @test !hasindex(l, 7)
        @test hasindex(l, 11)
        @test !hasindex(l, 12)
        @test hasindex(l, 21)
    end

    @testset "SparseLeaf update" begin
        l = SparseLeaf{Int}()

        @test hasindex(update(l, 1, 1)[1], 1)
        l = update(l, 10, 10)
        @test hasindex(l[1], 10) && !hasindex(l[1], 11)

        l = SparseLeaf{Int}([1, 5], 2^0 | 2^4)
        l, _ = update(l, 2, 2)
        @test arrayof(l) == [1, 2, 5]
        @test length(l) == 3
        @test hasindex(l, 1)
        @test hasindex(l, 2)
        @test !hasindex(l, 3)
        @test !hasindex(l, 4)
        @test hasindex(l, 5)
        @test !hasindex(l, 6)

        @test arrayof(update(l, 2, 100)[1]) == [1, 100, 5]
    end

    @testset "SparseNode update" begin
        n, _ = update(SparseNode(String), 1, "foo")
        @test length(arrayof(n)) == 1

        leaf = arrayof(n)[1].arr[1].arr[1].arr[1].arr[1].arr[1].arr[1]
        @test hasindex(leaf, 1)
        @test leaf.arr[1] == "foo"

        n2, _ = update(n, 33, "bar")
        leaf2 = n2.arr[1].arr[1].arr[1].arr[1].arr[1].arr[1].arr[2]
        @test hasindex(leaf2, 33)
        @test arrayof(leaf2)[1] == "bar"
    end

    @testset "SparseBitmappedTrie get" begin
        n, _ = update(SparseNode(Int), 33, 33)
        @test get(n, 33, "missing") == 33
        @test get(update(SparseNode(Int), 12345, 12345)[1], 12345, "missing") == 12345
        @test get(n, 12345, "missing") == "missing"
    end

    @testset "SparseBitmappedTrie length & items" begin
        n = SparseNode(Int)
        for i=1:1000
            n, _ = update(n, i, i)
        end
        @test lastindex(n) == length(n) == 1000
        @test [i for i=n] == collect(1:1000)
    end

end


text:
using MCMCTempering

struct Joint{Tℓprior, Tℓll}
    ℓprior      :: Tℓprior
    ℓlikelihood :: Tℓll
end

function (joint::Joint)(θ)
    return joint.ℓprior(θ) .+ joint.ℓlikelihood(θ)
end

struct TemperedJoint{Tℓprior, Tℓll, T<:Real} # <: Function
    ℓprior      :: Tℓprior
    ℓlikelihood :: Tℓll
    β           :: Real
end
function TemperedJoint(a,b,c)
    return TemperedJoint{typeof(a), typeof(b), typeof(c)}(a,b,c)
end

function (tj::TemperedJoint)(θ)
    return tj.ℓprior(θ) .+ (tj.ℓlikelihood(θ) .* tj.β)
end

function MCMCTempering.make_tempered_model(
    model::DifferentiableDensityModel,
    β::Real
)
    # @show typeof(model.ℓπ.ℓprior) typeof(model.ℓπ.ℓlikelihood) typeof(β)
    # error()
    ℓπ_β = TemperedJoint(model.ℓπ.ℓprior, model.ℓπ.ℓlikelihood, β)
    ∂ℓπ∂θ_β = TemperedJoint(model.∂ℓπ∂θ.ℓprior, model.∂ℓπ∂θ.ℓlikelihood, β)
    model_β = DifferentiableDensityModel(ℓπ_β, ∂ℓπ∂θ_β)
    return model_β
end

function MCMCTempering.make_tempered_loglikelihood(
    model::DifferentiableDensityModel,
    β::Real
)
    function logπ(z)
        return model.ℓπ.ℓlikelihood(z) * β
    end
    return logπ
end

function MCMCTempering.get_params(trans::AdvancedHMC.Transition)
    return trans.z.θ
end

# Fallback when no random number generator is provided (as is usually the case)
function temperedhmc(system::System, target_accept=0.85, ensemble=MCMCSerial(); kwargs...)
    return temperedhmc(Random.default_rng(), system, target_accept, ensemble; kwargs...)
end

function temperedhmc(
    rng::Random.AbstractRNG,
    system::System,
    target_accept=0.85,
    ensemble=MCMCSerial();
    num_temperatures=4,
    num_chains=1,
    adaptation,
    iterations,
    tree_depth=10,
    discard_initial=adaptation,
    initial_samples=50_000,
    initial_parameters=nothing,
    step_size=nothing,
    progress=true
)

    # Choose parameter dimensionality and initial parameter value
    initial_θ_0 = sample_priors(system)
    D = length(initial_θ_0)

    #ln_prior = make_ln_prior(system)
    ln_prior_transformed = make_ln_prior_transformed(system)
    arr2nt = DirectDetections.make_arr2nt(system)
    priors_vec = _list_priors(system)

    # Capture these variables in a let binding to improve performance
    # ln_post = let system=system, ln_prior=ln_prior, arr2nt=arr2nt
    #     function (θ)
    #         θ_res = arr2nt(θ)
    #         ll = ln_prior(θ) + ln_like(θ_res, system)
    #         return ll
    #     end
    # end
    ℓlikelihood = let system=system, arr2nt=arr2nt, priors_vec=priors_vec
        function (θ_t)
	    #return 0
            # Transform back from the unconstrained support to constrained support for the likelihood function
            θ = Bijectors.invlink.(priors_vec, θ_t)

	#    for i in eachindex(θ)
	#        if !insupport(priors_vec[i], θ[i])
	#	    @warn "likelihood call with value outside the support of the priors"
	#	    return -Inf
	#	end
            #end

	    θ_res = arr2nt(θ)
            return ln_like(θ_res, system)
        end
    end

    # model ii= AdvancedHMC.DifferentiableDensityModel(ℓπ, ForwardDiff)

    # Define the target distribution
    function ℓprior(θ_t)
        θ = Bijectors.invlink.(priors_vec, θ_t)
        ln_prior_transformed(θ)
    end
    # TODO: These can be more efficient by retrieving the value and gradients in a single shot.
    ∂ℓprior∂θ(θ) = (ℓprior(θ), ForwardDiff.gradient(ℓprior, θ))
    ∂ℓlikelihood∂θ(θ) = (ℓlikelihood(θ), ForwardDiff.gradient(ℓlikelihood, θ))

    model = DifferentiableDensityModel(
        Joint(ℓprior, ℓlikelihood),
        Joint(∂ℓprior∂θ, ∂ℓlikelihood∂θ)
    )

    if isnothing(initial_parameters)
        initial_θ = guess_starting_position(system,initial_samples)
    else
        initial_θ = initial_parameters
    end

    # Transform from constrained support to unconstrained support
    initial_θ_t = Bijectors.link.(priors_vec, initial_θ)

    # Define a Hamiltonian system
    metric = DenseEuclideanMetric(D)

    hamiltonian = Hamiltonian(metric, model.ℓπ, model.∂ℓπ∂θ)

    if !isnothing(step_size)
        initial_ϵ = step_size
    else
        initial_ϵ = find_good_stepsize(hamiltonian, initial_θ_t)
        @info "Found initial stepsize" initial_ϵ
    end

    integrator = Leapfrog(initial_ϵ)

    # # We have change some parameters when running with image data
    if !isnothing(system.images) && target_accept > 0.6 && isnothing(step_size)
        @warn "Sampling from images with target accept greater than 0.6. This may lead to insufficient exploration."
    end

    mma = MassMatrixAdaptor(metric)
    if isnothing(step_size)
        @info "Adapting step size and mass matrix"
        ssa = StepSizeAdaptor(target_accept, integrator)
        adaptor = StanHMCAdaptor(mma, ssa)
    else
        @info "Adapting mass matrix only"
        adaptor = MassMatrixAdaptor(metric)
    end

    # κ = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
    κ = NUTS(integrator, max_depth=tree_depth)
    spl = AdvancedHMC.HMCSampler(κ, metric, adaptor)
    tempered_sampler = tempered(
        spl,
        num_temperatures,
        swap_strategy=:nonrev,
        N_swap=5
        # see for more options: https://github.com/TuringLang/MCMCTempering.jl/blob/67ca10916fc9949014475ae2df93488f37dd9d1a/src/tempered.jl#L27
    )
    start_time = time()

    AbstractMCMC.setprogress!(true)
    # mc_samples = sample(
    #     model,
    #     tempered_sampler,
    #     iterations;
    #     discard_initial = adaptation,
    #     progress=progress,
    #     verbose=false
    # )

    mc_samples_all_chains = sample(
        rng,
        model,
        tempered_sampler,
        ensemble,
        iterations+adaptation, num_chains;
        nadapts = adaptation,
        init_params = initial_θ_t,
        discard_initial,
        progress=progress,
        verbose=false
    )
    stop_time = time()

    @info "Sampling compete. Building chains."
    # Go through each chain and repackage results
    chains = MCMCChains.Chains[]
    logposts = Vector{Float64}[]
    for mc_samples in mc_samples_all_chains
        stat = map(s->s.stat, mc_samples)
        logpost = map(s->s.z.ℓπ.value, mc_samples)

        mean_accept = mean(getproperty.(stat, :acceptance_rate))
        num_err_frac = mean(getproperty.(stat, :numerical_error))
        mean_tree_depth = mean(getproperty.(stat, :tree_depth))
        max_tree_depth_frac = mean(getproperty.(stat, :tree_depth) .== tree_depth)

        println("Sampling report:")
        @show mean_accept
        @show num_err_frac
        @show mean_tree_depth
        @show max_tree_depth_frac

        # Report some warnings if sampling did not work well
        if num_err_frac == 1.0
            @error "Numerical errors encountered in ALL iterations. Check model and priors."
        elseif num_err_frac > 0.1
            @warn "Numerical errors encountered in more than 10% of iterations" num_err_frac
        end
        if max_tree_depth_frac > 0.1
            @warn "Maximum tree depth hit in more than 10% of iterations (reduced efficiency)" max_tree_depth_frac
        end

        logpost = map(s->s.z.ℓπ.value, mc_samples)

        # Transform samples back to constrained support
        samples = map(mc_samples) do s
            θ_t = s.z.θ
            θ = Bijectors.invlink.(priors_vec, θ_t)
            return θ
        end
        chain_res = arr2nt.(samples)
        push!(chains, DirectDetections.result2mcmcchain(system, chain_res))
        push!(logposts, logpost)
    end
    # mcmcchains = cat(chains..., dims=3);
    mcmcchains = AbstractMCMC.chainscat(chains...)
    logposts_mat = reduce(hcat, logposts)
    mcmcchains_with_info = MCMCChains.setinfo(
        mcmcchains,
        (;
            start_time,
            stop_time,
            model=system,
            adaptor,
            mc_samples=mc_samples_all_chains,
            tempered_sampler,
            logpost=logposts_mat
        )
    )
    return mcmcchains_with_info
end

export MCMCSerial, MCMCThreads, MCMCDistributed


text:
module SphericalWrapper

using JavaCall

include("OnedWrapper/OnedWrapper.jl")
include("TwodWrapper/TwodWrapper.jl")

end


text:
@testset "Connected components" begin
    @testset "label_components" begin
        @testset "interface" begin
            for T in (Bool, N0f8, Float64, Gray{N0f8})
                A = rand(T, 11, 11)
                @test label_components(A) == label_components(A, strel_diamond(A))
                @test label_components(A; r=2, dims=1) == label_components(A, strel_diamond(A, 1; r=2))

                label = similar(A, Int)
                label_components!(label, A)
                @test label_components(A) == label
                label_components!(label, A; r=2, dims=1)
                @test label_components(A; r=2, dims=1) == label
                label_components!(label, A, strel_box((3, 3)))
                @test label_components(A, strel_box((3, 3))) == label

                Ac = complement.(A)
                @test label_components(Ac; bkg=oneunit(T)) == label_components(A; bkg=zero(T))
            end

            A = rand(Bool, 11, 11)
            se = strel(CartesianIndex, centered(Bool[1 1 1; 1 1 0; 0 0 0]))
            label = @suppress_err label_components(A, se) # deprecated usage in v0.4
            @test label == label_components(A, strel_box((3, 3)))

            msg = "Non-symmetric structuring element is not supported yet"
            se = centered(Bool[1 1 1; 1 1 1; 1 1 0])
            @test_throws ArgumentError(msg) label_components(rand(Bool, 5, 5), se)

            msg = "axes of input and output must match, got (Base.OneTo(11), Base.OneTo(11)) and (-1:1, -1:1). The second argument seems to be a structuring element, it is expected to be the input array."
            @test_throws DimensionMismatch(msg) label_components!(rand(Bool, 11, 11), strel_diamond((3, 3)))

            msg = "labels exhausted, use a larger integer type"
            @test_throws ErrorException(msg) label_components!(zeros(Bool, 4, 4), rand(4, 4))
        end

        A = [
            true  true  false true
            true  false true  true
        ]
        lbltarget = [
            1 1 0 2
            1 0 2 2
        ]
        lbltarget1 = [
            1 2 0 4
            1 0 3 4
        ]
        @test label_components(A) == lbltarget
        @test label_components(A; dims=:) == lbltarget
        @test label_components(A; dims=1) == lbltarget1
        connectivity = centered([
            false true  false
            true  false true
            false true  false
        ])
        @test label_components(A, connectivity) == lbltarget
        connectivity = centered(trues(3, 3))
        lbltarget2 = [
            1 1 0 1
            1 0 1 1
        ]
        @test label_components(A, connectivity) == lbltarget2
    end

    @testset "component_boxes" begin
        A = rand(0:5, 11, 11)
        label = label_components(A, strel_box((5, 5)))
        boxes = component_boxes(label)
        @test eltype(boxes) <: CartesianIndices
        @test axes(boxes) == (0:maximum(label),)

        # the bounding box is the smallest region that contains the label value
        @test all(eachindex(boxes)) do i
            n1 = count(isequal(i), label[boxes[i]])
            n2 = count(isequal(i), label)
            n1 == n2
        end

        A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
        label = label_components(A)
        boxes = component_boxes(label)
        @test boxes == OffsetArray(
            [
                CartesianIndices((2:5, 2:4)),
                CartesianIndices((1:1, 1:5)),
                CartesianIndices((2:5, 1:3)),
                CartesianIndices((3:5, 3:5)),
                CartesianIndices((2:3, 4:5)),
            ], -1)
        @test A[boxes[1]] == [2 2 2 2 2]
        @test A[boxes[4]] == [0 1; 1 1]

        label = zeros(Int, 5, 5)
        label[1] = -1
        msg = "The input labeled array should contain background label `0` as the minimum value"
        @test_throws ArgumentError(msg) component_boxes(label)

        boxes = component_boxes(ones(Int, 5, 5))
        @test length(boxes) == 1
    end

    @testset "component_indices" begin
        A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
        label = label_components(A)
        indices = component_indices(label)
        @test IndexStyle(label) == IndexLinear()
        @test eltype(indices) == Vector{Int}
        @test axes(indices) == (0:maximum(label),)
        @test indices == OffsetVector([
                [8, 10, 17], [1, 6, 11, 16, 21], [2, 3, 4, 5, 7, 9, 12],
                [13, 14, 15, 19, 20, 24, 25], [18, 22, 23],
            ], -1)

        @test indices == component_indices(Int, label) == component_indices(IndexLinear(), label)

        indices = component_indices(CartesianIndex, label)
        @test eltype(indices) == Vector{CartesianIndex{ndims(A)}}
        @test axes(indices) == (0:maximum(label),)
        @test indices == [CartesianIndices(label)[R] for R in component_indices(Int, label)]
        @test indices == component_indices(IndexCartesian(), label)

        plabel = PaddedView(0, label, (0:6, 0:6))
        indices = component_indices(CartesianIndex, label)
        pindices = component_indices(plabel)
        @test IndexStyle(plabel) == IndexCartesian()
        @test eltype(pindices) == Vector{CartesianIndex{ndims(A)}}
        @test pindices[1:end] == indices[1:end] # they only differ in background indices

        label = zeros(Int, 5, 5)
        label[1] = -1
        msg = "The input labeled array should contain background label `0` as the minimum value"
        @test_throws ArgumentError(msg) component_indices(label)

        indices = component_indices(ones(Int, 5, 5))
        @test axes(indices, 1) == 0:1
        @test isempty(indices[0])
    end

    @testset "component_lengths" begin
        A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
        label = label_components(A)
        counts = component_lengths(label)
        indices = component_indices(label)
        @test counts == length.(indices)
        @test axes(counts) == (0:maximum(label),)

        label = zeros(Int, 5, 5)
        label[1] = -1
        msg = "The input labeled array should contain background label `0` as the minimum value"
        @test_throws ArgumentError(msg) component_lengths(label)

        counts = component_lengths(ones(Int, 5, 5))
        @test axes(counts, 1) == 0:1
        @test counts[0] == 0
    end

    @testset "component_centroids" begin
        A = [2 2 2 2 2; 1 1 1 0 1; 1 0 2 1 1; 1 1 2 2 2; 1 0 2 2 2]
        label = label_components(A)
        centers = component_centroids(label)
        indices = component_indices(CartesianIndex, label)
        @test centers == [sum(X).I ./ length(X) for X in indices]
        @test axes(centers) == (0:maximum(label),)

        label = zeros(Int, 5, 5)
        label[1] = -1
        msg = "The input labeled array should contain background label `0` as the minimum value"
        @test_throws ArgumentError(msg) component_centroids(label)

        centers = component_centroids(ones(Int, 5, 5))
        @test axes(centers, 1) == 1:1
    end

    @testset "deprecations" begin
        A = [
            true  true  false true
            true  false true  true
        ]
        label = label_components(A)
        indices = @suppress_err component_subscripts(label)
        @test indices ==
            Array{Tuple}[[(2, 2), (1, 3)], [(1, 1), (2, 1), (1, 2)], [(2, 3), (1, 4), (2, 4)]]
    end

end


text:
__precompile__()
module POMCP

if VERSION >= v"0.6.0"
    warn("The POMCP package is deprecated. Please use BasicPOMCP (https://github.com/JuliaPOMDP/BasicPOMCP.jl) instead.")
end

import POMDPs

import POMDPs: action, solve
import Base: rand, mean
import POMDPs: update, updater, initialize_belief
import POMDPToolbox
import StatsBase: WeightVec, sample
import MCTS: RandomActionGenerator, RolloutEstimator, next_action, estimate_value, init_N, convert_to_policy
using ParticleFilters

using Compat

export
    AbstractPOMCPSolver,
    POMCPSolver,
    POMCPDPWSolver,
    POMCPPlanner,
    RootUpdater,
    BeliefNode,
    RootNode,
    ObsNode,
    AbstractActNode,
    ActNode,
    POMCPTreeVisualizer,

    solve,
    action,
    create_policy,
    update,
    updater,
    create_belief,
    init_V,
    init_N,
    sparse_actions,
    estimate_value,
    extract_belief,
    convert_estimator,
    add_N,

    RandomActionGenerator,
    next_action,

    RolloutEstimator,
    PORollout,
    FORollout,
    FOValue,

    ParticleCollection,
    ParticleReinvigorator,
    reinvigorate!,
    handle_unseen_observation,
    DefaultReinvigoratorStub,

    NoDecision,
    AllSamplesTerminal,
    ExceptionRethrow,
    default_action,

    # deprecated
    PORolloutEstimator

include("tree.jl")
include("particle_filter.jl")
include("exceptions.jl")

abstract type AbstractPOMCPSolver <: POMDPs.Solver end

"""
POMCP Solver type

Fields:

    eps::Float64
        Rollouts and tree expansion will stop when discount^depth is less than this.
        default: 0.01

    max_depth::Int
        Rollouts and tree expension will stop when this depth is reached.
        default: 10

    c::Float64
        UCB exploration constant - specifies how much the solver should explore.
        default: 1.0

    tree_queries::Int
        Number of iterations during each action() call.
        default: 100

    rng::AbstractRNG
        Random number generator.
        default: Base.GLOBAL_RNG

    node_belief_updater::Updater
        Calculates the belief for a new belief node (see notebooks/Belief_and_Particle_Filter_Options.ipynb for more info.)
        default: DefaultReinvigoratorStub() - this will simply keep the particles as described in the paper without doing any reinvigoration.

    estimate_value::Any (rollout policy can be specified by setting this to RolloutEstimator(policy))
        Function, object, or number used to estimate the value at the leaf nodes.
        If this is a function `f`, `f(pomdp, s, h::BeliefNode, steps)` will be called to estimate the value.
        If this is an object `o`, `estimate_value(o, pomdp, s, h::BeliefNode, steps)` will be called.
        If this is a number, the value will be set to that number
        default: RolloutEstimator(RandomSolver(rng))

    init_V::Any
        Function, object, or number used to set the initial V(h,a) value at a new node.
        If this is a function `f`, `f(pomdp, h, a)` will be called to set the value.
        If this is an object `o`, `init_V(o, pomdp, h, a)` will be called.
        If this is a number, V will be set to that number
        default: 0.0

    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(pomdp, h, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, pomdp, h, a)` will be called.
        If this is a number, N will be set to that number
        default: 0

    num_sparse_actions::Int
        Number of actions to be considered at each node.
        If <= 0, the entire action space will be considered.
        default: 0

    default_action::Any
        Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
        If this is a Function `f`, `f(belief, ex)` will be called.
        If this is a Policy `p`, `action(p, belief)` will be called.
        If it is an object `a`, `default_action(a, belief, ex) will be called, and
        if this method is not implemented, `a` will be returned directly.
"""
type POMCPSolver <: AbstractPOMCPSolver
    eps::Float64 # will stop simulations when discount^depth is less than this
    max_depth::Int
    c::Float64 # UCB exploration constant
    tree_queries::Int
    rng::AbstractRNG
    node_belief_updater::Union{POMDPs.Updater, DefaultReinvigoratorStub}

    estimate_value::Any

    init_V::Any
    init_N::Any

    num_sparse_actions::Int # = 0 or less if not used
    default_action::Any
end

"""
POMCP Solver type

Fields:

    eps::Float64
        Rollouts and tree expansion will stop when discount^depth is less than this.
        default: 0.01

    max_depth::Int
        Rollouts and tree expension will stop when this depth is reached.
        default: 10

    c::Float64
        UCB exploration constant - specifies how much the solver should explore.
        default: 1.0

    tree_queries::Int
        Number of iterations during each action() call.
        default: 100

    rng::AbstractRNG
        Random number generator.
        default: Base.GLOBAL_RNG

    node_belief_updater::Updater
        Calculates the belief for a new belief node (see notebooks/Belief_and_Particle_Filter_Options.ipynb for more info.)
        default: DefaultReinvigoratorStub() - this will simply keep the particles as described in the paper without doing any reinvigoration.

    estimate_value::Any (rollout policy can be specified by setting this to RolloutEstimator(policy))
        Function, object, or number used to estimate the value at the leaf nodes.
        If this is a function `f`, `f(pomdp, s, h::BeliefNode, steps)` will be called to estimate the value.
        If this is an object `o`, `estimate_value(o, pomdp, s, h::BeliefNode, steps)` will be called.
        If this is a number, the value will be set to that number
        default: RolloutEstimator(RandomSolver(rng))

    enable_action_pw::Bool
        Turn progressive widening of the number of actions considered on or off.
        If false, all actions will be considered.
        default: true

    k_action::Float64
    alpha_action::Float64
    k_observation::Float64
    alpha_observation::Float64
        These constants control the double progressive widening. A new observation
        or action will be added if the number of children is less than or equal to kN^alpha.
        defaults: k:10, alpha:0.5

    init_V::Any
        Function, object, or number used to set the initial V(h,a) value at a new node.
        If this is a function `f`, `f(pomdp, h, a)` will be called to set the value.
        If this is an object `o`, `init_V(o, pomdp, h, a)` will be called.
        If this is a number, V will be set to that number
        default: 0.0

    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(pomdp, h, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, pomdp, h, a)` will be called.
        If this is a number, N will be set to that number
        default: 0

    next_action::Any
        Function or object used to choose the next action to be considered for progressive widening.
        The next action is determined based on the POMDP, the belief, `b`, and the current `BeliefNode`, `h`.
        If this is a function `f`, `f(pomdp, b, h)` will be called to set the value.
        If this is an object `o`, `next_action(o, pomdp, b, h)` will be called.
        default: RandomActionGenerator(rng)

    default_action::Any
        Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
        If this is a Function `f`, `f(belief, ex)` will be called.
        If this is a Policy `p`, `action(p, belief)` will be called.
        If it is an object `a`, `default_action(a, belief, ex) will be called, and
        if this method is not implemented, `a` will be returned directly.

For more information on the k and alpha parameters, see Couëtoux, A., Hoock, J.-B., Sokolovska, N., Teytaud, O., & Bonnard, N. (2011). Continuous Upper Confidence Trees. In Learning and Intelligent Optimization. Rome, Italy. Retrieved from http://link.springer.com/chapter/10.1007/978-3-642-25566-3_32
"""
type POMCPDPWSolver <: AbstractPOMCPSolver
    eps::Float64 # will stop simulations when discount^depth is less than this
    max_depth::Int
    c::Float64
    tree_queries::Int
    rng::AbstractRNG
    node_belief_updater::Union{POMDPs.Updater, DefaultReinvigoratorStub}

    estimate_value::Any

    enable_action_pw::Bool

    alpha_observation::Float64
    k_observation::Float64
    alpha_action::Float64
    k_action::Float64
    init_V::Any
    init_N::Any
    next_action::Any
    default_action::Any
end

"""
Policy that builds a POMCP tree to determine an optimal next action.

Note, you should construct this using the create_policy function
"""
type POMCPPlanner{S, A, O, B, SolverType<:AbstractPOMCPSolver} <: POMDPs.Policy
    problem::POMDPs.POMDP{S,A,O}
    solver::SolverType
    node_belief_updater::POMDPs.Updater{B}
    solved_estimate::Any

    #XXX hack
    _tree_ref::Nullable{Any}
end

include("constructor.jl")
include("updater.jl")
include("actions.jl")
include("rollout.jl")
include("solver.jl")
include("visualization.jl")

end # module


text:
struct Entity{S,D,I} # state, definition, identification
    state::S
    def::D
    id::I
end
Entity(entity::Entity{S,D,I}, s::S) where {S,D,I} = Entity(s, entity.def, entity.id)

Base.:(==)(a::Entity{S,D,I}, b::Entity{S,D,I}) where {S,D,I} = a.state == b.state && a.def == b.def && a.id == b.id

max_stars_repo_name: tawheeler/Records

max_issues_repo_name: tawheeler/Records

max_forks_repo_name: tawheeler/Records


text:
@doc raw"""
    add_variable(psi_container::PSIContainer,
                      devices::D,
                      var_name::Symbol,
                      binary::Bool,
                      expression_name::Symbol,
                      sign::Float64)

Adds a variable to the optimization model and to the affine expressions contained
in the psi_container model according to the specified sign. Based on the inputs, the variable can
be specified as binary.

# Bounds

``` lb_value_function <= varstart[name, t] <= ub_value_function ```

If binary = true:

``` varstart[name, t] in {0,1} ```

# LaTeX

``  lb \ge x^{device}_t \le ub \forall t ``

``  x^{device}_t \in {0,1} \forall t iff \text{binary = true}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* devices : Vector or Iterator with the devices
* var_name::Symbol : Base Name for the variable
* binary::Bool : Select if the variable is binary
* expression_name::Symbol : Expression_name name stored in psi_container.expressions to add the variable
* sign::Float64 : sign of the addition of the variable to the expression_name. Default Value is 1.0

# Accepted Keyword Arguments
* ub_value_function : Provides the function over device to obtain the value for a upper_bound
* lb_value_function : Provides the function over device to obtain the value for a lower_bound. If the variable is meant to be positive define lb = x -> 0.0
* initial_value_function : Provides the function over device to obtain the warm start value

"""
function add_variable(
    psi_container::PSIContainer,
    devices::D,
    var_name::Symbol,
    binary::Bool,
    expression_name::Union{Nothing, Symbol} = nothing,
    sign::Float64 = 1.0;
    kwargs...,
) where {D <: Union{Vector{<:PSY.Device}, IS.FlattenIteratorWrapper{<:PSY.Device}}}
    time_steps = model_time_steps(psi_container)
    variable = add_var_container!(
        psi_container,
        var_name,
        (PSY.get_name(d) for d in devices),
        time_steps,
    )

    lb_f = get(kwargs, :lb_value, nothing)
    init_f = get(kwargs, :initial_value, nothing)
    ub_f = get(kwargs, :ub_value, nothing)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        variable[name, t] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "$(var_name)_{$(name), $(t)}",
            binary = binary
        )

        !isnothing(ub_f) && JuMP.set_upper_bound(variable[name, t], ub_f(d))
        !isnothing(lb_f) && !binary && JuMP.set_lower_bound(variable[name, t], lb_f(d))
        !isnothing(init_f) && JuMP.set_start_value(variable[name, t], init_f(d))

        if !(isnothing(expression_name))
            bus_number = PSY.get_number(PSY.get_bus(d))
            add_to_expression!(
                get_expression(psi_container, expression_name),
                bus_number,
                t,
                variable[name, t],
                sign,
            )
        end
    end

    return

end

@doc raw"""
    set_variable_bounds!(
        psi_container::PSIContainer,
        bounds::DeviceRange,
        var_type::AbstractString,
        device_type::Type{PSY.Device},
    )

Adds a bounds to a variable in the optimization model.

# Bounds

``` bounds.min <= varstart[name, t] <= bounds.max  ```

# LaTeX

``  x^{device}_t >= bound^{min;} \forall t ``

``  x^{device}_t <= bound^{max} \forall t ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* bounds::DeviceRange : contains names and vector of min / max
* var_type::AbstractString : type of the variable
* T: type of the device

"""
function set_variable_bounds!(
    psi_container::PSIContainer,
    bounds::Vector{DeviceRange},
    var_type::AbstractString,
    ::Type{T},
) where {T <: PSY.Device}
    var = get_variable(psi_container, var_type, T)
    for t in model_time_steps(psi_container), bound in bounds
        _var = var[bound.name, t]
        JuMP.set_upper_bound(_var, bound.limits.max)
        JuMP.set_lower_bound(_var, bound.limits.min)
    end
end