text:
abstract type AbstractMTable{I,Tw} <: MultiplicativeStructure{Tw,I} end

Base.size(mt::AbstractMTable) = size(mt.table)

_check(mt::AbstractMTable) = _check(mt.table, basis(mt), _istwisted(mt))

function _check(product_matrix, basis, twisted::Bool)
    idx = findfirst(iszero, product_matrix)
    if idx != nothing
        i, j = Tuple(idx)
        x = (twisted ? star(basis[i]) : basis[i])
        throw(
            ProductNotDefined(
                i,
                j,
                "$x · $(basis[j]) = $(_product(Val(twisted), x, basis[j]))",
            ),
        )
    end
    return true
end

_iscached(mt::AbstractMTable, i, j) = !iszero(mt.table[i, j])

## MTables

struct MTable{I,Twisted,M<:AbstractMatrix{I}} <: AbstractMTable{I,Twisted}
    table::M
end

MTable{Tw}(mt::AbstractMatrix{<:Integer}) where {Tw} = MTable{eltype(mt),Tw,typeof(mt)}(mt)

MTable(b::AbstractBasis; table_size) = MTable{false}(b; table_size = table_size)

function MTable{Tw}(basis::AbstractBasis; table_size) where {Tw}
    @assert length(table_size) == 2
    @assert 1 <= first(table_size) <= length(basis)
    @assert 1 <= last(table_size) <= length(basis)

    table = zeros(SparseArrays.indtype(basis), table_size)

    complete!(table, basis, Val(Tw))

    _check(table, basis, Tw)

    return MTable{Tw}(table)
end

for twisted in (:true, :false)
    @eval begin
        function complete!(table, basis, v::Val{$twisted})
            Threads.@threads for j in 1:size(table, 2)
                y = basis[j]
                for i in 1:size(table, 1)
                    table[i, j] = basis[_product(v, basis[i], y)]
                end
            end
            return table
        end
    end
end

basis(mt::MTable) = throw("No basis is defined for a simple $(typeof(mt))")

Base.@propagate_inbounds function Base.getindex(m::MTable, i::Integer, j::Integer)
    @boundscheck checkbounds(Bool, m, i, j) ||
                 throw(ProductNotDefined(i, j, "out of Mtable bounds"))
    @boundscheck iszero(m.table[i, j]) && throw(ProductNotDefined(i, j))
    return m.table[i, j]
end

## CachedMTables

struct CachedMTable{T,I,B<:Basis{T,I},M,Twisted} <: AbstractMTable{I,Twisted}
    basis::B
    table::M
end

CachedMTable(basis::AbstractBasis; table_size) =
    CachedMTable{false}(basis; table_size = table_size)

function CachedMTable{Tw}(basis::AbstractBasis{T,I}; table_size) where {Tw,T,I}
    return CachedMTable{Tw}(basis, zeros(I, table_size))
end

CachedMTable(basis::AbstractBasis, mt::AbstractMatrix{<:Integer}) =
    CachedMTable{false}(basis, mt)

function CachedMTable{Tw}(
    basis::AbstractBasis{T,I},
    mt::AbstractMatrix{<:Integer},
) where {Tw,T,I}
    return CachedMTable{T,I,typeof(basis),typeof(mt),Tw}(basis, mt)
end

basis(m::CachedMTable) = m.basis

Base.@propagate_inbounds function Base.getindex(cmt::CachedMTable, i::Integer, j::Integer)
    cache!(cmt, i, j)
    return cmt.table[i, j]
end

Base.@propagate_inbounds function cache!(cmt::CachedMTable, i::Integer, j::Integer)
    @boundscheck checkbounds(Bool, cmt, i, j) ||
                 throw(ProductNotDefined(i, j, "out of Mtable bounds"))
    if !_iscached(cmt, i, j)
        b = basis(cmt)
        g, h = b[i], b[j]
        gh = _product(cmt, g, h)
        gh in b || throw(ProductNotDefined(i, j, "$g · $h = $gh"))
        cmt.table[i, j] = b[gh]
    end
    return cmt
end

function cache!(
    cmt::CachedMTable,
    suppX::AbstractVector{<:Integer},
    suppY::AbstractVector{<:Integer},
)
    Threads.@threads for j in suppY
        for i in suppX
            if !_iscached(cmt, i, j)
                cache!(cmt, i, j)
            end
        end
    end
    return cmt
end

complete!(cmt::CachedMTable) = cache!(cmt, 1:size(cmt, 1), 1:size(cmt, 2))


text:
module GaussianDistributions
# there are fundamental operations with Gaussian distributions:
# 1. marginals
# 2. conditionals
# 3. linear functionals
# 4. sampling
# 5. optimizing hyper-parameters of Gaussian Processes
# this package provides convenient and efficient implementations for these operations
function marginal end
function conditional end

using LinearAlgebra
using LinearAlgebra: checksquare

using Base.Threads
using StaticArrays
using FillArrays
using Statistics
import StatsBase: sample

using CovarianceFunctions
using CovarianceFunctions: gramian, MercerKernel, MultiKernel, AbstractKernel
Statistics.cov(k::MercerKernel, x::AbstractVector) = gramian(k, x)

using ForwardDiff # necessary for value_gradient in util.jl
using DiffResults

export Gaussian, Univariate, Multivariate, VectorGaussian, 𝒩
export mean, cov, var, std, pdf, cdf, nld, nlml
export marginal, conditional, conditional!, sample

using LazyInverse: inverse, Inverse, pseudoinverse
using LazyLinearAlgebra
using LinearAlgebraExtensions: difference, AbstractMatOrFacOrUni, AbstractMatOrFac
using WoodburyIdentity

abstract type AbstractGaussian{M, S} end
Statistics.mean(G::AbstractGaussian) = G.μ
Statistics.cov(G::AbstractGaussian) = G.Σ # ensure Hermitian-ness if Σ is a matrix?
const tol = 1e-6 # default error tolerance for matrix factorizations and samples
const VecOfVec{T} = AbstractVector{<:AbstractVector{T}}

############################# Gaussian Distribution ############################
struct Gaussian{M, S} <: AbstractGaussian{M, S}
    μ::M
    Σ::S
	function Gaussian(μ::Real, σ²::Real)
		σ² ≥ 0 ? new{typeof(μ), typeof(σ²)}(float(μ), σ²) : error("σ² < 0")
	end
	function Gaussian(μ::AbstractVector, Σ::AbstractMatOrFacOrUni)
		length(μ) == checksquare(Σ) || error("dimensions of μ, Σ not compatible")
		new{typeof(μ), typeof(Σ)}(μ, Σ)
	end
	Gaussian(μ, Σ) = new{typeof(μ), typeof(Σ)}(μ, Σ)
end
const 𝒩 = Gaussian
# 1-argument input constructors, assuming mean zero
Gaussian() = Gaussian(0., 1.)
# passing variances only
Gaussian(σ²::Real) = Gaussian(zero(σ²), σ²)
Gaussian(σ²::AbstractVector) = Gaussian(Diagonal(σ²))
Gaussian(Σ::AbstractMatOrFacOrUni) = Gaussian(Zeros(size(Σ, 1)), Σ)
# passing mean and variances
Gaussian(μ::AbstractVector, σ²::Real) = Gaussian(μ, σ²*I(length(μ)))
Gaussian(μ::AbstractVector, σ²::AbstractVector) = Gaussian(μ, Diagonal(σ²))

zero_mean(x) = zero(eltype(x))
# TODO: make this a mutable static vector
zero_mean(x, d::Int) = zeros(eltype(x), d)
zero_mean(d::Int) = x->zero_mean(x, d)

Gaussian(k) = Gaussian(zero_mean, k)

function Base.length(G::Gaussian)
	if G isa Univariate || G isa Multivariate
		return length(G.μ)
	else
		return Inf # assumed to be process
	end
end

const UnivariateGaussian{M, S} = Gaussian{<:Real, <:Real}
const Univariate = UnivariateGaussian
const Uni = Univariate

# (G::Univariate)(i::Int) = marginal(G, i)
Base.getindex(G::Univariate, i) = marginal(G, i) # bracket-indexing

const MultivariateGaussian{M, S} = Gaussian{<:AbstractVector, <:AbstractMatOrFacOrUni}
const Multivariate = MultivariateGaussian
const Multi = Multivariate

# (G::Multivariate)(i::Int) = marginal(G, i)
Base.getindex(G::Multivariate, i) = marginal(G, i) # bracket-indexing

# Gaussian is assumed to be process-valued, if neither Uni, nor Multi
(G::Gaussian)(x) = Gaussian(mean(G)(x), cov(G)(x, x))

############################# convenience ######################################
Statistics.cov(G::Multivariate, i, j) = cov(G)[i, j] # view?
Statistics.var(G::Union{Uni, Multi}) = diag(cov(G))
Statistics.std(G::Union{Uni, Multi}) = sqrt.(max.(var(G), 0))

# variance (computes diagonal of covariance), should be specialized where
# more efficient computation is possible
Statistics.var(G::Gaussian) = x->max(G.Σ(x, x), 0)
Statistics.std(G::Gaussian) = x->sqrt(var(G)(x))
Statistics.mean(G::Gaussian, x::AbstractVector) = marginal(mean(G), x)
# IDEA: allow for randomized approximation of this function akin to "Scaling Gaussian Process Regression with Derivatives"
# Statistics.var(G::Gaussian, x::AbstractVector) = var(G).(x)
function Statistics.var(G::Gaussian, x::AbstractVector)
	σ² = fill(var(G)(x[1]), length(x))
	@threads for i in 2:length(x) # virtually perfect parallel scaling
		σ²[i] = var(G)(x[i])
	end
	return σ²
end
function Statistics.std(G::Gaussian, x::AbstractVector)
	σ = var(G, x)
	@. σ = sqrt(max(σ, 0))
end

Base.copy(G::Gaussian) = Gaussian(copy(G.μ), copy(G.Σ))
function Base.:(==)(G::Gaussian, H::Gaussian)
	return mean(G) == mean(H) && cov(G) == cov(H)
end
function Base.isapprox(G::Gaussian, H::Gaussian; atol = eps(), rtol = eps())
    isapprox(G.μ, H.μ, atol = atol, rtol = rtol) && isapprox(G.Σ, H.Σ, atol = atol, rtol = rtol)
end

function LinearAlgebra.cholesky(G::Multivariate; tol::Real = tol, check::Bool = false) # low rank via cholesky
	cov(G) isa AbstractMatrix || return G # to prevent calling cholesky on a factorization
	S = Symmetric(cov(G))
    C = cholesky(S, Val(true), tol = tol, check = check)
    if check
		issuccess(C) || throw("pivoted cholesky not successful")
	end
    return Gaussian(mean(G), C)
end
function LinearAlgebra.factorize(G::Multivariate; tol::Real = tol)
	cov(G) isa AbstractMatrix || return G # to prevent calling cholesky on a non-factorizable object
	return cholesky(G, tol = tol)
end

include("util.jl")
include("sample.jl")
include("marginal.jl")
include("conditional.jl")
include("multi.jl")
include("linear.jl")
include("density.jl")
include("optimization.jl")
include("plot.jl")

# non-essential extensions
include("inputnoise.jl")

end # GaussianDistributions

ext: jl


text:
using StatsBase

tests = ["means",
         "scalarstats",
         "counts",
         "ranking",
         "corr",
         "sampling",
         "empirical",
         "misc"]

println("Running tests:")

for t in tests
    tfile = joinpath("test", "$(t).jl")
    println(" * $(tfile) ...")
    include(tfile)
end


text:
"""
    struct ModifyConstraintNotAllowed{F<:AbstractFunction, S<:AbstractSet,
                                      C<:AbstractFunctionModification} <: NotAllowedError
        constraint_index::ConstraintIndex{F, S}
        change::C
        message::String
    end

An error indicating that the constraint modification `change` cannot be applied
to the constraint of index `ci`.
"""
struct ModifyConstraintNotAllowed{F<:AbstractFunction, S<:AbstractSet,
                                  C<:AbstractFunctionModification} <: NotAllowedError
    constraint_index::ConstraintIndex{F, S}
    change::C
    message::String
end
function ModifyConstraintNotAllowed(
    ci::ConstraintIndex{F, S},
    change::AbstractFunctionModification,
    message="") where {F<:AbstractFunction, S<:AbstractSet}
    return ModifyConstraintNotAllowed{F, S, typeof(change)}(ci, change, message)
end
throw_modify_not_allowed(ci::ConstraintIndex, args...) = throw(ModifyConstraintNotAllowed(ci, args...))

operation_name(err::ModifyConstraintNotAllowed{F, S}) where {F, S} = "Modifying the constraints $(err.constraint_index) with $(err.change)"

"""
    struct ModifyObjectiveNotAllowed{C<:AbstractFunctionModification} <: NotAllowedError
        change::C
        message::String
    end

An error indicating that the objective modification `change` cannot be applied
to the objective.
"""
struct ModifyObjectiveNotAllowed{C<:AbstractFunctionModification} <: NotAllowedError
    change::C
    message::String
end
function ModifyObjectiveNotAllowed(change::AbstractFunctionModification)
    return ModifyObjectiveNotAllowed(change, "")
end
throw_modify_not_allowed(::ObjectiveFunction, args...) = throw(ModifyObjectiveNotAllowed(args...))

operation_name(err::ModifyObjectiveNotAllowed) = "Modifying the objective function with $(err.change)"

"""
## Constraint Function

    modify(model::ModelLike, ci::ConstraintIndex, change::AbstractFunctionModification)

Apply the modification specified by `change` to the function of constraint `ci`.

An [`ModifyConstraintNotAllowed`](@ref) error is thrown if modifying
constraints is not supported by the model `model`.

### Examples

```julia
modify(model, ci, ScalarConstantChange(10.0))
```

## Objective Function

    modify(model::ModelLike, ::ObjectiveFunction, change::AbstractFunctionModification)

Apply the modification specified by `change` to the objective function of
`model`. To change the function completely, call `set` instead.

An [`ModifyObjectiveNotAllowed`](@ref) error is thrown if modifying
objectives is not supported by the model `model`.

### Examples

```julia
modify(model, ObjectiveFunction{ScalarAffineFunction{Float64}}(), ScalarConstantChange(10.0))
```
"""
function modify end

function modify(model::ModelLike, ci::ConstraintIndex,
                change::AbstractFunctionModification)
    throw_modify_not_allowed(ci, change)
end

function modify(model::ModelLike, attr::ObjectiveFunction,
                change::AbstractFunctionModification)
    throw_modify_not_allowed(attr, change)
end


text:
#TODO write MPI tests

#TODO ensure result of CSRMatrix(rowMap, colMap, localMatrix, plist) is fill complete

#TODO make testing version of checkInternalState to run during testing

using TypeStability

#### Type Stability Tests ####
method_sigs = [(CSRMatrix{D, G, P, L}, L)
                for (D, G, P, L) in Base.Iterators.product(
                    [Float64, ComplexF32], #Data
                    [UInt64, Int64, UInt32], #GID
                    [UInt8, Int8, UInt32], #PID
                    [UInt32, Int32]) #LID
                ]
report = check_function(JuliaPetra.getLocalRowViewPtr, method_sigs, RegexDict((r"rowValues", Union{Vector, SubArray})))
stability_warn(JuliaPetra.getLocalRowViewPtr, report)
@test is_stable(report)

method_sigs = [(DenseMultiVector{D, G, P, L}, CSRMatrix{D, G, P, L},
                DenseMultiVector{D, G, P, L}, TransposeMode, D, D)
                for (D, G, P, L) in Base.Iterators.product(
                    [Float64, ComplexF32], #Data
                    [UInt64, Int64, UInt32], #GID
                    [UInt8, Int8, UInt32], #PID
                    [UInt32, Int32]) #LID
            ]
report = check_function(localApply, method_sigs, RegexDict((r"(?:#temp#)?@_[0-9]+", Union{Nothing, Tuple{UInt32, UInt32}, Tuple{Int32, Int32}})))
stability_warn(localApply, report)
@test is_stable(report)

#### Serial Tests####s

n = 8
m = 6

Data = Float32
GID = UInt16
PID = Bool
LID = UInt8

commObj = SerialComm{GID, PID, LID}()
rowMap = BlockMap(n, n, commObj)

mat = CSRMatrix{Data}(rowMap, m, STATIC_PROFILE)
@test isa(mat, CSRMatrix{Data, GID, PID, LID})

mat = CSRMatrix{Data}(rowMap, m, STATIC_PROFILE, Dict{Symbol, Any}())
@test isa(mat, CSRMatrix{Data, GID, PID, LID})
@test STATIC_PROFILE == getProfileType(mat)
@test isFillActive(mat)
@test !isFillComplete(mat)
@test isGloballyIndexed(mat)
@test !isLocallyIndexed(mat)
@test rowMap == getRowMap(mat)
@test !hasColMap(mat)
@test n == getGlobalNumRows(mat)
@test n == getLocalNumRows(mat)

@test 0 == getNumEntriesInLocalRow(mat, 2)
@test 0 == getNumEntriesInGlobalRow(mat, 2)
@test 0 == getLocalNumEntries(mat)
@test 0 == getGlobalNumEntries(mat)
@test 0 == getGlobalNumDiags(mat)
@test 0 == getLocalNumDiags(mat)
@test 0 == getGlobalMaxNumRowEntries(mat)
@test 0 == getLocalMaxNumRowEntries(mat)
insertGlobalValues(mat, 2, GID[1, 3, 4], Data[2.5, 6.21, 77])
@test 3 == getNumEntriesInLocalRow(mat, 2)
@test 3 == getNumEntriesInGlobalRow(mat, 2)
@test 3 == getLocalNumEntries(mat)
@test 0 == getLocalNumDiags(mat)
rowInfo = JuliaPetra.getRowInfo(mat.myGraph, LID(2))
@test 3 == rowInfo.numEntries
JuliaPetra.recycleRowInfo(rowInfo)
#skipped many of the global methods because those require re-generating and may not be up to date

row = getGlobalRowCopy(mat, 2)
@test isa(row, Tuple{<: AbstractArray{GID, 1}, <: AbstractArray{Data, 1}})
@test GID[1, 3, 4] == row[1]
@test Data[2.5, 6.21, 77] == row[2]

row = getGlobalRowView(mat, 2)
@test isa(row, Tuple{<: AbstractArray{GID, 1}, <: AbstractArray{Data, 1}})
@test GID[1, 3, 4] == row[1]
@test Data[2.5, 6.21, 77] == row[2]

fillComplete(mat)

row = getLocalRowCopy(mat, 2)
@test isa(row, Tuple{<: AbstractArray{LID, 1}, <: AbstractArray{Data, 1}})
@test LID[1, 2, 3] == row[1]
@test Data[2.5, 6.21, 77] == row[2]

row = getLocalRowView(mat, 2)
@test isa(row, Tuple{<: AbstractArray{LID, 1}, <: AbstractArray{Data, 1}})
@test LID[1, 2, 3] == row[1]
@test Data[2.5, 6.21, 77] == row[2]

#=

getGlobalNumCols(mat::CSRMatrix) = -1#TODO figure out
getLocalNumCols(mat::CSRMatrix) = numCols(mat.localMatrix)
=#

map = BlockMap(2, 2, commObj)

mat = CSRMatrix{Data}(map, 2, STATIC_PROFILE)
insertGlobalValues(mat, 1, GID[1, 2], Data[2, 3])
insertGlobalValues(mat, 2, GID[1, 2], Data[5, 7])
fillComplete(mat)

@test [1, 2, 1, 2] == (mat.myGraph.localIndices1D)
@test (LID[1, 2], Data[2, 3]) == getLocalRowCopy(mat, 1)
@test (LID[1, 2], Data[5, 7]) == getLocalRowCopy(mat, 2)
@test (LID[1, 2], Data[2, 3]) == getLocalRowView(mat, 1)
@test (LID[1, 2], Data[5, 7]) == getLocalRowView(mat, 2)

Y = DenseMultiVector(map, Data[1 0; 0 2])
X = DenseMultiVector(map, Data[2 2; 2 2])

@test Y === apply!(Y, mat, X, NO_TRANS, Data(3), Data(.5))

@test Data[2 2; 2 2] == X.data #ensure X isn't mutated
@test Data[30.5 30; 72 73] == Y.data

Y = DenseMultiVector(map, Data[1 0; 0 2])
X = DenseMultiVector(map, Data[2 2; 2 2]) #ensure bugs in the previous test don't affect this test

@test Y === apply!(Y, mat, X, TRANS, Float32(3), Float32(.5))

@test Data[2 2; 2 2] == X.data #ensure X isn't mutated
@test [42.5 42; 60 61] == Y.data

# Julia's LA API

X = DenseMultiVector(map, Data[2 2; 2 2])

Y = mat * X
@test Data[2 2; 2 2] == X.data
@test Y isa DenseMultiVector
@test getMap(Y) === getMap(X)
@test Data[10 10;24 24] == Y.data

Y = DenseMultiVector(map, Data[1 0; 0 2])
Y = mul!(Y, mat, X)
@test Data[2 2; 2 2] == X.data
@test Data[10 10;24 24] == Y.data


text:
include("Graph.jl")
include("Tools.jl")

using SparseArrays
using Laplacians

function Exact(G, s) # Returns ci, d, p, idc
	T = time()
	L = getSparseL(G)
	W = getW(G)
	avgs = sum(s)/G.n
	s2 = zeros(G.n)
	for i = 1 : G.n
		s2[i] = s[i] - avgs
	end
	z = W * s
	z2 = W * s2
	# calculate C_I(G)
	lz = L * z
	ci = lz' * lz
	# calculate D(G)
	d = z2' * L * z2
	# calculate P(G)
	p = z2' * z2
	# calculate C(G)
	c = z' * z
	# calculate I_dc(G)
	idc = d + c
	# END CALCULATION
	T = time() - T
	return T, ci, d, p, idc
end

function Approx(G, s; eps = 1e-6) # Returns aci, ad, ap, aidc
	T = time()
	IpL = getSparseIpL(G)
	sL = getSparseL(G)
	f = approxchol_sddm(IpL, tol=0.1*eps)
	avgs = sum(s)/G.n
	s2 = zeros(G.n)
	for i = 1 : G.n
		s2[i] = s[i] - avgs
	end
	z = f(s)
	z2 = f(s2)

	# calculate aC_I(G)
	alz = sL * z
	aci = alz' * alz
	# calculate aD(G)
	ad = z2' * sL * z2
	# calculate aP(G)
	ap = z2' * z2
	# calculate aC(G)
	ac = z' * z
	# calculate aI_dc(G)
	aidc = ad + ac
	# END CALCULATION
	T = time() - T
	return T, aci, ad, ap, aidc
end

function doExp(G, s, lg)
	T, ci, d, p, idc = Exact(G, s)
	T2, aci, ad, ap, aidc = Approx(G, s)
	println(lg, "Exact  Time : ", T)
	println(lg, "Approx Time : ", T2)
	println(lg, "ERROR of ci : ", abs(aci-ci)/ci)
	println(lg, "ERROR of d : ", abs(ad-d)/d)
	println(lg, "ERROR of p  : ", abs(ap-p)/p)
	println(lg, "ERROR of idc : ", abs(aidc-idc)/idc)
	println(lg)
end

function doLarge(G, s, lg)
	T2, aci, ad, ap, aidc = Approx(G, s)
	println(lg, "Approx Time : ", T2)
	println(lg)
end


text:
using GenieBDD
using Test

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


text:
using ModelTES, PyPlot

"""`noise_power_demo()`

Plot a noise power spectral density for five TES model systems.
For one, plot the 4 components from the Irwin-Hilton model."""

function noise_power_demo(SI_amp=5e-22)
    models = Dict(
        :LowE => ModelTES.lowEpix(),
        :HighE => ModelTES.highEpix(),
        :Holmes12 => ModelTES.pholmes(12e-9),
        :Holmes24 => ModelTES.pholmes(24e-9),
        :Holmes48 => ModelTES.pholmes(48e-9))
    colors=Dict(:HighE=>"r", :LowE=>"gold",
        :Holmes12=>"purple",:Holmes24=>"blue",:Holmes48=>"green");

    clf()
    freq = logspace(1, 6, 50);
    for m in (:LowE, :HighE, :Holmes12, :Holmes24, :Holmes48)
        model = IrwinHiltonTES(models[m])
        psd,N1,N2,N3,N4 = ModelTES.noisePSD(model, freq, SI_amp)
        loglog(freq,psd, color=colors[m], label = string(m))
        if m==:HighE
            loglog(freq, N1, "--", color=colors[m], label="HighE TES")
            loglog(freq, N2, ":", color=colors[m], label="HighE Load")
            loglog(freq, N3, "-.", color=colors[m], label="HighE TFN")
            loglog(freq, N4, "--", color="orange", label="HighE ampl")
        end
    end
    xlabel("Frequency (Hz)")
    ylabel("Noise power (A\$^2\$/Hz)")
    title("Noise Power Spectrum for several TES models")
    ylim([.3*SI_amp,3e-20+SI_amp])
    legend(loc="best")
    grid(true)
end


text:
using FileIO, Gnuplot, ColorSchemes
using GeometryBasics
# based on https://matplotlib.org/matplotblog/posts/custom-3d-engine/
# download this file
# https://matplotlib.org/matplotblog/posts/custom-3d-engine/bunny.obj
bunny = load("bunny.obj")
ptstemp = bunny.position
npts = length(ptstemp)
points = zeros(npts, 3)
for i in 1:npts
    points[i,:] = ptstemp[i]
end
# normalize into a box
mnpts = minimum(points, dims=1)
mxpts = maximum(points, dims=1)
# all the data needed.
verpts = (points .- (mxpts .+ mnpts)./2)./maximum(mxpts .- mnpts)
triface = decompose(TriangleFace{Int32}, bunny)
# polygons to gnuplot
tmp = Vector{String}()
for i in 1:length(triface)
    polytri = verpts[triface[i],:]
    x, y = polytri[:,1], polytri[:,2] # polytri[:,3]
    push!(tmp, join(string.(x) .* " ".* string.(y), "\n") * "\n")
end

@gp "set size ratio -1" :-
@gp :- tmp "w l notit lc 'black' lw 0.2"
@gp :- "unset xtics; unset ytics; unset border"
save(term = "pngcairo size 800,600", output = "bunnyLines.png")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
__precompile__(true)
module SpatialEcology

using DataFrames
import DataFrames: aggregate
import DataFramesMeta: @with
using SparseArrays
using Statistics
using RandomBooleanMatrices
import RandomNumbers.Xorshifts: Xoroshiro128Plus
import RandomBooleanMatrices: matrixrandomizer, MatrixGenerator

import EcoBase
import EcoBase: asindices, nthings, nplaces, occupancy, richness, nrecords, placenames, thingnames,
        occurring, noccurring, occupied, noccupied, occurrences,
        placeoccurrences, thingoccurrences, cooccurring, places, things,
        indices, coordinates, xcells, ycells, cells, xmin, xmax, ymin, ymax,
        xrange, yrange, xcellsize, ycellsize, cellsize, getcoords

import RecipesBase
import Base: copy, getindex, setindex!, size, show, summary, view, Meta.isexpr
import StatsBase: nquantile
import Distances: pairwise, PreMetric

export SiteData, ComMatrix, Assemblage, coordtype #types and their constructors
export AbstractComMatrix
export nspecies, nsites, occupancy, richness, nrecords, sitenames, speciesnames, coordinates
export occurring, noccurring, occupied, noccupied, occurrences, cooccurring
export traits, sitestats, sitestatnames, traitnames, commatrix
export addtraits!, addsitestats!
export asquantiles, asquantiles!
export coordstype, subset
export xcells, ycells, cells, xmin, xmax, ymin, ymax, xrange, yrange, xcellsize, ycellsize, cellsize, boundingbox #it is possible that I will export none of these
export sitetotals, speciestotals, getspecies, getsite
export groupspecies, groupsites
export aggregate
export @with, @traits, @sitestats
export matrixrandomizer, matrixrandomizations
export dispersionfield
export pairwise

include("DataTypes.jl")
include("Constructor_helperfunctions.jl")
include("Constructors.jl")
include("Sparse_matrixfunctions.jl")
include("Conveniencefunctions.jl")
include("ComMatrix.jl")
include("Subsetting.jl")
include("GetandSetdata.jl")
include("Gridfunctions.jl")
include("Grouping.jl")
include("Operations.jl")
include("Randomizations.jl")
include("PlotRecipes.jl")

end # module


text:
# Bi-objective shortest path problem (bisp)
#
# Exercise 9.5 page 269 of
# Multicriteria Optimization (2nd edt), M. Ehrgott, Springer 2005.

# ---- Packages to use
using vOptGeneric, JuMP, GLPK

# ---- Values of the instance to solve
M  = 50
C1 =  [ M 4 5 M M M ;  # coefficients's vector for arc (i,j) of the objective 1
        M M 2 1 2 7 ;
        M M M 5 2 M ;
        M M 5 M M 3 ;
        M M M M M 4 ;
        M M M M M M  ]

C2 =  [ M 3 1 M M M ; # coefficients's vector for arc (i,j) of the objective 2
        M M 1 4 2 2 ;
        M M M 1 7 M ;
        M M 1 M M 2 ;
        M M M M M 2 ;
        M M M M M M  ]

n  = size(C2,1)       # number of nodes

# ---- setting the model
bisp = vModel( GLPK.Optimizer ) ; JuMP.set_silent( bisp )

@variable( bisp, x[1:n,1:n], Bin )

@addobjective( bisp, Min , sum(x[i,j]*C1[i,j] for i=1:n, j=1:n) )
@addobjective( bisp, Min , sum(x[i,j]*C2[i,j] for i=1:n, j=1:n) )

@constraint( bisp, node_s, sum(x[1,j] for j = 1:n) == 1 )
@constraint( bisp, node_t, sum(x[i,n] for i = 1:n) == 1 )
@constraint( bisp, cstr[i=2:n-1], sum(x[i,j] for j = 1:n) - sum(x[j,i] for j =1:n) == 0 )

# ---- Invoking the solver (epsilon constraint method)
vSolve( bisp, method=:epsilon, step=0.5, verbose=true )

# ---- Querying the results
Y_N = getY_N( bisp )

# ---- Displaying the results
for s = 1:length(Y_N)
    X = value.(x, s)
    print("Path: ")
    for ind in findall(val -> val ≈ 1, X)
        i,j = ind.I
        print(" $i->$j ")
    end
    println("| Z = ",Y_N[s])
end

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric


text:
using ERA5Reanalysis
using Documenter

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

makedocs(;
    modules  = [ERA5Reanalysis,GeoRegions],
    doctest  = false,
    format   = Documenter.HTML(;
        prettyurls=get(ENV,"CI","false") == "true",
        canonical="https://natgeo-wong.github.io/ERA5Reanalysis.jl",
        assets=String[],
    ),
    authors  = "Nathanael Wong <natgeo.wong@outlook.com>",
    sitename = "ERA5Reanalysis.jl",
    pages    = [
        "Home" => "index.md",
        "ERA5Modules" => [
            "CDS Datasets"   => "modules/index.md",
            "Usage Examples" => "modules/examples.md",
        ],
        "ERA5Variables" => [
            "Variable Types"                    => "variables/index.md",
            "Retrieving Existing ERA5Variables" => "variables/read.md",
            "List of Existing ERA5Variables"    => "variables/list.md",
            "Creating Custom ERA5Variables"     => [
                "Custom Pressure-Level Variables" => "variables/pressures.md",
                "Custom Single-Level Variables"   => "variables/singles.md",
                "Resetting ERA5Variables"         => "variables/reset.md",
            ],
        ],
        "ERA5Regions" => [
            "The Basics of an ERA5Region"       => "regions/index.md",
            "Is it in an ERA5Region"            => "regions/isin.md",
            "Extracting Data using ERA5Regions" => "regions/gridded.md",
        ],
        "Downloading ERA5 Data" => [
            "The CDSAPI"              => "downloads/cdsapi.md"
            "Downloading with Julia"  => "downloads/julia.md"
            "Downloading with Python" => "downloads/python.md"
        ]
    ],
)

deploydocs(
    repo = "github.com/natgeo-wong/ERA5Reanalysis.jl.git",
    devbranch = "main"
)


text:
function read_merafile(filename::String, datatype::Symbol;

                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,

                        printfields::Bool=false,
                        verbose::Bool=verbose_mode)

        return read_merafile(filename=filename,
                        datatype=datatype,
                        xrange=xrange,
                        yrange=yrange,
                        zrange=zrange,
                        center=center,
                        range_unit=range_unit,
                        printfields=printfields,
                        verbose=verbose)
end

function read_merafile(filename::String;
                        datatype::Symbol=:hydro,

                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,

                        printfields::Bool=false,
                        verbose::Bool=verbose_mode)

        return read_merafile(filename=filename,
                        datatype=datatype,
                        xrange=xrange,
                        yrange=yrange,
                        zrange=zrange,
                        center=center,
                        range_unit=range_unit,
                        printfields=printfields,
                        verbose=verbose)
end

# todo: select variables
function read_merafile(;
                        filename::String="",
                        datatype::Symbol=:hydro,

                        xrange::Array{<:Any,1}=[missing, missing],
                        yrange::Array{<:Any,1}=[missing, missing],
                        zrange::Array{<:Any,1}=[missing, missing],
                        center::Array{<:Any,1}=[0., 0., 0.],
                        range_unit::Symbol=:standard,

                        printfields::Bool=false,
                        verbose::Bool=verbose_mode)

    info = info_merafile(filename=filename,
                        datatype=datatype,
                        printfields=printfields,
                        verbose=false) # suppress screen print from info_merafile

    if datatype == :hydro
        dtype = HydroDataType()
        if verbose println() end
        printtime("Get hydro data: ", verbose)
    elseif datatype == :particles
        dtype = PartDataType()
        if verbose println() end
        printtime("Get particle data: ", verbose)
    elseif datatype == :clumps
        dtype = ClumpDataType()
        if verbose println() end
        printtime("Get clump data: ", verbose)
    end

    #------------------
    #default variable selection for gethydro, getparticles, getclumps, getgravity, ....
    vars=[:all]
    if datatype==:clumps
        # read clumps-data of the selected variables
        column_names, NColumns = getclumpvariables(info, vars, info.fnames)
    else
        lmax=info.levelmax
        # create variabe-list and vector-mask (nvarh_corr) for gethydrodata-function
        # print selected variables on screen
        nvarh_list, nvarh_i_list, nvarh_corr, read_cpu, used_descriptors = prepvariablelist(info, datatype, vars, lmax, verbose)
    end
    #------------------

    #------------------
    # convert given ranges and print overview on screen
    ranges = prepranges(info, range_unit, verbose, xrange, yrange, zrange, center)
    #------------------
    if datatype==:clumps
        if verbose
            println("Read $NColumns colums: ")
            println(column_names)
        end
    end

    # read only data for database
    file = h5open(filename, "r")

        group = file[string(datatype)]
        data = group["data"]
        merafile_version = read(group["merafile_version"])
        column_names = read(data["names"])

        if datatype == :hydro
            if in("level", column_names)
                pkey=[:level,:cx, :cy, :cz]
            else
                pkey=[:cx, :cy, :cz]
            end
            t = table([read(data[i]) for i in column_names]...,
            names=Symbol.(column_names), pkey=pkey, presorted = false)

        elseif datatype == :particles

            if info.levelmax != info.levelmin # if AMR
                if info.descriptor.pversion == 0
                    Nkeys = [:level, :x, :y, :z, :id]
                elseif info.descriptor.pversion > 0
                    Nkeys = [:level, :x, :y, :z, :id, :family, :tag]
                end
            else # if uniform grid
                if info.descriptor.pversion == 0
                    Nkeys = [:x, :y, :z, :id]
                elseif info.descriptor.pversion > 0
                    Nkeys = [:x, :y, :z, :id, :family, :tag]
                end
            end

            t = table([read(data[i]) for i in column_names]...,
            names=Symbol.(column_names), pkey=collect(Nkeys), presorted = false)
        elseif datatype == :clumps
            t = table([read(data[i]) for i in column_names]...,
            names=Symbol.(column_names), presorted = false)
        end

        # Fill datatype
        dtype.info = info
        dtype.data = t
        dtype.boxlen = read(group["boxlen"])
        dtype.ranges = read(group["ranges"])

        desc = read(group["used_descriptors"])
        if desc == 0
            dtype.used_descriptors = Dict()
        else
            dtype.used_descriptors = desc
        end
        dtype.scale = info.scale

        if datatype == :hydro # Fill specific HydroDataType fields
            dtype.lmin = read(group["lmin"])
            dtype.lmax = read(group["lmax"])
            dtype.selected_hydrovars = read(group["selected_hydrovars"])
            dtype.smallr = read(group["smallr"])
            dtype.smallc = read(group["smallc"])
        elseif datatype == :particles # Fill specific PartDataType fields
            dtype.lmin = read(group["lmin"])
            dtype.lmax = read(group["lmax"])
            dtype.selected_partvars = Symbol.(read(group["selected_partvars"]))
        elseif datatype == :clumps # Fill specific ClumpDataType fields
            dtype.selected_clumpvars = Symbol.(read(group["selected_clumpvars"]))
        end

        # filter selected data region
        dtype = subregion(dtype, :cuboid,
                         xrange=xrange,
                         yrange=yrange,
                         zrange=zrange,
                         center=center,
                         range_unit=range_unit,
                         verbose=false)

        printtablememory(dtype, verbose)

        """
        if verbose
            ctype = read(group["compression_type"])
            cn    = read(group["compression"])

            println("Directory: ", info.path )
            println("merafile_version: ", merafile_version)
            println("Simulation code: ", info.simcode)
            println("DataType: ", datatype)
            println("Data variables: ", column_names)
            #if info.descriptor.usehydro
            #     println("Descriptor: ", dtype.used_descriptors )
            #end
            #println("Total fields: ", count_fields)
            println("Compression type: ", ctype, " with: ", cn )
            mem = usedmemory(dtype, false)
            println("Uncompressed memory size: ", round(mem[1], digits=3)," ", mem[2])
            s = filesize(filename)
            svalue, sunit = humanize(Float64(s), 3, "memory")
            println("File size: ", svalue, " ", sunit)
        end
        """
    close(file)
    return dtype
end


text:
using ProgenyTestingTools
using QMSimFiles
using Test
using Random
using Statistics

@testset "parameters" begin
   # no repeatability
   var_a = 6.0
   var_g = 30.0
   var_u = var_a + var_g
   var_pe = 0.0
   var_e = 64.0

   # pissibility 1: rep = 0
   var_p = var_u + var_pe + var_e
   h2_poly = var_a/var_p
   h2_qtl = var_g/var_p
   rep = (var_u+var_pe)/var_p
   par = PTParameters(50, var_p, h2_poly, h2_qtl, rep, 1.0)
   @test isnothing( check_parameters(par) )
   @test get_var_poly(par) ≈ var_a
   @test get_var_qtl(par) ≈ var_g
   @test get_var_pe(par) ≈ var_pe
   @test get_var_error(par) ≈ var_e

   # pissibility 2: rep = heritability
   rep = h2_poly + h2_qtl
   par = PTParameters(50, var_p, h2_poly, h2_qtl, rep, 1.0)
   @test isnothing( check_parameters(par) )
   @test get_var_poly(par) ≈ var_a
   @test get_var_qtl(par) ≈ var_g
   @test get_var_pe(par) ≈ var_pe
   @test get_var_error(par) ≈ var_e

   # repeated records
   var_a = 6.0
   var_g = 30.0
   var_u = var_a + var_g
   var_pe = 9.0
   var_e = 55.0
   var_p = var_u + var_pe + var_e
   h2_poly = var_a/var_p
   h2_qtl = var_g/var_p
   rep = (var_u+var_pe)/var_p
   par = PTParameters(50, var_p, h2_poly, h2_qtl, rep, 1.0)
   par = PTParameters(50, var_p, h2_poly, h2_qtl, rep, 1.0)
   @test isnothing( check_parameters(par) )
   @test get_var_poly(par) ≈ var_a
   @test get_var_qtl(par) ≈ var_g
   @test get_var_pe(par) ≈ var_pe
   @test get_var_error(par) ≈ var_e

   # out of range
   @test_throws ErrorException check_parameters(PTParameters(50, -100, 0.1, 0.1, 0.5, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, 100, -0.1, 0.1, 0.5, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, 100, 1.1, 0.1, 0.5, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, -100, 0.1, -0.1, 0.5, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, -100, 0.1, 1.1, 0.5, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, -100, 0.1, 0.1, -0.1, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, -100, 0.1, 0.1, 0.1, 1.0))
   @test_throws ErrorException check_parameters(PTParameters(50, -100, 0.1, 0.1, 1.1, 1.0))
end

@testset "generate_population" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   @test hp.maxAnimal == 10
   @test sum(hp.df[!,:male] .== true) == 5
   @test sum(hp.df[!,:male] .== false) == 5
end

@testset "migrate_from_hp!" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   idlist = [3,8]
   newidlist = migrate_from_hp!(hp,pop,idlist)
   @test pop.df[1,:male]
   @test !pop.df[2,:male]
   @test hp.df[3,:tbv] ≈ pop.df[1,:tbv]
   @test hp.df[8,:tbv] ≈ pop.df[2,:tbv]
   @test pop.df[!,:id] ≈ [1,2]
   @test pop.maxAnimal == 2

   # duplicated ID
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   idlist = [3,8,8,3,8]
   newidlist = migrate_from_hp!(hp,pop,idlist)
   @test pop.df[1,:male]
   @test !pop.df[2,:male]
   @test hp.df[3,:tbv] ≈ pop.df[1,:tbv]
   @test hp.df[8,:tbv] ≈ pop.df[2,:tbv]
   @test pop.df[!,:id] ≈ [1,2]
   @test pop.maxAnimal == 2
end

@testset "random_sampling" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=10,nf=10)
   cull!(hp,[10,20])

   # male only
   idlist = random_sampling(hp, 10, male=true, aliveonly=true, allowempty=true)
   @test all( sort(idlist) .== [1,2,3,4,5,6,7,8,9] )
   @test_throws ErrorException random_sampling(hp, 10, male=true, aliveonly=true, allowempty=false)
   idlist = random_sampling(hp, 10, male=true, aliveonly=false, allowempty=true)
   @test all( sort(idlist) .== [1,2,3,4,5,6,7,8,9,10] )
   ntests = 100
   ntested = 0
   for i in 1:ntests
      idlist = random_sampling(hp, 5, male=true, aliveonly=true, allowempty=true)
      if sum(map(x->in(x,[1,2,3,4,5,6,7,8,9]),idlist))==5
         ntested = ntested + 1
      end
   end
   @test ntests == ntested

   # female only
   idlist = random_sampling(hp, 10, female=true, aliveonly=true, allowempty=true)
   @test all( sort(idlist) .== [11,12,13,14,15,16,17,18,19] )
   @test_throws ErrorException random_sampling(hp, 10, female=true, aliveonly=true, allowempty=false)
   idlist = random_sampling(hp, 10, female=true, aliveonly=false, allowempty=true)
   @test all( sort(idlist) .== [11,12,13,14,15,16,17,18,19,20] )
   ntests = 100
   ntested = 0
   for i in 1:ntests
      idlist = random_sampling(hp, 5, female=true, aliveonly=true, allowempty=true)
      if sum(map(x->in(x,[11,12,13,14,15,16,17,18,19]),idlist))==5
         ntested = ntested + 1
      end
   end
   @test ntests == ntested

   # male and female
   idlist = random_sampling(hp, 20, male=true, female=true, aliveonly=true, allowempty=true)
   @test all( sort(idlist) .== [1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,17,18,19] )
   @test_throws ErrorException random_sampling(hp, 20, male=true, female=true, aliveonly=true, allowempty=false)
   idlist = random_sampling(hp, 20, male=true, female=true, aliveonly=false, allowempty=true)
   @test all( sort(idlist) .== [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20] )
   ntests = 100
   ntested = 0
   for i in 1:ntests
      idlist = random_sampling(hp, 15, male=true, female=true, aliveonly=true, allowempty=true)
      if sum(map(x->in(x,[1,2,3,4,5,6,7,8,9,11,12,13,14,15,16,17,18,19]),idlist))==15
         ntested = ntested + 1
      end
   end
   @test ntests == ntested
end

@testset "assign_year" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   # test 1
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   idlist = [3,5,8,9]
   newidlist = migrate_from_hp!(hp,pop,idlist)
   ProgenyTestingTools.assign_year!(pop,newidlist,[10,20])
   @test pop.df[1,:year] == 10
   @test pop.df[2,:year] == 20
   @test pop.df[3,:year] == 10
   @test pop.df[4,:year] == 20
   # test 2
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   idlist = [3,5,8,9]
   newidlist = migrate_from_hp!(hp,pop,idlist)
   ProgenyTestingTools.assign_year!(pop,newidlist,10)
   @test pop.df[1,:year] == 10
   @test pop.df[2,:year] == 10
   @test pop.df[3,:year] == 10
   @test pop.df[4,:year] == 10
end

@testset "generate_group (from population)" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [2,4,5]
   hp_females = [7,8,10]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   group = generate_group(pop, sires=pop_bulls[1:2], dams=pop_dams[1:2])
   @test group.n == 4 && group.maxSire==2 && group.maxDam == 2
   @test all(group.sires .== [1,2]) && all(group.dams .== [4,5])
end

@testset "copy_group" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [2,4,5]
   hp_females = [7,8,10]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   group = generate_group(pop, sires=pop_bulls[1:2], dams=pop_dams[1:2])
   group.generation .= [1,2,3,4]
   c = copy_group(group)
   @test c.n==group.n && c.maxSire==group.maxSire && c.maxDam==group.maxDam && c.groupid==group.groupid+1 &&
        all(c.id .== group.id) && all(c.generation .== group.generation) &&
        all(c.sires .== group.sires) && all(c.sires .== group.sires)
end

@testset "selectid from population" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [3,2,4,5]
   hp_females = [8,7,10,9]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   pop.df[2,:alive] = false
   pop.df[7,:alive] = false
   pop.df[!,:tbv] = [2,4,3,1, 2,1,3,4] * 1.0
   pop.df[!,:gebv] = [2,4,3,missing, 2,1,missing,4] * 1.0
   # simple selection: default aliveonly=true
   id = selectid(:male => x -> x==true, pop)
   @test all( id .== [1,3,4] )
   id = selectid(:male => x -> x==true, pop, aliveonly=false)
   @test all( id .== [1,2,3,4] )
   id = selectid(:male => x -> x==false, pop)
   @test all( id .== [5,6,8] )
   id = selectid(:male => x -> x==false, pop, aliveonly=false)
   @test all( id .== [5,6,7,8] )
   # multiple conditions
   id = selectid([:male,:tbv] => (x,y) -> x==true && y>=3.0, pop)
   @test all( id .== [3] )
   id = selectid([:male,:tbv] => (x,y) -> x==true && y>=3.0, pop, aliveonly=false)
   @test all( id .== [2,3] )
   # missing values
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=2.0, pop)
   @test all( id .== [1,3] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=2.0, pop, aliveonly=false)
   @test all( id .== [1,2,3] )
   # allowempty
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, pop, allowempty=true)
   @test all( id .== [] )
   @test_throws ErrorException selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, pop, allowempty=false)
   # particular id
   id = selectid(:male => x -> x==true, pop, idlist=[1,4,6,8])
   @test all( id .== [1,4] )
   id = selectid(:male => x -> x==true, pop, idlist=[1,1,6,4,1,8,4])
   @test all( id .== [1,4] )
   # sort
   id = selectid(:male => x -> x==true, pop, sortby=:tbv)
   @test all( id .== [3,1,4] )
   id = selectid(:male => x -> x==true, pop, sortby=:tbv, rev=false)
   @test all( id .== [4,1,3] )
   id = selectid(:male => x -> x==true, pop, sortby=:tbv, aliveonly=false)
   @test all( id .== [2,3,1,4] )
   id = selectid(:male => x -> x==true, pop, sortby=:gebv)
   @test all( id .== [4,3,1] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y), pop, sortby=:gebv)
   @test all( id .== [3,1] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, pop, sortby=:gebv, allowempty=true)
   @test all( id .== [] )
   # limit
   id = selectid(:male => x -> x==false, pop, limit=2)
   @test all( id .== [5,6] )
   id = selectid(:male => x -> x==false, pop, limit=4)
   @test all( id .== [5,6,8] )
   id = selectid(:male => x -> x==false, pop, limit=1, aliveonly=false)
   @test all( id .== [5] )
   id = selectid(:male => x -> x==false, pop, limit=4, aliveonly=false)
   @test all( id .== [5,6,7,8] )
end

@testset "selectid from group" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [3,2,4,5]
   hp_females = [8,7,10,9]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   pop.df[2,:alive] = false
   pop.df[7,:alive] = false
   pop.df[!,:tbv] = [2,4,3,1, 2,1,3,4] * 1.0
   pop.df[!,:gebv] = [2,4,3,missing, 2,1,missing,4] * 1.0
   group = generate_group(pop, sires=pop_bulls[1:3], dams=pop_dams[1:3], aliveonly=false)
   group.generation = collect(1:length(group.generation))*10
   # simple selection: default aliveonly=true
   id = selectid(:male => x -> x==true, group)
   @test all( id .== [1,3] )
   id = selectid(:male => x -> x==true, group, aliveonly=false)
   @test all( id .== [1,2,3] )
   id = selectid(:male => x -> x==false, group)
   @test all( id .== [5,6] )
   id = selectid(:male => x -> x==false, group, aliveonly=false)
   @test all( id .== [5,6,7] )
   # multiple conditions
   id = selectid([:male,:tbv] => (x,y) -> x==true && y>=3.0, group)
   @test all( id .== [3] )
   id = selectid([:male,:tbv] => (x,y) -> x==true && y>=3.0, group, aliveonly=false)
   @test all( id .== [2,3] )
   # missing values
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=2.0, group)
   @test all( id .== [1,3] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=2.0, group, aliveonly=false)
   @test all( id .== [1,2,3] )
   # allowempty
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, group, allowempty=true)
   @test all( id .== [] )
   @test_throws ErrorException selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, group, allowempty=false)
   # sort
   id = selectid(:male => x -> x==true, group, sortby=:tbv)
   @test all( id .== [3,1] )
   id = selectid(:male => x -> x==true, group, sortby=:tbv, rev=false)
   @test all( id .== [1,3] )
   id = selectid(:male => x -> x==true, group, sortby=:tbv, aliveonly=false)
   @test all( id .== [2,3,1] )
   id = selectid(:male => x -> x==true, group, sortby=:gebv)
   @test all( id .== [3,1] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y), group, sortby=:gebv)
   @test all( id .== [3,1] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, group, sortby=:gebv, allowempty=true)
   @test all( id .== [] )
   # limit
   id = selectid(:male => x -> x==false, group, limit=1, aliveonly=false)
   @test all( id .== [5] )
   id = selectid(:male => x -> x==false, group, limit=2, aliveonly=false)
   @test all( id .== [5,6] )
   id = selectid(:male => x -> x==true, group, limit=2, sortby=:tbv, aliveonly=false)
   @test all( id .== [2,3] )
   id = selectid(:male => x -> x==true, group, limit=1, sortby=:tbv, aliveonly=false)
   @test all( id .== [2] )
   id = selectid(:male => x -> x==true, group, limit=4, sortby=:tbv, aliveonly=false)
   @test all( id .== [2,3,1] )

   # removed dead individuals
   group = generate_group(pop, sires=pop_bulls[1:3], dams=pop_dams[1:3], aliveonly=true)
   # simple selection: default aliveonly=true
   id = selectid(:male => x -> x==true, group)
   @test all( id .== [1,3] )
   id = selectid(:male => x -> x==true, group, aliveonly=false)
   @test all( id .== [1,3] )
   id = selectid(:male => x -> x==false, group)
   @test all( id .== [5,6] )
   id = selectid(:male => x -> x==false, group, aliveonly=false)
   @test all( id .== [5,6] )
   # multiple conditions
   id = selectid([:male,:tbv] => (x,y) -> x==true && y>=3.0, group)
   @test all( id .== [3] )
   id = selectid([:male,:tbv] => (x,y) -> x==true && y>=3.0, group, aliveonly=false)
   @test all( id .== [3] )
   # missing values
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=2.0, group)
   @test all( id .== [1,3] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=2.0, group, aliveonly=false)
   @test all( id .== [1,3] )
   # allowempty
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, group, allowempty=true)
   @test all( id .== [] )
   @test_throws ErrorException selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, group, allowempty=false)
   # sort
   id = selectid(:male => x -> x==true, group, sortby=:tbv)
   @test all( id .== [3,1] )
   id = selectid(:male => x -> x==true, group, sortby=:tbv, rev=false)
   @test all( id .== [1,3] )
   id = selectid(:male => x -> x==true, group, sortby=:tbv, aliveonly=false)
   @test all( id .== [3,1] )
   id = selectid(:male => x -> x==true, group, sortby=:gebv)
   @test all( id .== [3,1] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y), group, sortby=:gebv)
   @test all( id .== [3,1] )
   id = selectid([:male,:gebv] => (x,y) -> x==true && !ismissing(y) && y>=5.0, group, sortby=:gebv, allowempty=true)
   @test all( id .== [] )
end

@testset "selectid from multiple groups" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [3,2,4,5]
   hp_females = [8,7,10,9]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   pop.df[2,:alive] = false
   pop.df[7,:alive] = false
   pop.df[!,:tbv] = [2,4,3,1, 2,1,3,4] * 1.0
   pop.df[!,:gebv] = [2,4,3,missing, 2,1,missing,4] * 1.0
   group1 = generate_group(pop, sires=pop_bulls[1:2], dams=pop_dams[1:2], aliveonly=false)
   group1.generation = collect(1:length(group1.generation))*10
   group2 = generate_group(pop, sires=pop_bulls[3:4], dams=pop_dams[3:4], aliveonly=false)
   group2.generation = collect(1:length(group2.generation))*10
   @test all( selectid([:male]=>x->x==true,[group1,group2]) .== [1,3,4])
   @test all( selectid([:male]=>x->x==true,[group1,group2],sortby=:tbv) .== [3,1,4])
end

@testset "add_sires! and add_dams!" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [3,2,4,5]
   hp_females = [8,7,10,9]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   pop.df[2,:alive] = false
   pop.df[7,:alive] = false
   group = generate_group(pop, sires=pop_bulls[3:3], dams=pop_dams[2:2])

   # capacity full
   @test 0 == add_sires!(group,[1,2])
   # extended
   group.maxSire = 2
   @test 1 == add_sires!(group,[1,2])
   group.maxSire = 3
   @test 0 == add_sires!(group,[1,3,1,2,3])
   @test 0 == add_sires!(group,[5,6,7,8])
   @test all( sort(group.sires) .==  [1,3])
   @test vacancy_for_sires(group) == 1

   # capacity full
   @test 0 == add_dams!(group,[5,6,7])
   # extended
   group.maxDam = 2
   @test 1 == add_dams!(group,[5,6,7])
   group.maxDam = 3
   @test 0 == add_dams!(group,[5,6,7])
   group.maxDam = 4
   @test 0 == add_dams!(group,[5,6,5,7,6,7])
   @test 0 == add_dams!(group,[1,2,3,4])
   @test all( sort(group.dams) .==  [5,6])
   @test vacancy_for_dams(group) == 2
end

@testset "candidates and status" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [3,2,4,5]
   hp_females = [8,7,10,9]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   pop.df[2,:alive] = false
   pop.df[7,:alive] = false
   mark_candidate!(pop,[1,3])
   @test all( pop.df[[1,3],:candidate] .== true)
   @test_throws ErrorException mark_candidate!(pop,[2])
   @test isnothing( mark_candidate!(pop,[2], check=false) )

   change_status!(pop,[5,6,8],1)
   @test all( pop.df[[5,6,8],:status] .== 1)
   @test_throws ErrorException change_status!(pop,[9],2)
   @test isnothing( change_status!(pop,[7],2) )
end

@testset "cull!" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=5)
   pop = generate_population(par)
   hp_males = [3,2,4,5]
   hp_females = [8,7,10,9]
   pop_bulls = migrate_from_hp!(hp,pop,hp_males)
   pop_dams = migrate_from_hp!(hp,pop,hp_females)
   cull!(pop,[2,7])
   @test !pop.df[2,:alive] && !pop.df[7,:alive]
   group = generate_group(pop, sires=pop_bulls[3:3], dams=pop_dams[2:2])

   group.maxSire = 3
   @test 2 == add_sires!(group,[1,2,3,4])
   @test all( sort(group.sires) .==  [1,3,4])
   @test 1 == cull!(group,[1])
   @test all( sort(group.sires) .==  [3,4])
   @test vacancy_for_sires(group) == 1
   @test 1 == cull!(group,[2,3])
   @test all( sort(group.sires) .==  [4])
   @test vacancy_for_sires(group) == 2
   @test 1 == cull!(group,[1,2,3,4,3,2,1])
   @test vacancy_for_sires(group) == 3

   group.maxDam = 4
   @test 2 == add_dams!(group,[5,6,7,8])
   @test all( sort(group.dams) .==  [5,6,8])
   @test 0 == cull!(group,[1])
   @test 0 == cull!(group,[1,2])
   @test 2 == cull!(group,[1,2,5,6])
   @test all( sort(group.dams) .==  [8])
   @test vacancy_for_dams(group) == 3
end

@testset "pedigree list and inbreeding" begin
   sires = [0,0,1,1,3,1,5]
   dams = [0,0,0,2,2,4,6]
   pedlist = ProgenyTestingTools.get_pedigree_list(sires,dams)
   f = ProgenyTestingTools.kernel_meuwissen_and_luo!(Float64, pedlist)
   ref_f = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.25, 0.15625]
   @test f ≈ ref_f
   g = ProgenyTestingTools.get_inbreeding(sires,dams)
   @test g ≈ ref_f

   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=3,nf=4)
   hp.df[:,:male] = [true,false,true,false,true,false,true]
   hp.df[:,:sire] = sires
   hp.df[:,:dam] = dams
   update_inbreeding!(hp)
   @test hp.df[:,:inb] ≈ ref_f

   f = copy(ref_f)
   f[end] = 0.0
   ProgenyTestingTools.get_inbreeding!(sires,dams,f,first=length(f))
   @test f ≈ ref_f
end

@testset "mating between groups! (dairy_standard_ai)" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=10,nf=10)
   pop = generate_population(par)
   bulls1 = migrate_from_hp!(hp,pop,random_sampling(hp,4,male=true))
   bulls2 = migrate_from_hp!(hp,pop,random_sampling(hp,4,male=true))
   cows = migrate_from_hp!(hp,pop,random_sampling(hp,8,female=true))
   mgroup1 = generate_group(pop, sires=bulls1)
   mgroup2 = generate_group(pop, sires=bulls2)
   fgroup = generate_group(pop, dams=cows)
   for i=1:4
      mating!(mgroup1, fgroup, "fgroup", n=1, method="dairy_standard_ai", plan="once_per_female", calving=false)
   end
   mating!(mgroup2, fgroup, "fgroup", n="all", method="dairy_standard_ai", plan="once_per_female", calving=false)

   # find progeny from sires in mgroup
   # extract sire group 1 by x->x[1]==mgroup1.gid
   @test all(fgroup.generation[9:16] .== 1)
   prog_id = selectid([:siregroup] => x->x==mgroup1.groupid, fgroup)
   @test all( (17 .<= prog_id) .& (prog_id .<= 20) )
   prog_id = selectid([:siregroup] => x->x==mgroup2.groupid, fgroup)
   @test all( (21 .<= prog_id) .& (prog_id .<= 24) )

   @test length(selectid([:pregnant] => x->x==true, pop))==8
   calving!(fgroup)
   @test length(selectid([:pregnant] => x->x==true, pop))==0

   # once more
   mgroup3 = copy_group(mgroup1)
   mgroup4 = copy_group(mgroup2)
   fgroup.dams .= sort(shuffle(17:24)[1:8])
   for i=1:4
      mating!(mgroup3, fgroup, "fgroup", n=1, method="dairy_standard_ai", plan="once_per_female", calving=false)
   end
   mating!(mgroup4, fgroup, "fgroup", n="all", method="dairy_standard_ai", plan="once_per_female", calving=false)

   @test all(fgroup.generation[17:24] .== 2)
   prog_id = selectid([:siregroup] => x->x==mgroup3.groupid, fgroup)
   @test all( (25 .<= prog_id) .& (prog_id .<= 28) )
   prog_id = selectid([:siregroup] => x->x==mgroup4.groupid, fgroup)
   @test all( (29 .<= prog_id) .& (prog_id .<= 32) )

   @test length(selectid([:pregnant] => x->x==true, pop))==8
   calving!(fgroup)
   @test length(selectid([:pregnant] => x->x==true, pop))==0
end

@testset "mating between groups! (unqual use of sires)" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=5,nf=1000)
   pop = generate_population(par)
   bulls = migrate_from_hp!(hp,pop,random_sampling(hp,5,male=true))
   cows = migrate_from_hp!(hp,pop,random_sampling(hp,1000,female=true))
   mgroup = generate_group(pop, sires=bulls)
   fgroup = generate_group(pop, dams=cows)
   prob = [0.05,0.05,0.05,0.05,0.8]
   mating!(mgroup, fgroup, "fgroup", n="all", method="dairy_standard_ai", plan="once_per_female", calving=false, prob=prob)
   n1 = sum(pop.df[:,:sire] .== 1)
   n2 = sum(pop.df[:,:sire] .== 2)
   n3 = sum(pop.df[:,:sire] .== 3)
   n4 = sum(pop.df[:,:sire] .== 4)
   n5 = sum(pop.df[:,:sire] .== 5)
   @show prob,n1,n2,n3,n4,n5
   @test n5>n1*10 && n5>n2*10 && n5>n3*10 && n5>n4*10
end

@testset "assign_phenotype!" begin
   par = PTParameters(50, 100, 0.5, 0, 0, 1.0)
   hp = generate_population(par,nm=20,nf=20)
   pop = generate_population(par)
   bulls = migrate_from_hp!(hp,pop,random_sampling(hp,4,male=true))
   cows = migrate_from_hp!(hp,pop,random_sampling(hp,4,female=true))
   mgroup = generate_group(pop, sires=bulls)
   fgroup = generate_group(pop, dams=cows)
   mating!(mgroup, fgroup, "fgroup", n="all", method="dairy_standard_ai", plan="once_per_female", calving=false, pmp=0.0)
   calving!(fgroup)
   assign_phenotype!(fgroup, gen=nothing, repeated=false)
   @test sum( map(x->length(x.cg), pop.animal) ) == sum(.!pop.df[fgroup.id,:male] )

   assign_phenotype!(fgroup, gen=nothing, repeated=false)
   assign_phenotype!(fgroup, gen=nothing, repeated=false)
   @test maximum(pop.df[:,:nrec]) == 1

   bulls2 = migrate_from_hp!(hp,pop,random_sampling(hp,4,male=true))
   mgroup2 = generate_group(pop, sires=bulls)
   fgroup.dams = selectid([:male]=>x->x!=true, fgroup, idlist=fgroup.id[fgroup.generation.==1])
   mating!(mgroup2, fgroup, "fgroup", n="all", method="dairy_standard_ai", plan="once_per_female", calving=false)
   calving!(fgroup)
   assign_phenotype!(fgroup, gen=2, repeated=false)
   @test length(selectid([:siregroup,:male]=>(x,y)->x==mgroup2.groupid && y==false,fgroup)) == sum(.!pop.df[fgroup.id[fgroup.generation .== 2],:male] )
   selected_id = selectid([:siregroup,:male]=>(x,y)->x==mgroup2.groupid && y==false,fgroup)
   if length(selected_id)>0
      @test maximum(pop.df[selected_id,:nrec])==1
   else
      @test_skip "invalid simulated data"
   end

   popx = generate_population(par)
   bulls3 = migrate_from_hp!(hp,popx,random_sampling(hp,4,male=true))
   cows3 = migrate_from_hp!(hp,popx,random_sampling(hp,4,female=true))
   mgroup3 = generate_group(popx, sires=bulls3)
   fgroup3 = generate_group(popx, dams=cows3)
   mating!(mgroup3, fgroup3, "fgroup", n="all", method="dairy_standard_ai", plan="once_per_female", calving=false, pmp=0.0)
   calving!(fgroup3)
   assign_phenotype!(fgroup3, gen=nothing, repeated=false, idlist=[5,8,11])
   @test sum(popx.df[:,:nrec] .> 0)==3
   @test all( popx.df[findall(popx.df[:,:nrec] .> 0),:id] .== [5,8,11] )
end

@testset "genomic data" begin
   par = PTParameters(50, 1.0, 0.2, 0.1, 0, 0.2)
   qmsimfile="base_001.h5"
   hp = generate_population(par, qmsimfile, nm=10, nf=10, gfile=qmsimfile)
   pop = generate_population(par, qmsimfile, gfile=tempname())
   base_males = migrate_from_hp!(hp,pop,random_sampling(hp,4,male=true),year=[1,2,3,4,5])
   base_females = migrate_from_hp!(hp,pop,random_sampling(hp,2,female=true),year=[6])
   bull_group = generate_group(pop, sires=base_males)
   cow_group = generate_group(pop, dams=base_females)
   progeny_id = mating!(bull_group, cow_group, "fgroup", year=10, method="dairy_standard_ai", n=2)
   a1 = read_qmsim_individual_hdf5(pop.map, pop.gfile, 1)
   a2 = read_qmsim_individual_hdf5(pop.map, pop.gfile, 2)
   a3 = read_qmsim_individual_hdf5(pop.map, pop.gfile, 3)
   a4 = read_qmsim_individual_hdf5(pop.map, pop.gfile, 4)
   a5 = read_qmsim_individual_hdf5(pop.map, pop.gfile, 5)
   a6 = read_qmsim_individual_hdf5(pop.map, pop.gfile, 6)
   @test pop.df[1,:qbv] ≈ a1.tbv
   @test pop.df[2,:qbv] ≈ a2.tbv
   @test pop.df[3,:qbv] ≈ a3.tbv
   @test pop.df[4,:qbv] ≈ a4.tbv
   @test pop.df[5,:qbv] ≈ a5.tbv
   @test pop.df[6,:qbv] ≈ a6.tbv
end

@testset "polygenic Mendelian sampling" begin
    ntest = 100000
    samples = zeros(ntest)
    inbs = 0.0
    inbd = 0.0
    var_a = 4.0
    #k = 0.0
    for i=1:ntest
       samples[i] = ProgenyTestingTools.get_ms_deviation(inbs, inbd, var_a)
    end
    @test isapprox(mean(samples),0.0, atol=0.1)
    @test isapprox(std(samples),sqrt(0.5*var_a), atol=0.1)

    inbs = 0.125
    inbd = 0.250
    var_a = 4.0
    #k = 0.0
    for i=1:ntest
       samples[i] = ProgenyTestingTools.get_ms_deviation(inbs, inbd, var_a)
    end
    @test isapprox(mean(samples),0.0, atol=0.1)
    @test isapprox(std(samples),sqrt((0.5 - 0.25*(inbs+inbd))*var_a), atol=0.1)

    inbs = 0.125
    inbd = 0.250
    var_a = 4.0
    k = 0.36
    for i=1:ntest
       samples[i] = ProgenyTestingTools.get_ms_deviation(inbs, inbd, var_a, k)
    end
    @test isapprox(mean(samples),0.0, atol=0.1)
    @test isapprox(std(samples),sqrt((1-k)*(0.5 - 0.25*(inbs+inbd))*var_a), atol=0.1)
end


text:
module Estimator

using ..Concepts

import Distributions
import StatsBase

export MLE,
  ProfileLikelihood,
  MOM

abstract type EstimationProcedure{T<:Any} end

## not super important
struct MLE{T<:Any} <:Concepts.AbstractEstimator
    of::T
    function MLE{T}() where T<:Any
        return new{T}()
    end

    function MLE{T}(of::T) where T<:Any
        return new{T}(of)
    end
end

## not super important
struct MOM{T<:Any} <:Concepts.AbstractEstimator
    of::T
    function MOM{T}() where T<:Any
        return new{T}()
    end

    function MOM{T}(of::T) where T<:Any
        return new{T}(of)
    end
end

struct ProfileLikelihood <: EstimationProcedure{MLE} end

function Concepts.estimator(name::MOM{AbstractNegativeBinomial}, data::AutoboxedArray{T}) where T<:Real
  EX = StatsBase.mean(data)
  VarX = StatsBase.var(data)
  pₘₒₘ = EX / VarX
  rₘₒₘ = EX * (pₘₒₘ) / (1 - pₘₒₘ)
  return Dict{Symbol, Float64}(:p => pₘₒₘ, :r => rₘₒₘ)
end

# const MLE(of::Union{T,Symbol}) where T<:Any =
#     !isa(of,Symbol) ? MLE{T}(of) : MLE(convert(of))

const MLE(of::Union{T,Symbol}) where T<:Any =
    !isa(of,Symbol) ? MLE{T}(of) : MLE(type_conversion(of))

@overload
function Concepts.estimator(name::MLE{AbstractGaussian}, data::AutoboxedArray{T};
                            method::EstimationProcedure{MLE} = ProfileLikelihood()) where T<:Real

    if method == ProfileLikelihood()
        # for now, we use Distributions.jl. Will be replaced in the future.
        est = Distributions.fit_mle(Distributions.Gaussian,data)
        return Dict(:μ => est.μ, :σ => est.σ)
    end
end

@overload
function Concepts.estimator(name::MLE{AbstractGamma},data::AutoboxedArray{T};
                            method::EstimationProcedure{MLE} = ProfileLikelihood()) where T<:Real
    # default uses profile likelihood method
    if method == ProfileLikelihood()
        # for now, we use Distributions.jl. Will be replaced in the future.
        est = Distributions.fit_mle(Distributions.Gamma,data)
        return Dict(:α => est.α, :θ => est.θ)
    end
end

# data layout needs to be implemented properly

@overload
function Concepts.estimator(name::MLE{AbstractPoisson},data::AutoboxedArray{T};
                            method::EstimationProcedure{MLE} = ProfileLikelihood(),
                            data_layout::Symbol = :flatten) where T<:Real
    # default uses profile likelihood method
    if method == ProfileLikelihood()
           λ = sum(data)/length(data)
        return Dict(:λ => λ)
    end
end

@overload
function Concepts.estimator(name::MLE{AbstractBernoulli},data::AutoboxedArray{T};
                            method::EstimationProcedure{MLE} = ProfileLikelihood(),
                            data_layout::Symbol = :flatten) where T<:Real
    # default uses profile likelihood method
    if method == ProfileLikelihood()
           est = Distributions.fit_mle(Distributions.Bernoulli,data)
        return Dict(:p => est.p)
    end
end

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(BinaryProvider.safe_isfile), String})
    precompile(Tuple{typeof(BinaryProvider.info_onchange), String, String, Int64})
    precompile(Tuple{typeof(BinaryProvider.libdir), BinaryProvider.Prefix})
    precompile(Tuple{typeof(BinaryProvider.readuntil_many), Base.Pipe, Array{Char, 1}})
    precompile(Tuple{typeof(BinaryProvider.locate), BinaryProvider.ExecutableProduct})
    precompile(Tuple{typeof(BinaryProvider.package), BinaryProvider.Prefix, String, Base.VersionNumber})
    precompile(Tuple{typeof(BinaryProvider.libdir), BinaryProvider.Prefix, BinaryProvider.Linux})
    precompile(Tuple{typeof(BinaryProvider.libdir), BinaryProvider.Prefix, BinaryProvider.MacOS})
    precompile(Tuple{typeof(BinaryProvider.locate), BinaryProvider.LibraryProduct})
    precompile(Tuple{typeof(BinaryProvider.collect_stdout), BinaryProvider.OutputCollector})
    precompile(Tuple{typeof(BinaryProvider.collect_stderr), BinaryProvider.OutputCollector})
end


text:
module AlgebraTests

using Test

@testset "AlgebraInterfaces" begin include("AlgebraInterfacesTests.jl") end

@testset "SparseMatrixCSC" begin include("SparseMatrixCSCTests.jl") end

@testset "SparseMatrixCSR" begin include("SparseMatrixCSRTests.jl") end

@testset "SymSparseMatrixCSR" begin include("SymSparseMatrixCSRTests.jl") end

@testset "NonlinearOperators" begin include("NonlinearOperatorsTests.jl") end

@testset "NonlinearSolvers" begin include("NonlinearSolversTests.jl") end

@testset "NLSolvers" begin include("NLSolversTests.jl") end

@testset "LinearSolvers" begin include("LinearSolversTests.jl") end

end # module


text: ack2(m::Integer, n::Integer) = m == 0 ? n + 1 : n == 0 ? ack2(m - 1, 1) : ack2(m - 1, ack2(m, n - 1))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@with_kw mutable struct Clouds
    Xp::Vector{Float64}
    Yp::Vector{Float64}
    Vp::Vector{Float64}
    Fp::Vector{Float64}
    Sp::Vector{Float64}
    Ip::Vector{Float64}
    logs::Dict{String, <:Any}

    belongs::Vector{Int64}=[]
    I::Vector{Int64}=[]
    k::Vector{Int64}=[]

    Xc::Vector{Float64}=[]
    Yc::Vector{Float64}=[]
    Vc::Vector{Float64}=[]
    Fc::Vector{Float64}=[]
    Sc::Vector{Float64}=[]
    Ic::Vector{Float64}=[]

    r::Vector{Float64}=[]
    R::Vector{Float64}=[]
    Φ::Vector{Float64}=[]
    VC::Vector{Float64}=[]
    dV::Vector{Float64}=[]
    Z::Vector{Float64}=[]
    inc::Vector{Float64}=[]
    pa::Vector{Float64}=[]
    dT::Vector{Float64}=[]
    Tnet::Vector{Float64}=[]
    ϕ::Vector{Float64}=[]
    tilt::Vector{Float64}=[]
end

"""
chunk3(arr::Vector{Float64})

Splits a vector in chunks of size 3
## Examples
```jldoctest
    julia> chunk3([1.,2.,3.,4.,5.,6.])
    2-element Vector{Vector{Float64}}:
    [1.0, 2.0, 3.0]
    [4.0, 5.0, 6.0]
```
```jldoctest
    julia> chunk3([1.,2.,3.,4.,5.])
    2-element Vector{Vector{Float64}}:
    [1.0, 2.0, 3.0]
    [4.0, 5.0]
```
"""
function chunk3(arr::Vector{Float64})::Vector{Vector{Float64}}
    return [arr[i:min(i + 2, end)] for i in 1:3:length(arr)]
end

"""
gaus(v,vc,A,sv) = ``A\\\\exp(\\\\frac{-(v-v_c)^2}{2 s_v^2.0})```

Gaussian function with scale ``A`` center at ``v_c`` and dispersion ``s_v``
## Examples
```jldoctest
    julia> gaus(80.0,100.0,0.1,37.0)
    0.08640781715118678
```
"""
function gaus(v::Float64,vc::Float64,A::Float64,sv::Float64)::Float64
    return A *exp(-0.5*(v-vc)^2.0 /sv^2.0)
end

"""
gmodel!(s::Vector{Float64},v::Vector{Float64},p::Vector{Float64})::Vector{Float64}

Multiple gaussians one dimension model with coordinates ``v``.

The parameters of the gaussians are infered from splitted ``p`` to N*[vc,A,sv] where ``N`` is the number of gaussians (infered also from p)

## Examples
```jldoctest
julia>  V = collect(-350.0:20:350)
        s=gmodel(V,[0.1,120.0,23.0,0.12,-60.0,33.0])
    36-element Vector{Float64}: ....
```
"""
function gmodel(v::Vector{Float64},p::Vector{Float64})::Vector{Float64}
    s=zeros(size(v))
    pars=chunk3(p)
    for (A,vc,sv) in pars
        si = [gaus(vi,vc,A,sv) for vi in v]
        s = s .+ si
    end
    return s
end

"""
gmodel_residual!(s::Vector{Float64},v::Vector{Float64},p::Vector{Float64})::Vector{Float64}

Sum of squared Residual of vector ``s`` after the subtraction of multiple gaussians one dimensional model with coordinates ``v``.

The parameters of the gaussians are infered from splitted ``p`` to N*[vc,A,sv] where ``N`` is the number of gaussians (infered also from p)

## Examples
```jldoctest
julia>  V = collect(-350.0:20:350)
        s=gmodel(V,[0.1,120.0,23.0,0.12,-60.0,33.0])
        gmodel_residual!(s,V,[0.1,120.0,23.0,0.12,-60.0,33.0])
        isapprox(gmodel_residual!(s,V,[0.1,120.0,23.0,0.12,-60.0,33.0]), 0.0; atol=eps(Float64), rtol=0)
        true
```
"""
function gmodel_residual!(s::Vector{Float64},v::Vector{Float64},p::Vector{Float64})::Float64
    pars=chunk3(p)
    for (A,vc,sv) in pars
        si = [gaus(vi,vc,A,sv) for vi in v]
        s = s .- si
    end
    return sum(s .^2.0)
end

"""
fit_pixel(signal::Vector{Float},rms::Float64,vv::Vector{Float64},er::Float64;vmin=-400.0,vmax=400.0,svmin=15.0,svmax=30.0,smax=1000.0,maxsteps=100000,popsize=80,debug=false,N0=1,Nmax=8)

Fit a multiple gaussian model in ``signal`` by minimizing the Bayesian Information Criterion (BIC):
``\\mathrm{BIC}=k\\log(n)-2\\log(\\hat{L})`` where ``k3\\times N`` is the number of the parametes, ``n`` the signal size and ``\\hat{L}`` the maximum likelihood estimation for a ``N``-gaussians model.
The algorithm fits the signal with multiple gaussians, calculating

## Input
- ``signal::Vector{Float}`` The input signal
- ``rms::Float64`` The noise of the data (it used in scale boundaries).
- ``er::Float64`` The error of the data as it used in the likelihood.

## Parameters
### Model Parameters boundaries
- ``vmin`` The lower boundary on central velocities (-400)
- ``vmax`` The upper boundary on central velocities (400)
- ``svmin`` The lower boundary on velocity dispersion (15)
- ``svmax`` The upper boundary on velocity dispersion (30)
- ``smax`` The upper boundary on the scaling in units of rms (1000)
- ``smax`` The upper boundary on the scaling in units of rms (1000)
### Fitter Parameters
- ``maxsteps`` The steps of the fitter for every model fitting (30000)
- ``popsize`` Population size for the differential evolution global fitter
- ``N0`` The starting number of gaussians (1)
- ``Nmax`` Maximum number of gaussians (8)
- ``debug`` Make plots of the procedure (false)

## Output
Best Parameters, Minimum BIC

## Examples
```jldoctest
julia>  V = collect(-350.0:20:350)
        s=gmodel(V,[0.1,120.0,23.0,0.12,-60.0,33.0])
        bpars,minbic=fit_pixel(s,0.02,V,0.02;vmin=-400.0,vmax=400.0,svmin=15.0,svmax=30.0,smax=1000.0,maxsteps=20000,popsize=80,N0=1,Nmax=8)
    ([0.1255719888895172, -60.00001027466198, 30.0, 0.0999767960308078, 119.99487202319484, 23.012779807985947], -193.52458049526706)
```
## Notes
The likelihood function is:
    ``2\\log\\hat{L} = -\\Big(\\sum_{i=1}^{n}\\frac{\\big(F_i-S(V,\\mathbf{V}_c,\\mathbf{F},\\mathbf{\\sigma})\\big)^2}{\\sigma_F^2} +n\\log(2\\pi\\sigma_F^2)\\Big)=-\\Big(\\chi^2+n\\log(2\\pi\\sigma_F^2)\\Big)``
"""
function fit_pixel(signal::Vector{Float64},rms::Float64,vv::Vector{Float64},er::Float64;vmin=-400.0,vmax=400.0,svmin=15.0,smin=0.75,svmax=30.0,smax=1000.0,maxsteps=30000,popsize=40,N0=1,Nmax=8)
    lowerc = [smin*rms,vmin,svmin]
    upperc = [smax*rms,vmax,svmax]

    lower=Float64[]
    upper=Float64[]

    minbic=1e10
    minchi=1000.0
    bic=minbic-11
    N=N0
    if N>1
        for ni in 1:N-1 #populate the boundaries
            lower=vcat(lower,lowerc)
            upper=vcat(upper,upperc)
        end
    end
    n=length(signal)
    bpars=Float64[]

    #chi_ = [minchi]
    L(p) = gmodel_residual!(signal,vv,p)
    while N<Nmax
        #@show N
        lower=vcat(lower,lowerc)
        upper=vcat(upper,upperc)
        k=length(lower)
        opt = bbsetup(L; SearchRange = collect(zip(lower, upper)),PopulationSize=popsize,MaxSteps = maxsteps, TraceMode = :silent)
        res1 = bboptimize(opt)
        pbb=best_candidate(res1)
        chi=best_fitness(res1)
        #@show pbb,chi
        logL=-0.5*(chi/er^2.0+n*log(2.0*pi*er^2.0))
       # @show logL
        bic=k*log(n)-2.0*logL
        #@show bic

        if (bic>minbic)
            break
        else
            minbic=bic
            bpars=pbb
            N+=1
        end
    end

    return bpars,minbic
end

"""
fit_pixels(ipixels::Vector{CartesianIndex{2}},data::Dict{String, Any};x0=-2.0,x1=2.0,y0=-5.0,y1=5.0,sigma=3.0,sigmac=3.0,crit=3,N0=1,vmin=-400.0,vmax=400.0,svmin=15.0,svmax=28.0,maxsteps=8000,popsize=22)

Fit a multiple gaussian model for every spaxel in ipixels ``signal`` by minimizing the Bayesian Information Criterion (BIC):
``\\mathrm{BIC}=k\\log(n)-2\\log(\\hat{L})`` where ``k3\\times N`` is the number of the parametes, ``n`` the signal size and ``\\hat{L}`` the maximum likelihood estimation for a ``N``-gaussians model.
The algorithm fits the signal with multiple gaussians, calculating

## Input
- ``Data`` Dictionary Data
## Parameters
### Data space boundaries and masking
- ``x0``, ``x1``, ``y0``, ``y1`` x,y boundaries in kpc
- ``sigma`` sigma mask on 0th moment

### Model Parameters boundaries
- ``vmin`` The lower boundary on central velocities (-400)
- ``vmax`` The upper boundary on central velocities (400)
- ``svmin`` The lower boundary on velocity dispersion (15)
- ``svmax`` The upper boundary on velocity dispersion (30)
- ``smax`` The upper boundary on the scaling in units of rms (1000)
- ``smax`` The upper boundary on the scaling in units of rms (1000)
### Fitter Parameters
- ``maxsteps`` The steps of the fitter for every model fitting in units of number of parameter=3N. So, for the first gaussian it will run ``maxsteps``, for the second ``3X maxsteps`` etc (30000)
- ``popsize`` Population size for the differential evolution global fitter
- ``N0`` The starting number of gaussians (1)
- ``Nmax`` Maximum number of gaussians (8)
- ``debug`` Make plots of the procedure (false)

## Output
Best Parameters, Minimum BIC

## Examples
```jldoctest
julia>  V = collect(-350.0:20:350)
        s=gmodel(V,[0.1,120.0,23.0,0.12,-60.0,33.0])
        bpars,minbic=fit_pixel(s,0.02,V,0.02;vmin=-400.0,vmax=400.0,svmin=15.0,svmax=30.0,smax=1000.0,maxsteps=20000,popsize=80,N0=1,Nmax=8)
    ([0.1255719888895172, -60.00001027466198, 30.0, 0.0999767960308078, 119.99487202319484, 23.012779807985947], -193.52458049526706)
```
"""
function fit_pixelst(data::Dict{String, Any};x0=-2.0,x1=2.0,y0=-5.0,y1=5.0, sigma=4.0,ncrit=1,N0=1,Nmax=8,vmin=-400.0,vmax=400.0, svmin=15.0,svmax=28.0, maxsteps=3000,popsize=22)
    beam_area=data["beam"][1]*data["beam"][2]*pi/(4.0*log(2.0)) #kpc^2
    pixel_area=data["dx"]^2.0
    mom0_data = sum(data["data"] .> sigma*data["rms"], dims = 3)[:, :, 1];
    ipixels=findall(mom0_data .> ncrit)
    @info "Starting line fitter on $(size(ipixels)[1]) spaxels"
    fit_results=[]
    Threads.@threads for n in 1:length(ipixels)
        ix=ipixels[n][1]
        iy=ipixels[n][2]
        x=data["X"][ix]
        y=data["Y"][iy]
        if ((x > x0) & (x < x1)) & ((y > y0) & (y < y1))
            signal=data["data"][ix,iy,:]
            pbb,_=fit_pixel(signal,sigma*data["rms"],data["V"],data["rms"];vmin=vmin,vmax=vmax,svmin=svmin,svmax=svmax,maxsteps=maxsteps,popsize=popsize,N0=N0,Nmax=Nmax)
            for gp in chunk3(pbb)
                A,vc,sv = gp
                Ip=A .* sv .* sqrt(2.0*pi) .*(pixel_area/beam_area) #jy km/s
                push!(fit_results,[x,y,vc,A,sv,Ip])
            end
        end
    end
    @info "Found $(length(fit_results)) subpixel clouds, dont forget to save them."
    cdp=Dict{String, Any}("svmax"=>svmax,"svmin"=>svmin,"x0"=>x0,"x1"=>x1,"y0"=>y0,"y1"=>y1,"vmin"=>vmin,"vmax"=>vmax,"sigma"=>sigma,"ncrit"=>ncrit,"maxsteps"=>maxsteps,"popsize"=>popsize,"N0"=>N0,"Nmax"=>Nmax)

	Xp=fill(NaN,length(fit_results))
	Yp=fill(NaN,length(fit_results))
	Vp=fill(NaN,length(fit_results))
	Fp=fill(NaN,length(fit_results))
	Sp=fill(NaN,length(fit_results))
	Ip=fill(NaN,length(fit_results))
	n=0
	for i in 1:length(fit_results)
		if isassigned(fit_results,i)
			n+=1
		    Xp[n]=fit_results[i][1]
			Yp[n]=fit_results[i][2]
			Vp[n]=fit_results[i][3]
			Fp[n]=fit_results[i][4]
			Sp[n]=fit_results[i][5]
			Ip[n]=fit_results[i][6]
		end
	end

    return Clouds(Xp=Xp,Yp=Yp,Vp=Vp,Fp=Fp,Sp=Sp,Ip=Ip,logs=cdp)
end

# function merge_pixel_clouds(X,Y,V,A;dx=0.1,dy=0.1,dv=30.0)
#     belongs = zeros(Int,length(X))
#     for i in reverse(sortperm(A)) #[1] #debug
#         #@show i #debug
#         if belongs[i] == 0 #if it hasnt merged to another source
#             belongs[i]=i #merge it to itself
#             Xi=X[i]
#             Yi=Y[i]
#             Vi=V[i]
#             Fi=A[i]
#             #@show Xi,Yi,Vi,Fi #debug
#             # find the distance to all other sources
#             dist= (Xi .-X) .^2.0 ./dx^2.0 .+(Yi .-Y) .^2.0 ./dy^2.0 .+(Vi .-V).^2.0 ./dv^2.0 #.<1.0

#             # sort the distance
#             local_sources_ind=sortperm(dist)#[1:50]
#             local_max_flux=Fi
#             #@show local_sources_ind
#             n=0
#             for j in local_sources_ind[2:end] #
#                 #@show n,j #debug
#                 d= dist[j]
#                 if belongs[j] == 0 #if it hasnt merged to another source
#                     n+=1
#                     #@show d #debug
#                     flux = A[j]
#                     #@show flux debug
#                     if flux <= local_max_flux
#                         belongs[j]=i#belongs[i]
#                     end
#                 end
#                 if  d>1
#                     break
#                 end
#             end
#         end
#     end
#     return belongs
# end

movingaverage(g, n) = [i < n ? mean(g[begin:i]) : mean(g[i-n+1:i]) for i in 1:length(g)]

# function merge_pixel_clouds(X::Vector{Float64},Y::Vector{Float64},V::Vector{Float64},F::Vector{Float64},sV::Vector{Float64};dr=0.1,dv=30.0,mp=4,df=5000.0,roll_d=0.1,max_d=2.0,plots=false)
function merge_pixel_clouds!(clouds::Clouds;dr=0.1,dv=30.0,mp=4,df=5000.0,roll_d=0.1,max_d=2.0)

    X=clouds.Xp
    Y=clouds.Yp
    V=clouds.Vp
    F=clouds.Fp
    Ip=clouds.Ip
    sV=clouds.Sp

    @info "Merging $(length(X)) subpixel clouds"
    belongs = zeros(Int,length(X))
    for (n,i) in enumerate(reverse(sortperm(F))) #debug
        #@show i #debug
        #@show n,i
        if belongs[i] == 0 #if it hasnt merged to another source
            belongs[i]=i #merge it to itself
            Xi=X[i]
            Yi=Y[i]
            Vi=V[i]
            Fi=F[i]
            dist= (Xi .-X) .^2.0 ./dr^2.0 .+(Yi .-Y) .^2.0 ./dr^2.0 .+(Vi .-V).^2.0 ./dv^2.0
            pd=sortperm(dist)
            d=dist[pd]
            Flux=F[pd]
            if sum(d .< max_d)>mp
                ma=movingaverage(Flux[d .< max_d],roll_d)
                # if n==plot
                #     fig,ax,sc=lines(d[d .< max_d],Flux[d .< max_d])
                #     lines!(ax,d[d .< max_d],ma)
                #     #vlines!(ax,1.0)
                #     #display(fig)
                # end
                f0=Fi
                #f=f0

                for (maj,j) in zip(ma,pd[d .< max_d])

                    if (belongs[j]==0)
                        f=maj
                        if (f<=f0+df)
                            belongs[j]=i
                            f0=f
                            # if n==plot
                            #     scatter!(ax,[dist[d .< max_d][j]],[F[d .< max_d][j]],marker='X')
                            #     display(fig)
                            #     sleep(0.05)
                            # end
                        else
                            break
                        end
                    end
                end
            end

        end
    end
    I_pre = unique(belongs)
    I=Vector{Int}()
    sI=Vector{Float64}()
    sIc=Vector{Float64}()
    for i in I_pre
        if sum(belongs .==i) > mp
            push!(I,i)
            push!(sI,sqrt(mean(sV[belongs .==i])^2.0 +std(V[belongs .==i] .- V[i])^2.0))
            push!(sIc,sum(Ip[belongs .==i]))
        end
    end
    @info "Found $(length(I)) unique clouds"

    setfield!(clouds, :belongs, belongs)
    setfield!(clouds, :I, I)
    setfield!(clouds, :Sc, sI)
    setfield!(clouds, :Xc, X[I])
    setfield!(clouds, :Yc, Y[I])
    setfield!(clouds, :Vc, V[I])
    setfield!(clouds, :Fc, F[I])
    setfield!(clouds, :Ic, sIc)
    setfield!(clouds, :logs, merge(clouds.logs,Dict{String, Any}("dr"=>dr,"dv"=>dv,"mp"=>mp,"df"=>df,"roll_d"=>roll_d,"max_d"=>max_d)))
end

function filter_clouds(clouds,rect::Vector{Dict{String, Float64}};neg=false,x0=-10.0,x1=10.0,y0=-10.0,z0=-10.0,z1=10.0,y1=10.0,v0=-1500.0,v1=1500.0,s0=0.0,s1=500.0,f0=0.0,f1=100.0,dv0=-1000.0,dv1=1000.0,R0=0.0,R1=100.0,r0=0.0,r1=100.0)

    mask=clouds.Xp .>1000.0
    maskc=clouds.Xc .>1000.0

	@info "Filtering $(sum(.!mask)) clouds"
    for rd in rect
		@info "Applying mask $(rd)"
        x0i = haskey(rd,"x0") ? rd["x0"] : x0
        x1i = haskey(rd,"x1") ? rd["x1"] : x1

        y0i = haskey(rd,"y0") ? rd["y0"] : y0
        y1i = haskey(rd,"y1") ? rd["y1"] : y1

        z0i = haskey(rd,"z0") ? rd["z0"] : z0
        z1i = haskey(rd,"z1") ? rd["z1"] : z1

		R0i = haskey(rd,"R0") ? rd["R0"] : R0
        R1i = haskey(rd,"R1") ? rd["R1"] : R1

		r0i = haskey(rd,"r0") ? rd["r0"] : r0
        r1i = haskey(rd,"r1") ? rd["r1"] : r1

        v0i = haskey(rd,"v0") ? rd["v0"] : v0
        v1i = haskey(rd,"v1") ? rd["v1"] : v1

		vc0i = haskey(rd,"vc0") ? rd["vc0"] : v0
        vc1i = haskey(rd,"vc1") ? rd["vc1"] : v1

        s0i = haskey(rd,"s0") ? rd["s0"] : s0
        s1i = haskey(rd,"s1") ? rd["s1"] : s1

        f0i = haskey(rd,"f0") ? rd["f0"] : f0
        f1i = haskey(rd,"f1") ? rd["f1"] : f1

		mask_rec= ((clouds.Xp .>x0i) .&& (clouds.Xp .<x1i)) .&& ((clouds.Fp .>f0i) .&& (clouds.Fp .<f1i)) .&& ((clouds.Yp .>y0i) .&& (clouds.Yp .<y1i)) .&& ((clouds.Vp .>v0i) .&& (clouds.Vp .<v1i)) .&& ((clouds.Sp .>s0i) .&& (clouds.Sp .<s1i))
		#@show sum(mask_rec)
		if length(clouds.r)>0
            dv0i = haskey(rd,"dv0") ? rd["dv0"] : dv0
			dv1i = haskey(rd,"dv1") ? rd["dv1"] : dv1
            mask_rec=mask_rec .&& ((clouds.dV .>dv0i) .&& (clouds.dV .<dv1i) )  .&& ((clouds.VC .>vc0i) .&& (clouds.VC .<vc1i)) .&& ((clouds.R .>R0i) .&& (clouds.R .<R1i)) .&& ((clouds.r .>r0i) .&& (clouds.r .<r1i)) .&& ((clouds.Z .>z0i) .&& (clouds.Z .<z1i))
			#@show sum(mask_rec)
        end

		@info "$(sum(mask_rec)) cloudlets in mask"
		mask =mask .|| (mask_rec)
		#@show sum(mask)

		# mask_rec= ((clouds.Xc .>x0i) .&& (clouds.Xc .<x1i)) .&& ((clouds.Fc .>f0i) .&& (clouds.Fc .<f1i)) .&& ((clouds.Yc .>y0i) .&& (clouds.Yc .<y1i)) .&& ((clouds.Vc .>v0i) .&& (clouds.Vc .<v1i)) .&& ((clouds.Sc .>s0i) .&& (clouds.Sc .<s1i))

		# maskc=maskc .|| (.!mask_rec)
    end
	if neg
		mask=.! mask
		@show sum(mask)
		maskc=.! maskc
	end

    @info "Found $(sum(mask)) unique sub-clouds and $(sum(maskc)) merged clouds"

    if length(clouds.r)>0 && length(clouds.ϕ)>0
        return GalaxyWarp.Clouds(Xp=clouds.Xp[mask],Yp=clouds.Yp[mask],Vp=clouds.Vp[mask],Fp=clouds.Fp[mask],Sp=clouds.Sp[mask],Ip=clouds.Ip[mask],logs=merge(clouds.logs,Dict{String, Any}("filtered"=>true)),belongs=clouds.belongs[mask],I=clouds.I[maskc],k=clouds.k[mask],Xc=clouds.Xc[maskc],Yc=clouds.Yc[maskc],Vc=clouds.Vc[maskc],Fc=clouds.Fc[maskc],Sc=clouds.Sc[maskc],Ic=clouds.Ic[maskc],r=clouds.r[mask],R=clouds.R[mask],Φ=clouds.Φ[mask],VC=clouds.VC[mask],dV=clouds.dV[mask],Z=clouds.Z[mask],inc=clouds.inc[mask],pa=clouds.pa[mask],dT=clouds.dT[mask],Tnet=clouds.Tnet[mask],ϕ=clouds.ϕ[mask],tilt=clouds.tilt[mask])
    elseif length(clouds.r)>0
        return GalaxyWarp.Clouds(Xp=clouds.Xp[mask],Yp=clouds.Yp[mask],Vp=clouds.Vp[mask],Fp=clouds.Fp[mask],Sp=clouds.Sp[mask],Ip=clouds.Ip[mask],logs=merge(clouds.logs,Dict{String, Any}("filtered"=>true)),belongs=clouds.belongs[mask],I=clouds.I[maskc],k=clouds.k[mask],Xc=clouds.Xc[maskc],Yc=clouds.Yc[maskc],Vc=clouds.Vc[maskc],Fc=clouds.Fc[maskc],Sc=clouds.Sc[maskc],Ic=clouds.Ic[maskc],r=clouds.r[mask],R=clouds.R[mask],Φ=clouds.Φ[mask],VC=clouds.VC[mask],dV=clouds.dV[mask],Z=clouds.Z[mask],inc=clouds.inc[mask],pa=clouds.pa[mask],dT=clouds.dT[mask],Tnet=clouds.Tnet[mask])
    else
        return GalaxyWarp.Clouds(Xp=clouds.Xp[mask],Yp=clouds.Yp[mask],Vp=clouds.Vp[mask],Fp=clouds.Fp[mask],Sp=clouds.Sp[mask],Ip=clouds.Ip[mask],logs=merge(clouds.logs,Dict{String, Any}("filtered"=>true)))
    end
end
#


text:
module HMM2d

#=
2 Dimensional Hidden Markov Models
=#

#include("hmm2d.jl")
#include("viterbi.jl")
#include("mle.jl")

end # module


text:
##############################################################################
#
# Symmetric triangular distribution from Distributions Handbook
#
##############################################################################

immutable Triangular <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function Triangular(l::Real, s::Real)
        if s > 0.0
            new(float64(l), float64(s))
        else
            error("scale must be positive")
        end
    end
end

Triangular(location::Real) = Triangular(location, 1.0)
Triangular() = Triangular(0.0, 1.0)

function insupport(d::Triangular, x::Number)
    return isreal(x) && isfinite(x) &&
           d.location - d.scale <= x <= d.location + d.scale
end

kurtosis(d::Triangular) = d.scale^4 / 15.0

mean(d::Triangular) = d.location

median(d::Triangular) = d.location

modes(d::Triangular) = [d.location]

function pdf(d::Triangular, x::Real)
    if insupport(d, x)
        return -abs(x - d.location) / (d.scale^2) + 1.0 / d.scale
    else
        return 0.0
    end
end

function rand(d::Triangular)
    xi1, xi2 = rand(), rand()
    return d.location + (xi1 - xi2) * d.scale
end

skewness(d::Triangular) = 0.0

var(d::Triangular) = d.scale^2 / 6.0


text:
## ------------------------------------------------------------------------
cd(@__DIR__)
# Install unregistered packages
using Pkg
Pkg.activate(Base.current_project(@__DIR__))
try
    pkg"rm Chemostat"
    pkg"rm UtilsJL"
catch; end
pkg"add https://github.com/josePereiro/UtilsJL.git#master"
pkg"add https://github.com/josePereiro/Chemostat#9b2070a"
pkg"instantiate"
pkg"build"
pkg"test Chemostat"


text:
module GeneticBitArrays

  import Base.show,
         Base.length,
         Base.getindex,
         Base.setindex!,
         Base.==,
         Base.rand,
         StatsBase.Weights,
         StatsBase.sample,
         StaticArrays.SVector

  abstract type GeneticSeq end

  struct DNASeq <: GeneticSeq
    data::BitArray{2}

    function DNASeq(x::BitArray{2}; checkinput::Bool=true)
      if checkinput
        if size(x, 1) == 4
          return new(x)
        else
          throw(ErrorException("Invalid input, must be `BitArray{4,n}`"))
        end
      else
        return new(x)
      end
    end
  end

  struct RNASeq <: GeneticSeq
    data::BitArray{2}

    function RNASeq(x::BitArray{2}; checkinput::Bool=true)
      if checkinput
        if size(x, 1) == 4
          return new(x)
        else
          throw(ErrorException("Invalid input, must be `BitArray{4,n}`"))
        end
      else
        return new(x)
      end
    end
  end

  function length(x::GeneticSeq)
    return size(x.data, 2)
  end

  function getindex(x::T, i) where {T <: GeneticSeq}
    return T(x.data[:, i])
  end

  const _bitslookup = SVector{16, BitArray{1}}([1;1;1;1], [1;1;1;0],
                                               [1;1;0;1], [1;1;0;0],
                                               [1;0;1;1], [1;0;1;0],
                                               [1;0;0;1], [1;0;0;0],
                                               [0;1;1;1], [0;1;1;0],
                                               [0;1;0;1], [0;1;0;0],
                                               [0;0;1;1], [0;0;1;0],
                                               [0;0;0;1], [0;0;0;0])

  const _dnacharlookup = SVector{16, Char}('N', 'V', 'H', 'M',
                                           'D', 'R', 'W', 'A',
                                           'B', 'S', 'Y', 'C',
                                           'K', 'G', 'T', '-')

  const _rnacharlookup = SVector{16, Char}('N', 'V', 'H', 'M',
                                           'D', 'R', 'W', 'A',
                                           'B', 'S', 'Y', 'C',
                                           'K', 'G', 'U', '-')

  const _lookup(::Type{DNASeq}) = _dnacharlookup
  const _lookup(::Type{RNASeq}) = _rnacharlookup
  const _seq(::Type{DNASeq}) = "DNA"
  const _seq(::Type{RNASeq}) = "RNA"

  function _bitarray(::Type{T}, x::Union{String, Vector{Char}}) where {T <: GeneticSeq}
    l = _lookup(T)
    s = BitArray{2}(undef, (4, length(x)))
    for i = eachindex(x)
      ind = findfirst(l .== x[i])
      if ind === nothing
        throw(ErrorException("Unrecognized $(_seq(T)) `Char` $(x[i]) at index $i"))
      else
        s[:, i] = _bitslookup[ind]
      end
    end
    return s
  end

  function _bitarray(::Type{T}, x::Char) where {T <: GeneticSeq}
    ind = findfirst(_lookup(T) .== x)
    if ind === nothing
      throw(ErrorException("Unrecognized $(_seq(T)) `Char` $x"))
    end
    return _bitslookup[ind]
  end

  function _bitarray(::Type{T}, x::BitArray{1}) where {T <: GeneticSeq}
    if length(x) != 4
      throw(ErrorException("Invalid input, must have a length of 4"))
    end
    return reshape(x, (4,1))
  end

  DNASeq(x) = DNASeq(_bitarray(DNASeq, x), checkinput=false)
  RNASeq(x) = RNASeq(_bitarray(RNASeq, x), checkinput=false)

  function setindex!(x::T, a, i) where {T <: GeneticSeq}
    x.data[:,i] = _bitarray(T, a)
    return x
  end

  function convert(::Type{String}, x::T) where {T <: GeneticSeq}
    s = ""
    for i in 1:length(x)
      if x.data[1,i]
        if x.data[2,i]
          if x.data[3,i]
            if x.data[4,i]
              # A C G T/U
              s *= "N"
            else
              # A C G
              s *= "V"
            end
          else
            if x.data[4,i]
              # A C T/U
              s *= "H"
            else
              # A C
              s *= "M"
            end
          end
        else
          if x.data[3,i]
            if x.data[4,i]
              # A G T/U
              s *= "D"
            else
              # A G
              s *= "R"
            end
          else
            if x.data[4,i]
              # A T/U
              s *= "W"
            else
              # A
              s *= "A"
            end
          end
        end
      else
        if x.data[2,i]
          if x.data[3,i]
            if x.data[4,i]
              # C G T/U
              s *= "B"
            else
              # C G
              s *= "S"
            end
          else
            if x.data[4,i]
              # C T/U
              s *= "Y"
            else
              # C
              s *= "C"
            end
          end
        else
          if x.data[3,i]
            if x.data[4,i]
              # G T/U
              s *= "K"
            else
              # G
              s *= "G"
            end
          else
            if x.data[4,i]
              if T == DNASeq
                # T
                s *= "T"
              else
                # U
                s *= "U"
              end
            else
              # gap
              s *= "-"
            end
          end
        end
      end
    end
    return s
  end

  function show(io::IO, x::T) where {T <: GeneticSeq}
    len = length(x)
    println(io, "$(len)nt $(_seq(T)) sequence")
    if len <= 26
      print(io, convert(String, x))
    else
      print(io, convert(String, x[1:13]) *
                "..." *
                convert(String, x[len-13:len]))
    end
  end

  function ==(x::T, y::T) where {T <: GeneticSeq}
    return x.data == y.data
  end

  function rand(::Type{T}, w::W, n::Int64; checkinput::Bool=true) where {T <: GeneticSeq, W<: Weights}
    if checkinput && length(w) != 4
      throw(ErrorException("Invalid sampling weights for $T generation"))
    end
    x = BitArray(fill(0, (4, n)))
    @simd for i = 1:n
      @inbounds x[sample(1:4, w), i] = 1
    end
    return T(x, checkinput=false)
  end

  function rand(::Type{T}, wv::Vector{W}; checkinput::Bool=true) where {T <: GeneticSeq, W<: Weights}
    len = length(wv)
    x = BitArray(fill(0, (4, len)))
    @simd for i = 1:len
      if checkinput && length(wv[i]) != 4
        throw(ErrorException("Invalid sampling weights for $(i)th nt in $T generation"))
      end
      @inbounds x[sample(1:4, wv[i]), i] = 1
    end
    return T(x, checkinput=false)
  end

  export GeneticSeq, DNASeq, RNASeq, Weights

end # module


text:
"""
    Aqua.test_project_toml_formatting(packages)
"""
function test_project_toml_formatting(packages)
    @testset "$(result.label)" for result in analyze_project_toml_formatting(packages)
        @debug result.label result
        @test result ⊜ true
    end
end

analyze_project_toml_formatting(packages) =
    [_analyze_project_toml_formatting_1(path) for path in project_toml_files_in(packages)]

project_toml_files_in(path::AbstractString) = [path]
project_toml_files_in(m::Module) = project_toml_files_in(PkgId(m))
function project_toml_files_in(pkg::PkgId)
    srcpath = Base.locate_package(pkg)
    if srcpath === nothing
        # TODO: record this as a test failure?
        error("Package file and directory not found: $pkg")
    end
    dir = dirname(dirname(srcpath))
    paths = [project_toml_path(dir)[1]]
    p, found = project_toml_path(joinpath(dir, "test"))
    found && push!(paths, p)
    return paths
end

project_toml_files_in(iterable) =
    [path for x in iterable for path in project_toml_files_in(x)]

function _analyze_project_toml_formatting_1(path::AbstractString)
    label = path

    if !isfile(path)
        return LazyTestResult(label, "Path `$path` is not an existing file.", false)
    end

    original = read(path, String)
    return _analyze_project_toml_formatting_2(path, original)
end

function _analyze_project_toml_formatting_2(path::AbstractString, original)
    @debug "Checking TOML style: `$path`" Text(original)
    label = path

    prj = TOML.parse(original)
    formatted = sprint(print_project, prj)
    if splitlines(original) == splitlines(formatted)
        LazyTestResult(
            label,
            "Running `Pkg.resolve` on `$(path)` did not change the content.",
            true,
        )
    else
        diff = format_diff(
            "Original $(basename(path))" => original,
            "Pkg's output" => formatted,
        )
        LazyTestResult(
            label,
            """
            Running `Pkg.resolve` on `$(path)` will change the content.

            $diff
            """,
            false,
        )
    end
end


text:
using Neet, Test

@test 1 + 2 == 3
@test 2 + 2 != 5


text:
import Pkg

function main(repo_url::AbstractString;
              build_directory::AbstractString = pwd(),
              master_branch::AbstractString = "master",
              allow_unchanged_prerelease::Bool = true,
              allow_skipping_versions::Bool = false,
              gh_set_output::Bool = get(ENV, "GITHUB_ACTIONS", "") == "true",
              gh_set_output_io::IO = stdout)::Nothing
    head_version = get_head_version(repo_url; build_directory = build_directory)
    master_version = get_upstream_version(repo_url; master_branch = master_branch)
    @info("Master version: $(master_version)")
    @info("Head version: $(head_version)")
    _gh_set_output_println(gh_set_output, gh_set_output_io, "master_version", "v$(master_version)")
    _gh_set_output_println(gh_set_output, gh_set_output_io, "pr_version", "v$(head_version)")
    compare_versions(master_version,
                     head_version;
                     allow_unchanged_prerelease = allow_unchanged_prerelease,
                     allow_skipping_versions = allow_skipping_versions,
                     gh_set_output = gh_set_output,
                     gh_set_output_io = gh_set_output_io)
    return nothing
end


text:
import Term: make_logo
# using BenchmarkTools

print(make_logo())
make_logo()
# display(@benchmark make_logo())