text:
"""
    VTKPlotter(eq::Equation, scenario::Scenario, grid::Grid,
               filename::String)

Initialize a VTKPlotter for equation `eq` and scenario `scenario`,
defined on `grid`.
Output name is `filename`.
"""
mutable struct VTKPlotter
    eq::Equation
    scenario::Scenario
    grid::Grid
    filename::String
    collection::WriteVTK.CollectionFile
    plot_counter::Int64

    function VTKPlotter(eq::Equation, scenario::Scenario, grid::Grid,
                        filename::String)
        collection = paraview_collection(filename)
        plot_counter = 0

        new(eq, scenario, grid, filename, collection, plot_counter)
    end

end

function evaluate_dof(eq, grid, dofidx)
    eval_data = Array{Float64, 2}(undef, (length(grid.cells), 1))
    offset = 1
    for cell in grid.cells
        eval_data[offset] = evaluate_basis(grid.basis, grid.dofs[:, dofidx, cell.dataidx], [0.5, 0.5])
        offset += 1
    end
    eval_data
end

function evaluate_dof_subcell(eq, grid, dofidx)
    order = grid.basis.order
    eval_data = Array{Float64, 2}(undef, (length(grid.cells)*order*order, 1))
    offset = 0
    for cell in grid.cells
        offset_cell = cell.center - [cell.size[1]/2, cell.size[2]/2]
        lin_idx = LinearIndices((order,order))
        for i in CartesianIndices((order,order))
            cellcenter_global = offset_cell + [i[1], i[2]] .* cell.size./order - cell.size./order ./ 2
            cellcenter_reference = localposition(cell, cellcenter_global)
            eval_data[offset+lin_idx[i]] = evaluate_basis(grid.basis, grid.dofs[:, dofidx, cell.dataidx], cellcenter_reference)
        end
        offset += order*order
    end
    eval_data
end

"""
    plot(plotter::VTKPlotter)

Write output with `plotter` for timestep.
"""
function plot(plotter::VTKPlotter)
    grid = plotter.grid
    eq = plotter.eq
    scenario = plotter.scenario
    order = grid.basis.order
    vtkcells = Array{MeshCell,1}(undef, length(grid.cells)*order*order)
    vtkpoints = Array{Float64, 2}(undef, (2, length(grid.cells)*order*order*4))

    for (i,cell) in enumerate(grid.cells)
        offset = cell.center - [cell.size[1]/2, cell.size[2]/2]
        subgrid = make_mesh(eq, scenario, order, cell.size/order, offset)
        for (j,subcell) in enumerate(subgrid)
            subcell_idx = (i-1)*order*order + j
            start = (subcell_idx-1) * 4 + 1
            vtkpoints[:, start+0] = subcell.center + [ subcell.size[1]/2,  subcell.size[2]/2]
            vtkpoints[:, start+1] = subcell.center + [ subcell.size[1]/2, -subcell.size[2]/2]
            vtkpoints[:, start+2] = subcell.center + [-subcell.size[1]/2, -subcell.size[2]/2]
            vtkpoints[:, start+3] = subcell.center + [-subcell.size[1]/2,  subcell.size[2]/2]
            vtkcells[subcell_idx] = MeshCell(VTKCellTypes.VTK_QUAD, [start, start+1, start+2, start+3])
        end
    end
    currentfile = @sprintf("%s_%d", plotter.filename, plotter.plot_counter)
    plotter.plot_counter += 1
    vtkfile = vtk_grid(currentfile, vtkpoints, vtkcells)

    for var=1:get_ndofs(eq)
        evaluated_data = evaluate_dof_subcell(eq, grid, var)
        vtk_cell_data(vtkfile, evaluated_data, string(get_variable_name(eq, var)))
    end
    collection_add_timestep(plotter.collection, vtkfile, grid.time)
end

"""
    save(plotter::VTKPlotter)

Save final output file for `plotter`.
"""
function save(plotter::VTKPlotter)
    vtk_save(plotter.collection)
end

ext: jl


text:
# dbn_distributions.jl
# 2019-11-12
# David Merrell
#
#

include("dbn_distributions.jl")

######################################
# IMPLEMENTATION OF GRAPH PRIOR
######################################

"""
parent set prior parameterized by (1) reference graph and (2) lambda
P(G | G', lambda) \\propto exp(-lambda * sum( z_ij, ij not in reference graph ))
"""
@gen (static) function ref_parents_prior(V::Int64, ref_parents::Vector{Int64},
                                                lambda::Float64)
    parents = @trace(refparentprior(V, ref_parents, lambda), :parents)
    return parents
end
# A `Map` operation over vertices yields a graph prior distribution.
ref_graph_prior = Gen.Map(ref_parents_prior)

"""
parent set prior parameterized by (1) edge confidences and (2) lambda;
a "smooth" version of the other graph prior.
"""
@gen function conf_parents_prior(idx::Int64, V::Int64,
                                 parent_confs::Dict,
                                 lambda::Float64)
    parents = @trace(confparentprior(idx, V, parent_confs, lambda), :parents)
    return parents
end
# A `Map` operation over vertices yields a graph prior distribution.
conf_graph_prior = Gen.Map(conf_parents_prior)

"""
`generate_Xp` models the linear relationship between a variable and its parents
in the time series data (i.e., this is the marginal likelihood function)
"""
@gen (static) function generate_Xp(ind::Int64, Xminus::Matrix{Float64},
                                   parents::Vector{Int64}, regression_deg::Int64)
    return @trace(cpdmarginal(Xminus, ind, parents, regression_deg), :Xp)
end
generate_Xplus = Gen.Map(generate_Xp)

# This `SingletonVec` type helps us use the
# `Map` combinator more efficiently.
import Base: getindex, length
mutable struct SingletonVec{T}
    item::T
    len::Int64
end
getindex(sa::SingletonVec, idx) = sa.item
length(sa::SingletonVec) = sa.len

"""
Uniform prior distribution for the lambda variables
(inverse temperatures).
"""
@gen (static) function generate_lambda(l_min::Float64, l_max::Float64)
#@gen function generate_lambda(l_min, l_max)
    return @trace(Gen.uniform(l_min, l_max), :lambda)
end
generate_lambda_vec = Gen.Map(generate_lambda)

"""
Model for the generative process, i.e.
P(G, Lambda|X) propto P(X|G) * P(G|Lambda) * P(Lambda)
"""
@gen (static) function vertex_lambda_dbn_model(parent_confs::Vector{Dict},
                                               Xminus::Array{Float64,2},
                                               lambda_min::Float64,
                                               lambda_max::Float64,
                                               regression_deg::Int64)

    V = size(Xminus, 2)

    # We make SingletonVecs out of these arguments
    # only so they can be shoe-horned into the `Map` operations.
    Vvec = SingletonVec(V, V)
    min_vec = SingletonVec(lambda_min, V)
    max_vec = SingletonVec(lambda_max, V)
    Xminus_vec = SingletonVec(Xminus, V)
    regdeg_vec = SingletonVec(regression_deg, V)

    # This is the actual substance of the probabilistic program
    lambda_vec = @trace(generate_lambda_vec(min_vec, max_vec), :lambda_vec)
    parent_sets = @trace(conf_graph_prior(1:V, Vvec, parent_confs, lambda_vec), :parent_sets)
    Xplus = @trace(generate_Xplus(1:V, Xminus_vec, parent_sets, regdeg_vec), :Xplus)
    return Xplus

end

#import Gen: simulate, generate, update, regenerate, get_args, get_retval, get_choices, get_score, get_gen_fn, project, has_argument_grads, accepts_output_grad
#
#"""
#This trace type is designed for
#"""
#mutable struct ParallelTrace <: Trace
#
#end

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
#using Pkg
#Pkg.add("Flux")
using Flux

using Flux.Tracker

# executable math
f(x) = x^2+1

# f'(x) = 2x
df(x) = gradient(f,x,nest=true)[1] # df is a tuple, [1] gets the first coordinate

df(4)

# f''(x) = 2
ddf(x) = gradient(df,x,nest=true)[1]
ddf(0)

h(x) = -cos(x)^cos(x)

# h'(x) = tan(x)cos(x)^(cos(x)+1)(log(cos(x))+1) obviously
dh(x)=gradient(h,x)[1]
dh(pi/4)

f(x,y,z) = x^2 + y^2 + z^2

#grad(f) = (2x,2y,2z)
gradient(f,1,2,3)

# Quick Example to introduce Params(): Linear Regression

# random initial parameters
W = rand(5,10)
b = rand(5)

fhat(x) = W*x + b

function loss(x,y)
    yhat = fhat(x) # our prediction for y
    return sum((y-yhat).^2)
end

x = rand(10)
y = rand(5)

loss(x,y) # big loss with random parameters
# I have 50+ paramters, how do I pass them all at once?

W = param(W)
b = param(b)

grads = gradient(() -> loss(x, y), Params([W, b]))

# gradient descent

alpha = 0.01 # learning rate, step size, etc
gW = grads[W]
gb = grads[b]

Tracker.update!(W,-alpha*gW) # essentially W = W - alpha * gW. It does something else I don't understand
Tracker.update!(b,-alpha*gb);

loss(x,y)

function sigmoid(x)
    return 1/(1+exp(-x))
end

W1 = param(rand(7, 10))
b1 = param(rand(7))
layer1(x) = W1 * x .+ b1

W2 = param(rand(5, 7))
b2 = param(rand(5))
layer2(x) = W2 * x .+ b2

model(x) = layer2(sigmoid.(layer1(x)))

layer1(x) = Dense(10,7,sigmoid)
layer2(x) = Dense(7,5)

model(x) = layer2(layer1(x))

# or equivalently
model2(x) = Chain(layer1,layer2)

# cool thing about Chain is that it supports indexing
model2(x)[1]

# To train a model call somthing like

train!(objective, parameters, data, optimizer, cb = () -> println("still training..."))

# cb stands for callback. Its useful to updating you about training (e.g. what the loss is currently)
# By default, it is called after every batch. Use Flux.throttle() to change this

@epoch 5 train!(...)

max_stars_repo_name: DAspden/MLGUI

max_issues_repo_name: DAspden/MLGUI

max_forks_repo_name: DAspden/MLGUI


text:
"""
    oocpca(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1.0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="", colsumlist::AbstractString="", dim::Number=3, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false)

Out-of-core PCA.

Input Arguments
---------
- `input` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `initW` : The CSV file saving the initial values of eigenvectors.
- `initV` : The CSV file saving the initial values of loadings.
- `logdir` : The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration.
- `perm` : Whether the data matrix is shuffled at random.

Output Arguments
---------
- `V` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `U` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
"""
function oocpca(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1.0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="", colsumlist::AbstractString="", dim::Number=3, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false)
    # Initial Setting
    pca = OOCPCA()
    pseudocount, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar = init(input, pseudocount, dim, rowmeanlist, rowvarlist, colsumlist, initW, initV, logdir, pca, scale)
    # Perform PCA
    out = oocpca(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, logdir, pca, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, perm)
    # Output
    if outdir isa String
        writecsv(joinpath(outdir, "Eigen_vectors.csv"), out[1])
        writecsv(joinpath(outdir, "Eigen_values.csv"), out[2])
        writecsv(joinpath(outdir, "Loadings.csv"), out[3])
        writecsv(joinpath(outdir, "Scores.csv"), out[4])
        writecsv(joinpath(outdir, "ExpVar.csv"), out[5])
        writecsv(joinpath(outdir, "TotalVar.csv"), out[6])
    end
    return out
end

function oocpca(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, logdir, pca, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, perm)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    l = dim + 5
    its = 3
    @assert 0 < dim ≤ l ≤ min(N, M)
    Ω = rand(Float32, M, l)
    Y = rand(Float32, N, l)
    B = zeros(Float32, l, M)
    # If not 0 the calculation is converged
    n = 1
    # Each epoch s
    println("Random Projection : A*Ω")
    progress = Progress(N)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        # Each step n
        while(n <= N)
            next!(progress)
            # Row vector of data matrix
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            # Random Projection
            tmpY = normx'*Ω
            @inbounds for i in 1:size(tmpY)[2]
                Y[n,i] = tmpY[1,i]
            end
            n += 1
        end
        close(stream)
    end
    # Renormalize with LU factorization
    F = lu!(Y)

    for i in 1:its
        println("Renormalization : A^T * L")
        n = 1
        AtL = zeros(Float32, M, l)
        progress = Progress(N)
        open(input) do file
            stream = ZstdDecompressorStream(file)
            read!(stream, tmpN)
            read!(stream, tmpM)
            # Each step n
            while(n <= N)
                next!(progress)
                # Row vector of data matrix
                read!(stream, x)
                normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec)
                if perm
                    normx .= normx[randperm(length(normx))]
                end
                AtL = AtL .+ normx*F.L[n,:]'
                n += 1
            end
            close(stream)
        end

        println("Renormalization : A * A^T * L")
        n = 1
        progress = Progress(N)
        open(input) do file
            stream = ZstdDecompressorStream(file)
            read!(stream, tmpN)
            read!(stream, tmpM)
            # Each step n
            while(n <= N)
                next!(progress)
                # Row vector of data matrix
                read!(stream, x)
                normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec)
                if perm
                    normx .= normx[randperm(length(normx))]
                end
                @inbounds for i in 1:size(AtL)[2]
                    Y[n,i] = normx'*AtL[:,i]
                end
                n += 1
            end
            close(stream)
        end
        if i < its
            # Renormalize with LU factorization
            F = lu!(Y)
        else
            # Renormalize with QR factorization
            F = qr!(Y)
        end
    end

    println("Calculation of small matrix : Q'A")
    Q = Matrix(F.Q)
    n = 1
    progress = Progress(N)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        # Each step n
        while(n <= N)
            next!(progress)
            # Row vector of data matrix
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            B = B .+ Q[n,:]*normx'
            n += 1
        end
        close(stream)
    end
    # SVD with small matrix
    println("SVD with small matrix : svd(Q'A)")
    W, λ, V = svd(B)
    U = Q*W
    # PC scores, Explained Variance
    Scores = zeros(Float32, M, dim)
    for n = 1:dim
        Scores[:, n] .= λ[n] .* V[:, n]
    end
    ExpVar = sum(λ) / TotalVar
    # Return
    return (V[:,1:dim], λ[1:dim], U[:,1:dim], Scores[:,1:dim], ExpVar, TotalVar)
end


text:
using Revise
using BridgeLandmarks
const BL=BridgeLandmarks
using Random
using Distributions
using DataFrames
using DelimitedFiles
using CSV
using StaticArrays
using LinearAlgebra
using JLD2
using FileIO

Random.seed!(9)

################ set directories ##########################################
workdir = @__DIR__
cd(workdir)
include("../../outdirpath.jl")

dir1 = joinpath(outdirpath,"cardiac/ms")
outdir_ms = mkpath(dir1)

dir2 = joinpath(outdirpath,"cardiac/ahs")
outdir_ahs = mkpath(dir2)

################ read data ##########################################
dat = load("data_cardiac.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]
n = dat["n"]

################ settings and mcmc #################################
ups =  [:innov, :rmmala_pos, :parameter]
adaptskip = 100
skip_saveITER = 25
printskip = 25
ITER = 2500#2500#10_00#0

δpos = (d*n)^(-1/6) * [1e-6, 1e-5]
p_ms = Pars_ms(δpos=δpos,  γinit=.1/√n,  ## LET OP
                  aprior=Pareto(1.0, 0.01), η =  n -> 0.0, dt = 0.01,
                adaptskip=adaptskip, skip_saveITER=skip_saveITER, ρlowerbound=0.8, δa=0.04, σobs=0.001)
p_ahs = Pars_ahs(δpos=δpos,  db=[2.5, 2.5],stdev=.25,γinit=.1, aprior=Pareto(1.0, 0.01), η =  n -> 0.0, dt = 0.001,
                                adaptskip=adaptskip, skip_saveITER=skip_saveITER, ρlowerbound=0.9)

template_estimation(xobsT; xinitq= xobsT[1], pars = p_ms, ITER=ITER, updatescheme = ups, printskip=printskip, outdir=outdir_ms)

#template_estimation(xobsT; xinitq= xobsT[1], pars = p_ahs, ITER=ITER, updatescheme = ups, printskip=printskip, outdir=outdir_ahs)

# from 'old experiments'
#p_ms = Pars_ms(δmom=0.01, σobs = 0.01,δpos=1.0e-6, δa=0.01)
# for ahs adjust domain bounds
#p_ahs = Pars_ahs(δmom=0.01,db=[2.5,1.5],stdev=.25)


text:
"""
NQCModels define the potentials and derivatives that
govern the dynamics of the particles.
These can exist as analytic models or as interfaces to other codes.
"""
module NQCModels

using Reexport: @reexport

export potential
export derivative
export nstates
export ndofs

"""
Top-level type for models.

# Implementation
When adding new models, this should not be directly subtyped. Instead, depending on
the intended functionality of the model, one of the child abstract types should be
subtyped.
If an appropriate type is not already available, a new abstract subtype should be created.
"""
abstract type Model end

Base.broadcastable(model::Model) = Ref(model)

"""
    potential(model::Model, R::AbstractMatrix)

Evaluate the potential at position `R` for the given `model`.
"""
function potential(model::Model, R::AbstractMatrix)
    if ndofs(model) == 1
        if size(R, 2) == 1
            return potential(model, R[1])
        else
            return potential(model, @view R[1,:])
        end
    elseif size(R, 2) == 1
        return potential(model, @view R[:,1])
    else
        throw(MethodError(potential, (model, R)))
    end
end

"""
    potential!(model::Model, V, R::AbstractMatrix)

In-place version of `potential`, used only when mutable arrays are preferred.

Currently used only for `LargeDiabaticModels`, see `diabatic/DiabaticModels.jl`.
"""
function potential!(model::Model, V, R::AbstractMatrix)
    if ndofs(model) == 1
        if size(R, 2) == 1
            return potential!(model, V, R[1])
        else
            return potential!(model, V, view(R, 1, :))
        end
    elseif size(R, 2) == 1
        return potential!(model, V, view(R, :, 1))
    else
        throw(MethodError(potential!, (model, V, R)))
    end

end

"""
    derivative!(model::Model, D, R::AbstractMatrix)

Fill `D` with the derivative of the electronic potential as a function of the positions `R`.

This must be implemented for all models.
"""
function derivative!(model::Model, D, R::AbstractMatrix)
    if ndofs(model) == 1
        if size(R, 2) == 1
            D[1] = derivative(model, R[1])
            return D
        else
            derivative!(model, view(D, 1, :), view(R, 1, :))
            return D
        end
    elseif size(R, 2) == 1
        derivative!(model, view(D, :, 1), view(R, :, 1))
        return D
    else
        throw(MethodError(derivative!, (model, D, R)))
    end
end

"""
    derivative(model::Model, R)

Allocating version of `derivative!`, this definition should be suitable for all models.

Implement `zero_derivative` to allocate an appropriate array then implement `derivative!`
to fill the array.
"""
function derivative(model::Model, R)
    D = zero_derivative(model, R)
    derivative!(model, D, R)
    return D
end

"""
    zero_derivative(model::Model, R)

Create an zeroed array of the right size to match the derivative.
"""
function zero_derivative end

"""
    nstates(::Model)

Get the number of electronic states in the model.
"""
nstates(::Model) = error("This should return the total number of electronic states.")

"""
    ndofs(::Model)

Get the number of degrees of freedom for every atom in the model. Usually 1 or 3.
"""
ndofs(::Model) = error("This should return the number of degrees of freedom for each atom.")

dofs(model::Model) = Base.OneTo(ndofs(model))

state_independent_potential(model, r) = 0.0
state_independent_derivative!(model, derivative, r) = fill!(derivative, zero(eltype(r)))
nelectrons(::Model) = error("This should return the total number of electrons.")
fermilevel(model) = 0.0

eachelectron(model::Model) = Base.OneTo(nelectrons(model))
eachstate(model::Model) = Base.OneTo(nstates(model))

mobileatoms(::Model, n::Int) = Base.OneTo(n)

include("adiabatic/AdiabaticModels.jl")
@reexport using .AdiabaticModels

include("friction/FrictionModels.jl")
@reexport using .FrictionModels

include("diabatic/DiabaticModels.jl")
@reexport using .DiabaticModels

include("plot.jl")

end # module


text:
export AdaptiveMass

mutable struct AdaptiveMass{K} <: MassComponent{Nothing,K}
    "Smoothing kernel."
    W::K
    "Scaling for the kernel widths (typically between 1.2 and 1.5)."
    η::Float
    "Tolerance for Newton's method."
    ε::Float
    "Mass of the fluid particles."
    m::Float
    kernel_widths::Vector{Float}
    mass_densities::Vector{Float}
    correction_scales::Vector{Float}

    """
        AdaptiveMass(
            W::Kernel{N},
            mass_distribution::InitialMassDistribution{N};
            η::Number = 1.2,
            ε::Number = 1e-3,
        ) -> (AdaptiveMass, Vector{SVectorF{N}})
    """
    function AdaptiveMass(
        W::K,
        mass_distribution::InitialMassDistribution{N};
        η::Number = 1.2,
        ε::Number = 1e-3,
    ) where {N,K<:Kernel{N}}
        @check η > 0
        @check ε > 0
        data = mass_distribution()
        m = data.m
        mass_densities = data.mass_densities
        positions = data.positions
        @check data.m > 0
        @checkeq(length(mass_densities), length(positions))
        kernel_widths =
            estimate_initial_kernel_width.(mass_densities, m, η, Ref(Dim{N}()))
        correction_scales = fill(1.0, length(data.positions))
        update_particle_volumes!(
            kernel_widths,
            mass_densities,
            correction_scales,
            positions,
            m,
            W,
            η,
            ε,
        )
        new{K}(W, η, ε, m, kernel_widths, mass_densities, correction_scales),
        positions
    end
end

get_kernel(adaptive_mass::AdaptiveMass) = adaptive_mass.W
get_particle_mass(adaptive_mass::AdaptiveMass) = adaptive_mass.m
get_mass_densities(adaptive_mass::AdaptiveMass) = adaptive_mass.mass_densities
get_kernel_widths(adaptive_mass::AdaptiveMass) = adaptive_mass.kernel_widths

updatemasses!(
    adaptive_mass::AdaptiveMass{<:Kernel{N}},
    positions::AbstractVector{SVectorF{N}},
) where {N} = update_particle_volumes!(adaptive_mass, positions)

function compute_m_ddr_W_r_h(
    ::Union{Unsigned,Tuple{Unsigned,Unsigned}},
    ::Number,
    adaptive_mass::AdaptiveMass,
)
    throw(InvalidStateException)
end

function compute_m_ddr_W_r_h(
    (i, j)::Tuple{Unsigned,Unsigned},
    (ddr_W_r_hᵢ, ddr_W_r_hⱼ)::Tuple{Number,Number},
    adaptive_mass::AdaptiveMass,
)
    m = adaptive_mass.m
    Ωᵢ, Ωⱼ =
        adaptive_mass.correction_scales[i], adaptive_mass.correction_scales[j]
    m * ddr_W_r_hᵢ / Ωᵢ, m * ddr_W_r_hⱼ / Ωⱼ
end

"""
    update_particle_volumes!(
        adaptive_mass::AdaptiveMass,
        positions::AbstractVector{SVectorF{N}},
    )

Update the particle volumes of the mass component according to the given
particle positions.
"""
function update_particle_volumes!(
    adaptive_mass::AdaptiveMass{<:Kernel{N}},
    positions::AbstractVector{SVectorF{N}},
) where {N}
    update_particle_volumes!(
        adaptive_mass.kernel_widths,
        adaptive_mass.mass_densities,
        adaptive_mass.correction_scales,
        positions,
        adaptive_mass.m,
        adaptive_mass.W,
        adaptive_mass.η,
        adaptive_mass.ε,
    )
end

function update_particle_volumes!(
    kernel_widths::AbstractVector{Float},
    mass_densities::AbstractVector{Float},
    correction_scales::AbstractVector{Float},
    positions::AbstractVector{SVectorF{N}},
    m::Number,
    W::Kernel{N},
    η::Number,
    ε::Number,
) where {N}
    MAX_ITER::UInt = 10000

    @dbgasserteq(length(mass_densities), length(positions))

    for (i, (rᵢ, hᵢ)) in enumerate(zip(positions, kernel_widths))
        h, ρ, Ω, n_iter =
            solve_h_ρ_Ω_newton(hᵢ, rᵢ, positions, m, W, η, ε, MAX_ITER)
        if n_iter >= MAX_ITER
            println(
                stderr,
                "warning: particle volume computation reached max number of iterations ($(MAX_ITER))",
            )
        end
        kernel_widths[i] = h
        mass_densities[i] = ρ
        correction_scales[i] = Ω
    end
end

function solve_h_ρ_Ω_newton(
    hᵢ::Number,
    rᵢ::SVectorF{N},
    positions::AbstractVector{SVectorF{N}},
    m::Number,
    W::Kernel{N},
    η::Number,
    ε::Number,
    max_iter::Unsigned,
) where {N}
    @dbgassert hᵢ > 0
    h = hᵢ
    ρ = 0.0
    Ω = 0.0
    n_iter = 0
    while n_iter < max_iter
        ρ, Ω = compute_ρ_Ω_for_h(h, rᵢ, positions, m, W, η)
        h_perturb = clamp(h_perturbance(h, ρ, Ω, m, η, N), -0.99, 0.99)
        h *= 1 + h_perturb

        n_iter += 1
        if abs(h_perturb) < ε
            break
        end
    end
    h, ρ, Ω, n_iter
end

h_perturbance(
    h::Number,
    ρ::Number,
    Ω::Number,
    m::Number,
    η::Number,
    N::Integer,
) = (m * (η / h)^N - ρ) / (N * ρ * Ω)

function compute_ρ_Ω_for_h(
    h::Number,
    rᵢ::SVectorF{N},
    positions::AbstractVector{SVectorF{N}},
    m::Number,
    W::Kernel{N},
    η::Number,
) where {N}
    ρ = 0.0
    Ω = 0.0
    for rⱼ in positions
        q = norm(rⱼ - rᵢ) / h
        if nonzeroat(W, q)
            ρ += W(q, h)
            Ω += ddh(W, q, h)
        end
    end
    ρ *= m
    Ω = 1 + m * Ω * h / (N * ρ)
    ρ, Ω
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
include("../julia/TBIO.jl")

using .TBIO
using Formatting
using JSON

arrs = [ [1.0, 2.0, 3.0], [4.0, 5.0] ]

println(TBIO.calChecksum(arrs))

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Base.Test
using XSV

cases_path = splitdir(Base.source_path())[1]

@test readxsv("1,2,3\n4,5,6")[1][1]=="1"
@test readxsv("1,2,3\n4,5,6")[1][3]=="3"
@test readxsv("1,2,3\n4,5,6")[2][2]=="5"

@test readxsv("1;2;3\n4;5;6", delimiter=';')[1][3]=="3"
@test readxsv("1;2;3\n4;5;6", delimiter=';')[2][2]=="5"

@test readxsv("1;2;\"\"\"3\"\n4;5;6", delimiter=';')[1][3]=="\"3"
@test readxsv("1;2;\"\"\"3\"\n4;5;6", delimiter=';')[2][2]=="5"

@test readxsv("\"1\";2;\"3\"\n4;5;6",     delimiter=';')[1][3]=="3"
@test readxsv("1;2;3\n\"4\";\"5\";\"6\"", delimiter=';')[2][2]=="5"

@test readxsv("B1B;2;B3B\n4;5;6",   delimiter=';', quotechar='B')[1][3]=="3"
@test readxsv("1;2;3\nB4B;B5B;B6B", delimiter=';', quotechar='B')[2][2]=="5"

data = freadxsv(joinpath(cases_path, "cases/1.csv"))

@test data[2][2]=="dog"
@test data[3][2]=="cat"
@test data[3][3]=="11"
@test data[4][2]=="penguin"
@test data[4][3]=="12"
@test length(data)==4

xsv_stream = fiterxsv(joinpath(cases_path, "cases/1.csv"))
@test consume(xsv_stream)[1]=="id"
@test consume(xsv_stream)[2]=="dog"
@test consume(xsv_stream)[3]=="11"
@test consume(xsv_stream)[2]=="penguin"

xsv_stream = fiterxsvh(joinpath(cases_path, "cases/1.csv"))
@test consume(xsv_stream)["type"]=="dog"
@test consume(xsv_stream)["count"]=="11"
@test consume(xsv_stream)["type"]=="penguin"

xsv_stream = fiterxsvht(joinpath(cases_path, "cases/1.csv"))
@test consume(xsv_stream).type_0=="dog"
@test consume(xsv_stream).count=="11"
@test consume(xsv_stream).type_0=="penguin"


text:
using Powseeker
using Plots
using POMDPs
using POMDPToolbox
using ForwardDiff

mdp = PowseekerMDP()
rng = MersenneTwister(100)

tries = [
            SkierState(1, [-550, -1150], deg2rad(60)),
            SkierState(1, [-50, 1400], deg2rad(270)),
            SkierState(1, [1500, 0], deg2rad(180)),
            SkierState(1, [-100, -100], deg2rad(225)),
            SkierState(1, [-100, -100], deg2rad(270)),
            SkierState(1, [-1600, 200], deg2rad(300)),
            SkierState(1, [-1600, 200], deg2rad(10)),
#             SkierState(1, [0, 1550], deg2rad(225)),
#             SkierState(1, [-350, -750], deg2rad(150)),
#             SkierState(1, [-350, -750], deg2rad(315)),
#             SkierState(1, [200, -1700], 0)
        ]

results = []
for s in tries
    n, r = generate_sr(mdp, s, 0.0, rng)
    @show r
    push!(results, n)
    @show mdp.topology(n.pos) - mdp.topology(s.pos)
end

plot(mdp)
for i in 1:length(tries)
    plot!([tries[i], results[i]])
end
gui()


text:
using FileIO
using JLD2
using Optim, LineSearches
using LinearAlgebra: I, norm
using TimerOutputs
using OMEinsum: get_size_dict, optimize_greedy,  MinSpaceDiff
using Zygote

"""
    logoverlap(Au, Ad, M)
````
    ┌── Au──┐                                 a ────┬──── c
    │   │   │         ┌──  Ad ──┐             │     b     │
    FL─ M ──FR  - 1/2 FL_n │   FR_n           ├─ d ─┼─ e ─┤
    │   │   │         └──  Ad ──┘             │     g     │
    └── Ad──┘                                 f ────┴──── h
````
"""
function logoverlap(Au, Ad, M)
    # Ad /= norm(Ad)
    _, FLud = leftenv(Au, conj(Ad), M)
    _, FRud = rightenv(Au, conj(Ad), M)
    _, FLd_n = norm_FL(Ad, conj(Ad))
    _, FRd_n = norm_FR(Ad, conj(Ad))

    nd = ein"(ad,acb),(dce,be) ->"(FLd_n,Ad,conj(Ad),FRd_n)[]/ein"ab,ab ->"(FLd_n,FRd_n)[]
    Ad /= sqrt(nd)
    AuM = ein"(((adf,abc),dgeb),fgh),ceh -> "(FLud,Au,M,conj(Ad),FRud)[]/ein"abc,abc -> "(FLud,FRud)[]
    -log(abs2(AuM))
end

function init_mps(;infolder = "./data/", D::Int = 2, χ::Int = 5, tol::Real = 1e-10, verbose::Bool = true)
    in_chkp_file = infolder*"/D$(D)_chi$(χ)_tol$(tol).jld2"
    if isfile(in_chkp_file)
        mps = load(in_chkp_file)["mps"]
        verbose && println("load mps from $in_chkp_file")
    else
        mps = rand(ComplexF64,χ,D,χ)
        verbose && println("random initial mps $in_chkp_file")
    end
    _, FL_n = norm_FL(mps, conj(mps))
    _, FR_n = norm_FR(mps, conj(mps))
    n = ein"(ad,acb),(dce,be) ->"(FL_n,mps,conj(mps),FR_n)[]/ein"ab,ab ->"(FL_n,FR_n)[]
    mps /= sqrt(n)
    return mps
end

function onestep(M::AbstractArray; atype = Array, infolder = "./data/", outfolder = "./data/", χ::Int = 5, tol::Real = 1e-10, f_tol::Real = 1e-6, opiter::Int = 100, optimmethod = LBFGS(m = 20), verbose= true, savefile = true)
    D = size(M, 1)
    mps = init_mps(infolder = infolder, D = D, χ = χ, tol = tol, verbose = verbose)

    Au, Ad = mps, mps
    to = TimerOutput()
    f(Ad) = @timeit to "forward" logoverlap(Au, Ad, M)
    ff(Ad) = logoverlap(Au, Ad, M)
    g(Ad) = @timeit to "backward" Zygote.gradient(ff,atype(Ad))[1]
    if verbose
        message = "time  iter   loss           grad_norm\n"
        printstyled(message; bold=true, color=:blue)
        flush(stdout)
    end
    res = optimize(f, g,
        Ad, optimmethod,inplace = false,
        Optim.Options(f_tol=f_tol, iterations=opiter,
        extended_trace=true,
        callback=os->writelog(os, outfolder, D, χ, tol, savefile, verbose)),
        )
    Ad = Optim.minimizer(res)
    _, FLud_n = norm_FL(Au, conj(Ad))
    _, FRud_n = norm_FR(Au, conj(Ad))
    fidelity = norm(ein"(ad,acb),(dce,be) ->"(FLud_n,Au,conj(Ad),FRud_n)[]/ein"ab,ab ->"(FLud_n,FRud_n)[])
    if verbose
        message = "fidelity   = $(fidelity) \ndifference = $(norm(Au-Ad)) \n"
        printstyled(message; bold=true, color=:green)
        flush(stdout)
    end
    return Ad
end

"""
    writelog(os::OptimizationState, key=nothing)

return the optimise infomation of each step, including `time` `iteration` `energy` and `g_norm`, saved in `/data/model_D_chi_tol_maxiter.log`. Save the final `ipeps` in file `/data/model_D_chi_tol_maxiter.jid2`
"""
function writelog(os::OptimizationState, outfolder, D, χ, tol, savefile, verbose)
    message = "$(round(os.metadata["time"],digits=1))    $(os.iteration)    $(round(os.value,digits=8))    $(round(os.g_norm,digits=8))\n"

    if verbose
        printstyled(message; bold=true, color=:blue)
        flush(stdout)
    end

    if savefile
        !(isdir(outfolder)) && mkpath(outfolder)
        logfile = open(outfolder*"/D$(D)_chi$(χ)_tol$(tol).log", "a")
        write(logfile, message)
        close(logfile)
        save(outfolder*"/D$(D)_chi$(χ)_tol$(tol).jld2", "mps", os.metadata["x"])
    end
    return false
end

function optimisemps(M::AbstractArray; atype = Array, infolder = "./data/", outfolder = "./data/", χ::Int = 5, tol::Real = 1e-10, f_tol::Real = 1e-6, opiter::Int = 100, optimmethod = LBFGS(m = 20), verbose= true,  mapsteps = 10, updown = true, downfromup = false)
    Au = nothing
    Ad = nothing
    Md = permutedims(M,(1,4,3,2))
    for i in 1:mapsteps
        if verbose
            message = "iter: $i   \n"
            printstyled(message; bold=true, color=:red)
            flush(stdout)
        end
        direction = "↑"
        Au = onestep(M; atype = atype, infolder = infolder*"$(direction)/", outfolder = outfolder*"$(direction)/", χ = χ, tol = tol, f_tol = f_tol, opiter = opiter, optimmethod = optimmethod, verbose = verbose, savefile = true)
        if updown
            !downfromup && (direction = "↓")
            Ad = onestep(Md; atype = atype, infolder = infolder*"$(direction)/", outfolder = outfolder*"$(direction)/", χ = χ, tol = tol, f_tol = f_tol, opiter = opiter, optimmethod = optimmethod, verbose = verbose, savefile = true)
        else
            Ad = Au
        end
    end
    return Au, Ad
end


text:
"""
Calculate the plasma dispersion function

```
principal(n, x) = x.^n .* exp(-x.^2) / √π
integrande(z, n, x) = P(x, n) / (x - z)
residue(z, n, x) = (imag(z) < 0 ? 2 : imag(z) == 0 ? 1 : 0) * im * π * principal(n, z)
Z(z, n) = ∫integrand(z, n, x)dx + residue(z, n, x)
```

# Example
```jldoctest
  plasma_dispersion_function(0.0) ≈ im * √π
```
"""
module PlasmaDispersionFunctions

using SpecialFunctions

export plasma_dispersion_function

"""@docs
    plasma_dispersion_function(z::T,power::Unsigned=UInt64(0),Z₋₁=missing)

Return the value of the plasma dispersion function
This implementation includes the residue, which is easy to verify
because Z(0) = im sqrt(π).

...
# Arguments
- `z::T`: the argument to the plasma dispersion function
- `power::Unsigned=UInt64(0)`: the moment of the integral
- `Z₋₁=missing`: (optional) the value of the plasma dispersion function for
power-1. If already calculated, passing this in saves some time.
...

# Example
```julia
  plasma_dispersion_function(1.0)
```

# References
[1] D. Fried and S. D. Conte, "The Plasma Dispersion Relation", Elsevier, 1961
[2] S.D. Baalrud, Phys. Plasmas 20, 012118 (2013) and put ν = -Inf
[3] M. Sampoorna et al., Generalized Voigt functions and their derivatives,
  Journal of Quantitative Spectroscopy & Radiative Transfer (2006),
  doi:10.1016/j.jqsrt.2006.08.011
"""
function plasma_dispersion_function(z::T, power::Unsigned=UInt64(0), Z₋₁=missing
    ) where {T<:Number}
  function _const(i::Signed)
    i == 1 && return T(1) # make sure that 1 is handled quickly
    return iseven(i) ? T(0) : T(prod(1:2:(i - 2)) * T(2)^div(1 - i, 2))
  end
  if ismissing(Z₋₁)
    Z̄ = sqrt(π) * erfcx(Complex(imag(z), -real(z))) # hot loop; -im * z
    Z = Complex(-imag(Z̄), real(Z̄)) # hot loop; im * Z̄
    for i ∈ 1:power
      Z = z * Z + _const(Signed(i))
    end
    return Z
  else
    return z * Z₋₁ + _const(Signed(power))
  end
end

function plasma_dispersion_function(z::T, n::Integer, Z₋₁=missing
    ) where {T}
  n < 0 && throw(ArgumentError("n, $n, must be >= 0"))
  return plasma_dispersion_function(z::T, Unsigned(n), Z₋₁)
end

end


text:
module OptimDash

include("render_funcs.jl")
using Dash, DashBootstrapComponents
using DashHtmlComponents, DashCoreComponents
using PlotlyJS
using Optim

function eval_definition(strs...)
    return eval.(Meta.parse.(strs))
end

#Initialize channel to get data during optimization
chn = Base.Channel{Vector{Float64}}(Inf)
x_graph = []
y_graph = []
y_graph1 = []
itr = 0
eval_f = nothing
range_value = []

function make_app()
    global chn, x_graph, y_graph, itr, eval_f
    global range_value

    chn = Base.Channel{Vector{Float64}}(Inf)
    x_graph = []
	y_graph = []
    y_graph1 = []
    range_value = []
    itr = 0
    eval_f = nothing

    external_stylesheets = "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/css/font-awesome.min.css"
    app = dash(external_stylesheets=[dbc_themes.BOOTSTRAP,external_stylesheets],
	    suppress_callback_exceptions=true)

    app.layout = html_div([
        render_all()
    ])
    callback!(app,
        Output("dopt-params-plot-loc", "children"),
        Input("dopt-btn-optim", "n_clicks"),
        State("dopt-user-funcs", "value"),
        State("dopt-user-func-init", "value"),
    ) do n_clicks, user_funcs, funcs_init
        ctx = callback_context()
        if isempty(ctx.triggered) || (user_funcs isa Nothing) || (funcs_init isa Nothing)
            throw(PreventUpdate())
        end
        global chn, eval_f, x_graph, y_graph, y_graph1, itr
        global range_value

        itr = 0
        eval_f, init_f = eval_definition(
            user_funcs, funcs_init,
        )
        x_graph = [0]
        y_graph = [init_f[1]]
        y_graph1 = [init_f[2]]

        if length(init_f) != 2
            throw(PreventUpdate())
        end

        global prevx = copy(init_f)
        function g(f)
            function _g(x)
                global prevx
                if !isapprox(x, prevx, atol = 1e-4)
                    put!(chn, x)
                    prevx = copy(x)
                end
                return Base.invokelatest(f, x)
            end
            return _g
        end
        alg = BFGS()
        @async optimize(g(eval_f), init_f, alg, Optim.Options(x_abstol = 1e-3, f_abstol = 1e-3, g_abstol = 1e-3, f_calls_limit = 100))
        return render_liveupdates()
    end

    callback!(app,
        Output("dopt-params-state", "figure"),
        Output("interval-component", "disabled"),
        Input("interval-component", "n_intervals"),
        State("dopt-params-state", "figure"),
        prevent_initial_call=true) do n, state_fig
        st = false
        global x_graph, y_graph, y_graph1, chn, itr, eval_f
        global range_value
        if isready(chn)
            val = take!(chn)
            itr += 1
            append!(x_graph, itr)
            append!(y_graph, val[1])
			append!(y_graph1, val[2])
        else
            st = true
        end
        return Dict(
                "data" => [
                    Dict(
                        "x" => x_graph,
                        "y" => y_graph,
                        "mode" => "line",
                    ),
                    Dict(
                        "x" => x_graph,
                        "y" => y_graph1,
                        "mode" => "line",
                    )
                ]
            ), st
    end
    return app
end

end # module


text:
using InteractiveUtils: code_llvm, code_warntype
using Test

"""
    llvm_ir(f, args) :: String

Get LLVM IR of `f(args...)` as a string.
"""
llvm_ir(f, args) = sprint(code_llvm, f, Base.typesof(args...))

nmatches(r, s) = count(_ -> true, eachmatch(r, s))

function vdot(xs, ys)
    d = zero(eltype(xs)) + zero(eltype(ys))
    @simd for i in eachindex(xs, ys)
        x = @inbounds xs[i]
        y = @inbounds ys[i]
        d += x * y
    end
    return d
end

@testset begin
    v = Float64[]
    @test nmatches(r"fadd fast <[0-9]+ x double>", llvm_ir(vdot, (v, v))) >= 4
end


text:
function integrate(y0, tspan, ps::AbstractCR3BPIndirectParams, dynamicsFlag::CR3BP,
    integrationTypeFlag::InitializationFlag, homotopyFlag::HomotopyFlag;
    copyParams = false, termCallbacks = false, inPlace = false)

   # Instantiate problem
   prob = createCR3BPODEProb(y0, tspan, ps, integrationTypeFlag, homotopyFlag;
        copyParams = copyParams, termCallbacks = termCallbacks, inPlace = inPlace)

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        save_everystep = false,
        save_start = false,
        initialize_save = false,
        maxiters = 1e6
        )

    # Compute time to desired final time if terminated
    timeToFinalTime = 0.0
    if termCallbacks && sol.retcode != :Success
        timeToFinalTime = tspan[2] - sol.t[end]
    end

    # Return final states and co-states
    return sol[end], timeToFinalTime
end

function integrate(y0, tspan, ps::AbstractCR3BPIndirectParams, dynamicsFlag::CR3BP,
    integrationTypeFlag::Solving, homotopyFlag::HomotopyFlag;
    copyParams = false, termCallbacks = false, inPlace = false)

    # Instantiate problem
    prob = createCR3BPODEProb(y0, tspan, ps, integrationTypeFlag, homotopyFlag;
        copyParams = copyParams, termCallbacks = termCallbacks, inPlace = inPlace)

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        save_everystep = false,
        save_start = false,
        initialize_save = false,
        maxiters = 1e6
        )

    # Return final states and co-states
    return sol[end]
end

function integrate(y0, tspan, ps::AbstractCR3BPIndirectParams, dynamicsFlag::CR3BP,
    integrationTypeFlag::FullSolutionHistory, homotopyFlag::HomotopyFlag;
    copyParams = false, termCallbacks = false, inPlace = false)

    # Instantiate problem
    prob = createCR3BPODEProb(y0, tspan, ps, integrationTypeFlag, homotopyFlag;
        copyParams = copyParams,
    termCallbacks = termCallbacks, inPlace = inPlace, save_positions = (true, true))

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        maxiters = 1e6
        )

    # Return final states and co-states
    return sol
end

function integrateWithHomotopy(y0, tspan, ϵ, ps::AbstractCR3BPIndirectParams, dynamicsFlag::CR3BP,
    integrationTypeFlag::Solving, homotopyFlag::HomotopyFlag;
    copyParams = false, termCallbacks = false, inPlace = false)

    # Instantiate problem
    prob = createCR3BPODEProb(y0, tspan, ps, integrationTypeFlag, homotopyFlag; copyParams = copyParams,
        termCallbacks = termCallbacks, inPlace = inPlace, ϵ = ϵ)

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        save_everystep = false,
        save_start = false,
        initialize_save = false,
        maxiters = 1e6
        )

    # Return final states and co-states
    return sol[end]
end

function integrate(z0, tspan, ps::CR3BPIndirectWithSTMParams, dynamicsFlag::CR3BP,
    integrationTypeFlag::SolvingWithSTM, homotopyFlag::HomotopyFlag;
    copyParams = false)

   # Instantiate problem
   prob = createCR3BPODEWithSTMProb(z0, tspan, ps, homotopyFlag; copyParams)

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        save_everystep = false,
        save_start = false,
        initialize_save = false,
        maxiters = 1e10
        )

    # Return final states and co-states
    return sol[end]
end

function integrateWithHomotopy(z0, tspan, ϵ, ps::CR3BPIndirectWithSTMParams, dynamicsFlag::CR3BP,
    integrationTypeFlag::SolvingWithSTM, homotopyFlag::HomotopyFlag; copyParams = false)

    # Instantiate problem
    prob = createCR3BPODEWithSTMProb(z0, tspan, ps, homotopyFlag; copyParams = copyParams, ϵ = ϵ)

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        save_everystep = false,
        save_start = false,
        initialize_save = false,
        maxiters = 1e10
        )

    # Return final states and co-states
    return sol[end]
end

function integrate(x0, tspan, scenario::String, dynamicsFlag::CR3BP,
    integrationTypeFlag::FullSolutionHistoryNoControl; inPlace = false)

   # Instantiate problem
   prob = createCR3BPODENoControlProb(x0, tspan, scenario; inPlace = inPlace)

    # Solve ode
    sol = solve(
        prob,
        Vern9(),
        reltol = 1e-14,
        abstol = 1e-14,
        maxiters = 1e6
        )

    # Return final states and co-states
    return sol
end


text:
module RecursiveClusterAlgorithm

# not currently in project but can compile locally!
using VortexDistributions

using VortexDistributions:VortexGroup, Dipole, xpos, ypos

using LinearAlgebra, SparseArrays, LightGraphs, SimpleWeightedGraphs

mutable struct Cluster <: VortexGroup
    vortices::Array{PointVortex,1}
    tree::Array{LightGraphs.SimpleGraphs.SimpleEdge{Int64},1}
end

include("get_dipoles.jl")
include("seed_clusters.jl")
include("utils.jl")
Cluster(vort::Array{PointVortex,1}) = Cluster(vort,spanning_tree(xpos(vort),ypos(vort)))

export Cluster, get_dipoles, seed_clusters, distances, periodic_distances,
sparse_distances, spanning_tree, vortex_marker, plot_vortices!,
plot_cluster!, vortex_marker

end # module


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

name = "GibbsSeaWater"
version = v"3.5.1"

# Collection of sources required to build GibbsSeaWater
# note 3.05.0-4 is too old to be cross-compiled on FreeBSD
sources = [
    ArchiveSource(
        "https://github.com/TEOS-10/GSW-C/archive/d392e91cb63341f543ed1609c4eff613055ab3cb.zip", "454c46ec00468aeba048fa4f3cee095f927925ab699b1d1b3a25a838cef2d22c"),
]

# Bash recipe for building across all platforms
script = raw"""
if [[ "${target}" == aarch64-apple-* ]]; then
    # Linking libomp requires the function `__divdc3`, which is implemented in
    # `libclang_rt.osx.a` from LLVM compiler-rt.
    LDFLAGS="-L${libdir}/darwin -lclang_rt.osx"
fi
cd $WORKSPACE/srcdir/GSW-C-*
cc -fPIC -c -O3 -Wall gsw_oceanographic_toolbox.c gsw_saar.c
cc $LDFLAGS -fPIC -shared -o libgswteos.$dlext gsw_oceanographic_toolbox.o gsw_saar.o -lm
mkdir -p ${libdir}
cp libgswteos.$dlext ${libdir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libgswteos", :libgswteos)
]

llvm_version = v"13.0.1"
# Dependencies that must be installed before this package can be built
dependencies = [
    # We need libclang_rt.osx.a, because this library provides the
    # implementation of `__divdc3`.
    BuildDependency(PackageSpec(name="LLVMCompilerRT_jll", uuid="4e17d02c-6bf5-513e-be62-445f41c75a11", version=llvm_version); platforms=[Platform("aarch64", "macos")]),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_llvm_version=llvm_version)

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
function variable(m::Model,k;opt...)
    push!(m[:Vs][k],num_variables(m)+1)
    variable(m;opt...)
end

function Model(optimizer,K;opt...)
    m = Model(optimizer;opt...)
    m[:Vs]= Dict(k=>Int[] for k in K)
    return m
end


text:
function get_rating(rating)
    input = (split.(readlines("day3.txt"),""))
    input_mat = hcat(input...) .|> x -> parse(Int64,x)

    for i in 1:(size(input_mat)[1])
        dd = [[j,count(x->x==j,input_mat[i,:])] for j in u]
        ddd = hcat(dd...)
        if rating == "oxygen"
            if all(y->y==ddd[2,1],ddd[2,:])
                inds = 1
            else inds = ddd[1,argmax(ddd[2,:])]
            end
        elseif rating == "CO2"
            if all(y->y==ddd[2,1],ddd[2,:])
                inds = 0
            else inds = ddd[1,argmin(ddd[2,:])]
            end
        end
        input_mat = input_mat[:,input_mat[i,:] .== inds]
        if size(input_mat)[2] == 1
            break
        end
    end
    return "0b"*join(input_mat) |> x -> parse(Int64,x)
end

get_rating("oxygen")*get_rating("CO2")


text:
# to run tests, type in Julia: Pkg.test("StateSpaceModelsEM")

include("../src/utils/utils_test.jl")
include("../src/ACQR/ACQR_test.jl")
include("../src/ACQRS/ACQRS_test.jl")
include("../src/ABCDQR/ABCDQR_test.jl")
include("../src/acqr/acqr_test.jl")


text:
import Base: length, map, show, copy, cat, start, done, next

export List, Nil, Cons, cons, nil, head, tail, list

abstract List{T}

type Nil{T} <: List{T}
end

type Cons{T} <: List{T}
    head::T
    tail::List{T}
end

cons{T}(h, t::List{T}) = Cons{T}(h, t)

nil(T) = Nil{T}()
nil() = nil(Any)

head(x::Cons) = x.head
tail(x::Cons) = x.tail

function show{T}(io::IO, l::List{T})
    if isa(l,Nil)
        if is(T,Any)
            print(io, "nil()")
        else
            print(io, "nil(", T, ")")
        end
    else
        print(io, "list(")
        while true
            show(io, head(l))
            l = tail(l)
            if isa(l,Cons)
                print(io, ", ")
            else
                break
            end
        end
        print(io, ")")
    end
end

list() = nil()
function list(elts...)
    l = nil()
    for i=length(elts):-1:1
        l = cons(elts[i],l)
    end
    return l
end
function list{T}(elts::T...)
    l = nil(T)
    for i=length(elts):-1:1
        l = cons(elts[i],l)
    end
    return l
end
function list{T}(elts::List{T}...)
    l = nil(List{T})
    for i=length(elts):-1:1
        l = cons(elts[i],l)
    end
    return l
end

length(l::Nil) = 0
length(l::Cons) = 1 + length(tail(l))

map(f::Base.Callable, l::Nil) = l
map(f::Base.Callable, l::Cons) = cons(f(head(l)), map(f, tail(l)))

copy(l::Nil) = l
copy(l::Cons) = cons(head(l), copy(tail(l)))

function append2(a, b)
    if isa(a,Nil)
        b
    else
        cons(head(a), append2(tail(a), b))
    end
end

cat(lst::List) = lst

function cat(lst::List, lsts...)
    n = length(lsts)
    l = lsts[n]
    for i = (n-1):-1:1
        l = append2(lsts[i], l)
    end
    return append2(lst, l)
end

start{T}(l::Cons{T}) = l
done{T}(l::Cons{T}, state::Cons{T}) = false
done{T}(l::Cons{T}, state::Nil{T}) = true
next{T}(l::Cons{T}, state::Cons{T}) = (state.head, state.tail)


text:
"""
    normalize_par_angles!(angles)

In-place version of [`normalize_par_angles`](@ref)
"""
function normalize_par_angles!(angles)
    # convert to positive
    @. angles[angles < 0] += 360
    idxs = sortperm(angles)

    correct = findfirst(d -> abs(d) > 180, diff(angles[idxs]))
    if !isnothing(correct)
        @. angles[angles < 180] += 360
    end

    return angles
end

"""
    normalize_par_angles(angles)

Ensures parallactic angle list (in degrees) is positive monotonic with no jumps greater than 180°.

# Examples
```jldoctest
julia> normalize_par_angles([-10, 20, 190])
3-element Array{Int64,1}:
 350
  20
 190
```
"""
normalize_par_angles(angles) = normalize_par_angles!(deepcopy(angles))


text:
using ModelingToolkit
using SteadyStateDiffEq
using Test

@parameters t r
@variables x(t)
@derivatives D'~t
eqs = [D(x) ~ x^2-r]
de = ODESystem(eqs)

for factor in [1e-1, 1e0, 1e10], u0_p in [(2.34,2.676),(22.34,1.632),(.3,15.676),(0.3,0.006)]
    u0 = [x => factor*u0_p[1]]
    p = [r => factor*u0_p[2]]
    ss_prob = SteadyStateProblem(de,u0,p)
    sol = solve(ss_prob,SSRootfind()).u[1]
    @test abs(sol^2 - factor*u0_p[2]) < 1e-8
    ss_prob = SteadyStateProblemExpr(de,u0,p)
    sol_expr = solve(eval(ss_prob),SSRootfind()).u[1]
    @test all(x->x==0,sol-sol_expr)
end


text:
export scalablexcorr, nodalxcorr, autocorrelate

function scalablexcorr(A::AbstractArray,Ntau::Int,kthresh::AbstractFloat)
    V,S,U = svd(A)
    kind = findall(S .> maximum(S) * 0.05)
    kmax = maximum(kind)

    Nt,Ns = size(A)

    # check that Ntau < Nt / 2
    Ntau = min(Ntau,Nt ÷ 2)

    # preallocate view
    V0 = @view V[Ntau÷2:Nt-Ntau÷2-1,1:kmax]
    U0 = @view U[:,1:kmax]
    if isa(V,AbstractGPUArray)
        VT0 = transpose(V[Ntau÷2:Nt-Ntau÷2-1,1:kmax])
        UT0 = transpose(U[:,1:kmax])
    else
        VT0 = V0'
        UT0 = U0'
    end

    Cout = similar(A,Ns,Ns,Ntau)

    for itau = -Ntau ÷ 2 + 1 : Ntau ÷ 2
        W = VT0 * V[itau + Ntau÷2:Nt + itau - Ntau÷2-1,1:kmax]
        Cout[:,:,itau+Ntau÷2] .= U0 * W * UT0
    end
    return Cout
end

function scalablexcorr(N::NodalData,Ntau::Real;kthresh=0.05)
    if isa(Ntau,AbstractFloat)
        Ntau = convert(Int,round(Ntau * N.fs[1]))
    end
    return scalablexcorr(N.data,Ntau,kthresh)
end

function nodalxcorr(N::NodalData,maxlag::Real)
    if isa(maxlag,AbstractFloat)
        maxlag = convert(Int,round(maxlag * N.fs[1]))
    end
    return nodalxcorr(N.data,maxlag)

end

function nodalxcorr(A::AbstractArray,maxlag::Int)
    Nt,Ns = size(A)
    Ncorr = Ns * (Ns -1) ÷ 2
    Cout = similar(A,maxlag * 2 + 1,Ncorr)
    FFTs = rfft(A,1)
    cstart = 0
    cend = 0
    for ii = 1:Ns-1
            cstart = cend + 1
            cend = cstart + Ns - ii - 1
            Cout[:,cstart:cend] .= correlate(FFTs[:,ii],FFTs[:,ii+1:end],Nt,maxlag)
    end
    return Cout
end

function autocorrelate(N::NodalData,maxlag::Real)
    if isa(maxlag,AbstractFloat)
        maxlag = convert(Int,round(maxlag * N.fs[1]))
    end
    return autocorrelate(N.data,maxlag)
end

function autocorrelate(A::AbstractArray,maxlag::Int64)
    Nt,Ns = size(A)
    FFT = rfft(A,1)
    corrT = irfft(conj.(FFT) .* FFT,Nt,1)
    # return corr[-maxlag:maxlag]
    t = vcat(0:Int(Nt  / 2)-1, -Int(Nt  / 2):-1)
    ind = findall(abs.(t) .<= maxlag)
    newind = fftshift(ind,1)
    return corrT[newind,:]
end


text:
#################### PDMatDistribution ####################

function unlist(d::PDMatDistribution, X::AbstractArray)
    n = dim(d)
    y = [X[i, j] for i = 1:n for j in 1:n if i >= j]
    y
end

function getix(row, col)::Int
    #assume left/lower triangular
    if col > row
        row = col
        col = row
    end
    (row - 1) * row / 2 + col
end

function relistlength(d::PDMatDistribution, X::AbstractArray{T}) where {T}
    n = dim(d)
    Y = [i >= j ? X[getix(i, j)] : X[getix(j, i)] for j = 1:n, i = 1:n]
    k = Int(n * (n + 1) / 2)
    (Y, k)
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
#=
  From Joost-Pieter Katoen
  "Probabilistic Programming Quantitative Modeling for the Masses?"
  (MMB 2018 Conference, Erlangen)

  Slide 4:
  """
  How likely does a student end up with a bad mood after getting
  a bad grade for an easy exam, given that she is well prepared?
  """

Distributions of variable difficulty
easy       =>    6022  (0.602200)
hard       =>    3978  (0.397800)

Distributions of variable preparation
good       =>   10000  (1.000000)

Distributions of variable grade
good       =>    6821  (0.682100)
bad        =>    3179  (0.317900)

Distributions of variable mood
bad        =>    5004  (0.500400)
good       =>    4996  (0.499600)

  Cf ~/webppl/student_mood_after_exam.wppl

=#
using Turing, StatsPlots, Distributions
include("jl_utils.jl")

@model function student_mood_after_exam()
    bad = 1
    good = 2
    badGood = [bad,good]

    easy = 1
    hard = 2

    difficulty ~ Categorical([0.6,0.4]) # [easy,hard]

    preparation ~ Categorical([0.7,0.3]) # badGood

    grade ~
        (difficulty==easy && preparation==bad)  ? Categorical([0.95,0.05]) : # badGood
        (difficulty==easy && preparation==good) ? Categorical([0.5,0.5])   :
        (difficulty==hard && preparation==bad)  ? Categorical([0.6,0.4])   :
        (difficulty==hard && preparation==good) ? Categorical([0.05,0.95]) : Dirac(easy)

    mood ~ grade == bad ? # badGood
        Categorical([0.9,0.1]) :
        Categorical([0.3,0.7])

    # true ~ Dirac(difficulty == easy)
    true ~ Dirac(preparation == good)
    # true ~ Dirac(grade == bad)
    # true ~ Dirac(mood == bad)

end

model = student_mood_after_exam()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

bad_good = ["bad","good"]
show_var_dist_pct(chns, :difficulty, ["easy","hard"])

println()
show_var_dist_pct(chns, :preparation, bad_good)

println()
show_var_dist_pct(chns, :grade, bad_good)

println()
show_var_dist_pct(chns, :mood, bad_good)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
import Base.Dates: now, unix2datetime

"""
    DateTime(::ZonedDateTime) -> DateTime

Returns an equivalent `DateTime` without any `TimeZone` information.
"""
DateTime(zdt::ZonedDateTime) = localtime(zdt)
@vectorize_1arg ZonedDateTime DateTime

"""
    now(::TimeZone) -> ZonedDateTime

Returns a `ZonedDateTime` corresponding to the user's system time in the specified `TimeZone`.
"""
function now(tz::TimeZone)
    utc = unix2datetime(time())
    ZonedDateTime(utc, tz, from_utc=true)
end

"""
    astimezone(zdt::ZonedDateTime, tz::TimeZone) -> ZonedDateTime

Converts a `ZonedDateTime` from its current `TimeZone` into the specified `TimeZone`.
"""
function astimezone end

function astimezone(zdt::ZonedDateTime, tz::VariableTimeZone)
    i = searchsortedlast(
        tz.transitions, zdt.utc_datetime,
        by=v -> typeof(v) == Transition ? v.utc_datetime : v,
    )

    if i == 0
        throw(NonExistentTimeError(localtime(zdt), tz))
    end

    zone = tz.transitions[i].zone
    return ZonedDateTime(zdt.utc_datetime, tz, zone)
end

function astimezone(zdt::ZonedDateTime, tz::FixedTimeZone)
    return ZonedDateTime(zdt.utc_datetime, tz, tz)
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function drawpolyboundingbox(apoly)
    gsave()
    setline(0.3)
    setdash("dotted")
    poly(convert(Vector{Point}, BoundingBox(apoly)), :stroke, close=true)
    box(BoundingBox(apoly), :stroke)
    grestore()
end

function drawpoly(p, x, y, counter)
    gsave()
    translate(x, y)
    sethue("purple")

    # label old points
    @layer begin
        sethue("grey50")
        fontsize(8)
        for (n, pt) in enumerate(p)
            circle(pt, 1, :fill)
            label(string(n), :NW, pt, offset=5)
        end
    end

    psort = polysortbyangle(p)
    drawpolyboundingbox(psort)
    poly(psort, close=true, :stroke)

    # label sorted points
    for (n, pt) in enumerate(psort)
        circle(pt, 1, :fill)
        text(string(n), pt.x, pt.y)
    end

    sethue("green")
    cp = polycentroid(psort)
    circle(cp, 5, :fill)
    text(string(counter), cp.x, cp.y)

    # highlight cases where centroid isn't inside poly. Might be an error...
    if ! isinside(cp, psort)
        gsave()
        setdash("dotted")
        setopacity(0.2)
        sethue("red")
        circle(cp, 100, :fillstroke)
        grestore()
    end

    grestore()
end

function polycentroidsort(width, height)
    pagetiles = Tiler(width, height, 5, 5, margin=50)
    tilesize = pagetiles.tilewidth/2
    for (pos, n) in pagetiles

        p = randompointarray(rand(-tilesize:-tilesize), rand(-tilesize:-tilesize), rand(tilesize:tilesize), rand(tilesize:tilesize), rand(5:12))

        # check for axis-aligned edges too
        rand(1:10) > 7 && (p = Random.shuffle!(box(O, 250, 300, vertices=true)))
        drawpoly(p, pos.x, pos.y, n)
    end
end

function polycentroidtest(fname)
    width, height = 3000, 3000
    Drawing(width, height, fname)
    fontsize(20)
    background("white")
    origin()
    polycentroidsort(width, height)
    @test finish() == true
end

fname = "polycentroidsort.png"
polycentroidtest(fname)

println("...finished test: output in $(fname)")


text:
"""
# module laminart

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

LAMINART equations from

R. D. S. Raizada and S. Grossberg, ‘Context-sensitive binding by the laminar circuits of V1 and V2: A unified model of perceptual grouping, attention, and orientation contrast’, Visual Cognition, vol. 8, no. 3–5, p. 431, Jun. 2001, doi: 10/d8phhc.

Uses CuArray and NNlib.conv with zeros boarding for convolution on GPU.
"""

module LaminartEqConv

using NNlib

function conv!(
    out::AbstractArray,
    img::AbstractArray,
    kern::AbstractArray,
    p::NamedTuple,
)
    @inbounds NNlib.conv!(
        out,
        img,
        kern,
        NNlib.DenseConvDims(
            img,
            kern,
            padding = (
                size(kern)[1] >> 1,
                size(kern)[1] >> 1,
                size(kern)[2] >> 1,
                size(kern)[2] >> 1,
            ),
            flipkernel = true,
        ),
    )
    return nothing
end

function fun_f!(arr::AbstractArray, p::NamedTuple)
    @inbounds @. arr = (p.μ * arr^p.n) / (p.ν_pw_n + arr^p.n)
    return nothing
end

# retina

function I_u!(r::AbstractArray, I::AbstractArray, p::NamedTuple)
    conv!(r, I, p.k_gauss_1, p)
    @. r = I - r
end

# LGN
function fun_dv!(
    dv::AbstractArray,
    v::AbstractArray,
    u::AbstractArray,
    x_lgn::AbstractArray,
    dv_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(dv_temp, x_lgn, p.k_gauss_1, p)
        @. dv =
            p.δ_v * (
                -v + ((1f0 - v) * max(u, 0f0) * (1f0 + p.C_1 * x_lgn)) -
                ((1f0 + v) * p.C_2 * dv_temp)
            )
    end
    return nothing
end

function fun_x_lgn!(x_lgn::AbstractArray, x::AbstractArray, p::NamedTuple)
    @inbounds NNlib.conv!(
        x_lgn,
        x,
        p.k_x_lgn,
        NNlib.DenseConvDims(x, p.k_x_lgn, padding = 0, flipkernel = true),
    )
    return nothing
end

function fun_v_C!(
    v_C::AbstractArray,
    v_p::AbstractArray,
    v_m::AbstractArray,
    V_temp_1::AbstractArray,
    V_temp_2::AbstractArray,
    Q_temp::AbstractArray,
    P_temp::AbstractArray,
    p::NamedTuple,
)

    @inbounds begin
        @. V_temp_2 = exp(-1.0f0 / 8.0f0) * (max(v_p, 0f0) - max(v_m, 0f0))
        conv!(V_temp_1, V_temp_2, p.k_gauss_2, p)
        conv!(Q_temp, V_temp_1, p.k_C_d, p)
        conv!(P_temp, V_temp_1, p.k_C_b, p)
        @. P_temp = abs(P_temp)
        @. v_C = p.γ * (max(Q_temp - P_temp, 0f0) + max(-Q_temp - P_temp, 0f0))
    end
    return nothing
end

# L6
function fun_dx_v1!(
    dx::AbstractArray,
    x::AbstractArray,
    C::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        @. dx =
            p.δ_c * (
                -x + (
                    (1.0f0 - x) * (
                        (p.α * C) + (p.ϕ * max(z - p.Γ, 0f0)) .+
                        (p.v_21 * x_v2) + p.att
                    )
                )
            )
    end
    return nothing
end

#     L4 excit
function fun_dy!(
    dy::AbstractArray,
    y::AbstractArray,
    C::AbstractArray,
    x::AbstractArray,
    m::AbstractArray,
    dy_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(dy_temp, m, p.k_W_p, p)
        @. dy_temp = m * dy_temp
        fun_f!(dy_temp, p)
        @. dy =
            p.δ_c *
            (-y + ((1f0 - y) * (C + (p.η_p * x))) - ((1f0 + y) * dy_temp))
    end
    return nothing
end

#     l4 inhib
function fun_dm!(
    dm::AbstractArray,
    m::AbstractArray,
    x::AbstractArray,
    dm_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(dm_temp, m, p.k_W_m, p)
        fun_f!(dm_temp, p)
        @. dm = p.δ_m * (-m + (p.η_m * x) - (m * dm_temp))
    end
    return nothing
end

#     L2/3 excit
function fun_dz!(
    dz::AbstractArray,
    z::AbstractArray,
    y::AbstractArray,
    H_z::AbstractArray,
    s::AbstractArray,
    dz_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(dz_temp, s, p.k_T_p, p)
        @. dz =
            p.δ_z * (
                -z + (
                    (1f0 - z) *
                    ((p.λ * max(y, 0f0)) + H_z + (p.a_23_ex * p.att))
                ) - ((z + p.ψ) * dz_temp)
            )
    end
    return nothing
end

#     L2/3 inhib
function fun_ds!(
    ds::AbstractArray,
    s::AbstractArray,
    H_z::AbstractArray,
    ds_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(ds_temp, s, p.k_T_m, p)
        @. ds = p.δ_s * (-s + H_z + (p.a_23_in * p.att) - (s * ds_temp))
    end
    return nothing
end

function fun_H_z!(
    H_z::AbstractArray,
    z::AbstractArray,
    H_z_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        @. H_z_temp = max(z - p.Γ, 0f0)
        conv!(H_z, H_z_temp, p.k_H, p)
    end
    return nothing
end

#     V2 L6
function fun_dx_v2!(
    dx::AbstractArray,
    x_v2::AbstractArray,
    z_v2::AbstractArray,
    z::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        @. dx =
            p.δ_c * (
                -x_v2 + (
                    (1f0 - x_v2) * (
                        (p.v12_6 * max(z - p.Γ, 0f0)) +
                        (p.ϕ * max(z_v2 - p.Γ, 0f0)) +
                        p.att
                    )
                )
            )
    end
    return nothing
end

# V2 L4 excit
function fun_dy_v2!(
    dy::AbstractArray,
    y_v2::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    m_v2::AbstractArray,
    dy_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(dy_temp, m_v2, p.k_W_p, p)
        fun_f!(dy_temp, p)
        @. dy =
            p.δ_c * (
                -y_v2 + (
                    (1.0f0 - y_v2) *
                    ((p.v12_4 * max(z - p.Γ, 0f0)) + (p.η_p * x_v2))
                ) - ((1.0f0 + y_v2) * dy_temp)
            )
    end
    return nothing
end

# # # for equilabrium

# # l6 equilabrium
function fun_x_equ!(
    x::AbstractArray,
    C::AbstractArray,
    z::AbstractArray,
    x_v2::AbstractArray,
    x_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        @. x_temp =
            (p.α * C) + (p.ϕ * max(z - p.Γ, 0f0)) + (p.v_21 * x_v2) + p.att
        @. x = x_temp / (1f0 + x_temp)
    end
    return nothing
end

# # l4 excit equilabrium
function fun_y_equ!(
    y::AbstractArray,
    C::AbstractArray,
    x::AbstractArray,
    m::AbstractArray,
    dy_temp::AbstractArray,
    p::NamedTuple,
)
    @inbounds begin
        conv!(dy_temp, m, p.k_W_p, p)
        @. dy_temp = m * dy_temp
        fun_f!(dy_temp, p)
        @. y = (C + (p.η_p * x) - dy_temp) / (1f0 + C + (p.η_p * x) + dy_temp)
    end
    return nothing
end

end


text:
struct ScaledLogistic{T} <: AbstractLink
    λ::T
end

(l::ScaledLogistic)(f::Real) = l.λ * logistic(f)

const AugPoisson = PoissonLikelihood{<:ScaledLogistic}

logλ(l::AugPoisson) = log(l.invlink.λ)

aux_field(::AugPoisson, Ω::NamedTuple) = values(Ω)
aux_field(::AugPoisson, Ω::TupleVector) = zip(Ω.ω, Ω.n)

function init_aux_variables(rng::AbstractRNG, ::AugPoisson, ndata::Int)
    return TupleVector((;
        ω=rand(rng, PolyaGamma(1, 0), ndata), n=rand(rng, Poisson(), ndata)
    ))
end

function init_aux_posterior(T::DataType, ::AugPoisson, n::Int)
    return For(TupleVector(; y=zeros(Int, n), c=zeros(T, n), λ=zeros(T, n))) do q
        PolyaGammaPoisson(q.y, q.c, q.λ)
    end
end

function aux_full_conditional(lik::AugPoisson, y::Int, f::Real)
    return PolyaGammaPoisson(y, abs(f), lik.invlink(-f))
end

function aux_posterior!(
    qΩ, lik::AugPoisson, y::AbstractVector{<:Int}, qf::AbstractVector{<:Normal}
)
    λ = lik.invlink.λ
    φ = only(qΩ.inds)
    @. φ.c = sqrt(second_moment(qf))
    @. φ.y = y
    @. φ.λ = λ * approx_expected_logistic(-mean(qf), φ.c)
    return qΩ
end

function auglik_potential(::AugPoisson, Ω, y::AbstractVector)
    return ((y .- Ω.n) / 2,) # use
end

function auglik_precision(::AugPoisson, Ω, ::AbstractVector)
    return (Ω.ω,)
end

function expected_auglik_potential(::AugPoisson, qΩ, y::AbstractVector)
    return ((y .- tvmean(qΩ).n) / 2,) # Short cut to get the mean n
end

function expected_auglik_precision(::AugPoisson, qΩ, ::AbstractVector)
    return (tvmean(qΩ).ω,) # It is not possible to shorten this since n is needed to compute ω
end

function expected_auglik_potential_and_precision(::AugPoisson, qΩ, y::AbstractVector)
    θ = tvmean(qΩ)
    return ((y .- θ.n) / 2,), (θ.ω,)
end

function logtilt(lik::AugPoisson, (ω, n)::Tuple{<:Real,<:Integer}, y::Integer, f::Real)
    return y * logλ(lik) - (y + n) * logtwo - logfactorial(y) +
           ((y - n) * f - abs2(f) * ω) / 2
end

function aux_prior(lik::AugPoisson, y::AbstractVector{<:Integer})
    λ = lik.invlink.λ
    return For(y) do yᵢ
        PolyaGammaPoisson(yᵢ, 0, λ)
    end
end

aux_prior(lik::AugPoisson, y::Integer) = PolyaGammaPoisson(y, 0, lik.invlink.λ)

function expected_logtilt(lik::AugPoisson, qΩ, y, qf::AbstractVector{<:Normal})
    logλ = log(lik.invlink.λ)
    return mapreduce(+, y, qf, @ignore_derivatives marginals(qΩ)) do yᵢ, qfᵢ, qω
        θ = ntmean(qω)
        m = mean(qfᵢ)
        return -(yᵢ + θ.n) * logtwo +
               ((yᵢ - θ.n) * m - (abs2(m) + var(qfᵢ)) * θ.ω) / 2 +
               yᵢ * logλ - logfactorial(yᵢ)
    end
end


text:
@testset "Loop chunk distribution" begin
using DiskArrays: DiskArrays, GridChunks, RegularChunks, IrregularChunks, AbstractDiskArray
using YAXArrayBase: YAXArrayBase
struct LargeDiskArray{N,CT<:GridChunks{N}} <: AbstractDiskArray{Float64,N}
    size::NTuple{N,Int}
    chunks::CT
    compressed::Bool
end
Base.size(a::LargeDiskArray) = a.size
DiskArrays.eachchunk(a::LargeDiskArray) = a.chunks
DiskArrays.haschunks(::LargeDiskArray) = DiskArrays.Chunked()
YAXArrayBase.iscompressed(a::LargeDiskArray) = a.compressed
s = (4000,2000,1500)
cs = (100,100,700)
YAXArrays.YAXDefaults.max_cache[] = 1.0e8
a1 = YAXArray(LargeDiskArray(s, GridChunks(s,cs),true))

#Test case where chunk has to be split
dc = mapslices(sum, a1, dims="Dim_1", debug = true)
ch = YAXArrays.DAT.getloopchunks(dc)

@test length(ch) == 2
@test ch[1] == RegularChunks(4,0,2000)
@test ch[2] == RegularChunks(700,0,1500)
dc.outcubes[1].cube
# Test that the allocated buffer is close to what the prescribes size
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc);
@test 0.5 < (sum(sizeof,incubes) + sum(sizeof,outcubes))/YAXArrays.YAXDefaults.max_cache[] <= 1.0;

#Test subsets and offset
a2 = a1[Dim_3=200..1300.5]
dc = mapslices(sum, a2, dims="Dim_1", debug = true)
ch = YAXArrays.DAT.getloopchunks(dc)
@test ch == (RegularChunks(4,0,2000), RegularChunks(700,200,1100))
# Test that the allocated buffer is close to what the prescribes size
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc);
@test 0.5 < (sum(sizeof,incubes) + sum(sizeof,outcubes))/YAXArrays.YAXDefaults.max_cache[] <= 1.0;

#Test with different max_cache
YAXArrays.YAXDefaults.max_cache[] = 2.0e8
dc = mapslices(sum, a2, dims="Dim_1", debug = true)
ch = YAXArrays.DAT.getloopchunks(dc)
#Test loop chunk sizes
@test ch == (RegularChunks(5,0,2000), RegularChunks(700,200,1100))
# Test that the allocated buffer is close to what the prescribes size
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc);
@test 0.5 < (sum(sizeof,incubes) + sum(sizeof,outcubes))/YAXArrays.YAXDefaults.max_cache[] <= 1.0# Test that the allocated buffer is close to what the prescribes size

# Test case that is chunk-friendly
YAXArrays.YAXDefaults.max_cache[] = 1.5e8
dc = mapslices(sum, a1, dims="Dim_3", debug = true)
ch = YAXArrays.DAT.getloopchunks(dc)
@test ch == (RegularChunks(100,0,4000), RegularChunks(100,0,2000))
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc)
@test 0.5 < (sum(sizeof,incubes) + sum(sizeof,outcubes))/YAXArrays.YAXDefaults.max_cache[] <= 1.0# Test that the allocated buffer is close to what the prescribes size

#With offset
a2 = a1[Dim_1=50..3050.5]
dc = mapslices(sum, a2, dims="Dim_3", debug = true)
ch = YAXArrays.DAT.getloopchunks(dc)
@test ch == (RegularChunks(100,50,3000), RegularChunks(100,0,2000))
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc)
@test 0.5 < (sum(sizeof,incubes) + sum(sizeof,outcubes))/YAXArrays.YAXDefaults.max_cache[] <= 1.0# Test that the allocated buffer is close to what the prescribes size

#With more working memory
YAXArrays.YAXDefaults.max_cache[] = 4.5e8
a2 = a1[Dim_1=50..3050.5]
dc = mapslices(sum, a1, dims="Dim_3", debug = true);
ch = YAXArrays.DAT.getloopchunks(dc)
@test ch == (RegularChunks(300,0,4000), RegularChunks(100,0,2000))
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc);
@test 0.5 < (sum(sizeof,incubes) + sum(sizeof,outcubes))/YAXArrays.YAXDefaults.max_cache[] <= 1.0# Test that the allocated buffer is close to what the prescribes size

#Now a completely different test, this is more for DataFrame functionality
# and tests if we can proagate irregular axes into the cubes
using Dates
yearchunks = map(y->isleapyear(y) ? 366 : 365, 2001:2008)
c = GridChunks(RegularChunks(100,0,4000), RegularChunks(100,0,2000), IrregularChunks(chunksizes=yearchunks))
s = last.(last(c))
YAXArrays.YAXDefaults.max_cache[] = 1.0e8
a1 = YAXArray(LargeDiskArray(s, c,true))

dc = mapCube(identity, a1, indims = InDims(), outdims = (), debug=true)
ch = YAXArrays.DAT.getloopchunks(dc)
@test ch == (RegularChunks(300,0,4000), RegularChunks(100,0,2000), RegularChunks(365,0,2922))
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc);
@test 0.5 < sum(sizeof,incubes)/YAXArrays.YAXDefaults.max_cache[] <= 1.0#

dc = mapCube(identity, a1, indims = InDims(), outdims = (), debug=true,irregular_loopranges = true)
ch = YAXArrays.DAT.getloopchunks(dc)
@test ch == (RegularChunks(300,0,4000), RegularChunks(100,0,2000), IrregularChunks(chunksizes=yearchunks))
incubes, outcubes = YAXArrays.DAT.getCubeCache(dc);
@test 0.5 < sum(sizeof,incubes)/YAXArrays.YAXDefaults.max_cache[] <= 1.0#
end

max_stars_repo_name: meggart/YAXArrays

max_forks_repo_name: meggart/YAXArrays


text:
import Term.Renderables: Renderable, RenderableText

lorem = "Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum."

println("Renderable text no style")
@time RenderableText(lorem)
@time RenderableText(lorem; width = 22)

println("Renderable text with style")
@time RenderableText(lorem; style = "red")
@time RenderableText(lorem; width = 22, style = "blue")

println("Re-create RT")
rt = RenderableText(lorem; style = "red")
@time RenderableText(rt; style = "red", width = rt.measure.w)
@time RenderableText(rt; style = "blue", width = rt.measure.w)

print("renderable of renderable")
ren = RenderableText(lorem; style = "red")
@time Renderable(ren);


text:
using ForwardDiff: gradient!

export update!, optimize!, optimize

"""
    update!(x, x̄)

Updates the array `x` according to `x .-= x̄`. `x` is of type `AbstractArray`
"""
function update!(x::Vector{T}, x̄) where {T<:AbstractFloat}
    x .-= x̄
end

"""
    update!(opt, x, g)

Updates array `x` according to gradient `g`
with respect to optimizer `opt`.
"""
function update!(opt, x, g)
    x .-= apply!(opt, x, g)
end

"""
    optimize!(x, opt, f, iter)

Uses optimizer `opt` to update `x` to try and minimize `f(x)` in `iter` iterations. This is the mutating version, which changes the input vector.
"""
function optimize!(
    x::Vector{T}, f::Function, opt; iter::Int=10, chunk::Union{Int, Nothing}=2) where {T <: AbstractFloat}
    grad = similar(x)
    cfg = grad_cfg(f,x ; chunk=chunk)
    for i in 1:iter
        gradient!(grad, f, x, cfg)
        update!(opt, x, grad)
    end
end

function optimize(x::Vector{T}, f::Function, opt; iter::Int=10, chunk::Union{Int, Nothing}=2) where {T <: AbstractFloat}
   y =copy(x)
   optimize!(y, f, opt; iter=iter, chunk=chunk)
   return y
end

function optimize!(A::Matrix{T}, x::Vector{T}, f::Function, opt; iter::Int=10, chunk::Union{Int, Nothing}=2) where {T <: AbstractFloat}
    if size(A)!=(size(x)[1]+1, iter)
        throw(ArgumentError("size of the matrix is not $(size(x)[1]+1)x$(iter)"))
    end
    grad = similar(x)
    cfg = grad_cfg(f, x; chunk=chunk)
    A[:, 1] = x
    for i in 1:iter
        A[:, i+1] = @view A[:, i]
        gradient!(grad, f, A[:, i+1], cfg)
        update!(opt, A[:, i+1], grad)
    end
end


text:
mutable struct LockedMinHeap{Heap,GetPriority,Priority}
    @const heap::Heap
    @const getpriority::GetPriority
    @atomic priority::Priority
    @atomic locked::Bool
end

lockedheap(
    ::Type{Eltype},
    ::Type{Priority} = Eltype;
    by = identity,
) where {Eltype,Priority} = LockedMinHeap(
    MutableBinaryHeap{Eltype}(Base.Order.By(by), EMPTY_VECTOR),
    by,
    typemax(Priority),
    false,
)

Base.eltype(h::LockedMinHeap) = eltype(typeof(h))
Base.eltype(::Type{LockedMinHeap{Heap}}) where {Heap} = eltype(Heap)

prioritytype(x) = prioritytype(typeof(x))
prioritytype(T::Type) = error("not defined for type ", T)
prioritytype(::Type{<:LockedMinHeap{<:Any,<:Any,Priority}}) where {Priority} = Priority

function Base.trylock(h::LockedMinHeap)
    if @atomic :monotonic h.locked
        return false
    end
    old = @atomicswap h.locked = true
    return !old
end

function Base.unlock(h::LockedMinHeap)
    @atomic h.locked = false
    return
end

function Base.push!(h::LockedMinHeap, x)
    @DBG @check h.getpriority(x) != typemax(prioritytype(h))
    push!(h.heap, x)
    @atomic :monotonic h.priority = h.getpriority(first(h.heap))
    return h
end

function Base.pop!(h::LockedMinHeap)
    x = pop!(h.heap)
    if isempty(h.heap)
        p = typemax(prioritytype(h))
    else
        p = h.getpriority(first(h.heap))
    end
    @atomic :monotonic h.priority = p
    return x
end

struct MultiQueue{Heap<:LockedMinHeap}
    heaps::Vector{Heap}
end

Base.eltype(q::MultiQueue) = eltype(typeof(q))
Base.eltype(::Type{MultiQueue{Heap}}) where {Heap} = eltype(Heap)

prioritytype(::Type{MultiQueue{Heap}}) where {Heap} = prioritytype(Heap)

function Base.push!(multiq::MultiQueue, x)
    nspins = 0
    while true
        nspins += 1
        h = rand(multiq.heaps)
        if trylock(h)
            try
                push!(h, x)
                return h
            finally
                unlock(h)
            end
        end
        if nspins > 1000
            errorwith(; multiq = multiq, x = x) do io, (; multiq, x)
                print(io, "Failed to push an element `", x, "` to ")
                show(io, MIME"text/plain"(), multiq)
            end
        end
    end
end

struct Contention end
struct LikelyEmpty end

function trypop!(multiq::MultiQueue)
    nempties = 0
    while true
        h1 = rand(multiq.heaps)
        h2 = rand(multiq.heaps)
        p1 = @atomic :monotonic h1.priority
        p2 = @atomic :monotonic h2.priority
        if p2 < p1
            p1, p2 = p2, p1
            h1, h2 = h2, h1
        end
        if p1 == typemax(prioritytype(multiq))
            nempties += 1
            if nempties > length(multiq.heaps)
                break
            else
                continue
            end
        end
        if trylock(h1)
            try
                p1 = @atomic :monotonic h1.priority
                if p1 == typemax(prioritytype(multiq))
                    continue
                end
                return Some{eltype(h1)}(pop!(h1))
            finally
                unlock(h1)
            end
        end
        # nempties = 0
    end

    lockfailed = false
    for h in multiq.heaps
        p = @atomic :monotonic h.priority
        if p == typemax(prioritytype(multiq))
            continue
        end
        if trylock(h)
            try
                p = @atomic :monotonic h.priority
                if p == typemax(prioritytype(multiq))
                    continue
                end
                return Some{eltype(h)}(pop!(h))
            finally
                unlock(h)
            end
        else
            lockfailed = true
        end
    end
    if lockfailed
        return Contention()
    else
        return LikelyEmpty()
    end
end


text:
#######################################################################
# functions to makle lags and leads
#######################################################################

"""
    F!(df, v, horizons=1, fun=identity)

Forward operator.
Applies function `fun` to the column `v` of `df` and then makes the
leads specified by horizons and adds them to `df` as the column `F(i)(v)`.
Modifies `df`.
Return the list of new columns created.
"""
function F!(df, v, horizons=1, fun=identity)
    # Forward operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("F$(i)$(v)")
        df[:, Fiv] = lead( map(fun, df[!, v]), i)
        push!(newvars, Fiv)
    end
    return newvars
end

"""
    L!(df, v, horizons=1, fun=identity)

Lag operator.
Applies function `fun` to the column `v` of `df` and then makes the
laggs specified by horizons and adds them to `df` as the column `L(i)(v)`.
Modifies `df`.
Return the list of new columns created.
"""
function L!(df, v, horizons=1, fun=identity)
    # Lag operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("L$(i)$(v)")
        df[:, Fiv] = lag(map(fun, df[!, v]), i)
        push!(newvars, Fiv)
    end
    return newvars
end

"""
    F!(df, v, horizons=1, fun=identity, baseperiod=0, relativebase=false)

Cumlative forward difference operator
Applies function `fun` to the column `v` of `df` and then calculates
the growth of this column, saving the results to a new columns in
`df` as the column `FΔ(i)(v)`.
If `x = map(fun, df[!, v])`, `df.FΔ(i)(v) = lead(x, horizons) - lead(x, baseperiod)`
or `x = map(fun, df[!, v])`, `df.FΔ(i)(v) = lead(x, horizons) - lead(x, horizon + baseperiod)`
if relativebase == true
Modifies `df`.
Return the list of new columns created.
"""
function FΔ!(df, v, horizons=1, fun=identity, baseperiod=0, relativebase=false)
    # Forward difference operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("FΔ$(i)$(v)")
        if relativebase
            df[:, Fiv] = lead( map(fun, df[!, v]), i) .- lead(map(fun, df[!, v]), i+baseperiod)
        else
            df[:, Fiv] = lead( map(fun, df[!, v]), i) .- lead(map(fun, df[!, v]), baseperiod)
        end
        push!(newvars, Fiv)
    end
    return newvars
end

"""
    LΔ!(df, v, horizons=1, fun=identity, baseperiod=0, relativebase=false)
Cumlative lag difference operator

Applies function `fun` to the column `v` of `df` and then calculates
the lagged growth of this column, saving the results to a new columns in
`df` as the column `LΔ(i)(v)`.
If `x = map(fun, df[!, v])`, `df.FΔ(i)(v) = lag(x, baseperiod) - lag(x, horizons)`
or `x = map(fun, df[!, v])`, `df.FΔ(i)(v) = lag(x, i+baseperiod) - lag(x, horizons)`
if relativebase == true
Modifies `df`.
Return the list of new columns created.
"""
function LΔ!(df, v, horizons=1, fun=identity, baseperiod=0, relativebase=false)
    # Lag difference operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("LΔ$(i)$(v)")
        if relativebase
            df[:, Fiv] = lag(map(fun, df[!, v]), i+baseperiod) .-  lag( map(fun, df[!, v]), i)
        else
            df[:, Fiv] = lag(map(fun, df[!, v]), i+baseperiod) .-  lag( map(fun, df[!, v]), i)
        end
        push!(newvars, Fiv)
    end
    return newvars
end

"""
    Δ(x, horizons=1, fun=identity)

Single period growth function
    `map(fun, x) .-  lag(map(fun, x), horizons)`
"""
function Δ(x, horizons=1, fun=identity)
    return map(fun, x) .-  lag(map(fun, x), horizons)
end


text:
using NeuralAttentionlib
using Documenter

DocMeta.setdocmeta!(NeuralAttentionlib, :DocTestSetup, :(using NeuralAttentionlib); recursive=true)

makedocs(;
    modules=[NeuralAttentionlib],
    authors="chengchingwen <adgjl5645@hotmail.com> and contributors",
    repo="https://github.com/chengchingwen/NeuralAttentionlib.jl/blob/{commit}{path}#{line}",
    sitename="NeuralAttentionlib.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chengchingwen.github.io/NeuralAttentionlib.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Terminology" => "term.md",
        "Example" => "example.md",
        "API Reference" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/chengchingwen/NeuralAttentionlib.jl",
)


text:
struct MDOptions
    n_iters::Int
    n_save_iters::Int
    n_save::Int
end

"""
    MDOptions(;n_iters = 10^4, n_save_iters=1)

Set options for samplers.

### Fields

* n_iters       - Set the number of iterations of the sampler
* n_save_iters  - Set the frequency at which iterations are saved.  If
                  n_save_iters=1, every iteration is saved.  If n_save_iters=n_iters,
                  only the final iteration is saved.
"""
function MDOptions(; n_iters = 10^4, n_save_iters = 1)

    return MDOptions(n_iters, n_save_iters, floor(Int, n_iters / n_save_iters))
end

"""
    Boltzmann_likelihood(x, V, β)

Compute the unnormalized Boltzmann density, exp(-β V(x))
for potential `V` at inverse temperature `β`

### Fields

* x     - Value at which to evalute the density
* V     - Potential
* β     - Inverse temperature
"""
function Boltzmann_likelihood(x, V, β)
    w = exp(-β * V(x))
    return w
end

"""
`trivial_constraint!` - Trivial constraint function
### Fields
* `state` - Current state of the sampler
* `i`     - Index of current iterate
"""
function trivial_constraint!(state, i)
    state
end


text:
function open_plexoszip(zippath::String, xmlname::String=defaultxml(zippath))

    resultsarchive = _open_plexoszip(zippath)
    xml = parsexml(resultsarchive[xmlname])

    data = PLEXOSSolutionDataset(xml)
    resultvalues = perioddata(resultsarchive)

    return data, resultvalues

end

function _open_plexoszip(zippath::String)
    isfile(zippath) || error("$zippath does not exist")
    archive = open_zip(zippath)
    return archive
end

defaultxml(zippath::String) = replace(basename(zippath), r".zip$"=>".xml")

function perioddata(archive::Archive)
    results = Dict{Int,Vector{UInt8}}()
    for filename in keys(archive)
        rgx = match(r"t_data_(\d).BIN", filename)
        isnothing(rgx) && continue
        data = archive[filename]
        results[parse(Int, rgx[1])] = data
    end
    return results
end

function getchildtext(name::String, e::Node)
    resultnode = findfirst("x:" * name, e, ["x"=>namespace(e)])
    isnothing(resultnode) && error("$e does not have child $name")
    return nodecontent(resultnode)
end

getchildfloat(name::String, e::Node) = parse(Float64, getchildtext(name, e))
getchildint(name::String, e::Node) = parse(Int, getchildtext(name, e))
getchildbool(name::String, e::Node) = parse(Bool, getchildtext(name, e))