text:
# Used for objectives and solvers where the gradient and Hessian is available/exists
mutable struct TwiceDifferentiable{T,TDF,TH,TX} <: AbstractObjective
    f
    df
    fdf
    h
    F::T
    DF::TDF
    H::TH
    x_f::TX
    x_df::TX
    x_h::TX
    f_calls::Vector{Int}
    df_calls::Vector{Int}
    h_calls::Vector{Int}
end
# compatibility with old constructor
function TwiceDifferentiable(f, g, fg, h, x::TX, F::T = real(zero(eltype(x))), G::TG = alloc_DF(x, F), H::TH = alloc_H(x, F); inplace = true) where {T, TG, TH, TX}
    x_f, x_df, x_h = x_of_nans(x), x_of_nans(x), x_of_nans(x)

    g! = df!_from_df(g, F, inplace)
    fg! = fdf!_from_fdf(fg, F, inplace)
    h! = h!_from_h(h, F, inplace)

    TwiceDifferentiable{T,TG,TH,TX}(f, g!, fg!, h!,
                                        copy(F), copy(G), copy(H),
                                        x_f, x_df, x_h,
                                        [0,], [0,], [0,])
end

function TwiceDifferentiable(f, g, h,
                             x::AbstractVector{TX},
                             F::Real = real(zero(eltype(x))),
                             G = alloc_DF(x, F),
                             H = alloc_H(x, F); inplace = true) where {TX}
    g! = df!_from_df(g, F, inplace)
    h! = h!_from_h(h, F, inplace)

    fg! = make_fdf(x, F, f, g!)
    x_f, x_df, x_h = x_of_nans(x), x_of_nans(x), x_of_nans(x)

    return TwiceDifferentiable(f, g!, fg!, h!, F, G, H, x_f, x_df, x_h, [0,], [0,], [0,])
end

function TwiceDifferentiable(f, g,
                             x_seed::AbstractVector{T},
                             F::Real = real(zero(T)); autodiff = :finite, inplace = true) where T
    n_x = length(x_seed)

    g! = df!_from_df(g, F, inplace)
    fg! = make_fdf(x_seed, F, f, g!)

    if is_finitediff(autodiff)

        # Figure out which Val-type to use for FiniteDiff based on our
        # symbol interface.
        fdtype = finitediff_fdtype(autodiff)

        jcache = FiniteDiff.JacobianCache(x_seed, fdtype)
        function h!(storage, x)
            FiniteDiff.finite_difference_jacobian!(storage, g!, x, jcache)
            return
        end

    elseif is_forwarddiff(autodiff)
        hcfg = ForwardDiff.HessianConfig(f, copy(x_seed))
        h! = (out, x) -> ForwardDiff.hessian!(out, f, x, hcfg)
    else
        error("The autodiff value $(autodiff) is not supported. Use :finite or :forward.")
    end
    TwiceDifferentiable(f, g!, fg!, h!, x_seed, F)
end

TwiceDifferentiable(d::NonDifferentiable, x_seed::AbstractVector{T} = d.x_f, F::Real = real(zero(T)); autodiff = :finite) where {T<:Real} =
    TwiceDifferentiable(d.f, x_seed, F; autodiff = autodiff)

function TwiceDifferentiable(d::OnceDifferentiable, x_seed::AbstractVector{T} = d.x_f,
                             F::Real = real(zero(T)); autodiff = :finite) where T<:Real
    if is_finitediff(autodiff)

        # Figure out which Val-type to use for FiniteDiff based on our
        # symbol interface.
        fdtype = finitediff_fdtype(autodiff)

        jcache = FiniteDiff.JacobianCache(x_seed, fdtype)
        function h!(storage, x)
            FiniteDiff.finite_difference_jacobian!(storage, d.df, x, jcache)
            return
        end
    elseif is_forwarddiff(autodiff)
        hcfg = ForwardDiff.HessianConfig(d.f, copy(gradient(d)))
        h! = (out, x) -> ForwardDiff.hessian!(out, d.f, x, hcfg)
    else
        error("The autodiff value $(autodiff) is not supported. Use :finite or :forward.")
    end
    return TwiceDifferentiable(d.f, d.df, d.fdf, h!, x_seed, F, gradient(d))
end

function TwiceDifferentiable(f, x::AbstractVector, F::Real = real(zero(eltype(x)));
                             autodiff = :finite, inplace = true)
    if is_finitediff(autodiff)

        # Figure out which Val-type to use for FiniteDiff based on our
        # symbol interface.
        fdtype = finitediff_fdtype(autodiff)
        gcache = FiniteDiff.GradientCache(x, x, fdtype)

        function g!(storage, x)
            FiniteDiff.finite_difference_gradient!(storage, f, x, gcache)
            return
        end
        function fg!(storage::Vector, x::Vector)
            g!(storage, x)
            return f(x)
        end

        function h!(storage::Matrix, x::Vector)
            FiniteDiff.finite_difference_hessian!(storage, f, x)
            return
        end
    elseif is_forwarddiff(autodiff)

        gcfg = ForwardDiff.GradientConfig(f, x)
        g! = (out, x) -> ForwardDiff.gradient!(out, f, x, gcfg)

        fg! = (out, x) -> begin
            gr_res = DiffResults.DiffResult(zero(eltype(x)), out)
            ForwardDiff.gradient!(gr_res, f, x, gcfg)
            DiffResults.value(gr_res)
        end

        hcfg = ForwardDiff.HessianConfig(f, x)
        h! = (out, x) -> ForwardDiff.hessian!(out, f, x, hcfg)
    else
        error("The autodiff value $(autodiff) is not supported. Use :finite or :forward.")
    end
    TwiceDifferentiable(f, g!, fg!, h!, x, F)
end

function hv_product!(obj::TwiceDifferentiable, x, v)
    H = hessian!(obj, x)
    return H*v
end


text:
@doc raw"""
First twenty roots of Bessel funcion ``J_0`` in Float64.
https://mathworld.wolfram.com/BesselFunctionZeros.html

# Examples
```jldoctest
julia> zeros = besselj0.(FastGaussQuadrature.BESSELJ0_ROOTS);

julia> all(zeros .< 1e-14)
true
```
"""
const BESSELJ0_ROOTS = @SVector [
    2.4048255576957728,
    5.5200781102863106,
    8.6537279129110122,
    11.791534439014281,
    14.930917708487785,
    18.071063967910922,
    21.211636629879258,
    24.352471530749302,
    27.493479132040254,
    30.634606468431975,
    33.775820213573568,
    36.917098353664044,
    40.058425764628239,
    43.199791713176730,
    46.341188371661814,
    49.482609897397817,
    52.624051841114996,
    55.765510755019979,
    58.906983926080942,
    62.048469190227170,
]

@doc raw"""
Coefficients of Chebyshev series approximations for the zeros of the Bessel functions
```math
j_{\nu, s} \approx \sum_{k}^{n} C_{k,s} T_k(\frac{\nu-2}{3})
```
where $j_{\nu, s}$ is a ``s``-th zero of Bessel function ``J_\nu``, ``\{T_k\}`` are Chebyshev polynomials and ``C_{k, s}`` is the coefficients.
"""
const PIESSENS_C = @SMatrix [
       2.883975316228  8.263194332307 11.493871452173 14.689036505931 17.866882871378 21.034784308088
       0.767665211539  4.209200330779  4.317988625384  4.387437455306  4.435717974422  4.471319438161
      -0.086538804759 -0.164644722483 -0.130667664397 -0.109469595763 -0.094492317231 -0.083234240394
       0.020433979038  0.039764618826  0.023009510531  0.015359574754  0.011070071951  0.008388073020
      -0.006103761347 -0.011799527177 -0.004987164201 -0.002655024938 -0.001598668225 -0.001042443435
       0.002046841322  0.003893555229  0.001204453026  0.000511852711  0.000257620149  0.000144611721
      -0.000734476579 -0.001369989689 -0.000310786051 -0.000105522473 -0.000044416219 -0.000021469973
       0.000275336751  0.000503054700  0.000083834770  0.000022761626  0.000008016197  0.000003337753
      -0.000106375704 -0.000190381770 -0.000023343325 -0.000005071979 -0.000001495224 -0.000000536428
       0.000042003336  0.000073681222  0.000006655551  0.000001158094  0.000000285903  0.000000088402
      -0.000016858623 -0.000029010830 -0.000001932603 -0.000000269480 -0.000000055734 -0.000000014856
       0.000006852440  0.000011579131  0.000000569367  0.000000063657  0.000000011033  0.000000002536
      -0.000002813300 -0.000004672877 -0.000000169722 -0.000000015222 -0.000000002212 -0.000000000438
       0.000001164419  0.000001903082  0.000000051084  0.000000003677  0.000000000448  0.000000000077
      -0.000000485189 -0.000000781030 -0.000000015501 -0.000000000896 -0.000000000092 -0.000000000014
       0.000000203309  0.000000322648  0.000000004736  0.000000000220  0.000000000019  0.000000000002
      -0.000000085602 -0.000000134047 -0.000000001456 -0.000000000054 -0.000000000004               0
       0.000000036192  0.000000055969  0.000000000450  0.000000000013               0               0
      -0.000000015357 -0.000000023472 -0.000000000140 -0.000000000003               0               0
       0.000000006537  0.000000009882  0.000000000043  0.000000000001               0               0
      -0.000000002791 -0.000000004175 -0.000000000014               0               0               0
       0.000000001194  0.000000001770  0.000000000004               0               0               0
      -0.000000000512 -0.000000000752               0               0               0               0
       0.000000000220  0.000000000321               0               0               0               0
      -0.000000000095 -0.000000000137               0               0               0               0
       0.000000000041  0.000000000059               0               0               0               0
      -0.000000000018 -0.000000000025               0               0               0               0
       0.000000000008  0.000000000011               0               0               0               0
      -0.000000000003 -0.000000000005               0               0               0               0
       0.000000000001  0.000000000002               0               0               0               0
]

@doc raw"""
Values of Bessel function ``J_1`` on first ten roots of Bessel function `J_0`.

# Examples
```jldoctest
julia> roots = approx_besselroots(0,10);

julia> (besselj1.(roots)).^2 ≈ FastGaussQuadrature.BESSELJ1_ON_BESSELJ0_ROOTS
true
```
"""
const BESSELJ1_ON_BESSELJ0_ROOTS = @SVector [
    0.2695141239419169,
    0.1157801385822037,
    0.07368635113640822,
    0.05403757319811628,
    0.04266142901724309,
    0.03524210349099610,
    0.03002107010305467,
    0.02614739149530809,
    0.02315912182469139,
    0.02078382912226786,
]

@doc raw"""
The first 11 roots of the Airy function in Float64 precision
https://mathworld.wolfram.com/AiryFunctionZeros.html

# Examples
```jldoctest
julia> zeros = airy.(FastGaussQuadrature.AIRY_ROOTS);

julia> all(zeros .< 1e-14)
true
```
"""
const AIRY_ROOTS = @SVector [
    -2.338107410459767,
    -4.08794944413097,
    -5.520559828095551,
    -6.786708090071759,
    -7.944133587120853,
    -9.022650853340981,
    -10.04017434155809,
    -11.00852430373326,
    -11.93601556323626,
    -12.828776752865757,
    -13.69148903521072,
]

ext: jl


text:
# Code taken from FluxML/Optimisers.jl

const Param{T<:Number} = Union{AbstractArray{T},T}

_apply(opt, x, x̄, state) = apply(opt, x, x̄, state)
_apply(opt, x, x̄, ::Nothing) = apply(opt, x, x̄)

# Immutable updates
update!(opt, x, xb, state=nothing) = _update!(opt, x, xb, state)
update(opt, x, xb, state=nothing)  = _update(opt, x, xb, state)

function _update(opt, x::Param, x̄::Param, state = nothing)
  Δ, state = _apply(opt, x, x̄, state)
  return x .- Δ, state
end

# Mutable updates

# Figure out if we can do in-place
inplace(x, y) = false
inplace(x, y::Nothing) = true
inplace(x::AbstractArray, x̄::AbstractArray) = true
inplace(x, x̄::NamedTuple) = all(inplace(getfield(x, f), getfield(x̄, f)) for f in fieldnames(typeof(x̄)))

function _update!(opt, x::AbstractArray{<:Number}, x̄::AbstractArray, state = nothing)
  Δ, state = _apply(opt, x, x̄, state)
  #println("fr = ", x)
  #println("to = ", Δ)
  x .-= Δ
  #println("after = ", x)
  return state
end

function _update!(opt, x::Tuple, x̄::Tuple, state=nothing)
  for f in propertynames(x̄)
    f̄ = getfield(x̄, f)
#  println(f,"-->",f̄)
    f̄ === nothing || _update!(opt, getfield(x, f), f̄, state)
  end
end

function _update!(opt, x, x̄::NamedTuple, state=nothing)
  for f in propertynames(x̄)
    f == :tuple_all_weights && continue
    f̄ = getproperty(x̄, f)
#    println(f,"-->",f̄)
    # println(f,"-->", typeof(getfield(x, f)))

    f̄ === nothing || _update!(opt, getfield(x, f), f̄, state)
  end
end

function _update!(opt, x, x̄, state=nothing)
  for f in propertynames(x̄)
    f̄ = getproperty(x̄, f)
#    println(f,"-->",f̄)
    # println(f,"-->", typeof(getfield(x, f)))

    f̄ === nothing || _update!(opt, getfield(x, f), f̄, state)
  end
end

# Package Integration

using Requires

@init @require Colors="5ae59095-9a9b-59fe-a467-6f913c188581" begin
  function _update(opt, x::Colors.RGB{T}, x̄::NamedTuple) where T
    Colors.RGB{T}(clamp(_update(opt, x.r, x̄.r)[1], 0, 1),
                  clamp(_update(opt, x.g, x̄.g)[1], 0, 1),
                  clamp(_update(opt, x.b, x̄.b)[1], 0, 1)), nothing
  end
end


text:
module Diban

######################## begin:Taken from Parquet.jl
using Thrift
using ProtoBuf
using Snappy
using CodecZlib, CodecZstd
using MemPool

import Base: show, open, close, values
import Thrift: isfilled

# export is_par_file, ParFile, show, nrows, ncols, rowgroups, columns, pages, bytes, values, colname, colnames
# export SchemaConverter, schema, JuliaConverter, ThriftConverter, ProtoConverter
# export RowCursor, ColCursor, RecCursor
# export AbstractBuilder, JuliaBuilder

# package code goes here
include("PAR2/PAR2.jl")
using .PAR2
include("codec.jl")
include("schema.jl")
include("reader.jl")
include("cursor.jl")
include("show.jl")
######################## end:Taken from Parquet.jl

# package code goes here
include("encoding.jl")
# include("thrift/parquet.jl")
include("metadata.jl")

include("column_reader.jl")
include("read_parquet.jl")

# include("BitPackedIterator.jl")

include("writer_consts.jl")
include("writer.jl")

# from column_reader.jl
# export read_column

# from metadata.jl
# export metadata

# from read_parquet.jl
export read_parquet

# # from BitPackedIterator
# export BitPackedIterator
# export iterate, length

# from writer.jl
export write_parquet

end # module


text:
using CompScienceMeshes
using Test

using StaticArrays

tet = @SVector[:a,:b,:c,:d]

@test CompScienceMeshes.relorientation(@SVector[:a,:b,:c], tet) == 4
@test CompScienceMeshes.relorientation(@SVector[:c,:a,:b], tet) == 4
@test CompScienceMeshes.relorientation(@SVector[:b,:c,:a], tet) == 4
@test CompScienceMeshes.relorientation(@SVector[:a,:c,:b], tet) == -4
@test CompScienceMeshes.relorientation(@SVector[:b,:a,:c], tet) == -4
@test CompScienceMeshes.relorientation(@SVector[:c,:b,:a], tet) == -4

@test CompScienceMeshes.relorientation(@SVector[:d,:c,:b], tet) == 1
@test CompScienceMeshes.relorientation(@SVector[:b,:d,:c], tet) == 1
@test CompScienceMeshes.relorientation(@SVector[:c,:b,:d], tet) == 1
@test CompScienceMeshes.relorientation(@SVector[:c,:d,:b], tet) == -1
@test CompScienceMeshes.relorientation(@SVector[:b,:c,:d], tet) == -1
@test CompScienceMeshes.relorientation(@SVector[:d,:b,:c], tet) == -1

@test CompScienceMeshes.relorientation(@SVector[:a,:c,:d], tet) == 2
@test CompScienceMeshes.relorientation(@SVector[:d,:a,:c], tet) == 2
@test CompScienceMeshes.relorientation(@SVector[:c,:d,:a], tet) == 2
@test CompScienceMeshes.relorientation(@SVector[:c,:a,:d], tet) == -2
@test CompScienceMeshes.relorientation(@SVector[:d,:c,:a], tet) == -2
@test CompScienceMeshes.relorientation(@SVector[:a,:d,:c], tet) == -2

@test CompScienceMeshes.relorientation(@SVector[:a,:d,:b], tet) == 3
@test CompScienceMeshes.relorientation(@SVector[:b,:a,:d], tet) == 3
@test CompScienceMeshes.relorientation(@SVector[:d,:b,:a], tet) == 3
@test CompScienceMeshes.relorientation(@SVector[:d,:a,:b], tet) == -3
@test CompScienceMeshes.relorientation(@SVector[:b,:d,:a], tet) == -3
@test CompScienceMeshes.relorientation(@SVector[:a,:b,:d], tet) == -3


text:
@testset "Maps and Reductions" begin
    @testset "sum" begin
        @testset "Vector" begin
            rng, M = MersenneTwister(123456), 3
            frule_test(sum, (randn(rng, M), randn(rng, M)))
            rrule_test(sum, randn(rng), (randn(rng, M), randn(rng, M)))
        end
        @testset "Matrix" begin
            rng, M, N = MersenneTwister(123456), 3, 4
            frule_test(sum, (randn(rng, M, N), randn(rng, M, N)))
            rrule_test(sum, randn(rng), (randn(rng, M, N), randn(rng, M, N)))
        end
        @testset "Array{T, 3}" begin
            rng, M, N, P = MersenneTwister(123456), 3, 7, 11
            frule_test(sum, (randn(rng, M, N, P), randn(rng, M, N, P)))
            rrule_test(sum, randn(rng), (randn(rng, M, N, P), randn(rng, M, N, P)))
        end
        @testset "keyword arguments" begin
            rng = MersenneTwister(33)
            n = 4
            X = randn(rng, n, n+1)
            y, pullback = rrule(sum, X; dims=2)
            ȳ = randn(rng, size(y))
            _, x̄_ad = pullback(ȳ)
            x̄_fd = only(j′vp(central_fdm(5, 1), x->sum(x, dims=2), ȳ, X))
            @test x̄_ad ≈ x̄_fd atol=1e-9 rtol=1e-9
        end
    end  # sum
end


text:
import FileIO:save
using DataFrames, DataFramesMeta, PyPlot

# Functions to plot gridsearch (optimization) results within julia.
# ------ Not used eventually for the paper, only for inspection.
function plot_fitness_multi(df::DataFrame, changeprobability::Float64
                    ;savepath = homedir(),
                    col = :meanterrors,
                    xlabelstring = "Leaky - Multi",
                    ylabelstring = "Mean error",
                    ylims = [0., 2.],
                    xlims = [0., 1.],
                    titlestring = "Pc: "*string(changeprobability),
                    filestr = "t"
                    )
    fig = figure(); ax = gca();
    stochvalues = unique(df[:stoch])
    clist = [cmap(i) for i in range(0, 1, length = length(stochvalues))]
    for i in length(stochvalues):-1:1
        plot_fitness_core_multi(df, changeprobability, stochvalues[i], ax,
                        col=col, color = clist[i],
                        xlabelstring = xlabelstring,
                        ylabelstring = ylabelstring,
                        ylims = ylims, xlims = xlims,
                        titlestring = titlestring)
    end
    PyPlot.savefig(joinpath(savepath,
        "gridsearchMulti_"* string(changeprobability)* "_" * xlabelstring *
        "_" * filestr *".png"))
end
function plot_fitness_gaussian(df::DataFrame, changeprobability::Float64
                    ;savepath = homedir(),
                    col=:meanmerrors,
                    xlabelstring = "Leak - Gaussian",
                    ylabelstring = "Mean error",
                    ylims = [0., 10.],
                    xlims = [0., 1.],
                    filestr = "g",
                    titlestring = "Pc: "*string(changeprobability))

    fig = figure(); ax = gca();
    sigmavalues = unique(df[:sigma])
    clist = [cmap(i) for i in range(0, 1, length = length(sigmavalues))]
    for i in length(sigmavalues):-1:1
        plot_fitness_core_gaussian(df, changeprobability, sigmavalues[i], ax,
                            col = col,
                            color = clist[i],
                            xlabelstring = xlabelstring,
                            ylabelstring = ylabelstring,
                            ylims = ylims, xlims = xlims,
                            titlestring = titlestring)
    end
    PyPlot.savefig(joinpath(savepath,
        "gridsearchGaussian_"* string(changeprobability) * "_" * xlabelstring *
        "_" * filestr * ".png"))
end
function plot_fitness_core_multi(df::DataFrame, changeprobability::Float64,
                            stochvalue::Float64,
                            ax; col = :meanterrors,
                            color = "b",
                            xlabelstring = "",
                            ylabelstring = "Mean error",
                            ylims = [0., 2.], #[0., 10^6],
                            xlims = [-0.03, 1.],
                            titlestring = "Pc: "*string(changeprobability))

    dftrimmed = @where(df, :changeprob .== changeprobability,
                            :stoch .== stochvalue)
    ax.plot(dftrimmed[:param], dftrimmed[col], #color = color,
            label = string(stochvalue))
    ylabel(ylabelstring); xlabel(xlabelstring); title(titlestring)
    ax.grid(true); ax.legend()
    ax.set_ylim(ylims[1], ylims[2]); ax.set_xlim(xlims[1], xlims[2])
end
function plot_fitness_core_gaussian(df::DataFrame, changeprobability::Float64,
                            sigmavalue::Float64,
                            ax; col = :meanmerrors,
                            color = "b",
                            xlabelstring = " ",
                            ylabelstring = "Mean error",
                            ylims = [0., 10.],
                            xlims = [-0.03, 1.],
                            titlestring = "Pc: "*string(changeprobability))
    dftrimmed = @where(df, :changeprob .== changeprobability,
                            :sigma .== sigmavalue)
    ax.plot(dftrimmed[:param], dftrimmed[col], #color = color,
            label = string(sigmavalue))
    ylabel(ylabelstring); xlabel(xlabelstring); title(titlestring)
    ax.grid(true); ax.legend()
    ax.set_ylim(ylims[1], ylims[2]); ax.set_xlim(xlims[1], xlims[2])
end
function plot_fitness_multi_stoch(df::DataFrame, stochvalue::Float64
                    ;savepath = homedir(),
                    col = :meanterrors,
                    xlabelstring = "Leaky - Multi",
                    ylabelstring = "Mean error",
                    ylims = [0., 2.],
                    xlims = [0., 1.],
                    titlestring = "Stoch: "*string(stochvalue),
                    filestr = "t"
                    )
    fig = figure(); ax = gca();
    changeprobabilityvalues = unique(df[:changeprobability])
    clist = [cmap(i) for i in range(0, 1, length = length(changeprobabilityvalues))]
    for i in length(changeprobabilityvalues):-1:1
        plot_fitness_core_multi_stoch(df, changeprobabilityvalues[i], stochvalue, ax,
                        col=col, color = clist[i],
                        xlabelstring = xlabelstring,
                        ylabelstring = ylabelstring,
                        ylims = ylims, xlims = xlims,
                        titlestring = titlestring)
    end
    PyPlot.savefig(joinpath(savepath,
        "gridsearchMulti_stoch"* string(stochvalue)* "_" * xlabelstring *
        "_" * filestr *".png"))
end
function plot_fitness_core_multi_stoch(df::DataFrame, changeprobability::Float64,
                            stochvalue::Float64,
                            ax; col = :meanterrors,
                            color = "b",
                            xlabelstring = "",
                            ylabelstring = "Mean error",
                            ylims = [0., 2.], #[0., 10^6],
                            xlims = [-0.03, 1.],
                            titlestring = "Pc: "*string(changeprobability))

    dftrimmed = @where(df, :changeprobability .== changeprobability,
                            :stochasticity .== stochvalue)
    ax.plot(dftrimmed[:param], dftrimmed[col], #color = color,
            label = string(changeprobability))
    ylabel(ylabelstring); xlabel(xlabelstring); title(titlestring)
    ax.grid(true); ax.legend()
    ax.set_ylim(ylims[1], ylims[2]); ax.set_xlim(xlims[1], xlims[2])
end
function plot_fitness_gaussian_sigma(df::DataFrame, sigmavalue::Float64
                    ;savepath = homedir(),
                    col = :meanmerrors,
                    xlabelstring = "Leaky - Multi",
                    ylabelstring = "Mean error",
                    ylims = [0., 10.],
                    xlims = [-0.03, 1.],
                    titlestring = "Sigma: "*string(sigmavalue),
                    filestr = "g"
                    )
    fig = figure(); ax = gca();
    changeprobabilityvalues = unique(df[:changeprobability])
    clist = [cmap(i) for i in range(0, 1, length = length(changeprobabilityvalues))]
    for i in length(changeprobabilityvalues):-1:1
        plot_fitness_core_gaussian_sigma(df, changeprobabilityvalues[i], sigmavalue, ax,
                        col=col, color = clist[i],
                        xlabelstring = xlabelstring,
                        ylabelstring = ylabelstring,
                        ylims = ylims, xlims = xlims,
                        titlestring = titlestring)
    end
    PyPlot.savefig(joinpath(savepath,
        "gridsearchGauss_s"* string(sigmavalue)* "_" * xlabelstring *
        "_" * filestr *".png"))
end
function plot_fitness_core_gaussian_sigma(df::DataFrame, changeprobability::Float64,
                            sigmavalue::Float64,
                            ax; col = :meanmerrors,
                            color = "b",
                            xlabelstring = "",
                            ylabelstring = "Mean error",
                            ylims = [0., 10.],
                            xlims = [-0.03, 1.],
                            titlestring = "Pc: "*string(changeprobability))

    dftrimmed = @where(df, :changeprobability .== changeprobability,
                            :sigma .== sigmavalue)
    ax.plot(dftrimmed[:param], dftrimmed[col], #color = color,
            label = string(changeprobability))
    ylabel(ylabelstring); xlabel(xlabelstring); title(titlestring)
    ax.grid(true); ax.legend()
    ax.set_ylim(ylims[1], ylims[2]); ax.set_xlim(xlims[1], xlims[2])
end


text:
using MacroTools: prewalk, postwalk
using MacroTools
import MacroTools

constants = [:π, :ℯ, :Inf]

# ----- #
# types #
# ----- #

struct stringifyError <: Exception
    ex::Expr
    msg::String
end

function stringifyError(ex::Expr)
    msg = """Dolang does not know how to stringify
    \t$(ex)
    """
    stringifyError(ex, msg)
end

function Base.showerror(io::IO, ne::stringifyError)
    print(io, ne.msg)
end

struct UnknownFunctionError <: Exception
    func_name::Symbol
    msg::String
end

# --------- #
# stringify #
# --------- #

"""
    stringify(var::Union{String,Symbol}, n::Integer)

stringify the string or symbol in the following way:

- if `n >= 0` return `_var__n_`
- if `n < 0` return `_var_mn_`

"""
function stringify(var::Symbol, n::Integer)
    Symbol(var, "_", n >= 0 ? "_" : "m", abs(n), "_")
end

"""
    stringify(x::Tuple{Symbol,Integer})

Same as `stringify(x[1], x[2])`
"""
stringify(x::Tuple{Symbol,T}) where {T<:Integer} = stringify(x[1], x[2])

"""
    stringify(x::Symbol)

stringify the symbol by returning `_x_` if `x` doesn't alread have leading
and trailling `_` characters
"""
function stringify(x::Symbol)
    if x in constants
        return x
    end
    str_x = string(x)
    if str_x[end] == "_"
        return x
    end
    return Symbol(string(x, "_"))
end

# ---------- #
# time_shift #
# ---------- #

operators = [:+, :*, :^, :-, :/]

create_variable(sym, date::Int64) = date>0 ? :($sym[t+$date]) : date<0 ? :($sym[t-$(-date)]) : :($sym[t])

### sanitize expressions

# variable match
function match_var(expr)
    res = @match expr begin
        v_Symbol[t+d_Int64] => (v,d)
        v_Symbol[t-d_Int64] => (v,-d)
        v_Symbol[t] => (v,0)
        _ => return nothing
    end
    return res
end

# variable/symbol match
function match_var_flexible(expr)
    res = @match expr begin
        v_Symbol(d_Int64) => (v,d)
        v_Symbol[t+d_Int64] => (v,d)
        v_Symbol[t-d_Int64] => (v,-d)
        v_Symbol[t] => (v,0)
    end
    return res
end

function match_varfun(expr)
    res = @match expr begin
        v_Symbol(d_Int64) => (v,d)
        fun_Symbol(args__) => (fun,args)
        _ => return nothing
    end
    return res
end

function sanitize(expr::Expr; variables::Vector{Symbol}=Symbol[])
    m = match_var(expr)
    if m !== nothing
        return create_variable(m...)
    end
    m = match_varfun(expr)
    if (m !== nothing)
        if (typeof(m) <: Tuple{Symbol, Int64}) & (m[1] in variables)
            return create_variable(m[1], m[2])
        else
            if m[1] in variables
                error("Incorrect subscript for variable ", m[1], " : ", m[2])
            # else
            #     return expr
            end
        end
    end
    head = expr.head
    args = expr.args
    new_args = []
    for a in args
        if a in variables
            push!(new_args, create_variable(a,0))
        else
            dig = sanitize(a; variables=variables)
            push!(new_args, dig)
        end
    end
    return Expr(head, new_args...)
end

sanitize(s::Symbol; variables::Vector{Symbol}=Symbol[]) = s in variables ? create_variable(s,0) : s
sanitize(s; variables::Vector{Symbol}=Symbol[]) = s

compat_string(s::AbstractString) = replace(replace(replace(s, "|"=>"⟂"), "**"=>"^"), "inf"=>"Inf")

function parse_string(s::AbstractString; variables=Symbol[]) :: Union{Expr, Symbol}
    expr = sanitize(Meta.parse(compat_string(s)), variables=variables)
    # if typeof(expr) <: Symbol
    #     return :(+$expr)
    # else
        return expr
    # end
end

# function parse_equation(expr)
#     m = @match expr begin
#         # (lhs_ = rhs_  | lb_ <= x_ <= ub)_     => (lhs, rhs, lb, x, ub)
#         # (lhs_ == rhs_ | lb_ <= x_ <= ub)_     => (lhs, rhs, lb, x, ub)
#         # (lhs_ == rhs_ | lb_ <= x_ <= ub)_     => (lhs, rhs, lb, x, ub)
#         (lhs_ = rhs)                         => (lhs, rhs)
#         (lhs_ == rhs)                        => (lhs, rhs)
#         # (lhs_)                                => (lhs,)
#     end
#     return m
# end

function match_equation(expr)
    m = @match expr begin
        # lhs_ == rhs_                      => (lhs, rhs)
        lhs_ ⟂ lb_ <= x_ <= ub_     => (lhs, lb, x, ub)
        (lhs_ = rhs_)                       => (lhs, rhs)
        (lhs_ == rhs_)                       => (lhs, rhs)
        lhs_                                => (lhs,)
    end
    return m
end

function match_equality(expr)
    m = @match expr begin
        (lhs_ = rhs_)                       => (lhs, rhs)
        (lhs_ == rhs_)                       => (lhs, rhs)
        lhs_                                => (lhs,)
    end
    return m
end

###
###
###

# variable/symbol match
function match_expr(expr)
    res = @match expr begin
        # v_Symbol(d_Int64) => (v,d)
        v_Symbol[t+d_Int64] => (v,d)
        v_Symbol[t-d_Int64] => (v,-d)
        v_Symbol[t] => (v,0)
        v_Symbol => (v in operators ? nothing : v)
    end
    return res
end

###
### Basic symbolic operations
###

struct SymbolList
    variables:: Vector{Tuple{Symbol, Int64}}
    parameters:: Vector{Symbol}
    functions:: Vector{Symbol}
end

function list_symbols(expr)

    funs = Symbol[]
    prewalk(u->(x=match_varfun(u);if x != nothing push!(funs, x[1]) end; u), expr)
    l = []
    prewalk(u->(x=match_expr(u);if x != nothing push!(l, x) else u end), expr)

    vars = [e for e in l if typeof(e) <: Tuple{Symbol, Int64}]

    parms = [e for e in l if (typeof(e) <: Symbol) & !(e in funs) & !(e in constants)]

    sl = SymbolList(unique(vars), unique(parms), unique([f for f in funs if !(f in operators)]))

    return sl
end

list_variables(expr) = list_symbols(expr).variables
list_parameters(expr) = list_symbols(expr).parameters

stringify(s::AbstractString) = stringify(parse_string(compat_string(s)))

function stringify(tree, parameters)
    function fun(u)
        a = match_expr(u)
        if a===nothing
            return u
        end
        if typeof(a) <: Symbol
            if a in parameters
                return stringify(a)
            else
                return a
            end
        end
        return stringify(a)
    end
    prewalk(fun, tree)
end

stringify(tree) = stringify(tree, list_parameters(tree))

time_shift(tree, Δt) = prewalk(u->(a=match_var(u); a === nothing ? u : create_variable(a[1],a[2]+Δt)),tree)
steady_state(tree) = prewalk(u->(a=match_var(u); a === nothing ? u : create_variable(a[1],0)),tree)

# ---- #
# subs #
# ---- #

function subs(expr, d::AbstractDict{Tuple{Symbol, Int64}, T}) where T
    # replace variables
    fun  = u -> (a = match_var(u); a === nothing ? u : (a in keys(d)) ? d[a] : u)
    postwalk(fun, expr)
end

subs(expr, pair::Pair{Tuple{Symbol,Int64}, Any}) = subs(expr, Dict(pair))

subs(expr, from, to) =  subs(expr, from=>to)

# function subs(expr, from, to)
#     fun  = u -> (a = match_var(u); a === nothing ? u : (a==from) ? to : u)
#     postwalk(fun, expr)
# end

# function subs(expr, d::AbstractDict{Tuple{Symbol, Int64}, T}) where T
#     # replace variables
#     res = expr
#     for (from, to) in d
#         res = subs(res, from, to)
#     end
# end

function tsubs(expr, d::AbstractDict{Tuple{Symbol, Int64}, T}) where T
    # replace symbol, with lags
    if Set([h[2] for h in keys(d)]) != Set([0])
        error("All variables to be replaced must be given at date t.")
    end
    fun  = u -> (a = match_var(u); a === nothing ? u : ((a[1],0) in keys(d)) ? time_shift(d[(a[1],0)], a[2]) : u)
    postwalk(fun, expr)
end

tsubs(expr, pair::Pair{Tuple{Symbol,Int64}, Any}) = tsubs(expr, from=>to)
tsubs(expr, from, to) = tsubs(expr, from=>to)

###########
#         #
###########

function csubs(expr, d::AbstractDict)
    dd = reorder_triangular_block(d)
    subs(expr, dd)
end

# ---------#
# arg_name #
# ---------#

function arg_name(s::Symbol)::Symbol
    # is the symbol already stringified?
    str_s = string(s)
    if str_s[1] == str_s[end] == '_'
        parts = match(r"^_(.+)_[m_]\d+_$", str_s)
        if parts === nothing
            return s
        else
            return Symbol(parts[1])
        end
    else
        return s
    end
end
arg_name(s::Tuple{Symbol,Int}) = s[1]
function arg_name(ex::Expr)::Symbol
    if is_time_shift(ex)
        return ex.args[1]
    else
        m = string(
            "Expression $(ex) does not look like a variable. ",
            "The only allowed expressions must satisfy `is_time_shift(ex)`"
        )
        error(m)
    end
end

function arg_time(s::Symbol)::Int
    # is the symbol already stringified?
    str_s = string(s)
    if str_s[1] == str_s[end] == '_'
        # yep, we are already stringified, need to extract the
        r = r"(m)?(\d+)_$"
        parts = match(r, str_s)
        if parts === nothing
            # not a match
            return 0
        else
            shift = Base.parse(Int, parts[2])
            return parts[1] === nothing ? shift : -shift
        end
    else
        return 0
    end
end
arg_time(s::Tuple{Symbol,Int}) = s[2]
function arg_time(ex::Expr)::Int
    if is_time_shift(ex)
        return ex.args[2]
    else
        m = string(
            "Expression $(ex) does not look like a variable. ",
            "The only allowed expressions must satisfy `is_time_shift(ex)`"
        )
        error(m)
    end
end

arg_name_time(s) = (arg_name(s), arg_time(s))

arg_names(s::AbstractVector) = arg_name.(s)
arg_names(s::AbstractDict) = vcat([arg_names(v) for v in values(s)]...)

# --------- #
# Utilities #
# --------- #

"""
```julia
is_time_shift(ex::Expr)
```

Determine if the expression has the form `var(n::Integer)`.
"""
is_time_shift(ex::Expr) = ex.head == :call &&
                          length(ex.args) == 2 &&
                          isa(ex.args[1], Symbol) &&
                          isa(ex.args[2], Integer)

max_issues_repo_name: EconForge/Dolang


text:
@symbol_func function K_n(cur_reactor::AbstractReactor)
  cur_nu_n = cur_reactor.nu_n

  cur_K = 2.0

  cur_K *= cur_reactor.N_G

  cur_K /= sqrt(cur_reactor.pi) ^ 3

  cur_K /= cur_reactor.epsilon ^ 2

  cur_K *= gamma( cur_nu_n + 3/2 )

  cur_K /= gamma( cur_nu_n + 2.0 )

  cur_K
end


text:
export
	AbstractRobotType,
	DeliveryBot,
	DefaultRobotType

abstract type AbstractRobotType end
struct DeliveryBot <: AbstractRobotType end
const DefaultRobotType = DeliveryBot

export
	AbstractID,
	TemplatedID,
    ObjectID,
    BotID,
    RobotID,
    LocationID,
    ActionID,
    OperationID,
	AgentID,
	VtxID

abstract type AbstractID end

"""
    struct TemplatedID{T} <: AbstractID

A simple way to dispatch by Node type.
"""
struct TemplatedID{T} <: AbstractID
    id::Int
end

@with_kw struct ObjectID <: AbstractID
	id::Int = -1
end
@with_kw struct BotID{R<:AbstractRobotType} <: AbstractID
	id::Int = -1
end
const RobotID = BotID{DeliveryBot}
@with_kw struct LocationID <: AbstractID
	id::Int = -1
end
@with_kw struct ActionID <: AbstractID
	id::Int = -1
end
@with_kw struct OperationID <: AbstractID
	id::Int = -1
end
"""
	AgentID
Special helper for identifying agents.
"""
@with_kw struct AgentID <: AbstractID
	id::Int = -1
end
"""
	VtxID
Special helper for identifying schedule vertices.
"""
@with_kw struct VtxID <: AbstractID
	id::Int = -1
end

Base.summary(id::A) where {A<:AbstractID} = string(string(A),"(",get_id(id),")")
Base.string(id::A) where {A<:AbstractID} = summary(id)

export
	get_unique_id,
	get_unique_invalid_id,
	reset_id_counter!,
	reset_all_id_counters!,
	reset_invalid_id_counter!,
	reset_all_invalid_id_counters!,
    reset_operation_id_counter!,
    get_unique_operation_id,
    reset_action_id_counter!,
    get_unique_action_id

VALID_ID_COUNTERS = Dict{DataType,Int}()
INVALID_ID_COUNTERS = Dict{DataType,Int}()

function get_unique_id(::Type{T}) where {T}
	global VALID_ID_COUNTERS
	id = get!(VALID_ID_COUNTERS,T,1)
	VALID_ID_COUNTERS[T] += 1
	return T(id)
end
function reset_id_counter!(::Type{T}) where {T}
	global VALID_ID_COUNTERS
	VALID_ID_COUNTERS[T] = 1
end
function reset_all_id_counters!()
	global VALID_ID_COUNTERS
	for k in collect(keys(VALID_ID_COUNTERS))
		reset_id_counter!(k)
	end
end
function get_unique_invalid_id(::Type{T}) where {T}
	global INVALID_ID_COUNTERS
	id = get!(INVALID_ID_COUNTERS,T,-1)
	INVALID_ID_COUNTERS[T] -=1
	return T(id)
end
function reset_invalid_id_counter!(::Type{T}) where {T}
	global INVALID_ID_COUNTERS
	INVALID_ID_COUNTERS[T] = -1
end
function reset_all_invalid_id_counters!()
	global INVALID_ID_COUNTERS
	for k in collect(keys(INVALID_ID_COUNTERS))
		reset_invalid_id_counter!(k)
	end
end
get_unique_operation_id() = get_unique_id(OperationID)
get_unique_action_id() = get_unique_id(ActionID)
reset_action_id_counter!() = reset_id_counter!(ActionID)
reset_operation_id_counter!() = reset_id_counter!(OperationID)

# OPERATION_ID_COUNTER = 0
# get_unique_operation_id() = OperationID(Int(global OPERATION_ID_COUNTER += 1))
# function reset_operation_id_counter!()
#     global OPERATION_ID_COUNTER = 0
# end
# ACTION_ID_COUNTER = 0
# get_unique_action_id() = ActionID(Int(global ACTION_ID_COUNTER += 1))
# function reset_action_id_counter!()
#     global ACTION_ID_COUNTER = 0
# end

export
	get_id,
	valid_id

get_id(id::AbstractID) = id.id
get_id(id::Int) = id
Base.:+(id::A,i::Int) where {A<:AbstractID} = A(get_id(id)+i)
Base.:+(id::A,i::A) where {A<:AbstractID} = A(get_id(id)+get_id(i))
Base.:-(id::A,i::Int) where {A<:AbstractID} = A(get_id(id)-i)
Base.:-(id::A,i::A) where {A<:AbstractID} = A(get_id(id)-get_id(i))
Base.:(<)(id1::AbstractID,id2::AbstractID) = get_id(id1) < get_id(id2)
Base.:(>)(id1::AbstractID,id2::AbstractID) = get_id(id1) > get_id(id2)
Base.isless(id1::AbstractID,id2::AbstractID) = id1 < id2
Base.convert(::Type{ID},i::Int) where {ID<:AbstractID} = ID(i)
Base.copy(id::ID) where {ID<:AbstractID} = ID(get_id(id))

valid_id(id::AbstractID) = get_id(id) > -1


text:
push!(LOAD_PATH, "..")

import Pkg
Pkg.add("Documenter")

using Documenter, Kinetic
using Kinetic: KitBase, KitML

tutorial_page = [
    "Examples" => "tutorial.md",
    "Advection diffusion" => "eg_advection.md",
    "Burgers" => "eg_burgers.md",
    "Shock tube" => "eg_shock.md",
    "Lid-driven cavity" => "eg_cavity.md",
]

type_page = [
    "Configuration" => "type_config.md",
    "Setup" => "type_setup.md",
    "Domain" => "type_domain.md",
    "Velocity" => "type_velocity.md",
    "Property" => "type_property.md",
    "Condition" => "type_ib.md",
    "FVM" => "type_fvm.md",
]

solver_page = [
    "Framework" => "solver.md",
    "Preprocess" => "solver_pre.md",
    "Timestep" => "solver_timestep.md",
    "Reconstruction" => "solver_reconstruction.md",
    "Flux" => "solver_flux.md",
    "Update" => "solver_update.md",
    "Postprocess" => "solver_post.md",
]

utility_page = [
    "I / O" => "api_io.md",
    "Math" => "api_math.md",
    "Theory" => "api_theory.md",
    "Physical space" => "api_geo.md",
    "Phase space" => "api_phase.md",
    "Configuration" => "api_config.md",
    "Stepper" => "api_step.md",
]

parallel_page = [
    "General" => "parallel.md",
    "Threading" => "para_thread.md",
    "Distributed" => "para_dist.md",
    "CUDA" => "para_cuda.md",
]

ml_page = [
    "KitML" => "kitml1.md",
    "UBE" => "kitml2.md",
]

fortran_page = [
    "KitFort" => "fortran1.md",
    "Benchmark" => "fortran2.md",
]

format = Documenter.HTML(
    assets = ["assets/favicon.ico"],
    collapselevel = 1,
)

makedocs(
    sitename= "Kinetic.jl",
    modules = [Kinetic, KitBase, KitML],
    pages = Any[
        "Home" => "index.md",
        "Installation" => "install.md",
        "Physics" => "physics.md",
        "Type" => type_page,
        "Solver" => solver_page,
        "Tutorial" => tutorial_page,
        "Parallelization" => parallel_page,
        "Utility" => utility_page,
        "SciML" => ml_page,
        "Fortran" => fortran_page,
        "Index" => "function_index.md",
        "Python" => "python.md",
        "Contribution" => "contribution.md",
        "Reference" => "reference.md",
    ],
    format = format,
)

deploydocs(
    repo = "github.com/vavrines/Kinetic.jl.git",
)


text:
# ---
# title: 290. Word Pattern
# id: problem290
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/word-pattern/description/>
# hidden: true
# ---
#
# Given a `pattern` and a string `s`, find if `s` follows the same pattern.
#
# Here **follow** means a full match, such that there is a bijection between a
# letter in `pattern` and a **non-empty** word in `s`.
#
#
#
# **Example 1:**
#
#
#
#     Input: pattern = "abba", s = "dog cat cat dog"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: pattern = "abba", s = "dog cat cat fish"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: pattern = "aaaa", s = "dog cat cat dog"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: pattern = "abba", s = "dog dog dog dog"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= pattern.length <= 300`
#   * `pattern` contains only lower-case English letters.
#   * `1 <= s.length <= 3000`
#   * `s` contains only lower-case English letters and spaces `' '`.
#   * `s` **does not contain** any leading or trailing spaces.
#   * All the words in `s` are separated by a **single space**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function helio_rv(jd::T, t::T, P::T, V0::T, K::T, ecc::T, ω::T) where {T<:AbstractFloat}
    E = kepler_solver(2 * T(pi) * (jd - t) / P, ecc)
    ν = trueanom(E, ecc)
    ω = deg2rad(ω)
    return K*(cos(ν + ω) + (ecc*cos(ω))) + V0
end

"""
    helio_rv(jd, T, P, V_0, K[, e, ω]) -> rv

### Purpose ###

Return the heliocentric radial velocity of a spectroscopic binary.

### Explanation ###

This function will return the heliocentric radial velocity of a spectroscopic
binary star at a given heliocentric date given its orbit.

### Arguments ###

* `jd`: time of observation, as number of Julian days.
* `T`: time of periastron passage (max. +ve velocity for circular orbits), same
  time system as `jd`
* `P`: the orbital period in same units as `jd`
* `V_0`: systemic velocity
* `K`: velocity semi-amplitude in the same units as `V_0`
* `e`: eccentricity of the orbit.  It defaults to 0 if omitted
* `ω`: longitude of periastron in degrees.  It defaults to 0 if omitted

### Output ###

The predicted heliocentric radial velocity in the same units as Gamma for the
date(s) specified by `jd`.

### Example ###

(1) What was the heliocentric radial velocity of the primary component of HU Tau
at 1730 UT 25 Oct 1994?

```jldoctest
julia> using AstroLib

julia> jd = juldate(94, 10, 25, 17, 30); # Obtain Geocentric Julian days

julia> hjd = helio_jd(jd, ten(04, 38, 16) * 15, ten(20, 41, 05)); # Convert to HJD

julia> helio_rv(hjd, 46487.5303, 2.0563056, -6, 59.3)
-62.965570107789475
```

NB: the functions `juldate` and `helio_jd` return a reduced HJD (HJD - 2400000)
and so T and P must be specified in the same fashion.

(2) Plot two cycles of an eccentric orbit, \$e=0.6\$, \$\\omega=45\\degree\$ for
both components of a binary star.  Use
[PyPlot.jl](https://github.com/JuliaPy/PyPlot.jl) for plotting.

```julia
using PyPlot
φ = range(0, stop=2, length=1000); # Generate 1000 phase points
plot(φ ,helio_rv.(φ, 0, 1, 0, 100, 0.6, 45)) # Plot 1st component
plot(φ ,helio_rv.(φ, 0, 1, 0, 100, 0.6, 45+180)) # Plot 2nd component
```

### Notes ###

The user should ensure consistency with all time systems being used (i.e. `jd`
and `t` should be in the same units and time system).  Generally, users should
reduce large time values by subtracting a large constant offset, which may
improve numerical accuracy.

If using the the function `juldate` and `helio_jd`, the reduced HJD time system
must be used throughtout.

Code of this function is based on IDL Astronomy User's Library.
"""
helio_rv(jd::Real, t::Real, P::Real, V0::Real,
         K::Real, ecc::Real=0, ω::Real=0) =
             helio_rv(promote(float(jd), float(t), float(P),
                              float(V0), float(K), float(ecc), float(ω))...)


text:
function a_computer_vision_toolbox()
    mat"""$i = license('test','Video_and_Image_Blockset')"""
    @assert i == 1 "the Matlab Computer Vision System Toolbox is not available"
end
a_csvfile(file) = @assert isfile(file) "missing csv file"
function a_table(t)
    @assert length(t) ≠ 0 "csv file has no lines"
    for x in keys(csvfile_columns)
        @assert x ∈ propertynames(t) "missing column $x from csv file"
    end
    for x in propertynames(t)
        @assert haskey(csvfile_columns, x) "unkown column, $x, in csv file"
    end
end

c_resfile(io, ::Missing) = (println(io, "- res file missing in csv file"); return missing)
function c_resfile(io, resfile)
    if isfile(resfile)
        resfile
    else
        println(io, "- res file does not exist")
        missing
    end
end

c_poi_names(io, ::Missing) = (println(io, "- POI names missing in csv file"); return missing)
function c_poi_names(io, poi_names)
    good = true
    if !all(x -> !occursin(' ', String(x)), poi_names)
        println(io, "- POI name/s contains space/s")
        good = false
    end
    if !allunique(poi_names)
        println(io, "- POI names must be unique")
        good = false
    end
    if all(!occursin("track", String(name)) for name in poi_names)
        println(io, "- a track POI is missing")
        good = false
    end
    good ? poi_names : missing
end

c_expected_locations(io, ::Missing) = nothing
function c_expected_locations(io, expected_locations)
    if !allunique(values(expected_locations))
        println(io, "- at least two of the POI's expected locations are identical")
        return missing
    end
    expected_locations
end

a_expected_locations(io, ::Missing, poi_names, npoints) = nothing
function a_expected_locations(io, expected_locations, poi_names, npoints)
    for (k, v) in expected_locations
        i = findfirst(==(k), poi_names)
        n = npoints[i]
        n == 1 || println(io, "- you expected $k to be at $v, yet there are $n lines in its column (column #$i) in the res-file (instead of just one)")
    end
end

a_singular_pois(io, poi_names, npoints) = map(poi_names, npoints) do name, n
    if n ≠ 1 && name ∈ (:feeder, :pickup, :dropoff, :nest, :north, :south)
        println(io, "- $name is expected to have exactly one row (point) in its res file, yet it has $n rows")
    end
end

a_coords(io, x, y, d, e) = nothing
function a_coords(io, resfile::SystemPath, poi_names::Vector{Symbol}, duration::Float64, expected_locations)
    matopen(string(resfile)) do mio
        for field in ("xdata", "ydata", "status")
            if !MAT.exists(mio, field)
                println(io, "- resfile missing $field")
                return nothing
            end
        end
        haskey(read(mio, "status"), "FrameRate") || println(io,  "- resfile missing frame rate")
        xdata = read(mio, "xdata")
        n = size(xdata, 2)
        npoints = [length(nzrange(xdata, j)) for j in 1:n]
        a_expected_locations(io, expected_locations, poi_names, npoints)
        a_singular_pois(io, poi_names, npoints)
        !all(iszero, npoints) || println(io, "- res file is empty")
        nmax = maximum(npoints)
        nmax > 5 || println(io, "- res file missing a POI with more than 5 data points (e.g. a track)")
        npois = length(poi_names)
        npois == n || println(io, "- number of POIs, $npois, doesn't match the number of res columns, $n")
        rows = rowvals(xdata)
        fr = read(mio, "status")["FrameRate"]
        haskey(read(mio, "status"), "nFrames") || println(io,  "- resfile missing number of total frames")
        nf = read(mio, "status")["nFrames"]
        duration2 = nf/fr
        isapprox(duration2, duration, atol = 1) || println(io, "- the duration of the POI video file, $duration s, and the one reported in the res file, $duration2 s, are not the same")
        if npois > n
            println(io, "- there are more POIs ($npois) than there are columns in the res file ($n)")
            return nothing
        end
        for (j, name) in enumerate(poi_names)
            i = nzrange(xdata, j)
            if !isempty(i)
                t = rows[i[end]]/fr
                t ≤ duration || println(io, "- the time-stamp of the $name POI is not in the video")
            end
        end
    end
end

c_videofile(io, ::Missing, what) = (println(io, "- ", what, " video file missing in csv file"); return missing)
function c_videofile(io, videofile, what)
    if !isfile(videofile)
        println(io, "- ", what, " video file does not exist")
        return missing
    else
        return videofile
    end
end

a_turning_point(io, x, y, d) = nothing
function a_turning_point(io, poi_videofile::AbstractString, t::Float64, duration::Float64)
    0 ≤ t ≤ duration || println(io, "- turning point time-stamp is not in the video")
end

a_checker_size(io, ::Missing) = println(io, "- checker-size is missing in the csv file")
a_checker_size(io, checker_size) = checker_size > 0 || println(io, "- checkers must be larger than zero")

a_extrinsic(io, ::Missing, duration) = println(io, "- extrinsic time-stamp is missing in the csv file")
a_extrinsic(io, extrinsic, duration) = 0 ≤ extrinsic ≤ duration || println(io, "- extrinsic time-stamp is not in the video")

function a_intrinsic(io, intrinsic, duration)
    all(x -> 0 ≤ x ≤ duration, intrinsic) || println(io, "- intrinsic time-stamp is not in the video")
    issorted(intrinsic) || println(io, "- the intrinsic starting time cannot come after the ending time")
end
a_intrinsic(io, ::Missing, _) = nothing

a_calibration(io, ::Missing, extrinsic, intrinsic, checker_size) = nothing
function a_calibration(io, calib_videofile, extrinsic, intrinsic, checker_size)
    a_checker_size(io, checker_size)
    duration = get_duration(calib_videofile)
    a_extrinsic(io, extrinsic, duration)
    a_intrinsic(io, intrinsic, duration)
end

c_nest2feeder(io, ::Missing) = missing
function c_nest2feeder(io, nest2feeder)
    if nest2feeder < 0
        println(io, "- nest to feeder distance must be larger than zero")
        return missing
    else
        return nest2feeder
    end
end
a_azimuth(io, ::Missing) = nothing
a_azimuth(io, azimuth) = 0 < azimuth < 360 || println(io, "- azimuth must be between 0° and 360°")

__getfeeder(x) = get(x, :initialfeeder, get(x, :pickup, get(x, :feeder, missing)))
_get_expected_nest2feeder(::Missing) = missing
_get_expected_nest2feeder(x) = _get_expected_nest2feeder(get(x, :nest, missing), __getfeeder(x))
_get_expected_nest2feeder(nest, feeder) = norm(nest - feeder)
_get_expected_nest2feeder(::Missing, feeder) = missing
_get_expected_nest2feeder(nest, ::Missing) = missing
_get_expected_nest2feeder(::Missing, ::Missing) = missing
a_expected_nest2feeder(io, ::Missing, nest2feeder::Missing) = nothing
a_expected_nest2feeder(io, ::Real, nest2feeder::Missing) = println(io, "- it seems like you have expectations on the distance between the nest and feeder, but nest2feeder is missing")
a_expected_nest2feeder(io, ::Missing, nest2feeder::Real) = nothing#println(io, "- it seems like you should have expectations on the distance between the nest and feeder since nest2feeder is not missing")
a_expected_nest2feeder(io, expected_nest2feeder::Real, nest2feeder::Real) = expected_nest2feeder ≈ nest2feeder || println(io, "- nest2feeder, $nest2feeder, is not equal to the distance between the expected nest and feeder locations, $expected_nest2feeder")

a_extra_correction(io, ::Bool) = nothing
a_extra_correction(io, _) = println(io, "- extra correction is missing/wrong (use true/false)")

function check4errors(x)
    io = IOBuffer()
    resfile = c_resfile(io, x.resfile)
    poi_names = c_poi_names(io, x.poi_names)
    expected_locations = c_expected_locations(io, x.expected_locations)
    poi_videofile = c_videofile(io, x.poi_videofile, "POI")
    duration = get_duration(poi_videofile)
    a_coords(io, resfile, poi_names, duration, expected_locations)
    a_turning_point(io, poi_videofile, x.turning_point, duration)
    calib_videofile = c_videofile(io, x.calib_videofile, "Calibration")
    a_calibration(io, calib_videofile, x.extrinsic, x.intrinsic, x.checker_size)
    nest2feeder = c_nest2feeder(io, x.nest2feeder)
    a_azimuth(io, x.azimuth)
    expected_nest2feeder = _get_expected_nest2feeder(expected_locations)
    a_expected_nest2feeder(io, expected_nest2feeder, nest2feeder)
    a_extra_correction(io, x.extra_correction)
    String(take!(io))
end

get_duration(::Missing) = missing
function get_duration(file)
    p = FFMPEG.exe(`-i $file -show_entries format=duration -v quiet -of csv="p=0"`, command=FFMPEG.ffprobe, collect=true)
    parse(Float64, only(p))
end

# function get_dimentions(file)
#     p = FFMPEG.exe(`-i $file -show_entries stream=width,height -of csv=p=0:s=x`, command=FFMPEG.ffprobe, collect=true)
#     txt = split(only(p), 'x')
#     parse.(Int, txt)
# end


text:
function build_op_model!(op_model::OperationModel{M}; kwargs...) where M<:AbstractOperationModel
    sys = get_system(op_model)
    _build_canonical!(op_model.canonical, op_model.model_ref, sys; kwargs...)
    return
end

function _build_canonical!(canonical::CanonicalModel, ref::ModelReference, sys::PSY.System; kwargs...)

    verbose = get(kwargs, :verbose, true)
    transmission = ref.transmission

    # Build Injection devices
    for (_, device_model) in ref.devices
        verbose && @info "Building $(device_model.device) with $(device_model.formulation) formulation"
        construct_device!(canonical, sys, device_model, transmission; kwargs...)
    end

    # Build Network
    verbose && @info "Building $(transmission) network formulation"
    construct_network!(canonical, sys, transmission; kwargs...)

    # Build Branches
    for (_, branch_model) in ref.branches
        verbose && @info "Building $(branch_model.device) with $(branch_model.formulation) formulation"
        construct_device!(canonical, sys, branch_model, transmission; kwargs...)
    end

    #Build Service
    for (_, service_model) in ref.services
        #construct_service!(canonical, sys, service_model, transmission; kwargs...)
    end

    # Objective Function
    verbose && @info "Building Objective"
    JuMP.@objective(canonical.JuMPmodel,
                    MOI.MIN_SENSE,
                    canonical.cost_function)

    return

end


text:
import .DataFrames

function _sanitize!(df::T, data::Vector{Data}, elements::_DataElements; kwargs...)::T where T <: DataFrames.AbstractDataFrame
    zero!(df)
    return df
end

function zero!(df::T; kwargs...)::T where T <: DataFrames.AbstractDataFrame
    for column in names(df)
        df[:, column] = zero!(df[:, column])
    end
    return df
end

function _elements!(all_elements::Vector{Any}, df::DataFrames.AbstractDataFrame; kwargs...) where T
    # push!(all_elements, df)
    for column in names(df)
        _elements!(all_elements, df[column])
    end
    return all_elements
end


text:
using Random
using JSON

ENV["PYTHONPATH"] = "$(homedir())/lib_acetb/dftbplus/build/install/lib/python3.7/site-packages/"
ENV["JULIA_PROJECT"] = "$(@__DIR__)"

const DFTB_PATH = "$(homedir())/lib_acetb/dftbplus/build/install/bin/"
const OUTPUT_DIR = joinpath(@__DIR__, "fits/")
const DATABCC = "data/BCC/"
const DATAFCC = "data/FCC/"
const PARAM_FILE = "Al-Al.acetb.json"

# always include lattice expansion data?
# FCC: /home/eng/essswb/ownCloud/SK_AITB/FCC/without_HS/LATTICE_EXTENSION/SUPERCELL-9x9x9/SK-supercell-noHS-000.h5
# BCC: /home/eng/essswb/ownCloud/SK_AITB/BCC/without_HS/LATTICE_EXPANSION/SK-supercell-noHS-000.h5

using ACEtb.FitTools
fc = FitContext(@__DIR__, DFTB_PATH, OUTPUT_DIR, PARAM_FILE)

params = Dict(
    :order => 3,
    :deg => 12,
    :rcut => 10.0,
    :renv => 3.0,
    :envdeg => 4,
    :pcut => 1,
    :method => :RegRRQR,
    :rtol => 1e-8
)

bcc_files = readlines(pipeline(`find -L $DATABCC -name SK-supercell-\*.h5`))
fcc_files = readlines(pipeline(`find -L $DATAFCC -name SK-supercell-\*.h5`))
bcc_data = reshape([ (file, index) for file in bcc_files, index in 1:729 ], :)
fcc_data = reshape([ (file, index) for file in fcc_files, index in 1:729 ], :)
all_data = [bcc_data; fcc_data]

@info "Found $(length(bcc_data)) FCC and $(length(fcc_data)) FCC atomic environments"

if !isfile("train_test.json")
    @info "Generating train/test split..."
    perm = randperm(length(all_data))
    open("train_test.json", "w") do json_file
        JSON.print(json_file, perm)
    end
else
    @info "Loading train/test split"
    perm = JSON.parsefile("train_test.json")
end
train = perm[1:length(perm)÷2]
test = perm[length(perm)÷2+1:end]

training_data = all_data[train]
test_data = all_data[test]

cparams = copy(params)

N_train = parse(Int, ARGS[1])
N_test = 100

@show N_train, N_test
model, param_hash, Bint_table, cutf_func = fit_model(fc, cparams,
                                                     training_data[1:N_train],
                                                     test_data[1:N_test],
                                                     "data/exact-onsite-BCC.json")


text:
using TupleTools

NotImplementedError(name::Symbol) = NotImplementedError(name, ())

function Base.show(io::IO, e::NotImplementedError)
    str = join(map(typeof, e.args), ", ::")
    str = "::" * str
    print(
        io,
        "NotImplementedError: $(e.name) is not implemented for (",
        str,
        "), ",
        "please implement this method for your custom type",
    )
end

function Base.show(io::IO, e::NotImplementedError{Tuple{}})
    print(io, "$(e.name) is not implemented.")
end

function Base.show(io::IO, e::LocationConflictError)
    print(io, "LocationConflictError: ", e.msg)
    # print(io, "locations of $(e.blk1) and $(e.blk2) is conflict.")
end

function show(io::IO, e::QubitMismatchError)
    print(io, e.msg)
end

export islocs_inbounds, islocs_conflict

_sort(x::Vector; by = identity) = sort(x, by = by)
_sort(x::Tuple; by = identity) = TupleTools.sort(x, by = by)

# NOTE: this method assumes its input is not empty, it gets rid of errors
nonempty_minimum(x::UnitRange) = x.start
nonempty_minimum(x::Integer) = x
nonempty_maximum(x::UnitRange) = x.stop
nonempty_maximum(x::Integer) = x

const AddressVector{T} = Vector{T} where {T<:Union{Integer,UnitRange}}
const AddressNTuple{N,T} = NTuple{N,T} where {T<:Union{Integer,UnitRange}}
const AddressList{T} = Union{AddressVector{T},AddressNTuple{N,T} where N}

"""
    islocs_inbounds(n, locs) -> Bool

Check if the input locations are inside given bounds `n`.
"""
function islocs_inbounds(n::Int, locs::AddressList)
    length(locs) == 0 && return true
    locs = _sort(locs, by = x -> nonempty_minimum(x))
    (minimum(first(locs)) > 0 && maximum(last(locs)) <= n) || return false
    return true
end

"""
    islocs_conflict(locs) -> Bool

Check if the input locations has conflicts.
"""
function islocs_conflict(locs::AddressList)
    locs = _sort(locs)
    for i = 1:length(locs)-1
        nonempty_minimum(locs[i+1]) > nonempty_maximum(locs[i]) || return true
    end
    return false
end

function process_msgs(msgs...; default = "")
    msg = isempty(msgs) ? default : msgs[1]
    if isa(msg, AbstractString)
        msg = msg # pass-through
    elseif !isempty(msgs) && (isa(msg, Expr) || isa(msg, Symbol))
        # message is an expression needing evaluating
        msg = :(Main.Base.string($(esc(msg))))
    elseif isdefined(Main, :Base) &&
           isdefined(Main.Base, :string) &&
           applicable(Main.Base.string, msg)
        msg = Main.Base.string(msg)
    else
        # string() might not be defined during bootstrap
        msg = :(Main.Base.string($(Expr(:quote, msg))))
    end
    return msg
end

# NOTE: we may use @assert in the future
#       these macro will help us keep original APIs
export @assert_locs_safe, @assert_locs_inbounds

"""
    @assert_locs_inbounds <number of total qudits> <locations list> [<msg>]

Assert if all the locations are inbounds.
"""
macro assert_locs_inbounds(n, locs, msgs...)
    msg = process_msgs(msgs...; default = "locations is out of bounds!")

    return quote
        islocs_inbounds($(esc(n)), $(esc(locs))) || error($msg)
        nothing
    end
end

"""
    @assert_locs_safe <number of total qudits> <locations list> [<msg>]

Assert if all the locations are:
    - inbounds.
    - do not have any conflict.
"""
macro assert_locs_safe(n, locs, msgs...)
    msg = process_msgs(msgs...; default = "locations conflict.")
    return quote
        @assert_locs_inbounds $(esc(n)) $(esc(locs))
        islocs_conflict($(esc(locs))) && throw(LocationConflictError($msg))
        nothing
    end
end


text:
#!/usr/bin/env julia6
using ReadGlobal, InplaceRealFFT, LESFilter

function adddx!(out,inp,kim)
    nx,ny,nz = size(inp)
    Threads.@threads for l = 1:nz
        for j = 1:ny
            @simd for i = 1:nx
                @inbounds out[i,j,l] = muladd(inp[i,j,l],kim[i], out[i,j,l])
            end
        end
    end
end

function adddy!(out,inp,kim)
    nx,ny,nz = size(inp)
    Threads.@threads for l = 1:nz
        for j = 1:ny
            @simd for i = 1:nx
                @inbounds out[i,j,l] = muladd(inp[i,j,l],kim[j], out[i,j,l])
            end
        end
    end
end

function adddz!(out,inp,kim)
    nx,ny,nz = size(inp)
    Threads.@threads for l = 1:nz
        for j = 1:ny
            @simd for i = 1:nx
                @inbounds out[i,j,l] = muladd(inp[i,j,l],kim[l], out[i,j,l])
            end
        end
    end
end

function div(t11::AbstractString,t12::AbstractString,t13::AbstractString,t22::AbstractString,t23::AbstractString,t33::AbstractString,f1::AbstractString,f2::AbstractString,f3::AbstractString)
    nx,ny,nz,lx,ly,lz = getdimsize()
    s = (nx,ny,nz)
    kxim = (2π .* rfftfreq(nx,lx) .* im)
    kyim = (2π .* fftfreq(ny,ly) .* im)
    kzim = (2π .* fftfreq(nz,lz) .* im)

    FFTW.set_num_threads(Threads.nthreads())
    isfile("fftw_wisdom") && FFTW.import_wisdom("fftw_wisdom")

    dtype, padded = checkinput(t11,nx,ny,nz)

    aux1 = PaddedArray(t11,s,padded)
    aux2 = similar(aux1)

    begin #f1
        rfft!(aux1)
        adddx!(aux2.c,aux1.c,kxim)

        read!(t12,aux1,padded)
        rfft!(aux1)
        adddy!(aux2.c,aux1.c,kyim)

        read!(t13,aux1,padded)
        rfft!(aux1)
        adddz!(aux2.c,aux1.c,kzim)

        irfft!(aux2)
        write(f1,aux2.r)
    end

    begin #f3
        fill!(aux2.data,0.0)
        adddx!(aux2.c,aux1.c,kxim)

        read!(t33,aux1,padded)
        rfft!(aux1)
        adddz!(aux2.c,aux1.c,kzim)

        read!(t23,aux1,padded)
        rfft!(aux1)
        adddy!(aux2.c,aux1.c,kyim)

        irfft!(aux2)
        write(f3,aux2.r)
    end

    begin #f2
        fill!(aux2.data,0.0)
        adddz!(aux2.c,aux1.c,kzim)

        read!(t22,aux1,padded)
        rfft!(aux1)
        adddy!(aux2.c,aux1.c,kyim)

        read!(t12,aux1,padded)
        rfft!(aux1)
        adddx!(aux2.c,aux1.c,kxim)

        irfft!(aux2)
        write(f2,aux2.r)
    end

    return nothing
end

function div(u1::AbstractString,u2::AbstractString,u3::AbstractString,f::AbstractString)
    nx,ny,nz,lx,ly,lz = getdimsize()
    s = (nx,ny,nz)
    kxim = (2π .* rfftfreq(nx,lx) .* im)
    kyim = (2π .* fftfreq(ny,ly) .* im)
    kzim = (2π .* fftfreq(nz,lz) .* im)

    FFTW.set_num_threads(Threads.nthreads())
    isfile("fftw_wisdom") && FFTW.import_wisdom("fftw_wisdom")

    dtype, padded = checkinput(u1,nx,ny,nz)

    aux1 = PaddedArray(u1,s,padded)
    aux2 = similar(aux1)

    rfft!(aux1)
    adddx!(aux2.c,aux1.c,kxim)

    read!(u2,aux1,padded)
    rfft!(aux1)
    adddy!(aux2.c,aux1.c,kyim)

    read!(u3,aux1,padded)
    rfft!(aux1)
    adddz!(aux2.c,aux1.c,kzim)

    irfft!(aux2)
    write(f,aux2.r)

    return nothing
end

div(ARGS...)


text:
module TestReadWrite
using PhaseSpaceIO
using Test
using PhaseSpaceIO.Testing

end


text:
# Relevant parts of ODESystem interface
MTK.independent_variable(me::AbstractMeanfieldEquations) = me.iv
MTK.states(me::AbstractMeanfieldEquations) = me.states

function MTK.equations(me::AbstractMeanfieldEquations)
    # Get the MTK variables
    varmap = me.varmap
    vs = MTK.states(me)
    vhash = map(hash, vs)

    # Substitute conjugate variables by explicit conj
    vs′ = map(_conj, vs)
    vs′hash = map(hash, vs′)
    i = 1
    while i <= length(vs′)
        if vs′hash[i] ∈ vhash
            deleteat!(vs′, i)
            deleteat!(vs′hash, i)
        else
            i += 1
        end
    end
    rhs = [substitute_conj(eq.rhs, vs′, vs′hash) for eq∈me.equations]

    # Substitute to MTK variables on rhs
    subs = Dict(varmap)
    rhs = [substitute(r, subs) for r∈rhs]
    vs_mtk = getindex.(varmap, 2)

    # Return equations
    t = MTK.independent_variable(me)
    D = MTK.Differential(t)
    return [Symbolics.Equation(D(vs_mtk[i]), rhs[i]) for i=1:length(vs)]
end

# Substitute conjugate variables
function substitute_conj(t,vs′,vs′hash)
    if SymbolicUtils.istree(t)
        if t isa Average
            if hash(t)∈vs′hash
                t′ = _conj(t)
                return conj(t′)
            else
                return t
            end
        else
            _f = x->substitute_conj(x,vs′,vs′hash)
            args = map(_f, SymbolicUtils.arguments(t))
            return SymbolicUtils.similarterm(t, SymbolicUtils.operation(t), args)
        end
    else
        return t
    end
end

# Conversion to ODESystem
MTK.isparameter(::SymbolicUtils.Sym{<:Parameter}) = true

function MTK.ODESystem(me::AbstractMeanfieldEquations, iv=me.iv; kwargs...)
    eqs = MTK.equations(me)
    return MTK.ODESystem(eqs, iv; kwargs...)
end


text:
using Mamba
using Distributions

xx = randn(100)
yy = 0.5 .+ 1*xx + randn(100)

model = Model(

  y = Stochastic(1,
    (mu, s2) ->  MvNormal(mu, sqrt(s2)),
    false
  ),

  mu = Logical(1,
    (xmat, beta) -> xmat * beta,
    false
  ),

  beta1 = Stochastic(1,
    () -> Normal(0, sqrt(0.1)),false
  ),

  beta2 = Stochastic(1,
   () -> Normal(0, sqrt(1000)),false
  ),

  beta = Logical(1,
    (beta1,beta2) -> [beta1; beta2]
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )
)

scheme = [NUTS([:beta, :s2])]

setsamplers!(model, scheme)

line = Dict{Symbol, Any}(
  :x => xx,
  :y => yy
)

line[:xmat] = [ones(100) line[:x]]

inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta1 => rand(Normal(0, 1),1),
    :beta2 => rand(Normal(0, 1),1),
    :s2 => rand(Gamma(1, 1))
  ) for i in 1:3]

sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

describe(sim1)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
#ASUM: BLAS function for calculating sum of absolute values.
doc"""
  asum_naive(n::Integer, X::Array{Posit}, stride::Integer)

  calculates the sum of the first n strided values in X.
"""
function asum_naive{N,ES}(n::Integer, X::DenseArray{<:PositOrComplex{N,ES}}, incx::Integer)
    #see LAPACK sasum.f for code reference:
    # http://www.netlib.org/lapack/explore-html/df/d28/group__single__blas__level1_gafc5e1e8d9f26907c0a7cf878107f08cf.html#gafc5e1e8d9f26907c0a7cf878107f08cf

    #corner case checking.
    (n <= 0 || incx <= 0) && return zero(Posit{N,ES})
    #throwing a bounds error early allows us to use @inbounds optimization in the
    #for loop.
    (n*incx) > length(X) && throw(BoundsError(X, n*incx))

    #when the increment is 1, we can use the builtin mapreduce function which we
    #presume should be language-optimized.
    (incx == 1) && return mapreduce(abs, +, (n == length(X)) ? X : X[1:n])

    #for non-one increment, just do a for loop manually.
    accumulator = zero(Posit{N,ES})
    for idx in 1:incx:n*incx
        @inbounds accumulator += abs(X[idx])
    end
    accumulator
end

function BLAS.asum{N,ES}(n::Integer, X::DenseArray{<:PositOrComplex{N,ES}}, incx::Integer)

    #corner case checking.
    (n <= 0 || incx <= 0) && return zero(Posit{N,ES})
    (n*incx) > length(X) && throw(BoundsError(X, n*incx))

    zero!(__real_quire)
    if (incx == 1)
        foreach(p -> add!(__real_quire, abs(p)), n == length(X) ? X : X[1:n])
        return Posit{N,ES}(__real_quire)
    end

    #declare the local accumulator value.
    accumulator = zero(Posit{N,ES})
    for idx in 1:incx:n*incx
        accumulator = add!(__real_quire, abs(X[idx]))
    end
    accumulator
end

function axpy!{T<:Posit}(α, x::AbstractArray{T}, y::AbstractArray{T})
    n = _length(x)
    if n != _length(y)
        throw(DimensionMismatch("x has length $n, but y has length $(_length(y))"))
    end
    for (IY, IX) in zip(eachindex(y), eachindex(x))
        @inbounds y[IY] = fma(α, x[IX], y[IY])
    end
    y
end

function dot_naive{T<:Posit}(n::Integer, x::AbstractArray{T}, incx::Integer, y::AbstractArray{T}, incy::Integer)
    (n <= 0 || incx <= 0) && return zero(T)
    (n * incx) > length(x) && throw(BoundsError(y, n * incx))
    (n * incy) > length(y) && throw(BoundsError(y, n * incy))
    accumulator = zero(T)
    #TODO:  make this use the fused multiply add function.
    for (xidx,yidx) in zip(1:incx:n*incx, 1:incy:n*incy)
        @inbounds accumulator = accumulator + x[xidx] * y[yidx]
    end
    accumulator
end

function Base.dot{T<:Posit}(n::Integer, x::AbstractArray{T}, incx::Integer, y::AbstractArray{T}, incy::Integer)
    (n <= 0 || incx <= 0) && return zero(T)
    (n * incx) > length(x) && throw(BoundsError(y, n * incx))
    (n * incy) > length(y) && throw(BoundsError(y, n * incy))

    accumulator = zero(T)
    zero!(__real_quire)
    for (xidx,yidx) in zip(1:incx:n*incx, 1:incy:n*incy)
        @inbounds accumulator = fdp!(__real_quire, x[xidx], y[yidx])
    end
    accumulator
end

function dot_naive{T<:Posit}(n::Integer, x::AbstractArray{Complex{T}}, incx::Integer, y::AbstractArray{Complex{T}}, incy::Integer)
    (n <= 0 || incx <= 0) && return zero(T)
    (n * incx) > length(x) && throw(BoundsError(y, n * incx))
    (n * incy) > length(y) && throw(BoundsError(y, n * incy))

    real_accumulator = zero(T)
    imag_accumulator = zero(T)
    #TODO:  make this use the fused multiply add function.
    for (xidx,yidx) in zip(1:incx:n*incx, 1:incy:n*incy)
        @inbounds real_accumulator = real_accumulator + real(x[xidx]) * real(y[yidx])
        @inbounds real_accumulator = real_accumulator + imag(x[xidx]) * imag(y[yidx])
        @inbounds imag_accumulator = imag_accumulator + real(x[xidx]) * imag(y[yidx])
        @inbounds imag_accumulator = imag_accumulator - imag(x[xidx]) * real(y[yidx])
    end
    accumulator
end

function Base.dot{T<:Posit}(n::Integer, x::AbstractArray{Complex{T}}, incx::Integer, y::AbstractArray{Complex{T}}, incy::Integer)
    (n <= 0 || incx <= 0) && return zero(T)
    (n * incx) > length(x) && throw(BoundsError(y, n * incx))
    (n * incy) > length(y) && throw(BoundsError(y, n * incy))

    zero!(__real_quire)
    zero!(__imag_quire)
    real_accumulator = zero(T)
    imag_accumulator = zero(T)
    for (xidx,yidx) in zip(1:incx:n*incx, 1:incy:n*incy)
        @inbounds real_accumulator = fdp!(__real_quire, real(x[xidx]), real(y[yidx]))
        @inbounds real_accumulator = fdp!(__real_quire, imag(x[xidx]), imag(y[yidx]))
        @inbounds imag_accumulator = fdp!(__imag_quire, real(x[xidx]), imag(y[yidx]))
        @inbounds imag_accumulator = fdp!(__imag_quire, -imag(x[xidx]), real(y[yidx]))
    end
    real_accumulator + imag_accumulator * im
end

function dot_naive{T<:Posit}(x::AbstractArray{T}, y::AbstractArray{T})
    if length(x) != length(y)
        throw(length(x) > length(y) ? BoundsError(x, length(y) + 1) : BoundsError(y, length(x) + 1))
    end
    accumulator = zero(T)
    for idx in 1:length(x)
        @inbounds accumulator = accumulator + x[idx] * y[idx]
    end
    accumulator
end

function Base.dot{N,ES}(x::AbstractArray{<:PositOrComplex{N,ES}}, y::AbstractArray{<:PositOrComplex{N,ES}})
    if length(x) != length(y)
        throw(length(x) > length(y) ? BoundsError(x, length(y) + 1) : BoundsError(y, length(x) + 1))
    end
    naked_dot(x, y)
end

function dot_naive{T<:Posit}(x::AbstractArray{Complex{T}}, y::AbstractArray{Complex{T}})
    if length(x) != length(y)
        throw(length(x) > length(y) ? BoundsError(x, length(y) + 1) : BoundsError(y, length(x) + 1))
    end
    real_accumulator = zero(T)
    imag_accumulator = zero(T)
    for idx in 1:length(x)
        #note that for complex values, the first vector is conjugated.
        @inbounds real_accumulator = real_accumulator + real(x[xidx]) * real(y[yidx])
        @inbounds real_accumulator = real_accumulator + imag(x[xidx]) * imag(y[yidx])
        @inbounds imag_accumulator = imag_accumulator + real(x[xidx]) * imag(y[yidx])
        @inbounds imag_accumulator = imag_accumulator - imag(x[xidx]) * real(y[yidx])
    end
    real_accumulator + imag_accumulator * im
end

#to implement:  vector 2-norm
#unconjugated dot product

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
"""
    monotonicity(prob) -> Int

Whether the solution to `prob` must be decreasing (`-1`), constant (`0`) or increasing (`+1`) in `r`.
"""
function monotonicity end

"""
    DirichletProblem(eq; i=<initial value>, b=<boundary value>, ϕb=0)
    DirichletProblem(D; i=<initial value>, b=<boundary value>, ϕb=0)

Semi-infinite problem with a Dirichlet boundary condition.

# Arguments
- `eq::Equation`: governing equation.
- `D`: diffusivity function. Shortcut for `DirichletProblem(DiffusionEquation(D), ...)`.

# Keyword arguments
- `i`: initial value.
- `b`: imposed boundary value.
- `ϕb=0` (`\\phi<tab>b`): boundary constant for an optional moving boundary.
At time `t`, the boundary is located at `ϕb*√t`. Must be positive if `eq` is radial.

# Examples
```jldoctest; setup = :(using Fronts)
julia> D(θ) = θ^4
D (generic function with 1 method)

julia> prob = Fronts.DirichletProblem(D, i=1, b=2)
⎧ ∂θ/∂t = ∂(D(θ)*∂θ/∂r)/∂r, r>0,t>0
⎨ θ(r,0) = 1, r>0
⎩ θ(0,t) = 2, t>0
```

See also: [`Equation`](@ref)
"""
struct DirichletProblem{Teq,_T,_Tϕ}
    eq::Teq
    i::_T
    b::_T
    ϕb::_Tϕ
    function DirichletProblem(eq::Equation{1}; i, b, ϕb=0)
        new{typeof(eq),promote_type(typeof(i),typeof(b)),typeof(ϕb)}(eq, i, b, ϕb)
    end
    function DirichletProblem(eq::Equation; i, b, ϕb)
        @argcheck ϕb>0
        new{typeof(eq),promote_type(typeof(i),typeof(b)),typeof(ϕb)}(eq, i, b, ϕb)
    end
end

DirichletProblem(D; i, b, ϕb=0) = DirichletProblem(DiffusionEquation(D), i=i, b=b, ϕb=ϕb)

function Base.show(io::IO, prob::DirichletProblem)
    if iszero(prob.ϕb)
        println(io, "⎧ ", prob.eq, ", r>0,t>0")
    else
        println(io, "⎧ ", prob.eq, ", r>rb(t),t>0")
    end
    println(io, "⎨ ", prob.eq.symbol, "(r,0) = ", prob.i, ", r>0")
    if iszero(prob.ϕb)
        print(io, "⎩ ", prob.eq.symbol, "(0,t) = ", prob.b, ", t>0")
    else
        println(io, "⎩ ", prob.eq.symbol, "(rb(t),t) = ", prob.b, ", t>0")
        print(io, "with rb(t) = ", prob.ϕb, "*√t")
    end
end

monotonicity(prob::DirichletProblem)::Int = sign(prob.i - prob.b)

function d_dϕ(prob::DirichletProblem{<:DiffusionEquation}, symbol::Symbol)
    @argcheck symbol === :b_hint
    (prob.i - prob.b)/(2*√prob.eq.D(prob.b))
end

function d_dϕ(prob::DirichletProblem{<:RichardsEquation}, symbol::Symbol)
    @argcheck symbol === :b_hint
    C_ = prob.eq.C(prob.b)
    if C_≤zero(C_)
        return (prob.i - prob.b)
    end
    return (prob.i - prob.b)/(2*√(prob.eq.K(prob.b)/C_))
end

"""
    FlowrateProblem(eq; i=<initial value>, Qb=<boundary flowrate>, angle=2π, height=1, ϕb=0)

Semi-infinite radial (polar/cylindrical) problem with an imposed-flowrate boundary condition.

# Arguments
- `eq::Equation{2}`: governing equation.

# Keyword arguments
- `i`: initial value.
- `Qb`: imposed boundary flowrate.
- `angle=2π`: total angle covered by the domain.
- `height=1`: domain height.
- `ϕb=0` (`\\phi<tab>b`): boundary constant for an optional moving boundary.
At time `t`, the boundary is located at `ϕb*√t`.

# Examples
```jldoctest; setup = :(using Fronts)
julia> D(θ) = θ^4
D (generic function with 1 method)

julia> eq = Fronts.DiffusionEquation{2}(D)
∂θ/∂t = 1/r*∂(r*D(θ)*∂θ/∂r)/∂r

julia> prob = Fronts.FlowrateProblem(eq, i=1, Qb=1)
⎧ ∂θ/∂t = 1/r*∂(r*D(θ)*∂θ/∂r)/∂r, r>0,t>0
⎨ θ(r,0) = 1, r>0
⎩ Qb(0,t) = 1, t>0
```

See also: [`Equation`](@ref)
"""
struct FlowrateProblem{Teq,_Tθ,_Tϕ,_TQ,_Th}
    eq::Teq
    i::_Tθ
    Qb::_TQ
    _αh::_Th
    ϕb::_Tϕ

    function FlowrateProblem(eq::Equation{2}; i, Qb, angle=2π, height=1, ϕb=0)
        @argcheck 0<angle≤2π; @argcheck height>zero(height)
        @argcheck ϕb≥zero(ϕb)
        αh = angle*height
        new{typeof(eq),typeof(i),typeof(ϕb),typeof(Qb),typeof(αh)}(eq, i, Qb, αh, ϕb)
    end
end

function Base.show(io::IO, prob::FlowrateProblem)
    if iszero(prob.ϕb)
        println(io, "⎧ ", prob.eq, ", r>0,t>0")
    else
        println(io, "⎧ ", prob.eq, ", r>rb(t),t>0")
    end
    println(io, "⎨ ", prob.eq.symbol, "(r,0) = ", prob.i, ", r>0")
    if iszero(prob.ϕb)
        print(io, "⎩ Qb(0,t) = ", prob.Qb, ", t>0")
    else
        println(io, "⎩ Qb(rb(t),t) = ", prob.Qb, ", t>0")
        print(io, "with rb(t) = ", prob.ϕb, "*√t")
    end
end

function d_dϕ(prob::FlowrateProblem, symbol::Symbol; b, ϕb=prob.ϕb)
    @argcheck symbol === :b
    return d_dϕ(prob.eq, b, ϕb, flux_mul_r=prob.Qb/prob._αh)
end

monotonicity(prob::FlowrateProblem)::Int = -sign(prob.Qb)

"""
    CauchyProblem(eq; b=<boundary value>, d_dϕb=<boundary ϕ-derivative>, ϕb=0)
    CauchyProblem(D; b=<boundary value>, d_dϕb=<boundary ϕ-derivative>, ϕb=0)

Semi-infinite problem with a Cauchy boundary condition (and unknown initial condition).

# Arguments
- `eq::Equation`: governing equation.
- `D`: diffusivity function. Shortcut for `CauchyProblem(DiffusionEquation(D), ...)`.

# Keyword arguments
- `b`: imposed boundary value.
- `d_dϕb`: imposed value of the ϕ-derivative of the solution at the boundary, where ϕ is the Boltzmann variable.
This value is equivalent to `√t*∂_∂r(<solution>, :b, t)` at any time `t>0`.
- `ϕb=0` (`\\phi<tab>b`): boundary constant for an optional moving boundary.
At time `t`, the boundary is located at `ϕb*√t`. Must be positive if `eq` is radial.

# Examples
```jldoctest; setup = :(using Fronts)
julia> D(θ) = θ^4
D (generic function with 1 method)

julia> prob = Fronts.CauchyProblem(D, b=2, d_dϕb=-0.1)
⎧ ∂θ/∂t = ∂(D(θ)*∂θ/∂r)/∂r, r>0,t>0
⎨ θ(0,t) = 2, t>0
⎩ √t*∂θ/∂r(0,t) = -0.1, t>0
```

See also: [`Equation`](@ref)
"""
struct CauchyProblem{Teq,_T,_Tϕ,_Td_dϕ}
    eq::Teq
    b::_T
    d_dϕb::_Td_dϕ
    ϕb::_Tϕ
    function CauchyProblem(eq::Equation{1}; b, d_dϕb, ϕb=0)
        new{typeof(eq),typeof(b),typeof(ϕb),typeof(d_dϕb)}(eq, b, d_dϕb, ϕb)
    end
    function CauchyProblem(eq::Equation; b, d_dϕb, ϕb)
        @argcheck ϕb>zero(ϕb)
        new{typeof(eq),typeof(b),typeof(ϕb),typeof(d_dϕb)}(eq, b, d_dϕb, ϕb)
    end
end

CauchyProblem(D; b, d_dϕb, ϕb=0) = CauchyProblem(DiffusionEquation(D), b=b, d_dϕb=d_dϕb, ϕb=ϕb)

function Base.show(io::IO, prob::CauchyProblem)
    if iszero(prob.ϕb)
        println(io, "⎧ ", prob.eq, ", r>0,t>0")
        println(io, "⎨ ", prob.eq.symbol, "(0,t) = ", prob.b, ", t>0")
        print(io,   "⎩ √t*∂", prob.eq.symbol, "/∂r(0,t) = ", prob.d_dϕb, ", t>0")
    else
        println(io, "⎧ ", prob.eq, ",r>rb(t),t>0")
        println(io, "⎨ ", prob.eq.symbol, "(rb(t),t) = ", prob.b, ", t>0")
        println(io, "⎩ √t*∂", prob.eq.symbol, "/∂r(rb(t),t) = ", prob.d_dϕb, ", t>0")
        print(io, "with rb(t) = ", prob.ϕb, "*√t")
    end
end

monotonicity(prob::CauchyProblem)::Int = sign(prob.d_dϕb)


text:
using HDF5, JLD

include("CenterOfMass.jl")
include("setup_parameters.jl")

actuator_limit = 4.0
max_robots = 4

prior = initialize_prior(circle_ps, resolution, interior_q, masses)

data_file = "dict_file_$(max_robots)_$(actuator_limit)"

precomputed_normals = Array{CombinationMap}()

if isreadable(data_file)
  println("loading normals")
  @time @load data_file precomputed_normals
else
  println("precomputing normals")

  @time precomputed_normals = generate_combinations_maps(attachment_ws, max_robots, prior,
    circle_ws, actuator_limit, interior_q)

  println("saving normals")
  @time @save data_file precomputed_normals
end

println("succesfully got normals")

data = get_data(prior)[:]
@time for ii = 1:length(precomputed_normals)
  for key = keys(precomputed_normals[ii])
    normals = precomputed_normals[ii][key]
    csqmi = compute_mutual_information(data, normals)
  end
end

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
@testset "recv_reorder.jl" begin
    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        segment_arrives(receiver; data="abcd", seqno=isn + 10, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""

        @test unassembled_bytes(receiver) == 4
        @test assembled_bytes(receiver) == 0
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        segment_arrives(receiver; data="efgh", seqno=isn + 5, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""

        @test unassembled_bytes(receiver) == 4
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="abcd", seqno=isn + 1, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 9)
        @test read!(stream_out(receiver)) == "abcdefgh"

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 8

    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        segment_arrives(receiver; data="efgh", seqno=isn + 5, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""

        @test unassembled_bytes(receiver) == 4
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="ab", seqno=isn + 1, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 3)
        @test read!(stream_out(receiver)) == "ab"

        @test unassembled_bytes(receiver) == 4
        @test assembled_bytes(receiver) == 2

        segment_arrives(receiver; data="cd", seqno=isn + 3, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 9)
        @test read!(stream_out(receiver)) == "cdefgh"

        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 8
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)

        segment_arrives(receiver; data="e", seqno=isn + 5, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 1
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="g", seqno=isn + 7, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 2
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="c", seqno=isn + 3, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 3
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="ab", seqno=isn + 1, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 4)
        @test read!(stream_out(receiver)) == "abc"
        @test unassembled_bytes(receiver) == 2
        @test assembled_bytes(receiver) == 3

        segment_arrives(receiver; data="f", seqno=isn + 6, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 4)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 3
        @test assembled_bytes(receiver) == 3

        segment_arrives(receiver; data="d", seqno=isn + 4, result=SegmentArrives_Result_OK)
        @test read!(stream_out(receiver)) == "defg"
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 7
    end

    @testset begin receiver = TCPReceiver(2358)
        isn = rand(UInt32(0):typemax(UInt32))
        segment_arrives(receiver; with_syn=true, seqno=isn, result=SegmentArrives_Result_OK)

        @test ackno(receiver) == WrappingInt32(isn + 1)

        segment_arrives(receiver; data="e", seqno=isn + 5, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 1
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="g", seqno=isn + 7, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 2
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="c", seqno=isn + 3, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 1)
        @test read!(stream_out(receiver)) == ""
        @test unassembled_bytes(receiver) == 3
        @test assembled_bytes(receiver) == 0

        segment_arrives(receiver; data="abcdefgh", seqno=isn + 1, result=SegmentArrives_Result_OK)
        @test ackno(receiver) == WrappingInt32(isn + 9)
        @test read!(stream_out(receiver)) == "abcdefgh"
        @test unassembled_bytes(receiver) == 0
        @test assembled_bytes(receiver) == 8
    end

end

main()

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using Test, SparseArrays
import LuxurySparse: SparseMatrixCOO, allocated_coo, isdense

coo1 = SparseMatrixCOO(
    [1, 4, 2, 3, 3, 3],
    [1, 1, 2, 4, 3, 4],
    [0.1, 0.2, 0.4im, 0.5, 0.3, 0.5im],
    4,
    4,
)

@testset "basic" begin
    @test eltype(coo1) == ComplexF64
    @test size(coo1, 1) == size(coo1, 2) == 4

    coo2 = allocated_coo(ComplexF64, 4, 4, 6)
    @test eltype(coo2) == ComplexF64
    @test nnz(coo2) == 6
    @test size(coo2) == (4, 4)
    @test size(coo2, 1) == size(coo2, 2) == 4

    @test coo1 == copy(coo1)
    @test coo1 == copyto!(coo2, coo1)

    @test Matrix(coo1) == [0.1 0 0 0; 0 0.4im 0 0; 0 0 0.3 0.5+0.5im; 0.2 0 0 0]
    @test sparse(coo1) == sparse(findnz(coo1)...)
    @test isdense(coo1) == false

    A = sprand(50, 50, 0.2) |> Matrix
    @test SparseMatrixCOO(A) == A
end

@testset "sparse" begin
    p1 = copy(coo1)
    @test LuxurySparse.nonzeros(p1) == p1.vs
    p1.vs[2] = 0
    @test LuxurySparse.dropzeros!(p1) ==
          [0.1 0 0 0; 0 0.4im 0 0; 0 0 0.3 0.5+0.5im; 0 0 0 0]
end

@testset "basicmath" begin
    @test coo1 * 2 == Matrix(coo1) * 2
    @test 2 * coo1 == 2 * Matrix(coo1)
    @test coo1 / 2 == Matrix(coo1) / 2
    @test coo1 * 2 isa SparseMatrixCOO
    @test coo1 / 2 isa SparseMatrixCOO
    @test 2 * coo1 isa SparseMatrixCOO
end

@testset "setindex" begin
    coo[1, 1] = 0.2
    @test coo.is[end] == 1
    @test coo.js[end] == 1
    @test coo.vs[end] == 0.2

    S = SparseMatrixCOO{Float32, Cint}(Cint[], Cint[], Float32[], 5, 5)
    S[1, 1] = 3
    @test S[1, 1] == 3
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

function tau_model(events::AbstractVector, D₁::AbstractArray, Dₙ::AbstractArray)
  nevents = length(events)

  # trivial mass function
  nevents == 1 && return [1.0]

  # sources of information: primary + auxiliary₁ + auxiliary₂ + ...
  nsources = 1 + length(Dₙ)

  # distance matrix
  D = zeros(nevents, nsources)
  D[:,1] = D₁[events]
  for j=2:nsources
    D[:,j] = Dₙ[j-1][events]
  end

  # convert distances to ranks
  idx = mapslices(sortperm, D, dims=1)
  for j=1:nsources
    r = 0; prevdist = -Inf
    for i in view(idx,:,j)
      D[i,j] > prevdist && (r += 1)
      prevdist = D[i,j]
      D[i,j] = r
    end
  end

  # conditional probabilities
  P = nevents .- D .+ 1
  P = broadcast(/, P, sum(P, dims=1))

  # prior to data all events are equally probable
  x₀ = (1 - 1/nevents) / (1/nevents)

  # assume no redundancy
  X = (1 .- P) ./ P
  x = x₀ * prod(X/x₀, dims=2)

  p = 1 ./ (1 .+ x)
end