text:
struct IndexedCounts{M<:AbstractMer,C<:CountMode}
    mer_index::Vector{M}            # Ordered list of kmer that contain counts
    count_names::Vector{Symbol}     # Names of the count vectors
    counts::Vector{Vector{UInt16}}  # Count vectors, each contains an entry per mer in the index
    name::Symbol                    # The name of this MerCounts
end

# Construction of empty MerCounts with a name
@inline function IndexedCounts{M,C}(name::Symbol) where {M<:AbstractMer,C<:CountMode}
    return IndexedCounts{M,C}(Vector{M}(), Vector{Symbol}(), Vector{Vector{UInt16}}(), name)
end

@inline function IndexedCounts(::Type{M}, name::Symbol, mode::C = CANONICAL) where {M<:AbstractMer,C<:CountMode}
    return IndexedCounts{M,C}(name)
end

# Construction from a sequence distance graph
function IndexedCounts(name::Symbol, sequences, mode::C = CANONICAL) where {M<:AbstractMer,C<:CountMode}
    @info string("Creating an empty kmer count datastore called ", name)
    mc = IndexedCounts{M,C}(name, mode)
    index!(mc, sequences)
    return mc
end

function _count_and_collapse!(mers::Vector{M}, counts::Vector{UInt16}) where {M<:AbstractMer}
    wi = 1
    ri = 1
    stop = lastindex(mers) + 1
    @inbounds while ri < stop
        mers[wi] = mers[ri]
        ci = one(UInt16)
        while (ri += 1) < stop && mers[ri] == mers[wi]
            ci = ci + one(UInt16)
        end
        push!(counts, ci)
        wi = wi + 1
    end
    resize!(mers, wi - 1)
    return nothing
end

function index!(mc::IndexedCounts{M,C}, ref_sequences) where {M<:AbstractMer,C<:CountMode}
    @assert eltype(ref_sequences) <: LongSequence
    @info "Indexing kmer counts of reference sequences"
    @info "Populating index with reference sequence mers"
    # Add all Kmers from reference sequences.
    # t = expected number of kmers.
    t = sum(length(x) for x in ref_sequences) - BioSequences.ksize(M) + 1
    mcindex = mc.mer_index
    resize!(mcindex, t)
    vec_i = 1
    f = C()
    for seq in ref_sequences
        for mer in each(M, seq)
            mcindex[vec_i] = f(mer)
            vec_i = vec_i + 1
        end
    end
    resize!(mcindex, vec_i)
    @info "Sorting the index"
    # Sort the index
    sort!(mcindex)
    @info "Counting and collapsing mers in index"
    # Create a first count;
    # Collapse the kmer index, saving the coverage to the first count.
    first_count = Vector{UInt16}()
    sizehint!(first_count, length(mcindex))
    push!(mc.counts, first_count)
    push!(mc.count_names, :sdg)
    _count_and_collapse!(mcindex, first_count)
    return mc
end

function Base.summary(io::IO, mc::IndexedCounts{M}) where {M}
    print(io, BioSequences.ksize(M), "-mer Indexed Counts Datastore '",  mc.name, "': ", length(mc.counts), " stored ", BioSequences.ksize(M), "-mer counts")
end

function Base.show(io::IO, mc::IndexedCounts{M}) where {M}
    summary(io, mc)
end

"""
    add_count!(mc::IndexedCounts{M}, name::Symbol, reads::ReadDatastore) where {M}

Count the k-mers in a ReadDatastore and add them to the IndexedCounts{M}.

Currently IndexedCounts use a serial and in-memory counting method, as

!!! note
    Only the k-mers in the IndexedCounts' index shall be counted from the reads.
"""
function add_count!(mc::IndexedCounts{M,C}, name::Symbol, reads::ReadDatastore) where {M,C}
    mode = C()
    # TODO: Check if name already exists.
    push!(mc.count_names, name)
    c = zeros(UInt16, length(mc.mer_index))
    push!(mc.counts, c)
    # Populate lookup map.
    kmer_map = Dict{M, Int}()
    sizehint!(kmer_map, length(mc.mer_index))
    for (i, j) in enumerate(mc.mer_index)
        kmer_map[j] = i
    end
    for read in reads
        for mer in each(M, read)
            findk = get(kmer_map, mode(mer), 0)
            @inbounds if findk > 0
                old = c[findk]
                inc = old + one(UInt16)
                c[findk] = ifelse(old > inc, typemax(UInt16), inc)
            end
        end
    end
end

function project_count(mc::IndexedCounts{M,C}, name::Symbol, sequence::BioSequence) where {M,C}
    countidx = findfirst(isequal(name), mc.count_names)
    if isnothing(countidx)
        throw(ArgumentError(string("No count dataset called", string(name), "exists")))
    end
    counts = mc.counts[countidx]
    querymers = collect_mers(M, C(), sequence)
    projection = zeros(UInt16, length(querymers))
    indexsize = length(mc.mer_index)
    index = mc.mer_index
    for (i, mer) in enumerate(querymers)
        kidx = min(searchsortedfirst(index, mer), indexsize)
        if index[kidx] === mer
            projection[i] = counts[kidx]
        end
    end
    return projection
end


text:
using Distributions, Random

@testset "test kernel function" begin
    @test DualOptimalFiltering.dirichletkernel_oneval([0.2,0.3,0.5], [0.16,  0.64,  0.2], 1.2) ≈ 2.3904148155639433 atol=10.0^(-5)
    @test DualOptimalFiltering.dirichletkernel_oneval([0.0,0.5,0.5], [0.16,  0.64,  0.2], 1.2) == 0
    # Random.seed!(1)
    X  = rand(Dirichlet(ones(4)),10)'
    w = ones(size(X, 1))
    @test_nowarn dirichletkernel([0.2,0.3,0.4, 0.1], X, 1.2, w, size(X, 1))
    # res = [7.85877, 7.34332, 8.70397, 8.05658, 4.97756, 8.25494, 8.63194, 7.72478, 7.36807, 8.70078]
    # for i in eachindex(w)
    #     @test w[i] .≈ res[i] atol=10.0^(-5)
    # end
    X_list = [X[i,:] for i in 1:size(X,1)]
    @test_nowarn  dirichletkernel([0.2,0.3,0.4, 0.1], X, 1.2, w, size(X, 1))
    # for i in eachindex(w)
    #     @test w[i] .≈ res[i] atol=10.0^(-5)
    # end
end;

@testset "test bandwidth selection functions" begin
    Random.seed!(1)
    X  = rand(Dirichlet(ones(4)),10)'
    @test DualOptimalFiltering.midrange([0,1]) ≈ 0.5
    @test_nowarn DualOptimalFiltering.lcv(X, dirichletkernel, DualOptimalFiltering.midrange(X), ones(size(X, 1)), size(X, 1))
    @test_nowarn DualOptimalFiltering.bwlcv(X, dirichletkernel)
    @test_nowarn DualOptimalFiltering.bwlcv_large_bounds(X, dirichletkernel)
    @test_nowarn DualOptimalFiltering.minus_log_leaveoneout(X, dirichletkernel, DualOptimalFiltering.midrange(X), ones(size(X, 1)), size(X, 1))
    @test_nowarn DualOptimalFiltering.bwloo_large_bounds(X, dirichletkernel)
    X_list = [X[i,:] for i in 1:size(X,1)]
    @test_nowarn DualOptimalFiltering.midrange(X_list)
    @test_nowarn  DualOptimalFiltering.lcv(X_list, dirichletkernel, DualOptimalFiltering.midrange(X_list), ones(size(X, 1)), size(X, 1))
    @test_nowarn DualOptimalFiltering.bwlcv(X_list, dirichletkernel)
end;

ext: jl

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using Distributions
using EMpht
using LinearAlgebra: norm
using Random
using Test

# Fit a simple exponential distribution using
# both methods (uniformization and the ODE solver)
trueRate = 10
trueDist = Exponential(1/trueRate)

Random.seed!(1)
sObs = EMpht.Sample(obs=rand(trueDist, 1_000))
ph1unif = empht(sObs, p=1, method=:unif)
ph1ode = empht(sObs, p=1, method=:ode)

@test -ph1unif.T[1,1] ≈ trueRate atol=1
@test -ph1ode.T[1,1] ≈ trueRate atol=1

# Fit binned samples from a gamma distribution.
trueDist = Gamma(20, 2/10)
xGrid = range(0, 8, length=1_000)
truePDFs = pdf.(trueDist, xGrid)

# Random.seed!(1)
# data = rand(trueDist, 1_000)
# int, intweight = bin_observations(data, 15)
int = [1.5 2.0; 2.0 2.5; 2.5 3.0; 3.0 3.5; 3.5 4.0; 4.0 4.5; 4.5 5.0; 5.0 5.5;
        5.5 6.0; 6.0 6.5; 6.5 7.0; 7.0 7.5]
intweight = [4.0, 34.0, 107.0, 170.0, 202.0, 222.0, 140.0, 77.0, 24.0, 14.0,
        4.0, 2.0]
sInt = EMpht.Sample(int=int, intweight=intweight)

phCF1 = empht(sInt, p=100, ph_structure="CanonicalForm1")
fitPDFs = pdf.(phCF1, xGrid)

@test norm(truePDFs .- fitPDFs, Inf) <= 1e-1

# Test that reading settings from a file works.
ph100 = empht("Coxian100.json")

# Test that some of the internal functions work as expected.
π=[0.005632862236468221, 0.00898302181442919, 0.007418596691981033,
        0.020026201266915976, 0.04062880012941406, 0.05812175206666447,
        0.04260242040415281, 0.1654039105298598, 0.3888838419368729,
        0.2622985929232415]
T = [-1.7486539983379394 1.7486539983379394 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 -1.9875507782505892 1.9875507782505892 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 -1.9962263354076826 1.9962263354076826 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 -2.064044688471394 2.064044688471394 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 -2.1489226352793827 2.1489226352793827 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 -2.2809066095317565 2.2809066095317565 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 -2.5379722045000195 2.5379722045000195 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.763746562995887 2.763746562995887 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.082614321671955 3.082614321671955;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.910958867504505]
t=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.910958867504505]
ph = EMpht.PhaseType(π, T, t)

u = zeros(ph.p * ph.p)
du = similar(u)
EMpht.ode_observations!(du, u, ph, 10.0)

duCorrect = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.60666e-10, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.8095e-9, 0.0, 0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0, 2.46331e-8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        0.0, 1.15854e-7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        4.26538e-7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.23682e-6,
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.71623e-6, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.29954e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 8.87025e-6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
        1.093e-5]
@test norm(du .- duCorrect, Inf) <= 1e-10

EMpht.ode_censored!(du, u, ph, 10.0)

duCorrect = [1.43358e-10, 1.43358e-10, 1.43358e-10, 1.43358e-10, 1.43358e-10,
        1.43358e-10, 1.43358e-10, 1.43358e-10, 1.43358e-10, 1.43358e-10,
        9.74057e-10, 9.74057e-10, 9.74057e-10, 9.74057e-10, 9.74057e-10,
        9.74057e-10, 9.74057e-10, 9.74057e-10, 9.74057e-10, 9.74057e-10,
        6.29847e-9, 6.29847e-9, 6.29847e-9, 6.29847e-9, 6.29847e-9, 6.29847e-9,
        6.29847e-9, 6.29847e-9, 6.29847e-9, 6.29847e-9, 2.96229e-8, 2.96229e-8,
        2.96229e-8, 2.96229e-8, 2.96229e-8, 2.96229e-8, 2.96229e-8, 2.96229e-8,
        2.96229e-8, 2.96229e-8, 1.09062e-7, 1.09062e-7, 1.09062e-7, 1.09062e-7,
        1.09062e-7, 1.09062e-7, 1.09062e-7, 1.09062e-7, 1.09062e-7, 1.09062e-7,
        3.16245e-7, 3.16245e-7, 3.16245e-7, 3.16245e-7, 3.16245e-7, 3.16245e-7,
        3.16245e-7, 3.16245e-7, 3.16245e-7, 3.16245e-7, 6.94519e-7, 6.94519e-7,
        6.94519e-7, 6.94519e-7, 6.94519e-7, 6.94519e-7, 6.94519e-7, 6.94519e-7,
        6.94519e-7, 6.94519e-7, 1.35505e-6, 1.35505e-6, 1.35505e-6, 1.35505e-6,
        1.35505e-6, 1.35505e-6, 1.35505e-6, 1.35505e-6, 1.35505e-6, 1.35505e-6,
        2.26805e-6, 2.26805e-6, 2.26805e-6, 2.26805e-6, 2.26805e-6, 2.26805e-6,
        2.26805e-6, 2.26805e-6, 2.26805e-6, 2.26805e-6, 2.7947e-6, 2.7947e-6,
        2.7947e-6, 2.7947e-6, 2.7947e-6, 2.7947e-6, 2.7947e-6, 2.7947e-6,
        2.7947e-6, 2.7947e-6]

@test norm(du .- duCorrect, Inf) <= 1e-10

cIntCorrect = [1.95424e-7, 1.37498e-7, 8.25671e-8, 3.9341e-8, 1.42689e-8,
        3.6283e-9, 5.76897e-10, 5.8746e-11, 2.60493e-12, 1.13063e-14,
        1.09161e-6, 7.68042e-7, 4.61207e-7, 2.19753e-7, 7.97037e-8, 2.02671e-8,
        3.22246e-9, 3.28146e-10, 1.45507e-11, 6.3155e-14, 4.93283e-6,
        3.47068e-6, 2.08413e-6, 9.93034e-7, 3.6017e-7, 9.15845e-8, 1.45619e-8,
        1.48285e-9, 6.57528e-11, 2.85389e-13, 1.56191e-5, 1.09894e-5,
        6.59911e-6, 3.1443e-6, 1.14043e-6, 2.89989e-7, 4.6108e-8, 4.69522e-9,
        2.08197e-10, 9.03644e-13, 3.88217e-5, 2.73145e-5, 1.64023e-5,
        7.81523e-6, 2.83456e-6, 7.20775e-7, 1.14603e-7, 1.16701e-8, 5.17479e-10,
        2.24603e-12, 7.8674e-5, 5.53541e-5, 3.324e-5, 1.58379e-5, 5.74438e-6,
        1.46069e-6, 2.32248e-7, 2.365e-8, 1.0487e-9, 4.55169e-12, 0.000128921,
        9.07072e-5, 5.44694e-5, 2.59532e-5, 9.41315e-6, 2.39359e-6, 3.80578e-7,
        3.87547e-8, 1.71847e-9, 7.45873e-12, 0.000195269, 0.000137389,
        8.25018e-5, 3.93099e-5, 1.42576e-5, 3.62543e-6, 5.7644e-7, 5.86995e-8,
        2.60287e-9, 1.12973e-11, 0.000264903, 0.000186383, 0.000111922,
        5.33279e-5, 1.93419e-5, 4.91827e-6, 7.82001e-7, 7.9632e-8, 3.53106e-9,
        1.5326e-11, 0.000282308, 0.000198628, 0.000119276, 5.68317e-5,
        2.06127e-5, 5.24141e-6, 8.3338e-7, 8.4864e-8, 3.76306e-9, 1.63329e-11]

@test norm(EMpht.c_integrand(5.0, ph, 10.0) .- cIntCorrect, Inf) <= 1e-8

dIntCorrect = [2.49838e-7, 1.38082e-7, 6.89019e-8, 2.75403e-8, 8.48015e-9,
        1.84014e-9, 2.4941e-10, 2.21039e-11, 8.4793e-13, 2.89092e-15,
        1.39556e-6, 7.71302e-7, 3.84875e-7, 1.53836e-7, 4.73688e-8, 1.02787e-8,
        1.39317e-9, 1.23469e-10, 4.7364e-12, 1.61482e-14, 6.30634e-6,
        3.48541e-6, 1.7392e-6, 6.95164e-7, 2.14053e-7, 4.64482e-8, 6.29553e-9,
        5.57939e-10, 2.14032e-11, 7.29717e-14, 1.99681e-5, 1.1036e-5,
        5.50693e-6, 2.20114e-6, 6.77769e-7, 1.47071e-7, 1.99339e-8, 1.76663e-9,
        6.77701e-11, 2.31054e-13, 4.96313e-5, 2.74304e-5, 1.36876e-5,
        5.47098e-6, 1.68461e-6, 3.6555e-7, 4.95462e-8, 4.39101e-9, 1.68444e-10,
        5.74291e-13, 0.00010058, 5.5589e-5, 2.77386e-5, 1.10872e-5, 3.41395e-6,
        7.40805e-7, 1.00408e-7, 8.8986e-9, 3.41361e-10, 1.16383e-12,
        0.000164818, 9.10922e-5, 4.54545e-5, 1.81683e-5, 5.59434e-6, 1.21394e-6,
        1.64535e-7, 1.45819e-8, 5.59378e-10, 1.90714e-12, 0.000249641,
        0.000137972, 6.88474e-5, 2.75185e-5, 8.47344e-6, 1.83868e-6, 2.49213e-7,
        2.20864e-8, 8.47259e-10, 2.88863e-12, 0.000338663, 0.000187174,
        9.33986e-5, 3.73317e-5, 1.14951e-5, 2.49436e-6, 3.38083e-7, 2.99624e-8,
        1.14939e-9, 3.91873e-12, 0.000360914, 0.000199471, 9.95351e-5,
        3.97845e-5, 1.22503e-5, 2.65825e-6, 3.60296e-7, 3.1931e-8, 1.22491e-9,
        4.1762e-12]

@test norm(EMpht.d_integrand(5.0, ph, 10.0) .- dIntCorrect, Inf) <= 1e-8

p = ph.p; Bs = zeros(p); Zs = zeros(p); Ns = zeros(p, p+1)

EMpht.e_step_observed_ode!(sObs, ph, Bs, Zs, Ns)

BsCorrect = [9.02369e-6, 9.24901e-5, 0.00041208, 0.00568032, 0.0558114,
        0.377397, 1.30242, 23.3691, 251.208, 723.681]
ZsCorrect = [7.43742e-7, 8.80107e-6, 4.71973e-5, 0.000600675, 0.00613343,
        0.04273, 0.160008, 2.10905, 21.4158, 71.7403]
NsCorrect = [0.0 9.02399e-6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.000101516 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.000513599 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0061939 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.062005 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.439402 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.74183 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 25.1109 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 276.319 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1000.0]

@test norm(Bs .- BsCorrect, Inf) <= 1e-3
@test norm(Zs .- ZsCorrect, Inf) <= 1e-3
@test norm(Ns .- NsCorrect, Inf) <= 1e-3

p = ph.p; Bs = zeros(p); Zs = zeros(p); Ns = zeros(p, p+1)

EMpht.e_step_censored_ode!(sInt, ph, Bs, Zs, Ns)

BsCorrect = [130.738, 174.776, 108.273, 188.54, 206.27, 126.544, 29.5952,
        26.7001, 8.34855, 0.214886]
ZsCorrect = [79.0644, 173.84, 245.305, 376.65, 533.48, 617.446, 565.593,
        538.451, 475.918, 341.487]
NsCorrect = [0.0 130.737 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 305.513 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 413.787 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 602.329 0.0 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 808.599 0.0 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 935.148 0.0 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 964.712 0.0 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 991.49 0.0 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 999.738 0.0;
        0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1000.0]

@test norm(Bs .- BsCorrect, Inf) <= 1e-3
@test norm(Zs .- ZsCorrect, Inf) <= 1e-3
@test norm(Ns .- NsCorrect, Inf) <= 1e-3


text:
struct Population
  risks::AbstractDataFrame
  distances::Union{Nothing, AbstractArray}
  individuals::Int64

  Population(risks, distances) = new(risks, distances, size(risks, 1))
end

Population(risks) = Population(risks, nothing)

function Base.show(io::IO, x::Population)
  return print(io, "Population object (n=$(x.individuals))")
end


text:
# ---
# title: 1562. Find Latest Group of Size M
# id: problem1562
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/find-latest-group-of-size-m/description/>
# hidden: true
# ---
#
# Given an array `arr` that represents a permutation of numbers from `1` to `n`.
# You have a binary string of size `n` that initially has all its bits set to
# zero.
#
# At each step `i` (assuming both the binary string and `arr` are 1-indexed)
# from `1` to `n`, the bit at position `arr[i]` is set to `1`. You are given an
# integer `m` and you need to find the latest step at which there exists a group
# of ones of length `m`. A group of ones is a contiguous substring of 1s such
# that it cannot be extended in either direction.
#
# Return _the latest step at which there exists a group of ones of length
# **exactly**_  `m`. _If no such group exists, return_  `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,5,1,2,4], m = 1
#     Output: 4
#     Explanation: Step 1: "00 _1_ 00", groups: ["1"]
#     Step 2: "0010 _1_ ", groups: ["1", "1"]
#     Step 3: " _1_ 0101", groups: ["1", "1", "1"]
#     Step 4: "1 _1_ 101", groups: ["111", "1"]
#     Step 5: "111 _1_ 1", groups: ["11111"]
#     The latest step at which there exists a group of size 1 is step 4.
#
# **Example 2:**
#
#
#
#     Input: arr = [3,1,5,4,2], m = 2
#     Output: -1
#     Explanation: Step 1: "00 _1_ 00", groups: ["1"]
#     Step 2: " _1_ 0100", groups: ["1", "1"]
#     Step 3: "1010 _1_ ", groups: ["1", "1", "1"]
#     Step 4: "101 _1_ 1", groups: ["1", "111"]
#     Step 5: "1 _1_ 111", groups: ["11111"]
#     No group of size 2 exists during any step.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [1], m = 1
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: arr = [2,1], m = 2
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `n == arr.length`
#   * `1 <= n <= 10^5`
#   * `1 <= arr[i] <= n`
#   * All integers in `arr` are  **distinct**.
#   * `1 <= m <= arr.length`
#
#
## @lc code=start
using LeetCode

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


text:
# missing methods in TaylorSeries.jl:

degree(f::TaylorN) = f.order
setindex!(f::TaylorN, x, i) = f.coeffs[i+1] = x

doc"""
A `Taylor1Model` represents a polynomial approximation to a function $f(t)$ of a single variable.

The fields are:
- `n`: degree of the polynomial
- `x0`: expansion point as IntervalBox
- `I`: IntervalBox over which the Taylor model is defined / valid
- `p`: the polynomial, of type `TaylorSeries.TaylorN`
- `Δ`: the interval bound
"""
struct TaylorNModel{N,T,S}
    n::Int      # degree
    x0::IntervalBox{N,T}  # expansion point
    I::IntervalBox{N,T}   # interval over which the Taylor model is valid
    p::TaylorN{S}  # Taylor Taylor1nomial
    Δ::Interval{T}   # interval remainder bound
    order_bounds::Vector{Interval{T}}
end

# constructors that calculate order bounds automatically:
function TaylorNModel(n, x0, I::IntervalBox{N,T}, p, Δ) where {N,T}
    order_bounds = [homog([I...]) for homog in p.coeffs]  # bound homogeneous polynomials

    return TaylorNModel(n, x0, I, p, Δ, order_bounds)
end

TaylorNModel(n, x0, I::IntervalBox{N,T}, p) where {N,T} =  TaylorNModel(n, x0, I, p, Interval{T}(0))

# convenience constructors
TaylorNModel(f, p, Δ) = TaylorNModel(f.n, f.x0, f.I, p, Δ)
TaylorNModel(f, Δ) = TaylorNModel(f.n, f.x0, f.I, f.p, Δ)

include("arithmetic.jl")
# include("functions.jl")
# include("bound.jl")
#
include("integrate.jl")
# include("draw.jl")

import Base.copy
copy(f::TaylorNModel) = TaylorNModel(f.n, f.x0, f.I, copy(f.p), f.Δ, copy(f.order_bounds))

import Base.show
function show(io::IO, f::TaylorNModel)
    print(io,
    """Taylor model of degree $(f.n):
       - x0: $(f.x0)
       - I: $(f.I)
       - p: $(f.p)
       - Δ: $(f.Δ)
    """
    )
end


text:
##
x = 5
##
typeof(x)

##
"""
`add5`
======

This adds the value of 5 to the input variable.

Returns: the same type as the input variable + the integer 5.

Example
-------
        add5(5.0)
        typeof(add5(5.0)) # returns Float64
        add5(5)
        typeof(add5(5)) # returns Int64
"""
function add5(x)
    return x+5
end
@show typeof(add5(5.0))
@show add5(5.0)

##
@show typeof(add5(5))
@show add5(5)
##

## Types of Matrices and Vectors
# a vector in Julia looks like:

x = [5,6,7]
@show typeof(x)

##
x = [5.0,6,7]
@show typeof(x)

##
A = [5 6 7]
##
A = [5.0 6 7; 1 2 3; 4 5 6]
##
B = A .+ 5

##
C = B' # this is the transpose
@show typeof(C)
##
D = copy(C)
##
E = D+B

##


text:
import StatisticalRethinking: plotcoef

ProjDir = @__DIR__

models=[m5_1s, m5_2s, m5_3s];
pars=[:bA, :bM];
fig="$(ProjDir)/AMD_fork_1.png"
title="Particles (Normal) estimates"

function plotcoef(models::Vector{SampleModel}, pars::Vector{Symbol}, fig::AbstractString,
  title="", func=nothing)

  mnames = [models[i].name for i in 1:length(models)]
  levels = length(models) * (length(pars) + 1)
  colors = [:blue, :red, :green, :darkred, :black]

  s = Vector{NamedTuple}(undef, length(models))
  for (mindx, mdl) in enumerate(models)
    if isnothing(func)
      s[mindx] = read_samples(mdl; output_format=:particles)
    else
      dfs = read_samples(mdl; output_format=:dataframe)
      s[mindx] = func(dfs)
    end
  end

  xmin = 0; xmax = 0.0
  for i in 1:length(s)
    for par in pars
      syms = Symbol.(keys(s[i]))
      if Symbol(par) in syms
        mp = mean(s[i][Symbol(par)])
        sp = std(s[i][Symbol(par)])
        xmin = min(xmin, mp - sp)
        xmax = max(xmax, mp + sp)
      end
    end
  end

  ylabs = String[]
  for j in 1:length(models)
    for i in 1:length(pars)
      l = length(String(pars[i]))
      str = repeat(" ", 9-l) * String(pars[i])
      append!(ylabs, [str])
    end
    l = length(models[j].name)
    str = models[j].name * repeat(" ", 9-l)
    append!(ylabs, [str])
  end

  ys = [string(ylabs[i]) for i = 1:length(ylabs)]
  plot(xlims=(xmin, xmax), leg=false, framestyle=:grid)
  title!(title)
  yran = range(1, stop=length(ylabs), length=length(ys))
  yticks!(yran, ys)

  line = 0
  for (mindx, model) in enumerate(models)
    line += 1
    hline!([line] .+ length(pars), color=:darkgrey, line=(2, :dash))
    for (pindx, par) in enumerate(pars)
      line += 1
      syms = Symbol.(keys(s[mindx]))
      if Symbol(par) in syms
        ypos = (line - 1)
        mp = mean(s[mindx][Symbol(par)])
        sp = std(s[mindx][Symbol(par)])
        plot!([mp-sp, mp+sp], [ypos, ypos], leg=false, color=colors[pindx])
        scatter!([mp], [ypos], color=colors[pindx])
      end
    end
  end
  savefig(fig)
  s
end

plotcoef(models, pars, fig, title)


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using ArgParse

ArgParse.parse_item(::Type{Vector{String}}, x::AbstractString) = split(x, ",")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
        "--min_nv"
            help = "Minimum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemin(Int)
        "--max_nv"
            help = "Maximum number of vertices a network had to have to be processed."
            arg_type = Int
            default = typemax(Int)
        "--recompute"
            help = "Wether to recompute the basic features if they're already not NULL in the database."
            default = false
            action = :store_true
        "--feature_names"
            help = "Feature(s) name(s) to compute"
            required = true
            arg_type = Vector{String}
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    db = SQLite.DB(parsed_args["dbpath"])
    save_single_features_instances!(db;
                                    feature_name=parsed_args["feature_names"],
                                    min_nv=parsed_args["min_nv"],
                                    max_nv=parsed_args["max_nv"],
                                    recompute=parsed_args["recompute"])
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using AdventOfCode2018
using AdventOfCode2018.Day7

input = get_aoc_input(7)
parsed_input = parse_input(input)
dag = build_dag(parsed_input)

println("Solution 1: ", part_1(dag))
println("Solution 2: ", part_2(dag, workers=5, per_task=60))


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbrect()::NURBS1D

or

    nrbrect(size::F)::NURBS1D where {F<:AbstractFloat}

or

    nrbrect(width::F,height::F)::NURBS1D where {F<:AbstractFloat}

Construct NURBS representation of a rectangle.

Construct a rectangle or square in the x-y plane with the bottom corner at
(0,0,0). If there are no arguments provided the function constructs a unit
square.

# Arguments:
- `size`: Size of the square (width = height).
- `width`: Width of the rectangle (along x-axis).
- `height`: Height of the rectangle (along y-axis).

# Output:
- `crv`: NURBS curve for a straight line.

# Examples:
```julia
julia> crv = nrbrect()
julia> crv = nrbrect(size)
julia> crv = nrbrect(width, height)
```
"""
function nrbrect(w::F,h::F)::NURBS1D where {F<:AbstractFloat}
coefs  = [0 w w w w 0 0 0;
          0 0 0 h h h h 0;
          0 0 0 0 0 0 0 0;
          1 1 1 1 1 1 1 1];
knots  = vec([0 0 0.25 0.25 0.5 0.5 0.75 0.75 1 1]);
return nrbmak(coefs, knots);
end # nrbrect

# overloading of function nrbrect with one arguments
function nrbrect(size::F)::NURBS1D where {F<:AbstractFloat}
return nrbrect(size,size);
end # nrbrect

# overloading of function nrbrect with zero arguments
function nrbrect()::NURBS1D
return nrbrect(1.0,1.0);
end # nrbrect


text:
using MetacommunityDynamics
using DynamicGrids
using Plots
using Distributions
using Dispersal: Moore, DispersalKernel
using ColorSchemes
using EcologicalNetworks: mpnmodel, trophic_level, UnipartiteNetwork, degree, richness

number_of_species = 200
connectance = 0.05
forbiddenlinkprob = 0.5

foodweb = mpnmodel(number_of_species, connectance, forbiddenlinkprob)

sp = DiscreteUnipartiteSpeciesPool(Symbol.(foodweb.S), Matrix(foodweb.edges)) # move these type changes to a method
trophicdict = trophic_level(foodweb)  # returns a dictionary

resource = filter(s -> trophicdict[String(s)] == 1, sp.species)
consumers = filter(s -> trophicdict[String(s)] > 1, sp.species)

consumermodel =
    FoodWebEating(consumers, resource, LotkaVolterra(0.2), metaweb(sp)) +
    AdjacentBernoulliDispersal(consumers, DispersalKernel(radius=1), 0.1) +
    RandomExtinction(consumers, probability=0.1) +
    LinearMortality(consumers, 0.01);

plantmodel =
    LogisticGrowth(resource) +
    AdjacentBernoulliDispersal(resource, DispersalKernel(radius=2), 0.1) ;

fullmodel = consumermodel + plantmodel;

dim = (50,50)
ntimesteps = 300

init = NamedTuple(merge(
    rand(Biomass, resource, Exponential(10), dim...),
    rand(Biomass, consumers, Exponential(10), dim...)));

arrayout = ArrayOutput(init, tspan=1:ntimesteps)
@time sim!(arrayout, fullmodel);

"""
    plots :
"""

mat = zeros(number_of_species, ntimesteps)

spnames = []
for t = 1:ntimesteps
  ind = 1

  for (k,v) in zip(keys(arrayout[t]), values(arrayout[t]))
    mat[ind, t] = sum(v)
    push!(spnames, k)
    ind += 1
  end
end

plt = plot(legend=:none)
xaxis!(plt, xticks=[i for i in 0:20:ntimesteps],"timestep")
yaxis!(plt, "log2(biomass)")
for s in 1:number_of_species
    speciesname = spnames[s]
    thiscol = get(ColorSchemes.thermal, trophicdict[String(speciesname)]/max(values(trophicdict)...))

    @show thiscol, speciesname, trophicdict[String(speciesname)]
    plot!(plt, 1:ntimesteps, mat[s,:], la=0.7, lw=3, c=thiscol, colorbar_title="trophic level")
end
plt

savefig(plt, "manyspeciestimeseries.png")


text:
u = 10 + 21

high(x::Int) = 10
high(x::Float64) = bar(x)
bar(y) = 2

module SubMod
double(x) = x+x
end

include("dummydir/D.jl")


text:
struct ContainerKey{N<:Real} <: ResourceKey
  priority :: Real
  id :: UInt
  amount :: N
end

mutable struct Container{N<:Real} <: AbstractResource
  env :: Environment
  capacity :: N
  level :: N
  seid :: UInt
  put_queue :: DataStructures.PriorityQueue{Put, ContainerKey{N}}
  get_queue :: DataStructures.PriorityQueue{Get, ContainerKey{N}}
  function Container{N}(env::Environment, capacity::N=one(N); level::N=zero(N)) where {N<:Real}
    new(env, capacity, level, zero(UInt), DataStructures.PriorityQueue{Put, ContainerKey{N}}(), DataStructures.PriorityQueue{Get, ContainerKey{N}}())
  end
end

function Container(env::Environment, capacity::N=one(N); level::N=zero(N)) where N<:Real
  Container{N}(env, capacity, level=level)
end

const Resource = Container{Int}

function put(con::Container{N}, amount::N; priority::Real=0) where N<:Real
  put_ev = Put(con.env)
  con.put_queue[put_ev] = ContainerKey(priority, con.seid+=one(UInt), amount)
  @callback trigger_get(put_ev, con)
  trigger_put(put_ev, con)
  put_ev
end

request(res::Resource; priority::Real=0) = put(res, 1; priority=priority)

function get(con::Container{N}, amount::N; priority::Real=0) where N<:Real
  get_ev = Get(con.env)
  con.get_queue[get_ev] = ContainerKey(priority, con.seid+=one(UInt), amount)
  @callback trigger_put(get_ev, con)
  trigger_get(get_ev, con)
  get_ev
end

release(res::Resource; priority::Real=0) = get(res, 1; priority=priority)

function do_put(con::Container{N}, put_ev::Put, key::ContainerKey{N}) where N<:Real
  con.level + key.amount > con.capacity && return false
  schedule(put_ev)
  con.level += key.amount
  true
end

function do_get(con::Container{N}, get_ev::Get, key::ContainerKey{N}) where N<:Real
  con.level - key.amount < zero(N) && return false
  schedule(get_ev)
  con.level -= key.amount
  true
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# https://twitter.com/dolphin7473/status/1466376598853148680

# %%
module O

using OffsetArrays

mutable struct Variables
    u::OffsetArray{ComplexF64, 1}
end

function update_u!(var::Variables)
    rhs = OffsetArray{ComplexF64}(undef, -2:2)
    rhs[-2:2] = randn(ComplexF64, 5)
    var.u .+= rhs
end

function main()
    u = OffsetArray{ComplexF64}(undef, -2:2)
    var = Variables(u)
    var.u[-2:2] .= 10:10:50
    @time update_u!(var)
end

end

O.main()
O.main()
O.main()

# %%
@code_warntype O.main()

# %%
module P

using OffsetArrays

mutable struct Variables{U<:OffsetArray{ComplexF64, 1}}
    u::U
end

function update_u!(var::Variables)
    rhs = OffsetArray{ComplexF64}(undef, -2:2)
    rhs[-2:2] = randn(ComplexF64, 5)
    var.u .+= rhs
end

function main()
    u = OffsetArray{ComplexF64}(undef, -2:2)
    var = Variables(u)
    var.u[-2:2] .= 10:10:50
    @time update_u!(var)
end

end

P.main()
P.main()
P.main()

# %%
@code_warntype P.main()

# %%
module Q

using OffsetArrays
using Random

struct Variables{U}
    u::U
end

function update_u!(var::Variables, rhs)
    randn!(rhs)
    var.u .+= rhs
    var
end

function main()
    u = OffsetVector{ComplexF64}(undef, -2:2)
    var = Variables(u)
    parent(var.u) .= 10:10:50
    rhs = similar(u)
    @time update_u!(var, rhs)
end

end

Q.main()
Q.main()
var = Q.main()
@show var
var.u

# %%
@code_warntype Q.main()

# %%
u = Q.OffsetVector{ComplexF64}(undef, -2:2)
u

# %%
typeof(u)

# %%

max_stars_repo_path: 0025/OffsetArrays example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0025/OffsetArrays example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0025/OffsetArrays example.jl

max_forks_repo_name: genkuroki/public


text:
@testset "Quadratic Costs" begin
    # Quadratic Costs
    n, m = rand(10:20), rand(5:10)
    Q = Diagonal(@SVector fill(0.1, n))
    R = Diagonal(@SVector fill(0.01, m))
    H = rand(m, n)
    q = @MVector rand(n)
    r = @MVector rand(m)
    c = rand()
    Qf = Diagonal(@SVector fill(10.0, n))
    xf = @SVector ones(n)

    @testset "Constructors" begin
        qcost = QuadraticCost(Q, R)
        @test qcost.Q == Q
        @test qcost.q == zeros(n)
        @test qcost.R == R
        @test qcost.r == zeros(m)
        @test qcost.c == 0
        @test state_dim(qcost) == n
        @test control_dim(qcost) == m
        @test TO.is_blockdiag(qcost) == true

        qcost = QuadraticCost(Q, R, H = H)
        @test qcost.H ≈ H
        @test qcost.Q == Q
        @test qcost.q == zero(q)

        qcost = QuadraticCost(Q, R, r = r, H = H, q = q, c = c, terminal = true)
        @test qcost.H ≈ H
        @test qcost.Q == Q
        @test qcost.q == q
        @test qcost.r == r
        @test qcost.c == c
        @test qcost.zeroH == false
        @test qcost.terminal == true

        qcost = QuadraticCost(Q, R, H, q, r, c, checks = false, terminal = false)
        @test qcost.H ≈ H
        @test qcost.Q == Q
        @test qcost.q == q
        @test qcost.r == r
        @test qcost.c == c
        @test qcost.zeroH == false
        @test qcost.terminal == false
        @test TO.is_blockdiag(qcost) == false

        qcost = QuadraticCost{Float64}(n, m)
        @test eltype(qcost.Q) == Float64
        @test qcost.Q == I(n)
        @test qcost.R == I(m)
        @test qcost.q == zero(q)
        @test qcost.r == zero(r)
        @test qcost.c == zero(c)

        qcost = QuadraticCost{Float32}(n, m)
        @test eltype(qcost.Q) == Float32
        @test eltype(qcost.r) == Float32
        @test eltype(qcost.c) == Float32

        @test_logs (:warn, "R is not positive definite") QuadraticCost(Q, R * 0, H, q, r, c)
        @test_nowarn QuadraticCost(Q * 0, R, H, q, r, c)
        @test_logs (:warn, "Q is not positive semidefinite") QuadraticCost(
            Q .- 0.2,
            R,
            H,
            q,
            r,
            c,
        )
        @test_logs (:warn, "R is not positive definite") QuadraticCost(Q, R * 0)
        @test_logs (:warn, "Q is not positive semidefinite") QuadraticCost(
            Q .- 0.2,
            R,
            q = q,
        )
        @test_nowarn QuadraticCost(Q .- 0.2, R, q = q, checks = false)

        # Diagonal Cost functions
        dcost = DiagonalCost(Q, R)
        @test dcost.Q === Q
        @test dcost.R === R
        @test dcost.r == zero(r)
        @test dcost.q == zero(q)

        dcost = DiagonalCost(Q, R, q = q)
        @test dcost.Q === Q
        @test dcost.R === R
        @test dcost.r == zero(r)
        @test dcost.q === q

        dcost = DiagonalCost(Q.diag, R.diag, q = q)
        @test dcost.Q === Q
        @test dcost.R === R
        @test dcost.r == zero(r)
        @test dcost.q === q

        dcost = DiagonalCost(Vector(Q.diag), Vector(R.diag), q = q)
        @test dcost.Q === Q
        @test dcost.R === R
        @test dcost.r == zero(r)
        @test dcost.q === q
        @test dcost.Q isa Diagonal{Float64,<:SVector}

        # Test convenience constructors
        quadcost = LQRCost(Q, R, xf)
        @test quadcost.Q == Q
        @test quadcost.q == -Q * xf
        @test quadcost.R == R
        @test quadcost.r == zeros(m)
        @test quadcost.c ≈ 0.5 * xf'Q * xf
        @test quadcost isa DiagonalCost
        @test TO.is_blockdiag(quadcost)
        @test TO.is_diag(quadcost)
    end

    @testset "Math Operations" begin
        # Test adding cost functions
        qcost = QuadraticCost(Q, R, H, q, r, c, checks = false, terminal = false)
        dcost = DiagonalCost(Vector(Q.diag), Vector(R.diag), q = q)
        addcost = dcost + qcost
        @test addcost.Q ≈ 2Q
        @test addcost.R ≈ 2R
        @test addcost.q ≈ 2q
        @test addcost.r ≈ r
        @test addcost.c ≈ c
        @test addcost isa QuadraticCost

        dcost = DiagonalCost(Vector(Q.diag), Vector(R.diag), q, r, c)
        addcost = dcost + dcost
        @test addcost.Q ≈ 2Q
        @test addcost.R ≈ 2R
        @test addcost.q ≈ 2q
        @test addcost.r ≈ 2r
        @test addcost.c ≈ 2c

        dcost2 = DiagonalCost(Q.diag, Q.diag, q, q, c)
        @test_throws AssertionError dcost2 + dcost

        # Test inversion
        dinv = inv(dcost)
        @test dinv.Q ≈ inv(Q)
        @test dinv.R ≈ inv(R)
        @test dinv.q == q

        qinv = inv(qcost)
        G = inv(Matrix([Q H'; H R]))
        @test G[1:n, 1:n] ≈ qinv.Q
        @test G[1:n, n.+(1:m)] ≈ qinv.H'
        @test G[n.+(1:m), n.+(1:m)] ≈ qinv.R
        @test dinv.q == q
        @test dinv.r == r
        @test dinv.c == c
        @test TO.is_diag(qinv) == false
        @test TO.is_diag(dinv) == true
        @test TO.is_blockdiag(qinv) == false
        @test TO.is_blockdiag(dinv) == true

        qinv = inv(QuadraticCost(Q, R))
        @test qinv.Q ≈ inv(Q)
        @test qinv.R ≈ inv(R)
        @test TO.is_blockdiag(qinv) == true
        @test TO.is_diag(qinv) == true

    end

    @testset "Conversion and Promotion" begin
        # Test promotion rules
        qcost = QuadraticCost(10Q, zero(R), q = q, c = c, terminal = true)
        dcost = DiagonalCost(2Q, zero(R), q = q, c = c, terminal = true)
        c1, c2 = promote(qcost, dcost)
        @test c1 isa QuadraticCost
        @test c2 isa QuadraticCost
        @test c1.Q ≈ qcost.Q
        @test c1.R ≈ qcost.R
        @test c2.Q ≈ dcost.Q

        c1, c2 = promote(dcost, dcost)
        @test c1 isa DiagonalCost
        @test c2 isa DiagonalCost
        @test c1.Q ≈ dcost.Q
        @test c1.R ≈ dcost.R
        @test c2.Q ≈ dcost.Q

        c1 = QuadraticCost(Q, R)
        c2 = QuadraticCost(2 * Matrix(Q), 2 * R)
        p1, p2 = promote(c1, c2)
        @test p1.Q ≈ Q
        @test p2.Q ≈ 2Q
        @test p1.R ≈ R
        @test p2.R ≈ 2R
        @test p1.Q isa SizedMatrix
        @test c1.Q isa Diagonal
        @test c2.Q isa Matrix
    end

    @testset "Evaluation and expansion" begin
        # Test cost functions and expansions
        x = @SVector rand(n)
        u = @SVector rand(m)
        dt = 0.1
        z = KnotPoint(x, u, dt)
        zterm = KnotPoint(x, u, 0.0)
        @test TO.is_terminal(zterm)

        qcost = QuadraticCost(Q, R, H, q, r, c)
        @test TO.stage_cost(qcost, z) ≈ TO.stage_cost(qcost, x, u) * dt
        @test TO.stage_cost(qcost, x, u) ≈
              0.5 * (x'Q * x + u'R * u) + q'x + r'u + c + u'H * x

        E = QuadraticCost{Float64}(n, m)
        TO.gradient!(E, qcost, zterm)
        @test E.q ≈ Q * x + q
        @test E.r ≈ zero(r)
        TO.gradient!(E, qcost, z)
        @test E.q ≈ Q * x + q + H'u
        @test E.r ≈ R * u + r + H * x

        TO.hessian!(E, qcost, zterm)
        @test E.Q ≈ Q
        @test E.R ≈ I(m)
        TO.hessian!(E, qcost, z)
        @test E.R ≈ R
        @test E.H ≈ H
        @test (@allocated TO.gradient!(E, qcost, zterm)) == 0
        # @test (@allocated TO.gradient!(E, qcost, z)) == 0
        @test (@allocated TO.hessian!(E, qcost, zterm)) == 0
        @test (@allocated TO.hessian!(E, qcost, z)) == 0

        dcost = DiagonalCost(Q, R, q, r, c)
        @test TO.stage_cost(dcost, z) ≈ TO.stage_cost(dcost, x, u) * dt
        @test TO.stage_cost(dcost, x, u) ≈ 0.5 * (x'Q * x + u'R * u) + q'x + r'u + c

        E = QuadraticCost{Float64}(n, m)
        TO.gradient!(E, dcost, zterm)
        @test E.q ≈ Q * x + q
        @test E.r ≈ zero(r)
        TO.gradient!(E, dcost, z)
        @test E.q ≈ Q * x + q
        @test E.r ≈ R * u + r

        TO.hessian!(E, dcost, zterm)
        @test E.Q ≈ Q
        @test E.R ≈ I(m)
        TO.hessian!(E, dcost, z)
        @test E.R ≈ R
        @test E.H ≈ zero(H)
        @test (@allocated TO.gradient!(E, dcost, zterm)) == 0
        @test (@allocated TO.gradient!(E, dcost, z)) == 0
        @test (@allocated TO.hessian!(E, dcost, zterm)) == 0
        @test (@allocated TO.hessian!(E, dcost, z)) == 0
    end

end


text:
@testset "mvnormal.jl" begin
    @testset "consistency with Normal" begin
        nsamples = 1_000
        ntestsamples = 5

        # create predictions
        predictions_μ = randn(nsamples)
        predictions_σ = rand(nsamples)
        predictions_normal = map(Normal, predictions_μ, predictions_σ)
        predictions_mvnormal = map(predictions_μ, predictions_σ) do μ, σ
            MvNormal([μ], σ^2 * I)
        end

        # create targets
        targets_normal = randn(nsamples)
        targets_mvnormal = map(vcat, targets_normal)

        # create test locations
        testpredictions_μ = randn(ntestsamples)
        testpredictions_σ = rand(ntestsamples)
        testpredictions_normal = map(Normal, testpredictions_μ, testpredictions_σ)
        testpredictions_mvnormal = map(testpredictions_μ, testpredictions_σ) do μ, σ
            MvNormal([μ], σ^2 * I)
        end
        testtargets_normal = randn(ntestsamples)
        testtargets_mvnormal = map(vcat, testtargets_normal)

        for kernel in (
            ExponentialKernel(; metric=Wasserstein()) ⊗ SqExponentialKernel(),
            ExponentialKernel(; metric=Wasserstein()) ⊗
            (SqExponentialKernel() ∘ ScaleTransform(rand())),
            ExponentialKernel(; metric=Wasserstein()) ⊗
            (SqExponentialKernel() ∘ ARDTransform([rand()])),
        )
            for estimator in
                (BiasedSKCE(kernel), UnbiasedSKCE(kernel), BlockUnbiasedSKCE(kernel, 5))
                skce_mvnormal = estimator(predictions_mvnormal, targets_mvnormal)
                skce_normal = estimator(predictions_normal, targets_normal)
                @test skce_mvnormal ≈ skce_normal
            end

            ucme_mvnormal = UCME(kernel, testpredictions_mvnormal, testtargets_mvnormal)(
                predictions_mvnormal, targets_mvnormal
            )
            ucme_normal = UCME(kernel, testpredictions_normal, testtargets_normal)(
                predictions_normal, targets_normal
            )
            @test ucme_mvnormal ≈ ucme_normal
        end
    end

    @testset "kernels with input transformations" begin
        nsamples = 100
        ntestsamples = 5

        for dim in (1, 10)
            # create predictions and targets
            predictions = [MvNormal(randn(dim), rand() * I) for _ in 1:nsamples]
            targets = [randn(dim) for _ in 1:nsamples]

            # create random test locations
            testpredictions = [MvNormal(randn(dim), rand() * I) for _ in 1:ntestsamples]
            testtargets = [randn(dim) for _ in 1:ntestsamples]

            for γ in (1.0, rand())
                kernel1 =
                    ExponentialKernel(; metric=Wasserstein()) ⊗
                    (SqExponentialKernel() ∘ ScaleTransform(γ))
                kernel2 =
                    ExponentialKernel(; metric=Wasserstein()) ⊗
                    (SqExponentialKernel() ∘ ARDTransform(fill(γ, dim)))
                kernel3 =
                    ExponentialKernel(; metric=Wasserstein()) ⊗
                    (SqExponentialKernel() ∘ LinearTransform(diagm(fill(γ, dim))))

                # check evaluation of the first two observations
                p1 = predictions[1]
                p2 = predictions[2]
                t1 = targets[1]
                t2 = targets[2]
                for f in (
                    CalibrationErrors.unsafe_skce_eval_targets,
                    CalibrationErrors.unsafe_ucme_eval_targets,
                )
                    out1 = f(kernel1.kernels[2], p1, t1, p2, t2)
                    out2 = f(kernel2.kernels[2], p1, t1, p2, t2)
                    out3 = f(kernel3.kernels[2], p1, t1, p2, t2)
                    @test out2 ≈ out1
                    @test out3 ≈ out1
                    if isone(γ)
                        @test f(SqExponentialKernel(), p1, t1, p2, t2) ≈ out1
                    end
                end

                # check estimates
                for estimator in (UnbiasedSKCE, x -> UCME(x, testpredictions, testtargets))
                    estimate1 = estimator(kernel1)(predictions, targets)
                    estimate2 = estimator(kernel2)(predictions, targets)
                    estimate3 = estimator(kernel3)(predictions, targets)
                    @test estimate2 ≈ estimate1
                    @test estimate3 ≈ estimate1
                    if isone(γ)
                        @test estimator(
                            ExponentialKernel(; metric=Wasserstein()) ⊗
                            SqExponentialKernel(),
                        )(
                            predictions, targets
                        ) ≈ estimate1
                    end
                end
            end
        end
    end

    @testset "apply" begin
        dim = 10
        μ = randn(dim)
        A = randn(dim, dim)

        for d in (
            MvNormal(μ, rand() * I),
            MvNormal(μ, Diagonal(rand(dim))),
            MvNormal(μ, Symmetric(I + A' * A)),
        )
            # unscaled transformation
            for t in (
                ScaleTransform(1.0),
                ARDTransform(ones(dim)),
                LinearTransform(Diagonal(ones(dim))),
            )
                out = CalibrationErrorsDistributions.apply(t, d)
                @test mean(out) ≈ mean(d)
                @test cov(out) ≈ cov(d)
            end

            # scaling
            scale = rand()
            d_scaled = diagm(fill(scale, dim)) * d
            for t in (
                ScaleTransform(scale),
                ARDTransform(fill(scale, dim)),
                LinearTransform(Diagonal(fill(scale, dim))),
            )
                out = CalibrationErrorsDistributions.apply(t, d)
                @test mean(out) ≈ mean(d_scaled)
                @test cov(out) ≈ cov(d_scaled)
            end
        end
    end

    @testset "scale_cov" begin
        dim = 10
        v = rand(dim)
        γ = rand()
        X = diagm(γ .* v .^ 2)

        for A in (ScalMat(dim, γ), PDiagMat(fill(γ, dim)), PDMat(diagm(fill(γ, dim))))
            Y = CalibrationErrorsDistributions.scale_cov(A, v)
            @test Matrix(Y) ≈ X
        end
    end

    @testset "invquad_diff" begin
        dim = 10
        x = rand(dim)
        y = rand(dim)
        γ = rand()
        u = sum(abs2, x - y) / γ

        for A in (ScalMat(dim, γ), PDiagMat(fill(γ, dim)), PDMat(diagm(fill(γ, dim))))
            v = CalibrationErrorsDistributions.invquad_diff(A, x, y)
            @test v ≈ u
        end
    end
end


text:
@testset "BinaryHeap" begin
    data = [16, 2, 6, 5, 7, 1, 3, 4, 11]
    @testset "MinHeap" begin
        @testset "init heap/buildheap!" begin
            h = MinHeap(data)
            @test peek(h) == 1
            @test length(h) == 9
            @test !isempty(h)
            @test isequal(h.tree.data, [1, 2, 3, 4, 7, 6, 16, 5, 11])
        end
        @testset "insert! and remove!" begin
            h = MinHeap(Vector{Int}())
            @test_throws ArgumentError peek(h)
            @test length(h) == 0
            @test isempty(h)
            expected = [
                 [16],
                 [2, 16],
                 [2, 16, 6],
                 [2, 5, 6, 16],
                 [2, 5, 6, 16, 7],
                 [1, 5, 2, 16, 7, 6],
                 [1, 5, 2, 16, 7, 6, 3],
                 [1, 4, 2, 5, 7, 6, 3, 16],
                 [1, 4, 2, 5, 7, 6, 3, 16, 11]]
            for i = 1:length(data)
                insert!(h, data[i])
                @test length(h) == i
                @test !isempty(h)
                @test isequal(h.tree.data, expected[i])
            end
            removed = [1, 2, 3, 4, 5, 6, 7, 11, 16]
            while !isempty(h)
                @test remove!(h) == popfirst!(removed)
            end
        end
        @testset "remove! at position" begin
            h = MinHeap(data)
            @test_throws BoundsError(h, 10) remove!(h, 10)
            @test remove!(h, length(h)) == 11
            @test remove!(h, 3) == 3
            @test length(h) == 7
        end
    end
    @testset "MaxHeap" begin
        @testset "init heap/buildheap!" begin
            h = MaxHeap(data)
            @test peek(h) == 16
            @test length(h) == 9
            @test !isempty(h)
            @test isequal(h.tree.data, [16, 11, 6, 5, 7, 1, 3, 4, 2])
        end
        @testset "insert! and remove!" begin
            h = MaxHeap(Vector{Int}())
            @test_throws ArgumentError peek(h)
            @test length(h) == 0
            @test isempty(h)
            expected = [
                 [16],
                 [16, 2],
                 [16, 2, 6],
                 [16, 5, 6, 2],
                 [16, 7, 6, 2, 5],
                 [16, 7, 6, 2, 5, 1],
                 [16, 7, 6, 2, 5, 1, 3],
                 [16, 7, 6, 4, 5, 1, 3, 2],
                 [16, 11, 6, 7, 5, 1, 3, 2, 4]]
            for i = 1:length(data)
                insert!(h, data[i])
                @test length(h) == i
                @test !isempty(h)
                @test isequal(h.tree.data, expected[i])
            end
            removed = [1, 2, 3, 4, 5, 6, 7, 11, 16]
            while !isempty(h)
                @test remove!(h) == pop!(removed)
            end
        end
        @testset "remove! at position" begin
            h = MaxHeap(data)
            @test_throws BoundsError(h, 10) remove!(h, 10)
            @test remove!(h, length(h)) == 2
            @test remove!(h, 2) == 11
            @test length(h) == 7
        end
    end
end


text:
using GraphRecipes, Printf

export get_cooccurrence_matrix, get_cooccurrence_prob_matrix, plot_cooccurrence_graph
export get_informed_proposal_matrix

struct CondProbMatrix
    data::Matrix{Float64}
    function CondProbMatrix(data::Matrix{Float64})
        @assert prod(sum(data, dims=1) .≈ 1.0) "Input not valid. Must have column sums equal to one."
        new(data)
    end
end

struct CumCondProbMatrix
    data::Matrix{Float64}
    function CumCondProbMatrix(data::Matrix{Float64})
        @assert size(data)[1] == (size(data)[2]+1) "Incorrect dimensions. Must have columns as cumulative probabilities, with 0 first entry."
        @assert prod(data[1,:] .== 0.0) "Incorrect initial entries. First entry of each column must be 0.0"
        @assert prod(data[end,:] .≈ 1.0) "Incorrect final entries. Final entry of each column must be 1.0"
        new(data)
    end
end

Base.size(P::CondProbMatrix) = size(P.data)
Base.getindex(P::CondProbMatrix, ind1, ind2) = P.data[ind1, ind2]
# Base.show(io::IO, P::CondProbMatrix) = print(io, P.data)

function Base.show(io::IO, P::CondProbMatrix)
    (n,m) = size(P)
    println(io, "$(n)x$(m) Conditional Probability Matrix (Colwise)")
    for i in 1:n
        for j in 1:m
            print(io, @sprintf("%.2f ", P[i,j]))
        end
        print(io, "\n")
    end
end

Base.size(P::CumCondProbMatrix) = size(P.data)
Base.getindex(P::CumCondProbMatrix, ind1, ind2) = P.data[ind1, ind2]
function Base.show(io::IO, P::CumCondProbMatrix)
    (n,m) = size(P)
    println(io, "$(n)x$(m) Cumulative Conditional Probability Matrix (Colwise)")
    for i in 1:n
        for j in 1:m
            print(io, @sprintf("%.2f ", P[i,j]))
        end
        print(io, "\n")
    end
end
function CumCondProbMatrix(P::CondProbMatrix)
    P_cusum = cumsum(P.data, dims=1)
    P_cusum = [zeros(1, size(P)[1]); P_cusum]
    return CumCondProbMatrix(P_cusum)
end

function Base.cumsum(P::CondProbMatrix)
    P_cusum = cumsum(P.data, dims=1)
    P_cusum = [zeros(1, size(P)[1]); P_cusum]
    return CumCondProbMatrix(P_cusum)
end

function get_cooccurrence_matrix(posterior::Union{SisPosterior,SimPosterior})
    V = length(posterior.V)
    C = zeros(Int, V, V)
    for S in posterior.data
        C_tmp = zeros(Int, V, V)
        for I in S
            tmp_counts = counts(I, 1:V)
            vals = unique(I)
            # @show vals
            for j in 1:length(vals)
                for i in 1:(j-1)
                    C_tmp[vals[i],vals[j]] += 1
                    C_tmp[vals[j],vals[i]] += 1
                end
                if tmp_counts[vals[j]] > 1 # If vertex appears more than once add to it's diag
                    C_tmp[vals[j], vals[j]] += 1
                end
            end
        end
        C += Int.(C_tmp .> 0)
    end
    d = Dict(i => i for i in 1:V)
    return C, d, d
end

# function get_cooccurrence_matrix(posterior::Union{SisPosterior{String},SimPosterior{String}})

#     vertex_set = posterior.V
#     V = length(vertex_set)
#     C = zeros(Int, V, V)
#     ind_map = Dict{String, Int}(v => i for (i,v) in enumerate(vertex_set))
#     for S in posterior.data
#         C_tmp = zeros(Int, V, V)
#         for I in S
#             tmp_counts = countmap(I)
#             vals = unique(I)
#             for j in 1:length(vals)
#                 for i in 1:(j-1)
#                     C_tmp[ind_map[vals[i]],ind_map[vals[j]]] += 1 # only filling upper tri
#                     C_tmp[ind_map[vals[j]],ind_map[vals[i]]] += 1 # only filling upper tri
#                 end
#                 if tmp_counts[vals[j]] > 1 # If vertex appears more than once add to it's diag
#                     C_tmp[ind_map[vals[j]], ind_map[vals[j]]] += 1
#                 end
#             end
#         end
#         C += Int.(C_tmp .> 0)
#     end
#     ind_map_inv = Dict(val => key for (val, key) in zip(values(ind_map), keys(ind_map)))
#     return C, ind_map, ind_map_inv
# end

function get_cooccurrence_prob_matrix(C::Matrix{Int}, α::Real)

    @assert issymmetric(C) "Input matrix should be symmetric."

    V = size(C)[1]
    P = convert(Array{Float64,2}, C)

    for i in 1:V
        Z = sum(P[:,i])
        if  Z == 0.0
            P[:,i] = fill(1/V, V)
        else
            P[:,i] = ((P[:,i] ./ Z) .+ α) ./ (1 + V*α)
        end
    end
    return CondProbMatrix(P)
end

function get_informed_proposal_matrix(posterior::Union{SisPosterior,SimPosterior}, α::Real)
    C, vmap, vmap_inv = get_cooccurrence_matrix(posterior)
    P = get_cooccurrence_prob_matrix(C, α)
    return cumsum(P), vmap, vmap_inv
end

function plot_cooccurrence_graph(posterior::Union{SisPosterior,SimPosterior})
    C, vmap, vmap_inv = get_cooccurrence_matrix(posterior)
    C₀ = copy(C)
    C₀[diagind(C₀)] .= 0
    fig = graphplot(
        C₀, names=vmap,
        edge_width= (s,d,w) -> 2*C[s,d]/maximum(C),
        markersize=0.15, shape=:circle,
        fontsize=10,
        msw=0
    )
    display(fig)
end


text:
module TestEnsembleKombine

using Ensemble
using LinearAlgebra

using Statistics

using Test: @test, @testset

function testgaussian()
    ndim = 10
    mu = randn(ndim)
    sigma = randn(ndim, ndim)
    sigma = sigma'*sigma

    sigmafact = cholesky(sigma)

    function logpost(x)
        dx = x .- mu
        -0.5*sum(dx .* (sigmafact \ dx))
    end

    pts = randn(10, 128)

    prop = EnsembleKombine.build_proposal_kde(pts)

    lnprobs = [logpost(pts[:,j,1]) for j in 1:size(pts,2)]
    lnprops = [EnsembleKombine.logpdf(prop, pts[:,j]) for j in 1:size(pts,2)]

    # Burnin
    for j in 1:10
        pts, lnprobs, lnprops = EnsembleKombine.run_mcmc(pts, lnprobs, lnprops, logpost, prop, 128)

        pts = pts[:,:,end]
        lnprobs = lnprobs[:,end]
        lnprops = lnprops[:,end]
        prop = EnsembleKombine.build_proposal_kde(pts)
    end

    pts, lnprobs, lnprops = EnsembleKombine.run_mcmc(pts, lnprobs, lnprops, logpost, prop, 1024)

    @testset "Gaussian distribution tests" begin
        for i in 1:ndim
            @test isapprox(mean(pts[i,:,:]), mu[i], atol=sqrt(sigma[i,i])/10.0)
        end
    end
end

function testall()
    @testset "EnsembleKombine tests" begin
        testgaussian()
    end
end

end


text:
module MKLSparse

using MKL_jll

function __init__()
    ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint, (Cint,), Base.USE_BLAS64 ? 1 : 0)
end

include(joinpath("BLAS", "BLAS.jl"))

end # module


text:
struct SensorType
    createSensor_func::Function
    domain::Tuple
end

function get_transmission(problem, β; kwargs...)

    create_sensors_func = get(kwargs, :create_sensors_func, _create_sphere_sensor)
    θₘₐₓ = problem.atoms.N < 50 ? deg2rad(45) : deg2rad(15)
    integral_domain = get(kwargs, :domain, ((0.0, 0.0), (θₘₐₓ, 2π)))

    sensor_setup = SensorType(create_sensors_func, integral_domain)

    scattering = get(kwargs, :scattering, :farField)
    I_scattered = _get_intensity_laser_plus_scattering(problem,β; scattering = scattering, sensor_setup = sensor_setup)
    I_laser =     _get_intensity_laser(problem; sensor_setup)
    T = I_scattered / I_laser
    return T
end

function _get_intensity_laser_plus_scattering(
    problem,
    β;
    scattering::Symbol,
    sensor_setup::SensorType,
)
    if scattering == :nearField
        scattering_func = _scattering_nearField
    elseif scattering == :farField
        scattering_func = _scattering_farField
    else
        @error "Invalid scattering Value: should be `:nearField` or `:farField` "
    end

    view_of_atoms = view(problem.atoms.r, :, :)

    _toCall_I_total(x) = _I_total(
        x,
        problem,
        view_of_atoms,
        β,
        scattering_func,
        sensor_setup.createSensor_func,
    )
    (int_scatt, _e) = hcubature(_toCall_I_total, sensor_setup.domain...)
    return int_scatt
end
function _I_total(
    x,
    problem::LinearOptics,
    v_r,
    β::Vector{ComplexF64},
    scattering_func::Function,
    getSensor::Function,
)

    sensor = getSensor(x, problem)

    return _get_intensity_over_sensor(
        problem.atoms.shape,
        problem.laser,
        v_r,
        sensor,
        β;
        scattering_func = scattering_func,
    )
end

function _get_intensity_laser(problem; sensor_setup::SensorType)
    integral_args = sensor_setup.domain
    myargs = (:problem => problem, :createSensor_func => sensor_setup.createSensor_func)

    (int_laser, _e) = hcubature(x -> _I_laser(x; myargs...), integral_args...)
    return int_laser
end
function _I_laser(x; kwargs...)
    problem = kwargs[:problem]
    getSensor = kwargs[:createSensor_func]
    sensor = getSensor(x, problem)

    return abs2(apply_laser_over_oneSensor(problem.laser, sensor))
end

@inline function _create_sphere_sensor(x, problem)
    θ, ϕ = x[1], x[2]
    new_R = how_far_is_FarField(problem.atoms)
    spherical_coordinate = [θ, ϕ, new_R]
    sensor = sph2cart(spherical_coordinate)
    return sensor
end
@inline function _create_plane_sensor(x, problem)
    x, y = x[1], x[2]
    z = size(problem.atoms)
    new_z = z + 0.5 * z
    sensor = [x, y, new_z]
    return sensor
end

#=
    The value of "farField_factor*size(atoms)" is necessary to avoid
    transmission above 100% for small number of particles - numerical erros
    goes to zero, and we get invalid results, if I allow the Far Field
    be to far, i induce such mistakes.

    The "5.1*(size(atoms)^2)" expression is based on math arguments.
    See the explanation inside benchmarks folder ('benchmarks/benchmark_farField.jl')
=#
@inline function how_far_is_FarField(atoms)
    if atoms.N < 50
        return farField_factor * size(atoms)
    else
        return 5.01 * (size(atoms)^2)
    end
end


text:
# utility functions to run a batch of jobs

# assign the young bulls to random herd
function test_mating!(young_stud, herds, birthyear, n_test_mating_per_bull; pmp=0.5, method="random", updateinb=false, prob::Union{Nothing,Vector{Float64}}=nothing)
   n_herds = length(herds)
   n_pregnant = 0
   if method=="random"
      for bull in young_stud.sires
         # n mating per young bull
         for i in 1:n_test_mating_per_bull
            h = rand(1:n_herds)
            preg_dams = mating!(young_stud, herds[h], "fgroup", year=birthyear, method="dairy_standard_ai", n=1, pmp=pmp, updateinb=updateinb, prob=prob);
            n_pregnant = n_pregnant + length(preg_dams)
         end
      end
   else
      throw(ArgumentError("use `random` for mating method"))
   end
   return n_pregnant
end

function regular_mating!(bull_stud, herds, birthyear; method="random", updateinb=false, prob::Union{Nothing,Vector{Float64}}=nothing)
   n_herds = length(herds)
   n_pregnant = 0
   if method=="random"
      for h in 1:n_herds
         preg_dams = mating!(bull_stud, herds[h], "fgroup", year=birthyear, method="dairy_standard_ai", updateinb=updateinb, prob=prob)
         n_pregnant = n_pregnant + length(preg_dams)
      end
   else
      throw(ArgumentError("use `random` for mating method"))
   end
   return n_pregnant
end

function cull_old_bulls!(bull_stud,birth_year_of_culled_bulls)
   return cull!(bull_stud, selectid([:male,:year] => (x,y) -> x==true && y==birth_year_of_culled_bulls, bull_stud))
end

function cull_old_cows!(herds,birth_year_of_culled_cows)
   n_herds = length(herds)
   nculled = 0
   for h=1:n_herds
      nculled = nculled + cull!(herds[h], selectid([:male,:year] => (x,y) -> x==false && y==birth_year_of_culled_cows, herds[h]))
   end
   return nculled
end

function mean_inbreeding(pop; year::Int=-1)
   if year>=0
      return mean(pop.df[pop.df[:,:year].==year, :inb])
   else
      return mean(pop.df[:, :inb])
   end
end

function number_of_active_females(herd::PTGroup)
   females = selectid([:male]=>x->x==false, herd, aliveonly=true)
   nheifers = length(setdiff(females, herd.dams))
   ndams = length(herd.dams)
   return ndams,nheifers
end

function number_of_active_females(herds::Vector{PTGroup})
   n_herds = length(herds)
   ndams = 0
   nheifers = 0
   for h=1:n_herds
      (nd,nh) = number_of_active_females(herds[h])
      nheifers = nheifers + nh
      ndams = ndams + nd
   end
   return ndams,nheifers
end

function phenotyped_cows(herds::Vector{PTGroup}; year=-1, test_status=0, plan="limited", n_breeders=0, n_daughters=-1, verbose=false)
   n_herds = length(herds)
   pop = herds[1].pop
   if plan=="limited"
      # rule: All tested cows have phenotypes.
      # rule: All cows in breeder herds have phenotypes.
      # rule: Each bull has the upper limit on the number of daughters each year.
      tested_id = zeros(Int,0)
      breeder_cow_id = zeros(Int,0)
      standard_cow_id = zeros(Int,0)
      phenotyped_cow_id = zeros(Int,0)
      if year>0
         # tested daughters
         tested_id = selectid([:male,:year,:status] => (x,y,z) -> x==false && y==year && z==test_status, pop)
         # regular daughters
         for i=1:n_breeders
            breeder_cow_id = [breeder_cow_id; selectid([:male,:year] => (x,y) -> x==false && y==year, herds[i])]
         end
         for i=n_breeders+1:n_herds
            standard_cow_id = [standard_cow_id; selectid([:male,:year] => (x,y) -> x==false && y==year, herds[i])]
         end
         n_tested = length(tested_id)
         n_breeder = length(breeder_cow_id)
         n_standard = length(standard_cow_id)
         if n_daughters<0
            n_additional = n_standard
         else
            n_additional = min(max(0,n_daughters-n_breeder),n_standard)
         end

         if verbose
            println("  phenotyped cows:")
            println("    tested:$(n_tested) breeder:$(n_breeder) standard:$(n_additional); total:$(n_daughters)")
         end

         # tested daughters + breeder cows + additional cows
         phenotyped_cow_id = [phenotyped_cow_id; tested_id; breeder_cow_id]
         if n_additional>0
            shuffle!(standard_cow_id)
            phenotyped_cow_id = [phenotyped_cow_id; standard_cow_id[1:n_additional]]
         end
         sort!(phenotyped_cow_id)
         return phenotyped_cow_id
      else
         throw(ArgumentError("year not given"))
      end
   else
      throw(ArgumentError("plan not supported: $(plan)"))
   end
   return zeros(Int,0)
end

# To get N daughters on average with SD
# possible use:
#   Gamma(0.8,1.0); scale=65, shift=50; mean=100 std=55 upper=200 (lower=50)
#   Gamma(1.0,1.0); scale=40, shift=50; mean=90  std=40 upper=250 (lower=50)
#   or Gamma(1.0,1.0) upper=4.0 (98% percentile)
function expected_frequency_of_mating(mgroup::PTGroup,dist; by=:ebv, rev=true, scale=1.0, shift=0.0, upper=0.0)
   pop = mgroup.pop
   n = length(mgroup.sires)
   benchmark = pop.df[mgroup.sires,by]
   perm = sortperm(benchmark,rev=rev) # large to small by default
   iperm = invperm(perm)
   x = sort((rand(dist,n)*scale .+ shift),rev=true)  # large to small
   if upper>0
      x[x .> upper] .= upper
   end
   x .= x[iperm]/sum(x)
   return x
end


text:
__precompile__()

"""
Never Fail Another Interview™
"""
module FizzBuzz

export fizzbuzz, @fizzbuzz

"""
    fizzbuzz(iterable::AbstractArray)
    fizzbuzz(iterable::AbstractArray, io:IO)

FizzBuzz-ify an `iterable`. Optionally print to `io` instead of `STDOUT`.

```jldoctest fb_env
julia> using FizzBuzz

julia> fizzbuzz(1:5)
1
2
fizz
4
buzz

julia> io = IOBuffer();

julia> fizzbuzz(1:15, io)

julia> String(take!(io))
"1\\n2\\nfizz\\n4\\nbuzz\\nfizz\\n7\\n8\\nfizz\\nbuzz\\n11\\nfizz\\n13\\n14\\nfizz buzz\\n"
```
"""
function fizzbuzz(iterable::T, io::Union{Base.TTY, IO}=STDOUT) where T<:AbstractArray{<:Real}
    for i = iterable
        println(io, fizzbuzzer(i))
    end
end

"""
Private functions that compute return values for elements of `iterable`.
"""
function fizzbuzzer(v::T)::Union{T, String} where T<:Real
    if v % 15 == 0
        return "fizz buzz"
    elseif v % 5 == 0
        return "buzz"
    elseif v % 3 == 0
        return "fizz"
    else
        return v
    end
end

function fizzbuzzer(v::T)::Union{T, String} where T<:AbstractFloat
    if v % 15. == 0.
        return "fizz buzz"
    elseif v % 5. == 0.
        return "buzz"
    elseif v % 3. == 0.
        return "fizz"
    else
        return v
    end
end

"""
    @fizzbuzz iterable::AbstractArray
    @fizzbuzz iterable::AbstractArray io::IO

Macro version of `fizzbuzz()`.

```jldoctest fb_env
julia> @fizzbuzz 1:5
1
2
fizz
4
buzz
```
"""
macro fizzbuzz(iterable::Expr)
    return quote
        fizzbuzz($iterable)
    end
end

end  # module


text:
"""
Clenshaw-Curtis quadrature.
"""
function ClenshawCurtisQuadrature(::Type{T}, s::Integer; IT=BigFloat) where {T}
    if s == 1
        throw(ErrorException("Clenshaw-Curtins quadrature is not defined for one stage."))
    end

    mymapreduce(f, op, itr) = isempty(itr) ? zero(IT) : mapreduce(f, op, itr)

    b(j,n) = fld(n,2) == j == cld(n,2) ? IT(1) : IT(2)
    c(k,n) = mod(k,n) == 0 ? IT(1) : IT(2)
    ϑ(k,n) = @big 2π * k / n

    cctrm(k,j,n) = @big b(j,n) / ( 4 * j^2 - 1 ) * cos(j * ϑ(k,n))
    ccsum(k,n) = c(k,n) / IT(2n) * ( IT(1) - mymapreduce(j -> cctrm(k,j,n), +, 1:div(n,2)) )

    x = chebyshev_points(IT, s, Val(2))
    w = [ccsum(i-1,s-1) for i in 1:s]

    QuadratureRule(s^2, x, w, T)
end

ClenshawCurtisQuadrature(s) = ClenshawCurtisQuadrature(Float64, s)


text:
# Adapted from:
# SOSDEMO7 --- Chebyshev polynomials
# Section 3.7 of SOSTOOLS User's Manual

@testset "SOSDEMO7 with $(factory.constructor)" for factory in sdp_factories
    isscs(factory) && continue
    ndeg = 8   # Degree of Chebyshev polynomial

    @polyvar x

    Z = monomials((x,), 0:ndeg-1)

    m = SOSModel(factory)

    @variable m γ
    @variable m p1 Poly(Z)

    p = p1 + γ * x^ndeg # the leading coeff of p is γ

    dom = @set x >= -1 && x <= 1
    @constraint(m, p <= 1, domain = dom)
    @constraint(m, p >= -1, domain = dom)

    @objective m Max γ

    JuMP.optimize!(m)

    @test JuMP.primal_status(m) == MOI.FeasiblePoint

    @test isapprox(JuMP.value(p), 128x^8 - 256x^6 + 160x^4 - 32x^2 + 1, ztol=1e-7, atol=1e-7)
    @test isapprox(JuMP.value(γ), 128)
end


text:
@testset "NormalDDP" begin
    rng = MersenneTwister(1)
    N, G, v0, r0, u0, s0 = 10, 4, 2.0, 3.0, 3.0, 9.0
    m = NormalDDP(rng, N, G; v0, r0, u0, s0)
    @test m.G  == 4
    @test m.v0 == 2
    @test m.r0 == 3
    @test m.u0 == 3.0
    @test m.s0 == 9.0
    @test m.v1 == [2 * ones(G)]
    @test m.r1 == [3 * ones(G)]
    @test m.u1 == [3 * ones(G)]
    @test m.s1 == [9 * ones(G)]
    @test m.gamma  == ones(Bool, G)
end

@testset "NormalDDP inherited accessors" begin
    rng = MersenneTwister(1)
    N, G, K0, a0, b0, v0, r0, u0, s0 = 10, 4, 5, 2.0, 4.0, 2.0, 3.0, 3.0, 9.0
    m = NormalDDP(rng, N, G; K0, a0, b0, v0, r0, u0, s0)
    @test dp_mass(m) > 0.0
    @test dp_mass(m) < Inf
    @test n_clusters(m) == K0
    @test sum(cluster_sizes(m) .== 0) == 1
    @test sum(cluster_sizes(m)) == N
    @test active_clusters(m) == Set(1:K0)
    @test passive_clusters(m) == Set(K0 + 1)
    @test cluster_capacity(m) == K0 + 1
end

@testset "add_cluster!" begin
    rng = MersenneTwister(1)
    N, G, K0, a0, b0, v0, r0, u0, s0 = 10, 4, 5, 2.0, 4.0, 2.0, 3.0, 3.0, 9.0
    m = NormalDDP(rng, N, G; K0, a0, b0, v0, r0, u0, s0)
    ANOVADDPTest.add_cluster!(m)
    @test length(m.v1) == 2
    @test length(m.r1) == 2
    @test length(m.u1) == 2
    @test length(m.s1) == 2
end

@testset "update_suffstats! (1)" begin
    rng = MersenneTwister(1)
    data = NormalData([1], [1.0])
    N, G, K0, v0, r0, u0, s0 = 1, 1, 1, 1.0, 1.0, 0.0, 1.0
    m = NormalDDP(rng, N, G; K0, v0, r0, u0, s0)
    ANOVADDPTest.update_suffstats!(m, data)
    @test m.v1[1][1] ≈ 2.0
    @test m.v1[2][1] ≈ 1.0
    @test m.r1[1][1] ≈ 2.0
    @test m.r1[2][1] ≈ 1.0
    @test m.u1[1][1] ≈ 0.5
    @test m.u1[2][1] ≈ 0.0
    @test m.s1[1][1] ≈ 1.5
    @test m.s1[2][1] ≈ 1.0
end

@testset "update_suffstats! (2)" begin
    rng = MersenneTwister(1)
    data = NormalData([1], [1.0])
    N, G, K0, v0, r0, u0, s0 = 1, 1, 1, 1.0, 1.0, 0.0, 1.0
    m = NormalDDP(rng, N, G; K0, v0, r0, u0, s0)
    ANOVADDPTest.update_suffstats!(m, data)
    ANOVADDPTest.update_suffstats!(m, data, 1, 1, 2)
    @test m.v1[1][1] ≈ 1.0
    @test m.v1[2][1] ≈ 2.0
    @test m.r1[1][1] ≈ 1.0
    @test m.r1[2][1] ≈ 2.0
    @test m.u1[1][1] ≈ 0.0
    @test m.u1[2][1] ≈ 0.5
    @test m.s1[1][1] ≈ 1.0
    @test m.s1[2][1] ≈ 1.5
end

@testset "logpredlik (empty clusters)" begin
    rng = MersenneTwister(1)
    data = NormalData([1], [1.0])
    N, G, K0, v0, r0, u0, s0 = 1, 1, 1, 1.0, 1.0, 0.0, 1.0
    m = NormalDDP(rng, N, G; K0, v0, r0, u0, s0)
    ANOVADDPTest.update_suffstats!(m, data)
    @test ANOVADDPTest.logpredlik(m, data, 1, first(passive_clusters(m))) ≈ (
        0.5 * 1.0 * log(1.0) -
        0.5 * 2.0 * log(1.5) +
        loggamma(2.0 / 2) -
        loggamma(1.0 / 2) +
        0.5 * log(1.0 / 2.0) -
        0.5 * log(π)
    )
end

@testset "logpredlik (non-empty clusters)" begin
    rng = MersenneTwister(1)
    data = NormalData([1, 1], [1.0, 0.0])
    N, G, K0, v0, r0, u0, s0 = 2, 1, 1, 1.0, 1.0, 0.0, 1.0
    m = NormalDDP(rng, N, G; K0, v0, r0, u0, s0)
    ANOVADDPTest.update_suffstats!(m, data)
    @test m.v1[1][1] ≈ 3.0
    @test m.r1[1][1] ≈ 3.0
    @test m.u1[1][1] ≈ 1/3
    @test m.s1[1][1] ≈ 5/3

    @test ANOVADDPTest.logpredlik(m, data, 2, 1) ≈ (
        0.5 * 2 * log(1.5) -
        0.5 * 3 * log(5/3) +
        loggamma(3/2) -
        loggamma(2/2) +
        0.5 * log(2/3) -
        0.5 * log(π)
    )
end

@testset "update! (1)" begin
    rng = MersenneTwister(1)
    data = NormalData([1, 1], [1.0, 0.0])
    N, G, K0, v0, r0, u0, s0 = 2, 1, 1, 1.0, 1.0, 0.0, 1.0
    m = NormalDDP(rng, N, G; K0, v0, r0, u0, s0)
    update!(rng, m, data)
end

@testset "update! (2)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = randn(rng, N)
    data = NormalData(x, y)
    m = NormalDDP(rng, N, G; K0)
    update!(rng, m, data)
end

@testset "final_example" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = randn(rng, N)
    data = NormalData(x, y)
    m = NormalDDP(rng, N, G; K0)
    for t in 1:10
        update!(rng, m, data)
    end
end

@testset "train (1)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = randn(rng, N)
    m = NormalDDP(rng, N, G; K0)
    data = NormalData(x, y)
    predict =
        expandgrid(1:G, range(-2.0, stop = 2.0, length = 50)) |>
        x -> NormalData(x...)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, false, false])
end

@testset "train (2)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = randn(rng, N)
    for i = 1:N
        if x[i] == 2
            y[i] += 0.5
        end
    end
    data = NormalData(x, y)
    predict =
        expandgrid(1:G, range(-2.0, stop = 2.0, length = 50)) |>
        x -> NormalData(x...)
    m = NormalDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, true, false])
end

@testset "train (3)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = randn(rng, N)
    for i = 1:N
        if x[i] == 3
            y[i] += 0.5
        end
    end
    data = NormalData(x, y)
    predict =
        expandgrid(1:G, range(-2.0, stop = 2.0, length = 50)) |>
        x -> NormalData(x...)
    m = NormalDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, false, true])
end

@testset "train (4)" begin
    rng = MersenneTwister(1)
    N, G, K0 = 1000, 3, 1
    x = rand(rng, 1:G, N)
    y = randn(rng, N)
    for i = 1:N
        if x[i] != 1
            y[i] += 0.5
        end
    end
    data = NormalData(x, y)
    predict =
        expandgrid(1:G, range(-2.0, stop = 2.0, length = 50)) |>
        x -> NormalData(x...)
    m = NormalDDP(rng, N, G; K0)
    chain = train(rng, m, data, predict)
    @test all(mode(chain.gamma) .== [true, true, true])
end

# # Note:
# # Remember that you can use denserank(x) for converting
# # vector codes to numeric codes


text:
### Vector/Collections API

function vcat(x::RLEVector, y::RLEVector)
  RLEVector( vcat(x.runvalues, y.runvalues), vcat(x.runends, y.runends + length(x)))
end

function pop!(x::RLEVector)
  runcount = nrun(x)
  isempty(x) && throw(ArgumentError("array must be non-empty"))
  item = x.runvalues[end]
  x.runends[end] -= 1
  if x.runends[end] == 0
    deleteat!(x.runvalues,runcount)
    deleteat!(x.runends,runcount)
  end
  return(item)
end

function push!{T,T2 <: Integer}(x::RLEVector{T,T2},item)
  item = convert(T,item) # Copying how base does it for arrays
  if !isempty(x) && (item == x.runvalues[end])
    x.runends[end] += 1
  else
    push!(x.runvalues,item)
    push!(x.runends, (isempty(x) ? 0 : x.runends[end]) + 1)
  end
  return(x)
end

function shift!(x::RLEVector)
  isempty(x) && throw(ArgumentError("array must be non-empty"))
  item = x.runvalues[1]
  x.runends -= 1
  if x.runends[1] == 0
    deleteat!(x.runvalues,1)
    deleteat!(x.runends,1)
  end
  return(item)
end

function shove!{T,T2 <: Integer}(x::RLEVector{T,T2},item)
  item = convert(T,item) # Copying how base does it for arrays
  x.runends += 1
  if item != x.runvalues[1]
    unshift!(x.runvalues,item)
    unshift!(x.runends,1)
  end
  return(x)
end
unshift!{T,T2 <: Integer}(x::RLEVector{T,T2},item) = shove(x,item) # Does unshift come from perl? Isn't Larry Wall a linguist? C'mon!

function insert!{T,T2 <: Integer}(x::RLEVector{T,T2},i::Integer,item)
  if i == length(x) + 1
    push!(x,item)
    return(x)
  end
  if i < 1 || i > length(x)
    throw(BoundsError())
  end
  (run, index_in_run, run_remainder) = ind2runcontext(x,i)
  if item == x.runvalues[run]
    x.runends[run:end] += 1
  else
    if index_in_run == 1
      ins_vals = [item; x.runvalues[run]]
      ins_ends = [i; x.runends[run]+1]
    else
      ins_vals = [x.runvalues[run]; item; x.runvalues[run]]
      ins_ends = [i-1; i; x.runends[run]+1]
    end
    x.runvalues = vcat( x.runvalues[1:run-1], ins_vals, x.runvalues[run+1:end] )
    x.runends = vcat( x.runends[1:run-1], ins_ends, x.runends[run+1:end] + 1 )
  end
  x.runvalues,x.runends = ree(x.runvalues,x.runends)
  return(x)
end

function deleterun!(x::RLEVector,i::Integer)
  x.runends[i:end] -= rwidth(x,i)
  if (i > 1 && i < nrun(x) && x.runvalues[i-1] == x.runvalues[i+1])
    splice!(x.runvalues,(i-1):i)
    splice!(x.runends,(i-1):i)
  else
    deleteat!(x.runvalues,i)
    deleteat!(x.runends,i)
  end
  return(x)
end

function decrement_run!(x::RLEVector,run::Integer)
  if rwidth(x,run) == 1
    deleterun!(x,run)
  else
    x.runends[run:end] -= 1
  end
  return(x)
end

function deleteat!(x::RLEVector,i::Integer)
  run = ind2run(x,i)
  decrement_run!(x,run)
end

_default_splice = RLEVector(Union{}[],Int64[])
function splice!(x::RLEVector, i::Integer, ins::RLEVector=_default_splice)
  if i < 1 || i > length(x)
    throw(BoundsError())
  end
  if length(ins) == 0
    run = ind2run(x,i)
    current = x.runvalues[run]
    decrement_run!(x,run)
  else
    (run, index_in_run, run_remainder) = ind2runcontext(x,i)
    current = x.runvalues[run]
    right_shift = length(ins) - length(i)
    x.runends[run:end] += right_shift
    ins.runends += (i-1)
    if index_in_run == 1
      ins_vals = [ins.runvalues; x.runvalues[run]]
      ins_ends = [ins.runends; x.runends[run]]
    else
      ins_vals = [x.runvalues[run]; ins.runvalues; x.runvalues[run]]
      ins_ends = [i-1; ins.runends; x.runends[run]]
    end
    x.runvalues = vcat( x.runvalues[1:run-1], ins_vals, x.runvalues[run+1:end] )
    x.runends = vcat( x.runends[1:run-1], ins_ends, x.runends[run+1:end])
    x.runvalues,x.runends = ree(x.runvalues,x.runends)
  end
  return(current)
end

function splice!(x::RLEVector, index::Range, ins::RLEVector=_default_splice) # Can I do index::Union(Integer,UnitRange) here to have just one method?
  i_left = start(index)
  i_right = last(index)
  if i_left < 1 || i_right > length(x)
    throw(BoundsError())
  end
  if length(index) == 0
  current = similar(x,0)
  (run_right, index_in_run_right, run_remainder_right) = (run_left, index_in_run_left, run_remainder_left) = ind2runcontext(x,i_left)
  else
    current = x[index]
    (run_left, index_in_run_left, run_remainder_left) = ind2runcontext(x,i_left)
    (run_right, index_in_run_right, run_remainder_right) = ind2runcontext(x,i_right)
  end
  ins.runends += (i_left - 1)
  right_shift = nrun(ins) - length(index)
  x.runends[run_right:end] += right_shift
  if index_in_run_left == 1
    ins_vals = [ins.runvalues; x.runvalues[run_right]]
    ins_ends = [ins.runends; x.runends[run_right]]
  else
    ins_vals = [x.runvalues[run_left]; ins.runvalues; x.runvalues[run_right]]
    ins_ends = [i_left-1; ins.runends; x.runends[run_right]]
  end
  x.runvalues = vcat( x.runvalues[1:run_left-1], ins_vals, x.runvalues[run_right+1:end] )
  x.runends = vcat( x.runends[1:run_left-1], ins_ends, x.runends[run_right+1:end] )
  x.runvalues, x.runends = ree(x.runvalues,x.runends)
  return(current)
end

function splice!(x::RLEVector, i::Integer, ins::AbstractArray)
  splice!(x,i,RLEVector(ins))
end

function splice!(x::RLEVector, i::Range, ins::AbstractArray)
  splice!(x,i,RLEVector(ins))
end

# Appended space initialized with zero unlike base array
function resize!(x::RLEVector, nl::Integer) # Based on base version for array
  l = length(x)
  if nl > l
    push!(x.runends,nl)
    push!(x.runvalues,0)
  else
    nl < 0 && throw(ArgumentError("new length must be ≥ 0"))
    (run, index_in_run, run_remainder) = ind2runcontext(x,nl)
    x.runends = x.runends[1:run]
    x.runends[end] = x.runends[end] - run_remainder
    x.runvalues = x.runvalues[1:run]
  end
  return(x)
end

function empty!(x::RLEVector)
  empty!(x.runvalues)
  empty!(x.runends)
  return(x)
end

### Set API
# Mostly handled by convert and promote rules in RLEVector-type.jl

function intersect(x::RLEVector, sets...)
    ok = trues(length(x.runvalues))
    for (i, v) in enumerate(x.runvalues)
        for s in sets
            if !in(v, s)
                ok[i] = false
                break
            end
        end
    end
    RLEVector( x.runvalues[ ok ], cumsum( rwidth(x)[ ok ] ) )
end


text:
"""
    FEMMDeforLinearModule

Module for operations on interiors of domains to construct system matrices and
system vectors for linear deformation models.
"""
module FEMMDeforLinearModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: FESet, manifdim
import FinEtools.IntegDataModule: IntegData
import FinEtools.FEMMDeforLinearBaseModule: FEMMDeforLinearAbstract
import FinEtools.DeforModelRedModule: DeforModelRed, DeforModelRed2DAxisymm
import FinEtools.MatDeforModule: MatDefor
import FinEtools.CSysModule: CSys

"""
    FEMMDeforLinear{S<:FESet, F<:Function, P<:PropertyDeformationLinear}

Class for linear deformation finite element modeling machine.
"""
mutable struct FEMMDeforLinear{MR<:DeforModelRed,  S<:FESet, F<:Function, M<:MatDefor} <: FEMMDeforLinearAbstract
    mr::Type{MR}
    integdata::IntegData{S, F} # geometry data
    mcsys::CSys # updater of the material orientation matrix
    material::M # material object
end

function FEMMDeforLinear(mr::Type{MR}, integdata::IntegData{S, F}, material::M) where {MR<:DeforModelRed, S<:FESet, F<:Function, M<:MatDefor}
    # @show mr
    # @show material.mr
    @assert mr === material.mr "Model reduction is mismatched"
    @assert (integdata.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    return FEMMDeforLinear(mr, integdata, CSys(manifdim(integdata.fes)), material)
end

function FEMMDeforLinear(mr::Type{MR}, integdata::IntegData{S, F}, material::M, print) where {MR<:DeforModelRed, S<:FESet, F<:Function, M<:MatDefor}
    @show mr
    @show material.mr
    @assert mr === material.mr "Model reduction is mismatched"
    @assert (integdata.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    return FEMMDeforLinear(mr, integdata, CSys(manifdim(integdata.fes)), material)
end

end


text:
using MRIsim
using Test

@testset "Spinors×Magnetization" begin
    # Spinor 2x2 representation should be equivalent to a 3x1 vector rotation
    x = rand(3); x = x./sum(x)
    θ = rand() * π
    n = rand(3); n = n./sqrt(sum(n.^2))
    z = Mag(x[1]+1im*x[2], x[3])
    # General rotation
    xx1 = Q(θ,n[1]+1im*n[2],n[3])*z; #Spinor rot Q.(φ, B1./B, Bz./B)
    xx2 = Un(θ,n)*x; #3D rot matrix
    xx1 = [real(xx1.xy), imag(xx1.xy), xx1.z]
    @test xx1 ≈ xx2
    # Rot x
    nx = [1,0,0]
    xx1 = Rx(θ)*z; #Spinor rot
    xx2 = Un(θ,nx)*x; #3D rot matrix
    xx1 = [real(xx1.xy), imag(xx1.xy), xx1.z]
    @test xx1 ≈ xx2
    # Rot y
    nx = [0,1,0]
    xx1 = Ry(θ)*z; #Spinor rot
    xx2 = Un(θ,nx)*x; #3D rot matrix
    xx1 = [real(xx1.xy), imag(xx1.xy), xx1.z]
    @test xx1 ≈ xx2
    # Rot z
    nx = [0,0,1]
    xx1 = Rz(θ)*z; #Spinor rot
    xx2 = Un(θ,nx)*x; #3D rot matrix
    xx1 = [real(xx1.xy), imag(xx1.xy), xx1.z]
    @test xx1 ≈ xx2
end

@testset "Sequence" begin
    # Rotation 2D case
    A1, A2, T, t = rand(4)
    s = Sequence([Grad(A1,T);
                  Grad(A2,T)])
    θ = π*t
    R = MRIsim.rotz(θ)
    s2 = R*s
    GR2 = R*s.GR.A
    @test s2.GR.A ≈ GR2
    # Rotation 3D case
    T, t1, t2, t3 = rand(4)
    N = 100
    GR = [Grad(rand(),T) for i=1:3, j=1:N]
    s = Sequence(GR)
    α, β, γ = π*t1, π*t2, π*t3
    Rx = MRIsim.rotx(α)
    Ry = MRIsim.roty(β)
    Rz = MRIsim.rotz(γ)
    R = Rx*Ry*Rz
    s2 = R*s
    GR2 = R*GR.A
    @test s2.GR.A ≈ GR2
    # Concatenation
    A1, A2, A3, T1 = rand(4)
    s1 = Sequence([Grad(A1,T1);
                   Grad(A2,T1)],
                    [RF(A3,T1)])
    B1, B2, B3, T2 = rand(4)
    s2 = Sequence([Grad(B1,T2);
                   Grad(B2,T2)],
                    [RF(B3,T2)])
    s = s1 + s2
    @test s.GR.A ≈ [s1.GR.A s2.GR.A]
    @test s.RF.A ≈ [s1.RF.A s2.RF.A]
    @test s.DAC.N ≈ [s1.DAC.N ; s2.DAC.N]
end


text:
# ---
# title: 402. Remove K Digits
# id: problem402
# author: AquaIndigo
# date: 2020-11-15
# difficulty: Medium
# categories: Stack, Greedy
# link: <https://leetcode.com/problems/remove-k-digits/description/>
# hidden: true
# ---
#
# Given a non-negative integer _num_ represented as a string, remove _k_ digits
# from the number so that the new number is the smallest possible.
#
# **Note:**
#
#   * The length of _num_ is less than 10002 and will be ≥ _k_.
#   * The given _num_ does not contain any leading zero.
#
# **
#
# **Example 1:**
#
#
#
#     Input: num = "1432219", k = 3
#     Output: "1219"
#     Explanation: Remove the three digits 4, 3, and 2 to form the new number 1219 which is the smallest.
#
#
# **Example 2:**
#
#
#
#     Input: num = "10200", k = 1
#     Output: "200"
#     Explanation: Remove the leading 1 and the number is 200. Note that the output must not contain leading zeroes.
#
#
# **Example 3:**
#
#
#
#     Input: num = "10", k = 2
#     Output: "0"
#     Explanation: Remove all the digits from the number and it is left with nothing which is 0.
#
#
#
## @lc code=start
using LeetCode

function remove_kdigits(num::String, k::Int)::String
    num_code = codeunits(num)
    stk = UInt8[]
    for ch in num_code
        if k == 0 || isempty(stk)
            push!(stk, ch)
        else
            while !isempty(stk) && stk[end] > ch && k > 0
                pop!(stk)
                k -= 1
            end
            push!(stk, ch)
        end
    end
    idx = findfirst(x -> x != UInt8('0'), stk)
    if idx === nothing
        return "0"
    end
    return String(@view stk[idx:(end - k)])
end
## @lc code=end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# fast math

@testset "check expansions" begin
    @test macroexpand(Main, :(@fastmath 1+2)) == :(Base.FastMath.add_fast(1,2))
    @test macroexpand(Main, :(@fastmath +)) == :(Base.FastMath.add_fast)
    @test macroexpand(Main, :(@fastmath min(1))) == :(Base.FastMath.min_fast(1))
    @test macroexpand(Main, :(@fastmath min)) == :(Base.FastMath.min_fast)
    @test macroexpand(Main, :(@fastmath x.min)) == :(x.min)
    @test macroexpand(Main, :(@fastmath sincos(x))) == :(Base.FastMath.sincos_fast(x))
end
const one32 = one(Float32)
const eps32 = eps(Float32)
const eps32_2 = eps32/2
# Note: Cannot use local functions since these are not yet optimized
fm_ieee_32(x) = x + eps32_2 + eps32_2
fm_fast_32(x) = @fastmath x + eps32_2 + eps32_2
const one64 = one(Float64)
const eps64 = eps(Float64)
const eps64_2 = eps64/2
# Note: Cannot use local functions since these are not yet optimized
fm_ieee_64(x) = x + eps64_2 + eps64_2
fm_fast_64(x) = @fastmath x + eps64_2 + eps64_2
fm_ieee_64_upd(x) = (r=x; r+=eps64_2; r+=eps64_2)
fm_fast_64_upd(x) = @fastmath (r=x; r+=eps64_2; r+=eps64_2)
@testset "basic arithmetic" begin
    @test fm_ieee_32(one32) == one32
    @test (fm_fast_32(one32) == one32 ||
           fm_fast_32(one32) == one32 + eps32 > one32)
    @test fm_ieee_64(one64) == one64
    @test (fm_fast_64(one64) == one64 ||
           fm_fast_64(one64) == one64 + eps64 > one64)
    # check updating operators
    @test fm_ieee_64_upd(one64) == one64
    @test (fm_fast_64_upd(one64) == one64 ||
           fm_fast_64_upd(one64) == one64 + eps64 > one64)

    let epsf = 1.0f0/2^15, one_epsf = 1+epsf
        @test @fastmath(one_epsf * one_epsf - 1) ≈ Float32(65537/1073741824)
    end
    let eps = 1.0/2^30, one_eps = 1+eps
        @test @fastmath(one_eps * one_eps - 1) ≈ 2147483649/1152921504606846976
    end

    for T in (Float32, Float64, BigFloat)
        zero = convert(T, 0)
        one = convert(T, 1) + eps(T)
        two = convert(T, 2) + 1//10
        three = convert(T, 3) + 1//100

        @test @fastmath(+two) ≈ +two
        @test @fastmath(-two) ≈ -two
        @test @fastmath(zero+one+two) ≈ zero+one+two
        @test @fastmath(zero-one-two) ≈ zero-one-two
        @test @fastmath(one*two*three) ≈ one*two*three
        @test @fastmath(one/two/three) ≈ one/two/three
        @test @fastmath(rem(two,three)) ≈ rem(two,three)
        @test @fastmath(mod(two,three)) ≈ mod(two,three)
        @test @fastmath(cmp(two,two)) == cmp(two,two)
        @test @fastmath(cmp(two,three)) == cmp(two,three)
        @test @fastmath(cmp(three,two)) == cmp(three,two)
        @test @fastmath(one/zero) == convert(T,Inf)
        @test @fastmath(-one/zero) == -convert(T,Inf)
        @test isnan(@fastmath(zero/zero)) # must not throw

        for x in (zero, two, convert(T, Inf), convert(T, NaN))
            @test @fastmath(isfinite(x))
            @test !@fastmath(isinf(x))
            @test !@fastmath(isnan(x))
            @test !@fastmath(issubnormal(x))
        end
    end

    for T in (ComplexF32, ComplexF64, Complex{BigFloat})
        zero = convert(T,0)
        one = convert(T,1) + im*eps(real(convert(T,1)))
        two = convert(T,2) + im//10
        three = convert(T,3) + im//100

        @test @fastmath(+two) ≈ +two
        @test @fastmath(-two) ≈ -two
        @test @fastmath(zero+one+two) ≈ zero+one+two
        @test @fastmath(zero-one-two) ≈ zero-one-two
        @test @fastmath(one*two*three) ≈ one*two*three
        @test @fastmath(one/two/three) ≈ one/two/three
        @test @fastmath(three == two) == (three == two)
        @test @fastmath(three != two) == (three != two)
        @test isnan(@fastmath(one/zero))  # must not throw
        @test isnan(@fastmath(-one/zero)) # must not throw
        @test isnan(@fastmath(zero/zero)) # must not throw

        for x in (zero, two, convert(T, Inf), convert(T, NaN))
            @test @fastmath(isfinite(x))
            @test !@fastmath(isinf(x))
            @test !@fastmath(isnan(x))
            @test !@fastmath(issubnormal(x))
        end
    end
end

# math functions

@testset "real arithmetic" begin
    for T in (Float32, Float64, BigFloat)
        half = 1/convert(T,2)
        third = 1/convert(T,3)

        for f in (:+, :-, :abs, :abs2, :conj, :inv, :sign,
                  :acos, :asin, :asinh, :atan, :atanh, :cbrt, :cos, :cosh,
                  :exp10, :exp2, :exp, :expm1, :lgamma, :log10, :log1p,
                  :log2, :log, :sin, :sinh, :sqrt, :tan, :tanh)
            @eval begin
                @test @fastmath($f($half)) ≈ $f($half)
                @test @fastmath($f($third)) ≈ $f($third)
            end
        end
        for f in (:acosh,)
            @eval begin
                @test @fastmath($f(1+$half)) ≈ $f(1+$half)
                @test @fastmath($f(1+$third)) ≈ $f(1+$third)
            end
        end
        for f in (:+, :-, :*, :/, :%, :(==), :!=, :<, :<=, :>, :>=, :^,
                  :atan2, :hypot, :max, :min, :log)
            @eval begin
                @test @fastmath($f($half, $third)) ≈ $f($half, $third)
                @test @fastmath($f($third, $half)) ≈ $f($third, $half)
            end
        end
        for f in (:minmax,)
            @eval begin
                @test @fastmath($f($half, $third)[1]) ≈ $f($half, $third)[1]
                @test @fastmath($f($half, $third)[2]) ≈ $f($half, $third)[2]
                @test @fastmath($f($third, $half)[1]) ≈ $f($third, $half)[1]
                @test @fastmath($f($third, $half)[2]) ≈ $f($third, $half)[2]
            end
        end
    end
end
@testset "complex arithmetic" begin
    for T in (ComplexF32, ComplexF64, Complex{BigFloat})
        half = (1+1im)/T(2)
        third = (1-1im)/T(3)

        # some of these functions promote their result to double
        # precision, but we want to check equality at precision T
        rtol = Base.rtoldefault(real(T))

        for f in (:+, :-, :abs, :abs2, :conj, :inv, :sign,
                  :acos, :acosh, :asin, :asinh, :atan, :atanh, :cis, :cos,
                  :cosh, :exp10, :exp2, :exp, :expm1, :log10, :log1p,
                  :log2, :log, :sin, :sinh, :sqrt, :tan, :tanh)
            @eval begin
                @test @fastmath($f($half)) ≈ $f($half) rtol=$rtol
                @test @fastmath($f($third)) ≈ $f($third) rtol=$rtol
            end
        end
        for f in (:+, :-, :*, :/, :(==), :!=, :^, :log)
            @eval begin
                @test @fastmath($f($half, $third)) ≈ $f($half, $third) rtol=$rtol
                @test @fastmath($f($third, $half)) ≈ $f($third, $half) rtol=$rtol
            end
        end
    end
end
@testset "mixed real/complex arithmetic" begin
    for T in (Float32, Float64, BigFloat)
        CT = Complex{T}
        half = 1/T(2)
        third = 1/T(3)
        chalf = (1+1im)/CT(2)
        cthird = (1-1im)/CT(3)

        for f in (:+, :-, :*, :/, :(==), :!=, :^, :log)
            @eval begin
                @test @fastmath($f($chalf, $third)) ≈ $f($chalf, $third)
                @test @fastmath($f($half, $cthird)) ≈ $f($half, $cthird)
                @test @fastmath($f($cthird, $half)) ≈ $f($cthird, $half)
                @test @fastmath($f($third, $chalf)) ≈ $f($third, $chalf)
            end
        end

        @test @fastmath(third^3) ≈ third^3
        @test @fastmath(chalf/third) ≈ chalf/third
        @test @fastmath(chalf^3) ≈ chalf^3
        @test @fastmath(cis(third)) ≈ cis(third)
    end
end
@testset "issue #10544" begin
    a = fill(1.,2,2)
    b = fill(1.,2,2)
    @test @fastmath(a[1] += 2.0) ≈ (b[1] += 2.0)
    @test @fastmath(a[2] -= 2.0) ≈ (b[2] -= 2.0)
    @test @fastmath(a[1,1] *= 2.0) ≈ (b[1,1] *= 2.0)
    @test @fastmath(a[2,2] /= 2.0) ≈ (b[2,2] /= 2.0)
    @test @fastmath(a[1,2] ^= 2.0) ≈ (b[1,2] ^= 2.0)

    # test fallthrough for unsupported ops
    local c = 0
    @test @fastmath(c |= 1) == 1
end

@testset "issue #23218" begin
    a = zeros(1)
    b = [1.0]
    idx = (1,)
    @fastmath a[idx...] += b[idx...]
    @test a == b
end

@testset "literal powers" begin
    @test @fastmath(2^-2) == @fastmath(2.0^-2) == 0.25
end

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork