text:
using Random

const SEED = 999
const NROWS = 128
const N_IDX_TESTS = 16

arrow_tempname = tempname()
push!(temps, arrow_tempname)

Random.seed!(SEED)

randdate() = Date(rand(0:4000), rand(1:12), rand(1:27))
randtime() = Dates.Time(rand(0:23), rand(0:59), rand(0:59))
randdatetime() = randdate() + randtime()

randstrings() = String[[randstring(rand(0:20)) for i ∈ 1:(NROWS-1)]; "a"]
function randstrings(::Missing)
    Union{String,Missing}[[rand(Bool) ? missing : randstring(rand(0:20)) for i ∈ 1:(NROWS-1)]; "a"]
end

convstring(str::AbstractString) = String(str)
convstring(::Missing) = missing

@testset "ArrowTests" begin

cols = [rand(Int32,NROWS),
    rand(Float64,NROWS),
    Date[randdate() for i ∈ 1:NROWS],
    DateTime[randdatetime() for i ∈ 1:NROWS],
    Dates.Time[randtime() for i ∈ 1:NROWS],
    Union{Int64,Missing}[rand(Bool) ? missing : rand(Int64) for i ∈ 1:NROWS],
    randstrings(),
    randstrings(missing),
    CategoricalArrays.categorical(randstrings()),
    CategoricalArrays.categorical(randstrings(missing))]

colnames = [:ints,:floats,:dates,:datetimes,:times,:missingints,:strings,
    :missingstrings,:catstrings,:catstringsmissing]

featherwrite(arrow_tempname, cols, colnames)

ndf = featherread(arrow_tempname)

@test ndf.names == colnames

@test typeof(ndf.columns[1]) == Arrow.Primitive{Int32}
@test typeof(ndf.columns[2]) == Arrow.Primitive{Float64}
@test typeof(ndf.columns[3]) == Arrow.Primitive{Arrow.Datestamp}
@test typeof(ndf.columns[4]) == Arrow.Primitive{Arrow.Timestamp{Dates.Millisecond}}
@test typeof(ndf.columns[5]) == Arrow.Primitive{Arrow.TimeOfDay{Dates.Nanosecond,Int64}}
@test typeof(ndf.columns[6]) == Arrow.NullablePrimitive{Int64}
@test typeof(ndf.columns[7]) == Arrow.List{String,Arrow.DefaultOffset,Arrow.Primitive{UInt8}}
@test typeof(ndf.columns[8]) == Arrow.NullableList{String,Arrow.DefaultOffset,Arrow.Primitive{UInt8}}
@test typeof(ndf.columns[9]) == Arrow.DictEncoding{String,Arrow.Primitive{Int32},
    Arrow.List{String,Arrow.DefaultOffset,Arrow.Primitive{UInt8}}}
@test typeof(ndf.columns[10]) ==
        Arrow.DictEncoding{Union{String,Missing},Arrow.NullablePrimitive{Int32},Arrow.List{String,Arrow.DefaultOffset,
        Arrow.Primitive{UInt8}}}

for j ∈ 1:N_IDX_TESTS
    i = rand(1:NROWS)
    @test cols[1][i] == ndf.columns[1][i]
    @test cols[2][i] == ndf.columns[2][i]
    @test cols[3][i] == convert(Date, ndf.columns[3][i])
    @test cols[4][i] == convert(DateTime, ndf.columns[4][i])
    @test cols[5][i] == convert(Dates.Time, ndf.columns[5][i])
    @test isequal(cols[6][i], ndf.columns[6][i])
    @test cols[7][i] == ndf.columns[7][i]
    @test isequal(cols[8][i], ndf.columns[8][i])
    @test cols[9][i] == String(ndf.columns[9][i])
    @test isequal(cols[10][i], convstring(ndf.columns[10][i]))
end
for j ∈ 1:N_IDX_TESTS
    a, b = extrema(rand(1:NROWS, 2))
    i = a:b
    @test cols[1][i] == ndf.columns[1][i]
    @test cols[2][i] == ndf.columns[2][i]
    @test cols[3][i] == convert.(Date, ndf.columns[3][i])
    @test cols[4][i] == convert.(DateTime, ndf.columns[4][i])
    @test cols[5][i] == convert.(Dates.Time, ndf.columns[5][i])
    @test isequal(cols[6][i], ndf.columns[6][i])
    @test cols[7][i] == ndf.columns[7][i]
    @test isequal(cols[8][i], ndf.columns[8][i])
    @test cols[9][i] == String.(ndf.columns[9][i])
    @test isequal(cols[10][i], convstring.(ndf.columns[10][i]))
end

end


text:
using MPIMapReduce
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
sz = MPI.Comm_size(comm)

@time y1 = MPIMapReduce.pmapreduce(x -> ones(2, 2) * x, *, 1:4);

if MPI.Comm_rank(comm) == 0
    show(stdout, MIME"text/plain"(), y1)
    println()
end

@time y2 = MPIMapReduce.pmapreduce(x -> ones(2, 2) * x, MPI.PROD, 1:sz);

if MPI.Comm_rank(comm) == 0
    show(stdout, MIME"text/plain"(), y2)
    println()
end

ext: jl


text:
using DrWatson
@quickactivate "OptimalEstimationProject"
using MATLAB
using IndirectTrajOpt
using StaticArrays
using LinearAlgebra

function main()
# Data info
tf      = 13.0
folder  = "numParticles2000_tf"*replace(string(tf), '.' => 'p')
trial   = 10
hStep   = 2

# Load data
df = readTextData(datadir(folder))

# Init parameters
ps = initCR3BPIndirectParams("Low Thrust 10 CR3BP")
μ = ps.crp.μ
LU = ps.crp.LU
TU = ps.crp.TU

# Trajectory time span
tspan = (0.0, tf * 86400 / TU)

# Initial conditions
ics = [-0.019488511458668, -0.016033479812051, 0.0,
        8.918881923678198, -4.081793688818725, 0.0,
        1.0]

# Target final conditions
xf = @SVector [(1.0 - μ) + 6060.483/LU, 19452.284/LU, -34982.968/LU, 0.082677*(TU/LU), 0.006820*(TU/LU), -0.368434*(TU/LU)]

# Integrate target trajectory
solTarg = integrate(xf, (0.0, 1.6), "Low Thrust 10 CR3BP", CR3BP(), FullSolutionHistoryNoControl())
rxTarg = [u[1] for u in solTarg.u]
ryTarg = [u[2] for u in solTarg.u]
rzTarg = [u[3] for u in solTarg.u]

# Grab requirements from DataFrame
λ0 = df.HomotopySolutions[trial][hStep]
ps.ϵ = (hStep == 1 ? 1.0 : df.HomotopyParams[trial][hStep - 1])

# Initialize final integration inputs
y0 = @SVector [ics[1], ics[2], ics[3], ics[4], ics[5], ics[6], ics[7],
            λ0[1], λ0[2], λ0[3], λ0[4], λ0[5], λ0[6], λ0[7]]

# Integrate
sol = integrate(y0, tspan, ps, CR3BP(), FullSolutionHistory(), MEMF())

# Compute thrusting and coasting arcs
cSc = ps.sp.c * ps.crp.TU / (ps.crp.LU * 1000.0)
rxt = Vector{Float64}(undef, length(sol)); fill!(rxt, NaN)
ryt = Vector{Float64}(undef, length(sol)); fill!(ryt, NaN)
rzt = Vector{Float64}(undef, length(sol)); fill!(rzt, NaN)
rxc = Vector{Float64}(undef, length(sol)); fill!(rxc, NaN)
ryc = Vector{Float64}(undef, length(sol)); fill!(ryc, NaN)
rzc = Vector{Float64}(undef, length(sol)); fill!(rzc, NaN)
Sold = 0.0
for i in 1:length(sol)
    λv = norm(@view(sol.u[i][11:13]))
    S = IndirectTrajOpt.computeS(sol.u[i], λv, cSc)
    #println(S)
    i == 1 ? Sold = S : ()
    if S > ps.ϵ || Sold > ps.ϵ
        rxc[i] = sol.u[i][1]
        ryc[i] = sol.u[i][2]
        rzc[i] = sol.u[i][3]
    end
    if S <= ps.ϵ || Sold <= ps.ϵ
        rxt[i] = sol.u[i][1]
        ryt[i] = sol.u[i][2]
        rzt[i] = sol.u[i][3]
    end
    Sold = S
end

mat"""
figure
plot3($rxc,$ryc,$rzc, 'b')
hold on
plot3($rxt,$ryt,$rzt, 'r')
plot3($rxTarg, $ryTarg, $rzTarg, '--k')
axis equal
grid on
"""
end
main()

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
function testJackknife(numtest, MCS, obs)
  inform("Start testing Jackknife / $obs")

  nsigma = 3
  counts = zeros(Int, nsigma)

  for itest in 1:numtest
    sobs = obs()
    for mcs in 1:MCS
      sobs << randn()
    end
    jk = jackknife(sobs)

    @test_approx_eq_eps mean(sobs) mean(jk) 1.0e-10
    @test_approx_eq_eps stderror(sobs) stderror(jk) 1.0e-10
  end
  inform("Success testing Jackknife / $obs")
end

testJackknife(1_000, 8192, SimpleObservable)
testJackknife(1_000, 8192, BinningObservable)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SurfelsLib_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SurfelsLib")
JLLWrappers.@generate_main_file("SurfelsLib", UUID("63138847-3b8b-53b7-bc56-9d476e6e0bd2"))
end  # module SurfelsLib_jll


text:
using EchogramPlots
using Test

a = rand(100,100)
echogram(a)

max_stars_repo_name: RobBlackwell/EchogramPlots

max_issues_repo_name: RobBlackwell/EchogramPlots

max_forks_repo_name: RobBlackwell/EchogramPlots


text:
using TextParse
"""
    standartize(x...)

    returned standartized matrices with mean and standard deviation calculated jointly across all datasets

"""
function standartize(x...)
    xx = hcat(x...)
    mask = std(xx, dims = 2)[:] .> 1e-6
    x = map(z -> z[mask, :], x)

    xx = hcat(x...)
    mn = mean(xx, dims = 2);
    sd = std(xx, dims = 2)
    map(z -> (z .- mn)./sd, x)
end

"""
    catwithlabels(x...)

    horizontally concatenate matrices `x` while return label vector, starting with first matrix with one

"""
catwithlabels(x...) = hcat(x...), mapreduce(i -> i[1]*ones(Int,i[2]),vcat,enumerate(size.(x,2)))

"""
    typedread(filename,T,transposed = true)

    load CSV file such that all columns are initiated to type T
"""
typedread(filename,T) = T.(Matrix(hcat(TextParse.csvread(filename,' ',header_exists=false,spacedelim = true)[1]...)'))

loaddataset(name,difficulty,idir,T=Float32) = (typedread(joinpath(idir,name,"normal.txt"),T),typedread(joinpath(idir,name,difficulty*".txt"),T))

"""
    trData, tstData, clusterdness = makeset(dataset, alpha, frequency, variation, [normalize, seed])

Sample a given dataset, return training and testing subsets and a measure of clusterdness.
See Emmott, Andrew F., et al. "Systematic construction of anomaly detection benchmarks from
real data.", 2013 for details.

alpha - the ratio of training to all data\n
difficulty - easy/medium/hard/very_hard problem based on similarity of anomalous measurements to normal\n
variation - low/high - should anomalies be clustered or not\n
seed - random seed
"""
function makeset(normal, anomalous, alpha, variation, seed=time_ns())
    # test correct parameters size
    !(0 <= alpha <= 1) && error("alpha must be in the interval [0,1]")

    # problem dimensions
    n = size(normal, 2)
    trn_n = Int(round(n*alpha))
    tst_n = n - trn_n

    # how many anomalous points to be sampled
    a = size(anomalous, 2)
    trn_a = Int(round(a*alpha))
    tst_a = a - trn_a

    # set seed
    Random.seed!(seed)
    # normalize the data to zero mean and unit variance
    normal, anomalous = standartize(normal, anomalous)

    # randomly sample the training and testing normal data
    inds = randperm(n)[1:trn_n]
    trn_n_data = normal[:, inds]
    tst_n_data = normal[:, setdiff(1:n,inds)]

    # now sample the anomalous data
    if variation == "low"
        # in this setting, simply sample trn_a anomalous points
        inds = randperm(a)[1:trn_a]
    elseif variation == "high"
        error("not impleemnted yet")
        # # in this setting, randomly choose a point and then trn_a-1 nearest points to it as a cluster
        # x = anomalous[:, sample(1:a)]
        # x = reshape(x, : , 1) # promote the vector to a 2D array
        # # here maybe other metrics could be used?
        # dists = pairwise(Euclidean(), x, anomalous) # get the distance vector
        # inds = sortperm(reshape(dists, length(dists))) # get the sorted indices
        # inds = inds[1:trn_a] # get the nearest ones
        # inds = inds[randperm(length(inds))] # scramble them
    end
    trn_a_data = anomalous[:, inds]
    tst_a_data = anomalous[:,setdiff(1:a,inds)]

    # compute the clusterdness - sample variance of normal vs anomalous instances
    varN = mean(pairwise(Euclidean(), trn_n_data[:, sample(1:size(trn_n_data, 2), min(1000, size(trn_n_data, 2)), replace=false)]))/2
    varA = (trn_n > 0) ? mean(pairwise(Euclidean(), trn_a_data[:, sample(1:size(trn_a_data, 2), min(1000, size(trn_a_data, 2)), replace=false)]))/2 : 0.0

    clusterdness =  (varA>0) ? varN/varA : clusterdness = Inf
    catwithlabels(trn_n_data, trn_a_data), catwithlabels(tst_n_data, tst_a_data), clusterdness
end

"""
    subsampleanomalous(x, α, seed = time_ns())
    subsampleanomalous(x, n, seed = time_ns())

    removes all but α-fraction (or n) non-anomalous samples from dataset `x = (data,labels)`. Unless
    there is no negative sample and `α` is zero, at least positive (anomalous) sample will be returned
"""
subsampleanomalous(x,α::AbstractFloat,seed = time_ns()) = subsampleanomalous(x,Int(ceil(α*sum(x[2] .== 1)/0.95)),seed)
function subsampleanomalous(x::Tuple{A,B},n::Int,seed = time_ns()) where {A<:AbstractMatrix, B<:Vector}
    inds =  _subsampleanomalous(x[2], n, seed)
    x[1][:,inds], x[2][inds]
end

function subsampleanomalous(x::Tuple{A,B},n::Int,seed = time_ns()) where {A<:AbstractVector, B<:Vector}
    inds =  _subsampleanomalous(x[2], n, seed)
    x[1][inds], x[2][inds]
end

function _subsampleanomalous(labels, n, seed)
    n == 0 && findall(labels .== 1)
    a = findall(labels .> 1)
    inds = sample(a,min(n,length(a)),replace=false)
    vcat(findall(labels .== 1), inds)
end

"""
    subsampletrainvalidation(x, α, β, seed = time_ns())
    subsampletrainvalidation(x, na, nl, seed = time_ns())

    return training and validation dataset, such that training set is poluted with `α` anomalous data
    and `β` legitimate data.
"""
function subsampletrainvalidation(x::Tuple, na, nl, seed = time_ns())
    (trn_idxs, val_idxs) = _subsampletrainvalidation(x[2], na, nl, seed)
    return((x[1][:, trn_idxs], x[2][trn_idxs]),(x[1][:,val_idxs], x[2][val_idxs]))
end

function _subsampletrainvalidation(labels, na, nl, seed)
    ai = na == 0 ? Vector{Int}() :  sampleindices(findall(labels .> 0), na)
    li = sampleindices(findall(labels .== 0), nl)
    trn_idxs = vcat(ai,li)
    val_idxs = setdiff(1:length(labels), trn_idxs)
    return(trn_idxs, val_idxs)
end

sampleindices(indices, α::AbstractFloat) = sampleindices(indices, round(Int,α*length(indices)))
function sampleindices(indices, n::Int)
    n = min(n, length(indices))
    return(sample(indices, n ,replace=false))
end


text:
#common, except for Jaeckel
abstract type SRSolver end

function normalizePrice(isCall::Bool, price::T, f::T, strike::T, df::T) where {T}
    c = price / f / df
    ex = f / strike

    if !isCall
        if ex <= 1
            c = c + 1 - 1 / ex # put call parity
        else
            #else duality + put call parity
            c = ex * c
            ex = 1 / ex
        end
    else
        if ex > 1
            # use c(-x0, v)
            c = (f * (c - 1) + strike) / strike # in out duality, c = ex*c + 1 - ex  //not as good numericall
            ex = 1 / ex
        end
    end
    return c, ex
end

function impliedVolatilitySRGuess(isCall::Bool, price::Real, f::Real, strike::Real, tte::Real, df::Real)::Real
    c,ex = normalizePrice(isCall, price, f, strike, df)
    scaledVol = impliedVolatilitySRGuessUndiscountedCall(c, ex, log(ex))
    return scaledVol / sqrt(tte)
end

# f = 1.0, strike = 1/ex, ex=ey = 1.0/strike=f/strike, y = log(f/strike) = log(ey), alpha = price*ex
#polya is 2 / pi, polya_factor = sqrt(pi/8) is aludaat,
function impliedVolatilitySRGuessUndiscountedCall(price::Real, ey::Real, y::Real, polya_factor::Real = 2/pi)::Real
    alpha = price *ey
    r = 2 * alpha - ey + 1
    em2piy = exp(-polya_factor * y)
    A = (ey * em2piy - 1.0 / (ey * em2piy))^2
    r2 = r * r
    B = 4 * (em2piy + 1.0 / em2piy) - 2 / ey * (ey * em2piy + 1.0 / (ey * em2piy)) * ((ey * ey) + 1 - r2)
    beta = 0.0
    if abs(alpha) < SqrtEpsilon
        C = -16 * (1 - 1.0 / ey) * alpha - 16 * (1 - 3 / ey + 1 / (ey * ey)) * alpha * alpha
        beta = C / B
        #    println(" alpha ",alpha, " C ",C," beta ",beta)
    elseif abs(ey - alpha) < SqrtEpsilon
        C = -16 * (1 + 1.0 / ey) * (alpha - ey) - 16 * (1 + 3 / ey + 1 / (ey * ey)) * (alpha - ey) * (alpha - ey)
        if C == 0
            beta = B / A #this is wrong sign but gives order of magnitude
        else
            beta = C / B
        end
        #    println(" alpha ",alpha, " C ",C," beta ",beta)
    elseif abs(y) < SqrtEpsilon
        a2 = alpha * alpha
        a4 = a2 * a2
        C = 16 * ((a2 - a4) + (2 * a4 + 2 * alpha * a2 - a2 - alpha) * y)
        B = 16 * (a2 - y * (alpha + a2))
        beta = C / B
    else
        eym1 = ey - 1
        eyp1 = ey + 1
        C = 1.0 / (ey * ey) * (r2 - eym1 * eym1) * (eyp1 * eyp1 - r2)
        beta = 2 * C / (B + sqrt(B * B + 4 * A * C))
        # println("beta ",beta, " ",C, " ",B, " ",A)
    end
    gamma = -log(beta) / polya_factor
    if y >= 0
        Asqrty = 0.5 * (1 + sqrt(1 - em2piy * em2piy))
        c0 = (Asqrty - 0.5/ey)
        gmy = max(gamma - y,0)      #machine epsilon issues
        if price <= c0
            return sqrt(gamma + y) - sqrt(gmy)
        end
        return (sqrt(gamma + y) + sqrt(gmy))
    end
    Asqrty = 0.5 * (1 - sqrt(1 - 1.0 / (em2piy * em2piy)))
    c0 = (0.5 - Asqrty/ey)
    gpy = gamma + y
    if gpy < 0
        #machine epsilon issues
        gpy = 0.0
    end
    if price <= c0
        return (-sqrt(gpy) + sqrt(gamma - y))
    end
    return (sqrt(gpy) + sqrt(gamma - y))
end


text:
using DebuggingUtilities
using Test

function foo()
    x = 5
    @showln x
    x = 7
    @showln x
end

# The more laborious approach, in case one needs speed
function foofl()
    x = 5
    @showfl(@__FILE__, @__LINE__, x)
    x = 7
    @showfl(@__FILE__, @__LINE__, x)
end

io = IOBuffer()
DebuggingUtilities.showlnio[] = io

@test foo() == 7

str = chomp(String(take!(io)))
target = ("x = 5", "(in foo at", "x = 7", "(in foo at")
for (i,ln) in enumerate(split(str, '\n'))
    ln = lstrip(ln)
    @test startswith(ln, target[i])
end

io = IOBuffer()
DebuggingUtilities.showlnio[] = io
@test foofl() == 7
str = chomp(String(take!(io)))
target = ("x = 5", "(at file ", "x = 7", "(at file ")
for (i,ln) in enumerate(split(str, '\n'))
    ln = lstrip(ln)
    @test startswith(ln, target[i])
end

# Just make sure these run
io = IOBuffer()
DebuggingUtilities.showlnio[] = io
test_showline("noerror.jl")
@test_throws DomainError test_showline("error.jl")
time_showline("noerror.jl")

DebuggingUtilities.showlnio[] = stdout
nothing


text:
struct MWSL3DIntegrand{C,O,L}
    test_triangular_element::C
    trial_triangular_element::C
    op::O
    test_local_space::L
    trial_local_space::L
end

const i4pi = 1 / (4pi)
function (igd::MWSL3DIntegrand)(u,v)
        α = igd.op.α
        β = igd.op.β
        γ = igd.op.gamma

        x = neighborhood(igd.test_triangular_element,u)
        y = neighborhood(igd.trial_triangular_element,v)

        r = cartesian(x) - cartesian(y)
        R = norm(r)
        iR = 1 / R
        G = exp(-γ*R)*(i4pi*iR)

        f = igd.test_local_space(x)
        g = igd.trial_local_space(y)

        j = jacobian(x) * jacobian(y)

        jG = j*G
        αjG = α*jG
        βjG = β*jG

        G = @SVector [αjG*g[1].value, αjG*g[2].value, αjG*g[3].value]
        H = @SVector [βjG*g[1].divergence, βjG*g[2].divergence, βjG*g[3].divergence]

        SMatrix{3,3}((
            dot(f[1].value,G[1]) + f[1].divergence*H[1],
            dot(f[2].value,G[1]) + f[2].divergence*H[1],
            dot(f[3].value,G[1]) + f[3].divergence*H[1],
            dot(f[1].value,G[2]) + f[1].divergence*H[2],
            dot(f[2].value,G[2]) + f[2].divergence*H[2],
            dot(f[3].value,G[2]) + f[3].divergence*H[2],
            dot(f[1].value,G[3]) + f[1].divergence*H[3],
            dot(f[2].value,G[3]) + f[2].divergence*H[3],
            dot(f[3].value,G[3]) + f[3].divergence*H[3]))
end

struct MWDL3DIntegrand{C,O,L}
        test_triangular_element::C
        trial_triangular_element::C
        op::O
        test_local_space::L
        trial_local_space::L
end

function (igd::MWDL3DIntegrand)(u,v)

        γ = igd.op.gamma

        x = neighborhood(igd.test_triangular_element,u)
        y = neighborhood(igd.trial_triangular_element,v)

        r = cartesian(x) - cartesian(y)
        R = norm(r)
        iR = 1/R
        G = exp(-γ*R)*(iR*i4pi)

        f = igd.test_local_space(x)
        g = igd.trial_local_space(y)

        j = jacobian(x) * jacobian(y)

        GG = -(γ + iR) * G * (iR * r)

        G = @SVector [cross(GG,(j*g[1].value)), cross(GG,(j*g[2].value)), cross(GG,(j*g[3].value))]
        SMatrix{3,3}((
            dot(f[1].value, G[1]),
            dot(f[2].value, G[1]),
            dot(f[3].value, G[1]),
            dot(f[1].value, G[2]),
            dot(f[2].value, G[2]),
            dot(f[3].value, G[2]),
            dot(f[1].value, G[3]),
            dot(f[2].value, G[3]),
            dot(f[3].value, G[3])))
end

kernel_in_bary(op::MWSingleLayer3D,
    test_local_space::RTRefSpace, trial_local_space::RTRefSpace,
    test_chart, trial_chart) = MWSL3DIntegrand(
        test_chart, trial_chart, op, test_local_space, trial_local_space)

kernel_in_bary(op::MWDoubleLayer3D,
    test_local_space::RTRefSpace, trial_local_space::RTRefSpace,
    test_chart, trial_chart) = MWDL3DIntegrand(
        test_chart, trial_chart, op, test_local_space, trial_local_space)

const MWOperator3D = Union{MWSingleLayer3D, MWDoubleLayer3D}
function momintegrals!(op::MWOperator3D,
    test_local_space::RTRefSpace, trial_local_space::RTRefSpace,
    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

    I, J, K, L = SauterSchwabQuadrature.reorder(
        test_triangular_element.vertices,
        trial_triangular_element.vertices, strat)

    test_triangular_element  = simplex(
        test_triangular_element.vertices[I[1]],
        test_triangular_element.vertices[I[2]],
        test_triangular_element.vertices[I[3]])

    trial_triangular_element = simplex(
        trial_triangular_element.vertices[J[1]],
        trial_triangular_element.vertices[J[2]],
        trial_triangular_element.vertices[J[3]])

    # igd = MWSL3DIntegrand(test_triangular_element, trial_triangular_element,
    #     op, test_local_space, trial_local_space)
    igd = kernel_in_bary(op, test_local_space, trial_local_space,
        test_triangular_element, trial_triangular_element)
    G = sauterschwab_parameterized(igd, strat)
    for j ∈ 1:3, i ∈ 1:3
        out[i,j] += G[K[i],L[j]]
    end

    nothing
end

function momintegrals_nested!(op::MWOperator3D,
    test_local_space::RTRefSpace, trial_local_space::RTRefSpace,
    test_chart, trial_chart, out, strat::SauterSchwabStrategy)

    # 1. Refine the trial_chart
    p1, p2, p3 = trial_chart.vertices

    # TODO: generalise this to include more general refinements
    e1 = cartesian(neighborhood(trial_chart, (0,1/2)))
    e2 = cartesian(neighborhood(trial_chart, (1/2,0)))
    e3 = cartesian(neighborhood(trial_chart, (1/2,1/2)))

    ct = cartesian(center(trial_chart))

    refined_trial_chart = [
        simplex(ct, p1, e3),
        simplex(ct, e3, p2),
        simplex(ct, p2, e1),
        simplex(ct, e1, p3),
        simplex(ct, p3, e2),
        simplex(ct, e2, p1)]

    qs = defaultquadstrat(op, test_local_space, trial_local_space)
    qd = quaddata(op, test_local_space, trial_local_space,
        [test_chart], refined_trial_chart, qs)

    for (q,chart) in enumerate(refined_trial_chart)
        qr = quadrule(op, test_local_space, trial_local_space,
            1, test_chart, q ,chart, qd, qs)

        Q = restrict(trial_local_space, trial_chart, chart)
        zlocal = zero(out)
        momintegrals!(op, test_local_space, trial_local_space,
            test_chart, chart, zlocal, qr)

        for j in 1:3
            for i in 1:3
                for k in 1:3
                out[i,j] += zlocal[i,k] * Q[j,k]
        end end end
    end
end

function momintegrals_nested!(op::IntegralOperator,
    test_local_space::RefSpace, trial_local_space::RefSpace,
    test_chart, trial_chart, out, strat)

    momintegrals!(op, test_local_space, trial_local_space,
        test_chart, trial_chart, out, strat)
end


text:
"""
	Join partecipant common and extra dataset
"""
function join_partecipant_common_and_extra_datasets(common_dataset::String,extra_dataset::String)
	df_common = CSV.File(common_dataset,missingstring="NA") |> DataFrame
	df_extra = CSV.File(extra_dataset,missingstring="NA") |> DataFrame

	df_result = @from r1 in df_common begin
				@join r2 in df_extra on r1.part_id equals r2.part_id
				@select {r1.part_id,r1.hh_id,r1.part_age,r1.part_gender,r2.part_occupation}
				@collect DataFrame
	end

	replace!(df_result.part_occupation, JOB_SEEKING => UNEMPLOYED,  missing => 6) # We handle joob seeking people as the ones staying at home
	set_occupation_for_record_without(df_result)

	CSV.write("resources/processed_partecipant_extra.csv",df_result)

	return df_result
end

function join_contact_common_and_extra_datasets(common_dataset::String,extra_dataset::String)
	df_common = CSV.File(common_dataset,missingstring="NA") |> DataFrame
	df_extra = CSV.File(extra_dataset,missingstring="NA") |> DataFrame

	df_result = @from r1 in df_common begin
				@join r2 in df_extra on r1.cont_id equals r2.cont_id
				@select {r1.cont_id,r1.part_id,r1.cnt_age_exact,r1.cnt_age_est_min,r1.cnt_age_est_max,r1.cnt_gender,
				r1.cnt_home,r1.cnt_work,r1.cnt_school,r1.cnt_transport,r1.cnt_leisure,r1.cnt_otherplace,
				r1.frequency_multi,r1.phys_contact,r1.duration_multi,r2.cnt_otherplace_family,r2.cnt_otherplace_grandparents,r2.cnt_hh_member}
				@collect DataFrame
	end

	CSV.write("resources/processed_contact_extra.csv",df_result)

	return df_result
end

"""
	Join partecipant and contact datasets, and gets some metadata
"""
function analyze_contact_data(partecipand_dataset::String,partecipant_extra_dataset::String,contact_dataset::String,contact_extra_dataset::String)
	df_partecipant = join_partecipant_common_and_extra_datasets(partecipand_dataset,partecipant_extra_dataset)
	df_contact = join_contact_common_and_extra_datasets(contact_dataset,contact_extra_dataset)

	#TODO: collapse all this replaces into a single function
	replace!(df_contact.cnt_age_exact, missing => -1)
	replace!(df_contact.cnt_age_est_max, missing => -1)
	replace!(df_contact.cnt_age_est_min, missing => -1)
	replace!(df_contact.cnt_hh_member, missing => "N")
	replace!(df_partecipant.part_occupation, RETIRED => UNEMPLOYED) #Replace "retired" with "home"

	df_contact.cnt_home = convert.(Int,df_contact.cnt_home)
	df_contact.cnt_work = convert.(Int,df_contact.cnt_work)
	df_contact.cnt_school = convert.(Int,df_contact.cnt_school)
	df_contact.cnt_transport = convert.(Int,df_contact.cnt_transport)
	df_contact.cnt_leisure = convert.(Int,df_contact.cnt_leisure)
	df_contact.cnt_otherplace = convert.(Int,df_contact.cnt_otherplace)

	#Calculate age for the records that has at least the age interval
	transform!(df_contact, [:cnt_age_exact,:cnt_age_est_min,:cnt_age_est_max] => ByRow(calculate_age) => :age)
	#Remove the record for which is impossible to infer the age
	df_contact = df_contact |> @filter(_.age > 0) |> DataFrame

	df_result = @from r1 in df_partecipant begin
				@join r2 in df_contact on r1.part_id equals r2.part_id
				@select {r1.part_id,r1.hh_id,r2.cont_id,r1.part_age,r2.age,r1.part_gender,r1.part_occupation,r2.cnt_gender,
						 r2.cnt_home,r2.cnt_work,r2.cnt_school,r2.cnt_transport,r2.cnt_leisure,
						 r2.cnt_otherplace,r2.frequency_multi,r2.phys_contact,r2.duration_multi,
						 r2.cnt_otherplace_family,r2.cnt_otherplace_grandparents,r2.cnt_hh_member}
				@collect DataFrame
	end

	CSV.write("resources/processed_partecipant_contact.csv",df_result)

	#### Generate occupation and household for each contact
	occupation_distribution = get_occupation_distribution(df_partecipant)

	transform!(df_result,[:age] => (x -> get_occupation(x,occupation_distribution)) => :cnt_occupation,
						 [:part_id,:cont_id,:hh_id,:cnt_hh_member,:cnt_otherplace_family,:cnt_otherplace_grandparents] => ByRow(get_household) => :cnt_household_id)
	CSV.write("resources/processed_contact_with_occupation.csv",df_result)

	### Append the resulted contact dataset to the partecipant one
	df_contact_to_partecipant = df_result |> @select(:cont_id,:cnt_household_id,:cnt_occupation,:cnt_gender,:age) |>
										   	 @rename(:cont_id => :part_id, :cnt_household_id => :hh_id, :age => :part_age, :cnt_occupation => :part_occupation, :cnt_gender => :part_gender) |>
										   	 DataFrame
	replace!(df_contact_to_partecipant.part_gender, missing => "NA")

	append!(df_partecipant,df_contact_to_partecipant)
	CSV.write("resources/processed_partecipant_and_contact.csv",df_partecipant)

	### Normalize indicies
	count_id = -1
	df_partecipant.part_id = map(i -> begin count_id+=1; return count_id end, df_partecipant.part_id)
	CSV.write("resources/processed_partecipant_and_contact_normalized.csv",df_partecipant)

end

function get_household(part_id,cnt_id,household_id,hh_member,family_cnt,grandparents_cnt)
	if hh_member == "Y" || family_cnt || grandparents_cnt
		return household_id
	else
		return string(part_id,cnt_id) # The concatenation of the 2 ids is used as key
	end
end

function get_occupation(age,occupation_distribution)
	result = []

	for i in 1:length(age)
		push!(result,get_occupation_single_person(age[i],occupation_distribution))
	end
	return result
end

"""
	Get occupation based on the occupation vector
	Special value return for age < CHILD_AGE && age > ELDER_AGE
"""
function get_occupation_single_person(age::Int, occupation_distribution)
	# Handle children
	if age < CHILD_AGE
		if age < INFANT
			return UNEMPLOYED
		else
			return SCHOOL
		end
	end

	# Handle elder
	if age >= ELDER_AGE
		return UNEMPLOYED
	end

	index = floor(Int,age / 10) - 1
	calculated_prob_sum = rand()
	prob_sum = 0

	for i in 1:length(occupation_distribution[index])
		prob_sum += occupation_distribution[index][i][:num]
		#Handle cases where the float sum is not 1
		if i == length(occupation_distribution[index])
			return occupation_distribution[index][i][:Key]
		end
		if calculated_prob_sum <= prob_sum
			return occupation_distribution[index][i][:Key]
		end
	end
end

"""
	Returns a DataFrame containing the number of contacts for each class and for each age interval
"""
function get_contact_age_distribution(df::DataFrame)
	min_age, max_age = extrema(df.age)
	df_result = DataFrame(min_age = Int[],max_age = Int[],size = Int[],home_perc = Float32[],work_perc= Float32[],
				school_perc= Float32[],transport_perc= Float32[],leisure_perc= Float32[],otherplace_perc= Float32[])

	for i in range(0,_next_10_multiple(max_age) - 10,step=10)
		size,cnt_vector = get_contact_distribution_per_age(df,i,i+10)

		perc = 100/size
		cnt_vector = perc .* cnt_vector
		result = vcat(i,i+10,size,cnt_vector)

		push!(df_result,result)
	end

	return df_result
end

"""
	Gets, for the given age interval, the number of contacts for each class
"""
function get_contact_distribution_per_age(df::DataFrame, min_age::Int, max_age::Int)
	df_result = @from i in df begin
				@where i.age > min_age && i.age < max_age
				@select i
				@collect DataFrame
				end

	cnt_home = sum(df_result.cnt_home)
	cnt_work = sum(df_result.cnt_work)
	cnt_school = sum(df_result.cnt_school)
	cnt_transport = sum(df_result.cnt_transport)
	cnt_leisure = sum(df_result.cnt_leisure)
	cnt_otherplace = sum(df_result.cnt_otherplace)

	cnt_vector = [cnt_home,cnt_work,cnt_school,cnt_transport,cnt_leisure,cnt_otherplace]

	return sum(cnt_vector),cnt_vector
end

function get_occupation_distribution(df::DataFrame,occupation_filter::Function)
	age_distribution = Array{Array{Any},1}()

	for i in range(CHILD_AGE,_next_10_multiple(ELDER_AGE) - 10,step=10)
		push!(age_distribution, get_occupation_distribution_per_age(df,i,i+10,occupation_filter))
	end

	return age_distribution
end

"""
	Return a multidimensional vector, where the position i contains the occupation distribution in the range i*10 - (i*10 + 10)
"""
function get_occupation_distribution(df::DataFrame)
	return get_occupation_distribution(df, x -> true)
end

function get_occupation_distribution_per_age(df::DataFrame,min_age::Int,max_age::Int,occupation_filter::Function)
	df_result = df |> @filter(_.part_age >= min_age && _.part_age < max_age && occupation_filter(_.part_occupation)) |>
					  @groupby(_.part_occupation) |>
					  @map({ Key = key(_), num = length(_)}) |>
					  @orderby(_.Key) |>
					  DataFrame

	df_result[!,:num] = df_result[!,:num] ./ sum(df_result[!,:num])

	return Tables.rowtable(df_result)
end

"""
	Assign an occupation at each record with part_occupation == 6 (without one)
"""
function set_occupation_for_record_without(df::DataFrame)
	occupation_distribution = get_occupation_distribution(df, x -> x != 6 ? true : false)

	filter = (age,occupation) -> occupation == 6 ?  get_occupation_single_person(age,occupation_distribution) : occupation

	transform!(df,[:part_age,:part_occupation] => ByRow(filter) => :part_occupation)
end

function get_dataframe_subset(df::DataFrame, num_rows::Int, occupation_proportion::Array{Int},seed::Int)
	if sum(occupation_proportion) != 100
		throw(ArgumentError("occupation proportion must cover all the rows space"))
	end

	num_workers = occupation_proportion[1] * num_rows / 100
	num_students = occupation_proportion[2] * num_rows / 100
	num_unemployeds = occupation_proportion[3] * num_rows / 100

	shuffled_indices = randperm(MersenneTwister(seed),size(df,1))

	students_added = 0
	workers_added = 0
	unemployeds_added = 0
	subset_rows = []
	for i in shuffled_indices
		record = df[i,:]
		if length(subset_rows) == num_rows
			break
		end

		if is_student(record) && students_added < num_students
			push!(subset_rows,i)
			students_added += 1
		elseif is_worker(record) && workers_added < num_workers
			push!(subset_rows,i)
			workers_added += 1
		elseif is_home(record) && unemployeds_added < num_unemployeds
			push!(subset_rows,i)
			unemployeds_added += 1
		end
	end

	if length(subset_rows) != num_rows
		throw(ArgumentError("The size or the proportion of the records cannot be applied to this dataframe.
										Try to reduce the number of rows or to change the target proportions"))
	end

	return df[subset_rows,:]
end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
export dtwplot

"""
    dtwplot(seq1, seq2, [dist=SqEuclidean()])
    dtwplot(seq1, seq2, D, i1, i2)

Given two sequences, perform dynamic time warping and plot
the results. If alignment has already been computed, pass
the indices `i1` and `i2` to make the plot.
"""
function handleargs(seq1, seq2, dist::SemiMetric = SqEuclidean(); kwargs...)
    D = dtw_cost_matrix(seq1, seq2, dist; kwargs...)
    cost, i1, i2 = DynamicAxisWarping.trackback(D)
    seq1, seq2, D, i1, i2
end

handleargs(h; kwargs...) = handleargs(h.args...; kwargs...)

@userplot DTWPlot

@recipe function f(h::DTWPlot; transportcost=1)
    seq1, seq2, D, i1, i2 = handleargs(h; transportcost=transportcost)

    n1, n2 = lastlength(seq1), lastlength(seq2)

    all = ndims(seq1) == 1
    # set up the subplots
    legend --> false
    link := :both
    grid --> false
    if all
        layout --> @layout [
            a b{0.9w,0.9h}
            _ c
        ]
    else
        layout --> 1
    end

    left_margin --> 0mm
    bottom_margin --> 0mm
    top_margin --> 0mm
    right_margin --> 0mm
    clims --> (0, 3 * D[end, end])

    # heatmap
    @series begin
        seriestype := :heatmap
        formatter --> (z) -> ""
        subplot := (all ? 2 : 1)
        D
    end

    # the rest of the plots are paths
    seriestype := :path
    linecolor --> RGB(0, 0, 0)

    # main plot
    s1 = @series begin
        linewidth --> 3
        subplot := (all ? 2 : 1)
        formatter --> (z) -> ""
        i1, i2
    end

    if all
        # left line plot
        @series begin
            subplot := 1
            seq2, 1:n2
        end

        # bottom line plot
        @series begin
            subplot := 3
            1:n1, seq1
        end
    end
end

@userplot MatchPlot

znorm(x) = (x = x.- mean(x); x ./= std(x))
using Statistics
@recipe function f(h::MatchPlot; transportcost=1, separation=2, ds=1)
    x, y, D, i1, i2 = handleargs(h; transportcost=transportcost)
    x,y = znorm.((x,y))
    s1 = x .- separation
    s2 = y .+ separation

    @series begin
        s1
    end
    @series begin
        s2
    end
    @series begin
        primary := false
        linecolor --> :black
        alpha --> 0.2
        i = fill(Inf, 1, length(i1))
        vec([i1'; i2'; i][:,1:ds:end]), vec([s1[i1]'; s2[i2]'; i][:,1:ds:end])
    end
end


text:
module MiscTest

using ReTest
using MiniLoggers
using MiniLoggers: getmode, NoTransformations, Squash

using Dates

@testset "flush" begin
    logger = MiniLogger(flush_threshold = 1000)
    @test logger.flush_threshold[] == 1000

    logger = MiniLogger(flush_threshold = Second(1))
    @test logger.flush_threshold[] == 1000
end

mutable struct MockIORecord
    buf::String
    flushed::Bool
end
MockIORecord(buf) = MockIORecord(buf, false)

mutable struct MockIO <: IO
    recs::Vector{MockIORecord}
end
MockIO() = MockIO(MockIORecord[])
function Base.flush(io::MockIO)
    foreach(io.recs) do rec
        rec.flushed = true
    end
end
Base.isopen(io::MockIO) = true
Base.write(io::MockIO, a::Vector{UInt8}) = push!(io.recs, MockIORecord(String(copy(a))))

@testset "delayed flush" begin
    io = MockIO()
    logger = MiniLogger(io = io, flush_threshold = 100, format = "{message}")
    sleep(0.2)
    with_logger(logger) do
        @info "Foo"
        @info "Bar"
    end
    @test length(io.recs) == 2
    @test io.recs[1].buf == "Foo\n"
    @test io.recs[2].buf == "Bar\n"
    @test io.recs[1].flushed
    @test !io.recs[2].flushed

    sleep(0.2)
    with_logger(logger) do
        @info "Baz"
    end
    @test length(io.recs) == 3
    @test all(x -> x.flushed, io.recs)
    @test io.recs[3].buf == "Baz\n"
end

@testset "file logs" begin
    fpath = joinpath(@__DIR__, "logtest1.log")
    fpath2 = joinpath(@__DIR__, "logtest2.log")
    logger = MiniLogger(io = fpath, ioerr = fpath2, format = "{message}")
    with_logger(logger) do
        @info "Foo"
    end
    open(fpath, "r") do f
        @test read(f, String) == "Foo\n"
    end
    close(logger)
    @test !isopen(logger.io)
    @test !isopen(logger.ioerr)

    logger = MiniLogger(io = fpath, format = "{message}", append = true)
    with_logger(logger) do
        @info "Bar"
    end
    open(fpath, "r") do f
        @test read(f, String) == "Foo\nBar\n"
    end
    close(logger)
    @test !isopen(logger.io)
    @test isopen(logger.ioerr)

    logger = MiniLogger(io = fpath, ioerr = stdout, format = "{message}")
    with_logger(logger) do
        @info "Baz"
    end
    open(fpath, "r") do f
        @test read(f, String) == "Baz\n"
    end
    close(logger)
    @test !isopen(logger.io)
    @test isopen(logger.ioerr)

    try
        rm(fpath, force = true)
        rm(fpath2, force = true)
    catch err
        @error "" exception = (err, catch_backtrace())
    end
end

@testset "modes" begin
    @test getmode("notransformations") isa NoTransformations
    @test getmode("squash") isa Squash
    @test getmode(:notransformations) isa NoTransformations
    @test getmode(:squash) isa Squash
end

end # module


text:
module RecurrentNN
import Base.tanh

export Model, RNN, LSTM, GRU, GFLSTM, GFGRU
export NNMatrix, randNNMat, forwardprop, softmax, Solver, step
export Graph, backprop, rowpluck

include("recurrent.jl")
include("graph.jl")
include("solver.jl")
include("rnn.jl")
include("lstm.jl")
include("gru.jl")
include("gflstm.jl")
include("gfgru.jl")

end # module

max_issues_repo_name: Andy-P/RecurrentNN

max_forks_repo_name: Andy-P/RecurrentNN


text:
@testset "AbstractZero" begin
    @testset "iszero" begin
        @test iszero(ZeroTangent())
        @test iszero(NoTangent())
    end

    @testset "Linear operators" begin
        @test view(ZeroTangent(), 1) == ZeroTangent()
        @test view(NoTangent(), 1, 2) == NoTangent()

        @test sum(ZeroTangent()) == ZeroTangent()
        @test sum(NoTangent(); dims=2) == NoTangent()

        @test reshape(NoTangent(), (1, :)) === NoTangent()
        @test reshape(ZeroTangent(), 2, 3, 4, 5) === ZeroTangent()
    end

    @testset "ZeroTangent" begin
        z = ZeroTangent()
        @test z + z === z
        @test z + 1 === 1
        @test 1 + z === 1
        @test z - z === z
        @test z - 1 === -1
        @test 1 - z === 1
        @test -z === z
        @test z * z === z
        @test z * 11.1 === ZeroTangent()
        @test 12.3 * z === ZeroTangent()
        @test dot(z, z) === z
        @test dot(z, 1.8) === z
        @test dot(2.1, z) === z
        @test dot([1, 2], z) === z
        @test dot(z, [1, 2]) === z
        for x in z
            @test x === z
        end
        @test broadcastable(z) isa Ref{ZeroTangent}
        @test zero(@thunk(3)) === z
        @test zero(NoTangent()) === z
        @test zero(ZeroTangent) === z
        @test zero(NoTangent) === z
        @test zero(Tangent{Tuple{Int,Int}}((1, 2))) === z
        for f in (transpose, adjoint, conj)
            @test f(z) === z
        end
        @test z / 2 === z / [1, 3] === z

        @test eltype(z) === ZeroTangent
        @test eltype(ZeroTangent) === ZeroTangent

        # use mutable objects to test the strong `===` condition
        x = ones(2)
        @test muladd(ZeroTangent(), 2, x) === x
        @test muladd(2, ZeroTangent(), x) === x
        @test muladd(ZeroTangent(), ZeroTangent(), x) === x
        @test muladd(2, 2, ZeroTangent()) === 4
        @test muladd(x, ZeroTangent(), ZeroTangent()) === ZeroTangent()
        @test muladd(ZeroTangent(), x, ZeroTangent()) === ZeroTangent()
        @test muladd(ZeroTangent(), ZeroTangent(), ZeroTangent()) === ZeroTangent()

        @test reim(z) === (ZeroTangent(), ZeroTangent())
        @test real(z) === ZeroTangent()
        @test imag(z) === ZeroTangent()

        @test complex(z) === z
        @test complex(z, z) === z
        @test complex(z, 2.0) === Complex{Float64}(0.0, 2.0)
        @test complex(1.5, z) === Complex{Float64}(1.5, 0.0)
        @test Complex(z, 2.0) === Complex{Float64}(0.0, 2.0)
        @test Complex(1.5, z) === Complex{Float64}(1.5, 0.0)
        @test ComplexF64(z, 2.0) === Complex{Float64}(0.0, 2.0)
        @test ComplexF64(1.5, z) === Complex{Float64}(1.5, 0.0)

        @test convert(Bool, ZeroTangent()) === false
        @test convert(Int64, ZeroTangent()) === Int64(0)
        @test convert(Float32, ZeroTangent()) === 0.0f0
        @test convert(ComplexF64, ZeroTangent()) === 0.0 + 0.0im
    end

    @testset "NoTangent" begin
        dne = NoTangent()
        @test dne + dne == dne
        @test dne + 1 == 1
        @test 1 + dne == 1
        @test dne - dne == dne
        @test dne - 1 == -1
        @test 1 - dne == 1
        @test -dne == dne
        @test dne * dne == dne
        @test dne * 11.1 == dne
        @test 12.1 * dne == dne
        @test dot(dne, dne) == dne
        @test dot(dne, 17.2) == dne
        @test dot(11.9, dne) == dne

        @test ZeroTangent() + dne == dne
        @test dne + ZeroTangent() == dne
        @test ZeroTangent() - dne == dne
        @test dne - ZeroTangent() == dne

        @test ZeroTangent() * dne == ZeroTangent()
        @test dne * ZeroTangent() == ZeroTangent()
        @test dot(ZeroTangent(), dne) == ZeroTangent()
        @test dot(dne, ZeroTangent()) == ZeroTangent()

        for x in dne
            @test x === dne
        end
        @test broadcastable(dne) isa Ref{NoTangent}
        for f in (transpose, adjoint, conj)
            @test f(dne) === dne
        end
        @test dne / 2 === dne / [1, 3] === dne

        @test convert(Int64, NoTangent()) == 0
        @test convert(Float64, NoTangent()) == 0.0
    end
end


text:
# the purpose of this script is to fix the wrongly saved value of zdim model files
using DrWatson
@quickactivate
using BSON, FileIO
using Flux
using GenerativeModels
using DistributionsAD
using ValueHistories
using GenerativeAD
using ArgParse
using ProgressMeter

s = ArgParseSettings()
@add_arg_table! s begin
	"path"
		arg_type = String
		help = "path"
	"--force"
		action = :store_true
		help = "force overwriting"
end
parsed_args = parse_args(ARGS, s)
@unpack path, force = parsed_args

path = abspath(path)
files = GenerativeAD.Evaluation.collect_files(path)
mfiles = filter(f->occursin("model", f), files)
sfiles = filter(f->!occursin("model", f), files)
@info "storing modelfiles in $path"

function fix_modelfile(mf)
	# get parameters
	savepath = dirname(mf)
	parameters = parse_savename(mf)[2]
	parameters = NamedTuple{Tuple(Symbol.(keys(parameters)))}(values(parameters))

	# get additional fit info
	# zdim is empty since there was a bug that saved it incorrectly in model file
	#infopars = merge(parameters, (score="latent", zdim="",lr=""))
	sfs = filter(x->occursin("$(parameters.init_seed)", x), sfiles)
	if length(sfs) == 0
#		@info "data for $mf not found"
		return ""
	else
		sd = load(sfs[1]) # score data
	end

	# now create the new parameters
	outpars = merge(parameters, (zdim=sd[:parameters].zdim,lr=sd[:parameters].lr))

	# now save the fixed model data and delete the old model
	sn = joinpath(savepath, savename("model", outpars, "bson",digits=5))
	if ((sn != mf) || force) # only do all of this if the old and new modelfiles are different
		# get model data
		model_data = load(mf)
		# also add the additional fields
		model_data["history"] = sd[:history]
		model_data["fit_t"] = sd[:fit_t]
		model_data["parameters"] = Base.structdiff(sd[:parameters], (score=nothing,))

		rm(mf)
		save(sn, model_data)
	end

	return sn
end

newfiles = map(fix_modelfile, mfiles)

# print files that were not modified
@info "these files were not updated, probably missing the corresponding score files:"
for f in mfiles[newfiles .== ""]
	println(f)
end


text:
using SPTK
using FFTW
using Base.Test

import SPTK:
     # Adaptive mel-generalized analysis
     acep!,
     agcep!,
     amcep!,

     # Mel-generalized cepstrum analysis
     mcep!,
     mcep,
     gcep!,
     gcep,
     mgcep!,
     mgcep,
     uels!,
     uels,
     fftcep!,
     fftcep,
     lpc!,
     lpc,

     # LPC, LSP and PARCOR conversions
     lpc2c!,
     lpc2c,
     lpc2lsp!,
     lpc2lsp,
     lpc2par!,
     lpc2par,
     par2lpc!,
     par2lpc,
     lsp2sp!,
     lsp2sp,

     # Mel-generalized cepstrum conversions
     mc2b!,
     mc2b,
     b2mc!,
     b2mc,
     b2c!,
     b2c,
     c2acr!,
     c2acr,
     c2ir!,
     c2ir,
     ic2ir!,
     ic2ir,
     c2ndps!,
     c2ndps,
     ndps2c!,
     ndps2c,
     gc2gc!,
     gc2gc,
     gnorm!,
     gnorm,
     ignorm!,
     ignorm,
     freqt!,
     freqt,
     frqtr!,
     frqtr,
     mgc2mgc!,
     mgc2mgc,
     mgc2sp!,
     mgc2sp,
     mgclsp2sp!,
     mgclsp2sp,

     # F0 estimation
     swipe,
     rapt,

     # Excitation generation
     excite,

     # Waveform generation filters
     poledf,
     poledf_delay,
     lmadf,
     lmadf_delay,
     lspdf,
     lspdf_delay,
     ltcdf,
     ltcdf_delay,
     glsadf,
     glsadf_delay,
     mlsadf,
     mlsadf_delay,
     mglsadf,
     mglsadf_delay,

     # window functions
     blackman,
     hamming,
     hanning,
     bartlett,
     trapezoid,
     rectangular,

     # Utils
     lspcheck

for fname in [
              "common",
              "lib",
              "adaptive",
              "conversions",
              "excite",
              "f0",
              "mfcc",
              "mgcep",
              "synthesis_filters",
              "utils",
              "window",
              "extend"
    ]
    include(string(fname, ".jl"))
end


text:
#From Chopin's algorithm in N. Chopin, "Fast simulation of truncated Gaussian distributions",
#Stat Comput (2011) 21:275-288

#------------------------------------------------------------
#include("constants.jl")

const N_tn = 4001# Index of the right tail
const yl0 = 0.053513975472# y_l of the leftmost rectangle
const ylN = 0.000914116389555# y_l of the rightmost rectangle

#------------------------------------------------------------
# Compute y_l from y_k
function yl(k::Int64)
    if (k == 1)
        return yl0
    elseif (k == N_tn)
        return ylN
    elseif (k <= 1954)
        return YU[k-1]
    else
        return YU[k+1]
    end
end
#------------------------------------------------------------
# Rejection algorithm with lower truncated exponential proposal
function rtexp(rng::Distributions.AbstractRNG, a::Float64, b::Float64)
    stop = false
    TWOASQ = 2 * (a^2)
    expab = _exp_c(-a * (b - a)) - 1
    z = 0.0
    while !stop
        z = _log_c(1 + Distributions.rand(rng) * expab)
        e = -_log_c(Distributions.rand(rng))
        stop = (TWOASQ * e > z^2)
    end
    return a - z / a
end

# Design variables
const XMIN = -2.00443204036# Left bound
const XMAX = 3.48672170399# Right bound
const KMIN = 5# if kb-ka < kmin then use lower rejection algorithm
const INVH = 1631.73284006# = 1/h, h being the minimal interval range
const I0 = 3271# = - floor(X_TN(0)/h)
const ALPHA = 1.837877066409345 # = _log_c(2*__PI)
const XSIZE = size(X_TN, 1)
const SQ2 = 7.071067811865475e-1# = 1/sqrt(2)
const SQPI = 1.772453850905516# = sqrt(__PI)

function Distributions.std(
    d::Distributions.Truncated{Distributions.Normal{T},Distributions.Continuous},
) where {T<:Real}
    d.untruncated.σ
end

function Distributions.var(
    d::Distributions.Truncated{Distributions.Normal{T},Distributions.Continuous},
) where {T<:Real}
    d.untruncated.σ^2
end

"""
Distributions.rand(rng::Distributions.AbstractRNG, d::Distributions.Truncated{Distributions.Normal{T},Distributions.Continuous}) where {T<:Real}

# Description
Pseudorandom numbers from lower and/or upper truncated Gaussian distribution.
The Gaussian has parameters μ (default 0) and σ (default 1) and is truncated on the interval `[lower, upper]`.

# Output
Returns a sample of `d`.
"""

function Distributions.rand(
    rng::Distributions.AbstractRNG,
    d::Distributions.Truncated{Distributions.Normal{T},Distributions.Continuous},
) where {T<:Real}
    stop = false
    d0 = d.untruncated
    μ = Distributions.mean(d0)
    σ = Distributions.std(d0)
    lower = d.lower
    upper = d.upper
    # Scaling
    if (μ != 0 || σ != 1)
        lower = (lower - μ) / σ
        upper = (upper - μ) / σ
    end

    #-----------------------------

    # Check if lower < upper
    if lower >= upper
        error("upper must be greater than lower!")
        # Check if |lower| < |upper|
    elseif (abs(lower) > abs(upper))
        r = -rand(rng, Distributions.TruncatedNormal(μ, σ, -upper, -lower))# Pair (r,p)
    # If lower in the right tail (lower > XMAX), use rejection algorithm with lower truncated exponential proposal
    elseif (lower > XMAX)
        r = rtexp(rng, lower, upper)
        # If lower in the left tail (lower < XMIN), use rejection algorithm with lower Gaussian proposal
    elseif (lower > XMIN)
        stop = false
        while !stop
            r = Distributions.rand(Distributions.Normal(0, 1))
            stop = (r >= lower) && (r <= upper)
        end
        # In other cases (XMIN < lower < XMAX), use Chopin's algorithm
    else
        # Compute ka
        i = max(1, Int64(I0 + (floor(lower * INVH) + 1)))
        ka = NCELL[i]

        # Compute kb
        if (upper >= XMAX)
            kb = N_tn
        else
            i = max(1, Int64(I0 + (floor(upper * INVH) + 1)))
            kb = NCELL[i]
        end

        # If |upper-lower| is small, use rejection algorithm with lower truncated exponential proposal
        if (abs(kb - ka) < KMIN)
            r = rtexp(rng, lower, upper)
            stop = true
        end

        while (!stop)
            # Sample integer between ka and kb
            k = Int64(floor(Distributions.rand(rng) * (kb - ka + 1)) + ka + 1) #was floor
            if (k >= N_tn)
                # Right tail
                lbound = X_TN[XSIZE]
                z = -_log_c(Distributions.rand(rng))
                e = -_log_c(Distributions.rand(rng))
                z = z / lbound
                if ((z^2 <= 2 * e) && (z < upper - lbound))
                    # Accept this proposition, otherwise reject
                    r = lbound + z
                    stop = true
                end
            elseif (k <= ka + 1) || ((k >= kb - 1) && (upper < XMAX))

                # Two leftmost and rightmost regions
                sim = X_TN[k] + (X_TN[k+1] - X_TN[k]) * Distributions.rand(rng)

                if ((sim >= lower) && (sim <= upper))
                    # Accept this proposition, otherwise reject
                    simy = YU[k] * Distributions.rand(rng)
                    if ((simy < yl(k)) || (sim * sim + 2 * _log_c(simy) + ALPHA) < 0)
                        r = sim
                        stop = true
                    end
                end
            else # All the other boxes
                u = Distributions.rand(rng)
                simy = YU[k] * u
                d = X_TN[k+1] - X_TN[k]
                ylk = yl(k)
                if (simy < ylk)# That's what happens most of the time
                    r = X_TN[k] + u * d * YU[k] / ylk
                    stop = true
                else
                    sim = X_TN[k] + d * Distributions.rand(rng)
                    # Otherwise, check you're below the pdf curve
                    if ((sim * sim + 2 * _log_c(simy) + ALPHA) < 0)
                        r = sim
                        stop = true
                    end
                end
            end
        end
    end

    #-----------------------------

    # Scaling
    if (μ != 0 || σ != 1)
        r = r * σ + μ
    end
    # Compute the probability
    #Z = sqpi * sq2 * σ * ( gsl_sf_erf(upper *s q2) - gsl_sf_erf(lower*sq2) )
    #p = _exp_c(-(((r-μ)/σ)^2)/2) / Z

    return r
end


text:
using Documenter: doctest
using LightQuery
using LightQuery: my_flatten, my_all, partial_map, my_setdiff, named_tuple
using Test: @test, @test_throws

doctest(LightQuery)

# TODO: figure out a better way to trigger this
@test_throws StackOverflowError iterate(
    InnerJoin(By(Int[], identity), By(Int[], identity)),
    (nothing, nothing, false, false),
)

A = (true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, true, false)
@test !my_all(A)
@test A == partial_map(&, true, A)
@test A == partial_map(&, true, A, A)
@test my_setdiff(A, (true,)) == (false,)
@test name"k"(named_tuple((a = 1, b = 2, c = 3, d = 4, e = 5, f = 6, g = 7, h = 8, i = 9, j = 10, k = 11, l = 12, m = 13, n = 14, o = 15, p = 16, q = 17, r = 18))) == 11


text:
# Abstracts
export AbstractFormat, JSONFormat, YAMLFormat

# DataFile
export DataFile

# Interface
export obtain, filename, filepath, url
import Base: dirname, download, parse, show

using JSON

abstract type AbstractFormat end
abstract type JSONFormat <: AbstractFormat end
abstract type YAMLFormat <: AbstractFormat end

"""
    DataFile

Physical data files, it stores the constants.
"""
struct DataFile{Format<:AbstractFormat, ST<:AbstractString}
    url::ST
    localdir::ST
    filename::ST
end

DataFile(url::ST, dir::ST, name::ST) where ST =
    DataFile{JSONFormat, ST}(url, abspath(dir), name)

# DataFile Interface
url(file::DataFile) = file.url
dirname(file::DataFile) = file.localdir
filename(file::DataFile) = file.filename
filepath(file::DataFile) = joinpath(file.localdir, file.filename)

"""
    obtain(file)

obtain data file. Do nothing if this file already exists.
"""
function obtain(file::DataFile)
    if !ispath(filepath(file))
        info(filename(file), prefix="[Download]: ")
        download(file)
    end
end

"""
    download(file)

download data file.
"""
function download(file::DataFile)
    mkpath(dirname(file))
    download(url(file), filepath(file))
end

"""
    parse(file) -> Dict

parse data file to julia dict.
"""
function parse(file::DataFile{JSONFormat})
    raw = open(filepath(file), "r") do f
        Compat.read(f, String)
    end

    # workaround for JSON
    # Download again if JSON.parse fails
    # TODO: md5 check
    try
        JSON.parse(raw)
    catch
        download(file)
        parse(file)
    end
end

function show(io::IO, file::DataFile{F, S}) where {F, S}
    print(io, "DataFile{$F, $S}\n")
    print(io, "  url: $(url(file))\n")
    print(io, "  local: $(filepath(file))")
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@testset "1-Particle Basis"  begin

##

using ACE
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d
using JuLIP.Potentials: i2z, numz

##

maxdeg = 8
r0 = 1.0
rcut = 3.0

trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)

Nat = 15
P1 = ACE.BasicPSH1pBasis(Pr; species = :X)
Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, :X)
evaluate(P1, Rs, Zs, z0)
evaluate_d(P1, Rs[1], Zs[1], z0)

##

for species in (:X, :Si, [:C, :O, :H])
   @info("species = $species")
   Nat = 15
   P1 = ACE.BasicPSH1pBasis(Pr; species = species)
   @info("   test de-serialisation")
   println(@test(all(JuLIP.Testing.test_fio(P1))))

   @info("   test evaluation")
   for ntest = 1:10
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      A = evaluate(P1, Rs, Zs, z0)
      A_ = sum( evaluate(P1, R, Z, z0) for (R, Z) in zip(Rs, Zs) )
      print_tf(@test A ≈ A_)
   end
   println()
   # test that the specification is reproduced correctly!
   @info("    Check specification is correct")
   for iz0 in numz(P1)
      z0 = i2z(P1, iz0)
      P1_spec = ACE.get_basis_spec(P1, z0)
      P1_spec_2 = [ ACE.get_basis_spec(P1, z0, i) for i = 1:length(P1, z0) ]
      println(@test P1_spec == P1_spec_2)
   end
   # Check gradients
   @info("    Check gradients")
   for ntest = 1:30
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      R, Z = Rs[1], Zs[1]
      U = rand(JVecF) .- 0.5; U /= norm(U)
      A = evaluate(P1, R, Z, z0)
      dA = dot.(Ref(U), evaluate_d(P1, R, Z, z0))
      errs = []
      for p = 2:10
         h = 0.1^p
         Ah = evaluate(P1, R + h * U, Z, z0)
         dAh = (Ah - A) / h
         # @show norm(dA - dAh, Inf)
         push!(errs, norm(dA - dAh, Inf))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end
   println()
end

##

end


text:
using Mimi

# When evaluated in the __init__() function, the surrounding module
# is Main rather than Mimi.
@defcomp adder begin
    add    = Parameter(index=[time])
    input  = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.output[t] = p.input[t] + p.add[t]
    end
end


text:
"""
The first part of the Boltzmann Wealth Distribution example.

In the first part of the Boltzmann example, the experiment is ran without
a spatial structure. In the second part, a spatial structure is added,
and agents are required to only give money to agents who are on the
same node.

This example can be ran by navigating to the examples/ folder, starting
a julia REPL session and running:

```
julia> include("boltzmann_wealth_distribution.jl")
```

This will instantiate the model and create a `DataFrame` `data` that
contains the result of running the model 10 steps. After running the
model, the results can be visualized in DataVoyager like this:

```
julia> visualize_data(data);
```

...which should result in a pop-up window that displays graphs
depicting the results of the experiment.

"""
using Agents

"""
Defines the agent type.

The agent type must be a subtype of AbstractAgent.

Commonly, an agent type will require a field for location value in the form
`pos::Tuple{T, T}`. In the first part of this example we will not be using a spatial
structure, therefore we will not define a field for position.

"""
mutable struct MyAgent{T<:Integer} <: AbstractAgent
  "The identifier number of the agent."
  id::T
  "The agent's wealth."
  wealth::T
end

"Define the model type."
mutable struct MyModel{T<:AbstractVector} <: AbstractModel
  "An array of agents."
  agents::T
  "A field for the scheduler function."
  scheduler::Function
end

"Function to instantiate the model."
function instantiate_model(; numagents)
  # Create a list of agents, each with position (1,1) and one unit of
  # wealth.
  agents = [MyAgent(i, 1) for i in 1:numagents]

  # Instantiate and return the model.
  model = MyModel(agents, random_activation)
  return model
end

"""
Define the agent step function.

Defines what the agent should do at each step.
"""
function agent_step!(agent::AbstractAgent, model::AbstractModel)
  # If the agent's wealth is zero, then do nothing.
  if agent.wealth == 0
    return
  # Otherwise, choose a random agent, subtract one unit of own wealth
  # and add one unit of wealth to the randomly chosen agent.
  else
    random_agent = model.agents[rand(1:nagents(model))]
    agent.wealth -= 1
    random_agent.wealth += 1
  end
end

# Instantiate the model.
model = instantiate_model(numagents=100)

# Run the model multiple steps and collect data.
# An array of Symbols for the agent fields that are to be collected, in
# this case wealth is the only variable to be collected.
agent_properties = [:wealth]
# Specifies at which steps data should be collected.
steps_to_collect_data = collect(1:10)
# Use the step function to run the model 10 times and collect data at
# each step.
data = step!(agent_step!, model, 10, agent_properties, steps_to_collect_data)


text:
struct Bound{T}
    state_lower::Vector{T}
    state_upper::Vector{T}
    action_lower::Vector{T}
    action_upper::Vector{T}
end

function Bound(num_state::Int=0, num_action::Int=0;
    state_lower=-Inf * ones(num_state),
    state_upper=Inf * ones(num_state),
    action_lower=-Inf * ones(num_action),
    action_upper=Inf * ones(num_action))
    return Bound(state_lower, state_upper, action_lower, action_upper)
end

const Bounds{T} = Vector{Bound{T}}


text:
#md # # Advanced Tutorial
#md # Notebook [download](advanced_tutorial.ipynb) [nbviewer](@__NBVIEWER_ROOT_URL__/advanced_tutorial.ipynb) [source](advanced_tutorial.jl)

# ## Introduction

# In this tutorial, we will fit the [coal mining disaster dataset](coal_mining_data.tsv) with a Gaussian process modulated
# Poisson process.

# ## Prepare the environment

# We start by importing:
# * MGVI for the posterior fit
# * `Distributions.jl` and `FFTW.jl` to define the statistical model
# * `Optim.jl` to pass `Optim.Options` to MGVI and to find Maximum a posteriori fit that we will use for comparison
# * `StatsBase.jl` for histogram construction from the data and also for error bands visualization
# * `Plots.jl` for visualization

using MGVI

using FillArrays
using DelimitedFiles
using LinearAlgebra
using Random
using StatsBase
using Distributions
using Optim
using ForwardDiffPullbacks

using Plots
using Plots.PlotMeasures
#!jl Plots.default(legendfontsize=10, tickfontsize=10, grid=false, dpi=120, size=(500, 300))
#jl Plots.default(legendfontsize=24, tickfontsize=24, grid=false, dpi=100, size=(1300, 700))

using FFTW

import ForwardDiff
#-
Random.seed!(84612);
#jl mkpath(joinpath(@__DIR__, "plots"));

# ## Load data

# The dataset, which is included with this repository, contains intervals in days between
# disasters occuring at British coal mines between March 1851 and March 1962.
# We build a model by splitting the entire time range into intervals of 365 days.

function read_coal_mining_data(filepath, binsize)
    init_year = empty
    data = empty
    open(filepath) do io
        raw = readline(io)
        while ! occursin("init_date", raw)
            raw = readline(io)
        end

        init_year = parse(Float64, split(split(strip(raw[2:end]), "\t")[2], "-")[1])
        data = readdlm(io, '\t', Int, '\n', comments=true)[:]
    end
    dates_fract_years = init_year .+ cumsum(data)/365
    left_edge = dates_fract_years[1]
    num_bins = ((dates_fract_years[end] - left_edge) ÷ binsize)
    right_edge = left_edge + binsize*num_bins
    fit(Histogram, dates_fract_years, left_edge:binsize:right_edge).weights
end

coal_mine_disaster_data = read_coal_mining_data(joinpath(@__DIR__, "coal_mining_data.tsv"), 1);

# ## Global parameters and the grid

# Now we define several model properties:
# * `DATA_DIM` is the shape of the dataset
# * `DATA_XLIM` specifies the time range of the data
# * `GP_GRAIN_FACTOR` determines the numbers of finer bins which a data bin is split into.
#    This is useful when there are several datasets defined on different grids.
# * `GP_PADDING` adds empty paddings to the dataset. We use a Fourier transform to sample from the Gaussian process
#    with a finite correlation length. `GP_PADDING` helps us to ensure that periodic boundary conditions
#    imposed by a Fourier transform won't affect the data region.

DATA_DIM = size(coal_mine_disaster_data, 1);

data = coal_mine_disaster_data;

DATA_XLIM = [1851., 1962.];

GP_GRAIN_FACTOR = 3;
GP_PADDING = 80;
#-
function produce_bins()
    data_binsize = (DATA_XLIM[2] - DATA_XLIM[1])/DATA_DIM
    gp_binsize = data_binsize/GP_GRAIN_FACTOR
    gp_dim = Integer(((DATA_XLIM[2] - DATA_XLIM[1]) + 2*GP_PADDING) ÷ gp_binsize)
    gp_left_bin_offset = gp_right_bin_offset = (gp_dim - DATA_DIM) ÷ 2
    if (2*gp_left_bin_offset + DATA_DIM*GP_GRAIN_FACTOR) % 2 == 1
        gp_left_bin_offset += 1
    end
    gp_left_xlim = DATA_XLIM[1] - gp_left_bin_offset*gp_binsize
    gp_right_xlim = DATA_XLIM[2] + gp_right_bin_offset*gp_binsize
    gp_left_xs = collect(gp_left_xlim + gp_binsize/2:gp_binsize:DATA_XLIM[1])
    gp_right_xs = collect(DATA_XLIM[2] + gp_binsize/2:gp_binsize:gp_right_xlim)
    gp_data_xs = collect(DATA_XLIM[1] + gp_binsize/2:gp_binsize:DATA_XLIM[2])
    gp_xs = [gp_left_xs; gp_data_xs; gp_right_xs]
    data_idxs = collect(gp_left_bin_offset+1:GP_GRAIN_FACTOR:gp_left_bin_offset+DATA_DIM*GP_GRAIN_FACTOR)
    gp_xs, gp_binsize, data_idxs
end;

# Based on the defined model properties, we generate the grid. GP grid is the fine-grained grid
# with offsets added to the data range.
# * `_GP_XS` represent bin centers of such a fine-grained grid
# * `_GP_BINSIZE` is the width of the bin (that is 1/`GP_GRAIN_FACTOR` of data bin size)
# * `_DATA_IDXS` - integer indices of the left edges of the data bins

_GP_XS, _GP_BINSIZE, _DATA_IDXS = produce_bins();
_GP_DIM = length(_GP_XS);
_GP_HARMONIC_DIST = 1/_GP_DIM/_GP_BINSIZE;

# ## Model parameters

# The Gaussian process in this tutorial is modeled in the Fourier space with zero mean
# and two hyperparameters defining properties of its kernel. To sample from this
# Gaussian process, we also need a parameter per bin that will represent the particular
# realization of the GP in the bin.

function assemble_paridx(;kwargs...)
    pos = 0
    res = []
    for (k, v) in kwargs
        new_start, new_stop = v.start+pos, v.stop+pos
        push!(res, (k, (v.start+pos):(v.stop+pos)))
        pos = new_stop
    end
    (;res...)
end;

# MGVI is an iterative procedure, so we will need to introduce an initial guess for the state of the model.
# We create a vector with size equal to the count of all parameters' `starting_point` and a NamedTuple
# `PARDIX` that assigns names to the sub-regions in this vector. In the correct case:
# * `gp_hyper` is two hyperparameters of the Gaussian process stored in the first two cells of the parameter vector
# * `gp_latent` `_GP_DIM` are parameters used to define the particular realization of the gaussian process,
#    stored at indices between `3` to `2 + _GP_DIM`.
#
# Function `assemble_paridx` is responsible for constructing such a NamedTuple from the parameter specification.

PARIDX = assemble_paridx(gp_hyper=1:2, gp_latent=1:_GP_DIM);

starting_point = randn(last(PARIDX).stop);

# ## Model implementation

function map_idx(idx::Real, idx_range::AbstractUnitRange{<:Integer})
    i = idx - minimum(idx_range)
    n = length(eachindex(idx_range))
    n_2 = n >> 1
    ifelse(i <= n_2, i, i - n)
end

function dist_k(idx::CartesianIndex, ax::NTuple{N,<:AbstractUnitRange{<:Integer}}, harmonic_distances::NTuple{N,<:Real}) where N
    mapped_idx = map(map_idx, Tuple(idx), ax)
    norm(map(*, mapped_idx, harmonic_distances))
end

function dist_array(dims::NTuple{N,<:Real}, harmonic_distances::NTuple{N,<:Real}) where N
    cart_idxs = CartesianIndices(map(Base.OneTo, dims))
    dist_k.(cart_idxs, Ref(axes(cart_idxs)), Ref(harmonic_distances))
end;

# A Gaussian process's covariance in the Fourier space is represented with a diagonal matrix. Values
# on the diagonal follow a squared exponential function with parameters depending on priors.
# A kernel that is diagonal and mirrored around the center represents a periodic and translationally invariant function
# in the coordinate space. This property restricts covariance to have a finite correlation length in the coordinate
# space.
#
# The kernel in the Fourier space is defined on the domain of wave numbers `k`. We model the mirror-symmetrical kernel
# by imposing the mirror symmetry on the vector of the wave numbers. (See `map_idx` for the symmetry implementation)

k = dist_array((_GP_DIM,), (_GP_HARMONIC_DIST,));

# MGVI assumes that all priors are distributed as standard normals `N(0, 1)`; thus,
# to modify the shapes of the priors, we explicitly rescale them at the model implementation phase.
#
# We also exponentiate each prior before using it to tune the squared exponential shape. In doing so,
# we ensure only positive values for the kernel's hyperparameters.
#
# Actually, for the sake of numeric stability we model already square root of the covariance.
# This can be traced by missing `sqrt` in the next level, where we sample from the Gaussian process.

function amplitude_spectrum(d::Real, ampl::Real, corrlen::Real)
    ampl * sqrt(2 * π * corrlen) * exp( -π^2 * d^2 * corrlen^2)
end;

function sqrt_kernel(p)
    kernel_A_c, kernel_l_c = p[PARIDX.gp_hyper]
    kernel_A = 2*exp(kernel_A_c*0.9)*GP_GRAIN_FACTOR
    kernel_l = 12*exp(kernel_l_c/15)/(GP_GRAIN_FACTOR^0.3)
    amplitude_spectrum.(k, kernel_A, kernel_l)
end;

# As a Fourier transform we choose the Discrete Hartley Transform, which ensures that Fourier
# coefficients of the real valued function remain real valued.

ht = FFTW.plan_r2r(zeros(_GP_DIM), FFTW.DHT);

# Before we proceed, let's have a brief look at the kernel's shape. Below
# we plot the kernel in the coordinate space `K(r) = K(x2 - x1)` as a function of time in years
# between two points. As we go further along the `x`-axis, the time interval will increase, and
# the covariance will decrease.

function plot_kernel_model(p, width; plot_args=(;))
    xs = collect(1:Int(floor(width/_GP_BINSIZE)))
    plot!(xs .* _GP_BINSIZE, (ht * (sqrt_kernel(p) .^ 2))[xs] .* _GP_HARMONIC_DIST, label=nothing, linewidth=2.5; plot_args...)
end

plot()
plot_kernel_model(starting_point, 20)

# To make it even more visual, we also plot the structure of the covariance matrix as a heatmap.
# We see that the finite correlation length shows up as a band around the diagonal. We also
# see small artifacts in the antidiagonal corners. These come from the assumption that the
# kernel is periodic.

function plot_kernel_matrix(p)
    xkernel = ht * (sqrt_kernel(p) .^ 2) .* _GP_HARMONIC_DIST
    res = reduce(hcat, [circshift(xkernel, i) for i in 0:(_GP_DIM-1)])'
    heatmap!(_GP_XS, _GP_XS, res; yflip=true, xmirror=true, tick_direction=:out, top_margin=20px, right_margin=30px)
end

plot()
plot_kernel_matrix(starting_point)
#jl png(joinpath(@__DIR__, "plots/gp-covariance-matrix.png"))

# After we defined the square root of the kernel function (`sqrt_kernel`),
# we just follow the regular procedure of sampling from the normal distribution.
# Since the covariance matrix in the Fourier space is diagonal, Gaussian variables
# in each bin are independent of each other. Thus, sampling ends up rescaling
# the `gp_latent` part of the prior vector responsible for the Gaussian process state.
#
# After we produced a sample of Gaussian random values following the kernel model,
# we apply a Fourier transform to return back to the coordinate space.

function gp_sample(p)
    flat_gp = sqrt_kernel(p) .* p[PARIDX.gp_latent]
    (ht * flat_gp) .* _GP_HARMONIC_DIST
end;

# Together with the implementation of `gp_sample` we also need
# to define its version of the `Dual`s. This will allow our
# application of the Hartley transform to be differentiatiable.

function gp_sample(dp::Vector{ForwardDiff.Dual{T, V, N}}) where {T,V,N}
    flat_gp_duals = sqrt_kernel(dp) .* dp[PARIDX.gp_latent]
    val_res = ht*ForwardDiff.value.(flat_gp_duals) .* _GP_HARMONIC_DIST
    psize = size(ForwardDiff.partials(flat_gp_duals[1]), 1)
    ps = x -> ForwardDiff.partials.(flat_gp_duals, x)
    val_ps = map((x -> ht*ps(x) .* _GP_HARMONIC_DIST), 1:psize)
    ForwardDiff.Dual{T}.(val_res, val_ps...)
end;

# Gaussian process realization is meant to serve as a Poisson rate of the Poisson
# process. Since the Gaussian process is not restricted to positive values, we
# exponentiate its values to forcefully make the function positive.

function poisson_gp_link(fs)
    exp.(fs)
end;

# Now when we have a function representing the Poisson rate density,
# we have to integrate it over each data bin to define the Poisson rate in these bins.
# Function `agg_lambdas` does precisely that. When `GP_GRAIN_FACTOR = 1`, this function
# just multiplies the value of the Gaussian process in the bin by the `_GP_BINSIZE`.
# When we have more GP bins per data bin (`GP_GRAIN_FACTOR > 1`), then we apply
# rectangular quadrature to integrate over the bin.

function _forward_agg(data, idxs, steps_forward)
    [sum(data[i:i+steps_forward-1]) for i in idxs]
end;

function agg_lambdas(lambdas)
    gps = _forward_agg(lambdas, _DATA_IDXS, GP_GRAIN_FACTOR) .* _GP_BINSIZE
    xs = _GP_XS[_DATA_IDXS .+ (GP_GRAIN_FACTOR ÷ 2)]
    xs, gps
end;

# Finally, we define the model by using the building blocks defined above:
# * `gp_sample` sample from the Gaussian process with defined `sqrt_kernel` covariance
# * `poisson_gp_link` ensures Gaussian process is positive
# * `agg_lambdas` integrates Gaussian process over each data bin to turn it into a Poisson rate for each bin
# * `model` maps parameters into the product of the Poisson distribution's counting events in each bin.

function model(params)
    fs = gp_sample(params)
    fine_lambdas = poisson_gp_link(fs)
    _, lambdas = agg_lambdas(fine_lambdas)
    Product(fwddiff(Poisson).(lambdas))
end;

# ## Visualization utilities

function agg_full_lambdas(lambdas)
    left_idxs = 1:GP_GRAIN_FACTOR:(_DATA_IDXS[1]-GP_GRAIN_FACTOR)
    left_gp = _forward_agg(lambdas, left_idxs, GP_GRAIN_FACTOR) .* _GP_BINSIZE
    left_xs = _GP_XS[left_idxs .+ (GP_GRAIN_FACTOR ÷ 2)]
    right_idxs = (_DATA_IDXS[end]+1):GP_GRAIN_FACTOR:(size(lambdas, 1) - GP_GRAIN_FACTOR)
    right_gp = _forward_agg(lambdas, right_idxs, GP_GRAIN_FACTOR) .* _GP_BINSIZE
    right_xs = _GP_XS[right_idxs .+ (GP_GRAIN_FACTOR ÷ 2)]
    middle_xs, middle_gp = agg_lambdas(lambdas)
    full_xs = [left_xs; middle_xs; right_xs]
    full_gp = [left_gp; middle_gp; right_gp]
    full_xs, full_gp
end;

function _mean(p; full=false)
    agg_func = if (!full) agg_lambdas else agg_full_lambdas end
    xs, gps = agg_func(poisson_gp_link(gp_sample(p)))
    xs, gps
end;

function plot_mean(p, label="mean"; plot_args=(;), full=false)
    plot!(_mean(p; full=full)...; label=label, linewidth=3, plot_args...)
end;

function plot_prior_samples(num_samples; mean_plot_args=(;))
    for _ in 1:num_samples
        p = randn(last(PARIDX).stop)
        plot_mean(p, nothing; plot_args=mean_plot_args)
    end
end;

function plot_kernel_prior_samples(num_samples, width)
    for _ in 1:num_samples
        p = randn(last(PARIDX).stop)
        plot_kernel_model(p, width)
    end
    plot!()
end;

function plot_data(; scatter_args=(;), smooth_args=(;))
    bar!(_GP_XS[_DATA_IDXS .+ (GP_GRAIN_FACTOR ÷ 2)], data, color=:deepskyblue2, la=0, markersize=2., markerstrokewidth=0, alpha=0.4, label="data"; scatter_args...)
    smooth_step = 4
    smooth_xs = _GP_XS[_DATA_IDXS .+ (GP_GRAIN_FACTOR ÷ 2)][(smooth_step+1):(end-smooth_step)]
    smooth_data = [sum(data[i-smooth_step:i+smooth_step])/(2*smooth_step+1) for i in (smooth_step+1):(size(data, 1)-smooth_step)]
    plot!(smooth_xs, smooth_data, color=:deeppink3, linewidth=3, linealpha=1, ls=:dash, label="smooth data"; smooth_args...)
end;

function plot_mgvi_samples(samples)
    for sample in eachcol(samples)
        if any(isnan.(sample))
            print("nan found in samples", "\n")
            continue
        end
        plot!(_mean(Vector(sample))..., linealpha=0.5, linewidth=2, label=nothing)
    end
    plot!()
end;

function plot_kernel_mgvi_samples(samples, width)
    for sample in eachcol(samples)
        if any(isnan.(sample))
            print("nan found in samples", "\n")
            continue
        end
        plot_kernel_model(sample, width; plot_args=(linealpha=0.5, linewidth=2, label=nothing))
    end
    plot!()
end;

function produce_posterior_samples(p, num_residuals)
    batch_size = 10

    if num_residuals <= 2*batch_size
        batch_size = num_residuals ÷ 2
    end

    est_res_sampler = MGVI._create_residual_sampler(model, p;
                                                    residual_sampler=ImplicitResidualSampler,
                                                    jacobian_func=FwdRevADJacobianFunc,
                                                    residual_sampler_options=(;cg_params=(;abstol=1E-2)))
    batches = []
    for _ in 1:(num_residuals ÷ batch_size ÷ 2)
        batch_residual_samples = MGVI.rand(Random.GLOBAL_RNG, est_res_sampler, batch_size)
        push!(batches, p .+ batch_residual_samples)
        push!(batches, p .- batch_residual_samples)
    end
    reduce(hcat, batches)
end

function _extract_quantile(sorted_gp_realizations, p)
    map(s -> quantile(s, p; sorted=true), eachrow(sorted_gp_realizations))
end;

function plot_posterior_bands(p, num_samples; full=false)
    bands = [(0.997, :red), (0.955, :goldenrod1), (0.683, :green)]
    samples = produce_posterior_samples(p, num_samples)
    xs, first_gp = _mean(samples[1:end, 1]; full=full)
    gp_realizations = reduce(hcat, [_mean(Vector(sample); full=full)[2] for sample in eachcol(samples[1:end, 2:end])]; init=first_gp)
    for (i, one_x_sample) in enumerate(eachrow(gp_realizations))
        gp_realizations[i, 1:end] .= sort(Vector(one_x_sample))
    end
    for (band, color) in bands
        quant_l = _extract_quantile(gp_realizations, (1-band)/2)
        quant_u = _extract_quantile(gp_realizations, (1+band)/2)
        plot!(xs, quant_l; fillrange=quant_u, fillcolor=color, linealpha=0, label=band)
    end
    sample_median = _extract_quantile(gp_realizations, 0.5)
    plot!(xs, sample_median; linewidth=2, linecolor=:grey25, label="median")
end;

# ## Visualization and fitting

# We start by plotting the Gaussian process's dynamic range by
# sampling many possible realizations of it unconditionally on the data.
# We expect the set of lines to populate the regions where there are data.

plot()
plot_data(;scatter_args=(;alpha=0.7))
plot_prior_samples(200, mean_plot_args=(;alpha=0.5))
plot!(ylim=[0, 8])
#jl png(joinpath(@__DIR__, "plots/poisson-dynamic-range.png"))

# We also plot prior samples for the kernel in the coordinate space. The plot below
# shows that the kernel is flexible in the amplitude while the correlation length
# is quite strongly predefined:

plot()
plot_kernel_prior_samples(200, 20)
#jl png(joinpath(@__DIR__, "plots/gp-kernel-dynamic-range.png"))

# Now that we see that the Gaussian process is potentially able to fit the data,
# we plot the initial guess (`starting_point`) to see where we should start from.
# This plot shows:
# * data points
# * smoothed data with a moving average of 9 years
# * Poisson rate for each bin
# * MGVI samples around the mean. At the later stages they can be used to estimate MGVI's uncertainty

plot()
plot_data()
plot_mean(starting_point, "starting point"; plot_args=(;color=:darkorange2))
plot_mgvi_samples(produce_posterior_samples(starting_point, 6))
#jl png(joinpath(@__DIR__, "plots/res-starting-point.png"))

# We also want to introduce the `full` plot that shows not only the data region,
# but includes the region with the padding we added with `GP_PADDING`. We will use
# this plot to make sure that periodic boundary conditions don't interfere with
# the data.

plot()
plot_data()
plot_mean(starting_point, "full gp"; full=true, plot_args=(;color=:pink))
plot_mean(starting_point, "starting point"; plot_args=(;color=:darkorange2))

# Below we also plot the kernel and MGVI samples that represent
# the possible variation of the kernel shape around the mean:

plot()
plot_kernel_model(starting_point, 20; plot_args=(;label="kernel model"))
plot_kernel_mgvi_samples(produce_posterior_samples(starting_point, 6), 20)
#jl png(joinpath(@__DIR__, "plots/kernel-starting-point.png"))

# Let's make a first iteration of the MGVI. For purposes of displaying the convergence curve, we limit `Optim.option` to 1 iteration so that
# MGVI will coverge more slowly.

first_iteration = mgvi_kl_optimize_step(Random.GLOBAL_RNG,
                                        model, data,
                                        starting_point;
                                        num_residuals=3,
                                        jacobian_func=FwdRevADJacobianFunc,
                                        residual_sampler=ImplicitResidualSampler,
                                        optim_options=Optim.Options(iterations=1, show_trace=false),
                                        residual_sampler_options=(;cg_params=(;abstol=1E-2,verbose=false)));

# We again plot data and the Poisson rate. We again show the Gaussian process with padding.
# After one iteration the Poisson rate doesn't seem to get much closer to the data.

plot()
plot_data()
plot_mean(first_iteration.result, "first iteration"; plot_args=(;color=:darkorange2))
plot_mgvi_samples(first_iteration.samples)
#-
plot()
plot_data()
plot_mean(first_iteration.result, "full gp"; full=true, plot_args=(;color=:pink))
plot_mean(first_iteration.result, "first iteration"; plot_args=(;color=:darkorange2))

# Kernel and its MGVI samples changed significantly
# in comparison to the `starting_point` even after the first iteration:

plot()
plot_kernel_model(first_iteration.result, 20; plot_args=(;label="kernel model"))
plot_kernel_mgvi_samples(first_iteration.samples, 20)

# In order to visualize convergence we prepare a few functions to compute,
# store and plot the average posterior likelihood of.

function compute_avg_likelihood(model, samples, data)
    tot = 0
    for sample in eachcol(samples)
        tot += -MGVI.posterior_loglike(model, sample, data)
    end
    tot/size(samples, 2)
end;

function show_avg_likelihood(series)
    scatter!(1:size(series, 1), series, label="-loglike")
end;

# Now we do 30 more iterations of the MGVI and store the average likelihood after each step.
# We feed the fitted result of the previous step as an input to the next iteration.

next_iteration = first_iteration;
avg_likelihood_series = [];
push!(avg_likelihood_series, compute_avg_likelihood(model, next_iteration.samples, data));
for i in 1:30
    tmp_iteration = mgvi_kl_optimize_step(Random.GLOBAL_RNG,
                                          model, data,
                                          next_iteration.result;
                                          num_residuals=3,
                                          jacobian_func=FwdRevADJacobianFunc,
                                          residual_sampler=ImplicitResidualSampler,
                                          optim_options=Optim.Options(iterations=1, show_trace=false),
                                          residual_sampler_options=(;cg_params=(;abstol=1E-2,verbose=false)))
    global next_iteration = tmp_iteration
    push!(avg_likelihood_series, compute_avg_likelihood(model, next_iteration.samples, data))
end;

# First, let's have a look at the convergence plots. We see that MGVI converged after 10 iterations
# while being limited to very poor `Optim` performance.

plot(yscale=:log)
show_avg_likelihood(avg_likelihood_series)
#jl png(joinpath(@__DIR__, "plots/convergence.png"))

# Below we plot the result of the fit. Together with the data and Poisson rate, we also plot
# MGVI residuals. These are samples from the Gaussian posterior, sampled with respect to the posterior's
# covariance structure. Thus MGVI residual samples are deviations from the MGVI fit and represent
# how confident we are about the prediction.

plot(ylim=[0,8])
plot_data()
plot_mgvi_samples(next_iteration.samples)
plot_mean(next_iteration.result, "many iterations"; plot_args=(;color=:darkorange2))
#jl png(joinpath(@__DIR__, "plots/res-many-iter.png"))

# To present credibility intervals we also plot credibility bands. We sample 400 residual samples
# from MGVI and then plot quantiles for each data bin. This should give us a feeling of how compatible
# the MGVI fit is with the data.

plot(ylim=[0,8])
plot_posterior_bands(next_iteration.result, 400)
plot_data()
plot_mean(next_iteration.result, "many iterations"; plot_args=(;color=:darkorange2))
#jl png(joinpath(@__DIR__, "plots/res-bands.png"))

# We also make sure boundary conditions do not interfere with the data. Here is the Gaussian process
# plot with the paddings included:

plot()
plot_data()
plot_mean(next_iteration.result; full=true, plot_args=(;color=:pink))
plot_mean(next_iteration.result, "many iterations"; plot_args=(;color=:darkorange2))

# Let's have a look at the kernel again. We expect the variation of samples
# to become narrower:

plot()
plot_kernel_model(next_iteration.result, 20; plot_args=(;label="kernel model"))
plot_kernel_mgvi_samples(next_iteration.samples, 20)
#jl png(joinpath(@__DIR__, "plots/kernel-many-iter.png"))

# ## Maximum a posteriori estimation

# We build a MAP as a cross check of MGVI results. We simply optimize the posterior likelihood by using `Optim`
# without any particular tuning settings:

max_posterior = Optim.optimize(x -> -MGVI.posterior_loglike(model, x, data), starting_point, LBFGS(), Optim.Options(show_trace=false, g_tol=1E-10, iterations=300));

# We observe that the bump in the middle (around 1910) is caught by the MAP while it is less pronounced in the MGVI fit.
# MAP also has finer structure around 1875 and 1835.

plot()
plot_data()
plot_mean(next_iteration.result, "mgvi mean"; plot_args=(;color=:darkorange2))
plot_mean(Optim.minimizer(max_posterior), "map")
#jl png(joinpath(@__DIR__, "plots/map.png"))

# We also can see the difference at the left edge of the data region. While MGVI smoothed the data,
# the MAP predicted a consequent peak:

plot()
plot_data()
plot_mean(Optim.minimizer(max_posterior), "full gp"; full=true, plot_args=(;color=:darkorange2))
plot_mean(next_iteration.result, "mgvi full gp"; full=true)


text:
module ICPTrajectoryTest

using QPWalkingControl
using StaticArrays
using MathOptInterface
using OSQP
using OSQP.MathOptInterfaceOSQP
using Test
using Random
using ForwardDiff
using LinearAlgebra
using PlanarConvexHulls
using QPControl.Trajectories

const MOI = MathOptInterface

using QPWalkingControl: integrate_icp
import QPWalkingControl: ICPTrajectoryGenerator, push_segment!, solve!, initial_icp, cop, find_segment, initialize!
import QPWalkingControl: HRep, SHRep

const visualize = false
@static if visualize
    import Plots
    using Plots: plot, plot!, gui
    Plots.gr()
end

@testset "integrate_icp" begin
    ξ0 = 1.0
    p = BezierCurve(1, 2, 3, 4)
    ω = 3.0
    for t in range(0.0, 1.0, length=10)
        ξ = integrate_icp(ξ0, p, ω, t)
        ξd = ForwardDiff.derivative(t -> integrate_icp(ξ0, p, ω, t), t)
        @test ξd ≈ ω * (ξ - p(t))
    end
end

function optimizer()
    optimizer = OSQP.Optimizer()
    MOI.set(optimizer, OSQPSettings.Verbose(), false)
    MOI.set(optimizer, OSQPSettings.EpsAbs(), 1e-6)
    MOI.set(optimizer, OSQPSettings.EpsRel(), 1e-8)
    MOI.set(optimizer, OSQPSettings.MaxIter(), 20000)
    # MOI.set(optimizer, OSQPSettings.Polish(), true)
    MOI.set(optimizer, OSQPSettings.AdaptiveRhoInterval(), 25) # required for deterministic behavior
    optimizer
end

@testset "straight line walking" begin
    visualize = true

    g = 9.81
    z = 0.95
    ω = sqrt(g / z)
    num_segments = 8
    generator = ICPTrajectoryGenerator{Float64, 6}(optimizer(), num_segments, ω)

    foot_width = 0.15
    foot_length = 0.3
    step_length = 0.5
    step_width = 0.2
    foot_points = SVector(
        SVector(-foot_length / 2, -foot_width / 2),
        SVector(-foot_length / 2,  foot_width / 2),
        SVector( foot_length / 2,  foot_width / 2),
        SVector( foot_length / 2, -foot_width / 2)
    )

    foot_locations = [foot_points .+ Ref(SVector(0.0, -step_width)), foot_points .+ Ref(SVector(0.0, step_width))]
    convex_hulls = [DConvexHull{Float64}() for i = 1 : num_segments]
    Δts = fill(NaN, num_segments)
    support_indices = Set([1, 2])
    next_swing_index = 1
    first_step = true
    rng = MersenneTwister(1)
    for i = 1 : num_segments
        convex_hull = convex_hulls[i]
        support_points = collect(Iterators.flatten(foot_locations[i] for i in support_indices))
        jarvis_march!(convex_hull, support_points)
        if length(support_indices) == length(foot_locations)
            Δts[i] = first_step ? 1.0 + 0.2 * rand(rng) : 0.2 + 0.2 * rand(rng)
            delete!(support_indices, next_swing_index)
        else
            Δts[i] = 0.6 + 0.2 * rand(rng)
            for i in eachindex(foot_locations)
                push!(support_indices, i)
            end
            foot_displacement = if first_step
                first_step = false
                SVector(step_length, 0.0)
            else
                SVector(2 * step_length, 0.0)
            end
            foot_locations[next_swing_index] = foot_locations[next_swing_index] .+ Ref(foot_displacement)
            next_swing_index = next_swing_index == 1 ? 2 : 1
        end
    end
    support_polygon = Piecewise(Constant.(convex_hulls), cumsum([0.0; Δts]))

    for num_active_segments in 1 : num_segments
        tf = sum(Δts[1 : num_active_segments]) - 10 * eps()
        empty!(generator)
        p0 = SVector(0.02, 0.01)
        ξ0 = p0
        ξf = centroid(convex_hulls[num_active_segments])
        initialize!(generator, p0, ξ0, ξf)
        for i = 1 : num_active_segments
            push_segment!(generator, Δts[i], convex_hulls[i], centroid(convex_hulls[i]))
        end

        solve!(generator)

        # Initial/final ICP
        @test first(generator(0.0)) ≈ generator.initial_icp[] atol=1e-4
        @test first(generator(tf)) ≈ generator.final_icp[] atol=1e-4

        # CoP in support polygon
        for t in range(0.0, tf; length=100)
            @test cop(generator, t) in support_polygon(t)
        end

        # Continuity
        let t = 0.0
            for i in 1 : num_active_segments
                for t′ in (max(0.0, t - eps()), t, t + eps())
                    ξ, ξd = generator(t′)
                    @test ξ ≈ initial_icp(generator, i) atol=1e-4
                    i′, t0′ = find_segment(generator, t′)
                    if t′ > t
                        @test i′ == i
                        @test t0′ == t
                    elseif t′ < t
                        @test i′ == i - 1
                    end
                end
                t += generator.Δts[i]
            end

            for t′ in (t, t + eps())
                ξ, ξd = generator(t′)
                @test ξ ≈ generator.final_icp[] atol=1e-4
                @test ξd ≈ SVector(0.0, 0.0) atol=1e-4
            end
        end

        allocs = @allocated solve!(generator)
        @test allocs == 0

        @static if visualize
            ts = range(0.0, sum(Δts), length=200)
            xlim = [-0.5, 2.5]
            ylim = [-0.5, 0.5]
            width = 1600
            height = round(Int, width * diff(ylim)[1] / diff(xlim)[1])
            plt = plot(xlabel = "x", ylabel = "y", xlim = xlim, ylim = ylim, size = (width, height))
            cops = cop.(Ref(generator), ts)
            plot!(plt, getindex.(cops, 1), getindex.(cops, 2), label = "cop")
            icps = first.(generator.(ts))
            plot!(plt, getindex.(icps, 1), getindex.(icps, 2), label = "icp")
            for i in 1 : num_active_segments
                hull = convex_hulls[i]
                points = push!(copy(vertices(hull)), first(vertices(hull)))
                plot!(plt, getindex.(points, 1), getindex.(points, 2), label = "hull $i")
            end
            gui(plt)
            sleep(1)
        end
    end
end

end # module


text:
# Autogenerated wrapper script for SoXResampler_jll for i686-w64-mingw32
export libsoxr, libsoxr_lsr

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SoXResampler")
JLLWrappers.@declare_library_product(libsoxr, "libsoxr.dll")
JLLWrappers.@declare_library_product(libsoxr_lsr, "libsoxr-lsr.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsoxr,
        "bin\\libsoxr.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsoxr_lsr,
        "bin\\libsoxr-lsr.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
let pointtypes = (GDAL.wkbPoint, GDAL.wkbPoint25D, GDAL.wkbPointM,
        GDAL.wkbPointZM),
    multipointtypes = (GDAL.wkbMultiPoint, GDAL.wkbMultiPoint25D,
        GDAL.wkbMultiPointM, GDAL.wkbMultiPointZM),
    linetypes = (GDAL.wkbLineString, GDAL.wkbLineString25D, GDAL.wkbLineStringM,
        GDAL.wkbLineStringZM),
    multilinetypes = (GDAL.wkbMultiLineString, GDAL.wkbMultiLineString25D,
        GDAL.wkbMultiLineStringM, GDAL.wkbMultiLineStringZM),
    polygontypes = (GDAL.wkbPolygon, GDAL.wkbPolygon25D, GDAL.wkbPolygonM,
        GDAL.wkbPolygonZM),
    multipolygontypes = (GDAL.wkbMultiPolygon, GDAL.wkbMultiPolygon25D,
        GDAL.wkbMultiPolygonM, GDAL.wkbMultiPolygonZM),
    collectiontypes = (GDAL.wkbGeometryCollection,
        GDAL.wkbGeometryCollection25D, GDAL.wkbGeometryCollectionM,
        GDAL.wkbGeometryCollectionZM)

    function GeoInterface.geotype(g::AbstractGeometry)
        gtype = getgeomtype(g)
        if gtype in pointtypes
            return :Point
        elseif gtype in multipointtypes
            return :MultiPoint
        elseif gtype in linetypes
            return :LineString
        elseif gtype == GDAL.wkbLinearRing
            return :LinearRing
        elseif gtype in multilinetypes
            return :MultiLineString
        elseif gtype in polygontypes
            return :Polygon
        elseif gtype in multipolygontypes
            return :MultiPolygon
        elseif gtype in collectiontypes
            return :GeometryCollection
        else
            @warn "unknown geometry type" gtype
            return :Unknown
        end
    end

    function GeoInterface.coordinates(g::AbstractGeometry)
        gtype = getgeomtype(g)
        ndim = getcoorddim(g)
        if gtype in pointtypes
            if ndim == 2
                return Float64[getx(g,0), gety(g,0)]
            elseif ndim == 3
                return Float64[getx(g,0), gety(g,0), getz(g,0)]
            else
                @assert ndim == 0
                @warn("Empty Point")
            end
        elseif gtype in multipointtypes
            return Vector{Float64}[
                GeoInterface.coordinates(getgeom(g,i-1)) for i in 1:ngeom(g)
            ]
        elseif gtype in linetypes || gtype == GDAL.wkbLinearRing
            return Vector{Float64}[
                collect(getpoint(g,i-1)[1:ndim]) for i in 1:ngeom(g)
            ]
        elseif gtype in multilinetypes || gtype in polygontypes
            return Vector{Vector{Float64}}[
                GeoInterface.coordinates(getgeom(g,i-1)) for i in 1:ngeom(g)
            ]
        elseif gtype in multipolygontypes
            return Vector{Vector{Vector{Float64}}}[
                GeoInterface.coordinates(getgeom(g,i-1)) for i in 1:ngeom(g)
            ]
        end
    end

end


text:
#
# METHOD TEMPLATES
# basic method template: "This method has priority over the more general T:<Union{....} template"
#set_wfm(fgen::INSTR{:Agilent33500B},func="sinusoid",ch=1) = write(fgen.obj,"SOURce$ch:FUNC $func")
# more general method template, superseded by basic method template
#set_wfm(fgen::T;func="sinusoid",ch=1) where {T<:Union{F332x1,F335x2}} = write(fgen.obj,"SOURce$ch:FUNC $func")

##
# Reset
reset_instr(obj::T) where {T<:KE200x} = write(obj,"*RST")
# Clear
clear_instr(obj::T) where {T<:KE200x} = write(obj,"*CLS;TRACe:CLEar")
# IDN
get_idn(obj::T) where {T<:KE200x} = query(obj,"*IDN?")

# Function
set_fc(obj::T;fc="VOLT:DC") where {T<:KE200x} = write(obj,"SENS:FUNC '$fc'")
# Function and range
set_fc_range(obj::T;fc="VOLT",rn=10) where {T<:KE200x} = write(obj,"SENS:$fc:RANG $rn")
set_fc_range_auto(obj::T;fc="VOLT") where {T<:KE200x} = write(obj,"SENS:$fc:RANG:AUTO ON")
# Function and NPLC
set_fc_nplc(obj::T;fc="VOLT",nplc=1) where {T<:KE200x} = write(obj,"SENS:$fc:NPLC $nplc")
# Init once
init_meas(obj::T) where {T<:KE200x} = write(obj,":INIT")
# Init continuous
init_meas_continuous(obj::T;st="ON") where {T<:KE200x} = write(obj,":INIT:CONT $st")
# Abort
abort_meas(obj::T) where {T<:KE200x} = write(obj,":INIT")
# Read measurement
function read_meas(obj::T) where {T<:KE200x} # this is broken
	raw  = query(obj,":READ?")
	aux  = split(raw,"N")
	y    = parse(Float64,aux[1])
	return y
end


text:
using Base.Test

include("zeckendorf.jl")

"Tests the zeckendorf function"
function test_zeckendorf()
    @test(zeckendorf(0) == [false])
    @test(zeckendorf(0) == [true])
    @test(zeckendorf(2) == [false, false, true])
    @test(zeckendorf(4) == [false, true, false, true])
    @test(zeckendorf(6) == [false, true, false, false, true])
end

zeckendorf_handler(r::Test.Success) = println("zeckendorf: test passed on $(r.expr)")
zeckendorf_handler(r::Test.Failure) = println("zeckendorf: failure on $(r.expr)")
zeckendorf_handler(r::Test.Error) = rethrow(r)

Test.with_handler(test_zeckendorf, zeckendorf_handler)

"Tests the binary_string function"
function test_binary_string()
    @test(binary_string(Bool[]) == "")
    @test(binary_string(Bool[false]) == "0")
    @test(binary_string(Bool[true]) == "1")
    @test(binary_string(Bool[false, true]) == "10")
    @test(binary_string(Bool[false, true, false]) == "010")
end

binary_string_handler(r::Test.Success) = println("binary_string: test passed on $(r.expr)")
binary_string_handler(r::Test.Failure) = println("binary_string: failure on $(r.expr)")
binary_string_handler(r::Test.Error) = rethrow(r)

Test.with_handler(test_binary_string, binary_string_handler)

max_stars_repo_name: mooxmirror/julia-experiments

max_issues_repo_name: mooxmirror/julia-experiments

max_forks_repo_name: mooxmirror/julia-experiments