text:
using DelimitedFiles
using Logging
using Random
# using StatsBase
using Statistics

"""
    DataSplit

A basic struct for encapsulating the four components of supervised training.
"""
mutable struct DataSplit
    train_x::Array
    test_x::Array
    train_y::Array
    test_y::Array
    DataSplit(train_x, test_x, train_y, test_y) = new(train_x, test_x, train_y, test_y)
end

"""
    DataSplit(data_x::Array, data_y::Array, ratio::Float)

Return a DataSplit struct that is split by the ratio (e.g. 0.8).
"""
function DataSplit(data_x::Array, data_y::Array, ratio::Real)
    dim, n_data = size(data_x)
    split_ind = Int(floor(n_data*ratio))

    train_x = data_x[:, 1:split_ind]
    test_x = data_x[:, split_ind+1:end]
    train_y = data_y[1:split_ind]
    test_y = data_y[split_ind+1:end]

    return DataSplit(train_x, test_x, train_y, test_y)
end

"""
    DataSplit(data_x::Array, data_y::Array, ratio::Real, seq_ind::Array)

Sequential loading and ratio split of the data.
"""
function DataSplit(data_x::Array, data_y::Array, ratio::Real, seq_ind::Array)
    dim, n_data = size(data_x)
    n_splits = length(seq_ind)

    train_x = Array{Float64}(undef, dim, 0)
    train_y = Array{Float64}(undef, 0)
    test_x = Array{Float64}(undef, dim, 0)
    test_y = Array{Float64}(undef, 0)

    # Iterate over all splits
    for ind in seq_ind
        local_x = data_x[:, ind[1]:ind[2]]
        local_y = data_y[ind[1]:ind[2]]
        # n_data = ind[2] - ind[1] + 1
        n_data = size(local_x)[2]
        split_ind = Int(floor(n_data*ratio))

        train_x = [train_x local_x[:, 1:split_ind]]
        test_x = [test_x local_x[:, split_ind+1:end]]
        train_y = [train_y; local_y[1:split_ind]]
        test_y = [test_y; local_y[split_ind+1:end]]
    end
    return DataSplit(train_x, test_x, train_y, test_y)
end # DataSplit(data_x::Array, data_y::Array, ratio::Real, seq_ind::Array)

"""
    load_iris(data_path::String ; split_ratio::Real = 0.8)

Loads the iris dataset for testing and examples.
"""
function load_iris(data_path::String ; split_ratio::Real = 0.8)
    raw_data = readdlm(data_path,',')
    labels = ["Iris-setosa", "Iris-versicolor", "Iris-virginica"]
    raw_x = convert(Array{Float64, 2}, raw_data[2:end, 2:5])
    raw_y_labels = raw_data[2:end, 6]
    raw_y = Array{Integer}(undef, 0)
    for ix = 1:length(raw_y_labels)
        for jx = 1:length(labels)
            if raw_y_labels[ix] == labels[jx]
                push!(raw_y, jx)
            end
        end
    end
    n_samples, n_features = size(raw_x)

    # Julia is column-major, so use columns for features
    raw_x = permutedims(raw_x)

    # Shuffle the data and targets
    ind_shuffle = Random.randperm(n_samples)
    x = raw_x[:, ind_shuffle]
    y = raw_y[ind_shuffle]

    data = DataSplit(x, y, split_ratio)

    return data
end # load_iris(data_path::String ; split_ratio::Real = 0.8)

"""
    get_cvi_data(data_file::String)

Get the cvi data specified by the data_file path.
"""
function get_cvi_data(data_file::String)
    # Parse the data
    data = readdlm(data_file, ',')
    data = permutedims(data)
    train_x = convert(Matrix{Float64}, data[1:2, :])
    train_y = convert(Vector{Integer}, data[3, :])

    return train_x, train_y
end # get_cvi_data(data_file::String)

ext: jl


text:
module initialiser
using population, _deme_, core, fitness, representation, individual
export Initialiser, initialise!

"""
The base type used by all initialisers.
"""
abstract Initialiser

"""
Initialises all members within a given population.
"""
initialise!(init::Initialiser, p::Population) = for d in p.demes
  initialise!(init, d)
end

"""
Initialises all members within a given deme.
"""
initialise!(init::Initialiser, d::Deme) = for i in 1:d.capacity
  initialise!(init, i)
end

"""
The default initialiser used by all evolutionary algorithms.
"""
type DefaultInitialiser <: Initialiser; end

"""
Initialises a given deme by generating its contents at pseudo-random according
to its default generator method.
"""
function initialise!(i::DefaultInitialiser, d::Deme)
  rep = d.species.genotype.representation
  chromo = chromosome(rep)
  d.members.fitnesses = d.offspring.fitnesses = Array(uses(d.species.fitness), d.capacity)
  d.members.stages[d.species.genotype.label] =
    d.offspring.stages[d.species.genotype.label] =
    IndividualStage{chromo}[IndividualStage{chromo}(generate(rep)) for i in 1:d.capacity]
end

end


text:
using SumBassSolver
using Plots

p=fill(0.001,2)
q=[0.1,0.11]
m=fill(100,2)
c0=fill(0,2)
alpha = fill(1,2)
F = fill(0,2)
T = fill(100,2)
Pi=90

sumBassSolver(p,q,m,c0,alpha,F,T,Pi)
sumBassSolver(p,q,m,c0,alpha,F,T,Pi,localopt=true)
sumBassSolver(p,q,m,c0,alpha,F,T,Pi,alg="SB")
sumBassSolver(p,q,m,c0,alpha,F,T,Pi,alg="SB",localopt=true)

Pis = 50:100
v = Vector{Float64}(undef,length(Pis))
x = Array{Float64}(undef,2,length(Pis))

@time for i in 1:length(Pis)
    x[:,i],v[i] = sumBassSolver(p,q,m,c0,alpha,F,T,Pis[i],localopt=true)
end

pl1 = plot(Pis,v);
pl2 = plot(Pis,x[1,:],label = "plot1");
plot!(Pis,x[2,:],label = "plot2");
plot!(Pis,fill(Bass_inflection(p[1],q[1],m[1],c0[1]),length(Pis)),label="inflection1");
plot!(Pis,fill(Bass_inflection(p[2],q[2],m[2],c0[2]),length(Pis)),label="inflection2");
pl2

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
function dimerdecay(;x = [1_000, 0, 0], c = [1.0, 0.002, 0.5, 0.04])
  m = Network("dimer-decay / Rathiam, Petzold, Cao, Gillespie 2003")

  m <= Species("S1", x[1])
  m <= Species("S2", x[2])
  m <= Species("S3", x[3])

  m <= Reaction("R1", c[1], "S1 --> 0")
  m <= Reaction("R2", c[2], "2 * S1 --> S2")
  m <= Reaction("R3", c[3], "S2 --> 2 * S1")
  m <= Reaction("R4", c[4], "S2 --> S3")

  return m
end

model = dimerdecay()


text:
"""
    macro sk_reg(expr)

Helper macro for defining interfaces of ScikitLearn regression models.
Struct fields require a type annotation and a default value as in the example
below. Constraints for parameters (fields) are introduced as for `field3`.
The constraint must refer to the parameter as `arg` or `_`. If the used
parameter does not meet a constraint the default value will be used.

@sk_reg mutable struct SomeRegression <: MMI.Deterministic
    field1::Int = 1
    field2::Any = nothing
    field3::Float64 = 0.5::(0 < arg < 0.8)
end

MMI.fit and MMI.predict methods are also produced. See also [`@sk_clf`](@ref)
"""
macro sk_reg(expr)
    modelname, params, clean_expr, expr = _sk_constructor(expr)
    fit_expr = _skmodel_fit_reg(modelname, params)
    _sk_finalize(modelname, clean_expr, fit_expr, expr)
end

"""
    macro sk_clf(expr)

Same as [`@sk_reg`](@ref) but for classifiers.
"""
macro sk_clf(expr)
    modelname, params, clean_expr, expr = _sk_constructor(expr)
    fit_expr = _skmodel_fit_clf(modelname, params)
    _sk_finalize(modelname, clean_expr, fit_expr, expr)
end

"""
    macro sk_uns(expr)

Same as [`@sk_reg`](@ref) but for unsupervised models.
"""
macro sk_uns(expr)
    modelname, params, clean_expr, expr = _sk_constructor(expr)
    fit_expr = _skmodel_fit_uns(modelname, params)
    _sk_finalize(modelname, clean_expr, fit_expr, expr)
end

# =================
# HELPER FUNCTIONS
# =================

"""
    _sk_constructor(expr)

Extracts the relevant information from the expr and build the expression
corresponding to the model constructor (see [`_model_constructor`](@ref)).
"""
function _sk_constructor(expr)
    # similar to @mlj_model
    expr, modelname, params, defaults, constraints = _process_model_def(@__MODULE__, expr)
    # keyword constructor
    const_expr = _model_constructor(modelname, params, defaults)
    # associate the constructor with the definition of the struct
    push!(expr.args[3].args, const_expr)
    # cleaner
    clean_expr = _model_cleaner(modelname, defaults, constraints)
    # return
    return modelname, params, clean_expr, expr
end

"""
    _sk_finalize(expr)

Call the different bits and pieces of code and assemble the final expression
which will:
* generate a constructor
* generate a `clean!` method
* generate the fit/predict/transform method as appropriate
* add traits.
"""
function _sk_finalize(m, clean_expr, fit_expr, expr)
    # call a different predict based on whether the model is  probabilistic,
    # deteterministic or unsupervised
    if expr.args[2].args[2] == :(MMI.Deterministic)
        predict_expr = _skmodel_predict(m)
    elseif expr.args[2].args[2] == :(MMI.Probabilistic)
        predict_expr = _skmodel_predict_prob(m)
    else
        predict_expr = nothing
    end
    expr = quote
        export $m       # make the model available
        $expr           # effectively the constructor
        $clean_expr     # clean! method
        $fit_expr       # fit method
        $predict_expr   # predict method if relevant
        MMI.load_path(::Type{<:$m})       = "$API_PKG_NAME.$(MMI.name($m))"
        MMI.package_name(::Type{<:$m})    = PKG_NAME
        MMI.is_pure_julia(::Type{<:$m})   = false
        MMI.package_license(::Type{<:$m}) = SK_LIC
        MMI.package_uuid(::Type{<:$m})    = SK_UUID
        MMI.package_url(::Type{<:$m})     = SK_URL
        MMI.is_wrapper(::Type{<:$m})      = false
    end
    esc(expr)
end

# =================================
# Specifics for SUPERVISED MODELS
# =================================

"""
    _skmodel_fit_reg(modelname, params)

Called as part of [`@sk_reg`](@ref), returns the expression corresponing to the
`fit` method for a ScikitLearn regression model.
"""
function _skmodel_fit_reg(modelname, params)
    expr = quote
        function MMI.fit(model::$modelname, verbosity::Int, X, y)
            # set X and y into a format that can be processed by sklearn
            Xmatrix   = MMI.matrix(X)
            yplain    = y
            targnames = nothing
            # check if it's a multi-target regression case, in that case keep
            # track of the names of the target columns so that the prediction
            # can be named accordingly
            if MMI.istable(y)
               yplain    = MMI.matrix(y)
               targnames = MMI.schema(y).names
            end
            # --------------------------------------------------------------
            # NOTE: be very careful before modifying the next few lines,
            # in particular, consider the following recommendations:
            # JuliaPy/PyCall.jl#using-pycall-from-julia-modules
            # --------------------------------------------------------------
            # sksym is the namespace (e.g. :SKLM)
            # skmod is the submodule (e.g. :linear_model)
            # mdl is the actual model (e.g. :ARDRegression)
            sksym, skmod, mdl = $(Symbol(modelname, "_"))
            # retrieve the namespace, if it's not there yet, import it
            parent = eval(sksym)
            ispynull(parent) && ski!(parent, skmod)
            # retrieve the effective ScikitLearn constructor
            skconstr = getproperty(parent, mdl)
            # build the scikitlearn model passing all the parameters
            skmodel = skconstr(
                        $((Expr(:kw, p, :(model.$p)) for p in params)...))
            # --------------------------------------------------------------
            # fit and organise results
            fitres = SK.fit!(skmodel, Xmatrix, yplain)
            # TODO: we may want to use the report later on
            report = NamedTuple()
            # the first nothing is so that we can use the same predict for
            # regressors and classifiers
            return ((fitres, nothing, targnames), nothing, report)
        end
    end
end

"""
    _skmodel_fit_clf(modelname, params)

Called as part of [`@sk_clf`](@ref), returns the expression corresponing to the
`fit` method for a ScikitLearn classifier model.
"""
function _skmodel_fit_clf(modelname, params)
    quote
        function MMI.fit(model::$modelname, verbosity::Int, X, y)
            Xmatrix = MMI.matrix(X)
            yplain  = MMI.int(y)
            # See _skmodel_fit_reg, same story
            sksym, skmod, mdl = $(Symbol(modelname, "_"))
            parent = eval(sksym)
            ispynull(parent) && ski!(parent, skmod)
            skconstr = getproperty(parent, mdl)
            skmodel = skconstr(
                        $((Expr(:kw, p, :(model.$p)) for p in params)...))
            fitres  = SK.fit!(skmodel, Xmatrix, yplain)
            # TODO: we may want to use the report later on
            report  = NamedTuple()
            # pass y[1] for decoding in predict method, first nothing
            # is targnames
            return ((fitres, y[1], nothing), nothing, report)
        end
    end
end

"""
    _skmodel_predict(modelname)

Called as part of [`@sk_model`](@ref), returns the expression corresponing to
the `predict` method for the ScikitLearn model (for a deterministic model).
"""
function _skmodel_predict(modelname)
    quote
        function MMI.predict(model::$modelname, (fitres, y1, targnames), Xnew)
            Xmatrix = MMI.matrix(Xnew)
            preds   = SK.predict(fitres, Xmatrix)
            if isa(preds, Matrix)
                # only regressors are possibly multitarget;
                # build a table with the appropriate column names
                return MMI.table(preds, names=targnames)
            end
            if y1 !== nothing
                # if it's a classifier
                return preds |> MMI.decoder(y1)
            end
            return preds
        end
    end
end

"""
    _skmodel_predict_prob(modelname)

Same as `_skmodel_predict` but with probabilities. Note that only classifiers
are probabilistic in sklearn so that we always decode.
"""
function _skmodel_predict_prob(modelname)
    quote
        # there are no multi-task classifiers in sklearn
        function MMI.predict(model::$modelname, (fitres, y1, _), Xnew)
            Xmatrix = MMI.matrix(Xnew)
            # this is an array of size n x c with rows that sum to 1
            preds   = SK.predict_proba(fitres, Xmatrix)
            classes = MMI.classes(y1)
            return MMI.UnivariateFinite(classes, preds)
        end
    end
end

# ============================================================================
# Specifics for UNSUPERVISED MODELS
# ============================================================================
# Depending on the model there may be
# * a transform
# * a inverse_transform
# * a predict

"""
    _skmodel_fit_uns(modelname, params)

Called as part of [`@sk_uns`](@ref), returns the expression corresponing to the
`fit` method for a ScikitLearn unsupervised model.
"""
function _skmodel_fit_uns(modelname, params)
    quote
        function MMI.fit(model::$modelname, verbosity::Int, X)
            Xmatrix = MMI.matrix(X)
           # See _skmodel_fit_reg, same story
            sksym, skmod, mdl = $(Symbol(modelname, "_"))
            parent = eval(sksym)
            ispynull(parent) && ski!(parent, skmod)
            skconstr = getproperty(parent, mdl)
            skmodel = skconstr(
                        $((Expr(:kw, p, :(model.$p)) for p in params)...))
            fitres  = SK.fit!(skmodel, Xmatrix)
            # TODO: we may want to use the report later on
            report  = NamedTuple()
            return (fitres, nothing, report)
        end
    end
end

"""
    macro sku_transform(modelname)

Adds a `transform` method to a declared scikit unsupervised model if
there is one supported.
"""
macro sku_transform(modelname)
    quote
        function MMI.transform(::$modelname, fitres, X)
            X = SK.transform(fitres, MMI.matrix(X))
            MMI.table(X)
        end
    end
end

"""
    macro sku_inverse_transform(modelname)

Adds an `inverse_transform` method to a declared scikit unsupervised model if
there is one supported.
"""
macro sku_inverse_transform(modelname)
    quote
        function MMI.inverse_transform(::$modelname, fitres, X)
            X = SK.inverse_transform(fitres, MMI.matrix(X))
            MMI.table(X)
        end
    end
end

"""
    macro sku_predict

Adds a `predict` method to a declared scikit unsupervised model if there is one
supported. Returns a MMI.categorical vector.
Only `AffinityPropagation`, `Birch`, `KMeans`, `MiniBatchKMeans` and
`MeanShift` support a `predict` method.

Note: for models that offer a `fit_predict`, the encoding is done in the
`fitted_params`.
"""
macro sku_predict(modelname)
    quote
        function MMI.predict(m::$modelname, fitres, X)
            # this is due to the fact that we have nested modules
            # so potentially have to extract the leaf node...
            sm = Symbol($modelname)
            ss = string(sm)
            sm = Symbol(split(ss, ".")[end])
            if sm in (:Birch, :KMeans, :MiniBatchKMeans)
                catv = MMI.categorical(1:m.n_clusters)
            elseif sm == :AffinityPropagation
                nc   = length(fitres.cluster_centers_indices_)
                catv = MMI.categorical(1:nc)
            elseif sm == :MeanShift
                nc   = size(fitres.cluster_centers_, 1)
                catv = MMI.categorical(1:nc)
            else
                throw(ArgumentError("Model $sm does not support `predict`."))
            end
            preds  = broadcast(+, SK.predict(fitres, MMI.matrix(X)), 1)
            return catv[preds]
        end
    end
end


text:
```
Reg(args...)  -->Plot

Estimate Linear, Exponential and Power regression type
``Linear\quad \Rightarrow y=mx+b\\ Exponential\quad \Rightarrow y=b{ e }^{ mx }\\ Power\quad \Rightarrow y=b{ x }^{ m }``

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|x|Mandatory|--|Array{Number,1}| x Data|
|y|Mandatory|--|Array{Number,1}| y Data|
|TypeReg|Optional|TypeReg="Linear"|TypeReg= "linear" or "Exponential"  or "Power"|Type of regression|
|pr|Optional|pr=false|pr=Bool|Print the results of regression|
```
function Reg(x,y; TypeReg="Linear", pr=false)

    if TypeReg=="Linear"
        m=(mean(x)*mean(y)-mean(x.*y))/(mean(x)^2-mean(x.^2))
        b=mean(y)-m*mean(x)
        SELine=sum(map((X,Y)->((X*m+b)-Y)^2,x,y))
        SEYmean=sum(map(Y->(Y-mean(y))^2,y))
        R2=1-(SELine/SEYmean)

        pr ? println(" **Linear Regression \n form y = mx+b \n m=$m \n b=$b \n r2=$R2") : :none
        return m,b,R2
    elseif TypeReg=="Exponential"
        m=(mean(x)*mean(log.(y))-mean(x.*log.(y) ) )/(mean(x)^2-mean(x.^2))
        Logb=mean(log.(y))-m*mean(x)
        b=exp(Logb)
        SELine=sum(map((X,Y)->((X*m+Logb)-Y)^2,x,log.(y)))
        SEYmean=sum(map(Y->(Y-mean(log.(y)))^2,log.(y)))
        R2=1-(SELine/SEYmean)

        pr ? println(" **Exponential Regression \n form y = b*exp(mx) \n m=$m \n b=$b \n r2=$R2") : :none
        return m, b,R2
    elseif TypeReg=="Power"
        m=(mean(log.(x))*mean(log.(y))-mean(log.(x).*log.(y) ) )/(mean(log.(x))^2-mean(log.(x).^2))
        Logb=mean(log.(y))-m*mean(log.(x))
        b=exp(Logb)
        SELine=sum(map((X,Y)->((X*m+Logb)-Y)^2,log.(x),log.(y)))
        SEYmean=sum(map(Y->(Y-mean(log.(y)))^2,log.(y)))
        R2=1-(SELine/SEYmean)

        pr ? println(" **Power Regression \n form y = b x^(m) \n m=$m \n b=$b \n r2=$R2") : :none
        return m, b,R2
    end

end


text:
const TupleN{T,N} = NTuple{N,T}
const AbstractTensor{N,T} = AbstractArray{T,N}

"""
    reshape_maybe(x, shape)

Like `reshape(x, shape)`, except that zero-dimensional outputs are returned as scalars.
"""
reshape_maybe(x::Number, ::Tuple{}) = x
reshape_maybe(x::AbstractArray, ::Tuple{}) = only(x)
reshape_maybe(x::AbstractArray, sz::Tuple{Int,Vararg{Int}}) = reshape(x, sz)
reshape_maybe(x::Number) = x
reshape_maybe(x::AbstractArray, sz::Int...) = reshape_maybe(x, sz)

"""
    samepad(kernel_size, dilation = 1)

Computes padding such that the input and output spatial sizes of the convolution
are the same.

If stride > 1, then `size(output, d) * stride[d] == size(input, d)` for
each spatial dimension `d`.
"""
function samepad end

function samepad(kernel_size::NTuple{K,Int}, dilation::NTuple{K,Int}) where {K}
    sz = (kernel_size .- 1) .* dilation
    return ntuple(Val(2K)) do k
        isodd(k) ? cld(sz[(k + 1) ÷ 2], 2) : fld(sz[k ÷ 2], 2)
    end
end

function samepad(kernel_size::NTuple{K,Int}, dilation::Int = 1) where {K}
    dilation_tuple = ntuple(Returns(dilation), Val(K))
    return samepad(kernel_size, dilation_tuple)
end

function parsepad(kernel_size::NTuple{N,Int}, dilation, pad) where {N}
    if pad === :same
        return samepad(kernel_size, dilation)
    elseif pad isa Int
        return ntuple(Returns(pad), Val(2N))
    elseif pad isa NTuple{N,Int}
        return ntuple(d -> pad[d ÷ 2 + 1], Val(2N))
    else
        return pad::NTuple{2N,Int}
    end
end

"""
    splitpad(pad)

Given `pad = (x_lo, x_hi, y_lo, y_hi, ...)`, returns
`lo = (x_lo, y_lo, ...), hi = (x_hi, y_hi, ...)`.
"""
function splitpad(pad::NTuple{N,Int}) where {N}
    @assert iseven(N)
    pad_lo = ntuple(d -> pad[2d - 1], Val(N ÷ 2))
    pad_hi = ntuple(d -> pad[2d], Val(N ÷ 2))
    return (lo = pad_lo, hi = pad_hi)
end

"""
    expand_tuple(Val(N), t)

Expands `t` into a tuple of length `N`.
"""
expand_tuple(::Val{N}, t::Int) where {N} = ntuple(Returns(t), Val(N))
expand_tuple(::Val{N}, t::NTuple{N,Int}) where {N} = t

"""
    output_size(kernel_size, input_size; stride = 1, pad = 0, dilation = 1)

Output size of the convolution.
"""
output_size(k::NTuple{N,Int}, i::NTuple{N,Int}; stride = 1, pad = 0, dilation = 1) where {N} = output_size(
    k, i,
    expand_tuple(Val(N), stride),
    parsepad(k, dilation, pad),
    expand_tuple(Val(N), dilation)
)

function output_size(
    kernel::NTuple{N,Int}, input::NTuple{N,Int},
    stride::NTuple{N,Int}, pad::NTuple{N2,Int}, dilation::NTuple{N,Int}
) where {N,N2}
    @assert N2 == 2N
    return ntuple(Val(N)) do i
        (input[i] + pad[2i - 1] + pad[2i] - (kernel[i] - 1) * dilation[i] - 1) ÷ stride[i] + 1
    end
end

"""
    out2in(j, k; stride = 1, dilation = 1)

Given kernel index `j` and output index `k`, gets the corresponding input index `i`.
"""
function out2in(j::CartesianIndex{N}, k::CartesianIndex{N}; stride = 1, dilation = 1) where {N}
    i = (Tuple(j) .- 1) .* dilation .+ (Tuple(k) .- 1) .* stride .+ 1
    return CartesianIndex(i)
end

"""
    translate(v, channel_size, input_size, Δ; mode = :pad, fillvalue = 0)

Translates an image along input dimensions by `Δ`.
"""
function translate(
    v::AbstractArray,
    channel_size::NTuple{C,Int},
    input_size::NTuple{N,Int},
    Δ::CartesianIndex{N};
    mode::Symbol = :pad,
    fillvalue = zero(eltype(v))
) where {C,N}
    batch_size = size(v)[(C + N + 1):end]
    @assert size(v) == (channel_size..., input_size..., batch_size...)
    if mode === :pad
        Δv = similar(v)
        for i in CartesianIndices(input_size)
            if i + Δ ∈ CartesianIndices(input_size)
                for n in CartesianIndices(batch_size), c in CartesianIndices(channel_size)
                    Δv[c,i,n] = v[c, i + Δ, n]
                end
            else
                for n in CartesianIndices(batch_size), c in CartesianIndices(channel_size)
                    Δv[c,i,n] = fillvalue
                end
            end
        end
        return Δv
    elseif mode === :periodic
        Δv = similar(v)
        for i in CartesianIndices(input_size)
            p = mod1.(Tuple(i + Δ), input_size)
            for n in CartesianIndices(batch_size), c in CartesianIndices(channel_size)
                Δv[c,i,n] = v[c, p, n]
            end
        end
        return Δv
    end
end

function translate(
    v::AbstractArray,
    channel_size::NTuple{C,Int},
    input_size::NTuple{N,Int},
    Δ::NTuple{N,Int};
    mode::Symbol = :pad,
    fillvalue = zero(eltype(v))
) where {C,N}
    return translate(v, channel_size, input_size, CartesianIndex(Δ); mode, fillvalue)
end


text:
using Distributed

"""
    integrate(ϕ, grid::CubicElementGrid, ψ)
    integrate(grid::CubicElementGrid, ρ)
    integrate(ϕ::QuantumState, ψ::QuantumState)

Low lever integration routines. (users should not use these, as they can change)
"""
function integrate(ϕ, grid::CubicElementGrid, ψ)
    ω = ω_tensor(grid)
    c = ϕ.*ψ
    return  @tensor ω[α,I]*ω[β,J]*ω[γ,K]*c[α,β,γ,I,J,K]
end

function integrate(grid::CubicElementGrid, ρ)
    ω = ω_tensor(grid)
    return  @tensor ω[α,I]*ω[β,J]*ω[γ,K]*ρ[α,β,γ,I,J,K]
end

function integrate(ϕ::QuantumState, ψ::QuantumState)
    @assert size(ϕ) == size(ψ)
    integrate(ϕ.psi, ψ.elementgrid, ψ.psi)
end

function integrate(ϕ::QuantumState{Any, Any, Complex}, ψ::QuantumState)
    @assert ϕ.elementgrid == ψ.elementgrid
    integrate(conj(ϕ).psi, ψ.elementgrid, ψ.psi)
end

"""
    bracket(ϕ::QuantumState, ψ::QuantumState)
    bracket(ϕ::QuantumState, op::AbstractOperator, ψ::QuantumState)

Equivalent to Dirac bracket notation <ϕ|ψ> and <ϕ|op|ψ>.
Returns expectation value.
"""
function bracket(ϕ::QuantumState, ψ::QuantumState)
    @assert size(ϕ) == size(ψ)
    return integrate(ϕ.elementgrid, ϕ ⋆ ψ)*unit(ϕ)*unit(ψ)
end

function bracket(ϕ::QuantumState, op::AbstractOperator{1}, ψ::QuantumState)
    @assert size(ϕ) == size(ψ) == size(op)
    return integrate(ϕ.elementgrid, ϕ ⋆ (op*ψ))*unit(ϕ)*unit(ψ)*unit(op)
end

function bracket(ϕ::QuantumState, op::AbstractOperator, ψ::QuantumState)
    @assert size(ϕ) == size(ψ) == size(op)
    return pmap( O->bracket(ϕ, O, ψ),  op)
end

function bracket(op::AbstractOperator, sd::SlaterDeterminant)
    # One electron operator, so operate on each orbital
    S = map(ψ->bracket(ψ,ψ), sd.orbitals )
    val = @distributed (+) for i in axes(S, 1)
        tmp = 2*bracket(sd.orbitals[i], op, sd.orbitals[i])
        p1 = i > 1 ? prod( S[1:i-1] ) : 1
        p2 = i < length(S) ? prod( S[i+1:end] ) : 1
        tmp .* (p1 * p2)
    end
    return val
end

"""
    magnetic_current(ψ::QuantumState, H::HamiltonOperator)
    magnetic_current(ψ::QuantumState, H::HamiltonOperatorMagneticField)

Return magnetic current for systen with given Hamiltonian.

# Returns
`Vector{AbstractArray{Float,6}}`
"""
function magnetic_current(ψ::QuantumState, H::HamiltonOperator)
    p = momentum_operator(H) * (-1u"e_au"/(2H.T.m))
    return _magnetic_current(p, ψ)
end

function magnetic_current(ψ::QuantumState, H::HamiltonOperatorMagneticField)
    p = momentum_operator(H) * (H.q/(2H.T.m))
    return _magnetic_current(p, ψ)
end

"""
    para_magnetic_current(ψ::QuantumState)

Calculates para magnetic current.

# Returns
`Vector{AbstractArray{Float,6}}`
"""
function para_magnetic_current(ψ::QuantumState)
    p = momentum_operator(ψ)
    return _magnetic_current(p, ψ)
end

function _magnetic_current(p, ψ)
    ϕ = conj(ψ)
    return pmap( x->real.(ψ⋆(x*ψ) .+ ϕ⋆(x*ϕ)), p)
end

## Coulomb integral / Poisson equation

function coulomb_tensor(ρ::AbstractArray, transtensor::AbtractTransformationTensor;
                        tmax=nothing, showprogress=false)
    V = similar(ρ)
    V .= 0
    v = similar(ρ)
    ptime = showprogress ? 1 : Inf
    @showprogress ptime "Calculating v-tensor..." for p in Iterators.reverse(eachindex(transtensor.wt))
        T = transtensor[:,:,:,:,p]
        @tensoropt v[α,β,γ,I,J,K] = T[α,α',I,I']*T[β,β',J,J']*T[γ,γ',K,K']*ρ[α',β',γ',I',J',K']
        V = V .+ transtensor.wt[p].*v
    end
    if tmax != nothing
        return muladd.(2/sqrt(π), V, coulomb_correction(ρ, tmax))
    end
    return 2/sqrt(π).*V
end

function coulomb_tensor(ρ::AbstractArray, transtensor::AbstractArray, wt::AbstractVector;
                        tmax=nothing, showprogress=false)
    @assert size(transtensor)[end] == length(wt)
    V = similar(ρ)
    V .= 0
    v = similar(ρ)
    ptime = showprogress ? 1 : Inf
    @showprogress ptime "Calculating v-tensor..." for p in Iterators.reverse(eachindex(wt))
        T = @view transtensor[:,:,:,:,p]
        @tensoropt v[α,β,γ,I,J,K] = T[α,α',I,I']*T[β,β',J,J']*T[γ,γ',K,K']*ρ[α',β',γ',I',J',K']
        V = V .+ wt[p].*v
    end
    if tmax != nothing
        return muladd.(2/sqrt(π), V, coulomb_correction(ρ, tmax))
    end
    return  2/sqrt(π).*V
end

coulomb_correction(ρ, tmax) = 2/sqrt(π)*(π/tmax^2).*ρ

function poisson_equation(ρ::AbstractArray,
                          transtensor::AbtractTransformationTensor,
                          I, J, K;
                          tmax=nothing)
    V = similar(ρ, size(ρ)[1:3]...)
    V .= 0
    v = similar(V)
    @showprogress for p in Iterators.reverse(eachindex(transtensor.wt))
        Tx = transtensor[:,:,:,I,p]
        Ty = transtensor[:,:,:,J,p]
        Tz = transtensor[:,:,:,K,p]
        @tensoropt v[α,β,γ] = Tx[α,α',I'] * Ty[β,β',J'] * Tz[γ,γ',K'] * ρ[α',β',γ',I',J',K']
        #@tensoropt v[α',β',γ',I',J',K'] = Tx[α,α',I'] * Ty[β,β',J'] * Tz[γ,γ',K'] * vρ[α,β,γ]
        V = V .+ transtensor.wt[p].*v
    end
    if tmax != nothing
        return muladd.(2/sqrt(π), V, coulomb_correction(ρ[:,:,:,I,J,K], tmax))
    end
    return  2/sqrt(π).*V
end

function poisson_equation(ρ::AbstractArray{<:Any,6},
                          transtensor::AbtractTransformationTensor,
                          t)
    V = similar(ρ)
    T = similar(ρ, size(transtensor)[1:end-1])
    T .= transtensor[:,:,:,:,t]
    @tensoropt V[α,β,γ,I,J,K] = T[α,α',I,I'] * T[β,β',J,J'] * T[γ,γ',K,K'] * ρ[α',β',γ',I',J',K']
    return  (2/sqrt(π)*transtensor.wt[t]).*V
end

function poisson_equation!(V::AbstractArray{<:Any,6},
                          ρ::AbstractArray{<:Any,6},
                          transtensor::AbtractTransformationTensor,
                          t)
    T = transtensor[:,:,:,:,t]
    @tensoropt V[α,β,γ,I,J,K] = T[α,α',I,I'] * T[β,β',J,J'] * T[γ,γ',K,K'] * ρ[α',β',γ',I',J',K']
    return  (2/sqrt(π)*transtensor.wt[t]).*V
end

function poisson_equation!(V::AbstractArray{<:Any,6},
                          ρ::AbstractArray{<:Any,6},
                          T::AbstractArray{<:Any,4},
                          wt)
    @assert size(V) == size(ρ)
    @tensoropt V[α,β,γ,I,J,K] = T[α,α',I,I'] * T[β,β',J,J'] * T[γ,γ',K,K'] * ρ[α',β',γ',I',J',K']
    V .*= (2/sqrt(π)*wt)
    return V
end

function poisson_equation(ρ::AbstractArray, transtensor::AbtractTransformationTensor;
                          tmax=nothing, showprogress=false)

    V = ρ.*transtensor.wt[1]
    nt = size(transtensor)[end]
    @debug "nt=$nt"
    @debug "V type = $(typeof(V))"
    ptime = showprogress ? 0.3 : Inf
    if nworkers() > 1
        tmp = similar(V) # Yes it is a hack
        V = @showprogress ptime "Poisson equation... " @distributed (+) for t in 1:nt
            poisson_equation!(tmp, ρ, transtensor, t)
        end
    else
        tmp = similar(V)
        @showprogress ptime "Poisson equation... " for t in 1:nt
            V .+= poisson_equation!(tmp, ρ, transtensor[:,:,:,:,t], transtensor.wt[t])
        end
    end
    if tmax != nothing
        return V .+ coulomb_correction(ρ, tmax)
    end
    return V
end

function poisson_equation(ψ::QuantumState, transtensor::AbtractTransformationTensor;
                          tmax=nothing, showprogress=false)
    @assert dimension(ψ) == dimension(u"bohr^-2")
    ψ = uconvert(u"bohr^-2", ψ)
    V = poisson_equation(ψ.psi, transtensor, tmax=tmax, showprogress=showprogress)
    return QuantumState(ψ.elementgrid, V)
end


text:
cd("/home/albyback/git/dotlabelproject/OuterProject/SparseInverseProblems.jl/examples/simple_example/")
include("simple_example.jl")
using Compose
evaluation_points = -5.0:0.25:5.0
model = SimpleExample(evaluation_points, -10.0:0.5:10.0)
means = [-2.0,1.0,3.5]
k = length(means)
weights = ones(k)+randn(k)*0.3
weights[1] = 0.2
target = max.(phi(model,Matrix(means'),weights) + randn(length(evaluation_points))*0.1,0.0);

# using Plots
# using PlutoUI
# using Images, FileIO
# heatmap(model.grid, aspect_ratio=1, c=:grays)
# plot(phi(model,Matrix(means'),weights) + randn(length(evaluation_points))*0.1)
# plot(target)

function callback(old_thetas,thetas, weights,output,old_obj_val)
  #evalute current OV
  new_obj_val,t = loss(LSLoss(), output - target)
  println("gap = $(old_obj_val - new_obj_val)")
  if old_obj_val - new_obj_val < 1E-1
    return true
  end
  return false
end

# ADCG(model, LSLoss(), target,10000.0;  callback=callback)

(means_est,weights_est) = ADCG(model, LSLoss(), target,100.0;  callback=callback)
print(means_est)
print(means)

using Plots
using PlutoUI
using Images, FileIO

# using PyPlot
# pygui(true)
# #draw the observations alone
plot(evaluation_points, target)

#draw the observations with the locations and weights of the true blurs
# anno = Guide.annotation(
#        compose(context(), circle([means,means], [zeros(k),weights], [2mm]), fill(nothing),
#        stroke("blue")))
# draw(SVG("truth.svg", 5inch, 2inch), plot(x=evaluation_points,y=target,anno))

#draw the estimation along with the predicted means/observations
anno = Guide.annotation(
       compose(context(), circle([means_est',means_est'], [zeros(k),weights_est], [2mm]), fill(nothing),
       stroke("red")))
output_est = phi(model, means_est, weights_est)
draw(SVG("est.svg", 5inch, 2inch), plot(x=evaluation_points, y=output_est,anno))


text:
q = Quaternion(1, 0, 0, 0)
r = λmap(q)
p = compressedλmap(q)

@test isapprox(r, ℝ³(1, 0, 0))
@test isapprox(p, ℝ³(0.7615941559557649, 0, 0))

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using EstimHill
using Documenter

makedocs(;
    modules=[EstimHill],
    authors="Shao-Ting Steven Chiu <r07945001@ntu.edu.tw>",
    repo="https://github.com/stevengogogo/EstimHill.jl/blob/{commit}{path}#L{line}",
    sitename="EstimHill.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://stevengogogo.github.io/EstimHill.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/stevengogogo/EstimHill.jl.git",
)


text:
# Autogenerated wrapper script for Tcl_jll for x86_64-apple-darwin
export libtcl

using Zlib_jll
JLLWrappers.@generate_wrapper_header("Tcl")
JLLWrappers.@declare_library_product(libtcl, "@rpath/libtcl8.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll)
    JLLWrappers.@init_library_product(
        libtcl,
        "lib/libtcl8.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export crossover

function onepoint_crossover(parent1::GRN, parent2::GRN)
    point = rand(1:max(length(parent1), length(parent2)))
    order = rand()
    if rand(Bool)
        p1 = parent1
        p2 = parent2
    else
        p1 = parent2
        p2 = parent1
    end
    point = min(min(point, length(p1)), length(p2)-1)
    ids = [p1.ids[1:point]; p2.ids[point+1:end]]
    enh = [p1.enh[1:point]; p2.enh[point+1:end]]
    inh = [p1.inh[1:point]; p2.inh[point+1:end]]
    nreg = length(ids) - p1.nin - p1.nout
    beta = parent1.beta
    if rand(Bool)
        beta = parent2.beta
    end
    delta = parent1.delta
    if rand(Bool)
        delta = parent2.delta
    end
    GRN(p1.nin, p1.nout, nreg, ids, enh, inh, beta, delta)
end

function crossover(parent1::GRN, parent2::GRN, config::Config)
    eval(config.crossover)(parent1, parent2)
end


text:
try
    using StopWatchesTests
    true
catch
    false
end || begin
    push!(LOAD_PATH, joinpath(@__DIR__, "StopWatchesTests"))
    using StopWatchesTests
end


text:
include("CommonToAll.jl")
using .CommonToAll
include("FooPhysics.jl") # example for API
using. FooPhysics # example for API
include("FooPhysicsInversion.jl") # example for API
using .FooPhysicsInversion # example for API
include("LineRegression.jl")
using .LineRegression
include("SurfaceRegression.jl")
using .SurfaceRegression
include("CSEM1DEr.jl") # add this before using HED for CSEM
using .CSEM1DEr
include("CSEM1DInversion.jl")
using .CSEM1DInversion
include("ImageRegression.jl")
using .ImageRegression
include("AEM_VMD_HMD.jl") # add this before using VMD or HMD for land/AEM geophysics
using .AEM_VMD_HMD
include("SkyTEM1DInversion.jl")
using .SkyTEM1DInversion
include("TEMPEST1DInversion.jl")
using .TEMPEST1DInversion
include("MT1D.jl")
using .MT1D
include("MT1DInversion.jl")
export Operator, get_misfit, Sounding, makeoperator, make_tdgp_opt, getresidual

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
using LibCURL

function callSOAP(url, infilename, outfilename)
    rfp = open(infilename, "r")
    wfp = open(outfilename, "w+")

    header = curl_slist_append(header, "Content-Type:text/xml")
    header = curl_slist_append(header, "SOAPAction: rsc");
    header = curl_slist_append(header, "Transfer-Encoding: chunked")
    header = curl_slist_append(header, "Expect:")

    curl = curl_easy_init();
    curl_easy_setopt(curl, CURLOPT_URL, URL)
    curl_easy_setopt(curl, CURLOPT_POST, 1L)
    curl_easy_setopt(curl, CURLOPT_READFUNCTION, read_data)
    curl_easy_setopt(curl, CURLOPT_READDATA, rfp)
    curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, write_data)
    curl_easy_setopt(curl, CURLOPT_WRITEDATA, wfp)
    curl_easy_setopt(curl, CURLOPT_HTTPHEADER, header)
    curl_easy_setopt(curl, CURLOPT_POSTFIELDSIZE_LARGE, (curl_off_t)-1)
    curl_easy_setopt(curl, CURLOPT_VERBOSE,1L)
    curl_easy_perform(curl)

    curl_easy_cleanup(curl)
end

try
    callSOAP(ARGS[1], ARGS[2], ARGS[3])
catch y
    println("Usage : $(@__FILE__) <URL of WSDL> <Input file path> <Output File Path>")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export distance

"""
    distance(H1::AbstractHyperrectangle{N}, H2::AbstractHyperrectangle{N};
             p::Real=2) where {N<:Real}

Compute the standard distance between two hyperrectangular sets, defined as

```math
    \\inf_{x \\in H_1, y \\in H_2} \\{ d(x, y) \\}.
```

### Input

- `H1` -- hyperrectangular set
- `H2` -- hyperrectangular set
- `p`  -- (optional; default: `2`) value of the ``p``-norm

### Output

The distance, which is zero if the sets intersect and otherwise the ``p``-norm
of the shortest line segment between any pair of points.

### Notes

See also [`hausdorff_distance`](@ref) for an alternative distance notion.
"""
function distance(H1::AbstractHyperrectangle{N},
                  H2::AbstractHyperrectangle{N};
                  p::Real=2) where {N<:Real}
    n = dim(H1)
    @assert n == dim(H2) "incompatible set dimensions $n and $(dim(H2))"

    d = Vector{N}(undef, n)
    @inbounds for i in 1:n
        # find relative position in dimension i
        # (if c1 == c2, the results are equivalent independent of the branch)
        if center(H1, i) >= center(H2, i)
            lhs = low(H1, i)
            rhs = high(H2, i)
        else
            lhs = low(H2, i)
            rhs = high(H1, i)
        end
        if _leq(lhs, rhs)
            d[i] = zero(N)
        else
            d[i] = rhs - lhs
        end
    end
    return norm(d, p)
end


text:
module ApproxMinimumDegree
const DEBUG = false

include("NodeWeightPQs/NWLLPriorityQueues.jl")
import .NWLLPriorityQueues
include("QuotientGraphs/QuotientGraphs.jl")
import .QuotientGraphs

function amd{Ti,Tv}(C::SparseMatrixCSC{Tv,Ti})
    matorder = C.n

    ### DATA STRUCTURE SETUP ###

    # Quotient Graph (elimination and supervariable quotient)
    #
    # As with most local preordering codes, this code maintains a quotient graph model
    # of the elimination graph rather than an explicit representation of the elimination
    # graph. Here by quotient graph we mean quotient in both elimination and
    # supervariable senses.
    #
    # The quotient graph data structure we've abstracted into the module QuotientGraphs.
    # In this method we primarily manipulate it via a relatively opaque, high-level
    # interface. For more on the QuotientGraph data structure, see Quotientgraphs.
    #
    # Build the initial quotient graph from input matrix C. The method below currently
    # assumes that C is symmetric and conatins all entries, not merely an upper or lower
    # triangle.
    qg = QuotientGraphs.QuotientGraph(C)
    # FUTURE Allow input matrix C to be asymmetric or contain only the upper or lower
    # triangle of a symmetric matrix, and build the initial quotient graph for (C + C^T)

    # Node weight tracking data structure
    #
    # We track the weight (minimum degree, minimum fill, approximate minimum degree, or
    # what have you) of each live supervariable (quotient graph node) in a specialized
    # priority-queue-like data structure. Here we primarily manipulate this data structure
    # via a relatively opaque, high-level interface. See the NWLLPriorityQueues module
    # for more information on the data structure.
    #
    # Initialize the node weight priority queue and then populate it
    nodeweightpq = NWLLPriorityQueues.NWLLPriorityQueue{Ti}(matorder)
    for k in 1:matorder
        NWLLPriorityQueues.enqueue!(nodeweightpq, k, qg.setstorage[qg.setpointers[k]])
    end
    # FUTURE Initializing the priority queue in one shot from the quotient graph may be
    # more efficient than building it incrementally as we do just above.

    # Allocate various workspaces
    #
    # Generic workspaces for various operations (for example by the QuotientGraph)
    workpoolA = Vector{Ti}(matorder)
    workpoolB = Vector{Ti}(matorder)
    #
    # Workspaces used for specific steps in approximate external degree calculation
    # W = Vector{Ti}(n)
    # We now attach W to workpoolA as operations on either are temporally segregated.
    W = workpoolA
    # lL = Vector{Ti}(n)
    # We now attached lL to nodeweightpq.weights as entries go into lL only after the
    # corresponding weight no longer matters.
    lL = nodeweightpq.weights

    ### MAIN ELIMINATION LOOP ###
    elimstepn = 1
    while elimstepn <= matorder
        # Identify a minimum-weight supervariable / quotient graph node by querying the
        # the node weight priority queue. This minimum-weight supervariable / node,
        # commonly the "pivot" p in the literature, we eliminate from the quotient
        # graph in this iteration of the main loop and then add to the ordering.
        p, pw = NWLLPriorityQueues.dequeue!(nodeweightpq)
        DEBUG && println("Beginning elimination step $elimstepn. Pivot $p of weight $pw.")

        # Transform the quotient graph, eliminating supervariable p. Specifically, we
        # transform supervariable / quotient graph node p into quotient graph element p.
        QuotientGraphs.elementize!(qg, nodeweightpq, p, workpoolA, workpoolB)

        # Update the remaining supervariables' / quotient graph nodes' weights.
        # Specifically, we need udpate the approximate external degree of each supervariable
        # i in [the set of supervariables to which now-element p links, L_p].
        #
        # First we compute the number of simple variables that now-element p links to, i.e.
        # the simple-variable count of set L_p. We calculate this value once and cache it
        # as we potentially use it many times and it does not change during its lifetime
        # (as L_p does not change the set of _simple_ variables it contains, though the
        # set of _super_variables it contains may change, till another element absorbs
        # element p and we no longer need this information). For efficiency we've fused
        # the loop over L_p calculating this quantity into similar loops over L_p for
        # Alg. 2 described below. The old, unfused code read:
        # lLp = 0
        # for kLp in QuotientGraphs.rangeL(qg, p)
        #     iLp = qg.setstorage[kLp]
        #     lLp += QuotientGraphs.numsimplevars(qg, iLp)
        # end
        # lL[p] = lLp
        #
        # Second, for each supervariable i in L_p, and for each element e in [the set of
        # elements to which supervariable i links, E_i], we must calculate
        # length(setdiff(L_e, L_p)). We accomplish this via Algorithm 2 from the seminal
        # AMD paper by Amestoy, Davis, and Duff. Notes:
        #
        # Algorithm 2 calculates the relevant quantity for all live elements,
        # but it does so implicitly. That is, it explicitly calculate the relevant quantity
        # for all (e in E[i] for i in L[p]). Calculating the relevant quantity
        # for live elements outside that set is trivial, as for element e outside that set
        # L[e] and L[p] are necessarily disjoint, not being connected to any of the same
        # live supervariables. Which...
        #
        # .. is fine, as updating the approximate external degree of supervariables i
        # in L[p] only requires the relevant quantity for (e in E[i] for i in
        # L[p]), i.e. only where the calculation is nontrivial. We restrict our calculation
        # to only the required subset.
        #
        # We begin Algorithm 2 with the implicit marking step, setting each W[e] = -1
        for kLp in QuotientGraphs.rangeL(qg, p)
            iLp = qg.setstorage[kLp]
            for kEi in QuotientGraphs.rangeE(qg, iLp)
                eEi = qg.setstorage[kEi]
                W[eEi] = -1
            end
        end
        # FUTURE We should be able to eliminate this marking step by introducing slight
        # additional complexity into the latter part of Alg. 2 and using some assumptions
        # on the contents of the workspace in which we place W.
        #
        # Then we perform Algorithm 2 proper, calculating the set of values simultaneously.
        #
        # NOTE We only access lL[e = p] in this loop to set W[e = p], which hypothetically
        # becomes simplevarcount(setdiff(L_{e = p}, L_p)), but of course we know that value
        # trivially (zero) and, though we might use it below, we don't actually need it.
        # This allowed us to fuse the lL[p] calculation discussed above with the following
        # loop. This way we avoid evaluating QuotientGraphs.numsimplevars(qg, i) twice
        # for each i in L_p, we may geet better cache behavior once we fuse lL with
        # qg.nodestatus instead of nodeweightpq.weights, and we are free to get rid of
        # the above marking loop if we find a nice way to do so.
        lLp = zero(Ti)
        for kLp in QuotientGraphs.rangeL(qg, p)
            i = qg.setstorage[kLp]
            lenIi = QuotientGraphs.numsimplevars(qg, i)
            lLp += lenIi
            for kEi in QuotientGraphs.rangeE(qg, i)
                e = qg.setstorage[kEi]
                We = W[e]
                if We < 0
                    W[e] = lL[e] - lenIi
                else
                    W[e] = We - lenIi
                end
            end
        end
        W[p] = 0
        lL[p] = lLp
        #
        # Third, with the quantities calculated above on hand, we calculate and update
        # the approximate external degrees of the supervariables i in L_p.
        for kLp in QuotientGraphs.rangeL(qg, p)
            iLp = qg.setstorage[kLp]
            # Three terms contribute to the approximate external degree for supervariable
            # i. We calculate each in turn, then the approximate external degree proper.
            #
            # Term 1, the simple-variable count of setdiff(L_p, I_i). We know that
            # supervariable i is an entry of set L_p, and hence in terms of simple
            # variables all entries in I_i are entries of L_p. Hence, having calculated
            # the simple-variable count of L_p above, we can calculate the simple-var count
            # of setdiff(L_p, I_i) trivially as simplevarcount(L_p) - simplevarcount(I_i):
            lengthLpLessIi = lLp - QuotientGraphs.numsimplevars(qg, iLp)
            #
            # Term 2, the simple-variable count of setdiff(A_i, I_i). If we are
            # maintaining A_i properly, we know that supervariable i should not be an
            # entry of set A_i; that is, when we build our QuotientGraph, we exclude
            # i from A_i, and if we maintain A_i properly, we should not reintroduce i
            # into A_i, and when merging supervariables etcetera we carefully strip
            # entries in I_i from A_i. So setdiff(A_i, I_i) should just be A_i in our
            # case. Hence the simple-variable count of that setdiff is just the
            # simple-variable count of A_i, which we calculate readily.
            #
            # Term 3, sum_{e in E_i} simplevarcount(setdiff(L_e, L_p)). We precalculated
            # each term in the sum above via Algorithm 2. Hence here we need simply sum
            # the appropriate precalculated terms.
            #
            # Given that A_i and E_i are stored contiguously in qg.setstorage and in that
            # order, we calculate terms 2 and 3 fused as follows for efficiency.
            #
            rangeAi, rangeEi = QuotientGraphs.rangesAandE(qg, iLp)
            # FUTURE This range-pair-retrieval operation is slow. Fix.
            lengthAiLessIi = zero(Ti)
            for kAi in rangeAi
                jAi = qg.setstorage[kAi]
                lengthAiLessIi += QuotientGraphs.numsimplevars(qg, jAi)
            end
            sumEi_lengthLeLessLp = zero(Ti)
            for kEi in rangeEi
                eEi = qg.setstorage[kEi]
                # eEi == p && continue # W[p] == 0 from above
                sumEi_lengthLeLessLp += W[eEi]
            end
            # We combine the three terms calculated above into the approximate
            # external degree proper.
            oldapproxextdeg = NWLLPriorityQueues.getnodeweight(nodeweightpq, iLp)
            bound1_maxextdegree = matorder - elimstepn
            bound2_crudemaxfill = oldapproxextdeg + lengthLpLessIi
            bound3_tightextdegest = lengthLpLessIi + lengthAiLessIi + sumEi_lengthLeLessLp
            newapproxextdeg = min(bound1_maxextdegree, bound2_crudemaxfill, bound3_tightextdegest)
            #
            # Finally, we update this supervariable's / quotien graph node's
            # approximate external degree / weight in the node weight priority queue.
            NWLLPriorityQueues.update!(nodeweightpq, iLp, newapproxextdeg)
        end
        elimstepn += QuotientGraphs.numsimplevars(qg, p)
    end

    postorder!(qg, workpoolA, workpoolB)
    return workpoolB
end

"""
    postorder!{T<:Signed}(qg::QuotientGraph{T}, perm::Vector{T}, iperm::Vector{T})

Given a `QuotientGraph` on which we completed elimination/ordering, returns an assembly-tree
postordering as a permutation vector in `perm` and as an inverse permutation vector in
`iperm`.
"""
function postorder!{T<:Signed}(
        qg::QuotientGraphs.QuotientGraph,
        perm::Vector{T},
        iperm::Vector{T})
    supervarcards = qg.nodestatus
    rootwardlinks = qg.setpointers
    leafwardlinks = perm
    # On entry, rootwardlinks and supervarcards together define the assembly tree generated
    # during the ordering process. Specifically, for index i into supervarcards,
    # superacards[i] indicates one of two possible things:
    #
    # (1) supervarcards[i] == 0 indicates that supervariable i was amalgmated into another
    # supervariable prior to (mass) elimination during the ordering process.
    #
    # (2) supervarcards[i] (== k) > 0 indicates that supervariable i was representative
    # when it was eliminated, and hence became an element. k is the number of simple
    # variables that supervariable i represented on (mass) elimination.
    #
    # For index i into rootwardlinks, rootwardlinks[i] indicates one of three things:
    #
    # (1) rootwardlinks[i] >= 0 indicates that as an element i was never absorbed into
    # another element, and hence constitutes a root of the assembly tree.
    #
    # (2) If supervarcards[i] == 0 and hence as above i was a supervariable amalgamated
    # into another supervariable prior to (mass) elimination during the ordering process,
    # then rootwardslinks[i] (== -j) < 0 indicates that supervariable i was amalgamated
    # into supervariable j.
    #
    # (3) If supervarcards[i] > 0 and hence as above supervariable i was representative
    # when it was eliminated and so became an element, then rootwardlinks[i] (== -j) < 0
    # indicates that element i was absorbed into elemnt j during the ordering process.
    #
    # In this fashion rootwardlinks and supervarcards together define the assembly tree
    # generated during the ordering process by way of child-to-parent links.
    #
    # Our objective is to generate a postordering of this assembly tree as a permutation
    # vector `perm` and an inverse permutation vector `iperm`. We do this via a
    # depth-first-search of the assembly tree (which can be a forest, and we handle that
    # case as well). To perform an efficient DFS, we need parent-to-child links in addition
    # to the child-to-parent links we have. Hence this method involves two stages: In
    # the first stage, we build such parent-to-child links in leafwardlinks. In the second
    # stage, we perform the depth-first-search of the assembly tree.
    #
    # Moreover, we need not only parent-to-child links, but we need an ordering among the
    # children of a given parent. Specifically, we need that those children corresponding to
    # assembly operations be 'elder' to those corresponding to supervariable amalgamations.
    # Hence, for nodes with younger siblings, we replace the link to the node's parent in
    # rootwardlinks to a link to the node's next-younger sibling. Having this information
    # allows us to efficiently walk the tree ordering elder siblings (assembly operations)
    # before younger siblings (supervariable amalgamations). To effect this, we build
    # leafwardlinks in two stages, first sweeping over amalgamated supervariables, and then
    # over elements (assembly operations), ensuring the ordering we need.
    #
    # For more information on this process, see the paper and technical report on MA27
    # from Duff and Reid.
    #
    # Additionally, as a side effect of stage one we build a linked list of assembly tree
    # roots. nextroot points to the first root. rootwardlinks[nextroot] then points to the
    # second root if it exists, rootwardslinks[rootwardlinks[nextroot]] then points
    # to the third root if it exists, and so on. A zero in the chain indicates its end.

    # Stage one / initialization.
    nextroot = 0
    fill!(leafwardlinks, 0)
    # Stage one / substage one, building parent-to-child and child-to-ysibling links
    # for amalgamated supervariables.
    for child in 1:length(supervarcards)
        supervarcards[child] == 0 || continue # Skip all but amalgamated supervariables
        nextroot = parentchildlink!(child, rootwardlinks, leafwardlinks, nextroot)
    end
    # Stage one / substage two, building parent-to-child and child-to-ysibling links
    # for absorbed elements.
    for child in 1:length(supervarcards)
        supervarcards[child] != 0 || continue # Skip all but elements
        nextroot = parentchildlink!(child, rootwardlinks, leafwardlinks, nextroot)
    end

    # Stage two, depth-first-search through the assembly tree to generate the postorder.
    korder = 1
    while nextroot != 0
        # We begin the DFS from the assembly-tree root identified in nextroot. Presently
        # rootwardlinks[nextroot] either contains 0 indicating nextroot is the last root
        # in the rootchain, or contains a positive integer identifying the next root in
        # the root chain. During the DFS, we need rootwardlinks[root] = 0 such that we
        # can terminate DFS once we've returned to the root. Hence we need set
        # rootwardlinks[nextroot] = 0, but we need the present value after we've completed
        # this DFS. Hence, retrieve rootwardlinks[nextroot], store it in nextroot after
        # seeding the present DFS, and then set rootwardlinks[nextroot] = 0.
        presnode = nextroot
        nextroot = rootwardlinks[presnode]
        rootwardlinks[presnode] = 0
        while true
            child = leafwardlinks[presnode]
            if child != 0
                # If the present node has a child, step to the child and then burn the
                # bridge so that we cannot walk down this path again.
                leafwardlinks[presnode] = 0
                presnode = child
                continue
            else
                # If the present node has no child, add this node to both representations
                # of the postordering.
                leafwardlinks[presnode] = korder # permutation representation
                iperm[korder] = presnode # inverse permutation represeentation
                korder += 1
                # Retrieve the present node's parent or younger sibling.
                ysiborpar = rootwardlinks[presnode]
                # If this node has no parent or younger sibling, we've returned to the root,
                # in which case terminate the search.
                ysiborpar >= 0 && break
                # Otherwise, step to this node's parent or sibling and continue the search.
                presnode = -ysiborpar
            end
        end
    end
end

@inline function parentchildlink!(child, rootwardlinks, leafwardlinks, nextroot)
    parent = rootwardlinks[child]
    if parent >= 0
        # rootwardlinks[child] >= 0 indicates child was not an amalgamated supervariable or
        # absorbed element, but rather an ultimately unabsorbed element and hence assembly
        # tree root node. Add this node to the chain of roots.
        rootwardlinks[child] = nextroot
        return child # nextroot = child
    else # parent < 0
        # If leafwardlinks[parent] already contains a link (i.e. != 0 and so parent must
        # already have at least one child attached), then child must have a younger sibling.
        # In this case, we replace child's child-to-parent link in rootwardlinks with a
        # child-to-ysibling link, and then place a parent-to-child link to this child in
        # leafwardlinks over the younger sibling's parent-to-child link.
        ysibling = leafwardlinks[-parent]
        leafwardlinks[-parent] = child
        ysibling != 0 && (rootwardlinks[child] = -ysibling)
    end
    return nextroot
end

end


text:
"""
    calculate_electric_potential(det::SolidStateDetector{T}; <keyword arguments>) where {T}

Compute the electric potential for the given Detector `det` on an adaptive grid
through successive over relaxation. It returns a collection struct `PotentialSimulationSetup{T}` which stores
the potential, the charge density, the dielectric distribution, pointtypes and the final grid.

There are serveral `<keyword arguments>` which can be used to tune the computation:

# Keywords
- `convergence_limit::Real`: `convergence_limit` times the bias voltage sets the convergence limit of the relaxation. The convergence value is the absolute maximum difference of the potential between two iterations of all grid points. Default of `convergence_limit` is `5e-6` (times bias voltage).
- `max_refinements::Int`: number of maximum refinements. Default is `2`. Set it to `0` to switch off refinement.
- `refinement_limits::Vector{Real}`: vector of refinement limits for each dimension (in case of cylindrical coordinates the order is `r`, `φ`, `z`). A refinement limit (e.g. `refinement_limits[1]`) times the bias voltage of the detector `det` is the maximum allowed voltage difference between two neighbouring grid points in the respective dimension. When the difference is larger, new points are created inbetween. Default is `[1e-4, 1e-4, 1e-4]`.
- `init_grid_spacing::Vector{Real}`: vector of the initial distances between two grid points for each dimension. For normal coordinates the unit is meter. For angular coordinates, the unit is radiance. It prevents the refinement to make the grid to fine. Default is `[0.005, 10.0, 0.005]``.
- `min_grid_spacing::Vector{Real}`: vector of the mimimum allowed distance between two grid points for each dimension. For normal coordinates the unit is meter. For angular coordinates, the unit is radiance. It prevents the refinement to make the grid to fine. Default is [`1e-4`, `1e-2`, `1e-4`].
- `grid::Grid{T, N, S}`: Initial grid used to start the simulation. Default is `Grid(detector, init_grid_spacing=init_grid_spacing)`.
- `depletion_handling::Bool`: enables the handling of undepleted regions. Default is false.
- `use_nthreads::Int`: Number of threads to use in the computation. Default is `Base.Threads.nthreads()`. The environment variable `JULIA_NUM_THREADS` must be set appropriately before the Julia session was started (e.g. `export JULIA_NUM_THREADS=8` in case of bash).
- `sor_consts::Vector{<:Real}`: Two element array. First element contains the SOR constant for `r` = 0. Second contains the constant at the outer most grid point in `r`. A linear scaling is applied in between. First element should be smaller than the second one and both should be ∈ [1.0, 2.0]. Default is [1.4, 1.85].
- `max_n_iterations::Int`: Set the maximum number of iterations which are performed after each grid refinement. Default is `10000`. If set to `-1` there will be no limit.
- `verbose::Bool=true`: Boolean whether info output is produced or not.
- `init_grid_spacing::Vector{<:Real}`: Initial spacing of the grid. Default is [2e-3, 5, 2e-3] <=> [2mm, 5 degree, 2mm ]
"""
function calculate_electric_potential(  detector::SolidStateDetector{T, :cylindrical};
                                        init_grid_spacing::Vector{<:Real} = [0.005, deg2rad(5.0), 0.005], # 2mm, 10 degree, 2 mm
                                        grid::Grid{T, N, S} = Grid(detector, init_grid_spacing=init_grid_spacing),
                                        convergence_limit::Real = 5e-6,
                                        max_refinements::Int = 3,
                                        refinement_limits::Vector{<:Real} = [1e-4, 1e-4, 1e-4],
                                        min_grid_spacing::Vector{<:Real} = [1e-4, 1e-2, 1e-4],  # mm, degree, mm
                                        depletion_handling::Bool = false,
                                        use_nthreads::Int = Base.Threads.nthreads(),
                                        sor_consts::Vector{<:Real}=[1.4, 1.85],
                                        max_n_iterations::Int=10000,
                                        verbose::Bool=true,
                                        ) where {T, N, S}
    refinement_limits::Vector{T} = T.(refinement_limits)
    min_grid_spacing::Vector{T} = T.(min_grid_spacing)
    sor_consts::Vector{T} = T.(sor_consts)
    refine::Bool = max_refinements > 0 ? true : false
    only_2d::Bool = length(grid.axes[2]) == 1 ? true : false
    check_grid(grid)
    cyclic::T = grid.axes[2].interval.right - grid.axes[2].interval.left
    n_φ_sym::Int = only_2d ? 1 : Int(round(T(2π) / cyclic, digits = 3))
    if use_nthreads > Base.Threads.nthreads()
        use_nthreads = Base.Threads.nthreads();
        @warn "`use_nthreads` was set to `1`. The environment variable `JULIA_NUM_THREADS` must be set appropriately before the julia session is started."
    end

    n_φ_sym_info_txt = if only_2d
        "φ symmetry: Detector is φ-symmetric -> 2D computation."
    else
        "φ symmetry: calculating just 1/$(n_φ_sym) in φ of the detector."
    end

    fssrb::PotentialSimulationSetupRB{T, 3, 4, :cylindrical} = PotentialSimulationSetupRB(detector, grid);

    if verbose
        println("Electric Potential Calculation\n",
                #"Bulk type: $(detector.bulk_type)\n",
                "Bias voltage: $(fssrb.bias_voltage) V\n",
                "$n_φ_sym_info_txt\n",
                "Precision: $T\n",
                "Convergence limit: $convergence_limit => $(round(abs(fssrb.bias_voltage * convergence_limit), sigdigits=2)) V\n",
                "Threads: $use_nthreads\n",
                "Coordinate system: Cylindrical\n",
                "Initial grid dimension: $(size(grid))\n",
                "Refine? -> $refine\n",
                "Refinement parameters:\n",
                "\tmaximum number of refinements: $(max_refinements)\n",
                "\tminimum grid spacing:\n",
                "\t\tr: $(min_grid_spacing[1]) m\n",
                "\t\tφ: $(min_grid_spacing[2]) rad\n",
                "\t\tz: $(min_grid_spacing[3]) m\n",
                "\tRefinement limits:\n",
                "\t\tr: $(refinement_limits[1]) -> $(round(abs(fssrb.bias_voltage * refinement_limits[1]), sigdigits=2)) V\n",
                "\t\tφ: $(refinement_limits[2]) -> $(round(abs(fssrb.bias_voltage * refinement_limits[2]), sigdigits=2)) V\n",
                "\t\tz: $(refinement_limits[3]) -> $(round(abs(fssrb.bias_voltage * refinement_limits[3]), sigdigits=2)) V\n",
                ""
        )
    end

    n_iterations_between_checks::Int = 500
    update_till_convergence!(   fssrb, convergence_limit, fssrb.bias_voltage,
                                depletion_handling = Val{depletion_handling}(), only2d = Val{only_2d}(),
                                use_nthreads=use_nthreads, max_n_iterations=max_n_iterations, n_iterations_between_checks = n_iterations_between_checks )

    potential::Array{T, 3} = ElectricPotentialArray(fssrb)

    refinement_counter::Int = 0
    if refine
        new_inds = check_for_refinement(potential, grid, abs.(refinement_limits .* fssrb.bias_voltage), min_grid_spacing)
        while !isempty(new_inds[1]) || !isempty(new_inds[2]) || !isempty(new_inds[3])
            refinement_counter += 1
            if refinement_counter <= max_refinements
                potential, grid = add_points_and_interpolate(potential, grid, new_inds...)
                if verbose println("Refinement $refinement_counter:\tNew grid size: $(size(potential))") end
                fssrb = PotentialSimulationSetupRB(detector, grid, potential, sor_consts = refinement_counter == max_refinements ? T[1, 1] : sor_consts)
                n_iterations_between_checks = div(n_iterations_between_checks, 2)
                if n_iterations_between_checks < 50 n_iterations_between_checks = 50 end
                update_till_convergence!(fssrb, convergence_limit, fssrb.bias_voltage, depletion_handling = Val{depletion_handling}(), only2d = Val{only_2d}(), use_nthreads=use_nthreads, max_n_iterations=max_n_iterations, n_iterations_between_checks = n_iterations_between_checks )
                potential = ElectricPotentialArray(fssrb)
                new_inds = check_for_refinement(potential, grid, abs.(refinement_limits .* fssrb.bias_voltage), min_grid_spacing)
            else
                break
            end
        end
    end

    pointtypes::Array{PointType, 3} = PointTypeArray(fssrb)
    # min_or_max_potential_for_bubble_marking::T = fssrb.bulk_is_ptype ? fssrb.minimum_applied_potential + 10 : fssrb.maximum_applied_potential - 10
    @inbounds for i in eachindex(pointtypes)
        if (pointtypes[i] & update_bit == 0)
            pointtypes[i] = PointType(0)
        else
            # if (pointtypes[i] & undepleted_bit > 0) # Deactivated bubble feature for now
            #     if fssrb.bulk_is_ptype
            #         if (min_or_max_potential_for_bubble_marking < potential[i] < fssrb.maximum_applied_potential)
            #             pointtypes[i] += bubble_bit
            #         end
            #     else
            #         if (fssrb.minimum_applied_potential < potential[i] < min_or_max_potential_for_bubble_marking)
            #             pointtypes[i] += bubble_bit
            #         end
            #     end
            # end
            if (pointtypes[i] & pn_junction_bit == 0)
                if pointtypes[i] & undepleted_bit > 0 pointtypes[i] -= undepleted_bit end
                # if pointtypes[i] & bubble_bit > 0 pointtypes[i] -= bubble_bit end
            end
        end
    end

    # Checks
    max::T = maximum(potential)
    min::T = minimum(potential)
    if max > fssrb.maximum_applied_potential
        @warn "Detector not fully depleted at this bias voltage ($(fssrb.bias_voltage) V). At least on grid point has a higher potential value ($(max) V) than the maximum applied potential ($(fssrb.maximum_applied_potential) V). This should not be. However, small overshoots might be due to over relaxation and/or not full convergence."
    end
    if min < fssrb.minimum_applied_potential
        @warn "Detector not fully depleted at this bias voltage ($(fssrb.bias_voltage) V). At least on grid point has a smaller potential value ($(min) V) than the minimum applied potential ($(fssrb.minimum_applied_potential) V). This should not be. However, small overshoots might be due to over relaxation and/or not full convergence."
    end

    return PotentialSimulationSetup{T, N, S}( Grid(fssrb), potential, pointtypes, ChargeDensityArray(fssrb), DielektrikumDistributionArray(fssrb)  )
end

"""
    ElectricPotential(detector::SolidStateDetector{T}; kwargs...) where {T}

Calls `calculate_electric_potential(detector; kwargs...)` end extracts only the electric potential and returns it.
"""
function ElectricPotential(detector::SolidStateDetector{T}; kwargs...) where {T}
    setup::PotentialSimulationSetup{T} = calculate_electric_potential(detector; kwargs...)
    return ElectricPotential(setup)
end

function calculate_electric_potential(  detector::SolidStateDetector{T, :cartesian};
                                        init_grid_spacing::Vector{<:Real} = [0.001, 0.001, 0.001], # 5mm each
                                        grid::Grid{T, N, S} = Grid(detector, init_grid_spacing=init_grid_spacing),
                                        convergence_limit::Real = 5e-6,
                                        max_refinements::Int = 3,
                                        refinement_limits::Vector{<:Real} = [1e-4, 1e-4, 1e-4],
                                        min_grid_spacing::Vector{<:Real} = [1e-5, 1e-5, 1e-5],  # mm, mm, mm
                                        depletion_handling::Bool = false,
                                        use_nthreads::Int = Base.Threads.nthreads(),
                                        sor_consts::Vector{<:Real}=[1.4],
                                        max_n_iterations::Int=10000,
                                        verbose::Bool=true,
                                        ) where {T, N, S}
    refinement_limits::Vector{T} = T.(refinement_limits)
    min_grid_spacing::Vector{T} = T.(min_grid_spacing)
    sor_consts::Vector{T} = T.(sor_consts)
    refine::Bool = max_refinements > 0 ? true : false

    if use_nthreads > Base.Threads.nthreads()
        use_nthreads = Base.Threads.nthreads();
        @warn "`use_nthreads` was set to `1`. The environment variable `JULIA_NUM_THREADS` must be set appropriately before the julia session is started."
    end

    fssrb::PotentialSimulationSetupRB{T, 3, 4, :cartesian} = PotentialSimulationSetupRB(detector, grid);

    if verbose
        println("Electric Potential Calculation\n",
                #"Bulk type: $(detector.bulk_type)\n",
                "Bias voltage: $(fssrb.bias_voltage) V\n",
                "Precision: $T\n",
                "Convergence limit: $convergence_limit => $(round(abs(fssrb.bias_voltage * convergence_limit), sigdigits=2)) V\n",
                "Threads: $use_nthreads\n",
                "Coordinate system: Cartesian\n",
                "Initial grid dimension: $(size(grid))\n",
                "Refine? -> $refine\n",
                "Refinement parameters:\n",
                "\tmaximum number of refinements: $(max_refinements)\n",
                "\tminimum grid spacing:\n",
                "\t\tx: $(min_grid_spacing[1]) m\n",
                "\t\ty: $(min_grid_spacing[2]) m\n",
                "\t\tz: $(min_grid_spacing[3]) m\n",
                "\tRefinement limits:\n",
                "\t\tx: $(refinement_limits[1]) -> $(round(abs(fssrb.bias_voltage * refinement_limits[1]), sigdigits=2)) V\n",
                "\t\ty: $(refinement_limits[2]) -> $(round(abs(fssrb.bias_voltage * refinement_limits[2]), sigdigits=2)) V\n",
                "\t\tz: $(refinement_limits[3]) -> $(round(abs(fssrb.bias_voltage * refinement_limits[3]), sigdigits=2)) V\n",
                ""
        )
    end

    n_iterations_between_checks::Int = 500
    update_till_convergence!(   fssrb, convergence_limit, fssrb.bias_voltage,
                                depletion_handling = Val{depletion_handling}(),
                                use_nthreads=use_nthreads, max_n_iterations=max_n_iterations, n_iterations_between_checks = n_iterations_between_checks )

    potential::Array{T, 3} = ElectricPotentialArray(fssrb)

    refinement_counter::Int = 0
    if refine
        new_inds = check_for_refinement(potential, grid, abs.(refinement_limits .* fssrb.bias_voltage), min_grid_spacing)
        while !isempty(new_inds[1]) || !isempty(new_inds[2]) || !isempty(new_inds[3])
            refinement_counter += 1
            if refinement_counter <= max_refinements
                potential, grid = add_points_and_interpolate(potential, grid, new_inds...)
                if verbose println("Refinement $refinement_counter:\tNew grid size: $(size(potential))") end
                fssrb = PotentialSimulationSetupRB(detector, grid, potential)
                n_iterations_between_checks = div(n_iterations_between_checks, 2)
                if n_iterations_between_checks < 50 n_iterations_between_checks = 50 end
                update_till_convergence!(fssrb, convergence_limit, fssrb.bias_voltage, depletion_handling = Val{depletion_handling}(), use_nthreads=use_nthreads, max_n_iterations=max_n_iterations, n_iterations_between_checks = n_iterations_between_checks )
                potential = ElectricPotentialArray(fssrb)
                new_inds = check_for_refinement(potential, grid, abs.(refinement_limits .* fssrb.bias_voltage), min_grid_spacing)
            else
                break
            end
        end
    end

    pointtypes::Array{PointType, 3} = PointTypeArray(fssrb)
    # min_or_max_potential_for_bubble_marking::T = fssrb.bulk_is_ptype ? fssrb.minimum_applied_potential + 10 : fssrb.maximum_applied_potential - 10
    @inbounds for i in eachindex(pointtypes)
        if (pointtypes[i] & update_bit == 0)
            pointtypes[i] = PointType(0)
        else
            # if (pointtypes[i] & undepleted_bit > 0) # Deactivated bubble feature for now
            #     if fssrb.bulk_is_ptype
            #         if (min_or_max_potential_for_bubble_marking < potential[i] < fssrb.maximum_applied_potential)
            #             pointtypes[i] += bubble_bit
            #         end
            #     else
            #         if (fssrb.minimum_applied_potential < potential[i] < min_or_max_potential_for_bubble_marking)
            #             pointtypes[i] += bubble_bit
            #         end
            #     end
            # end
            if (pointtypes[i] & pn_junction_bit == 0)
                if pointtypes[i] & undepleted_bit > 0 pointtypes[i] -= undepleted_bit end
                # if pointtypes[i] & bubble_bit > 0 pointtypes[i] -= bubble_bit end
            end
        end
    end

    # Checks
    max::T = maximum(potential)
    min::T = minimum(potential)
    if max > fssrb.maximum_applied_potential
        @warn "Detector not fully depleted at this bias voltage ($(fssrb.bias_voltage) V). At least on grid point has a higher potential value ($(max) V) than the maximum applied potential ($(fssrb.maximum_applied_potential) V). This should not be. However, small overshoots might be due to over relaxation and/or not full convergence."
    end
    if min < fssrb.minimum_applied_potential
        @warn "Detector not fully depleted at this bias voltage ($(fssrb.bias_voltage) V). At least on grid point has a smaller potential value ($(min) V) than the minimum applied potential ($(fssrb.minimum_applied_potential) V). This should not be. However, small overshoots might be due to over relaxation and/or not full convergence."
    end

    return PotentialSimulationSetup{T, N, S}( Grid(fssrb), potential, pointtypes, ChargeDensityArray(fssrb), DielektrikumDistributionArray(fssrb)  )
end


text:
abstract type PositionRecord end

mutable struct FullJob
    history::Vector{PositionRecord}
end

struct InTransit <: PositionRecord
    time::Float64 #time transit started
    end_time::Float64 #time transit will end
end

struct InQueue <: PositionRecord
    time::Float64 #time entered queue
    queue::Int #queue index
end

struct LeaveSystem <: PositionRecord
    time::Float64 #time left system
end

FullJob() = FullJob([])
FullJob(entry::PositionRecord) = FullJob([entry])

""" Returns the time the job entered the system """
entry_time(job::FullJob)::Float64 = first(job.history).time

""" Returns the time the job exited the system or time of last event if still in the system """
exit_time(job::FullJob)::Float64 = last(job.history).time

""" Returns the duration of time the job has spent in the system """
duration(job::FullJob)::Float64 = exit_time(job) - entry_time(job)

function push!(job::FullJob, pos::PositionRecord)
    push!(job.history, pos)
end

isless(p1::PositionRecord, p2::PositionRecord) = p1.time < p2.time
isless(p1::InTransit, p2::InTransit) = p1.end_time < p2.end_time
isless(j1::FullJob, j2::FullJob) = last(j1.history) < last(j2.history)

function show(io::IO, job::FullJob)
    print(io, "FullJob(Entry: $(entry_time(job)), Pos: $(last(job.history)), Duration: $(duration(job)))")
end

function show(io::IO, pos::InTransit)
    print(io, "(T: $(round(pos.time, digits = 4)))")
end

function show(io::IO, pos::InQueue)
    print(io, "(Q$(pos.queue): $(round(pos.time, digits = 4)))")
end

function show(io::IO, pos::LeaveSystem)
    print(io, "(E: $(round(pos.time, digits = 4)))")
end

function show(io::IO, history::Vector{PositionRecord})
    length(history) == 0 && print(io, "[]") && return
    print(io, "[$(first(history))")
    for i in 2:length(history)
        print(io, "-> $(history[i])")
    end
    print(io, "]")
end


text:
module FluxUtil

using ..Reproduce
import Flux

function opt_settings!(as::ArgParseSettings)
    @add_arg_table as begin
        "--opt"
        arg_type=String
        required=true
        "--optparams"
        arg_type=Float64
        required=true
    end
end

function get_optimizer(parsed::Dict)
    kt = keytype(parsed)
    get_optimizer(parsed[kt("opt")], parsed[kt("optparams")])
end

function get_optimizer(opt_string::AbstractString, params)
    opt_func = getproperty(Flux, Symbol(opt_string))
    return opt_func(params...)
end

end


text:
#-------------monte_carlo.jl---------------------------------------------------#
#
# Purpose: A simple monte carlo script to integrate x^2 in a given range
#
#------------------------------------------------------------------------------#

# A function to integrate x^2 in a given range with num number of points
function monte_carlo(range::Float64, num::Int64)
    count = 0
    for i = 1:num
        x_pos = rand() * range - range * 0.5
        y_pos = rand() * 0.25 * range * range

        if (y_pos < x_pos*x_pos)
            count += 1
        end
    end

    #println("The integral of x^2 between ", -range*0.5, " and ",
            #range*0.5, " is ", (count / num) *0.25*range*range*range)
end

monte_carlo(6.0, 10000)

max_stars_repo_name: vchuravy/skillpill-julia

max_issues_repo_name: oist/skillpill-julia

max_forks_repo_name: oist/skillpill-julia


text:
"""
Star log likelihood maker.  Creates a star log prob function that is
parameterized by a flat vector

  - th : [log_fluxes; unconstrained_pos]

Args:
  imgs: Array of observed data Images with the .pixel field
  pos0: Initial location of the source in ra/dec
  pos_delta: (optional) determines how much ra/dec we allow the sampler
    to drift away from pos0
"""
function make_star_loglike(imgs::Array{Image},
                           pos0::Array{Float64, 1};
                           pos_delta::Array{Float64, 1}=[1., 1.])

    # create a function that constrains the pixel location to be within
    # pos0 +- [pos_delta]
    constrain_pos, unconstrain_pos =
        make_position_transformations(pos0, pos_delta)

    function star_loglike(th::Array{Float64, 1})

      # unpack location and log fluxes (smushed)
      lnfluxes, unc_pos = th[1:5], th[6:end]
      pos = constrain_pos(unc_pos)

      ll = 0.
      for img in imgs

          # sky pixel intensity (sky image)
          epsilon    = img.sky[1, 1]
          iota       = img.nelec_per_nmgy[1]
          sky_pixels = [epsilon * iota for h=1:img.H, w=1:img.W]

          # create and cache unit flux src image
          #src_pixels = [0. for h=1:img.H, w=1:img.W]
          src_pixels = zeros(img.H, img.W)
          write_star_unit_flux(img, pos, src_pixels)

          # compute model image
          rates = sky_pixels .+ (exp(lnfluxes[img.b])*src_pixels)

          # sum per-pixel likelihood contribution
          H, W = size(rates)
          for h in 1:H, w in 1:W
              ll += poisson_lnpdf(img.pixels[h,w], rates[h,w])
          end

      end
    return ll
    end

    return star_loglike, constrain_pos, unconstrain_pos
end

"""
Galaxy log likelihood maker.  Creates a galaxy log prob function that is
parameterized by a flat vector

  - th : [log_fluxes, unconstrained_pos, unconstrained_gal_shape]

Args:
  imgs: Array of observed data Images with the .pixel field
  pos0: Initial location of the source in ra/dec
  pos_delta: (optional) determines how much ra/dec we allow the sampler
    to drift away from pos0
"""
function make_gal_loglike(imgs::Array{Image},
                          pos0::Array{Float64, 1};
                          pos_delta::Array{Float64, 1}=[1., 1.])

    # create a function that constrains the pixel location to be within
    # pos0 +- [pos_delta]
    constrain_pos, unconstrain_pos =
        make_position_transformations(pos0, pos_delta)

    function gal_loglike(th::Array{Float64, 1})

        # unpack location and log fluxes (smushed)
        lnfluxes, unc_pos, ushape = th[1:5], th[6:7], th[8:end]
        pos   = constrain_pos(unc_pos)
        shape = Model.constrain_gal_shape(ushape)

        gal_frac_dev, gal_ab, gal_angle, gal_scale =
            Model.constrain_gal_shape(ushape)

        ll = 0.
        for img in imgs
            # sky pixel intensity (sky image)
            epsilon    = img.sky[1, 1]
            iota       = img.nelec_per_nmgy[1]
            sky_pixels = [epsilon * iota for h=1:img.H, w=1:img.W]

            # create and cache unit flux src image
            src_pixels = [0. for h=1:img.H, w=1:img.W]
            write_galaxy_unit_flux(img, pos, gal_frac_dev,
                               gal_ab, gal_angle, gal_scale, src_pixels)

            # compute model image
            rates = sky_pixels .+ exp(lnfluxes[img.b])*src_pixels

            # sum per-pixel likelihood contribution
            H, W = size(rates)
            for h in 1:H, w in 1:W
                ll += poisson_lnpdf(img.pixels[h,w], rates[h,w])
            end
        end

    return ll
    end

    return gal_loglike, constrain_pos, unconstrain_pos
end

function make_position_transformations(pos0::Array{Float64, 1},
                                       pos_delta::Array{Float64, 1})

    # create a function that constrains the pixel location to be within
    # pos0 +- [pos_delta]
    pos_lo = pos0 .- pos_delta
    pos_hi = pos0 .+ pos_delta
    pos_range = pos_hi .- pos_lo

    function constrain_pos(upos::Array{Float64, 1})
        unit = 1 ./ (1. + exp.(-upos))
        cpos  = (unit.*pos_range) .+ pos_lo
        return cpos
    end

    function unconstrain_pos(cpos::Array{Float64, 1})
        unit = (cpos .- pos_lo) ./ pos_range
        upos = log.(unit) .- log.(1 .- unit)
        return upos
    end

    return constrain_pos, unconstrain_pos
end

#####################################
# simple log prob implementation    #
#####################################

"""
given a single image, create a function that returns log likelihood of image
as a function of the log flux = [lnu, lnr, lng, lni, lnz]
"""
function make_single_image_logflux_loglike(img0::Image,
                                           pos::Array{Float64, 1};
                                           is_star               = true,
                                           gal_frac_dev::Float64 = .5,
                                           gal_ab::Float64       = .7,
                                           gal_angle::Float64    = .79,
                                           gal_scale::Float64    = 4.)
    # sky pixel intensity
    epsilon    = img0.sky[1, 1]
    iota       = img0.nelec_per_nmgy[1]
    sky_pixels = [epsilon * iota for h=1:img0.H, w=1:img0.W]

    # create and cache unit flux src image
    src_pixels = [0. for h=1:img0.H, w=1:img0.W]
    if is_star
        write_star_unit_flux(img0, pos, src_pixels)
    else
        write_galaxy_unit_flux(img0, pos, gal_frac_dev,
                               gal_ab, gal_angle, gal_scale, src_pixels)
    end

    # create log like handle and return
    function lnflux_loglike(lnflux)

        # compute model image
        rates = sky_pixels .+ exp(lnflux)*src_pixels

        # sum per-pixel likelihood contribution
        ll = 0.
        H, W = size(rates)
        for h in 1:H, w in 1:W
            ll += poisson_lnpdf(img0.pixels[h,w], rates[h,w])
        end

        return ll
    end

    return lnflux_loglike
end

"""
create vectorized star log flux log like function with multiple images
"""
function make_logflux_loglike(imgs::Array{Image},
                              pos::Array{Float64, 1};
                              is_star = true,
                              gal_frac_dev::Float64 = .5,
                              gal_ab::Float64       = .7,
                              gal_angle::Float64    = .79,
                              gal_scale::Float64    = 4.)

    # create per image scalar log flux function (caches unit flux pixel image)
    if is_star
        scalar_funs = [make_single_image_logflux_loglike(img, pos;
                          is_star=true)
                       for img in imgs]
    else
        scalar_funs = [make_single_image_logflux_loglike(img, pos;
                          is_star=false,
                          gal_frac_dev=gal_frac_dev,
                          gal_ab=gal_ab,
                          gal_angle=gal_angle,
                          gal_scale=gal_scale)
                       for img in imgs]
    end

    function lnflux_loglike(lnflux::Vector{Float64})
        # for each of those image/function pairs, compute flux-specific ll
        ll = 0.
        for i in 1:length(imgs)
            b   = imgs[i].b
            ll += scalar_funs[i](lnflux[b])
        end
        return ll
    end

    return lnflux_loglike
end

#############################
# Priors and prior samplers #
#############################

const pp = Model.load_prior()

"""
prior over log fluxes (for stars and galaxies)
"""
function logflux_logprior(lnfluxes::Vector{Float64}; is_star::Bool=true)
    # toggle star / galaxy
    type_i = is_star ? 1 : 2

    # convert to color space for prior calc
    lnr, colors = logfluxes_to_colors(lnfluxes)

    # compute brightness distribution
    llr = logpdf(Normal(pp.flux_mean[type_i], 2*pp.flux_var[type_i]), lnr)

    # compute color likelihoods (mixture model --- computes all component lls)
    nc, nk, ns = size(pp.color_mean)
    llk = [logpdf(MvNormal(pp.color_mean[:,k,type_i],
                           pp.color_cov[:,:,k,type_i]), colors)
           for k in 1:nk]
    lnpik = log.(pp.k[:, type_i])
    llc   = Model.logsumexp(llk .+ lnpik)

    # add and return
    return llr + llc
end

"""
More numerically stable converter from log fluxes to colors
"""
function logfluxes_to_colors(lnfluxes::Vector{Float64})
    lnr = lnfluxes[3]
    colors = Vector{Float64}(4)
    colors[1] = lnfluxes[2] - lnfluxes[1]
    colors[2] = lnfluxes[3] - lnfluxes[2]
    colors[3] = lnfluxes[4] - lnfluxes[3]
    colors[4] = lnfluxes[5] - lnfluxes[4]
    return lnr, colors
end

"""
sample fluxes from prior
"""
function sample_logfluxes(; is_star=true, lnr=nothing)
    # toggle star / galaxy
    type_i = is_star ? 1 : 2

    # sample log r
    if lnr == nothing
        lnr = rand(Normal(pp.flux_mean[type_i], 2*pp.flux_var[type_i]))
    end

    # sample colors
    k = rand(Distributions.Categorical(pp.k[:, type_i]))
    c = rand(MvNormal(pp.color_mean[:,k,type_i], pp.color_cov[:,:,k,type_i]))

    # convert to fluxes
    lnfluxes = log.(Model.colors_to_fluxes(lnr, c))
    return lnfluxes
end

"""
sample colors from prior
"""
function sample_colors(; is_star=true)
    # toggle star / galaxy
    type_i = is_star ? 1 : 2
    # sample colors
    k = rand(Distributions.Categorical(pp.k[:, type_i]))
    c = rand(MvNormal(pp.color_mean[:,k,type_i], pp.color_cov[:,:,k,type_i]))
    return c
end

function sample_logr(; is_star=true)
    type_i = is_star ? 1 : 2
    lnr = rand(Normal(pp.flux_mean[type_i], 2*pp.flux_var[type_i]))
    return lnr
end

function sample_fluxes(i::Int, r_s)
    k_s = rand(Distributions.Categorical(pp.k[i]))
    c_s = rand(Distributions.MvNormal(pp.color[i][:, k_s], pp.color[i][:, :, k_s]))

    l_s = Vector{Float64}(5)
    l_s[3] = r_s
    l_s[4] = l_s[3] * exp(c_s[3])
    l_s[5] = l_s[4] * exp(c_s[4])
    l_s[2] = l_s[3] / exp(c_s[2])
    l_s[1] = l_s[2] / exp(c_s[1])
    l_s
end


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/" "${BASH_SOURCE[0]}" "$@" -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

println("\033[1;34m===>\033[0;38m\033[1;38m\tLoading required libraries (it will take a moment to precompile if it is your first time doing this)...\033[0;38m")

include(joinpath(dirname(dirname(@__FILE__)), "src", "FaceDetection.jl"))

using .FaceDetection
const FD = FaceDetection
using Images: imresize
using LightXML

println("...done")

function main()
    xml_data = parse_file(joinpath(dirname(dirname(@__FILE__)), "data", "haarcascades", "haarcascade_frontalface_default.xml")) # ../data/haarcascades/haarcascade_frontalface_default.xml

	for c in child_nodes(root(xml_data))  # c is an instance of XMLNode
    if is_elementnode(c)
        e = XMLElement(c)  # this makes an XMLElement instance
        println(name(e))
    end
end
end

main()


text:
module structs

    include("io.jl")

    #using Logging
    #Logging.configure(level = Logging.CRITICAL)
    using JDBC, JavaCall, DataFrames

    #init JDBC
    for driver_family in io.driver_families()
        JavaCall.addClassPath(joinpath(dirname(@__FILE__), "drivers", "$driver_family.jar"))
    end
    JDBC.init()

    mutable struct con
        #values
        result::Dict

        #hidden values
        _driver::Any
        _status::String
        _con::Any
        _cursor::Any
        _con_info::Dict

        #functions
        execute::Function
        close::Function

        function con(con_info)
            instance = new(con_info)

            #set con object
            instance._driver = con_info["driver"] ; con_info = delete!(con_info, "driver")
            try
                instance._con = JDBC.DriverManager.getConnection(con_info["url"]) ; con_info = delete!(con_info, "url")
                instance._cursor = JDBC.createStatement(instance._con)
                instance._status = "success"
                instance._con_info = con_info
            catch error
                instance._con = nothing
                instance._status = "failed"
                instance._con_info = ""
            end

            instance.execute = function(query)
                datas, columns, result = [], [], []
                column_set = false

                rs = JDBC.executeQuery(instance._cursor, query)
                for row in JDBC.JDBCRowIterator(rs)
                    append!(datas, [[Any(item) for item in row]])
                end

                meta_data = JDBC.getTableMetaData(JDBC.executeQuery(instance._cursor, query))
                if !column_set
                    columns = [item[1] for item in meta_data]
                    column_set = true
                end

                for data in datas
                    item = Dict()
                    for idx in eachindex(columns)
                        item[columns[idx]] = data[idx]
                    end

                    append!(result, [item])
                end

                return result
            end

            instance.close = function()
                JDBC.close(instance._con)
            end

            return instance
        end
    end
    function Base.show(io::Base.IO, instance::con)
        args_string = join([join([sss for sss in split(ss, "=")], " = ") for ss in split(instance._con_info["args"], "&") if ss != ""], "\n")
        print(io, string(
            "driver : ", instance._driver, "\n",
            "status : ", instance._status,
            args_string != "" ? string("\n", args_string) : ""
        ))
    end
end


text:
"""
    NonparametricBootSimulator

The NonparametricBootSimulator type holds the parameters and storages for simulating the multinomial counts
"""
struct NonparametricBootSimulator{T <: LinearAlgebra.BlasReal}
    # estimated obtained from the subset
    β_subset::Vector{T}
    Σ_subset::Matrix{T}
    σ²_subset::Vector{T}
    # for simulating multinomial counts
    ns::Vector{Int64}
    mult_prob::Vector{T}
    mult_dist::Distributions.Multinomial
end

"""
    NonparametricBootSimulator(m)

The constructor for the `NonparametricBootSimulator` type.
"""
function NonparametricBootSimulator(
    m::blblmmModel{T},
    ) where T <: BlasReal
    # *_subset will not change throughout bootstrap iterations
    β_subset = similar(m.β)
    copyto!(β_subset, m.β)
    Σ_subset = similar(m.Σ)
    copyto!(Σ_subset, m.Σ)
    σ²_subset = similar(m.σ²)
    copyto!(σ²_subset, m.σ²)

    # initialize a vector for storing multinomial counts
    ns = zeros(Int64, m.b)
    mult_prob = ones(m.b) / m.b
    mult_dist = Multinomial(m.N, mult_prob)
    NonparametricBootSimulator(β_subset, Σ_subset, σ²_subset, ns, mult_prob, mult_dist)
end

"""
    ParametricBootSimulator

The `ParametricBootSimulator` type holds the parameters and storages for simulating the response and the multinomial counts
"""
struct ParametricBootSimulator{T <: LinearAlgebra.BlasReal}
    # estimated obtained from the subset
    β_subset::Vector{T}
    Σ_subset::Matrix{T}
    ΣL_subset::Matrix{T}
    σ²_subset::Vector{T}
    # for simulating the response y
    Xβ::Vector{Vector{T}}
    storage_q::Vector{T}
    re_storage::Vector{T}
    # re_dist::MvNormal
    # for simulating multinomial counts
    ns::Vector{Int64}
    mult_prob::Vector{T}
    mult_dist::Distributions.Multinomial
end

"""
    ParametricBootSimulator(m)

The constructor for the `ParametricBootSimulator` type.
"""
function ParametricBootSimulator(
    m::blblmmModel{T},
    ) where T <: BlasReal
    # *_subset will not change throughout bootstrap iterations
    β_subset = similar(m.β)
    copyto!(β_subset, m.β)
    Σ_subset = similar(m.Σ)
    copyto!(Σ_subset, m.Σ)
    ΣL_subset = similar(m.ΣL) #Matrix(cholesky(Symmetric(Σ_subset)).L)
    copyto!(ΣL_subset, m.ΣL)
    σ²_subset = similar(m.σ²)
    copyto!(σ²_subset, m.σ²)

    # Since xtβ is constant throughout simulation, we pre-compute it
    Xβ = Vector{Vector{T}}(undef, m.b)
    @inbounds @views for i in 1:m.b
        Xβ[i] = Vector{T}(undef, m.data[i].n)
        BLAS.gemv!('N', T(1), m.data[i].X, β_subset, T(0), Xβ[i])
    end

    # initialize a vector for storing α
    re_storage = Vector{T}(undef, m.q)
    storage_q = Vector{T}(undef, m.q)
    # distribution of the random effect
    # re_dist = MvNormal(zeros(m.q), Σ_subset)
    # initialize a vector for storing multinomial counts
    ns = zeros(Int64, m.b)
    mult_prob = ones(m.b) / m.b
    mult_dist = Multinomial(m.N, mult_prob)
    ParametricBootSimulator(
        β_subset, Σ_subset, ΣL_subset, σ²_subset,
        Xβ, storage_q, re_storage,
        ns, mult_prob, mult_dist
        )
end

function simulate!(
    rng::Random.AbstractRNG,
    m::MixedModelsBLB.blblmmModel{T},
    simulator::ParametricBootSimulator{T}
    ) where T<: LinearAlgebra.BlasReal
    σ = sqrt(simulator.σ²_subset[1])
    @inbounds @views for bidx = 1:m.b
        randn!(rng, m.data[bidx].y) # y = standard normal error
        BLAS.axpby!(T(1), simulator.Xβ[bidx], σ, m.data[bidx].y) # y = Xβ + σ * standard normal error
        # simulate random effect: ΣL * standard normal
        randn!(rng, simulator.storage_q)
        BLAS.gemv!('N', 1., simulator.ΣL_subset, simulator.storage_q, 0., simulator.re_storage)
        BLAS.gemv!('N', T(1), m.data[bidx].Z, simulator.re_storage, T(1), m.data[bidx].y) # y = Xβ + Zα + error

        # Update quantities related to y in blblmmObs
        m.data[bidx].yty[1] = dot(m.data[bidx].y, m.data[bidx].y)
        BLAS.gemv!('T', T(1), m.data[bidx].X, m.data[bidx].y, T(0), m.data[bidx].xty)
        BLAS.gemv!('T', T(1), m.data[bidx].Z, m.data[bidx].y, T(0), m.data[bidx].zty)
        # copyto!(m.data[bidx].xty, transpose(m.data[bidx].X) * m.data[bidx].y)
        # copyto!(m.data[bidx].zty, transpose(m.data[bidx].Z) * m.data[bidx].y)
    end
end