text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/PDEAssembler.jl/blob/master/LICENSE

module PDEAssembler

using Reexport
@reexport using FEMBase
using TimerOutputs

include("poisson.jl")

"""
    get_global_matrices(problems, time)

Assemble a set of problems and return assembles global matrices M for mass,
C for damping, K for stiffness, and right hand side vector f. Semi-discrete
equation of motion is then

    Mu'' + Cu' + Ku = f

Problems can describe a PDE in whole domain or in part of it. It can also describe
a boundary of the domain.
"""

function get_global_matrices(problems, time)

    @timeit "assemble problems" for problem in problems
        empty!(problem.assembly)
        assemble!(problem, time)
        assemble_mass_matrix!(problem, time)
    end

    @timeit "combine field problems" begin
        M = SparseMatrixCOO()
        K = SparseMatrixCOO()
        f = SparseMatrixCOO()
        C1 = SparseMatrixCOO()
        C2 = SparseMatrixCOO()
        D = SparseMatrixCOO()
        g = SparseMatrixCOO()
        for problem in problems
            A = problem.assembly
            append!(M, A.M)
            append!(K, A.K)
            append!(K, A.Kg)
            append!(f, A.f)
            append!(C1, A.C1)
            append!(C2, A.C2)
            append!(D, A.D)
            append!(g, A.g)
        end
    end

    @timeit "create LinearSystem(dim)" begin

        # Based on the nodal/dof mapping of PDEs, there might be empty rows in
        # stiffness matrix giving singular solution. For that reason we have
        # to find nonzero rows and add 1 to diagonal of those to indicate that
        # this particular dof is "disabled"

        dim = size(K, 1)
        nonzero_rows = zeros(dim)
        for i in K.I
            nonzero_rows[i] = 1.0
        end
        zero_indices = find(nonzero_rows .== 0.0)
        append!(K.I, zero_indices)
        append!(K.J, zero_indices)
        append!(K.V, ones(zero_indices))

        ls = LinearSystem(dim)
        ls.M = sparse(M, dim, dim)
        ls.K = sparse(K, dim, dim)
        ls.f = sparse(f, dim, 1)
        ls.C1 = sparse(C1, dim, dim)
        ls.C2 = sparse(C2, dim, dim)
        ls.D = sparse(D, dim, dim)
        ls.g = sparse(g, dim, 1)
    end

    # Now we have a linear system of equations to solve problem
    #
    # M*u'' + C*u' + K*u + C1*la = f
    # C2*u + D*la = g
    #
    # for u and la. In a case C1 == C2 and D = 0, we have a saddle
    # point problem where it is possible to either eliminate Lagrange
    # multipliers or solve constrained optimization problem using penalty
    # method. In simples kind of boundary conditions, constraint matrix is
    # a diagonal matrix having 1.0 in those diagonal entries that are
    # constrained.

    if ls.C1 == ls.C2' && nnz(ls.D) == 0
        ls.K += 1e36*ls.C1'*ls.C1
        ls.f += 1e36*ls.C1'*ls.g
    else
        error("Cannot handle boundary conditions using penalty method.")
    end

    # TODO: we are missing damping matrix ?!
    return ls.M, spzeros(dim,dim), ls.K, ls.f
end

export get_global_matrices

end


text:
function solve_poisson_FFT( Gvec::GvectorsT, rhoR )
    #
    G2 = Gvec.G2
    Ns = Gvec.Ns
    Npoints = prod(Ns)
    #
    ctmp = 4.0*pi*R_to_G( Ns, rhoR )
    #
    for ip = 2:Npoints
        ctmp[ip] = ctmp[ip]/G2[ip]
    end
    ctmp[1] = 0.0
    return real( G_to_R(Ns,ctmp) )
end

# In case we forget to convert the input, we convert it in this version
function R_to_G( Ns::Array{Int,1}, fR_::Array{Float64,1} )
    fR = convert(Array{ComplexF64,1},fR_)
    out = reshape( fft( reshape(fR,Ns[1],Ns[2],Ns[3]) ), size(fR) )
end

function G_to_R( Ns::Array{Int,1}, fG::Array{ComplexF64,1} )
    out = reshape( ifft( reshape(fG,Ns[1],Ns[2],Ns[3]) ), size(fG) )
end

ext: jl


text:
using FastBEAST
using StaticArrays
using LinearAlgebra

using Printf

function OneoverRkernel(testpoint::SVector{3,T}, sourcepoint::SVector{3,T}) where T
    if isapprox(testpoint, sourcepoint, rtol=eps()*1e1)
        return 0.0
    else
        return 1.0 / (norm(testpoint - sourcepoint))
    end
end

function assembler(kernel, testpoints, sourcepoints)
    kernelmatrix = zeros(promote_type(eltype(testpoints[1]),eltype(sourcepoints[1])),
                length(testpoints), length(sourcepoints))

    for i = 1:length(testpoints)
        for j = 1:length(sourcepoints)
            kernelmatrix[i,j] = kernel(testpoints[i], sourcepoints[j])
        end
    end
    return kernelmatrix
end

function assembler(kernel, matrix, testpoints, sourcepoints)
    for i = 1:length(testpoints)
        for j = 1:length(sourcepoints)
            matrix[i,j] = kernel(testpoints[i], sourcepoints[j])
        end
    end
end

##

N =  100
NT = N

spoints = [@SVector rand(3) for i = 1:N]
tpoints = 0.1*[@SVector rand(3) for i = 1:NT] + [1.0*SVector(3.5, 3.5, 3.5) for i = 1:NT]

@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(OneoverRkernel, matrix, tpoints[tdata], spoints[sdata])
stree = create_tree(spoints, treeoptions=BoxTreeOptions(nmin=50))
ttree = create_tree(tpoints, treeoptions=BoxTreeOptions(nmin=50))
kmat = assembler(OneoverRkernel, tpoints, spoints)
hmat = HMatrix(OneoverRkernelassembler, ttree, stree, compressor=:naive, T=Float64)

@printf("Accuracy test: %.2e\n", estimate_reldifference(hmat,kmat))
@printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)

##
N = 3000
NT = N

spoints = [@SVector rand(3) for i = 1:N]
##
@views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(OneoverRkernel, matrix, spoints[tdata], spoints[sdata])
stree = create_tree(spoints, treeoptions=BoxTreeOptions(nmin=100))
kmat = assembler(OneoverRkernel, spoints, spoints)
hmat = HMatrix(OneoverRkernelassembler, stree, stree, compressor=:aca, T=Float64)

@printf("Accuracy test: %.2e\n", estimate_reldifference(hmat,kmat))
@printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)

##
totalspoints = 0
for child in stree.children
    totalspoints += length(child.data)
end

@assert totalspoints == N

##

function hmatrix3d_benchmark(N)
    spoints = [@SVector rand(3) for i = 1:N]

    @views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(OneoverRkernel, matrix, spoints[tdata], spoints[sdata])
    stree = create_tree(spoints, nmin=400)

    @time hmat = HMatrix(OneoverRkernelassembler, stree, stree, compressor=:aca, tol=1e-4, T=Float64)

    @printf("Memory usage: %.2f GiB\n", nnz(hmat)*8/(1024*1024*1024.0))
    @printf("Compression rate: %.2f %%\n", compressionrate(hmat)*100)
    return nothing #stree, hmat
end

max_stars_repo_name: JoshuaTetzner/FastBEAST

max_issues_repo_name: JoshuaTetzner/FastBEAST

max_forks_repo_name: JoshuaTetzner/FastBEAST


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct AggregateDistributedGenerationA <: DynamicInjection
        name::String
        Pf_Flag::Int
        Freq_Flag::Int
        PQ_Flag::Int
        Gen_Flag::Int
        Vtrip_Flag::Int
        Ftrip_Flag::Int
        T_rv::Float64
        Trf::Float64
        dbd_pnts::Tuple{Float64, Float64}
        K_qv::Float64
        Tp::Float64
        T_iq::Float64
        D_dn::Float64
        D_up::Float64
        fdbd_pnts::Tuple{Float64, Float64}
        fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Tpord::Float64
        Kpg::Float64
        Kig::Float64
        I_max::Float64
        vl_pnts::Vector{Tuple{Float64,Float64}}
        vh_pnts::Vector{Tuple{Float64,Float64}}
        Vrfrac::Float64
        fl::Float64
        fh::Float64
        tfl::Float64
        tfh::Float64
        Tg::Float64
        rrpwr::Float64
        Tv::Float64
        Vpr::Float64
        Iq_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        V_ref::Float64
        Pfa_ref::Float64
        Q_ref::Float64
        P_ref::Float64
        base_power::Float64
        states::Vector{Symbol}
        n_states::Int
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

Parameters of the DERA1 model in PSS/E

# Arguments
- `name::String`
- `Pf_Flag::Int`: Flag for Power Factor Control, validation range: `(0, 1)`
- `Freq_Flag::Int`: Flag to enable/disable frequency control, validation range: `(0, 1)`
- `PQ_Flag::Int`: Flag used to enforce maximum current, validation range: `(0, 1)`
- `Gen_Flag::Int`: Flag to specify generator or storage, validation range: `(0, 1)`
- `Vtrip_Flag::Int`: Flag to enable/disable voltage trip logic, validation range: `(0, 1)`
- `Ftrip_Flag::Int`: Flag to enable/disable frequency trip logic, validation range: `(0, 1)`
- `T_rv::Float64`: Voltage measurement transducer time constant, validation range: `(0, nothing)`
- `Trf::Float64`: Frequency measurement transducer time constant, validation range: `(0, nothing)`
- `dbd_pnts::Tuple{Float64, Float64}`: Voltage deadband thresholds `(dbd1, dbd2)`
- `K_qv::Float64`: Proportional voltage control gain (pu), validation range: `(0, nothing)`
- `Tp::Float64`: Power measurement transducer time constant, validation range: `(0, nothing)`
- `T_iq::Float64`: Time constant for low-pass filter for state q_V when QFlag = 0, validation range: `(0, nothing)`
- `D_dn::Float64`: Reciprocal of droop for over-frequency conditions (>0) (pu), validation range: `(0, nothing)`
- `D_up::Float64`: Reciprocal of droop for under-frequency conditions <=0) (pu), validation range: `(0, nothing)`
- `fdbd_pnts::Tuple{Float64, Float64}`: Frequency control deadband thresholds `(fdbd1, fdbd2)`
- `fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Frequency error limits (femin, femax)
- `P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Power limits (Pmin, Pmax)
- `dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Power reference ramp rate limits (dPmin, dPmax)
- `Tpord::Float64`: Power filter time constant, validation range: `(0, nothing)`
- `Kpg::Float64`: PI controller proportional gain (pu), validation range: `(0, nothing)`
- `Kig::Float64`: PI controller integral gain (pu), validation range: `(0, nothing)`
- `I_max::Float64`: Maximum limit on total converter current (pu), validation range: `(0, nothing)`
- `vl_pnts::Vector{Tuple{Float64,Float64}}`: Low voltage cutout points `[(tv10, vl0), (tv11, vl1)]`
- `vh_pnts::Vector{Tuple{Float64,Float64}}`: High voltage cutout points `[(tvh0, vh0), (tvh1, vh1)]`
- `Vrfrac::Float64`: Fraction of device that recovers after voltage comes back to within vl1 < V < vh1 (0 <= Vrfrac <= 1), validation range: `(0, 1)`
- `fl::Float64`: Inverter frequency break-point for low frequency cut-out (Hz), validation range: `(0, nothing)`
- `fh::Float64`: Inverter frequency break-point for high frequency cut-out (Hz), validation range: `(0, nothing)`
- `tfl::Float64`: Low frequency cut-out timer corresponding to frequency fl (s), validation range: `(0, nothing)`
- `tfh::Float64`: High frequency cut-out timer corresponding to frequency fh (s), validation range: `(0, nothing)`
- `Tg::Float64`: Current control time constant (to represent behavior of inner control loops) (> 0) (s), validation range: `(0, nothing)`
- `rrpwr::Float64`: Ramp rate for real power increase following a fault (pu/s), validation range: `(0, nothing)`
- `Tv::Float64`: Time constant on the output of the multiplier (s), validation range: `(0, nothing)`
- `Vpr::Float64`: Voltage below which frequency tripping is disabled (pu), validation range: `(0, nothing)`
- `Iq_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Reactive current injection limits (Iqll, Iqhl)
- `V_ref::Float64`: User defined voltage reference. If 0, PSID initializes to initial terminal voltage, validation range: `(0, nothing)`
- `Pfa_ref::Float64`: Reference power factor, validation range: `(0, nothing)`
- `Q_ref::Float64`: Reference reactive power, in pu, validation range: `(0, nothing)`
- `P_ref::Float64`: Reference active power, in pu, validation range: `(0, nothing)`
- `base_power::Float64`: Base power
- `states::Vector{Symbol}`: The states of AggregateDistributedGenerationA depends on the Flags
- `n_states::Int`: The states of AggregateDistributedGenerationA depends on the Flags
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AggregateDistributedGenerationA <: DynamicInjection
    name::String
    "Flag for Power Factor Control"
    Pf_Flag::Int
    "Flag to enable/disable frequency control"
    Freq_Flag::Int
    "Flag used to enforce maximum current"
    PQ_Flag::Int
    "Flag to specify generator or storage"
    Gen_Flag::Int
    "Flag to enable/disable voltage trip logic"
    Vtrip_Flag::Int
    "Flag to enable/disable frequency trip logic"
    Ftrip_Flag::Int
    "Voltage measurement transducer time constant"
    T_rv::Float64
    "Frequency measurement transducer time constant"
    Trf::Float64
    "Voltage deadband thresholds `(dbd1, dbd2)`"
    dbd_pnts::Tuple{Float64, Float64}
    "Proportional voltage control gain (pu)"
    K_qv::Float64
    "Power measurement transducer time constant"
    Tp::Float64
    "Time constant for low-pass filter for state q_V when QFlag = 0"
    T_iq::Float64
    "Reciprocal of droop for over-frequency conditions (>0) (pu)"
    D_dn::Float64
    "Reciprocal of droop for under-frequency conditions <=0) (pu)"
    D_up::Float64
    "Frequency control deadband thresholds `(fdbd1, fdbd2)`"
    fdbd_pnts::Tuple{Float64, Float64}
    "Frequency error limits (femin, femax)"
    fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power limits (Pmin, Pmax)"
    P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power reference ramp rate limits (dPmin, dPmax)"
    dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power filter time constant"
    Tpord::Float64
    "PI controller proportional gain (pu)"
    Kpg::Float64
    "PI controller integral gain (pu)"
    Kig::Float64
    "Maximum limit on total converter current (pu)"
    I_max::Float64
    "Low voltage cutout points `[(tv10, vl0), (tv11, vl1)]`"
    vl_pnts::Vector{Tuple{Float64,Float64}}
    "High voltage cutout points `[(tvh0, vh0), (tvh1, vh1)]`"
    vh_pnts::Vector{Tuple{Float64,Float64}}
    "Fraction of device that recovers after voltage comes back to within vl1 < V < vh1 (0 <= Vrfrac <= 1)"
    Vrfrac::Float64
    "Inverter frequency break-point for low frequency cut-out (Hz)"
    fl::Float64
    "Inverter frequency break-point for high frequency cut-out (Hz)"
    fh::Float64
    "Low frequency cut-out timer corresponding to frequency fl (s)"
    tfl::Float64
    "High frequency cut-out timer corresponding to frequency fh (s)"
    tfh::Float64
    "Current control time constant (to represent behavior of inner control loops) (> 0) (s)"
    Tg::Float64
    "Ramp rate for real power increase following a fault (pu/s)"
    rrpwr::Float64
    "Time constant on the output of the multiplier (s)"
    Tv::Float64
    "Voltage below which frequency tripping is disabled (pu)"
    Vpr::Float64
    "Reactive current injection limits (Iqll, Iqhl)"
    Iq_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "User defined voltage reference. If 0, PSID initializes to initial terminal voltage"
    V_ref::Float64
    "Reference power factor"
    Pfa_ref::Float64
    "Reference reactive power, in pu"
    Q_ref::Float64
    "Reference active power, in pu"
    P_ref::Float64
    "Base power"
    base_power::Float64
    "The states of AggregateDistributedGenerationA depends on the Flags"
    states::Vector{Symbol}
    "The states of AggregateDistributedGenerationA depends on the Flags"
    n_states::Int
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AggregateDistributedGenerationA(name, Pf_Flag, Freq_Flag, PQ_Flag, Gen_Flag, Vtrip_Flag, Ftrip_Flag, T_rv, Trf, dbd_pnts, K_qv, Tp, T_iq, D_dn, D_up, fdbd_pnts, fe_lim, P_lim, dP_lim, Tpord, Kpg, Kig, I_max, vl_pnts, vh_pnts, Vrfrac, fl, fh, tfl, tfh, Tg, rrpwr, Tv, Vpr, Iq_lim, V_ref=1.0, Pfa_ref=0.0, Q_ref=0.0, P_ref=1.0, base_power=100.0, ext=Dict{String, Any}(), )
    AggregateDistributedGenerationA(name, Pf_Flag, Freq_Flag, PQ_Flag, Gen_Flag, Vtrip_Flag, Ftrip_Flag, T_rv, Trf, dbd_pnts, K_qv, Tp, T_iq, D_dn, D_up, fdbd_pnts, fe_lim, P_lim, dP_lim, Tpord, Kpg, Kig, I_max, vl_pnts, vh_pnts, Vrfrac, fl, fh, tfl, tfh, Tg, rrpwr, Tv, Vpr, Iq_lim, V_ref, Pfa_ref, Q_ref, P_ref, base_power, ext, PowerSystems.get_AggregateDistributedGenerationA_states(Freq_Flag)[1], PowerSystems.get_AggregateDistributedGenerationA_states(Freq_Flag)[2], InfrastructureSystemsInternal(), )
end

function AggregateDistributedGenerationA(; name, Pf_Flag, Freq_Flag, PQ_Flag, Gen_Flag, Vtrip_Flag, Ftrip_Flag, T_rv, Trf, dbd_pnts, K_qv, Tp, T_iq, D_dn, D_up, fdbd_pnts, fe_lim, P_lim, dP_lim, Tpord, Kpg, Kig, I_max, vl_pnts, vh_pnts, Vrfrac, fl, fh, tfl, tfh, Tg, rrpwr, Tv, Vpr, Iq_lim, V_ref=1.0, Pfa_ref=0.0, Q_ref=0.0, P_ref=1.0, base_power=100.0, states=PowerSystems.get_AggregateDistributedGenerationA_states(Freq_Flag)[1], n_states=PowerSystems.get_AggregateDistributedGenerationA_states(Freq_Flag)[2], ext=Dict{String, Any}(), internal=InfrastructureSystemsInternal(), )
    AggregateDistributedGenerationA(name, Pf_Flag, Freq_Flag, PQ_Flag, Gen_Flag, Vtrip_Flag, Ftrip_Flag, T_rv, Trf, dbd_pnts, K_qv, Tp, T_iq, D_dn, D_up, fdbd_pnts, fe_lim, P_lim, dP_lim, Tpord, Kpg, Kig, I_max, vl_pnts, vh_pnts, Vrfrac, fl, fh, tfl, tfh, Tg, rrpwr, Tv, Vpr, Iq_lim, V_ref, Pfa_ref, Q_ref, P_ref, base_power, states, n_states, ext, internal, )
end

# Constructor for demo purposes; non-functional.
function AggregateDistributedGenerationA(::Nothing)
    AggregateDistributedGenerationA(;
        name="init",
        Pf_Flag=0,
        Freq_Flag=0,
        PQ_Flag=0,
        Gen_Flag=0,
        Vtrip_Flag=0,
        Ftrip_Flag=0,
        T_rv=0,
        Trf=0,
        dbd_pnts=(0.0, 0.0),
        K_qv=0,
        Tp=0,
        T_iq=0,
        D_dn=0,
        D_up=0,
        fdbd_pnts=(0.0, 0.0),
        fe_lim=(min=0.0, max=0.0),
        P_lim=(min=0.0, max=0.0),
        dP_lim=(min=0.0, max=0.0),
        Tpord=0,
        Kpg=0,
        Kig=0,
        I_max=0,
        vl_pnts=[(0.0, 0.0), (0.0, 0.0)],
        vh_pnts=[(0.0, 0.0), (0.0, 0.0)],
        Vrfrac=0,
        fl=0,
        fh=0,
        tfl=0,
        tfh=0,
        Tg=0,
        rrpwr=0,
        Tv=0,
        Vpr=0,
        Iq_lim=(min=0.0, max=0.0),
        V_ref=0,
        Pfa_ref=0,
        Q_ref=0,
        P_ref=0,
        base_power=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`AggregateDistributedGenerationA`](@ref) `name`."""
get_name(value::AggregateDistributedGenerationA) = value.name
"""Get [`AggregateDistributedGenerationA`](@ref) `Pf_Flag`."""
get_Pf_Flag(value::AggregateDistributedGenerationA) = value.Pf_Flag
"""Get [`AggregateDistributedGenerationA`](@ref) `Freq_Flag`."""
get_Freq_Flag(value::AggregateDistributedGenerationA) = value.Freq_Flag
"""Get [`AggregateDistributedGenerationA`](@ref) `PQ_Flag`."""
get_PQ_Flag(value::AggregateDistributedGenerationA) = value.PQ_Flag
"""Get [`AggregateDistributedGenerationA`](@ref) `Gen_Flag`."""
get_Gen_Flag(value::AggregateDistributedGenerationA) = value.Gen_Flag
"""Get [`AggregateDistributedGenerationA`](@ref) `Vtrip_Flag`."""
get_Vtrip_Flag(value::AggregateDistributedGenerationA) = value.Vtrip_Flag
"""Get [`AggregateDistributedGenerationA`](@ref) `Ftrip_Flag`."""
get_Ftrip_Flag(value::AggregateDistributedGenerationA) = value.Ftrip_Flag
"""Get [`AggregateDistributedGenerationA`](@ref) `T_rv`."""
get_T_rv(value::AggregateDistributedGenerationA) = value.T_rv
"""Get [`AggregateDistributedGenerationA`](@ref) `Trf`."""
get_Trf(value::AggregateDistributedGenerationA) = value.Trf
"""Get [`AggregateDistributedGenerationA`](@ref) `dbd_pnts`."""
get_dbd_pnts(value::AggregateDistributedGenerationA) = value.dbd_pnts
"""Get [`AggregateDistributedGenerationA`](@ref) `K_qv`."""
get_K_qv(value::AggregateDistributedGenerationA) = value.K_qv
"""Get [`AggregateDistributedGenerationA`](@ref) `Tp`."""
get_Tp(value::AggregateDistributedGenerationA) = value.Tp
"""Get [`AggregateDistributedGenerationA`](@ref) `T_iq`."""
get_T_iq(value::AggregateDistributedGenerationA) = value.T_iq
"""Get [`AggregateDistributedGenerationA`](@ref) `D_dn`."""
get_D_dn(value::AggregateDistributedGenerationA) = value.D_dn
"""Get [`AggregateDistributedGenerationA`](@ref) `D_up`."""
get_D_up(value::AggregateDistributedGenerationA) = value.D_up
"""Get [`AggregateDistributedGenerationA`](@ref) `fdbd_pnts`."""
get_fdbd_pnts(value::AggregateDistributedGenerationA) = value.fdbd_pnts
"""Get [`AggregateDistributedGenerationA`](@ref) `fe_lim`."""
get_fe_lim(value::AggregateDistributedGenerationA) = value.fe_lim
"""Get [`AggregateDistributedGenerationA`](@ref) `P_lim`."""
get_P_lim(value::AggregateDistributedGenerationA) = value.P_lim
"""Get [`AggregateDistributedGenerationA`](@ref) `dP_lim`."""
get_dP_lim(value::AggregateDistributedGenerationA) = value.dP_lim
"""Get [`AggregateDistributedGenerationA`](@ref) `Tpord`."""
get_Tpord(value::AggregateDistributedGenerationA) = value.Tpord
"""Get [`AggregateDistributedGenerationA`](@ref) `Kpg`."""
get_Kpg(value::AggregateDistributedGenerationA) = value.Kpg
"""Get [`AggregateDistributedGenerationA`](@ref) `Kig`."""
get_Kig(value::AggregateDistributedGenerationA) = value.Kig
"""Get [`AggregateDistributedGenerationA`](@ref) `I_max`."""
get_I_max(value::AggregateDistributedGenerationA) = value.I_max
"""Get [`AggregateDistributedGenerationA`](@ref) `vl_pnts`."""
get_vl_pnts(value::AggregateDistributedGenerationA) = value.vl_pnts
"""Get [`AggregateDistributedGenerationA`](@ref) `vh_pnts`."""
get_vh_pnts(value::AggregateDistributedGenerationA) = value.vh_pnts
"""Get [`AggregateDistributedGenerationA`](@ref) `Vrfrac`."""
get_Vrfrac(value::AggregateDistributedGenerationA) = value.Vrfrac
"""Get [`AggregateDistributedGenerationA`](@ref) `fl`."""
get_fl(value::AggregateDistributedGenerationA) = value.fl
"""Get [`AggregateDistributedGenerationA`](@ref) `fh`."""
get_fh(value::AggregateDistributedGenerationA) = value.fh
"""Get [`AggregateDistributedGenerationA`](@ref) `tfl`."""
get_tfl(value::AggregateDistributedGenerationA) = value.tfl
"""Get [`AggregateDistributedGenerationA`](@ref) `tfh`."""
get_tfh(value::AggregateDistributedGenerationA) = value.tfh
"""Get [`AggregateDistributedGenerationA`](@ref) `Tg`."""
get_Tg(value::AggregateDistributedGenerationA) = value.Tg
"""Get [`AggregateDistributedGenerationA`](@ref) `rrpwr`."""
get_rrpwr(value::AggregateDistributedGenerationA) = value.rrpwr
"""Get [`AggregateDistributedGenerationA`](@ref) `Tv`."""
get_Tv(value::AggregateDistributedGenerationA) = value.Tv
"""Get [`AggregateDistributedGenerationA`](@ref) `Vpr`."""
get_Vpr(value::AggregateDistributedGenerationA) = value.Vpr
"""Get [`AggregateDistributedGenerationA`](@ref) `Iq_lim`."""
get_Iq_lim(value::AggregateDistributedGenerationA) = value.Iq_lim
"""Get [`AggregateDistributedGenerationA`](@ref) `V_ref`."""
get_V_ref(value::AggregateDistributedGenerationA) = value.V_ref
"""Get [`AggregateDistributedGenerationA`](@ref) `Pfa_ref`."""
get_Pfa_ref(value::AggregateDistributedGenerationA) = value.Pfa_ref
"""Get [`AggregateDistributedGenerationA`](@ref) `Q_ref`."""
get_Q_ref(value::AggregateDistributedGenerationA) = value.Q_ref
"""Get [`AggregateDistributedGenerationA`](@ref) `P_ref`."""
get_P_ref(value::AggregateDistributedGenerationA) = value.P_ref
"""Get [`AggregateDistributedGenerationA`](@ref) `base_power`."""
get_base_power(value::AggregateDistributedGenerationA) = value.base_power
"""Get [`AggregateDistributedGenerationA`](@ref) `states`."""
get_states(value::AggregateDistributedGenerationA) = value.states
"""Get [`AggregateDistributedGenerationA`](@ref) `n_states`."""
get_n_states(value::AggregateDistributedGenerationA) = value.n_states
"""Get [`AggregateDistributedGenerationA`](@ref) `ext`."""
get_ext(value::AggregateDistributedGenerationA) = value.ext
"""Get [`AggregateDistributedGenerationA`](@ref) `internal`."""
get_internal(value::AggregateDistributedGenerationA) = value.internal

"""Set [`AggregateDistributedGenerationA`](@ref) `Pf_Flag`."""
set_Pf_Flag!(value::AggregateDistributedGenerationA, val) = value.Pf_Flag = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Freq_Flag`."""
set_Freq_Flag!(value::AggregateDistributedGenerationA, val) = value.Freq_Flag = val
"""Set [`AggregateDistributedGenerationA`](@ref) `PQ_Flag`."""
set_PQ_Flag!(value::AggregateDistributedGenerationA, val) = value.PQ_Flag = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Gen_Flag`."""
set_Gen_Flag!(value::AggregateDistributedGenerationA, val) = value.Gen_Flag = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Vtrip_Flag`."""
set_Vtrip_Flag!(value::AggregateDistributedGenerationA, val) = value.Vtrip_Flag = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Ftrip_Flag`."""
set_Ftrip_Flag!(value::AggregateDistributedGenerationA, val) = value.Ftrip_Flag = val
"""Set [`AggregateDistributedGenerationA`](@ref) `T_rv`."""
set_T_rv!(value::AggregateDistributedGenerationA, val) = value.T_rv = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Trf`."""
set_Trf!(value::AggregateDistributedGenerationA, val) = value.Trf = val
"""Set [`AggregateDistributedGenerationA`](@ref) `dbd_pnts`."""
set_dbd_pnts!(value::AggregateDistributedGenerationA, val) = value.dbd_pnts = val
"""Set [`AggregateDistributedGenerationA`](@ref) `K_qv`."""
set_K_qv!(value::AggregateDistributedGenerationA, val) = value.K_qv = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Tp`."""
set_Tp!(value::AggregateDistributedGenerationA, val) = value.Tp = val
"""Set [`AggregateDistributedGenerationA`](@ref) `T_iq`."""
set_T_iq!(value::AggregateDistributedGenerationA, val) = value.T_iq = val
"""Set [`AggregateDistributedGenerationA`](@ref) `D_dn`."""
set_D_dn!(value::AggregateDistributedGenerationA, val) = value.D_dn = val
"""Set [`AggregateDistributedGenerationA`](@ref) `D_up`."""
set_D_up!(value::AggregateDistributedGenerationA, val) = value.D_up = val
"""Set [`AggregateDistributedGenerationA`](@ref) `fdbd_pnts`."""
set_fdbd_pnts!(value::AggregateDistributedGenerationA, val) = value.fdbd_pnts = val
"""Set [`AggregateDistributedGenerationA`](@ref) `fe_lim`."""
set_fe_lim!(value::AggregateDistributedGenerationA, val) = value.fe_lim = val
"""Set [`AggregateDistributedGenerationA`](@ref) `P_lim`."""
set_P_lim!(value::AggregateDistributedGenerationA, val) = value.P_lim = val
"""Set [`AggregateDistributedGenerationA`](@ref) `dP_lim`."""
set_dP_lim!(value::AggregateDistributedGenerationA, val) = value.dP_lim = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Tpord`."""
set_Tpord!(value::AggregateDistributedGenerationA, val) = value.Tpord = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Kpg`."""
set_Kpg!(value::AggregateDistributedGenerationA, val) = value.Kpg = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Kig`."""
set_Kig!(value::AggregateDistributedGenerationA, val) = value.Kig = val
"""Set [`AggregateDistributedGenerationA`](@ref) `I_max`."""
set_I_max!(value::AggregateDistributedGenerationA, val) = value.I_max = val
"""Set [`AggregateDistributedGenerationA`](@ref) `vl_pnts`."""
set_vl_pnts!(value::AggregateDistributedGenerationA, val) = value.vl_pnts = val
"""Set [`AggregateDistributedGenerationA`](@ref) `vh_pnts`."""
set_vh_pnts!(value::AggregateDistributedGenerationA, val) = value.vh_pnts = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Vrfrac`."""
set_Vrfrac!(value::AggregateDistributedGenerationA, val) = value.Vrfrac = val
"""Set [`AggregateDistributedGenerationA`](@ref) `fl`."""
set_fl!(value::AggregateDistributedGenerationA, val) = value.fl = val
"""Set [`AggregateDistributedGenerationA`](@ref) `fh`."""
set_fh!(value::AggregateDistributedGenerationA, val) = value.fh = val
"""Set [`AggregateDistributedGenerationA`](@ref) `tfl`."""
set_tfl!(value::AggregateDistributedGenerationA, val) = value.tfl = val
"""Set [`AggregateDistributedGenerationA`](@ref) `tfh`."""
set_tfh!(value::AggregateDistributedGenerationA, val) = value.tfh = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Tg`."""
set_Tg!(value::AggregateDistributedGenerationA, val) = value.Tg = val
"""Set [`AggregateDistributedGenerationA`](@ref) `rrpwr`."""
set_rrpwr!(value::AggregateDistributedGenerationA, val) = value.rrpwr = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Tv`."""
set_Tv!(value::AggregateDistributedGenerationA, val) = value.Tv = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Vpr`."""
set_Vpr!(value::AggregateDistributedGenerationA, val) = value.Vpr = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Iq_lim`."""
set_Iq_lim!(value::AggregateDistributedGenerationA, val) = value.Iq_lim = val
"""Set [`AggregateDistributedGenerationA`](@ref) `V_ref`."""
set_V_ref!(value::AggregateDistributedGenerationA, val) = value.V_ref = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Pfa_ref`."""
set_Pfa_ref!(value::AggregateDistributedGenerationA, val) = value.Pfa_ref = val
"""Set [`AggregateDistributedGenerationA`](@ref) `Q_ref`."""
set_Q_ref!(value::AggregateDistributedGenerationA, val) = value.Q_ref = val
"""Set [`AggregateDistributedGenerationA`](@ref) `P_ref`."""
set_P_ref!(value::AggregateDistributedGenerationA, val) = value.P_ref = val
"""Set [`AggregateDistributedGenerationA`](@ref) `base_power`."""
set_base_power!(value::AggregateDistributedGenerationA, val) = value.base_power = val
"""Set [`AggregateDistributedGenerationA`](@ref) `ext`."""
set_ext!(value::AggregateDistributedGenerationA, val) = value.ext = val


text:
using Logging
Logging.configure(level=DEBUG)
include("../calibrate_params.jl")
using CalibrateParameters

if !haskey(parameters, :S)
	parameters[:S] = 101
end

parameters[:numerical_zero] = 1e-12

parameters[:bp_weights] = [0.774074394803123; -0.201004684236153; -0.135080548288772; -0.0509519648766360]

# set random seed so that all scenarios comparable
srand(9181)
CalibrateParameters.calibrate_parameters!(parameters)

# adaptive step size. large lambda means large steps
parameters[:inner_step_size] = exp(-0.06*(parameters[:J]-1)^0.75)
# large substitution needs more dampening
parameters[:middle_step_size] = exp(-0.8*max(1,parameters[:sigma]))
parameters[:adjustment_step_size] = 0.25
# any deviation from sigma=1 needs more dampening
parameters[:outer_step_size] = exp(-0.5*(1+abs(log(parameters[:sigma]))))
# this is log points of average input price differences
parameters[:inner_tolerance]  = 0.0002
parameters[:middle_tolerance] = 0.0003
parameters[:adjustment_tolerance] = 0.0004
parameters[:outer_tolerance] = 0.001

# maximum number of iterations in each loop
parameters[:max_iter_inner] = 200
parameters[:max_iter_middle] = 100
parameters[:max_iter_adjustment] = 100
parameters[:max_iter_outer] = 100


text:
using Skewlia
using Random
using Test

@testset "SkewNode" begin
    h = Skewlia.SkewHeap()
    a = Skewlia.SkewNode(1)
    b = Skewlia.SkewNode(4)

    @test nothing + nothing === nothing
    @test a + nothing === a
    @test nothing + b === b
    @test a + b == Skewlia.SkewNode(1, b, nothing) # in order
    @test b + a == Skewlia.SkewNode(1, b, nothing) # out of order
end

@testset "SkewHeap" begin
    ordered = 0:9
    shuffled = Random.shuffle(ordered)

    h = Skewlia.SkewHeap()
    @test Skewlia.take!(h) === nothing
    @test Skewlia.peek(h) === nothing
    @test Skewlia.length(h) == 0
    @test Skewlia.is_empty(h)

    @test Skewlia.put!(h, shuffled...) == length(shuffled)
    @test Skewlia.length(h) == length(shuffled)

    size = length(shuffled)

    for i in ordered
        @test Skewlia.length(h) == size

        size -= 1

        @test Skewlia.peek(h) == i
        @test Skewlia.take!(h) == i
        @test Skewlia.peek(h) != i
        @test Skewlia.length(h) == size
    end

    @test Skewlia.take!(h) === nothing
    @test Skewlia.peek(h) === nothing
    @test Skewlia.length(h) == 0
    @test Skewlia.is_empty(h)
end

@testset "Iteration" begin
    ordered = 0:9
    shuffled = Random.shuffle(ordered)

    h = Skewlia.SkewHeap()
    Skewlia.put!(h, shuffled...)

    expect = 0
    for i in h
       @test i == expect
       expect += 1
    end
end

max_stars_repo_name: sysread/Skewlia

max_issues_repo_name: sysread/Skewlia

max_forks_repo_name: sysread/Skewlia


text:
mutable struct TestReward <: SeaPearl.AbstractReward
    value::Float32
end

TestReward(model::SeaPearl.CPModel) = TestReward(0)

function SeaPearl.set_reward!(::Type{SeaPearl.InitializingPhase}, lh::SeaPearl.LearnedHeuristic{SR, TestReward, A}, model::SeaPearl.CPModel) where {
    SR <: SeaPearl.AbstractStateRepresentation,
    A <: SeaPearl.ActionOutput
}
    lh.reward.value += 1
    nothing
end
function SeaPearl.set_reward!(::Type{SeaPearl.StepPhase}, lh::SeaPearl.LearnedHeuristic{SR, TestReward, A}, model::SeaPearl.CPModel) where {
    SR <: SeaPearl.AbstractStateRepresentation,
    A <: SeaPearl.ActionOutput
}
    lh.reward.value += 2
    nothing
end
function SeaPearl.set_reward!(::Type{SeaPearl.DecisionPhase}, lh::SeaPearl.LearnedHeuristic{SR, TestReward, A}, model::SeaPearl.CPModel) where {
    SR <: SeaPearl.AbstractStateRepresentation,
    A <: SeaPearl.ActionOutput
}
    lh.reward.value += 3
    nothing
end
function SeaPearl.set_reward!(::Type{SeaPearl.EndingPhase}, lh::SeaPearl.LearnedHeuristic{SR, TestReward, A}, model::SeaPearl.CPModel, symbol::Union{Nothing, Symbol}) where {
    SR <: SeaPearl.AbstractStateRepresentation,
    A <: SeaPearl.ActionOutput
}
    lh.reward.value += -5
    nothing
end

@testset "rewards.jl" begin

    include("defaultreward.jl")
    include("tsptwreward.jl")

    @testset "Custom reward" begin
        @testset "set_reward(StepPhase)" begin
            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            lh = SeaPearl.LearnedHeuristic(agent)
            SeaPearl.update_with_cpmodel!(lh, model)

            lh.reward.value = 0
            model.statistics.numberOfNodes = 1
            SeaPearl.set_reward!(SeaPearl.StepPhase, lh, model,:Infeasible)
            @test lh.reward.value == -10

            lh.reward.value = 0
            SeaPearl.set_reward!(SeaPearl.StepPhase, lh, model,:FoundSolution)
            @test lh.reward.value == -1

            lh.reward.value = 0
            SeaPearl.set_reward!(SeaPearl.StepPhase, lh, model,:Feasible)
            @test lh.reward.value == -1

            lh.reward.value = 0
            SeaPearl.set_reward!(SeaPearl.StepPhase, lh, model,:BackTracking)
            @test lh.reward.value == -1

        end

        @testset "set_reward!(InitializingPhase)" begin
            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, TestReward, SeaPearl.FixedOutput}(agent)
            SeaPearl.update_with_cpmodel!(lh, model)

            lh.reward.value = 0
            SeaPearl.set_reward!(SeaPearl.InitializingPhase, lh, model)
            @test lh.reward.value == 1
        end

        @testset "set_reward!(StepPhase)" begin
            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, TestReward, SeaPearl.FixedOutput}(agent)
            SeaPearl.update_with_cpmodel!(lh, model)

            lh.reward.value = 0
            SeaPearl.set_reward!(SeaPearl.StepPhase, lh, model)
            @test lh.reward.value == 2
        end

        @testset "set_reward!(DecisionPhase)" begin
            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, TestReward, SeaPearl.FixedOutput}(agent)
            SeaPearl.update_with_cpmodel!(lh, model)

            lh.reward.value = 0
            SeaPearl.set_reward!(SeaPearl.DecisionPhase, lh, model)
            @test lh.reward.value == 3
        end

        @testset "set_reward!(EndingPhase)" begin
            trailer = SeaPearl.Trailer()
            model = SeaPearl.CPModel(trailer)

            lh = SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, TestReward, SeaPearl.FixedOutput}(agent)
            SeaPearl.update_with_cpmodel!(lh, model)

            lh.reward.value = 6
            SeaPearl.set_reward!(SeaPearl.EndingPhase, lh, model, nothing)
            @test lh.reward.value == 1
        end
    end

end


text:
export iauTttdb
"""
Time scale transformation:  Terrestrial Time, TT, to Barycentric
Dynamical Time, TDB.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tt1,tt2    double    TT as a 2-part Julian Date
   dtr        double    TDB-TT in seconds

Returned:
   tdb1,tdb2  double    TDB as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Notes:

   1. tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tdb1,tdb2 follow
   suit.

   2. The argument dtr represents the quasi-periodic component of the
   GR transformation between TT and TCB.  It is dependent upon the
   adopted solar-system ephemeris, and can be obtained by numerical
   integration, by interrogating a precomputed time ephemeris or by
   evaluating a model such as that implemented in the SOFA function
   iauDtdb.   The quantity is dominated by an annual term of 1.7 ms
   amplitude.

   3. TDB is essentially the same as Teph, the time argument for the JPL
   solar system ephemerides.

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   IAU 2006 Resolution 3

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTttdb(double tt1, double tt2, double dtr,
#              double *tdb1, double *tdb2)

function iauTttdb(tt1::Real, tt2::Real, dtr::Real)
   ref_tdb1 = Ref{Float64}(0.0)
   ref_tdb2 = Ref{Float64}(0.0)

   status = ccall((:iauTttdb, libsofa_c), Cint,
                   (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tt1), convert(Float64, tt2), convert(Float64, dtr),
                   ref_tdb1, ref_tdb2)

   return status, ref_tdb1[], ref_tdb2[]
end


text:
struct Reaction
    product_stoich_coeffs::SparseMatrixCSC{Float64,Int64}
    reactant_stoich_coeffs::SparseMatrixCSC{Float64,Int64}
    reactant_orders::SparseMatrixCSC{Float64,Int64}
    is_reversible::Array{Bool,1}
    Arrhenius_coeffs::Array{Float64,2}
    Arrhenius_0::Array{Float64,2}
    Troe_::Array{Float64,2}
    index_three_body::Array{Int64,1}
    index_falloff::Array{Int64,1}
    index_falloff_Troe::Array{Int64,1}
    efficiencies_coeffs::SparseMatrixCSC{Float64,Int64}
    i_reactant::Array{Array{Int64,1},1}
    i_product::Array{Array{Int64,1},1}
    n_reactions::Int64
    vk::SparseMatrixCSC{Float64,Int64}
    vk_sum::Array{Float64,1}
end

abstract type Thermo end

struct Transport
    poly_order::Int64
    species_viscosities_poly::Array{Float64,2}
    thermal_conductivity_poly::Array{Float64,2}
    binary_diff_coeffs_poly::Array{Float64,2}
end

struct Solution
    n_species::Int64
    n_reactions::Int64
    MW::Vector{Float64}
    species_names::Vector{String}
    elements::Vector{String}
    ele_matrix::Array{Float64,2}
    thermo::Thermo
    trans::Transport
    reaction::Reaction
end


text:
using  HorizonSideRobots

function Obhod(r,HS)
    x = 0
    Flag = true
    while isborder(r,HorizonSide(HS)) && !isborder(r,HorizonSide((HS + 1) % 4))
        move!(r,HorizonSide((HS + 1) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 1) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(HS))
        move!(r,HorizonSide(HS))
    end
    while !isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 3) % 4))
        move!(r,HorizonSide(HS))
    end
    for i = 1:x
        move!(r,HorizonSide((HS+3) % 4))
    end
    return Flag
end

function Obhod3(r,HS,y)
    x = 0
    Flag = true
    while isborder(r,HorizonSide(HS)) && !isborder(r,HorizonSide((HS + 3) % 4))
        move!(r,HorizonSide((HS + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 3) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(HS))
        move!(r,HorizonSide(HS))
        y-=1
    end
    while !isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 1) % 4))
        move!(r,HorizonSide(HS))
        y-= 1
    end
    for i = 1:x
        move!(r,HorizonSide((HS+1) % 4))
    end
    return Flag,y
end

function Vozvrat(r,y,z)
    Flag = true
    for i = 1:z
        move!(r,West)
    end
    while y > 0
        if !isborder(r,Nord)
            move!(r,Nord)
            y -= 1
        else
            Flag, y = Obhod3(r,0,y)
        end
    end

end
function ugl(r)
    Flag = true
    y = 0
    z = 0
    while !isborder(r,Sud) || Flag
        if !isborder(r,Sud)
            move!(r,Sud)
            y += 1
        end
        if isborder(r,Sud)
            Flag,y =  Obhod2(r,2,y)
        end
    end
    while !isborder(r,Ost)
        move!(r,Ost)
        z += 1
    end
    return y,z
end
function Obhod2(r,HS,y)

    x = 0
    Flag = true
    while isborder(r,HorizonSide(HS)) && !isborder(r,HorizonSide((HS + 3) % 4))
        move!(r,HorizonSide((HS + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 3) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(HS))
        move!(r,HorizonSide(HS))
        y+=1
    end
    while !isborder(r,HorizonSide(HS)) && isborder(r,HorizonSide((HS + 1) % 4))
        move!(r,HorizonSide(HS))
        y+= 1
    end
    for i = 1:x
        move!(r,HorizonSide((HS+1) % 4))
    end
    return Flag,y
end

function markers(r)
    y = 0
    z = 0
    y, z = ugl(r)
    putmarker!(r)
    for HS = 0:3
        while !isborder(r,HorizonSide(HS))
        move!(r,HorizonSide(HS))
        putmarker!(r)
        end
    end
    Vozvrat(r,y,z)
end

markers(r)


text:
using TimerOutputs

@params struct LinearElasticityResult
    comp::Any
    u::Any
end
function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::LinearElasticityResult)
    return println("TopOpt linear elasticity result")
end

function simulate(
    problem::StiffnessTopOptProblem,
    topology=ones(getncells(TopOptProblems.getdh(problem).grid));
    round=true,
    hard=true,
    xmin=0.001,
)
    if round
        if hard
            solver = FEASolver(Direct, problem; xmin=0.0)
        else
            solver = FEASolver(Direct, problem; xmin=xmin)
        end
    else
        solver = FEASolver(Direct, problem; xmin=xmin)
    end
    vars = solver.vars
    fill_vars!(vars, topology; round=round)
    solver(Val{true})
    comp = dot(solver.u, solver.globalinfo.f)
    return LinearElasticityResult(comp, copy(solver.u))
end

function fill_vars!(vars::Array, topology; round)
    if round
        vars .= Base.round.(topology)
    else
        copyto!(vars, topology)
    end
    return vars
end


text:
# Thanks to load order issues, we have to do the documentation seperately

const README, compatability =
    let readme=read("$(dirname(dirname(@__DIR__)))/README.org", String)
        prolouge=match(r"\n([^#\n][\s\S]*)\n\* Progress", readme).captures[1]
        progress=match(r"\n\* Progress\s*((?:\|[^\n]+\n)+)", readme).captures[1]
        Base.Docs.doc!(@__MODULE__,
                       Base.Docs.Binding(@__MODULE__, Symbol(@__MODULE__)),
                       Base.Docs.docstr(parse(Org, prolouge),
                                        Dict(:path => joinpath(@__DIR__, @__FILE__),
                                             :linenumber => @__LINE__,
                                             :module => @__MODULE__)),
                       Union{})
        (parse(Org, readme),
         parse(Org,
               replace(replace(replace(progress,
                                       "| X " => "| =✓="),
                               r"\| +\|" => "| ~⋅~ |"),
                       r"\| +\|" => "| ~⋅~ |")))
    end

@doc org"""
#+begin_src julia
orgmatcher(::Type{C}) where {C <: OrgComponent}
#+end_src

Return a /matcher/ for components of type ~C~.
This will either be:
+ a regular expression which matcher the entire component
+ a function which takes a string and returns either
  - nothing, if the string does not start with the component
  - the substring which has been identified as an instance of the component
  - a tuple of the substring instance of the component, and the component data structure
""" matcher

@doc org"""
#+begin_src julia
consume(component::Type{<:OrgComponent}, text::AbstractString)
#+end_src
Try to /consume/ a ~component~ from the start of ~text~.

Returns a tuple of the consumed text and the resulting component
or =nothing= if this is not possible.
""" consume


text:
# Configuration for running spatial simulation
export simtype
# simtype=2
@everywhere simtype = 2
const num_trials = 1000
@everywhere const N_list = [6400,3200,1600,800,400,200,100,50,25,10,5,2]        # population size
#N_mut_list = [0.125, 0.25, 0.5, 1.0]
mutation_stddev_list = [0.002, 0.005,0.01,0.02,0.04]
const num_subpops = 1                    # Number of subpopulations
const num_attributes_list = [1,5]        # number attributes for quantitative representation
const ngens = 1				# Generations after burn-in
const burn_in= 800    # if an Int, generations of burn_in. If a Float, int_burn_in = burn_in*N+50
const ideal = 1.0
const neutral = true
const additive_error=false
const wrap_attributes=false

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
module BioCCP

using Base: Integer
using Distributions
export expectation_minsamplesize, std_minsamplesize, success_probability,
        expectation_fraction_collected, prob_occurrence_module

"""
    exp_ccdf(n, T; p=ones(n)/n, m=1, r=1, normalize=true)

Calculates `1 - F(t)`, which is the complement of the success probability
`F(t)=P(T ≤ t)` (= probability that the expected minimum
  number of designs T is smaller than `t` in order to
  see each module at least `m` times). This function
  serves as the integrand for calculating `E[T]`.

- `n`: number of modules in the design space
- `p`: vector with the probabilities/abundances of the different modules in the design space during library generation
- `T`: number of designs
- `m`: number of times each module has to observed in the sampled set of designs
- `r`: number of modules per design
- normalize: if true, normalize `p`

References:
- Doumas, A. V., & Papanicolaou, V. G. (2016). The coupon collector’s problem revisited:
        generalizing the double Dixie cup problem of Newman and Shepp. ESAIM: Probability and Statistics, 20, 367-399.
- Boneh, A., & Hofri, M. (1997). The coupon-collector problem revisited—a survey of
        engineering problems and computational methods. Stochastic Models, 13(1), 39-66.

## Examples

```julia-repl
julia> n = 100
julia> t = 500
julia> exp_ccdf(n, t; p=ones(n)/n, m=1, r=1, normalize=true)
0.4913906004535237
```
"""
function exp_ccdf(n, t; p=ones(n)/n, m=1, r=1, normalize=true)
    @assert length(p) == n

    # Normalize probabilities
    if normalize
        p=p ./ sum(p)
    end
    # Initialize probability P
    P_cdf = 1
    for i in 1:n
          Sm = 0
        for j in 1:m
            Sm += ((p[i]*r*t)^(j-1))/factorial(big(j-1)) #formulas see paper reference [1]
        end
        P_cdf *= (1 - Sm*exp(-p[i]*r*t))
    end
    P = 1 - P_cdf
    return P
end

"""
    approximate_moment(n, fun; p=ones(n)/n, q=1, m=1, r=1,
steps=1000, normalize=true)

Calculates the q-th rising moment of `T[N]` (number of designs that are needed to collect
all modules `m` times). Integral is approximated by the Riemann sum.

Reference:
- Doumas, A. V., & Papanicolaou, V. G. (2016). The coupon collector’s problem revisited:
        generalizing the double Dixie cup problem of Newman and Shepp. ESAIM: Probability
        and Statistics, 20, 367-399.

## Examples

```julia-repl
julia> n = 100
julia> fun = exp_ccdf
julia> approximate_moment(n, fun; p=ones(n)/n, q=1, m=1, r=1,
steps=10000, normalize=true)
518.8175339489885
```
"""
function approximate_moment(n, fun; p=ones(n)/n, q=1, m=1, r=1,
	        steps=500, normalize=true)
    @assert length(p) == n
    a = 0; b = n*log(n)
    ϵ = 0.001 # error tolerance
    while fun(n, b; p=p, m=m, r=r, normalize=normalize) > ϵ
        b += n
    end

    # integration exp_ccdf, see paper References [1]
    # build in adaptive integration (exp_ccdf is a very steep function): minimize function evaluation at constant function value, only evaluate function at steep part
    a = deepcopy(b)
    while fun(n, a; p=p, m=m, r=r, normalize=normalize) < 1 - ϵ
	a += -n/10
    end
    δ = (b-a)/steps; t = a:δ:b
    qth_moment = q * sum(δ .* 1 .* (0:δ:a-δ).^[q-1])  + q * sum(δ .* fun.(n, t; p=p, m=m, r=r, normalize=normalize) .* t.^[q-1])
    return qth_moment
end

"""
    expectation_minsamplesize(n; p=ones(n)/n, m=1, r=1, normalize=true)

Calculates the expected minimum number of designs  `E[T]` to observe each module at least `m` times.

- `n`: number of modules in the design space
- `p`: vector with the probabilities or abundances of the different modules
- `m`: number of times each module has to be observed in the sampled set of designs
- `r`: number of modules per design
- normalize: if true, normalize `p`

References:
- Doumas, A. V., & Papanicolaou, V. G. (2016). The coupon collector’s problem revisited:
    generalizing the double Dixie cup problem of Newman and Shepp. ESAIM: Probability
    and Statistics, 20, 367-399.
- Boneh, A., & Hofri, M. (1997). The coupon-collector problem revisited—a survey of
    engineering problems and computational methods. Stochastic Models, 13(1), 39-66.

## Examples

```julia-repl
julia> n = 100
julia> expectation_minsamplesize(n; p=ones(n)/n, m=1, r=1, normalize=true)
519
```
"""
function expectation_minsamplesize(n::Integer; p=ones(n)/n, m::Integer=1, r=1, normalize=true)
    @assert length(p) == n
    @assert n > 0
    @assert all(p .>= 0)
    @assert m > 0
    @assert r > 0
    E = approximate_moment(n, exp_ccdf; p=p, q=1, m=m, r=r, normalize=normalize)
    return Int(ceil(E))
end

"""
    std_minsamplesize(n::Integer; p=ones(n)/n, m::Integer=1, r=1, normalize=true)

Calculates the standard deviation on the minimum number of designs to observe each module at least `m` times.

- `n`: number of modules in the design space
- `p`: vector with the probabilities or abundances of the different modules
- `m`: number of complete sets of modules that need to be collected
- `r`: number of modules per design
- normalize: if true, normalize `p`

## Examples

```julia-repl
julia> n = 100
julia> std_minsamplesize(n; p=ones(n)/n, m=1, r=1, normalize=true)
126
```
"""
function std_minsamplesize(n::Integer; p=ones(n)/n, m::Integer=1, r=1, normalize=true)
    @assert length(p) == n
    @assert n > 0
    @assert all(p .>= 0)
    @assert m > 0
    @assert r > 0
    M1 = approximate_moment(n, exp_ccdf; p=p, q=1, m=m, r=r,  normalize=normalize)
    M2 = approximate_moment(n, exp_ccdf; p=p, q=2, m=m, r=r, normalize=normalize)
    var = M2 - M1 - M1^2
    return Int(ceil(sqrt(var)))
end

"""
    success_probability(n::Integer, t::Integer; p=ones(n)/n, m::Integer=1, r=1, normalize=true)

Calculates the success probability `F(t) = P(T ≤ t)` or the probability that
the minimum number of designs `T` to see each module at least `m` times
is smaller than `t`.

- `n`: number of modules in design space
- `t`: sample size/number of designs for which to calculate the success probability
- `p`: vector with the probabilities or abundances of the different modules
- `m`: number of complete sets of modules that need to be collected
- `r`: number of modules per design
- normalize: if true, normalize `p`

References:
- Boneh, A., & Hofri, M. (1997). The coupon-collector problem revisited—a survey of engineering problems and computational methods. Stochastic Models, 13(1), 39-66.

## Examples

```julia-repl
julia> n = 100
julia> t = 600
julia> success_probability(n, t; p=ones(n)/n, m=1, r=1, normalize=true)
0.7802171997092149
```
"""
function success_probability(n::Integer, t::Integer; p=ones(n)/n, m::Integer=1, r=1, normalize=true)
    @assert length(p) == n
    @assert n > 0
    @assert all(p .>= 0)
    @assert t >= 0
    @assert m > 0
    @assert r > 0
    P_success = 1 - exp_ccdf(n, t; p=p, m=m, r=r, normalize=normalize)
    return P_success
end

"""
    expectation_fraction_collected(n::Integer, t::Integer; p=ones(n)/n, r=1, normalize=true)

Calculates the fraction of all modules that is expected to be observed
after collecting `t` designs.

- `n`: number of modules in design space
- `t`: sample size/number of designs
- `p`: vector with the probabilities or abundances of the different modules
- `r`: number of modules per design
- normalize: if true, normalize `p`

References:
- Boneh, A., & Hofri, M. (1997). The coupon-collector problem revisited—a survey of engineering problems and computational methods. Stochastic Models, 13(1), 39-66.

## Examples

```julia-repl
julia> n = 100
julia> t = 200
julia> expectation_fraction_collected(n, t; p=ones(n)/n, r=1, normalize=true)
0.8660203251420364
```
"""
function expectation_fraction_collected(n::Integer, t::Integer; p=ones(n)/n, r=1, normalize=true)
    @assert length(p) == n
    @assert n > 0
    @assert all(p .>= 0)
    @assert t >= 0
    @assert r > 0
    if normalize
        p = p./sum(p)
    end
    frac = sum( (1-(1-p[i])^(t*r)) for i in 1:n )/n
    return frac
end

"""
    prob_occurrence_module(pᵢ, t::Integer, r, k::Integer)

Calculates probability that specific module with module probability `pᵢ`
has occurred `k` times after collecting `t` designs.

Sampling processes of individual modules are assumed to be independent Poisson processes.

- `pᵢ`: module probability
- `t`: sample size/number of designs
- `k`: number of occurrence

References:
- Boneh, A., & Hofri, M. (1997). The coupon-collector problem revisited—a survey of engineering problems and computational methods. Stochastic Models, 13(1), 39-66.

## Examples

```julia-repl
julia> pᵢ = 0.005
julia> t = 500
julia> k = 2
julia> r = 1
julia> prob_occurrence_module(pᵢ, t, r, k)
0.25651562069968376
```
"""
function prob_occurrence_module(pᵢ, t::Integer, r, k::Integer)
    @assert pᵢ > 0 && pᵢ <= 1
    @assert t >= 0
    @assert r >= 0
    @assert k >= 0
    poisson = Poisson(pᵢ * t * r)
	return pdf(poisson, k)
end

end

max_forks_repo_name: kirstvh/BioCCP


text:
# This file is a part of ArraysOfArrays.jl, licensed under the MIT License (MIT).

using ArraysOfArrays
using Test

using ElasticArrays
using UnsafeArrays
using Adapt

using Statistics
using StatsBase: cov2cor

@testset "array_of_similar_arrays" begin
    function rand_flat_array(Val_N::Val{N}) where {N}
        sz_max = (2,3,2,4,5)
        sz = ntuple(i -> sz_max[i], Val_N)
        x = rand(sz...)
    end

    function rand_nested_similar_arrays(Val_M::Val{M}, Val_N::Val{N}) where {M,N}
        sz_max = (2,3,2,4,5)
        sz_inner = ntuple(i -> sz_max[i], Val_M)
        sz_outer = ntuple(i -> sz_max[i + M], Val_N)

        A = Array{Array{Float64, M}, N}(undef, sz_outer...)
        for i in eachindex(A)
            A[i] = rand(sz_inner...)
        end
        A
    end

    @inline function test_from_flat(::Type{TT}, ::Type{RT}, Val_L::Val{L}) where {
        TT<:ArrayOfSimilarArrays,
        RT<:ArrayOfSimilarArrays,
        L
    }
        @testset "$TT from Array{Float64,$L}" begin
            A = rand_flat_array(Val_L)
            @test typeof(@inferred TT(A)) == RT
            @test typeof(@inferred convert(TT, A)) == RT
        end
    end

    @inline function test_from_nested(::Type{TT}, ::Type{RT}, Val_M::Val{M}, Val_N::Val{N}) where {
        TT<:ArrayOfSimilarArrays,
        RT<:ArrayOfSimilarArrays,
        M, N
    }
        @testset "$TT from Array{Array{Float64,$M},$N}" begin
            A = rand_nested_similar_arrays(Val_M, Val_N)

            A2_ctor = @inferred TT(A)
            U = eltype(flatview(A2_ctor))

            A_U = Array{Array{U,M},N}(A)

            @test typeof(A2_ctor) == RT
            @test A2_ctor == A_U

            A2_conv = @inferred convert(TT, A)
            @test typeof(A2_conv) == RT
            @test A2_conv == A2_ctor

            U = eltype(flatview(A2_ctor))
            A3 = @inferred Array(A2_ctor)
            @test typeof(A3) == Array{Array{U,M},N}
            @test A3 == A_U
        end
    end

    @testset "construct/convert from flat array" begin
        test_from_flat(ArrayOfSimilarArrays{Float64,2,3}, ArrayOfSimilarArrays{Float64,2,3,5,Array{Float64,5}}, Val(5))
        test_from_flat(ArrayOfSimilarArrays{Float64,2}, ArrayOfSimilarArrays{Float64,2,3,5,Array{Float64,5}}, Val(5))
        test_from_flat(ArrayOfSimilarVectors{Float64}, ArrayOfSimilarVectors{Float64,2,3,Array{Float64,3}}, Val(3))
        test_from_flat(VectorOfSimilarArrays{Float64}, VectorOfSimilarArrays{Float64,2,3,Array{Float64,3}}, Val(3))
        test_from_flat(VectorOfSimilarVectors{Float64}, VectorOfSimilarVectors{Float64,Array{Float64,2}}, Val(2))
        test_from_flat(VectorOfSimilarVectors, VectorOfSimilarVectors{Float64,Array{Float64,2}}, Val(2))

        test_from_flat(ArrayOfSimilarArrays{Float32,2,3}, ArrayOfSimilarArrays{Float32,2,3,5,Array{Float32,5}}, Val(5))
        test_from_flat(ArrayOfSimilarArrays{Float32,2}, ArrayOfSimilarArrays{Float32,2,3,5,Array{Float32,5}}, Val(5))
        test_from_flat(ArrayOfSimilarVectors{Float32}, ArrayOfSimilarVectors{Float32,2,3,Array{Float32,3}}, Val(3))
        test_from_flat(VectorOfSimilarArrays{Float32}, VectorOfSimilarArrays{Float32,2,3,Array{Float32,3}}, Val(3))
        test_from_flat(VectorOfSimilarVectors{Float32}, VectorOfSimilarVectors{Float32,Array{Float32,2}}, Val(2))
        test_from_flat(VectorOfSimilarVectors{Float32}, VectorOfSimilarVectors{Float32,Array{Float32,2}}, Val(2))
    end

    @testset "construct/convert from nested arrays" begin
        test_from_nested(ArrayOfSimilarArrays, ArrayOfSimilarArrays{Float64,2,3,5,Array{Float64,5}}, Val(2), Val(3))
        test_from_nested(ArrayOfSimilarArrays{Float64,2,3}, ArrayOfSimilarArrays{Float64,2,3,5,Array{Float64,5}}, Val(2), Val(3))
        test_from_nested(ArrayOfSimilarArrays{Float64}, ArrayOfSimilarArrays{Float64,2,3,5,Array{Float64,5}}, Val(2), Val(3))
        test_from_nested(ArrayOfSimilarArrays{Float32}, ArrayOfSimilarArrays{Float32,2,3,5,Array{Float32,5}}, Val(2), Val(3))

        test_from_nested(ArrayOfSimilarArrays, VectorOfSimilarArrays{Float64,4,5,Array{Float64,5}}, Val(4), Val(1))
        test_from_nested(ArrayOfSimilarArrays, ArrayOfSimilarVectors{Float64,4,5,Array{Float64,5}}, Val(1), Val(4))
        test_from_nested(ArrayOfSimilarArrays, VectorOfSimilarVectors{Float64,Array{Float64,2}}, Val(1), Val(1))

        test_from_nested(VectorOfSimilarArrays{Float64,4}, VectorOfSimilarArrays{Float64,4,5,Array{Float64,5}}, Val(4), Val(1))
        test_from_nested(VectorOfSimilarArrays{Float64}, VectorOfSimilarArrays{Float64,4,5,Array{Float64,5}}, Val(4), Val(1))
        test_from_nested(VectorOfSimilarArrays{Float32}, VectorOfSimilarArrays{Float32,4,5,Array{Float32,5}}, Val(4), Val(1))
        test_from_nested(VectorOfSimilarArrays, VectorOfSimilarArrays{Float64,4,5,Array{Float64,5}}, Val(4), Val(1))

        test_from_nested(ArrayOfSimilarVectors{Float64,4}, ArrayOfSimilarVectors{Float64,4,5,Array{Float64,5}}, Val(1), Val(4))
        test_from_nested(ArrayOfSimilarVectors{Float64}, ArrayOfSimilarVectors{Float64,4,5,Array{Float64,5}}, Val(1), Val(4))
        test_from_nested(ArrayOfSimilarVectors{Float32}, ArrayOfSimilarVectors{Float32,4,5,Array{Float32,5}}, Val(1), Val(4))
        test_from_nested(ArrayOfSimilarVectors, ArrayOfSimilarVectors{Float64,4,5,Array{Float64,5}}, Val(1), Val(4))

        test_from_nested(VectorOfSimilarVectors{Float64}, VectorOfSimilarVectors{Float64,Array{Float64,2}}, Val(1), Val(1))
        test_from_nested(VectorOfSimilarVectors{Float32}, VectorOfSimilarVectors{Float32,Array{Float32,2}}, Val(1), Val(1))
        test_from_nested(VectorOfSimilarVectors, VectorOfSimilarVectors{Float64,Array{Float64,2}}, Val(1), Val(1))

        r = @inferred(rand(5,5))
        @test @inferred(flatview(ArrayOfSimilarVectors(r))) == r
    end

    @testset "add remove" begin
        EA_ref1 = rand_flat_array(Val(3))
        EA_ref2 = rand_flat_array(Val(3))
        EA1 = ElasticArray{Float64, 3}(EA_ref1)
        EA2 = ElasticArray{Float64, 3}(EA_ref2)
        AEA1 = ArrayOfSimilarArrays{Float64, 3}(EA1)
        AEA2 = ArrayOfSimilarArrays{Float64, 3}(EA2)
        AEA1_ref = copy(AEA1)
        AEA2_ref = copy(AEA2)
        append!(AEA1, AEA2)
        prepend!(AEA2, AEA1_ref)
        N = size(AEA1.data)[end]
        for i in 1:N
            @test @inferred(reverse(AEA2.data, dims=3)[:,:,i]) == @inferred(AEA1.data[:,:,N+1-i])
        end

        A1 = ArrayOfSimilarArrays{Float64,1}(rand_flat_array(Val(1)))
        A2 = ArrayOfSimilarArrays{Float64,1}(rand_flat_array(Val(1)))
        A1_data = copy(A1.data)
        A2_data = copy(A2.data)
        append!(A1, A2)
        @test A1.data == vcat(A1_data, A2_data)
        prepend!(A1, A1)
        len = @inferred(length(A1.data))
        @test A1.data[1:Int(len/2)] == A1.data[Int(len/2 + 1):end]
    end

    @testset "similar and copyto!" begin
        A = ArrayOfSimilarArrays{Float64,1}(rand_flat_array(Val(1)))
        @test (@inferred copyto!((@inferred similar(A)), A)) == A

        A = ArrayOfSimilarArrays{Float64,2}(rand_flat_array(Val(5)))
        @test (@inferred copyto!((@inferred similar(A)), A)) == A

        A_data = rand_flat_array(Val(4))
        A = ArrayOfSimilarArrays{Float64, 2}(A_data)
        A_similar = similar(A, Array{Float64, 2}, size(A))
        @test @inferred(size(A)) == @inferred(size(A_similar))
        @test @inferred(size(A.data)) == @inferred(size(A_similar.data))
        @test typeof(A_similar.data) == typeof(A_data)
        @test typeof(A_similar) == typeof(A)
    end

    @testset "adapt" begin
        A_flat = rand(2,3,4,5,6)
        A_nested = nestedview(A_flat, 2)
        @test @inferred(adapt(identity, A_nested)) == A_nested
        @test typeof(adapt(identity, A_nested)) == typeof(A_nested)
    end

    @testset "deepcopy" begin
        A = ArrayOfSimilarArrays{Float64,1}(rand_flat_array(Val(1)))
        @test (@inferred deepcopy(A)) == A
        @test typeof(deepcopy(A)) == typeof(A)
    end

    @testset "flatview" begin
        A = rand_nested_similar_arrays(Val(3), Val(2))
        B = ArrayOfSimilarArrays(A)
        @inferred(flatview(B))[:] == collect(flatview(A))
    end

    @testset "deepgetindex" begin
        A = rand_nested_similar_arrays(Val(3), Val(2))
        B = ArrayOfSimilarArrays(A)

        @test deepgetindex(A, 3, 4, 2, 1, 2) == @inferred deepgetindex(B, 3, 4, 2, 1, 2)
        @test deepgetindex(A, 2:3, 4, 2, 1, 2) == @inferred deepgetindex(B, 2:3, 4, 2, 1, 2)
        @test deepgetindex(A, 2:3, 2:4, 2, 1, 2) == @inferred deepgetindex(B, 2:3, 2:4, 2, 1, 2)
        @test deepgetindex(A, 2, 4, :, 1, 2) == @inferred deepgetindex(B, 2, 4, :, 1, 2)
        @test deepgetindex(A, 2, 4, :, 1, 1:2) == @inferred deepgetindex(B, 2, 4, :, 1, 1:2)
        @test deepgetindex(A, 2:3, 4, :, 1, 2) == @inferred deepgetindex(B, 2:3, 4, :, 1, 2)
        @test deepgetindex(A, 2:3, 4, :, 1, 1:2) == @inferred deepgetindex(B, 2:3, 4, :, 1, 1:2)
    end

    @testset "deepsetindex!" begin
        function testdata()
            A = rand_nested_similar_arrays(Val(3), Val(2))
            B = ArrayOfSimilarArrays(A)
            A, B
        end

        A, B = testdata()
        @test deepsetindex!(A, 42, 3, 4, 2, 1, 2) == @inferred deepsetindex!(B, 42, 3, 4, 2, 1, 2)
        @test deepgetindex(B, 3, 4, 2, 1, 2) == 42

        A, B = testdata()
        X1 = rand(2)
        @test deepsetindex!(A, X1, 2:3, 4, 2, 1, 2) == @inferred deepsetindex!(B, X1, 2:3, 4, 2, 1, 2)
        @test deepgetindex(B, 2:3, 4, 2, 1, 2) == X1

        A, B = testdata()
        X2 = rand(2,2)
        @test deepsetindex!(A, X2, 2, 4, :, 1, 1:2) == @inferred deepsetindex!(B, X2, 2, 4, :, 1, 1:2)
        @test deepgetindex(B, 2, 4, :, 1, 1:2) == X2

        A, B = testdata()
        X3 = [rand(2,2), rand(2,2)]
        @test deepsetindex!(A, X3, 2:3, 4, :, 1, 1:2) == @inferred deepsetindex!(B, X3, 2:3, 4, :, 1, 1:2)
        @test deepgetindex(B, 2:3, 4, :, 1, 1:2) == X3
    end

    @testset "deepview" begin
        A = rand_nested_similar_arrays(Val(3), Val(2))
        B = ArrayOfSimilarArrays(A)

        @test deepview(A, 3, 4, 2, 1, 2) == @inferred deepview(B, 3, 4, 2, 1, 2)
        @test deepgetindex(A, 2:3, 4, 2, 1, 2) == @inferred deepview(B, 2:3, 4, 2, 1, 2)
        @test deepgetindex(A, 2:3, 2:4, 2, 1, 2) == @inferred deepview(B, 2:3, 2:4, 2, 1, 2)
        @test deepview(A, 2, 4, :, 1, 2) == @inferred deepview(B, 2, 4, :, 1, 2)
        @test deepview(A, 2, 4, :, 1, 1:2) == @inferred deepview(B, 2, 4, :, 1, 1:2)
        @test deepview(A, 2:3, 4, :, 1, 2) == @inferred deepview(B, 2:3, 4, :, 1, 2)
        @test deepview(A, 2:3, 4, :, 1, 1:2) == @inferred deepview(B, 2:3, 4, :, 1, 1:2)
    end

    @testset "empty" begin
        A = [rand(2,3), rand(2,3), rand(2,3)]
        B = ArrayOfSimilarArrays(A)
        @test typeof(@inferred empty(B)) == typeof(B)
        @test empty(A) == empty(B)
    end

    @testset "stats" begin
        VV = [rand(3) for i in 1:10]
        VV_aosa = ArrayOfSimilarArrays(VV)

        VA = [rand(2,3,3) for i in 1:10]
        VA_aosa = ArrayOfSimilarArrays(VA)

        array_cmp(A, B) = (A ≈ B) && (size(A) == size(B))

        function test_statistics_op(op::Function)
            @testset "$op" begin
                @test array_cmp(@inferred(op(VV_aosa)), op(VV))

                if op in (var, cov)
                    @test array_cmp(@inferred(op(VV_aosa, corrected = false)), op(VV, corrected = false))
                end

                if (op in (sum, mean, var))
                    @test array_cmp(@inferred(op(VA_aosa)), op(VA))
                end
            end
        end

        test_statistics_op(sum)
        test_statistics_op(mean)
        test_statistics_op(var)
        test_statistics_op(std)
        test_statistics_op(cov)

        @testset "cor" begin
            # Statistics.cor currently results in an error for Vector{Vector},
            # this should be considered a bug, though, since Statistics.cov
            # works fine.
            @test array_cmp(@inferred(cor(VV_aosa)), cov2cor(cov(VV), std(VV)))
        end
    end

    @testset "examples" begin
        A_flat = rand(2,3,4,5,6)
        A_nested = nestedview(A_flat, 2)

        @test A_nested isa AbstractArray{<:AbstractArray{T,2},3} where T
        @test flatview(A_nested) === A_flat

        A_flat = rand(4,4)
        A_nested = @inferred(nestedview(A_flat))

        @test A_nested.data == A_flat
        @test @inferred(size(A_nested))[1] == @inferred(size(A_flat))[1]
        @test @inferred(innersize(A_nested,1)) == @inferred(size(A_flat))[1]

        # -------------------------------------------------------------------

        A_flat = rand(2,3,4,5)
        ASA = @inferred(ArrayOfSimilarArrays{Float64,2,2}(A_flat))
        @test ASA.data == A_flat

        # -------------------------------------------------------------------
        A_nested = nestedview(ElasticArray{Float64}(undef, 2, 3, 0), 2)
        A_nested_copy = deepcopy(A_nested)

        for i in 1:4
            push!(A_nested, rand(2, 3))
        end
        @test size(flatview(A_nested)) == (2, 3, 4)

        resize!(A_nested, 6)
        @test size(flatview(A_nested)) == (2, 3, 6)

        for i in 1:4
            pushfirst!(A_nested, rand(2,3))
        end
        @test size(flatview(A_nested)) == (2, 3, 10)

        for i in 1:4
            pop!(A_nested)
        end
        @test size(flatview(A_nested)) == (2, 3, 6)

        for i in 1:size(A_nested)[1]
            pop!(A_nested)
        end
        @test_throws ArgumentError pop!(A_nested)

    end

    @testset "misc" begin
        N = 4
        r1 = rand(1,4); r2 = rand(1,4); r3 = rand(1,4); r4 = rand(1,4)
        r = vcat(r1,r2,r3,r4)
        VSV = VectorOfSimilarVectors(r)
        VSA = VectorOfSimilarArrays(r)
        ASA = ArrayOfSimilarArrays([r1,r2,r3,r4])

        f = x -> x.*2

        @test @inferred(IndexStyle(VSV)) == IndexLinear()
        @test @inferred(IndexStyle(VSA)) == IndexLinear()
        @test @inferred(IndexStyle(ASA)) == IndexLinear()
        @test VSA == VSV
        @test flatview(VSA) == flatview(VSV)

        @test @inferred(deepmap(f, ASA)).data == ASA.data.*2
        @test @inferred(innermap(f, ASA)).data == ASA.data.*2

        @test @inferred(ArraysOfArrays._innerlength(VSV)) == N
    end

    @testset "map and broadcast" begin
        A_flat = rand(2,3,4,5,6)
        A = nestedview(A_flat, 2)

        for do_map in (map, broadcast)
            @test @inferred(do_map(identity, A)) === A
        end
    end
end


text:
"""
    randomshuffle(ts::AbstractArray{T, 1} where T)

Generate a random constrained surrogate for `ts`. Destroys any linear
correlation in the signal, but preserves its amplitude distribution.

## Arguments
- **`ts`**: the time series for which to generate the surrogate realization.
"""
function randomshuffle(ts::AbstractArray{T, 1} where T)
    n = length(ts)
    ts[sample(1:n, n, replace = false)]
end


text:
include("base_extension_structs.jl")
include("single_qubit_gates.jl")
include("multi_qubit_gates.jl")


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at

#   http://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using ImageView
using Images
using FileIO
using Colors

# Run include("showimage.jl")

# 1. Read an image (could be replaced by one of the
# sample images)

# Test images
using TestImages

## XXX Should must load the raw images, should also load them directly into
##     a tuple or an array, using a loop.
scorpion1 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_1.png")))
scorpion2 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_2.png")))
scorpion3 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_3.png")))
scorpion4 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_4.png")))
scorpion5 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_5.png")))
scorpion6 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_6.png")))
scorpion7 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_7.png")))
scorpion8 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_8.png")))
scorpion9 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_9.png")))
scorpion10 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_10.png")))
scorpion11 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_11.png")))
scorpion12 = channelview(Gray.(load("data/scorpion-lapse/scorpion_lapse_12.png")))
scorpionlist = (scorpion1, scorpion2, scorpion3, scorpion4, scorpion5, scorpion6, scorpion7, scorpion8, scorpion9, scorpion10, scorpion11, scorpion12)

mandrilimg= testimage("mandrill")
lenaimg=testimage("lena")

testimgs = (mandrilimg, lenaimg)
testarrays = map(x -> channelview(Gray.(x)), testimgs)
mandg = channelview(Gray.(mandrilimg))

img = mandrilimg

imgg = Gray.(img)
fftimg = fft(channelview(imgg))
cimg = channelview(imgg)

zor = log.(abs.(fft(cimg)))
normfft = 255 .*(zor./maximum(zor))

#imshow(img)

# 2. Block algorithm: Make an algorithm that will pick out a block
#    centered around a pixel coordinate, do a calculation on that
#    block, and put the result into a matrix that has the size of the
#    image (could/should be interpreted as an image to evaluate
#    reasonableness).  Be able to adjust block size.

"""
Determine the blurriness of a picture using the blurriness
estimator from XXX.  Returns a number between 0 and 1.
0 means very sharp, 1 means very blurry.

*examples:

     julia> foo = blurriness(channelview(imgg))
     0.0693817138671875
"""
function blurriness(image)
    # Find M and N
    local   f = fft(image)
    # local  Fc = fft(shiftOrigin(image))
    local fc = f # Just taking a chance

    # local  AF = abs(Fc) # Absolute value of centered fourier tranformof image I
    local af = abs.(fc)

    # local  Mf = max(AF) # Maximum value of frequency component in F
    local mf = maximum(af)

    # local  TH = getNoOfPixelsAboveTreshold(Mf/1000)
    local th = length(find(x-> (x > mf/1000), af))
    # local  measure = TH / (M*N)
    local (m, n) = size(image)
    local measure = th / (m * n)

    # Fm = Frequency Domain Image Blur Measure
    return measure
end

"""
Apply the blurriness to a neigbourhood around every picture around a pixel in
the image, and use that estimated blurriness, as an estimate of the blurriness
of that particular pixel. Returns as an array

*Example

   tbd
"""
function blurrimap(img)

    ## XXX This thing does not correctly traverse all of the
    ##     pixels in the correct way so the blurriness estimate
    ##     is wrong.
    local (ydim, xdim) = size(img)
    local blocksize=5 # not all numbers works, that's bad.

    #  XXX the offset thing is not getting it right yet.
    local offset= blocksize

    local count=0
    local xlim=trunc(Int, (round(xdim/blocksize)*(blocksize - 1)))
    local ylim=trunc(Int, (round(ydim/blocksize)*(blocksize - 1)))

    local cview=channelview(img)
    local result = Array{Float64}(ydim, xdim)
    for y in 1:ylim, x in 1:xlim
        segment = cview[y:y+blocksize, x:x+blocksize]
        result[y + offset, x + offset] = blurriness(segment)
        count = count + 1
    end
    print("Number of segments in image  = $count\n")
    return result
end

function blurrimapstack(images)
    map(blurrimap, images)
end

blurritests = blurrimapstack(testarrays)

# TODO
# Find the index of the minimum blurriness value for every pixel, put that
# index in an array.   It will be a "depth" matrix.
# Then use this array to generate a new image consisting of the least blurry
# pixel for every position in the image stack.

bmm = blurrimap(imgg)

# 2.1 Make a sharpness estimator (based on the paper refered to
#     elsewhere) (done)

# 3. Sharpness stack generation: Make an algorithm that can process a
#    stack of images, and for each image, run the block algorithm and
#    generate a representation of that image, and store it alongside
#    the original.

# 4. Index picture generation: For the generated images, for each
#    pixel position, pick the index of the image that has the highest
#    value, put that value into a separate image.  This image will now
#    represent a "height matrix" into the focus-layer space.  It can
#    be interpreted as a conventional height matrix just for fun.
# 4.1 This algorithm could be extended into generating a 3D density
#     map generating algorithm, for picking out the sharp parts of the
#     images, assuming that these will outline actual features in a 3D
#     structure.

# 5. Generate a "focus stacked image" where original pixels are
#    collected, using the index picture to pick the pixel from the
#    sharpest image generating a "focused stacked" picture that
#    usually will be sharper than any of the picures being input
#    into the stacking algorithm.

typeof(img)


text:
function PolygonalFieldOfView(arg0::SphericalPolygonsSet, arg1::jdouble)
    return PolygonalFieldOfView((SphericalPolygonsSet, jdouble), arg0, arg1)
end

function PolygonalFieldOfView(arg0::Vector3D, arg1::PolygonalFieldOfView_DefiningConeType, arg2::Vector3D, arg3::jdouble, arg4::jint, arg5::jdouble)
    return PolygonalFieldOfView((Vector3D, PolygonalFieldOfView_DefiningConeType, Vector3D, jdouble, jint, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end


text:
"""
    coarsen(L::LT, f::FT, d::Tuple{IT,IT}) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Integer}

This function will *aggregate* a layer by merging cells. The function `f` passed
as its second argument is expected to take an array as input, and return a
single value of any type (but it is sort of a social contract here that this
will be a number of some sort, and if it is not, that you really know what
you're doing).

The size of the cells to aggregate is given by the tuple, so that `(2,2)` will
make groups of two cells vertically and two cells horizontally, for a total of
four cells. By default, the cells containing `nothing` will be *ignored*, so
that `f` is only applied to non-nothing values.

In cases where the number of cells to aggregate is not matching with the size of
the grid, and `ArgumentError` will be thrown. Note that there is no expectation
that the two values in `d` will be the same.
"""
function coarsen(L::LT, f::FT, d::Tuple{IT,IT}) where {LT <: SimpleSDMLayer, FT <: Function, IT <: Integer}
    cx, cy = d
    CX, CY = size(L)

    # Check that the arguments are compatible with the size of the grid
    mod(CX, cx) == 0 || throw(ArgumentError("The number of cells to aggregate ($(cx)) must be compatible with the number of rows ($(CX))"))
    mod(CY, cy) == 0 || throw(ArgumentError("The number of cells to aggregate ($(cy)) must be compatible with the number of columns ($(CY))"))

    # Then we create a new grid, full of undefined values of the same type as
    # the elements of L. This is not the best type here, of course, but I'm not
    # sure how to get a better idea _before_. This will all be converted to the
    # more compact type possible at the end, before returning a layer.
    nx = convert(Int64, CX/cx)
    ny = convert(Int64, CY/cy)
    newgrid = Array{Any}(undef, (nx, ny))

    # At this point, we do not need to create the new SimpleSDMLayer object, as
    # it can be a SimpleSDMPredictor which is not mutable. Instead, it is better
    # to work directly on the grid. We will iterate directly on the new grid,
    # and copy information from the old grid as needed.
    for i in 1:size(newgrid, 1)
        old_i = ((i-1)*cx+1):(i*cx)
        for j in 1:size(newgrid, 2)
            old_j = ((j-1)*cy+1):(j*cy)

            # This operation is super-duper important because it is the only
            # thing that makes the function work for some types - why it would
            # be the case is beyond me, but it turned out that being explicit
            # about converting things to a vector keep Julia happy, and my sense
            # of self worth is directly tied to how little error messages I get.
            V = convert(Vector, vec(L.grid[old_i, old_j]))

            # We remove the nothing values from the grid
            filter!(!isnothing, V)

            # If there is nothing left in V, we return nothing -- if not, we return f(V)
            newgrid[i,j] = length(V) == 0 ? nothing : f(V)

        end
    end

    # Now we change the type
    internal_types = unique(typeof.(newgrid))
    newgrid = convert(Matrix{Union{internal_types...}}, newgrid)

    # Now that everything is done, we can return an object of the correct type
    NT = LT <: SimpleSDMPredictor ? SimpleSDMPredictor : SimpleSDMResponse
    return NT(newgrid, L.left, L.right, L.bottom, L.top)
end


text:
# This file requires a GPU in order to run. It gets tested as part of the
# JuliaGPU CI on Gitlab. Contact @odow for more details.

# CUDA_jll is not in our Project.toml, so we need to install it on GITLAB_CI.
import Pkg
Pkg.add(Pkg.PackageSpec(name = "CUDA_jll", version = "10.1"))

using CUDA_jll  # CUDA_jll must be loaded _before_ SCS!
using SCS

using Test

@test SCS.GpuIndirectSolver in SCS.available_solvers

include("test_problems.jl")
feasible_basic_problems(SCS.GpuIndirectSolver)

# TODO(odow): consider re-enabling these
# include("MOI_wrapper.jl")
# moi_tests(SCS.GpuIndirectSolver)


text:
import Base: copy, getindex, sqrt, length, push!, setindex!

import ITensors:
  # types
  MPO,
  # circuits/gates.jl
  space,
  state,
  noise,
  dag

import LinearAlgebra: normalize!, tr, norm

import SCS

import Convex

import Optimisers.state as optimizerstate

import Observers: update!


text:
if VERSION < v"0.7.0"
    if Pkg.installed("DifferentialEquations") == nothing
        Pkg.add("DifferentialEquations")
    end
end

using DifferentialEquations


text:
mutable struct WebIOPlot
    p::PlotlyBase.Plot
    scope::Scope
end

Base.getindex(p::WebIOPlot, key) = p.scope[key] # look up Observables

WebIO.render(p::WebIOPlot) = WebIO.render(p.scope)
Base.show(io::IO, mm::MIME"text/plain", p::WebIOPlot) = show(io, mm, p.p)
Base.show(io::IO, mm::MIME"text/html", p::WebIOPlot) = show(io, mm, p.scope)

function WebIOPlot(
        args...;
        events::AbstractDict=Dict(),
        options::AbstractDict=Dict("showLink"=> false),
        kwargs...
    )
    # build plot, get json, setup options
    p = Plot(args...; kwargs...)
    lowered = JSON.lower(p)
    id = string("#plot-", p.divid)

    # setup scope
    deps = [
        "Plotly" => joinpath(@__DIR__, "..", "assets", "plotly-latest.min.js"),
        joinpath(@__DIR__, "..", "assets", "plotly_webio_bundle.js")
    ]
    scope = Scope(imports=deps)
    scope.dom = dom"div"(id=string("plot-", p.divid), events=events)

    # INPUT: Observables for plot events
    scope["svg"] = Observable("")
    scope["hover"] = Observable(Dict())
    scope["selected"] = Observable(Dict())
    scope["click"] = Observable(Dict())
    scope["relayout"] = Observable(Dict())

    # OUTPUT: setup an observable which sends modify commands
    scope["_commands"] = Observable{Any}([])

    # Do the respective action when _commands is triggered
    onjs(scope["_commands"], @js function (args)
        @var fn = args.shift()
        @var elem = this.plotElem
        @var Plotly = this.Plotly
        args.unshift(elem) # use div as first argument

        Plotly[fn].apply(this, args).then(function(gd)
            Plotly.toImage(elem, Dict("format" => "svg"))
        end).then(function(data)
            # TODO: make this optional
            @var svg_data = data.replace("data:image/svg+xml,", "")
            $(scope["svg"])[] = decodeURIComponent(svg_data)
        end)
    end)

    onimport(scope, JSExpr.@js function (Plotly)

        # set up container element
        @var gd = this.dom.querySelector($id);

        # save some vars for later
        this.plotElem = gd
        this.Plotly = Plotly
        if (window.Blink !== undefined)
            # set css style for auto-resize
            gd.style.width = "100%";
            gd.style.height = "100vh";
            gd.style.marginLeft = "0%";
            gd.style.marginTop = "0vh";
        end

        window.onresize = () -> Plotly.Plots.resize(gd)

        function save_svg()
            Plotly.toImage(gd, $(Dict("format" => "svg"))).then(function(data)
                @var svg_data = data.replace("data:image/svg+xml,", "")
                $(scope["svg"])[] = decodeURIComponent(svg_data)
            end)
        end

        # Draw plot in container
        Plotly.newPlot(
            gd, $(lowered[:data]), $(lowered[:layout]), $(options)
        ).then(save_svg);

        # hook into plotly events
        gd.on("plotly_hover", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "hover");
            if !(filtered_data.isnil)
                $(scope["hover"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_unhover", () -> $(scope["hover"])[] = Dict())

        gd.on("plotly_selected", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "selected");
            if !(filtered_data.isnil)
                $(scope["selected"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_deselect", () -> $(scope["selected"])[] = Dict())

        gd.on("plotly_relayout", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "relayout");
            if !(filtered_data.isnil)
                $(scope["relayout"])[] = filtered_data.out
            end
        end)

        gd.on("plotly_click", function (data)
            @var filtered_data = WebIO.CommandSets.Plotly.filterEventData(gd, data, "click");
            if !(filtered_data.isnil)
                $(scope["click"])[] = filtered_data.out
            end
        end)
    end)

    # create no-op `on` callback for svg so it is _always_ sent
    # to us
    on(scope["svg"]) do x end

    WebIOPlot(p, scope)
end


text:
using VegaLite

function relate_heat(data::DataFrame, method::Symbol)
    mth = String(method)
    data2 = select(data, :sample_2 => :sample_1, :sample_1 => :sample_2, :n_loci, :Lynch)
    append!(data2, data)
    sort!(data2, [:sample_1, :sample_2])
    transform!(data2, method => (i -> round.(i, digits = 4)) => method )
    tt_mthd = mth * ": '+datum."*mth
    data2 |> @vlplot(
        transform = [
            {calculate="''+datum.sample_1+ ' '+' '+datum.sample_2+'  '+' $tt_mthd", as="tt"}
        ],
        :rect,
        tooltip = "tt:n",
        :sample_1,
        :sample_2,
        color=method
    )
end

function relate_heat2(data::DataFrame, method::Symbol)
    mth = String(method)
    data2 = select(data, :sample_2 => :sample_1, :sample_1 => :sample_2, :n_loci, :Lynch)
    append!(data2, data)
    sort!(data2, [:sample_1, :sample_2])
    transform!(data2, method => (i -> round.(i, digits = 4)) => method )
    data2 |> @vlplot(
        hconcat=[
          {
            height = 700,
            width = :container,
            selection={
              brush={
                type="interval"
              }
            },
            encoding={
              x={
                field="sample_1",
                type="nominal"
              },
              color={
                field=mth
              },
              y={
                field="sample_2",
                type="nominal"
              }
            },
            mark="rect"
          },
          {
            hconcat=[
              {
                encoding={
                  text={
                    field="sample_1",
                    type="nominal"
                  },
                  y={
                    axis=nothing,
                    field="row_number",
                    type="ordinal"
                  }
                },
                title="Sample 1",
                width=50,
                mark="text"
              },
              {
                encoding={
                  text={
                    field="sample_2",
                    type="nominal"
                  },
                  y={
                    axis=nothing,
                    field="row_number",
                    type="ordinal"
                  }
                },
                title="Sample 2",
                width=50,
                mark="text"
              },
              {
                encoding={
                  text={
                    field=mth,
                    type="nominal"
                  },
                  y={
                    axis=nothing,
                    field="row_number",
                    type="ordinal"
                  }
                },
                title="Est",
                width=50,
                mark="text"
              }
            ],
            transform=[
              {
                filter={
                  selection="brush"
                }
              },
              {
                window=[
                  {
                    as="rank",
                    op="rank"
                  }
                ]
              },
              {
                filter={
                  field="rank",
                  lt=20
                }
              }
            ]
          }
        ],
        width=300,
        resolve={
          legend={
            color="independent"
          }
        },
        transform=[
          {
            window=[
              {
                as="row_number",
                op="row_number"
              }
            ]
          }
        ]
      )
end

c = relate_heat2(a, :Lynch)

function relate_hist(data::DataFrame, method::Symbol)
    mth = String(method)
    a |> @vlplot(
        width = :container,
        height = 350,
        :bar, #tooltip = "true",
        x={method, bin=true, step = 5, title = mth},
        y={"count()", title = "Count"}
    )
end

@vlplot(
  height=300,
  hconcat=[
    {
      selection={
        brush={
          type="interval"
        }
      },
      encoding={
        x={
          field="sample_1",
          type="nominal"
        },
        color={
          field="Lynch"
        },
        y={
          field="sample_2",
          type="nominal"
        }
      },
      mark="rect"
    },
    {
      hconcat=[
        {
          encoding={
            text={
              field="sample_1",
              type="nominal"
            },
            y={
              axis=nothing,
              field="row_number",
              type="ordinal"
            }
          },
          title="Sample 1",
          width=50,
          mark="text"
        },
        {
          encoding={
            text={
              field="sample_2",
              type="nominal"
            },
            y={
              axis=nothing,
              field="row_number",
              type="ordinal"
            }
          },
          title="Sample 2",
          width=50,
          mark="text"
        },
        {
          encoding={
            text={
              field="Lynch",
              type="nominal"
            },
            y={
              axis=nothing,
              field="row_number",
              type="ordinal"
            }
          },
          title="Est",
          width=50,
          mark="text"
        }
      ],
      transform=[
        {
          filter={
            selection="brush"
          }
        },
        {
          window=[
            {
              as="rank",
              op="rank"
            }
          ]
        },
        {
          filter={
            field="rank",
            lt=20
          }
        }
      ]
    }
  ],
  width=300,
  resolve={
    legend={
      color="independent"
    }
  },
  transform=[
    {
      window=[
        {
          as="row_number",
          op="row_number"
        }
      ]
    }
  ]
)


text:
############################################################################################
# struct
############################################################################################
"""
    SemFit

Fitted structural equation model.

# Interfaces
- `minimum(::SemFit)` -> minimum objective value
- `solution(::SemFit)` -> parameter estimates
- `start_val(::SemFit)` -> starting values
- `model(::SemFit)`
- `optimization_result(::SemFit)`

- `optimizer(::SemFit)` -> optimization algorithm
- `n_iterations(::SemFit)` -> number of iterations
- `convergence(::SemFit)` -> convergence properties
"""
mutable struct SemFit{Mi, So, St, Mo, O}
    minimum::Mi
    solution::So
    start_val::St
    model::Mo
    optimization_result::O
end

############################################################################################
# pretty printing
############################################################################################

function Base.show(io::IO, semfit::SemFit)
    print(io, "Fitted Structural Equation Model \n")
    print(io, "=============================================== \n")
    print(io, "--------------------- Model ------------------- \n")
    print(io, "\n")
    print(io, semfit.model)
    print(io, "\n")
    #print(io, "Objective value: $(round(semfit.minimum, digits = 4)) \n")
    print(io, "------------- Optimization result ------------- \n")
    print(io, "\n")
    print(io, semfit.optimization_result)
end

############################################################################################
# additional methods
############################################################################################

# access fields
minimum(sem_fit::SemFit) = sem_fit.minimum
solution(sem_fit::SemFit) = sem_fit.solution
start_val(sem_fit::SemFit) = sem_fit.start_val
model(sem_fit::SemFit) = sem_fit.model
optimization_result(sem_fit::SemFit) = sem_fit.optimization_result

# optimizer properties
optimizer(sem_fit::SemFit) = optimizer(optimization_result(sem_fit))
n_iterations(sem_fit::SemFit) = n_iterations(optimization_result(sem_fit))
convergence(sem_fit::SemFit) = convergence(optimization_result(sem_fit))


text:
using Printf
using LinearAlgebra
using Random

# These must be set manually :-(
include("../common/PWGrid_v02.jl")
#include("../common/PWGrid_v03.jl")

include("driver_harm.jl")

@time test_main( Ns_in=[30, 30, 30], solution_method="diag_lobpcg" )  # default
@time test_main( Ns_in=[30, 30, 30], solution_method="Emin" )
#@time test_main( ecutwfc_Ry=54.0 )


text:
abstract type DDFType end

abstract type AbstractDDF end

struct DDF{C <: DDFType,OVERDX} <: AbstractDDF end

"""
    DDF([ddftype=Fields.Roma],[dx=1.0])

Construct a discrete delta function operator. This is generally only needed
internally by the `Regularize` operator, so the user doesn't have much need
for accessing this directly. The default DDF is the `Roma` function, which
has a support of 3 grid cells. Other choices are the `Goza` operator, which is
a truncated Gaussian with 28 cells support, and the `Witchhat`, which has 2
cells support. The resulting operator is evaluated with one, two or three
coordinate arguments, producing, respectively, 1-d, 2-d, or 3-d smeared delta functions.
It can also be called with the usual Julia vectorized dot notation with arrays of
arguments.
The optional cell spacing argument `dx` rescales the coordinates by this spacing,
and the result is also rescaled by this spacing (raised to the number of dimensions).
This spacing argument defaults to 1.0.

```jldoctest
julia> ddf = DDF(ddftype=Fields.Roma)
Discrete delta function operator of type ViscousFlow.Fields.Roma, with spacing 1.0

julia> ddf(1)
0.16666666666666666

julia> ddf(-1)
0.16666666666666666

julia> ddf.([-1,0,1])
3-element Array{Float64,1}:
 0.16666666666666666
 0.6666666666666666
 0.16666666666666666
```
"""
function DDF(;dx::Real=1.0,ddftype=Yang3)
  DDF{ddftype,1.0/dx}()
end

struct GradDDF{C <: Fields.DDFType,OVERDX,D} <: AbstractDDF end

function GradDDF(dir::Int;dx::Real=1.0,ddftype=Yang3)
  GradDDF{ddftype,1.0/dx,dir}()
end

# This macro creates an abstract type in the DDFType system for the named ddftype
# It also creates functions to evaluate the ddf in 1-d,2-d or 3-d.
macro ddffunc(ddftype)
    fname = Symbol("ddf_",lowercase(string(ddftype)))
    dfname = Symbol("ddf_d",lowercase(string(ddftype)))
    return esc(quote
            abstract type $ddftype <: DDFType end
            @inline (::DDF{$ddftype,OVERDX})(x::T) where {OVERDX,T <: Real} =
                  $fname(abs(x)*OVERDX)*OVERDX
            @inline (::DDF{$ddftype,OVERDX})(x::T,y::T) where {OVERDX,T <: Real} =
                          $fname(abs(x)*OVERDX)*OVERDX*$fname(abs(y)*OVERDX)*OVERDX
            @inline (::DDF{$ddftype,OVERDX})(x::T,y::T,z::T) where {OVERDX,T <: Real} =
                          $fname(abs(x)*OVERDX)*OVERDX*$fname(abs(y)*OVERDX)*OVERDX*$fname(abs(z)*OVERDX)*OVERDX

            @inline (::GradDDF{$ddftype,OVERDX,1})(x::T) where {OVERDX,T <: Real} =
                   $dfname(x*OVERDX)*OVERDX*OVERDX

            @inline (::GradDDF{$ddftype,OVERDX,1})(x::T,y::T) where {OVERDX,T <: Real} =
                   $dfname(x*OVERDX)*OVERDX*OVERDX*$fname(abs(y)*OVERDX)*OVERDX
            @inline (::GradDDF{$ddftype,OVERDX,2})(x::T,y::T) where {OVERDX,T <: Real} =
                   $fname(abs(x)*OVERDX)*OVERDX*$dfname(y*OVERDX)*OVERDX*OVERDX

            @inline (::GradDDF{$ddftype,OVERDX,1})(x::T,y::T,z::T) where {OVERDX,T <: Real} =
                   $dfname(x*OVERDX)*OVERDX*OVERDX*$fname(abs(y)*OVERDX)*OVERDX*$fname(abs(z)*OVERDX)*OVERDX
            @inline (::GradDDF{$ddftype,OVERDX,2})(x::T,y::T,z::T) where {OVERDX,T <: Real} =
                   $fname(abs(x)*OVERDX)*OVERDX*$dfname(y*OVERDX)*OVERDX*OVERDX*$fname(abs(z)*OVERDX)*OVERDX
            @inline (::GradDDF{$ddftype,OVERDX,3})(x::T,y::T,z::T) where {OVERDX,T <: Real} =
                   $fname(abs(x)*OVERDX)*OVERDX*$fname(abs(y)*OVERDX)*OVERDX*$dfname(z*OVERDX)*OVERDX*OVERDX
            end)
end

#==== ROMA ====#

@ddffunc Roma

roma1(r) = (1+sqrt(-3r^2+1))/3
roma2(r) = (5-3r-sqrt(1-3(1-r)^2))/6

@inline ddf_roma(r::Real) = r > 1.5 ? 0.0 : r <= 0.5 ? roma1(r) : roma2(r)

droma1(r) = -r/sqrt(-3r^2+1)
droma2(r) = -0.5*(1+(1-r)/sqrt(1-3(1-r)^2))

@inline ddf_droma(r::Real) = abs(r) > 1.5 ? 0.0 : abs(r) <= 0.5 ? sign(r)*droma1(abs(r)) : sign(r)*droma2(abs(r))

#==== GOZA ====#

@ddffunc Goza

@inline ddf_goza(r::Real) = r > 14 ? 0.0 : exp(-π^2/36 * r^2)*sqrt(π/36)

@inline ddf_dgoza(r::Real) = abs(r) > 14 ? 0.0 : (-π^2/18)*r*exp(-π^2/36 * r^2)*sqrt(π/36)

#==== WITCHHAT ====#

@ddffunc Witchhat

@inline ddf_witchhat(r::Real) = r > 1 ? 0.0 : 1-r

@inline ddf_dwitchhat(r::Real) = abs(r) > 1 ? 0.0 : -sign(r)

#==== M3 ====#

@ddffunc M3

M31(r) = 0.75-r^2
M32(r) = 0.5*(1.5-r)^2
@inline ddf_m3(r::Real) = r > 1.5 ? 0.0 : r <= 0.5 ? M31(r) : M32(r)

dM31(r) = -2*r
dM32(r) = r-1.5
@inline ddf_dm3(r::Real) = abs(r) > 1.5 ? 0.0 : abs(r) <= 0.5 ? dM31(r) : sign(r)*dM32(abs(r))

#==== YANG3 ====#

@ddffunc Yang3

Yang31(r) = 17/48+sqrt(3)π/108+r/4-r^2/4+(1-2*r)/16*sqrt(-12*r^2+12*r+1)-sqrt(3)/12*asin(sqrt(3)/2*(2*r-1))
Yang32(r) = 55/48-sqrt(3)π/108-13*r/12+r^2/4+(2*r-3)/48*sqrt(-12*r^2+36*r-23)+sqrt(3)/36*asin(sqrt(3)/2*(2*r-3))

@inline ddf_yang3(r::Real) = r > 2.0 ? 0.0 : r <= 1.0 ? Yang31(r) : Yang32(r)

dYang31(r) =    1/4-r/2-sqrt(-12*r^2+12*r+1)/4
dYang32(r) = -13/12+r/2+sqrt(-12*r^2+36*r-23)/12

@inline ddf_dyang3(r::Real) = abs(r) > 2.0 ? 0.0 : abs(r) <= 1.0 ? sign(r)*dYang31(abs(r)) : sign(r)*dYang32(abs(r))

function Base.show(io::IO, ddf::DDF{C,OVERDX}) where {C<:DDFType,OVERDX}
    print(io, "Discrete delta function operator of type $C, with spacing $(1.0/OVERDX)")
end