text:
"""
    one_fold(x,y,path_length,folds,fold) -> Vector{Float}

For a regularization path given by `models`,
this function computes an out-of-sample error for the fold given by `fold` in a `q`-fold crossvalidation scheme.
The folds are indexed by the `Int` vector `folds`.

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.
- `models` is an `Int` vector containing the regularization path to compute.
- `folds` is the `Int` array that indicates which data to hold out for testing.
- `fold` indexes the current fold.

Optional Arguments:

- `tol` is the convergence tolerance to pass to the path computations. Defaults to `1e-6`.
- `max_iter` caps the number of permissible iterations in the algorithm. Defaults to `100`.
- `window` is an `Int` to dictate the dimension of the search window for potentially exchanging predictors.
   Defaults to `max(20, min(maximum(models), size(x,2)))` or the greater of `20` and `r`, where `r` is the lesser of the maximum model size and the number of predictors.
   Decreasing this quantity tells the algorithm to search through fewer current active predictors,
   which can decrease compute time but can also degrade model recovery performance.
- `quiet` is a Boolean to activate output. Defaults to `true` (no output).
"""
function one_fold{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    models   :: DenseVector{Int},
    folds    :: DenseVector{Int},
    fold     :: Int;
    tol      :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    window   :: Int  = max(20, min(maximum(models), size(x,2))),
    quiet    :: Bool = true
)

    # find testing indices and size of test set
    test_idx = folds .== fold
    test_size = sum(test_idx)

    # train_idx is the vector that indexes the TRAINING set
    train_idx = !test_idx

    # allocate the arrays for the training set
    x_train = x[train_idx,:]
    y_train = y[train_idx]

    # will return a sparse matrix of betas
    betas = exlstsq(x_train, y_train, models=models, window=window, max_iter=max_iter, tol=tol, quiet=quiet)

    # compute the mean out-of-sample error for the TEST set
    errors = vec(sumabs2(broadcast(-, y[test_idx], x[test_idx,:] * betas), 1)) ./ (2*test_size)

    return errors :: Vector{T}
end

"""
    pfold(x, y, path, folds, nfolds) -> Vector

This function is the parallel execution kernel in `cv_exlstsq`. It is not meant to be called outside of `cv_exlstsq`.
It will distribute `nfolds` crossvalidation folds across the processes supplied by the optional argument `pids` and call `one_fold` for each fold.
Each fold will compute a regularization path given by the `Int` vector `models`.
`pfold` collects the vectors of MSEs from each process, sum-reduces them, and returns their average across all folds.
"""
function pfold{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T},
    models   :: DenseVector{Int},
    folds    :: DenseVector{Int},
    nfolds   :: Int;
    pids     :: DenseVector{Int} = procs(x),
    tol      :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    window   :: Int  = max(20, min(maximum(models), size(x,2))),
    quiet    :: Bool = true,
)
    # how many CPU processes can pfold use?
    np = length(pids)

    # report on CPU processes
    quiet || println("pfold: np = ", np)
    quiet || println("pids = ", pids)

    # set up function to share state (indices of folds)
    i = 1
    nextidx() = (idx=i; i+=1; idx)

    # preallocate cell array for results
    results = Array{Any}(nfolds)

    # master process will distribute tasks to workers
    # master synchronizes results at end before returning
    @sync begin

        # loop over all workers
        for worker in pids

            # exclude process that launched pfold, unless only one process is available
            if worker != myid() || np == 1

                # asynchronously distribute tasks
                @async begin
                    while true

                        # grab next fold
                        current_fold = nextidx()

                        # if current fold exceeds total number of folds then exit loop
                        current_fold > nfolds && break

                        # report distribution of fold to worker and device
                        quiet || print_with_color(:blue, "Computing fold $current_fold on worker $worker.\n\n")

                        # launch job on worker
                        # worker loads data from file paths and then computes the errors in one fold
                        results[current_fold] = remotecall_fetch(worker) do
                            one_fold(x, y, models, folds, current_fold, tol=tol, max_iter=max_iter, window=window, quiet=quiet)
                        end # end remotecall_fetch()
                    end # end while
                end # end @async
            end # end if
        end # end for
    end # end @sync

    # return reduction (row-wise sum) over results
    return (reduce(+, results[1], results) ./ nfolds) :: Vector{T}
end

"""
    cv_exlstsq(x, y, models, q) -> ELSQResults

This function will perform `q`-fold crossvalidation for the best model size in the exchange algorithm for least squares regression.
Each path is asynchronously spawned using any available indexed processor.
The function to compute each path, `one_fold`, will return a vector of mean out-of-sample errors (MSEs).
`cv_exlstsq` reduces across the `q` folds yield averaged MSEs for each model size.

Arguments:

- `x` is the `n` x `p` design matrix.
- `y` is the `n`-vector of responses.

Optional Arguments:

- `models` is an integer vector to specify the the regularization path to compute.
- `q` is the number of folds to compute. Defaults to `max(3, min(CPU_CORES, 5))`.
- `folds` is the `Int` array that indicates which data to hold out for testing. Defaults to `cv_get_folds(sdata(y),q)`.
- `folds` is the partition of the data. Defaults to a random partition into `q` disjoint sets.
- `tol` is the convergence tolerance to pass to the path computations. Defaults to `1e-6`.
- `max_iter` caps the number of permissible iterations in the exchange algorithm. Defaults to `100`.
- `quiet` is a `Bool` to activate output. Defaults to `true` (no output).
   *NOTA BENE*: each processor outputs feed to the console without regard to the others,
   so setting `quiet=false` can yield very messy output!

Output:

An `ELSQResults` container object with the following fields:
- `mses` is a vector of averaged MSEs across all folds, with one component per model computed.
- `k` is the best crossvalidated model size. Here model size refers to the number of nonzero regression coefficients.
- `b`, a vector of `k` components with the computed effect sizes of the best model β.
- `bidx`, an `Int` vector of `k` components indicating the support of `b` at `k`.
"""
function cv_exlstsq{T <: Float}(
    x        :: DenseMatrix{T},
    y        :: DenseVector{T};
    models   :: DenseVector{Int} = collect(1:min(20,size(x,2))),
    q        :: Int              = set_cv_fold_num(3, 5),
    pids     :: DenseVector{Int} = procs(),
    folds    :: DenseVector{Int} = cv_get_folds(sdata(y),q),
    tol      :: T    = convert(T, 1e-6),
    max_iter :: Int  = 100,
    window   :: Int  = max(20, min(maximum(models), size(x,2))),
    quiet    :: Bool = true,
)

#    1 <= minimum(models) <= maximum(models) <= size(x,2) || throw(ArgumentError("Model sizes must be positive and cannot exceed number of predictors"))

    mses = pfold(x, y, models, folds, q, pids=pids, tol=tol, max_iter=max_iter, quiet=quiet, window=window)

    # what is the best model size?
    k = convert(Int, floor(mean(models[mses .== minimum(mses)])))

    # print results
    quiet || print_cv_results(mses, models, k)

    # refit coefficients
    b, bidx = refit_exlstsq(x, y, k, models=models, tol=tol, max_iter=max_iter, quiet=quiet)

    return ELSQCrossvalidationResults(mses, b, bidx, k, sdata(models))
end


text:
"""
Main module for `ParameterEstimocean.jl` -- a Julia software for parameter estimation
for Oceananigans models using Ensemble Kalman Inversion.
"""
module ParameterEstimocean

export
    SyntheticObservations,
    BatchedSyntheticObservations,
    InverseProblem,
    FreeParameters,
    Transformation,
    SpaceIndices,
    TimeIndices,
    RescaledZScore,
    ZScore,
    forward_map,
    forward_run!,
    observation_map,
    observation_times,
    observation_map_variance_across_time,
    ensemble_column_model_simulation,
    ConcatenatedOutputMap,
    ConcatenatedVectorNormMap,
    eki,
    lognormal,
    ScaledLogitNormal,
    iterate!,
    EnsembleKalmanInversion,
    Resampler,
    NormExceedsMedian,
    FullEnsembleDistribution,
    SuccessfulEnsembleDistribution,
    ConstantConvergence

include("Utils.jl")
include("Transformations.jl")
include("Observations.jl")
include("Parameters.jl")
include("EnsembleSimulations.jl")
include("InverseProblems.jl")
include("EnsembleKalmanInversions.jl")
include("PseudoSteppingSchemes.jl")

using .Utils
using .Transformations
using .Observations
using .EnsembleSimulations
using .Parameters
using .InverseProblems
using .EnsembleKalmanInversions
using .PseudoSteppingSchemes

#####
##### Data!
#####

using DataDeps

function __init__()
    # Register LESbrary data
    lesbrary_url = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/glw/lesbrary2/LESbrary/idealized/"

    cases = ["free_convection",
             "weak_wind_strong_cooling",
             "strong_wind_weak_cooling",
             "strong_wind",
             "strong_wind_no_rotation"]

    two_day_suite_url = lesbrary_url * "two_day_suite/"

    glom_url(suite, resolution, case) = string(lesbrary_url,
                                               suite, "/", resolution, "_resolution/",
                                               case, "_instantaneous_statistics.jld2")

    two_day_suite_1m_paths  = [glom_url( "two_day_suite", "2m_2m_1m", case) for case in cases]
    two_day_suite_2m_paths  = [glom_url( "two_day_suite", "4m_4m_2m", case) for case in cases]
    two_day_suite_4m_paths  = [glom_url( "two_day_suite", "8m_8m_4m", case) for case in cases]
    four_day_suite_1m_paths = [glom_url("four_day_suite", "2m_2m_1m", case) for case in cases]
    four_day_suite_2m_paths = [glom_url("four_day_suite", "4m_4m_2m", case) for case in cases]
    four_day_suite_4m_paths = [glom_url("four_day_suite", "8m_8m_4m", case) for case in cases]
    six_day_suite_1m_paths  = [glom_url( "six_day_suite", "2m_2m_1m", case) for case in cases]
    six_day_suite_2m_paths  = [glom_url( "six_day_suite", "4m_4m_2m", case) for case in cases]
    six_day_suite_4m_paths  = [glom_url( "six_day_suite", "8m_8m_4m", case) for case in cases]

    dep = DataDep("two_day_suite_1m", "Idealized 2 day simulation data with 1m vertical resolution", two_day_suite_1m_paths)
    DataDeps.register(dep)

    dep = DataDep("two_day_suite_2m", "Idealized 2 day simulation data with 2m vertical resolution", two_day_suite_2m_paths)
    DataDeps.register(dep)

    dep = DataDep("two_day_suite_4m", "Idealized 2 day simulation data with 4m vertical resolution", two_day_suite_4m_paths)
    DataDeps.register(dep)

    dep = DataDep("four_day_suite_1m", "Idealized 4 day simulation data with 1m vertical resolution", four_day_suite_1m_paths)
    DataDeps.register(dep)

    dep = DataDep("four_day_suite_2m", "Idealized 4 day simulation data with 2m vertical resolution", four_day_suite_2m_paths)
    DataDeps.register(dep)

    dep = DataDep("four_day_suite_4m", "Idealized 4 day simulation data with 4m vertical resolution", four_day_suite_4m_paths)
    DataDeps.register(dep)

    dep = DataDep("six_day_suite_1m", "Idealized 6 day simulation data with 1m vertical resolution", six_day_suite_1m_paths)
    DataDeps.register(dep)

    dep = DataDep("six_day_suite_2m", "Idealized 6 day simulation data with 2m vertical resolution", six_day_suite_2m_paths)
    DataDeps.register(dep)

    dep = DataDep("six_day_suite_4m", "Idealized 6 day simulation data with 4m vertical resolution", six_day_suite_4m_paths)
    DataDeps.register(dep)
end

end # module

ext: jl


text:
"""
    VariationalConvBE(filter, in => out, rank,
                    ensemble_size, σ = identity;
                    stride = 1, pad = 0, dilation = 1,
                    groups = 1, [bias, weight, init])
    VariationalConvBE(layer, alpha_sampler, gamma_sampler,
                    ensemble_bias, ensemble_act, rank)

Creates a bayesian conv BatchEnsemble layer.
Batch ensemble is a memory efficient alternative for deep ensembles. In deep ensembles,
if the ensemble size is N, N different models are trained, making the time and memory
complexity O(N * complexity of one network).
BatchEnsemble generates weight matrices for each member in the ensemble using a
couple of rank 1 vectors R (alpha), S (gamma), RS' and multiplying the result with
weight matrix W element wise. We also call R and S as fast weights. In the bayesian
version of batch ensemble, instead of having point estimates of the fast weights, we
sample them form a trainable parameterized distribution.

Reference - https://arxiv.org/abs/2005.07186

During both training and testing, we repeat the samples along the batch dimension
N times, where N is the ensemble_size. For example, if each mini batch has 10 samples
and our ensemble size is 4, then the actual input to the layer has 40 samples.
The output of the layer has 40 samples as well, and each 10 samples can be considered
as the output of an esnemble member.

# Fields
- `layer`: The conv layer which transforms the pertubed input to output
- `alpha_sampler`: Sampler for the first Fast weight of size (in_dim, ensemble_size)
- `gamma_sampler`: Sampler for the second Fast weight of size (out_dim, ensemble_size)
- `ensemble_bias`: Bias added to the ensemble output, separate from conv layer bias
- `ensemble_act`: The activation function to be applied on ensemble output
- `rank`: Rank of the fast weights (rank > 1 doesn't work on GPU for now)

# Arguments
- `filter::NTuple{N,Integer}`: Kernel dimensions, eg, (5, 5)
- `ch::Pair{<:Integer,<:Integer}`: Input channels => output channels
- `rank::Integer`: Rank of the fast weights
- `ensemble_size::Integer`: Number of models in the ensemble
- `σ::F=identity`: Activation of the dense layer, defaults to identity
- `init=glorot_normal`: Initialization function, defaults to glorot_normal
- `bias::Bool=true`: Toggle the usage of bias in the dense layer
- `ensemble_bias::Bool=true`: Toggle the usage of ensemble bias
- `ensemble_act::F=identity`: Activation function for enseble outputs
- `alpha_init=TrainableMvNormal`: Initialization function for the alpha fast weight,
                        defaults to TrainableMvNormal
- `gamma_init=TrainableMvNormal`: Initialization function for the gamma fast weight,
                        defaults to TrainableMvNormal
"""
struct VariationalConvBE{L,F,M,B}
    layer::L
    alpha_sampler::M
    gamma_sampler::M
    ensemble_bias::B
    ensemble_act::F
    rank
end

function VariationalConvBE(
    k::NTuple{N,Integer},
    ch::Pair{<:Integer,<:Integer},
    rank::Integer,
    ensemble_size::Integer,
    σ = identity;
    init = glorot_normal,
    stride = 1,
    pad = 0,
    dilation = 1,
    groups = 1,
    bias = true,
    ensemble_bias = true,
    ensemble_act = identity,
    alpha_init = TrainableMvNormal,
    gamma_init = TrainableMvNormal,
    device = cpu,
) where {N}
    layer = Flux.Conv(
        k,
        ch,
        σ;
        stride = stride,
        pad = pad,
        dilation = dilation,
        init = init,
        groups = groups,
        bias = bias,
    )
    in_dim = ch[1]
    out_dim = ch[2]
    if rank >= 1
        alpha_shape = (in_dim, ensemble_size)
        gamma_shape = (out_dim, ensemble_size)
    else
        error("Rank must be >= 1.")
    end
    # Alpha fast weight sampler
    alpha_sampler = alpha_init(alpha_shape, init = init, device = device)
    # Gamma fast weight sampler
    gamma_sampler = gamma_init(gamma_shape, init = init, device = device)
    gamma = gamma_sampler()
    ensemble_bias = create_bias(gamma, ensemble_bias, out_dim, ensemble_size)

    return VariationalConvBE(
        layer,
        alpha_sampler,
        gamma_sampler,
        ensemble_bias,
        ensemble_act,
        rank,
    )
end

@functor VariationalConvBE

function (be::VariationalConvBE)(x)
    # Conv Batch Ensemble params
    layer = be.layer
    alpha = be.alpha_sampler()
    gamma = be.gamma_sampler()
    e_b = be.ensemble_bias
    e_σ = be.ensemble_act

    batch_size = size(x)[end]
    in_size = size(alpha)[1]
    out_size = size(gamma)[1]
    ensemble_size = size(alpha)[2]
    samples_per_model = batch_size ÷ ensemble_size

    # TODO: Do we really need repeat? Or can we use broadcast?
    # Alpha, gamma shapes - [units, ensembles, rank]
    e_b = repeat(e_b, samples_per_model)
    alpha = repeat(alpha, samples_per_model)
    gamma = repeat(gamma, samples_per_model)
    # Reshape alpha, gamma to [units, batch_size, rank]
    e_b = reshape(e_b, (1, 1, out_size, batch_size))
    alpha = reshape(alpha, (1, 1, in_size, batch_size))
    gamma = reshape(gamma, (1, 1, out_size, batch_size))

    perturbed_x = x .* alpha
    output = layer(perturbed_x) .* gamma
    output = e_σ.(output .+ e_b)

    return output
end


text:
include("./script_resnet.jl")
##

using BenchmarkTools
using Flux: cpu, gpu

@profview self_play_step!(envsim)
@enter self_play_step!(envres)
@timed self_play_step!(envsim)
@timed self_play_step!(envres)
@timed learning_step!(envres)
@profview learning_step!(envres)
@profview self_play_step!(envres)

game = GI.init(gspec)
nns = envres.bestnns
player = MuPlayer((InitialOracle(nns), RecurrentOracle(nns)), envres.params.self_play.mcts)

@enter AlphaZero.think(player, game)
AlphaZero.think(player, game)
@profview for i in 1:100; AlphaZero.think(player, game); end;

batchsize = 8
actionbatch = [rand(1:5) for _ in 1:batchsize]

obsbatch = []
for a in actionbatch
  GI.play!(game, a)
  obs = GI.current_state(game)
  push!(obsbatch, obs)
end

obsbatch

statebatchres = evaluate_batch(envres.bestnns.h, obsbatch)
statebatchsim = evaluate_batch(envsim.bestnns.h, obsbatch)

pbatch_vbatchres = evaluate_batch(envres.bestnns.f, statebatchres)
pbatch_vbatchsim = evaluate_batch(envsim.bestnns.f, statebatchsim)
envsim.bestnns.f

rbatch_snextbatch = evaluate_batch(envres.bestnns.g, collect(zip(statebatchres, actionbatch)))
rbatch_snextbatch = evaluate_batch(envsim.bestnns.g, collect(zip(statebatchsim, actionbatch)))
size(rbatch_snextbatch[1])

envres.bestnns.f(hiddenstate)

r, s =envres.bestnns.g((hiddenstate,3))

res = Vector{Any}(undef,8)
res[1:8] = collect(rbatch_snextbatch)

##
# for batchifierbatchsize in 1:4
batchifierbatchsize = 32
obsbatch = (rand(rand(envsim.memory).states) for _ in 1:batchifierbatchsize)

statebatchres = (rand(Float32,5,4,128) for _ in 1:batchifierbatchsize)
actionbatch = (rand(1:5) for _ in 1:batchifierbatchsize)
gresbatch = zip(statebatchres, actionbatch)

statebatchsim = (rand(Float32, 64) for _ in 1:batchifierbatchsize)
gsimbatch = zip(statebatchsim, actionbatch)

benchmarks = Dict()

#f
benchmarks["f_res_batch_cpu"] = @benchmark evaluate_batch($(envres.bestnns.f),       $collect(statebatchres))
benchmarks["f_res_batch_gpu"] = @benchmark evaluate_batch($(envres.bestnns.f|>gpu),  $collect(statebatchres))

benchmarks["f_sim_batch_cpu"] = @benchmark evaluate_batch($(envsim.bestnns.f),       $collect(statebatchsim))
benchmarks["f_sim_batch_gpu"] = @benchmark evaluate_batch($(envsim.bestnns.f|>gpu),  $collect(statebatchsim))

#g
benchmarks["g_res_batch_cpu"] = @benchmark evaluate_batch($(envres.bestnns.g),       $collect(gresbatch))
benchmarks["g_res_batch_gpu"] = @benchmark evaluate_batch($(envres.bestnns.g|>gpu),  $collect(gresbatch))

benchmarks["g_sim_batch_cpu"] = @benchmark evaluate_batch($(envsim.bestnns.g),       $collect(gsimbatch))
benchmarks["g_sim_batch_gpu"] = @benchmark evaluate_batch($(envsim.bestnns.g|>gpu),  $collect(gsimbatch))

#h
benchmarks["h_res_batch_cpu"] = @benchmark evaluate_batch($(envres.bestnns.h),       $collect(obsbatch))
benchmarks["h_res_batch_gpu"] = @benchmark evaluate_batch($(envres.bestnns.h|>gpu),  $collect(obsbatch))

benchmarks["h_sim_batch_cpu"] = @benchmark evaluate_batch($(envsim.bestnns.h),       $collect(obsbatch))
benchmarks["h_sim_batch_gpu"] = @benchmark evaluate_batch($(envsim.bestnns.h|>gpu),  $collect(obsbatch))

bench_batchsize[batchifierbatchsize] = benchmarks
# end
##
@profview for _ in 1: 1000; evaluate_batch(g_net, collect(gbatch)); end;

# FileIO.save("bench_batchsize2.jld2", "bench_batchsize", bench_batchsize)
bench_batchsize = FileIO.load("bench_batchsize2.jld2", "bench_batchsize")
using BenchmarkPlots, StatsPlots

X = 4:80
# for x in X
  ts = bench_batchsize[x]["h_res_batch_gpu"].times
  t̄ = mean(ts)
  t̃ = median(ts)
  σ = std(ts)

scatter([x],[t̃])

plot!([x],[t̄], ribbon=[σ])
plotly()
bench_batchsize[x]["h_res_batch_gpu"]

evaluate_batch((envres.bestnns.f),       collect(statebatchres))
evaluate_batch((envres.bestnns.f|>gpu),  collect(statebatchres))

evaluate_batch((envsim.bestnns.f),       collect(statebatchsim))
evaluate_batch((envsim.bestnns.f|>gpu),  collect(statebatchsim))

#g
evaluate_batch((envres.bestnns.g),       collect(gresbatch))
typeof(evaluate_batch((player.recurrent_inference),       collect(gresbatch))[1][3])
evaluate_batch((envres.bestnns.g|>gpu),  collect(gresbatch))

evaluate_batch((envsim.bestnns.g),       collect(gsimbatch))
evaluate_batch((envsim.bestnns.g|>gpu),  collect(gsimbatch))

#h
evaluate_batch((envres.bestnns.h),       collect(obsbatch))
evaluate_batch((init),       collect(obsbatch))
evaluate_batch((envres.bestnns.h|>gpu),  collect(obsbatch))

evaluate_batch((envsim.bestnns.h),       collect(obsbatch))
evaluate_batch((envsim.bestnns.h|>gpu),  collect(obsbatch))

evaluate(init,collect(obsbatch)[1])

# state = GI.current_state(game)
# obsbatch = [state, state]

nns = deepcopy(envres.bestnns)|>Flux.gpu|>Flux.testmode!
player = MuPlayer((InitialOracle(nns), RecurrentOracle(nns)), envres.params.self_play.mcts)

collected_obsbatch = collect(obsbatch)
@profview for _ in 1:1000; initres = evaluate_batch((player.initial_inference), collected_obsbatch); end;
@timed for _ in 1:1000; initres = evaluate_batch((player.initial_inference), collected_obsbatch); end;
@profview for _ in 1:1000; recurres = evaluate_batch((player.recurrent_inference), collect(gresbatch)); end;
@timed for _ in 1:1000; recurres = evaluate_batch((player.recurrent_inference), collect(gresbatch)); end;

collect(obsbatch)
typeof(recurres)
typeof(initres)

p = rand(Float32,5,32)
v = rand(Float32,1,32)
s = rand(Float32,5,4,128,32)

cp = cp|>Flux.gpu
cv = cv|>Flux.gpu
cs = cs|>Flux.gpu
@benchmark (p,v,s)|>Flux.gpu
@benchmark (cp,cv,cs)|>Flux.cpu
@benchmark map(Flux.cpu, (cp,cv,cs))
@benchmark map(x->Flux.unstack(x,ndims(x)), (p,v,s))
@benchmark map(x->Flux.unstack(x,ndims(x)), (cp,cv,cs))
@benchmark P,V,S = map($convert_output, $(cp,cv,cs))

@benchmark evaluate_batch($(player.initial_inference), $collect(obsbatch));

P, V, S = map(convert_output, (cp,cv,cs))
@benchmark V = [v[1] for v in $V]
@code_typed evaluate_batch((player.initial_inference), collect(obsbatch))

# only forward
X = Flux.batch(GI.vectorize_state(gspec, b) for b in obsbatch) #obsrvation
Xnet = to_nndevice(player.initial_inference.f, X)
function initforwardonly(init, Xnet)
  S = forward(init.h, Xnet)
  P, V = forward(init.f, S)
  return P, V, S
end
function recurforwardonly(recur, S_A_net)
  R, S = forward(recur.g, S_A_net)
  P, V = forward(recur.f, S)
  return P, V, S, R
end
S = Flux.batch(b[1] for b in gresbatch)
A = Flux.batch(encode_a(gspec, b[2]; batchdim=4) for b in gresbatch)
S_A = cat(S,A,dims=3)
S_A_net = S_A|>gpu

@timed for _ in 1:1000; recurforwardonly(player.recurrent_inference, S_A_net); end;

P,V,S = initforwardonly((player.initial_inference), Xnet)
@benchmark P|>Flux.cpu
@benchmark V|>Flux.cpu
@benchmark Flux.cpu(S)
@benchmark

@benchmark initforwardonly($(player.initial_inference), Xnet)
@profview for _ in 1:1000; P,V,S = initforwardonly((player.initial_inference), Xnet); end;
@timed for _ in 1:1000; P,V,S = initforwardonly((player.initial_inference), Xnet); end;

function newhash(x)
  x = vec(x)
  h = hash(x)
  # x = reshape(x, 5,4,3)
  return h, x
end

@benchmark hash($rand(Float32,5,4,128))
@benchmark hash($rand(Float32,5,4,3))
@benchmark hash($rand(Float32,64))
@benchmark hash($rand(Float32,32))
@benchmark hash($GI.current_state(game))
@benchmark flatten($rand(Float32,5,4,3))
@benchmark reshape($rand(Float32,60), 5,4,3)
@benchmark newhash($rand(Float32, 5,4,3))

@enter hash(rand(Float32,5,4,128))

@benchmark hash($rand(Int))
@benchmark hash($rand(Int8))
@benchmark hash(($rand(Float32,5,4,3),$rand(Float32)))
@benchmark hash(($rand(Float32,5,4,3), $rand(Int)))
@benchmark hash($rand(Float32,5,12))

@benchmark hash($rand(128,128,128))
@benchmark hash($vec(rand(128,128,128)))
@benchmark hash($rand(128^3))

Base.summarysize(rand(128^3))
Base.summarysize(rand(128,128,128))


text:
using NMRTools
using Plots
using LsqFit
using Measurements
using ColorSchemes
using LightGraphs

## input files and settings

# filenames
peaklistfilename = "peak.list"
spectrafilename = "999/test.ft2"

# set size of ROI for fitting
stripwidthH = 2 * 0.025; # ppm
stripwidthC = 2 * 1.2; # ppm

# relaxation times (in s)
τ = [0.]

# load data
spectra = loadnmr(spectrafilename)
noise = spectra[:noise]

windowC = spectra[Y, :window]
windowH = spectra[X, :window]

## includes

include("multiplet_calculations.jl")

## utility functions

findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2]

## load peak list

methyls = []
peakH = Vector{Float64}()
peakC = Vector{Float64}()
println("Reading peak list ($peaklistfilename):")
open(peaklistfilename) do file
    for ln in eachline(file)
        fields = split(ln)
        println("\t$(fields[1])\t$(fields[2])\t$(fields[3])")
        push!(methyls, fields[1])
        push!(peakH, parse(Float64, fields[3]))
        push!(peakC, parse(Float64, fields[2]))
    end
end
npeaks = length(methyls)
println("$npeaks peak positions loaded.")
println()

## find non-overlapping clusters
overlaps(i, j) = (abs(peakH[i]-peakH[j]) < stripwidthH) && (abs(peakC[i]-peakC[j]) < stripwidthC)
adjacencymatrix = [overlaps(i,j) for i=1:npeaks, j=1:npeaks]
clusters = connected_components(SimpleGraph(adjacencymatrix))
#clusters = [[1], [2, 3]]
nclusters = length(clusters)
println("$nclusters clusters of peaks found:")
for cluster in clusters
    println("\t$(foldr((a,b)->a*" "*b, methyls[cluster]))")
end
println()

## loop over clusters and fit
fitresults = []
for cluster in clusters

    ncluster = length(cluster)
    peaknames = methyls[cluster]
    peakgroup = foldr((a,b)->a*"-"*b, peaknames)
    peakx = peakH[cluster]
    peaky = peakC[cluster]

    # get spectrum region covering selected peaks
    xmin = minimum(peakx) - stripwidthH/2
    xmax = maximum(peakx) + stripwidthH/2
    ymin = minimum(peaky) - stripwidthC/2
    ymax = maximum(peaky) + stripwidthC/2

    # extract data ROI into arrays
    δobsH = xval(spectra[NMRTools.Between(xmin,xmax),NMRTools.Between(ymin,ymax)])
    δobsC = yval(spectra[NMRTools.Between(xmin,xmax),NMRTools.Between(ymin,ymax)])
    dat = data(spectra[NMRTools.Between(xmin,xmax),NMRTools.Between(ymin,ymax)]) / noise
    dat = reshape(dat, size(dat)..., 1)
    nx, ny, nt = size(dat)

    # form mask for ROIs
    mask = zeros(Bool, (nx, ny))
    for i=1:ncluster
        ix = (peakx[i] - stripwidthH/2) .<= δobsH .<= (peakx[i] + stripwidthH/2)
        iy = (peaky[i] - stripwidthC/2) .<= δobsC .<= (peaky[i] + stripwidthC/2)
        mask[ix, iy] .= true
    end
    mask3 = cat([mask for i=1:nt]..., dims=3)
    maskeddat = dat[mask3]

    simmultiplet(dH, dC, R2H, R2C, S2tc, csa, J) = multiplet([dH, dC, R2H, R2C, S2tc, csa, J], δobsH, δobsC, τ, spectra[X,:bf], spectra[Y,:bf], windowH, windowC)
    # p = (ncluster x 4), p[i, :] = [R2H, R2C, S2tc, csa]
    function resid(p)
        simmultiplets = zeros(Float64, (nx, ny, nt, ncluster))
        maskedmultiplets = zeros(Float64, (length(maskeddat), ncluster))
        ysim = zeros(Float64, (nx, ny, nt))
        for i = 1:ncluster
            m = simmultiplet(peakx[i], peaky[i], p[i, 1], p[i, 2], p[i, 3], p[i, 4], 124.5)
            simmultiplets[:,:,:,i] = m
            maskedmultiplets[:, i] = m[mask3]
        end
        A = maskedmultiplets \ maskeddat
        return maskeddat - maskedmultiplets*A
    end

    # second round of fitting, including peak positions and J coupling
    function sim(p)
        simmultiplets = zeros(Float64, (nx, ny, nt, ncluster))
        maskedmultiplets = zeros(Float64, (length(maskeddat), ncluster))
        ysim = zeros(Float64, (nx, ny, nt))
        for i = 1:ncluster
            m = simmultiplet(p[i, 1], p[i, 2], p[i, 3], p[i, 4], p[i, 5], p[i, 6], p[i, 7])
            simmultiplets[:,:,:,i] = m
            maskedmultiplets[:, i] = m[mask3]
        end
        A = maskedmultiplets \ maskeddat
        for i=1:ncluster
            simmultiplets[:,:,:,i] = A[i] * simmultiplets[:,:,:,i]
        end
        return sum(simmultiplets, dims=4)[:,:,:,1]
    end
    function resid2(p)
        simmultiplets = zeros(Float64, (nx, ny, nt, ncluster))
        maskedmultiplets = zeros(Float64, (length(maskeddat), ncluster))
        ysim = zeros(Float64, (nx, ny, nt))
        for i = 1:ncluster
            m = simmultiplet(p[i, 1], p[i, 2], p[i, 3], p[i, 4], p[i, 5], p[i, 6], p[i, 7])
            simmultiplets[:,:,:,i] = m
            maskedmultiplets[:, i] = m[mask3]
        end
        A = maskedmultiplets \ maskeddat
        return maskeddat - maskedmultiplets*A
    end
    p1 = zeros(Float64, ncluster, 7)
    p1[:, 1] = peakx
    p1[:, 2] = peaky
    p1[:, 3:6] .= [50.0 40.0 10.0 30.0]
    p1[:, 7] .= 124.5

    fit = LsqFit.lmfit(p -> resid2(reshape(p, :, 7)), vec(p1), Float64[], show_trace=true, maxIter=100, x_tol=1e-3, g_tol=1e-6, autodiff=:finiteforward)

    pfit = reshape(fit.param, :, 7)
    pfiterr = reshape(stderror(fit), :, 7)
    p = pfit .± pfiterr
    ysim = sim(pfit)

    # store results
    for i=1:ncluster
        push!(fitresults, (peaknames[i], p[i,:]))
    end

    # plot results - overlay of contour plots
    clev = 10 .^LinRange(0.7, 3, 15)
    plots = []
    for i=1:nt
        # 2D overlay
        plt = heatmap(reverse(δobsH), reverse(δobsC), mask[end:-1:1,end:-1:1]')
        contour!(plt, reverse(δobsH), reverse(δobsC), dat[end:-1:1,end:-1:1,i]', c=:dodgerblue, levels=clev,
          colorbar=false, xflip=true, yflip=true, xguide="δH / ppm", yguide="δC / ppm", aspect_ratio=0.1)
        contour!(plt, reverse(δobsH), reverse(δobsC), ysim[end:-1:1,end:-1:1,i]', c=:orangered, levels=clev)
        xlims!(xmin, xmax)
        ylims!(ymin, ymax)
        push!(plots, plt)
    end
    plt = plot(plots..., size=(200,600))
    display(plt)
    savefig(plt, "contour-$(peakgroup).pdf")

end

## print results
println()
println("methyl\tdH / ppm\tdC / ppm \tR2H / s-1\tR2C / s-1\tS2tc / ns\tcsa / ppm\tJ / Hz")
println("methyl\tdH / ppm\tdC / ppm \tR2H / s-1\tR2C / s-1\tS2tc / ns\tcsa / ppm\tJ / Hz")
for fr in fitresults
    println("$(fr[1])\t$(fr[2][1])\t$(fr[2][2])\t$(fr[2][3])\t$(fr[2][4])\t$(fr[2][5])\t$(fr[2][6])\t$(fr[2][7])")
end
println()


text:
using ReachabilityAnalysis

r₀ = Interval(1.19, 1.21)
θ₀ = Interval(0.49, 0.51)
vr₀ = Singleton([0.0])
vθ₀ = Singleton([0.0])

X₀ = r₀ × θ₀ × vr₀ × vθ₀;

Dict([
    "opts" => Dict(:T=>5.0, :alg=>TMJets(orderT=8, orderQ=2, maxsteps=1_000)),
    "X0"   => X₀,
    "info" => Dict(
        "name" => "spring_pendulum",
        "dim" => 4,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
import Distributions.Categorical

function resample(log_weights::Array{Float64,1}, n::Int)
    dist = exp.(log_weights - logsumexp(log_weights))
    rand(Categorical(dist), n)
end

immutable SMCScheme
    initializer
    incrementers::Array
    num_particles::Int
end

function smc(scheme::SMCScheme)
    num_steps = length(scheme.incrementers) + 1
    num_particles = scheme.num_particles
    log_weights = Array{Float64,1}(num_particles)
    log_ml_estimate::Float64 = 0.0

    # a particle is an array of components of arbitrary type
    particles = Array{Any,1}(num_particles)
    new_particles = Array{Any,1}(num_particles)

    # initialize particles and log_weights
    for i=1:num_particles
        particles[i], log_weights[i] = forward(scheme.initializer)
    end
    log_ml_estimate += (logsumexp(log_weights) - log(num_particles))

    # iterate through steps
    for t=2:num_steps
        parents = resample(log_weights, num_particles)
        for i=1:num_particles
            parent_particle = particles[parents[i]]
            x, l = forward(scheme.incrementers[t-1], particles[parents[i]])
            new_particles[i], log_weights[i] = (x, l)
        end
        log_ml_estimate += (logsumexp(log_weights) - log(num_particles))
        tmp = particles
        particles = new_particles
        new_particles = tmp
    end

    # pick output particle
    output_index = resample(log_weights, 1)[1]
    output = particles[output_index]
    (output, log_ml_estimate)
end

function conditional_smc(scheme::SMCScheme, output::Any)
    num_steps = length(scheme.incrementers) + 1
    num_particles = scheme.num_particles
    log_weights = Array{Float64,1}(num_particles)
    log_ml_estimate::Float64 = 0.0

    # a particle is an array of components of arbitrary type
    particles = Array{Any,1}(num_particles)
    new_particles = Array{Any,1}(num_particles)

    # select indices of ancestors of the output particle
    ancestry = rand(Categorical(num_particles), num_steps)

    # sample values for the ancestors of the output particle
    ancest_particles = Array{Any,1}(num_steps)
    ancest_particles[end] = output
    ancest_log_weights = Array{Float64,1}(num_steps);
    for t=num_steps:-1:2
        x, l = backward(scheme.incrementers[t-1], ancest_particles[t])
        ancest_particles[t-1], ancest_log_weights[t] = (x, l)
    end
    ancest_log_weights[1] = backward(scheme.initializer, ancest_particles[1])

    # initialize particles and log_weights
    for i=1:num_particles
        if ancestry[1] == i
            x, l = (ancest_particles[1], ancest_log_weights[1])
        else
            x, l = forward(scheme.initializer)
        end
        particles[i], log_weights[i] = (x, l)
    end
    log_ml_estimate += (logsumexp(log_weights) - log(num_particles))

    # iterate through steps
    for t=2:num_steps
        parents = resample(log_weights, num_particles)
        for i=1:num_particles
            if ancestry[t] == i
                x, l = ancest_particles[t], ancest_log_weights[t]
            else
                parent_particle = particles[parents[i]]
                x, l = forward(scheme.incrementers[t-1], parent_particle)
            end
            new_particles[i], log_weights[i] = (x, l)
        end
        log_ml_estimate += (logsumexp(log_weights) - log(num_particles))
        tmp = particles
        particles = new_particles
        new_particles = tmp
    end
    log_ml_estimate
end


text:
"""
Implements configuration B, i.e. BLG in Bernal stacking with self energy on both the ZZ and
the AC edges. There is a different phase difference between left and right sides of the
weak link (-ϕ/2 for x<0 ϕ/2 for x>0).
Neither smoothness nor leads are implemented. Also disorder and rotation are not considered.
"""

function rectangle_weaklink(p, selfy = true, hole = false; Δx_mask = 0, Δy_mask = 0)
    (; Ls, Δ, Ws, B, Ln, W) = p
    (; model0, field!, modelinter) = modelS(p)
    lat_top, lat_bot = latBLG(p, mask = hole, Δx_mask = Δx_mask, Δy_mask = Δy_mask) #0 angle rotation
    #old lat_top, lat_bot = latBLG(p,0)
    println("Flux: ",B[3]*Ln*W/Φ0)
    # PEIERLS PHASES
    #model paper the phase at the sc contacts is not correctly defined
    # diagphi(φ) = Diagonal(SA[cis(φ), cis(φ), cis(-φ), cis(-φ)])
    # piecewise(y) = clamp(y, -W/2, W/2) # clamp(x, -Ln/2 - a0/(2*√3) , Ln/2 + a0/(2*√3)))
    # A(r, B) = [-B[3]/ħoec, 0, 0] * piecewise(r[2])
    # eφ(r, dr, B) = diagphi(dot(A(r, B), dr))
    # peierls! =@hopping!((t, r, dr) -> t * eφ(r, dr, B))

    diagphi(φ) = Diagonal(SA[cis(φ), cis(φ), cis(-φ), cis(-φ)])
    piecewise(x) =  clamp(x, -Ln+Ls, Ln-Ls)
    A(r, B) = [0, B[3]/ħoec, 0] * piecewise(r[1])
    eφ(r, dr, B) = diagphi(dot(A(r, B), dr))
    peierls! =@hopping!((t, r, dr) -> t * eφ(r, dr, B))

    # LOCAL SELF-ENERGY MODEL
    xmin, xmax = extrema(r->r[1], sitepositions(Quantica.combine(lat_top, lat_bot)))
    ymin, ymax = extrema(r->r[2], sitepositions(Quantica.combine(lat_top, lat_bot)))
    self_region(r) = ifelse(selfy == false, !(xmin+Ls <= r[1] <= xmax - Ls),
        !(xmin+Ls <= r[1] <= xmax - Ls) || !(ymin+Ws <= r[2] <= ymax - Ws))

    sCself! = @onsite!((o, r; ϕ) -> o + self_region(r) * Δ *
        eφ(r, r + ifelse(sign(r[1])<0, [-Ln+Ls,-W/2,0], [Ln-Ls,-W/2,0]), B) *
        diagphi(sign(r[1])* ϕ/4) * σyτy *
            diagphi(sign(r[1])*ϕ/4)' * eφ(r, r+ifelse(sign(r[1])<0, [-Ln+Ls,-W/2,0], [Ln-Ls,-W/2,0]), B)')
        # eφ(r, r+ifelse(sign(r[1][0,0,0], B) * diagphi(sign(r[1])* ϕ/4) * σyτy *
        #     diagphi(sign(r[1])*ϕ/4)' * eφ(r, r+[0,0,0], B)')

    # HAMILTONIAN BUILD
    h_top = lat_top |> hamiltonian(model0; orbitals = Val(4)) |>
        unitcell(mincoordination = 5)
    h_bot = lat_bot |> hamiltonian(model0; orbitals = Val(4)) |>
        unitcell(mincoordination = 5)
    ph = Quantica.combine(h_top, h_bot; coupling = modelinter) |>
        parametric(peierls!, sCself!)
            return ph
end

"""
    spectrumvsphase(philist, p; nev = 16, kw...)
computes a number nev of eigenvalues as a function of the phase difference across the weak
link.
    see: rectangle_weakling(::Params()), Params(), modelS()
"""
function spectrumvsphase(philist, p; nev = 16, kw...)
    ph = rectangle_weaklink(p; kw...)
    elist = zeros(Float64, nev, length(philist))
    for i in 1:length(philist)
        s = spectrum(ph(phi = philist[i]), method = ArpackPackage(sigma = 1e-7, nev = nev))
        elist[:,i] = s.energies
    end
    return philist, elist
end

"flux test"
function quantumflux(p)
    Φ0 = 2067.8338 # T*nm^2
    println("Set Bz to: ",Φ0/(p.Ln*p.W))
    return Φ0/(p.Ln*p.W)
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libx264"], :libx264),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/jpsamaroo/x264Builder/releases/download/v2018.2.12-noyasm"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.aarch64-linux-gnu.tar.gz", "2a0ab4f4b1ae142abf2b8cadb3acb1c176e1d29921ee52749f01c46e6c5051d8"),
    Linux(:aarch64, :musl) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.aarch64-linux-musl.tar.gz", "85e9a923faa455dedb6b4d06829799f49b6942f8eb8c8d5957c79b9e06cc33ae"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.arm-linux-gnueabihf.tar.gz", "e1173cb209ce8e9b0f4c9f623044b0a9fec03ec76d43e5e39e0824ef3070eb66"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.arm-linux-musleabihf.tar.gz", "610cf15612ca054537cc8a0845f55304562163662c75eb8a5c57d41c4f44196c"),
    Linux(:i686, :glibc) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.i686-linux-gnu.tar.gz", "e4ac5eba78397f85e3968664651d897f68c4ba9f0a84677569158eb652295ca1"),
    Linux(:i686, :musl) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.i686-linux-musl.tar.gz", "a731e80af58cf3f3203a598331b41bb1dc8f4d0f16696609ff29bf85f24426d4"),
    Windows(:i686) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.i686-w64-mingw32.tar.gz", "f449c7e177a1a45be63417487996166131d4d857b0d26e4ab0527b19a6ca3ff9"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.powerpc64le-linux-gnu.tar.gz", "d95e0eaa336b0977071d4a5212def6c354f02219c2c699c9917c47b238497f78"),
    MacOS(:x86_64) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.x86_64-apple-darwin14.tar.gz", "c51834f3bc6f95601a4e8434b24d809f61683133feccaf61670e6974c7c4ed82"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.x86_64-linux-gnu.tar.gz", "65daed583973ae54787bf528f48ca6ce23375aee19a428a6045d6ea7436b1f3a"),
    Linux(:x86_64, :musl) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.x86_64-linux-musl.tar.gz", "99c67163be9810476f39478ef3b64d82bf5d237592557eeac3690ab8aeb3a692"),
    FreeBSD(:x86_64) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.x86_64-unknown-freebsd11.1.tar.gz", "a437b7ada02d263470da92a273cdc34e447bfa56cae7d79764c9b6120692e6dd"),
    Windows(:x86_64) => ("$bin_prefix/x264Builder.v2018.2.12-pre-noyasm.x86_64-w64-mingw32.tar.gz", "70b2e53a82cca36e4aa9774a2be39443b690696da6b15534c49338bd78928972"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose, ignore_platform = true)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps_codec.jl"), products)


text:
using MeshBasisFunctions
using Documenter

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

makedocs(;
    modules=[MeshBasisFunctions],
    authors="Sam Urmy <oceanographerschoice@gmail.com>, John K Best <isposdef@gmail.com>",
    repo="https://github.com/jkbest2/MeshBasisFunctions.jl/blob/{commit}{path}#{line}",
    sitename="MeshBasisFunctions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jkbest2.github.io/MeshBasisFunctions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jkbest2/MeshBasisFunctions.jl",
    devbranch="main",
)


text:
"""
Operator to compute the far field of a current distribution. In particular, given the current distribution ``j`` this operator allows for the computation of

```math
A j = n × ∫_Ω j e^{γ ̂x ⋅ y} dy
```

where ``̂x`` is the unit vector in the direction of observation. Note that the assembly routing expects the observation directions to be normalised by the caller.
"""
struct VIEFarField3D{K,P}
  gamma::K
  tau::P
end

function VIEFarField3D(;wavenumber=error("wavenumber is a required argument"), tau=nothing)
  gamma = nothing

  if iszero(real(wavenumber))
    gamma = -imag(wavenumber)
  else
    gamma = wavenumber*im
  end

  tau == nothing && (tau = x->1.0)

  VIEFarField3D(gamma, tau)
end

quaddata(op::VIEFarField3D,rs,els) = quadpoints(rs,els,(3,))
quadrule(op::VIEFarField3D,refspace,p,y,q,el,qdata) = qdata[1,q]

kernelvals(op::VIEFarField3D,y,p) = exp(op.gamma*dot(y,cartesian(p)))
function integrand(op::VIEFarField3D,krn,y,f,p)
    (y × (op.tau(cartesian(p)) * krn * f[1])) × y
end


text:
using DataFrames, Query

df1 = DataFrame(a=[1,2,3], b=[1.,2.,3.])
df2 = DataFrame(c=[2.,4.,2.], d=["John", "Jim","Sally"])

x = @from i in df1 begin
    @join j in df2 on i.a equals convert(Int,j.c)
    @select {i.a,i.b,j.c,j.d,e="Name: $(j.d)"}
    @collect DataFrame
end

println(x)


text:
@testset "Table" begin
    v1 = [1,2]
    v2 = [2.0, 4.0]
    @test @inferred(Table(l"a" => v1, l"b" => v2)) isa Table{(:a, :b), Tuple{Vector{Int}, Vector{Float64}}}

    table = Table(l"a" => v1, l"b" => v2)
    @test @inferred(table[1]) === Row(l"a" => 1, l"b" => 2.0)
    @test @inferred(table[2]) === Row(l"a" => 2, l"b" => 4.0)
    @test_throws BoundsError  table[3]

    @test names(table) == (:a, :b)
    @test @inferred(labels(table)) == (l"a", l"b")
end


text:
# for fully-observable case
struct MoveTowards <: Policy end

function POMDPs.action(p::MoveTowards, s::LTState)
    # try to sneak up diagonally
    diff = s.opponent-s.robot
    dx = sign(diff[1])
    dy = sign(diff[2])
    if abs(dx) != 1 && abs(dy) == 1
        return DIR_TO_ACTION[SVector(dx, dy)]
    else
        if abs(dx) == 1
            dx = 0
        end
        if abs(dy) == 1
            dy = 0
        end
        return DIR_TO_ACTION[SVector(dx, dy)]
    end
end

struct MoveTowardsSampled{RNG} <: Policy
    rng::RNG
end
MoveTowardsSampled() = MoveTowardsSampled(Random.GLOBAL_RNG)

function POMDPs.action(p::MoveTowardsSampled, b)
    s = rand(p.rng, b)
    return action(MoveTowards(), s)
end

struct OptimalMLSolver <: Solver
    solver::Solver
end

struct OptimalML{P<:Policy} <: Policy
    fo_policy::P
end

solve(sol::OptimalMLSolver, p::Union{MDP,POMDP}) = OptimalML(solve(sol.solver, p))
POMDPs.action(pol::OptimalML, b) = action(pol.fo_policy, mode(b))

struct BestExpectedSolver <: Solver
    solver::Solver
end

struct BestExpected{P<:Policy} <: Policy
    fo_policy::P
end

solve(sol::BestExpectedSolver, p::Union{MDP,POMDP}) = BestExpected(solve(sol.solver, p))
function POMDPs.action(pol::BestExpected, b)
    best_eu = -Inf
    best_s = first(iterator(b))
    for s in iterator(b)
        eu = pdf(b, s)*value(pol.fo_policy, s)
        if eu > best_eu
            best_eu = eu
            best_s = s
        end
    end
    return action(pol.fo_policy, best_s)
end


text:
function ma1kalman_variance(T, V0, F, W, H, A)
  # Note that this not the usual Kalman filter because the noise in
  # dlogk is an MA(1) instead of independent. We use
  # https://arxiv.org/pdf/1909.10582.pdf
  # and adopt their notation
  varR = Vector{typeof(V0)}(undef,T)
  varRprior = Vector{typeof(V0)}(undef,T)
  K = Vector{typeof(V0)}(undef,T)
  varz = Vector{typeof(V0)}(undef, T)
  zcoef = Matrix{typeof(V0)}(undef,2, T)
  varZ = MMatrix{2,2,typeof(V0)}(zeros(4))
  covzZ = MVector{2, typeof(V0)}(zeros(2))
  for t in 1:T
    varZ .= 0
    if t == 1
      varRprior[t] = V0
    else
      varRprior[t] = F*varR[t-1]*F' + W
      varZ[2,2] = H*varRprior[t-1]*H'
    end
    varZ[1,1] = H*varRprior[t]*H'
    varZ .+= A
    varz[t] = varZ[1,1]
    covzZ[1] = varZ[1,2]
    ivarZ = inv(varZ)
    zcoef[:,t] .= vec(covzZ'*ivarZ)
    iL = 1/(varZ[1,1] + dot(zcoef[:,t], covzZ))
    J = varRprior[t]*H
    K[t] = J*iL
    varR[t] =  varRprior[t]-K[t]*J' # P_t
  end
  return(varR, varRprior, K, varz, zcoef)
end

function ma1kalman_mean!(x, zhat, z, F, H, x0, K, zcoef)
  T = length(z)
  elag = @MVector zeros(eltype(zcoef[1,1]*z[1]), 2)
  for t in 1:T
    xminus = F*(t==1 ? x0 : x[t-1])
    ez = H*xminus
    zhat[t] = ez
    zhat[t] += zcoef[:,t]'*elag
    x[t] = xminus + K[t]*(z[t] - zhat[t])
    elag[2] = elag[1]
    elag[1] = z[t] - H*xminus
  end
end

function ma1kalman_mean!(x, zhat, z, F, H, x0, K, zcoef, xshift)
  T = length(z)
  elag = @MVector zeros(eltype(zcoef[1,1]*z[1]), 2)
  for t in 1:T
    xminus = F*(t==1 ? x0 : x[t-1] ) + xshift[t]
    ez = H*xminus
    zhat[t] = ez
    #zhat[t] += zcoef[:,t]'*elag
    x[t] = xminus + K[t]*(z[t] - zhat[t])
    elag[2] = elag[1]
    elag[1] = z[t] - H*xminus
  end
end

function kalman(dlogk, σR, σk, σR0, μR0, γ, ρ=1, Xα=zeros(length(dlogk)); MA1=false)
  T = length(dlogk)
  A = MA1 ?  SMatrix{2,2}([2*σk^2  -σk^2; -σk^2  2*σk^2]) :  SMatrix{2,2}([2*σk^2  0; 0  2*σk^2])
  meanR = Vector{typeof(σR)}(undef, T)
  zhat = Vector{typeof(σR)}(undef, T)
  varR, varRprior, K, varz, zcoef = ma1kalman_variance(T, σR0^2, ρ, σR^2, γ, A)
  z = dlogk .+ γ .- γ* (Xα)
  ma1kalman_mean!(meanR, zhat, z, ρ, γ, μR0, K, zcoef)
  meanR .+= Xα
  zhat .+= -γ .+γ*(Xα)
  return(meanR, varR, zhat)
end


text:
 cmd = "eventcreate /T INFORMATION /ID 123 /D \"Rosetta Code Write to Windows event log task example\""
 Base.run(`$cmd`)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    # tr neglog
    ((3, MatNegLog()),),
    ((3, MatNegLogEigOrd()),),
    ((3, MatNegLogDirect()),),
    # tr negentropy
    ((3, MatNegEntropy()),),
    ((3, MatNegEntropyEigOrd()),),
    # tr power12
    ((3, MatPower12(1.5)),),
    ((3, MatPower12EigOrd(1.5)),),
    # tr neg2sqrt
    ((3, MatNeg2Sqrt()),),
    ((3, MatNeg2SqrtEigOrd()),),
    ]
insts["fast"] = [
    ((40, MatNegGeom()),),
    ((15, MatNegGeomEFExp()),),
    ((15, MatNegGeomEFPow()),),
    ((40, MatInv()),),
    ((8, MatInvEigOrd()),),
    ((20, MatInvDirect()),),
    ((30, MatNegLog()),),
    ((12, MatNegLogEigOrd()),),
    ((15, MatNegLogDirect()),),
    ((30, MatNegEntropy()),),
    ((8, MatNegEntropyEigOrd()),),
    ((30, MatPower12(1.5)),),
    ((10, MatPower12EigOrd(1.5)),),
    ]
insts["various"] = [
    ((100, MatNegGeom()),),
    ((50, MatNegGeomEFExp()),),
    ((30, MatNegGeomEFPow()),),
    ((100, MatInv()),),
    ((12, MatInvEigOrd()),),
    ((50, MatInvDirect()),),
    ((100, MatNegLog()),),
    ((18, MatNegLogEigOrd()), nothing, (default_tol_relax = 1000,)),
    ((80, MatNegLogDirect()),),
    ((75, MatNegEntropy()),),
    ((14, MatNegEntropyEigOrd()),),
    ((30, MatPower12(1.5)),),
    ((6, MatPower12EigOrd(1.5)),),
    ]
return (CovarianceEstJuMP, insts)


text:
using StatGeochem
using Test, Statistics, StatsBase

# Utilities
@testset "Import" begin include("testImport.jl") end
@testset "Resampling" begin include("testResampling.jl") end
@testset "Changepoint" begin include("testChangepoint.jl") end
@testset "Geochemistry" begin include("testGeochemistry.jl") end
@testset "Other Utilities" begin include("testUtilities.jl") end

# Resources
@testset "Crust 1.0" begin include("testCrust1.jl") end
@testset "Other Resources" begin include("testResources.jl") end


text:
#using Pkg
#Pkg.activate(".")

push!(LOAD_PATH, "Simulation")
using Base.Threads
using JSON
using ArgParse
using Random
using FileIO
using JLD2
using DataFrames
using ProgressMeter

import Simulation

function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table! s begin
		"--JSON"
		help = "path to a JSON file with the parameter settings. If parameters are given sepeartely, the parameters of the JSON are overwritten in the program."
		"--output_prefix"
		help = "path to the output file"
	end
return parse_args(s)
end

function read_params(json, rng::AbstractRNG)
  constant_kernel_param = json["transmission_probabilities"]["constant"]  |> float
  household_kernel_param = json["transmission_probabilities"]["household"] |> float
  hospital_kernel_param = get(json["transmission_probabilities"], "hospital", 0.0) |> float
  friendship_kernel_param = get(json["transmission_probabilities"], "friendship", 0.0) |> float

  mild_detection_prob = json["detection_mild_proba"]  |> float

  tracking_prob = json["contact_tracking"]["probability"]  |> float
  tracking_backward_delay = json["contact_tracking"]["backward_detection_delay"]  |> float
  tracking_forward_delay = json["contact_tracking"]["forward_detection_delay"]  |> float
  testing_time = json["contact_tracking"]["testing_time"]  |> float

  phone_tracking_usage = json["phone_tracking"]["usage"] |> float
  phone_tracking_testing_delay = json["phone_tracking"]["detection_delay"] |> float

  population_path = json["population_path"] # <= JSON
  population_path::AbstractString # checks if it was indeed a string

  individuals_df = load(population_path)["individuals_df"]

  infection_modulation_name, infection_modulation_params = if !haskey(json, "modulation")
    nothing, NamedTuple{}()
  else
    modulation = json["modulation"]
    params = get(modulation, "params", Dict{String,Any}())
    modulation["function"], NamedTuple{Tuple(Symbol.(keys(params)))}(values(params))
  end

  spreading = get(json, "spreading", nothing)
  spreading_alpha = isnothing(spreading) ? nothing : spreading["alpha"]
  spreading_x0 = isnothing(spreading) ? 1 : get(spreading, "x0", 1)
  spreading_truncation = isnothing(spreading) ? Inf : get(spreading, "truncation", Inf)

  Simulation.load_params(
    rng,
    population = individuals_df,

    mild_detection_prob = mild_detection_prob,

    constant_kernel_param = constant_kernel_param,
    household_kernel_param = household_kernel_param,
    hospital_kernel_param = hospital_kernel_param,
    friendship_kernel_param = friendship_kernel_param,

    backward_tracking_prob = tracking_prob,
    backward_detection_delay = tracking_backward_delay,

    forward_tracking_prob = tracking_prob,
    forward_detection_delay = tracking_forward_delay,

    testing_time = testing_time,

    phone_tracking_usage = phone_tracking_usage,
    phone_detection_delay = phone_tracking_testing_delay,

    infection_modulation_name=infection_modulation_name,
    infection_modulation_params=infection_modulation_params,

    spreading_alpha=spreading_alpha,
    spreading_x0=spreading_x0,
    spreading_truncation=spreading_truncation
)
end
const OptTimePoint = Union{Missing, Simulation.TimePoint}

struct DetectionCallback
    detection_times::Vector{OptTimePoint}
    detection_types::Vector{UInt8}

    tracking_times::Vector{OptTimePoint}
    tracking_sources::Vector{UInt32}
    tracking_types::Vector{UInt8}

    max_num_infected::UInt32
end
DetectionCallback(sz::Integer, max_num_infected::Integer=10^8) = DetectionCallback(
    Vector{OptTimePoint}(missing, sz),
    fill(UInt8(0), sz),
    Vector{OptTimePoint}(missing, sz),
    fill(UInt32(0), sz),
    fill(UInt8(0), sz),
    max_num_infected
)

optreal2float(optreal::Union{Missing,T} where T<:Real) = ismissing(optreal) ? NaN32 : Float32(optreal)

function saveparams(dict, cb::DetectionCallback, prefix::AbstractString="")
  dict[prefix*"detection_times"] = optreal2float.(cb.detection_times)
  dict[prefix*"detection_types"] = cb.detection_types

  dict[prefix*"tracking_times"] = optreal2float.(cb.tracking_times)
  dict[prefix*"tracking_sources"] = cb.tracking_sources
  dict[prefix*"tracking_types"] = cb.tracking_types
end

function reset!(cb::DetectionCallback)
  fill!(cb.detection_times, missing)
  fill!(cb.detection_types, 0)
  fill!(cb.tracking_sources, 0)
  fill!(cb.tracking_types, 0)
end

function (cb::DetectionCallback)(event::Simulation.Event, state::Simulation.SimState, params::Simulation.SimParams)
  eventkind = Simulation.kind(event)
  contactkind = Simulation.contactkind(event)
  subject = Simulation.subject(event)
  if Simulation.isdetection(eventkind)
    cb.detection_times[subject] = Simulation.time(event)
    cb.detection_types[subject] = Simulation.detectionkind(event) |> UInt8
  elseif Simulation.istracking(eventkind)
    cb.tracking_times[subject] = Simulation.time(event)
    cb.tracking_sources[subject] = Simulation.source(event)
    cb.tracking_types[subject] = Simulation.trackingkind(event) |> UInt8
  end
  return Simulation.numinfected(state.stats) < cb.max_num_infected
end

function save_infections_and_detections(path::AbstractString, simstate::Simulation.SimState, callback::DetectionCallback)
  f = jldopen(path, "w", compress=true)
  try
    Simulation.saveparams(f, simstate)
    saveparams(f, callback)
  finally
    close(f)
  end
  nothing
end

function main()
  # set more threads by setting JULIA_NUM_THREADS environment variable
  @info "Stated" nthreads()

  parsed_args = parse_commandline()
  # parse the JSON file if provided
  if parsed_args["JSON"] === nothing
    println("give me sim_params, option: --JSON path")
    exit()
  end

  json = JSON.parsefile(parsed_args["JSON"])
  output_prefix = ""
  if parsed_args["output_prefix"] !== nothing
    output_prefix = parsed_args["output_prefix"]
  end

  max_num_infected = json["stop_simulation_threshold"] |> Int
  num_trajectories = json["num_trajectories"] |> Int
  num_initial_infected = json["initial_conditions"]["cardinalities"]["infectious"] |> Int
  params_seed = get(json, "params_seed", 0)

  @info "loading population and setting up parameters" params_seed
  rng = MersenneTwister(params_seed)
  params = read_params(json, rng)
  num_individuals =  Simulation.numindividuals(params)
  param_save_path = output_prefix*"run_params.jld2"

  @info "saving parameters" num_individuals param_save_path
  jldopen(param_save_path, "w", compress=false) do dict
    Simulation.saveparams(dict, params)
  end

  states = [Simulation.SimState(num_individuals, seed=123) for _ in 1:nthreads()]
  callbacks = [DetectionCallback(num_individuals, max_num_infected) for _ in 1:nthreads()]

  @info "starting simulation" num_trajectories
  writelock = ReentrantLock()
  progress = ProgressMeter.Progress(num_trajectories)
  GC.gc()
  #@threads
  for trajectory_id in 1:num_trajectories
    state = states[threadid()]
    Simulation.reset!(state, trajectory_id)
    Simulation.initialfeed!(state, num_initial_infected)

    callback = callbacks[threadid()]
    reset!(callback)
    #try
      Simulation.simulate!(state, params, callback)
      try
        lock(writelock) # JLD2 is not thread-safe, not even when files are separate
        save_infections_and_detections(output_prefix*"run_$trajectory_id.jld2", state, callback)
      finally
        GC.gc()
        unlock(writelock)
      end
    #catch err
    #  println(stderr, "Failed on thread ", threadid(), " iteration ", trajectory_id, " failed: ", err)
    #  foreach(x -> println(stderr, x), stacktrace(catch_backtrace()))
    #end

    ProgressMeter.next!(progress) # is thread-safe
  end
end

main()


text:
module NCTiles

using NetCDF, MeshArrays, Printf, NCDatasets

#include("utilities.jl")
#export read_nctiles

import Base: getindex

export NCvar, BinData, NCData, TileData, readbin
export createfile, addDim, addVar, addData, addDimData,
        readncfile, parsemeta, readAvailDiagnosticsLog

include("writenctiles.jl")

end # module


text:
using DiffEqBase, Test

condtion= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  t - 2.95
end

affect! = function (integrator)
  integrator.u = integrator.u + 2
end

rootfind = true
save_positions = (true,true)
callback = ContinuousCallback(condtion,affect!;save_positions=save_positions)

cbs= CallbackSet(nothing)
@test typeof(cbs.discrete_callbacks) <: Tuple
@test typeof(cbs.continuous_callbacks) <: Tuple
cbs = CallbackSet(callback,nothing)
@test typeof(cbs.discrete_callbacks) <: Tuple
@test typeof(cbs.continuous_callbacks) <: Tuple
cbs = CallbackSet(callback,CallbackSet())
@test typeof(cbs.discrete_callbacks) <: Tuple
@test typeof(cbs.continuous_callbacks) <: Tuple

condtion = function (integrator)
  true
end
affect! = function (integrator) end
save_positions = (true,false)
saving_callback = DiscreteCallback(condtion,affect!;save_positions=save_positions)

cbs1 = CallbackSet(callback,saving_callback)

@test length(cbs1.discrete_callbacks) == 1
@test length(cbs1.continuous_callbacks) == 1

cbs2 = CallbackSet(callback)
@test length(cbs2.continuous_callbacks) == 1
@test length(cbs2.discrete_callbacks) == 0

cbs3 = CallbackSet(saving_callback)
@test length(cbs3.discrete_callbacks) == 1
@test length(cbs3.continuous_callbacks) == 0

cbs4 = CallbackSet()
@test length(cbs4.discrete_callbacks) == 0
@test length(cbs4.continuous_callbacks) == 0

cbs5 = CallbackSet(cbs1,cbs2)

@test length(cbs5.discrete_callbacks) == 1
@test length(cbs5.continuous_callbacks) == 2


text:
using LinuxPerf
using Test

import LinuxPerf: make_bench, enable!, disable!, reset!, reasonable_defaults, counters

@testset "LinuxPerf" begin

@testset "simple benchmark" begin
    @test begin
        bench = make_bench(reasonable_defaults);
        @noinline function g(a)
            enable!(bench)
            c = 0
            for x in a
                if x > 0
                    c += 1
                end
            end
            disable!(bench)
            c
        end
        g(zeros(10000))

        counters(bench)

        true  # Succeeded without any exceptions...
    end
end

@testset "@measure" begin
    c1 = @measure 2 + 2
    v, b1 = @measured 2 + 2

    @test v === 4
    @test typeof(c1) == typeof(counters(b1))
end

end


text:
@testset "set fluents" begin

# Load domain and problem
path = joinpath(dirname(pathof(PDDL)), "..", "test", "sets")
domain = load_domain(joinpath(path, "domain.pddl"))
problem = load_problem(joinpath(path, "problem.pddl"))

# Make sure function declarations have the right output type
@test PDDL.get_function(domain, :heard).type == :set
@test PDDL.get_function(domain, :known).type == :set

# State initialization should error if set functions are not registered yet
@test_throws Exception state = initstate(domain, problem)

# Register set theory
PDDL.Sets.register!()

state = initstate(domain, problem)
implementations = [
    "concrete interpreter" => (domain, state),
    "ground interpreter" => (ground(domain, state), state),
    "concrete compiler" => compiled(domain, state),
]

@testset "set fluents ($name)" for (name, (domain, _)) in implementations
    # Initialize state, test set membership and goal
    state = initstate(domain, problem)
    @test domain[state => pddl"(member (heard hanau) rumpelstiltskin)"] == 1
    @test domain[state => pddl"(member (heard steinau) cinderella)"] == 1
    @test domain[state => pddl"(subset (known jacob) story-set)"] == true
    @test domain[state => pddl"(subset (known wilhelm) (heard steinau))"] == false
    @test satisfy(domain, state, problem.goal) == false

    # Jacob tells stories at Steinau, Wilhem at Hanau
    state = execute(domain, state, pddl"(entertain jacob steinau)")
    state = execute(domain, state, pddl"(entertain wilhelm hanau)")
    @test domain[state => pddl"(cardinality (heard steinau))"] == 3
    @test domain[state => pddl"(cardinality (heard hanau))"] == 3

    # Both tell stories at Marburg
    state = execute(domain, state, pddl"(entertain jacob marburg)")
    state = execute(domain, state, pddl"(entertain wilhelm marburg)")
    @test domain[state => pddl"(cardinality (heard marburg))"] == 4

    # Check that goal is achieved
    @test satisfy(domain, state, problem.goal) == true
end

# Deregister set theory
PDDL.Sets.deregister!()
end # set fluents


text:
module PubChemCrawler

using HTTP
using LightXML

export atomregex, parse_formula, get_cid, get_for_cids, query_substructure, query_substructure_pug

include("utils.jl")
include("query.jl")
include("pugxml.jl")

end


text:
################################################################################
#
#   ENERGY MANIFOLD PLOTTING
#
################################################################################

# function to plot an energy manifold (2D)
function plotEnergyManifold(
            em :: AEM
            ;
            new_figure :: Bool = true,
            figsize :: Tuple = (6,6),
            color :: Vector{<:Integer} = [100,120,255],
            plot_bz :: Bool = true,
            ms :: Real = 5,
            BZ_corners :: Bool = false,
            kwargs...
        ) where {LS,S<:AbstractSite{LS,2},L,B,UC<:AbstractUnitcell{S,B},HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # create a new figure
    if plot_bz && new_figure
        # plot the brillouin zone
        plotBrillouinZone(getBrillouinZone(unitcell(hamiltonian(em))), figsize=figsize, BZ_corners=BZ_corners)
        # get the figure
        fig = gcf()
    elseif new_figure
        # only new figure
        fig = figure(figsize=figsize)
    else
        fig = gcf()
    end

    ###########################
    #   PLOT POINTS
    ###########################

    # obtain the points
    k_points = kpoints(em)
    kx = [k[1] for k in k_points]
    ky = [k[2] for k in k_points]

    # scatter all points
    scatter(kx,ky, color=color./255, s=ms)

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # equal axis
    gca().set_aspect("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # return the figure object
    return fig
end

# function to plot an energy manifold (3D)
function plotEnergyManifold(
            em :: AEM
            ;
            new_figure :: Bool = true,
            figsize :: Tuple = (6,6),
            color :: Vector{<:Integer} = [100,120,255],
            plot_bz :: Bool = true,
            ms :: Real = 5,
            BZ_corners :: Bool = false,
            kwargs...
        ) where {LS,S<:AbstractSite{LS,3},L,B,UC<:AbstractUnitcell{S,B},HB,H<:AbstractHamiltonian{L,UC,HB}, AEM <: AbstractEnergyManifold{H}}

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # create a new figure
    if plot_bz && new_figure
        # plot the brillouin zone
        plotBrillouinZone(getBrillouinZone(unitcell(hamiltonian(em))), figsize=figsize, BZ_corners=BZ_corners)
        # get the figure
        fig = gcf()
    elseif new_figure
        # only new figure
        fig = figure(figsize=figsize)
    else
        fig = gcf()
    end

    ###########################
    #   PLOT POINTS
    ###########################

    # obtain the points
    k_points = kpoints(em)
    kx = [k[1] for k in k_points]
    ky = [k[2] for k in k_points]
    kz = [k[3] for k in k_points]

    # scatter all points
    scatter3D(kx,ky,kz, color=color./255, s=ms)

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # equal axis
    gca().set_aspect("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # return the figure object
    return fig
end

# export plotting function
export plotEnergyManifold

# pass unknown arguments directly to construction of band structure
function plotEnergyManifold(
        args...
        ;
        kwargs...
    )

    # create and plot a bandstructure
    plotEnergyManifold(getEnergyManifold(args...); kwargs...)
end


text:
using Test
import PartialRejectionSampling
const PRS = PartialRejectionSampling

# using Documenter
# doctest(PartialRejectionSampling, fix=true)

@testset "PartialRejectionSampling.jl" begin
    # Write your tests here.
end


text:
module SPH

using Piecewise, Plots

export differentiate

include("kernels.jl")

const PARTICLE_MASS = 1.0
const H_INITIAL = 0.1
const ρ_INITIAL = 0.5
const T_INITIAL = 1

# TODO: Generalize for n dimensions
struct Particle{S}
    x::S
    y::S
    ρ::S
    h::S
    vx::S
    vy::S
    T::S
    Particle(x::S, y::S, ρ::S, h::S, vx::S, vy::S, T::S) where S = new{S}(x, y, ρ, h, vx, vy, T)
end

Particle(x::T, y::T) where T = Particle(x, y, ρ_INITIAL * one(T), H_INITIAL)

Particle(x::T, y::T, ρ::T, h::T) where T = Particle(x, y, ρ, h, zero(T), zero(T), one(T))

particles_in_box(n, w=1, h=1) = [Particle(w*rand(), h*rand()) for i in 1:n]

function plot_particles(particles; kwargs...)
    xs = [p.x for p in particles]
    ys = [p.y for p in particles]
    hs = [p.h for p in particles]
    ρs = [p.ρ for p in particles]

    p = plot()
    #draw_circles!(xs, ys, hs, lw = 0.5, fillalpha = 0)
    scatter!(xs, ys; mz = ρs, markerstrokewidth=0, markersize = 20, kwargs...)
    #surface!(xs, ys, ρs, camera=(0, 90); kwargs...)
    return p
end

function draw_circles!(xs, ys, hs; kwargs...)
    θ = LinRange(0, 2π, 50)
    for (x, y, h) in zip(xs, ys, hs)
        circle_x = x .+ h*sin.(θ)
        circle_y = y .+ h*cos.(θ)
        plot!(circle_x, circle_y; seriestype = [:shape], kwargs...)
    end
end

const MAX_TIME_STEPS = 100

# TODO: Use sparse upper triangular matrix for this
compute_distances(particles) = [hypot(p.x - q.x, p.y - q.y) for p in particles, q in particles]

function update_variable(particles, kernel = M6; binary = false, η = 1, dimension = 2)
    distances = compute_distances(particles)
    n_particles = length(particles)

    #new_particles = Vector{Particle}(undef, n_particles)
    new_particles = copy(particles)

    w  = kernel[dimension]
    ∂w = differentiate(kernel[dimension])

    d = dimension
    m = PARTICLE_MASS

    #binary = true

    for i in 1:n_particles
        xᵢ, yᵢ  = particles[i].x, particles[i].y
        ρᵢ, hᵢ  = particles[i].ρ, particles[i].h

        @views rs = distances[:, i]
        res = 1
        niter = 1

        if binary
            maxh = 1000
            minh = 0
            hᵢ = (maxh + minh)/2
        end

        q = zeros(n_particles)
        Wᵢ = zeros(n_particles)
        ∂Wᵢ = zeros(n_particles)

        #while abs(res) > 1e-8 && niter < MAX_TIME_STEPS
        #    q = rs / hᵢ
        #    Wᵢ = w.(q) / hᵢ^d
        #    ∂Wᵢ = ∂w.(q) / hᵢ^d
        #
        #    f  = m*((η/hᵢ)^d - sum(Wᵢ))
        #
        #    if binary
        #        if signbit(f)
        #            maxh = hᵢ
        #        else
        #            minh = hᵢ
        #        end
        #        newh = (maxh + minh)/2
        #        res = newh - hᵢ
        #        hᵢ = newh
        #    else
        #        f′ = -m*(d * η^d / hᵢ^(1 + d) + sum(q/hᵢ .* ∂Wᵢ .+ d/hᵢ * Wᵢ))
        #        res = - f/f′
        #        hᵢ += res
        #    end
        #end

        #ρᵢ = m*(η/hᵢ)^d
        qs = rs / particles[i].h
        ρᵢ = sum(m)

        Wᵢ = w.(q) / hᵢ^d
        #    ∂Wᵢ = ∂w.(q) / hᵢ^d

        #dhdρ = -hᵢ/ρᵢ/d

        dWdh = q/h * ∂Wᵢ - d/h * Wᵢ

        Ωᵢ = 1 - dhdρ * m * sum(dWdh)

        Pᵢ = ρᵢ * particles[i].T

        dvdt = 0

        #@show hᵢ, ρᵢ

        new_particles[i] = Particle(xᵢ, yᵢ, ρᵢ, hᵢ, vx, vy, T)
    end
    return new_particles
end

function update(particles; dt = 0.01, dimension = 2, kernel = M6, η = 1)
    distances = compute_distances(particles)
    n_particles = length(particles)

    #new_particles = Vector{Particle}(undef, n_particles)
    new_particles = copy(particles)

    w  = kernel[dimension]
    ∂w = differentiate(kernel[dimension])

    d = dimension
    m = PARTICLE_MASS

    h = 0.1
    q = distances / h

    # Update density
    new_ρ = zeros(n_particles)

    for i = 1:n_particles
        ρᵢ = 0.0
        for j in 1:n_particles
            if i != j
                Wᵢⱼ = 1/h^d * w(q[i, j])
                ρᵢ += Wᵢⱼ
            end
        end

        new_ρ[i] = ρᵢ
    end

    for i in 1:n_particles
        ρᵢ = 0.0
        dvx_dt = 0.0
        dvy_dt = 0.0

        particle_i = particles[i]
        xᵢ = particle_i.x
        yᵢ = particle_i.y

        for j in 1:n_particles
            if i != j
                ∂Wᵢⱼ = 1/h^d * ∂w(q[i, j])

                particle_j = particles[j]
                xⱼ = particle_j.x
                yⱼ = particle_j.y

                ρᵢ = new_ρ[i]
                ρⱼ = new_ρ[j]

                pᵢ = ρᵢ * particle_i.T
                pⱼ = ρⱼ * particle_j.T

                term = ∂Wᵢⱼ/distances[i, j]/h

                ∇Wx = term * (xᵢ - xⱼ)
                ∇Wy = term * (yᵢ - yⱼ)

                term = (pᵢ/ρᵢ^2 + pⱼ/ρⱼ^2)

                dvx_dt -= term * ∇Wx
                dvy_dt -= term * ∇Wy

                #isnan(dvx_dt) && println("dvx_dt nan at $i, $j")
                #isnan(dvy_dt) && println("dvy_dt nan at $i, $j")
            end
        end
        ρᵢ = m*ρᵢ
        dvx_dt = m*dvx_dt
        dvy_dt = m*dvy_dt
        vx = particle_i.vx + dt * dvx_dt
        vy = particle_i.vy + dt * dvy_dt

        x = particle_i.x + vx*dt
        y = particle_i.y + vy*dt

        if x >= 1
            x -= 1
        elseif x <= 0
            x += 1
        elseif y > 1
            y -= 1
        elseif y <= 0
            y += 1
        end

        new_particles[i] = Particle(x, y, ρᵢ, h, vx, vy, particle_i.T)
    end

    return new_particles
end

end


text:
using Combinatorics

function degeneracy(cluster)
    idx = Dict()
    multiple = Dict()
    for (i, a) in enumerate(cluster)
        if get(multiple, a, nothing) == nothing
            multiple[a] = 0
            idx[a] = []
        end
        multiple[a] += 1
        push!(idx[a], i)
    end

    total_mul = 1
    for (k, v) in pairs(multiple)
        total_mul *= factorial(v)
    end

    permutedidx = Dict()
    for s in keys(idx)
        permutedidx[s] = permutations(idx[s])
    end

    allpermutations = []
    orgidx = 1:length(cluster) |> collect
    for combinations in Iterators.product(values(permutedidx)...)
        copied = copy(orgidx)
        for (s, combo) in zip(keys(permutedidx), combinations)
            copied[idx[s]] = combo
        end

        push!(allpermutations, copied)
    end

    # return multiple, total_mul, allpermutations
    return allpermutations
end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
using Test, MathOptInterface

const MOI  = MathOptInterface
const MOIB = MathOptInterface.Bridges
const MOIT = MathOptInterface.Test
const MOIU = MOI.Utilities

import Clp
const OPTIMIZER = Clp.Optimizer()
MOI.set(OPTIMIZER, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "Clp"
end

@testset "supports_default_copy_to" begin
    @test !MOIU.supports_allocate_load(OPTIMIZER, false)
    @test !MOIU.supports_allocate_load(OPTIMIZER, true)
    @test !MOIU.supports_default_copy_to(OPTIMIZER, false)
    @test !MOIU.supports_default_copy_to(OPTIMIZER, true)
end

const CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const CACHED = MOIU.CachingOptimizer(CACHE, OPTIMIZER)
MOI.set(CACHED, MOI.Silent(), true)

const BRIDGED = MOI.Bridges.full_bridge_optimizer(CACHED, Float64)
MOI.set(BRIDGED, MOI.Silent(), true)

const CONFIG = MOIT.TestConfig(
    dual_objective_value=false
)

@testset "basic_constraint_tests" begin
    MOIT.basic_constraint_tests(CACHED, CONFIG)
end

@testset "Unit Tests" begin
    MOIT.unittest(BRIDGED, CONFIG, [
        # Unsupported attributes
        "number_threads",  # not supported by Clp
        "time_limit_sec",  # Weird behaviour of Clp
        "solve_time",  # not supported by Clp
        # Tests that require integer variables
        "solve_integer_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3",
        "solve_objbound_edge_cases",
        # Tests that require quadratic objective / constraints
        "solve_qcp_edge_cases",
        "solve_qp_edge_cases",
        # Tests that require SOC
        "delete_soc_variables",
    ])
end

@testset "Linear tests" begin
    MOIT.contlineartest(BRIDGED, CONFIG, [
        # Queries an infeasibility certificate
        "linear8a",
        # linear12 test requires the InfeasibilityCertificate for variable
        # bounds. These are available through C++, but not via the C interface.
        "linear12",
        # partial_start requires VariablePrimalStart to be implemented by the
        # solver.
        "partial_start"
    ])
end

@testset "ModelLike" begin
    # solver = Clp.Optimizer(LogLevel = 0)
    @testset "nametest" begin
        MOIT.nametest(BRIDGED)
    end
    @testset "validtest" begin
        MOIT.validtest(BRIDGED)
    end
    @testset "emptytest" begin
        MOIT.emptytest(BRIDGED)
    end
    # @testset "copytest" begin
    #     solver2 = Clp.Optimizer(LogLevel = 0)
    #     MOIT.copytest(solver,solver2)
    # end
    # Clp returns C_NULL when queried for the infeasibility ray in this case.
    @testset "Inexistant unbounded ray" begin
        # o = Clp.Optimizer(LogLevel = 0)
        MOI.empty!(BRIDGED)

        x = MOI.add_variables(BRIDGED, 5)
        MOI.set(BRIDGED, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                   MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0))
        MOI.set(BRIDGED, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        MOI.optimize!(BRIDGED)
        status = MOI.get(BRIDGED, MOI.TerminationStatus())
        @test status == MOI.DUAL_INFEASIBLE
    end
end


text:
function SQ._collect(src, q::SQ.FilterNode)
    res = default(src)
    flds = fields(src)
    for (field, T) in zip(flds, map(eltype, eltypes(src, flds...)))
        res[field] = NullableVector{T}()
    end
    for h in q.helpers
        apply!(res, h, src)
    end
    return res
end

function apply!(res, h::SQ.FilterHelper, src)::Void
    f, arg_fields, = SQ.parts(h)
    row_itr = eachrow(src)
    _apply!(res, h, row_itr)
    return
end

function _apply!(res, h::SQ.FilterHelper, row_itr)::Void
    f, arg_fields = SQ.parts(h)
    cols = columns(res)
    idx = index(res)
    indices = [ idx[field] for field in arg_fields ]
    # TODO: Make sure following is safe (in case of error in push!ing)
    for whole_row in row_itr
        args = whole_row[indices]
        if hasnulls(args)
            continue
        elseif f(map(unsafe_get, args))
            for (j, v) in enumerate(whole_row)
                push!(cols[j], v)
            end
        else
            continue
        end
    end
    return
end


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 = "seatd"
version = v"0.5.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://git.sr.ht/~kennylevinsen/seatd/archive/$version.tar.gz", "274b56324fc81ca6002bc1cdd387668dee34a6e1063e5f3896805c3770948988"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/seatd-*/
atomic_patch ../patches/fix_crosscompile.patch
mkdir build
cd build
apk add meson scdoc

PKG_CONFIG_SYSROOT_DIR="" meson ../ --cross-file="${MESON_TARGET_TOOLCHAIN}"
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p ->Sys.islinux(p) || Sys.isfreebsd(p), supported_platforms(; experimental=true))

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

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
mutable struct ZipCodeStr
  val::AbstractString
end

"""
```julia
mutable struct ZipCodeCleaner
  name::AbstractString
  pattern::Regex
  cleanfun::Function
end
```

Holds a textual identifier for clarity purposes, a regex
pattern to match a potential ZIP code to, and a function to
apply when a match is found.
"""
struct ZipCodeCleaner
  name::AbstractString
  pattern::Regex
  cleanfun::Function
end

"Returns a ZipCodeCleaner for stripping whitespace"
function removewhitespace()::ZipCodeCleaner
  ZipCodeCleaner("Whitespace", r"[ ]*", x -> strip(x))
end

"Returns a ZipCodeCleaner for removing the suffix on ZIP+4 patterns"
function removezipsuffix()::ZipCodeCleaner
  ZipCodeCleaner("Z+4 Suffix", r"-[0-9]{1,4}", x -> string(split(x,"-")[1]))
end

"Returns a ZipCodeCleaner for padding left zeros to ensure ZIP code is 5 digits"
function padleftzeros()::ZipCodeCleaner
  ZipCodeCleaner("Left Zeros", r"^[0-9]{1,4}$", x -> lpad(x, 5, "0"))
end

"""
```julia
cleanzipcode(
  Zip::ZipCodeStr,
  cleaner::ZipCodeCleaner
)
```

Executes `cleaner.cleanfun` on `Zip` when `Zip` matches `cleaner.pattern`.
"""
function cleanzipcode!(Zip::ZipCodeStr, cleaner::ZipCodeCleaner)::AbstractString
  Zip.val = ismatch(cleaner.pattern, Zip.val) ? cleaner.cleanfun(Zip.val) : Zip.val
end

"""
```julia
cleanzipcode(
    Zip::AbstractString;
    whitespace::Bool=true,      # Remove leading/trailing whitespace
    suffix::Bool=true,          # Remove "Z-{4}" suffix
    padzeros::Bool=true,        # Pad left zeros
    returnNA::Bool=true,        # Return unrecognized values as NA
    enforcestring::Bool=false   # Ensure return values are Strings
)
```

Description
===========
Returns `Zip` as an AbstractString with the following default corrections:
* Leading and trailing whitespace removed
* Z+4 suffix removed (e.g. 99999-1234 → 99999)
* Padded for left zeros(e.g. 245 → 00245)

If `Zip` does not match the "12345" pattern even after the above
corrections an NA value will be returned by default.

Usage
=====
```julia
cleanzipcode("01234 ")
cleanzipcode("1234")
cleanzipcode("1234-0000")
cleanzipcode("notaZip")
cleanzipcode("notaZip", enforcestring=true)
```
"""
function cleanzipcode(
            Zip::AbstractString;
            whitespace::Bool=true,      # Remove leading/trailing whitespace
            suffix::Bool=true,          # Remove "Z-{4}" suffix
            padzeros::Bool=true,        # Pad left zeros
            returnNA::Bool=true,        # Return unrecognized values as NA
            enforcestring::Bool=false   # Ensure return type is a String
          )::StringNA
  if ismatch(PATTERN_ZIPCODE, Zip)
    return Zip
  else
    # Clean routine
    Zip = ZipCodeStr(Zip)
    if whitespace cleanzipcode!(Zip, removewhitespace()) end
    if suffix cleanzipcode!(Zip, removezipsuffix()) end
    if padzeros cleanzipcode!(Zip, padleftzeros()) end

    # Return routine
    if ismatch(PATTERN_ZIPCODE, Zip.val)
      return Zip.val
    else
      return returnNA ? enforcestring ? "NA" : NA : Zip.val
    end

  end
end

function cleanzipcode(Zip::Any; args...)::StringNA
  cleanzipcode(string.(Zip); args...)
end

function cleanzipcode(Zip::AbstractVector; args...)::Vector{StringNA}
  [cleanzipcode(x) for x in Zip]
end


text:
# Automatically generated using Clang.jl

export bpConverterTypesC_Float
export bpConverterTypesC_UInt8, bpConverterTypesC_UInt16, bpConverterTypesC_UInt32, bpConverterTypesC_UInt64
const bpConverterTypesC_UInt8 = Cuchar # UInt8
const bpConverterTypesC_UInt16 = UInt16
const bpConverterTypesC_UInt32 = UInt32
const bpConverterTypesC_UInt64 = Culonglong # UInt64
const bpConverterTypesC_Float = Cfloat # Float32

const bpImageConverterC = Cvoid # Nothing
const bpImageConverterCPtr = Ptr{bpImageConverterC} # Ptr{Nothing}
const bpConverterTypesC_ProgressCallback = Ptr{Cvoid} # Ptr{Nothing}

export bpConverterTypesC_DataType
export bpConverterTypesC_UInt8Type, bpConverterTypesC_UInt16Type, bpConverterTypesC_UInt32Type, bpConverterTypesC_FloatType
@enum bpConverterTypesC_DataType::UInt32 begin
    bpConverterTypesC_UInt8Type = 0
    bpConverterTypesC_UInt16Type = 1
    bpConverterTypesC_UInt32Type = 2
    bpConverterTypesC_FloatType = 3
end
Base.convert(::Type{bpConverterTypesC_DataType}, ::Type{UInt8}) = bpConverterTypesC_UInt8Type
Base.convert(::Type{bpConverterTypesC_DataType}, ::Type{UInt16}) = bpConverterTypesC_UInt16Type
Base.convert(::Type{bpConverterTypesC_DataType}, ::Type{UInt32}) = bpConverterTypesC_UInt32Type
Base.convert(::Type{bpConverterTypesC_DataType}, ::Type{Float32}) = bpConverterTypesC_Float32Type

export bpConverterTypesC_ImageExtent
struct bpConverterTypesC_ImageExtent
    mExtentMinX::Cfloat
    mExtentMinY::Cfloat
    mExtentMinZ::Cfloat
    mExtentMaxX::Cfloat
    mExtentMaxY::Cfloat
    mExtentMaxZ::Cfloat
end
bpConverterTypesC_ImageExtent(; X = (0,1), Y = (0,1), Z = (0,1)) =
    bpConverterTypesC_ImageExtent(X[1], Y[1], Z[1], X[2], Y[2], Z[2])

const ImageExtent = bpConverterTypesC_ImageExtent

const bpConverterTypesC_ImageExtentPtr = Ptr{bpConverterTypesC_ImageExtent}
export bpConverterTypesC_ImageExtentPtr

export bpConverterTypesC_Dimension
export bpConverterTypesC_DimensionX, bpConverterTypesC_DimensionY, bpConverterTypesC_DimensionZ, bpConverterTypesC_DimensionC, bpConverterTypesC_DimensionT
@enum bpConverterTypesC_Dimension::UInt32 begin
    bpConverterTypesC_DimensionX = 0
    bpConverterTypesC_DimensionY = 1
    bpConverterTypesC_DimensionZ = 2
    bpConverterTypesC_DimensionC = 3
    bpConverterTypesC_DimensionT = 4
end
function bpConverterTypesC_Dimension(s::Symbol)
    if s == :X
        bpConverterTypesC_DimensionX
    elseif s == :Y
        bpConverterTypesC_DimensionY
    elseif s == :Z
        bpConverterTypesC_DimensionZ
    elseif s == :C
        bpConverterTypesC_DimensionC
    elseif s == :T
        bpConverterTypesC_DimensionT
    else
        error("Symbol cannot be mapped to a known dimension")
    end
end
bpConverterTypesC_Dimension(str::String) = bpConverterTypesC_Dimension(Symbol(str))
Base.convert(::Type{bpConverterTypesC_Dimension},s::Union{Symbol,String}) = bpConverterTypesC_Dimension(s)

export bpConverterTypesC_DimensionSequence5D
struct bpConverterTypesC_DimensionSequence5D
    mDimension0::bpConverterTypesC_Dimension
    mDimension1::bpConverterTypesC_Dimension
    mDimension2::bpConverterTypesC_Dimension
    mDimension3::bpConverterTypesC_Dimension
    mDimension4::bpConverterTypesC_Dimension
end
const DimensionSequence5D = bpConverterTypesC_DimensionSequence5D

const bpConverterTypesC_DimensionSequence5DPtr = Ptr{bpConverterTypesC_DimensionSequence5D}

export bpConverterTypesC_Size5D
struct bpConverterTypesC_Size5D
    mValueX::UInt32
    mValueY::UInt32
    mValueZ::UInt32
    mValueC::UInt32
    mValueT::UInt32
    bpConverterTypesC_Size5D(mValueX, mValueY = 1, mValueZ = 1, mValueC = 1, mValueT = 1) = new(mValueX, mValueY, mValueZ, mValueC, mValueT)
end
const Size5D = bpConverterTypesC_Size5D
bpConverterTypesC_Size5D(; X = 1, Y = 1, Z = 1, C = 1, T = 1) = bpConverterTypesC_Size5D(X, Y, Z, C, T)

const bpConverterTypesC_Size5DPtr = Ptr{bpConverterTypesC_Size5D}

export bpConverterTypesC_Index5D
mutable struct bpConverterTypesC_Index5D
    mValueX::UInt32
    mValueY::UInt32
    mValueZ::UInt32
    mValueC::UInt32
    mValueT::UInt32
    bpConverterTypesC_Index5D(mValueX, mValueY = 0, mValueZ = 0, mValueC = 0, mValueT = 0) = new(mValueX, mValueY, mValueZ, mValueC, mValueT)
end
const Index5D = bpConverterTypesC_Index5D

bpConverterTypesC_Index5D(; X = 0, Y = 0, Z = 0, C = 0, T =0) = bpConverterTypesC_Index5D(X, Y, Z, C, T)

const bpConverterTypesC_Index5DPtr = Ptr{bpConverterTypesC_Index5D}
const bpConverterTypesC_String = Cstring

export tCompressionAlgorithmType
@enum tCompressionAlgorithmType::UInt32 begin
    eCompressionAlgorithmNone = 0
    eCompressionAlgorithmGzipLevel1 = 1
    eCompressionAlgorithmGzipLevel2 = 2
    eCompressionAlgorithmGzipLevel3 = 3
    eCompressionAlgorithmGzipLevel4 = 4
    eCompressionAlgorithmGzipLevel5 = 5
    eCompressionAlgorithmGzipLevel6 = 6
    eCompressionAlgorithmGzipLevel7 = 7
    eCompressionAlgorithmGzipLevel8 = 8
    eCompressionAlgorithmGzipLevel9 = 9
    eCompressionAlgorithmShuffleGzipLevel1 = 11
    eCompressionAlgorithmShuffleGzipLevel2 = 12
    eCompressionAlgorithmShuffleGzipLevel3 = 13
    eCompressionAlgorithmShuffleGzipLevel4 = 14
    eCompressionAlgorithmShuffleGzipLevel5 = 15
    eCompressionAlgorithmShuffleGzipLevel6 = 16
    eCompressionAlgorithmShuffleGzipLevel7 = 17
    eCompressionAlgorithmShuffleGzipLevel8 = 18
    eCompressionAlgorithmShuffleGzipLevel9 = 19
    eCompressionAlgorithmLZ4 = 21
    eCompressionAlgorithmLShuffleLZ4 = 31
end
for i in instances(tCompressionAlgorithmType)
    @eval export $(Symbol(i))
end

export bpConverterTypesC_Options
struct bpConverterTypesC_Options
    mThumbnailSizeXY::UInt32
    mFlipDimensionX::Bool
    mFlipDimensionY::Bool
    mFlipDimensionZ::Bool
    mForceFileBlockSizeZ1::Bool
    mEnableLogProgress::Bool
    mNumberOfThreads::UInt32
    mCompressionAlgorithmType::tCompressionAlgorithmType
    mDisablePyramid::Bool
    bpConverterTypesC_Options(mThumbnailSizeXY, mFlipDimensionX, mFlipDimensionY, mFlipDimensionZ, mForceFileBlockSizeZ1, mEnableLogProgress, mNumberOfThreads, mCompressionAlgorithmType, mDisablePyramid = false) =
        new(mThumbnailSizeXY, mFlipDimensionX, mFlipDimensionY, mFlipDimensionZ, mForceFileBlockSizeZ1, mEnableLogProgress, mNumberOfThreads, mCompressionAlgorithmType, mDisablePyramid)
end
const Options = bpConverterTypesC_Options

# All keyword version with reasonable defaults
function bpConverterTypesC_Options(;
    mThumbnailSizeXY = 256,
    mFlipDimensionX = false,
    mFlipDimensionY = false,
    mFlipDimensionZ = false,
    mForceFileBlockSizeZ1 = false,
    mEnableLogProgress = true,
    mNumberOfThreads = Threads.nthreads(),
    mCompressionAlgorithmType = eCompressionAlgorithmLShuffleLZ4,
    mDisablePyramid = false
    )
    Options(
        mThumbnailSizeXY,
        mFlipDimensionX,
        mFlipDimensionY,
        mFlipDimensionZ,
        mForceFileBlockSizeZ1,
        mEnableLogProgress,
        mNumberOfThreads,
        mCompressionAlgorithmType,
        mDisablePyramid
    )
end

const bpConverterTypesC_OptionsPtr = Ptr{bpConverterTypesC_Options}

export bpConverterTypesC_Parameter
struct bpConverterTypesC_Parameter
    mName::bpConverterTypesC_String
    mValue::bpConverterTypesC_String
end
function bpConverterTypesC_Parameter(mName::String, mValue::String)
    bpConverterTypesC_Parameter(Cstring(pointer(mName)),
                                Cstring(pointer(mValue)))
end

struct Parameter
    mName::String
    mValue::String
end
bpConverterTypesC_Parameter(p::Parameter) = bpConverterTypesC_Parameter(p.mName, p.mValue)
Base.unsafe_convert(::Type{bpConverterTypesC_Parameter}, p::Parameter) =
    bpConverterTypesC_Parameter(p)

export bpConverterTypesC_ParameterSection
struct bpConverterTypesC_ParameterSection
    mName::bpConverterTypesC_String
    mValues::Ptr{bpConverterTypesC_Parameter}
    mValuesCount::UInt32
end

struct ParameterSection
    mName::String
    mValues::Vector{bpConverterTypesC_Parameter}
end
ParameterSection(mName::String, mValues::Vector{Parameter}) =
    ParameterSection(mName, bpConverterTypesC_Parameter.(mValues) )

ParameterSection(mName::String, p::Parameter) =
    ParameterSection(mName, Parameter[p])

Base.unsafe_convert(::Type{bpConverterTypesC_ParameterSection}, ps::ParameterSection) =
    bpConverterTypesC_ParameterSection(ps)

function bpConverterTypesC_ParameterSection(ps::ParameterSection)
    bpConverterTypesC_ParameterSection(
        Base.unsafe_convert(Cstring, ps.mName),
        pointer( ps.mValues ),
        length( ps.mValues )
    )
end

export bpConverterTypesC_Parameters
struct bpConverterTypesC_Parameters
    mValues::Ptr{bpConverterTypesC_ParameterSection}
    mValuesCount::UInt32
end
const Parameters = Vector{bpConverterTypesC_ParameterSection}

const bpConverterTypesC_ParametersPtr = Ptr{bpConverterTypesC_Parameters}

export bpConverterTypesC_TimeInfo
struct bpConverterTypesC_TimeInfo
    mJulianDay::UInt32
    mNanosecondsOfDay::Culonglong
end
const TimeInfo = bpConverterTypesC_TimeInfo

export bpConverterTypesC_TimeInfos
struct bpConverterTypesC_TimeInfos
    mValues::Ptr{bpConverterTypesC_TimeInfo}
    mValuesCount::UInt32
end
const TimeInfos = Vector{TimeInfo}
Base.unsafe_convert(::Type{bpConverterTypesC_TimeInfos}, v::TimeInfos) =
    bpConverterTypesC_TimeInfos(v)

bpConverterTypesC_TimeInfos(v::TimeInfos) =
    bpConverterTypesC_TimeInfos(
        pointer( v ),
        length( v )
    )

const bpConverterTypesC_TimeInfoVector = Ptr{bpConverterTypesC_TimeInfos}

export bpConverterTypesC_Color
struct bpConverterTypesC_Color
    mRed::Cfloat
    mGreen::Cfloat
    mBlue::Cfloat
    mAlpha::Cfloat
end
const Color = bpConverterTypesC_Color

export bpConverterTypesC_ColorInfo
struct bpConverterTypesC_ColorInfo
    mIsBaseColorMode::Bool
    mBaseColor::bpConverterTypesC_Color
    mColorTable::Ptr{bpConverterTypesC_Color}
    mColorTableSize::UInt32
    mOpacity::Cfloat
    mRangeMin::Cfloat
    mRangeMax::Cfloat
    mGammaCorrection::Cfloat
end

struct ColorInfo
    mIsBaseColorMode::Bool
    mBaseColor::bpConverterTypesC_Color
    mColorTable::Vector{bpConverterTypesC_Color}
    mOpacity::Cfloat
    mRangeMin::Cfloat
    mRangeMax::Cfloat
    mGammaCorrection::Cfloat
end

ColorInfo(
    mColorTable::Vector{bpConverterTypesC_Color};
    mOpacity = 0,
    mRangeMin = 0,
    mRangeMax = 255,
    mGammaCorrection = 1
) = ColorInfo(
    false,
    bpConverterTypesC_Color(0,0,0,0),
    mColorTable,
    mOpacity, mRangeMin, mRangeMax, mGammaCorrection)

ColorInfo(
    mBaseColor::bpConverterTypesC_Color;
    mOpacity = 0,
    mRangeMin = 0,
    mRangeMax = 255,
    mGammaCorrection = 1
) = ColorInfo(
    true,
    mBaseColor,
    bpConverterTypesC_Color[],
    mOpacity, mRangeMin, mRangeMax, mGammaCorrection
)

Base.unsafe_convert(::Type{bpConverterTypesC_ColorInfo}, ci::ColorInfo) =
    bpConverterTypesC_ColorInfo(ci)

bpConverterTypesC_ColorInfo(ci::ColorInfo) =
    bpConverterTypesC_ColorInfo(
        ci.mIsBaseColorMode,
        ci.mBaseColor,
        pointer( ci.mColorTable ),
        length( ci.mColorTable ),
        ci.mOpacity,
        ci.mRangeMin,
        ci.mRangeMax,
        ci.mGammaCorrection
    )

export bpConverterTypesC_ColorInfos
struct bpConverterTypesC_ColorInfos
    mValues::Ptr{bpConverterTypesC_ColorInfo}
    mValuesCount::UInt32
end
const bpConverterTypesC_ColorInfoVector = Ptr{bpConverterTypesC_ColorInfos}

const ColorInfos = Vector{bpConverterTypesC_ColorInfos}
function Base.unsafe_convert(::Type{bpConverterTypesC_ColorInfos}, cis::ColorInfos)
    bpConverterTypesC_ColorInfos(
        pointer( cis ),
        length( cis )
    )
end

# Skipping MacroDefinition: BP_IMARISWRITER_DLL_EXPORT __declspec ( dllexport )
# Skipping MacroDefinition: BP_IMARISWRITER_DLL_IMPORT __declspec ( dllimport )
# const BP_IMARISWRITER_DLL_API = BP_IMARISWRITER_DLL_IMPORT