text:
mutable struct RowType{T} <: NamedType
    type::String
    namespace::String
    name::String
    doc::Union{String, Nothing}
    aliases::Union{Vector{String}, Nothing}
    fields::Vector{FieldType}
    schema::T # Tables.Schema
end

==(x::RowType, y::RowType) = x.name == y.name && x.fields == y.fields
RowType(name, fields::Vector{FieldType}, sch::Tables.Schema, namespace="", doc=nothing, aliases=nothing) =
    RowType("record", namespace, name, doc, aliases, fields, sch)
RowType() = RowType("record", "", "", nothing, nothing, FieldType[])
StructTypes.StructType(::Type{<:RowType}) = StructTypes.Mutable()
StructTypes.excludes(::Type{<:RowType}) = (:schema,)

Tables.Schema(::Type{Record{names, types, N}}) where {names, types, N} = Tables.Schema{names, types}()

schematype(::Type{T}) where {T <: Tables.Schema} = schematype(T())
schematype(sch::Tables.Schema{names, types}, name=autoname(Record{names, types, fieldcount(types)})) where {names, types} =
    RowType(name, FieldType[FieldType(String(names[i]), schematype(fieldtype(types, i))) for i = 1:length(names)], sch)

mutable struct RowWriteClosure{B, S, T, KW}
    type::RowType{S}
    buf::T
    pos::Int
    len::Int
    opts::KW
end

function (f::RowWriteClosure{B, S, T, KW})(val, i, nm) where {B, S, T, KW}
    # @debug 2 "writing $typ at pos = $(f.pos)"
    f.pos = writevalue(B(), f.type.fields[i].type, val, f.buf, f.pos, f.len, f.opts)
end

function writevalue(B::Binary, T::RowType{S}, row, buf, pos, len, opts) where {S}
    c = RowWriteClosure{Binary, S, typeof(buf), typeof(opts)}(T, buf, pos, len, opts)
    Tables.eachcolumn(c, T.schema, row)
    return c.pos
end

mutable struct NBytesRowClosure{S}
    RT::RowType{S}
    n::Int
end

function (f::NBytesRowClosure{S})(val, i, nm) where {S}
    f.n += nbytes(f.RT.fields[i].type, val)
end

function nbytes(RT::RowType, row)
    c = NBytesRowClosure(RT, 0)
    Tables.eachcolumn(c, RT.schema, row)
    return c.n
end


text:
# This file contains functions to handle the symetries.
# The type SymOp is defined in Symop.jl

# A symmetry (W, w) (or (S, τ)) induces a symmetry in the Brillouin zone
# that the Hamiltonian at S k is unitary equivalent to that at k, which we exploit to
# reduce computations. The relationship is
# S = W'
# τ = -W^-1 w
# (valid both in reduced and cartesian coordinates). In our notation
# the rotation matrix W and translation w are such that, for each atom of
# type A at position a, W a + w is also an atom of type A.

# The full (reducible) Brillouin zone is implicitly represented by
# a set of (irreducible) kpoints (see explanation in docs). Each
# irreducible k-point k comes with a list of symmetry operations
# (S, τ) (containing at least the trivial operation (I, 0)), where S
# is a unitary matrix (/!\ in cartesian but not in reduced coordinates)
# and τ a translation vector. The k-point is then used to represent
# implicitly the information at all the kpoints Sk. The relationship
# between the Hamiltonians is
# H_{Sk} = U H_k U*, with
# (Uu)(x) = u(W x + w)
# or in Fourier space
# (Uu)(G) = e^{-i G τ} u(S^-1 G)
# In particular, we can choose the eigenvectors at Sk as u_{Sk} = U u_k

# We represent then the BZ as a set of irreducible points `kpoints`,
# and a set of weights `kweights` (summing to 1). The value of
# observables is given by a weighted sum over the irreducible kpoints,
# plus a symmetrization operation (which depends on the particular way
# the observable transforms under the symmetry).

# There is by decreasing cardinality
# - The group of symmetry operations of the lattice
# - The group of symmetry operations of the crystal (model.symmetries)
# - The group of symmetry operations of the crystal that preserves the BZ mesh (basis.symmetries)

# See https://juliamolsim.github.io/DFTK.jl/stable/advanced/symmetries for details.

@doc raw"""
Return the ``k``-point symmetry operations associated to a lattice and atoms.
"""
function symmetry_operations(lattice, atoms, positions, magnetic_moments=[];
        tol_symmetry=SYMMETRY_TOLERANCE)
    @assert length(atoms) == length(positions)
    atom_groups = [findall(Ref(pot) .== atoms) for pot in Set(atoms)]
    Ws, ws = spglib_get_symmetry(lattice, atom_groups, positions, magnetic_moments; tol_symmetry)
    [SymOp(W, w) for (W, w) in zip(Ws, ws)]
end

"""
Filter out the symmetry operations that respect the symmetries of the discrete BZ grid
"""
function symmetries_preserving_kgrid(symmetries, kcoords)
    kcoords_normalized = normalize_kpoint_coordinate.(kcoords)
    T = eltype(kcoords[1])
    atol = T <: Rational ? 0 : sqrt(eps(T))
    is_approx_in(x, X) = any(y -> isapprox(x, y; atol), X)
    function preserves_grid(symop)
        all(is_approx_in(normalize_kpoint_coordinate(symop.S * k), kcoords_normalized)
            for k in kcoords_normalized)
    end
    filter(preserves_grid, symmetries)
end

@doc raw"""
Apply various standardisations to a lattice and a list of atoms. It uses spglib to detect
symmetries (within `tol_symmetry`), then cleans up the lattice according to the symmetries
(unless `correct_symmetry` is `false`) and returns the resulting standard lattice
and atoms. If `primitive` is `true` (default) the primitive unit cell is returned, else
the conventional unit cell is returned.
"""
function standardize_atoms(lattice, atoms, positions, magnetic_moments=[]; kwargs...)
    @assert length(atoms) == length(positions)
    @assert isempty(magnetic_moments) || (length(atoms) == length(magnetic_moments))
    atom_groups = [findall(Ref(pot) .== atoms) for pot in Set(atoms)]
    ret = spglib_standardize_cell(lattice, atom_groups, positions, magnetic_moments; kwargs...)
    (; ret.lattice, atoms, ret.positions, ret.magnetic_moments)
end

"""
Apply a symmetry operation to eigenvectors `ψk` at a given `kpoint` to obtain an
equivalent point in [-0.5, 0.5)^3 and associated eigenvectors (expressed in the
basis of the new ``k``-point).
"""
function apply_symop(symop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
    S, τ = symop.S, symop.τ
    symop == one(SymOp) && return kpoint, ψk

    # Apply S and reduce coordinates to interval [-0.5, 0.5)
    # Doing this reduction is important because
    # of the finite kinetic energy basis cutoff
    @assert all(-0.5 .≤ kpoint.coordinate .< 0.5)
    Sk_raw = S * kpoint.coordinate
    Sk = normalize_kpoint_coordinate(Sk_raw)
    kshift = convert.(Int, Sk - Sk_raw)
    @assert all(-0.5 .≤ Sk .< 0.5)

    # Check whether the resulting k-point is in the basis:
    ikfull = findfirst(1:length(basis.kpoints)) do idx
        all(isinteger, basis.kpoints[idx].coordinate - Sk)
    end
    if isnothing(ikfull)
        # Build a new k-point datastructure:
        Skpoint = build_kpoints(basis, [Sk])[1]
    else
        Skpoint = basis.kpoints[ikfull]
        @assert Skpoint.coordinate ≈ Sk
    end

    # Since the eigenfunctions of the Hamiltonian at k and Sk satisfy
    #      u_{Sk}(x) = u_{k}(S^{-1} (x - τ))
    # their Fourier transform is related as
    #      ̂u_{Sk}(G) = e^{-i G \cdot τ} ̂u_k(S^{-1} G)
    invS = Mat3{Int}(inv(S))
    Gs_full = [G + kshift for G in G_vectors(basis, Skpoint)]
    ψSk = zero(ψk)
    for iband in 1:size(ψk, 2)
        for (ig, G_full) in enumerate(Gs_full)
            igired = index_G_vectors(basis, kpoint, invS * G_full)
            @assert igired !== nothing
            ψSk[ig, iband] = cis2pi(-dot(G_full, τ)) * ψk[igired, iband]
        end
    end

    Skpoint, ψSk
end

"""
Apply a symmetry operation to a density.
"""
function apply_symop(symop::SymOp, basis, ρin)
    symop == one(SymOp) && return ρin
    symmetrize_ρ(basis, ρin; symmetries=[symop])
end

# Accumulates the symmetrized versions of the density ρin into ρout (in Fourier space).
# No normalization is performed
function accumulate_over_symmetries!(ρaccu, ρin, basis, symmetries)
    T = eltype(basis)
    for symop in symmetries
        # Common special case, where ρin does not need to be processed
        if symop == one(SymOp)
            ρaccu .+= ρin
            continue
        end

        # Transform ρin -> to the partial density at S * k.
        #
        # Since the eigenfunctions of the Hamiltonian at k and Sk satisfy
        #      u_{Sk}(x) = u_{k}(S^{-1} (x - τ))
        # with Fourier transform
        #      ̂u_{Sk}(G) = e^{-i G \cdot τ} ̂u_k(S^{-1} G)
        # equivalently
        #     ρ ̂_{Sk}(G) = e^{-i G \cdot τ} ̂ρ_k(S^{-1} G)
        invS = Mat3{Int}(inv(symop.S))
        for (ig, G) in enumerate(G_vectors_generator(basis.fft_size))
            igired = index_G_vectors(basis, invS * G)
            if igired !== nothing
                @inbounds ρaccu[ig] += cis2pi(-T(dot(G, symop.τ))) * ρin[igired]
            end
        end
    end  # symop
    ρaccu
end

# Low-pass filters ρ (in Fourier) so that symmetry operations acting on it stay in the grid
function lowpass_for_symmetry!(ρ, basis; symmetries=basis.symmetries)
    for symop in symmetries
        symop == one(SymOp) && continue
        for (ig, G) in enumerate(G_vectors_generator(basis.fft_size))
            if index_G_vectors(basis, symop.S * G) === nothing
                ρ[ig] = 0
            end
        end
    end
    ρ
end

"""
Symmetrize a density by applying all the basis (by default) symmetries and forming the average.
"""
@views @timing function symmetrize_ρ(basis, ρ; symmetries=basis.symmetries)
    ρin_fourier = r_to_G(basis, ρ)
    ρout_fourier = zero(ρin_fourier)
    for σ = 1:size(ρ, 4)
        accumulate_over_symmetries!(ρout_fourier[:, :, :, σ],
                                    ρin_fourier[:, :, :, σ], basis, symmetries)
        lowpass_for_symmetry!(ρout_fourier[:, :, :, σ], basis; symmetries)
    end
    G_to_r(basis, ρout_fourier ./ length(symmetries))
end

"""
Symmetrize the stress tensor, given as a Matrix in cartesian coordinates
"""
function symmetrize_stresses(model::Model, stresses; symmetries)
    # see (A.28) of https://arxiv.org/pdf/0906.2569.pdf
    stresses_symmetrized = zero(stresses)
    for symop in symmetries
        W_cart = matrix_red_to_cart(model, symop.W)
        stresses_symmetrized += W_cart * stresses / W_cart
    end
    stresses_symmetrized /= length(symmetries)
    stresses_symmetrized
end
function symmetrize_stresses(basis::PlaneWaveBasis, stresses)
    symmetrize_stresses(basis.model, stresses; basis.symmetries)
end

"""
Symmetrize the forces in *reduced coordinates*, forces given as an
array forces[iel][α,i]
"""
function symmetrize_forces(model::Model, forces; symmetries)
    symmetrized_forces = zero(forces)
    for group in model.atom_groups, symop in symmetries
        positions_group = model.positions[group]
        W, w = symop.W, symop.w
        for (idx, position) in enumerate(positions_group)
            # see (A.27) of https://arxiv.org/pdf/0906.2569.pdf
            # (but careful that our symmetries are r -> Wr+w, not R(r+f))
            other_at = W \ (position - w)
            is_approx_integer(r) = all(ri -> abs(ri - round(ri)) ≤ SYMMETRY_TOLERANCE, r)
            i_other_at = findfirst(a -> is_approx_integer(a - other_at), positions_group)
            symmetrized_forces[idx] += W * forces[group[i_other_at]]
        end
    end
    symmetrized_forces / length(symmetries)
end
function symmetrize_forces(basis::PlaneWaveBasis, forces)
    symmetrize_forces(basis.model, forces; basis.symmetries)
end

""""
Convert a `basis` into one that doesn't use BZ symmetry.
This is mainly useful for debug purposes (e.g. in cases we don't want to
bother thinking about symmetries).
"""
function unfold_bz(basis::PlaneWaveBasis)
    if length(basis.symmetries) == 1
        return basis
    else
        kcoords = unfold_kcoords(basis.kcoords_global, basis.symmetries)
        new_basis = PlaneWaveBasis(basis.model,
                                   basis.Ecut, basis.fft_size, basis.variational,
                                   kcoords, [1/length(kcoords) for _ in kcoords],
                                   basis.kgrid, basis.kshift, basis.symmetries, basis.comm_kpts)
    end
end

# find where in the irreducible basis `basis_irred` the k-point `kpt_unfolded` is handled
function unfold_mapping(basis_irred, kpt_unfolded)
    for ik_irred = 1:length(basis_irred.kpoints)
        kpt_irred = basis_irred.kpoints[ik_irred]
        for symop in basis_irred.symmetries
            Sk_irred = normalize_kpoint_coordinate(symop.S * kpt_irred.coordinate)
            k_unfolded = normalize_kpoint_coordinate(kpt_unfolded.coordinate)
            if (Sk_irred ≈ k_unfolded) && (kpt_unfolded.spin == kpt_irred.spin)
                return ik_irred, symop
            end
        end
    end
    error("Invalid unfolding of BZ")
end

function unfold_array_(basis_irred, basis_unfolded, data, is_ψ)
    if basis_irred == basis_unfolded
        return data
    end
    if !(basis_irred.comm_kpts == basis_irred.comm_kpts == MPI.COMM_WORLD)
        error("Brillouin zone symmetry unfolding not supported with MPI yet")
    end
    data_unfolded = similar(data, length(basis_unfolded.kpoints))
    for ik_unfolded in 1:length(basis_unfolded.kpoints)
        kpt_unfolded = basis_unfolded.kpoints[ik_unfolded]
        ik_irred, symop = unfold_mapping(basis_irred, kpt_unfolded)
        if is_ψ
            # transform ψ_k from data into ψ_Sk in data_unfolded
            kunfold_coord = kpt_unfolded.coordinate
            @assert normalize_kpoint_coordinate(kunfold_coord) ≈ kunfold_coord
            _, ψSk = apply_symop(symop, basis_irred,
                                  basis_irred.kpoints[ik_irred], data[ik_irred])
            data_unfolded[ik_unfolded] = ψSk
        else
            # simple copy
            data_unfolded[ik_unfolded] = data[ik_irred]
        end
    end
    data_unfolded
end

function unfold_bz(scfres)
    basis_unfolded = unfold_bz(scfres.basis)
    ψ = unfold_array_(scfres.basis, basis_unfolded, scfres.ψ, true)
    eigenvalues = unfold_array_(scfres.basis, basis_unfolded, scfres.eigenvalues, false)
    occupation = unfold_array_(scfres.basis, basis_unfolded, scfres.occupation, false)
    E, ham = energy_hamiltonian(basis_unfolded, ψ, occupation;
                                scfres.ρ, eigenvalues, scfres.εF)
    @assert E.total ≈ scfres.energies.total
    new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation)
    merge(scfres, new_scfres)
end

function unfold_kcoords(kcoords, symmetries)
    all_kcoords = [normalize_kpoint_coordinate(symop.S * kcoord)
                   for kcoord in kcoords, symop in symmetries]

    # the above multiplications introduce an error
    unique(all_kcoords) do k
        digits = ceil(Int, -log10(SYMMETRY_TOLERANCE))
        normalize_kpoint_coordinate(round.(k; digits))
    end
end

ext: jl


text:
"Holds the tableau of a general linear method."
struct TableauGLM{T} <: AbstractTableau{T}
    name::Symbol
    o::Int
    s::Int
    r::Int

    a::Matrix{T}
    b::Matrix{T}
    u::Matrix{T}
    v::Matrix{T}
    c::Vector{T}

    function TableauGLM{T}(name, o, s, r, a, b, u, v, c) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)
        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"
        @assert s==length(c)
        @assert s==size(a,1)==size(a,2)==length(u,1)==length(b,2)
        @assert r==size(v,1)==size(v,2)==length(u,2)==length(b,1)
        new(name, o, s, r, a, b, u, v, c)
    end
end

# TODO Add external constructor for TableauGLM.

# TODO function readTableauGLMFromFile(dir::AbstractString, name::AbstractString)

# TODO function writeTableauToFile{Name, T}(dir::AbstractString, tab::TableauGLM{Name, T})


text:
using LinearAlgebra, StatsModels, ExperimentalDesign, Distributions, DataFrames

"""
$(TYPEDSIGNATURES)

```jldoctest
julia> 1 + 1
2
```
"""
function d_criterion(model_matrix;
                     tolerance = 1e-9)
    det(((model_matrix' * model_matrix) + (I * tolerance)) /
        size(model_matrix, 1)) ^
        (1 / size(model_matrix, 2))
end

"""
$(TYPEDSIGNATURES)

Builds optimum designs using the KL exchange algorithm, as described by Atkinson
*et al.*.   The ideia is  to iteratively swap  ``K`` design elements  with ``L``
candidate elements,  in order to  *maximize* an optimality  criterion.  Although
any criteria based  on information matrices could be used  here, the KL exchange
algorithm     leverages    determinant     properties     to    optimize     the
[`d_criterion`](@ref).

> Atkinson, A., Donev, A., & Tobias, R. (2007). Optimum experimental designs, with
> SAS (Vol. 34). Oxford University Press, Chapter 12.

```julia
julia> candidates = FullFactorial(fill([-1, 0, 1], 5));

julia> candidates_formula = ConstantTerm(0) ~ sum(Term.(Symbol.(names(candidates.matrix))));

julia> candidates_formula = FormulaTerm(candidates_formula.lhs,
                                        candidates_formula.rhs +
                                        (@formula 0 ~ factor3 ^ 2).rhs);

julia> selected_rows = kl_exchange(candidates_formula,
                                   candidates.matrix,
                                   seed_design_size = 2,
                                   experiments = 11,
                                   design_k = 11,
                                   candidates_l = size(candidates.matrix, 1) - 11);

julia> d_criterion(selected_rows)
0.730476820204009

julia> candidates.matrix[selected_rows.indices[1], :]
11×5 DataFrames.DataFrame
│ Row │ factor1 │ factor2 │ factor3 │ factor4 │ factor5 │
│     │ Int64   │ Int64   │ Int64   │ Int64   │ Int64   │
├─────┼─────────┼─────────┼─────────┼─────────┼─────────┤
│ 1   │ -1      │ 1       │ -1      │ 1       │ 1       │
│ 2   │ 1       │ -1      │ 1       │ 1       │ 1       │
│ 3   │ 1       │ 1       │ 0       │ -1      │ 1       │
│ 4   │ -1      │ 1       │ 1       │ -1      │ 1       │
│ 5   │ -1      │ -1      │ 0       │ 1       │ 1       │
│ 6   │ 1       │ 1       │ 1       │ 1       │ -1      │
│ 7   │ -1      │ -1      │ 1       │ -1      │ -1      │
│ 8   │ 1       │ -1      │ 0       │ 1       │ -1      │
│ 9   │ -1      │ 1       │ -1      │ 1       │ -1      │
│ 10  │ -1      │ 1       │ 0       │ -1      │ -1      │
│ 11  │ 1       │ -1      │ -1      │ -1      │ -1      │
```

```julia
julia> n_candidates = 3000;

julia> n_factors = 30;

julia> design_generator = DesignDistribution(Distributions.Uniform(0, 1), n_factors);

julia> candidates = rand(design_generator, n_candidates);

julia> selected_rows = kl_exchange(candidates.formula,
                                   candidates.matrix,
                                   experiments = n_factors + 2,
                                   design_k = n_factors,
                                   candidates_l = round(Int, (n_candidates - n_factors) / 2));

julia> d_criterion(selected_rows)
0.07904918814259465

julia> candidates.matrix[selected_rows.indices[1], :]
32×30 DataFrame. Omitted printing of 24 columns
│ Row │ factor1    │ factor2   │ factor3   │ factor4   │ factor5   │ factor6   │
│     │ Float64    │ Float64   │ Float64   │ Float64   │ Float64   │ Float64   │
├─────┼────────────┼───────────┼───────────┼───────────┼───────────┼───────────┤
│ 1   │ 0.832487   │ 0.875035  │ 0.433788  │ 0.805738  │ 0.481517  │ 0.382364  │
│ 2   │ 0.761498   │ 0.0398811 │ 0.998741  │ 0.632286  │ 0.29121   │ 0.55497   │
│ 3   │ 0.639793   │ 0.49432   │ 0.678795  │ 0.685274  │ 0.269203  │ 0.896397  │
│ 4   │ 0.00370418 │ 0.5139    │ 0.0630861 │ 0.175531  │ 0.570244  │ 0.60512   │
│ 5   │ 0.698514   │ 0.207593  │ 0.766977  │ 0.0719112 │ 0.665055  │ 0.499193  │
│ 6   │ 0.412744   │ 0.785263  │ 0.432861  │ 0.909976  │ 0.642545  │ 0.832013  │
│ 7   │ 0.353518   │ 0.116685  │ 0.140352  │ 0.929095  │ 0.0484366 │ 0.838524  │
⋮
│ 25  │ 0.577343   │ 0.790529  │ 0.763994  │ 0.548131  │ 0.402066  │ 0.293336  │
│ 26  │ 0.0941214  │ 0.52891   │ 0.293733  │ 0.92224   │ 0.982713  │ 0.929744  │
│ 27  │ 0.0256128  │ 0.723072  │ 0.23927   │ 0.993221  │ 0.521335  │ 0.407319  │
│ 28  │ 0.221631   │ 0.375189  │ 0.67251   │ 0.658613  │ 0.272448  │ 0.230501  │
│ 29  │ 0.843885   │ 0.375497  │ 0.116069  │ 0.989937  │ 0.935444  │ 0.0466424 │
│ 30  │ 0.181611   │ 0.956062  │ 0.98997   │ 0.80352   │ 0.994557  │ 0.203797  │
│ 31  │ 0.738407   │ 0.25753   │ 0.929704  │ 0.407498  │ 0.0934688 │ 0.99359   │
│ 32  │ 0.228975   │ 0.81272   │ 0.803107  │ 0.931826  │ 0.598371  │ 0.588823  │
```
"""
function kl_exchange(formula::FormulaTerm,
                     candidate_set::DataFrame;
                     tolerance::Float64 = 1e-9,
                     seed_design_size::Int = 2,
                     max_iterations::Int = 1000,
                     experiments::Int,
                     design_k::Int,
                     candidates_l::Int)
    model_matrix = ModelMatrix(ModelFrame(formula, candidate_set)).m

    seed = view(model_matrix,
                rand(1:size(model_matrix, 1), seed_design_size),
                :)

    l_variances = zeros(Int, candidates_l)
    k_variances = zeros(Int, design_k)

    design_information = zeros(size(model_matrix, 2),
                               size(model_matrix, 2))

    variances = zeros(size(model_matrix, 1) - size(seed, 1))

    # Seed Iterations
    for i = 1:(experiments - seed_design_size)
        seed_information = view(design_information, 1:size(seed, 2), 1:size(seed, 2))
        seed_information .= inv((seed' * seed) + (I * tolerance))

        candidates = view(model_matrix,
                          setdiff(1:size(model_matrix, 1), seed.indices[1]),
                          :)

        variances_candidates = view(variances, 1:size(candidates, 1))
        variances_candidates .= dot.(eachrow(candidates * seed_information),
                                     eachrow(candidates))

        l_variances .= partialsortperm(variances_candidates, 1:candidates_l, rev = true)

        seed = view(model_matrix,
                    vcat(candidates.indices[1][findmax(l_variances)[2]],
                         seed.indices[1]),
                    :)

    end

    design = seed
    variances_design = zeros(size(design, 1))
    variances_candidates = view(variances, 1:(size(model_matrix, 1) - size(design, 1)))

    kl_deltas = zeros(design_k, candidates_l)

    # Exchange Steps
    for i = 1:max_iterations
        design_information .= inv((design' * design) + (I * tolerance))
        candidates = view(model_matrix,
                          setdiff(1:size(model_matrix, 1), design.indices[1]),
                          :)

        variances_design .= dot.(eachrow(design * design_information),
                                 eachrow(design))
        k_variances .= partialsortperm(variances_design, 1:design_k, rev = false)

        variances_candidates .= dot.(eachrow(candidates * design_information),
                                     eachrow(candidates))
        l_variances .= partialsortperm(variances_candidates, 1:candidates_l, rev = true)

        k_swappable = view(model_matrix,
                           [design.indices[1][k] for k in k_variances],
                           :)

        l_swappable = view(model_matrix,
                           [candidates.indices[1][l] for l in l_variances],
                           :)

        kl_deltas .= k_swappable * design_information * l_swappable'

        swap = ((-1, -1), -1)

        for k = 1:design_k
            for l = 1:candidates_l
                delta = ((1 - variances_design[k_variances[k]]) *
                         (1 + variances_candidates[l_variances[l]])) +
                         (kl_deltas[k, l] ^ 2)

                if delta > swap[2]
                    swap = ((k, l), delta)
                end
            end
        end

        if swap[2] - 1.0 <= tolerance
            break
        end

        design = view(model_matrix,
                      vcat(deleteat!(copy(design.indices[1]),
                                     k_variances[swap[1][1]]),
                           candidates.indices[1][l_variances[swap[1][2]]]),
                      :)
    end

    return design
end

function test()
    c = 300000
    n = 50

    println("Allocating memory")
    design_generator = DesignDistribution(Distributions.Uniform(0, 1), n)
    @time candidates = rand(design_generator, c)

    println("Design created, calling exchange")
    println(design_generator.formula)

    @time b = kl_exchange(design_generator.formula,
                          candidates,
                          experiments = n + 2,
                          design_k = n - 3,
                          candidates_l = round(Int, (c - n) / 2))
    return b
end


text:
using GPUArrays, Colors, FileIO, ImageFiltering
using CLArrays

img = RGB{Float32}.(load(joinpath(@__DIR__, "..", "Poincare", "poincare.png")))

a = CLArray(img);
out = similar(a);
k = CLArray(Float32.(collect(Kernel.gaussian(7))));
imgc = similar(img)

GPUArrays.convolution!(a, out, k)
copy!(imgc, out)
save(joinpath(@__DIR__, "gaussian.png"), imgc)


text:
# approximation algorithm from ``Algorithm 917: Complex Double-Precision
# Evaluation of the Wright ω Function'' by Piers W. Laurence, Robert M. Corless
# and David J. Jeffrey
# Series approximations can be found in the book chapter ``The Wright ω Function''
# by Corless, R. M. and Jeffrey, D. J.
#
# see also https://github.com/emsr/tr29124_test

module WrightOmega

export wrightomega

# TODO handle infs
# TODO don't use two methods
# FIXME test failures around zero and regularization

function wrightomega(z::T) where {T <: Real}
    if z < -2
        t = exp(z)
        w = t * (1 + t * (-1 + t * (T(3) / 2 + t * (-T(8) / 3 + T(125) / 24 * t))))
    elseif z < 1 + T(π)
        z1 = z - 1
        w = 1 + z1 / 2 * (1 + z1 / 8 * (1 + z1 / 12 * (-1 + z1 / 16 *
            (-1 + z1 * T(13) / 20))))
    else
        l = log(z)
        w = z + l * (-1 + (1 + (l / 2 - 1 + (l * (l / 3 - T(3) / 2) + 1) / z) / z) / z)
    end
    # refinement steps
    for _ in 1:5
        r = z - w - log(w)
        w1 = w + 1
        t = w1 * (w1 + T(2) / 3 * r)
        w *= 1 + r / w1 * (t - r / 2) / (t - r)
        fscn = abs(r^4 * (2 * w * (w - 4) - 1))
        if fscn < eps(float(T)) * 72 * w1^6
            break
        end
    end
    return w::float(T)
end

function wrightomega(z::Complex{T}) where {T <: Real}
    x = real(z)
    y = imag(z)
    Tπ = float(T)(π)
    if x == -1 && abs(y) == Tπ
        return zero(Complex{T})
    end
    # region 1
    if -2 < x <= 1 && 1 < y < 2Tπ
        imp = im * conj(sqrt(2 * conj(z + 1 - im * Tπ)))
        w = -1 + imp * (1 + imp / 3 * (-1 + imp / 6 * (inv(T(2)) + imp / 15 *
            (1 + imp / 16))))
    # region 2
    elseif -2 < x <= 1 && -2Tπ < y < -1
        imp = im * conj(sqrt(2 * conj(z + 1 - im * Tπ)))
        w = -1 + imp * (-1 + imp / 3 * (-1 + imp / 6 * (-inv(T(2)) + imp / 15 *
            (1 - imp / 16))))
    # region 3
    elseif x < -2 && -Tπ < y <= Tπ
        t = exp(z)
        w = t * (1 + t * (-1 + t * (T(3) / 2 + t * (-T(8) / 3 + T(125) / 24 * t))))
    # region 4
    elseif (-2 < x <= 1 && -1 <= y <= 1) || (-1 < x && (x - 1)^2 + y^2 <= Tπ^2)
        z1 = z - 1
        w = 1 + z1 / 2 * (1 + z1 / 8 * (1 + z1 / 12 * (-1 + z1 / 16 *
            (-1 + z1 * T(13) / 20))))
    # region 5
    elseif x <= -2 && Tπ < y && y - Tπ <= -3 / T(4) * (x + 1)
        t = z - im * Tπ
        l = log(-t)
        w = t - l + l * (1 + (l / 2 - 1 + (l^2 / 3 - 3 * l / 2 + 1) / t) / t) / t
    # region 6
    elseif x <= -2 && y <= -Tπ && y + Tπ >= T(3) / 4 * (x + 1)
        t = z + im * Tπ
        l = log(-t)
        w = t - l + l * (1 + (l / 2 - 1 + (l^2 / 3 - 3 * l / 2 + 1) / t) / t) / t
    # region 7
    else
        l = log(z)
        w = z + l * (-1 + (1 + (l / 2 - 1 + (l * (l / 3 - T(3) / 2) + 1) / z) / z) / z)
    end
    if x <= -0.99 && (abs(y - Tπ) <= 0.01 || abs(y + Tπ) <= 0.01)
        s = -1
        if abs(y - Tπ) <= 0.01
            # TODO setrounding doesn't work for non-BigFloat, decide what to do
            # T == BigFloat && setrounding(T, RoundUp)
            z = x + im * (y - Tπ)
        else
            z = x + im * (y + Tπ)
        end
    else
        s = 1
    end
    w *= s
    # refinement steps
    for _ in 1:5
        r = z - s * w - log(Complex(w))
        w1 = s * w + 1
        t = w1 * (w1 + T(2) / 3 * r)
        w *= 1 + r / w1 * (t - r / 2) / (t - r)
        fscn = abs(r^4 * (2 * w * (w - 4) - 1))
        if fscn < eps(float(T)) * 72 * abs(w1)^6
            break
        end
    end
    w *= s
    return w::Complex{float(T)}
end

end


text:
"""
    mutable struct Event{T <: SSDFloat}

Collection struct for individual events.
This (mutable) struct is meant to be used to look at individual events,
not to process a huge amount of events.

## Fields
* `locations::VectorOfArrays{CartesianPoint{T}}`: Vector of the positions of all hits of the event.
* `energies::VectorOfArrays{T}`: Vector of energies corresponding to the hits of the event.
* `drift_paths::Union{Vector{EHDriftPath{T}}, Missing}`: Calculated drift paths of each hit position.
* `waveforms::Union{Vector{<:Any}, Missing}`: Generated signals (waveforms) of the event.
"""
mutable struct Event{T <: SSDFloat}
    locations::VectorOfArrays{CartesianPoint{T}}
    energies::VectorOfArrays{T}
    drift_paths::Union{Vector{EHDriftPath{T}}, Missing}
    waveforms::Union{Vector{<:Any}, Missing}
    Event{T}() where {T <: SSDFloat} = new{T}(VectorOfArrays(Vector{CartesianPoint{T}}[]), VectorOfArrays(Vector{T}[]), missing, missing)
end

function Event(location::AbstractCoordinatePoint{T}, energy::RealQuantity = one(T))::Event{T} where {T <: SSDFloat}
    evt = Event{T}()
    evt.locations = VectorOfArrays([[CartesianPoint(location)]])
    evt.energies = VectorOfArrays([[T(to_internal_units(energy))]])
    return evt
end

function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity} = ones(T, length(locations)))::Event{T} where {T <: SSDFloat}
    evt = Event{T}()
    evt.locations = VectorOfArrays(broadcast(pt -> [CartesianPoint(pt)], locations))
    evt.energies = VectorOfArrays(broadcast(E -> [T(to_internal_units(E))], energies))
    return evt
end

function Event(locations::Vector{<:Vector{<:AbstractCoordinatePoint{T}}}, energies::Vector{<:Vector{<:RealQuantity}} = [[one(T) for i in j] for j in locations])::Event{T} where {T <: SSDFloat}
    evt = Event{T}()
    evt.locations = VectorOfArrays(broadcast(pts -> CartesianPoint{T}.(pts), locations))
    evt.energies = VectorOfArrays(Vector{T}.(to_internal_units.(energies)))
    return evt
end

function Event(nbcc::NBodyChargeCloud{T})::Event{T} where {T <: SSDFloat}
    evt = Event{T}()
    evt.locations = VectorOfArrays([nbcc.locations])
    evt.energies = VectorOfArrays([nbcc.energies])
    return evt
end

function Event(nbccs::Vector{<:NBodyChargeCloud{T}})::Event{T} where {T <: SSDFloat}
    evt = Event{T}()
    evt.locations = VectorOfArrays(broadcast(nbcc -> nbcc.locations, nbccs))
    evt.energies = VectorOfArrays(broadcast(nbcc -> nbcc.energies, nbccs))
    return evt
end

function Event(locations::Vector{<:AbstractCoordinatePoint{T}}, energies::Vector{<:RealQuantity}, N::Int;
               particle_type::Type{PT} = Gamma, number_of_shells::Int = 2,
               radius::Vector{<:RealQuantity{T}} = radius_guess.(T.(to_internal_units.(energies)), particle_type)
              )::Event{T} where {T <: SSDFloat, PT <: ParticleType}

    return Event(broadcast(i ->
                NBodyChargeCloud(locations[i], energies[i], N, particle_type,
                radius = radius[i], number_of_shells = number_of_shells),
           eachindex(locations)))
end

function Event(evt::NamedTuple{(:evtno, :detno, :thit, :edep, :pos),
         <:Tuple{
            Union{Integer, AbstractVector{<:Integer}},
            Union{Integer, AbstractVector{<:Integer}},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:RealQuantity},
            AbstractVector{<:AbstractVector{<:RealQuantity}}
        }}, T = missing)
    if ismissing(T) T = eltype(to_internal_units(evt[:edep][:])) end

    evt = Event(
        [CartesianPoint{T}.(to_internal_units.(evt[:pos][:]))],
        [T.(to_internal_units.(evt[:edep][:]))]
    )

    return evt
end

in(evt::Event, det::SolidStateDetector) = all( pt -> pt in det, flatview(evt.locations))
in(evt::Event, sim::Simulation) = all( pt -> pt in sim.detector, flatview(evt.locations))

function move_charges_inside_semiconductor!(evt::Event{T}, det::SolidStateDetector{T}; fraction::T = T(0.2))::Event{T} where {T <: SSDFloat}
    move_charges_inside_semiconductor!(evt.locations, evt.energies, det; fraction)
    evt
end

function move_charges_inside_semiconductor!(
        locations::AbstractVector{<:AbstractVector{CartesianPoint{T}}}, energies::AbstractVector{<:AbstractVector{T}},
        det::SolidStateDetector{T}; fraction::T = T(0.2)) where {T <: SSDFloat}
    for n in eachindex(locations)
        idx_in = broadcast( pt -> pt in det.semiconductor, locations[n]);
        if !all(idx_in)
            charge_center = sum(locations[n] .* energies[n]) / sum(energies[n])
            @assert charge_center in det.semiconductor "The center of the charge cloud ($(charge_center)) is not inside the semiconductor."
            surf = ConstructiveSolidGeometry.surfaces(det.semiconductor.geometry)
            for (k,m) in enumerate(findall(.!idx_in))
                l = locations[n][m]
                line = ConstructiveSolidGeometry.Line(charge_center, CartesianVector{T}(l - charge_center))
                for s in surf
                    int = ConstructiveSolidGeometry.intersection(s, line)
                    for i in int
                        proj::T = (i - charge_center) ⋅ (l - charge_center) / sum((l - charge_center).^2)
                        if 0 ≤ proj ≤ 1
                            if i in det.semiconductor
                                locations[n][m] = (1 - fraction * proj) * i + fraction * proj * charge_center
                            end
                        end
                    end
                end
            end
            charge_center_new = sum(locations[n] .* energies[n]) / sum(energies[n])
            @warn "$(sum(.!idx_in)) charges of the charge cloud at $(round.(charge_center, digits = (T == Float64 ? 12 : 6)))"*
            " are outside. Moving them inside...\nThe new charge center is at $(round.(charge_center_new, digits = (T == Float64 ? 12 : 6))).\n"
        end
    end
    nothing
end

"""
    drift_charges!(evt::Event{T}, sim::Simulation{T}; kwargs...)::Nothing where {T <: SSDFloat}

Calculates the electron and hole drift paths for the given [`Event`](@ref) and [`Simulation`](@ref)
    and stores them in `evt.drift_paths`.

## Arguments
* `evt::Event{T}`: [`Event`](@ref) for which the charges should be drifted.
* `sim::Simulation{T}`: [`Simulation`](@ref) which defines the setup in which the charges in `evt` should drift.

## Keywords
* `max_nsteps::Int = 1000`: Maximum number of steps in the drift of each hit.
* `Δt::RealQuantity = 5u"ns"`: Time step used for the drift.
* `diffusion::Bool = false`: Activate or deactive diffusion of charge carriers via random walk.
* `self_repulsion::Bool = false`: Activate or deactive self-repulsion of charge carriers of the same type.
* `verbose = true`: Activate or deactivate additional info output.

## Example
```julia
drift_charges!(evt, sim, Δt = 1u"ns", verbose = false)
```

!!! note
    Using values with units for `Δt` requires the package [Unitful.jl](https://github.com/PainterQubits/Unitful.jl).
"""
function drift_charges!(evt::Event{T}, sim::Simulation{T}; max_nsteps::Int = 1000, Δt::RealQuantity = 5u"ns", diffusion::Bool = false, self_repulsion::Bool = false, verbose::Bool = true)::Nothing where {T <: SSDFloat}
    !in(evt, sim.detector) && move_charges_inside_semiconductor!(evt, sim.detector)
    evt.drift_paths = drift_charges(sim, evt.locations, evt.energies, Δt = Δt, max_nsteps = max_nsteps, diffusion = diffusion, self_repulsion = self_repulsion, verbose = verbose)
    nothing
end
function get_signal!(evt::Event{T}, sim::Simulation{T}, contact_id::Int; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}
    @assert !ismissing(evt.drift_paths) "No drift path for this event. Use `drift_charges(evt::Event, sim::Simulation)` first."
    @assert !ismissing(sim.weighting_potentials[contact_id]) "No weighting potential for contact $(contact_id). Use `calculate_weighting_potential!(sim::Simulation, contact_id::Int)` first."
    if ismissing(evt.waveforms)
        evt.waveforms = Union{Missing, RadiationDetectorSignals.RDWaveform}[missing for i in eachindex(sim.detector.contacts)]
    end
    evt.waveforms[contact_id] = get_signal(sim, evt.drift_paths, flatview(evt.energies), contact_id, Δt = Δt)
    nothing
end

"""
    get_signals!(evt::Event{T}, sim::Simulation{T}; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}

Generates the signals/waveforms from the drift paths of an [`Event`](@ref) for each [`Contact`](@ref),
for which a [`WeightingPotential`](@ref) is specified in `sim.weighting_potentials`.

The output is stored in `evt.waveforms`.

## Arguments
* `evt::Event{T}`: [`Event`](@ref) for which the waveforms should be generated.
* `sim::Simulation{T}`: [`Simulation`](@ref) which defines the setup in which the waveforms are generated.

## Keywords
* `Δt::RealQuantity = 5u"ns"`: Time steps with which the drift paths were calculated.

## Example
```julia
get_signals!(evt, sim, Δt = 1u"ns") # if evt.drift_paths were calculated in time steps of 1ns
```

!!! note
    This method only works if `evt.drift_paths` has already been calculated and is not `missing`.
"""
function get_signals!(evt::Event{T}, sim::Simulation{T}; Δt::RealQuantity = 5u"ns")::Nothing where {T <: SSDFloat}
    @assert !ismissing(evt.drift_paths) "No drift path for this event. Use `drift_charges!(evt::Event, sim::Simulation)` first."
    if ismissing(evt.waveforms)
        evt.waveforms = Union{Missing, RadiationDetectorSignals.RDWaveform}[missing for i in eachindex(sim.detector.contacts)]
    end
    for contact in sim.detector.contacts
        if any(ismissing, sim.weighting_potentials) "No weighting potential(s) for some contact(s).." end
        if !ismissing(sim.weighting_potentials[contact.id])
            evt.waveforms[contact.id] = get_signal(sim, evt.drift_paths, flatview(evt.energies), contact.id, Δt = Δt)
        end
    end
    nothing
end

"""
    simulate!(evt::Event{T}, sim::Simulation{T}; kwargs...)::Nothing where {T <: SSDFloat}

Simulates the waveforms for the [`Event`](@ref) for a given [`Simulation`](@ref) by

1. calculating the drift paths of all energy hits, at `evt.locations` and
2. generating the waveforms for each [`Contact`](@ref), for which a [`WeightingPotential`](@ref) is specified in `sim.weighting_potentials`.

The output is stored in `evt.drift_paths` and `evt.waveforms`.

## Arguments
* `evt::Event{T}`: [`Event`](@ref) for which the charges should be drifted.
* `sim::Simulation{T}`: [`Simulation`](@ref) which defines the setup in which the charges in `evt` should drift.

## Keywords
* `max_nsteps::Int = 1000`: Maximum number of steps in the drift of each hit.
* `Δt::RealQuantity = 5u"ns"`: Time step used for the drift.
* `diffusion::Bool = false`: Activate or deactive diffusion of charge carriers via random walk.
* `self_repulsion::Bool = false`: Activate or deactive self-repulsion of charge carriers of the same type.
* `verbose = true`: Activate or deactivate additional info output.

## Example
```julia
simulate!(evt, sim, Δt = 1u"ns", verbose = false)
```

See also [`drift_charges!`](@ref) and [`get_signals!`](@ref).
"""
function simulate!(evt::Event{T}, sim::Simulation{T}; max_nsteps::Int = 1000, Δt::RealQuantity = 5u"ns", diffusion::Bool = false, self_repulsion::Bool = false, verbose::Bool = true)::Nothing where {T <: SSDFloat}
    drift_charges!(evt, sim, max_nsteps = max_nsteps, Δt = Δt, diffusion = diffusion, self_repulsion = self_repulsion, verbose = verbose)
    get_signals!(evt, sim, Δt = Δt)
    nothing
end

"""
    add_baseline_and_extend_tail(wv::RadiationDetectorSignals.RDWaveform{T,U,TV,UV}, n_baseline_samples::Int, total_waveform_length::Int) where {T,U,TV,UV}

Adds a zero-valued baseline in front of the waveform `wv` and extends (or cuts off) the waveform at the end with the last value of `wv`.
A waveform of length `total_waveform_length` is returned.

## Arguments
* `wv::RadiationDetectorSignals.RDWaveform{T,U,TV,UV}`: A waveform (signal over time).
* `n_baseline_samples::Int`: Number of samples added in front of the waveform with values 0.
* `total_waveform_length::Int`: Number of samples of the extended waveform which is returned.

## Examples
    add_baseline_and_extend_tail(wv, 1000, 5000)

!!! note
    This functions assumes that the time steps between the samples of the input waveform `wv` are the same. Thus, that the input waveform is
    sampled with a fixed frequency.
"""
function add_baseline_and_extend_tail(wv::RadiationDetectorSignals.RDWaveform{T,U,TV,UV}, n_baseline_samples::Int, total_waveform_length::Int) where {T,U,TV,UV}
    new_signal::Vector{eltype(UV)} = Vector{eltype(UV)}(undef, total_waveform_length)
    new_signal[1:n_baseline_samples] .= zero(eltype(UV))
    if length(wv.value) <= total_waveform_length - n_baseline_samples
        new_signal[n_baseline_samples+1:n_baseline_samples+length(wv.value)] = wv.value
        new_signal[n_baseline_samples+length(wv.value)+1:end] .= wv.value[end]
    else
        new_signal[n_baseline_samples+1:end] = wv.value[1:total_waveform_length - n_baseline_samples]
    end
    new_times = if TV <: AbstractRange
        range( zero(first(wv.time)), step = step(wv.time), length = total_waveform_length )
    else
        error("Not yet definted for timestamps of type `$(TV)`")
    end
    return RDWaveform( new_times, new_signal )
end

"""
    get_electron_and_hole_contribution(evt::Event{T}, sim::Simulation{T}, contact_id::Int)

Returns the electron and hole contribution to the waveform of a [`Contact`](@ref) with a given
`contact_id` of an [`Event`](@ref) as a `NamedTuple` with two entries:
`electron_contribution` and `hole_contribution`.

## Arguments
* `evt::Event{T}`: [`Event`](@ref) in which the charges have already been drifted.
* `sim::Simulation{T}`: [`Simulation`](@ref) which defines the setup in which the charges in `evt` were drifted.
* `contact_id::Int`: The `id` of the [`Contact`](@ref) for which the waveform should be split into electron and hole contribution.

## Example
```julia
using Plots
using SolidStateDetector
T = Float32

simulation = Simulation{T}(SSD_examples[:InvertedCoax])
simulate!(simulation)
event = Event([CartesianPoint{T}(0.02,0.01,0.05)])
simulate!(event, simulation)

contact_id = 1
wf = get_electron_and_hole_contribution(evt, sim, contact_id)
```

!!! note
    The drift paths in `evt` need to be calculated using [`drift_charges!`](@ref) before calling this function.

See also [`plot_electron_and_hole_contribution`](@ref).
"""
function get_electron_and_hole_contribution(evt::Event{T}, sim::Simulation{T, S}, contact_id::Int
            )::NamedTuple{(:electron_contribution, :hole_contribution), <:Tuple{RDWaveform, RDWaveform}} where {T <: SSDFloat, S}

    @assert !ismissing(evt.drift_paths) "The charge drift is not yet simulated. Please use `drift_charges!(evt, sim)`!"

    dt::T = T(ustrip(uconvert(u"ns", diff(evt.drift_paths[1].timestamps_e)[1]*u"s")))
    wp::Interpolations.Extrapolation{T, 3} = interpolated_scalarfield(sim.weighting_potentials[contact_id])
    signal_e::Vector{T} = zeros(T, length(maximum(map(p -> p.timestamps_e, evt.drift_paths))))
    signal_h::Vector{T} = zeros(T, length(maximum(map(p -> p.timestamps_h, evt.drift_paths))))

    for i in eachindex(evt.drift_paths)
        energy = flatview(evt.energies)[i]

        dp_e::Vector{CartesianPoint{T}} = evt.drift_paths[i].e_path
        dp_e_t::Vector{T} = evt.drift_paths[i].timestamps_e
        add_signal!(signal_e, dp_e_t, dp_e, dp_e_t, T(-1)*energy, wp, S)

        dp_h::Vector{CartesianPoint{T}} = evt.drift_paths[i].h_path
        dp_h_t::Vector{T} = evt.drift_paths[i].timestamps_h
        add_signal!(signal_h, dp_h_t, dp_h, dp_h_t, energy, wp, S)
    end
    unitless_energy_to_charge = _convert_internal_energy_to_external_charge(sim.detector.semiconductor.material)
    return (electron_contribution = RDWaveform(range(zero(T) * u"ns", step = dt * u"ns", length = length(signal_e)), signal_e * unitless_energy_to_charge),
            hole_contribution = RDWaveform(range(zero(T) * u"ns", step = dt * u"ns", length = length(signal_h)), signal_h * unitless_energy_to_charge))
end

export get_electron_and_hole_contribution

"""
    plot_electron_and_hole_contribution(evt::Event{T}, sim::Simulation{T}, contact_id::Int)

Plots the waveform as well as the electron and hole contribution to the waveform
of a [`Contact`](@ref) with a given `contact_id` of an [`Event`](@ref).

## Arguments
* `evt::Event{T}`: [`Event`](@ref) in which the charges have already been drifted.
* `sim::Simulation{T}`: [`Simulation`](@ref) which defines the setup in which the charges in `evt` were drifted.
* `contact_id::Int`: The `id` of the [`Contact`](@ref) for which the waveform should be split into electron and hole contribution.

## Keywords
* `n_samples::Int`: Number of samples with which the waveforms will be plotted. The default is the number of samples of the original waveform.

## Example
```julia
using Plots
using SolidStateDetector
T = Float32

simulation = Simulation{T}(SSD_examples[:InvertedCoax])
simulate!(simulation)
event = Event([CartesianPoint{T}(0.02,0.01,0.05)])
simulate!(event, simulation)

contact_id = 1
plot_electron_and_hole_contribution(evt, sim, contact_id, n_samples = 300)
```

!!! note
    The drift paths in `evt` need to be calculated using [`drift_charges!`](@ref) before calling this function.

!!! note
    This method requires to load the package [Plots.jl](https://github.com/JuliaPlots/Plots.jl).

See also [`get_electron_and_hole_contribution`](@ref).
"""
function plot_electron_and_hole_contribution end

@userplot Plot_electron_and_hole_contribution
@recipe function f(gdd::Plot_electron_and_hole_contribution; linewidth = 2, n_samples = missing)

    sim::Simulation = gdd.args[2]
    T = get_precision_type(sim)
    evt::Event{T} = gdd.args[1]
    contact_id::Int = gdd.args[3]

    ismissing(n_samples) ? n_samples = length(evt.waveforms[contact_id].value) : nothing

    wf::NamedTuple{(:electron_contribution, :hole_contribution), <:Tuple{RDWaveform, RDWaveform}} = get_electron_and_hole_contribution(evt, sim, contact_id)

    unitformat --> :slash
    yguide --> "Charge"

    @series begin
        linecolor := :red
        linewidth --> linewidth
        label --> "Electron contribution"
        add_baseline_and_extend_tail(wf.electron_contribution, 0, n_samples)
    end

    @series begin
        linecolor := :green
        linewidth --> linewidth
        label --> "Hole contribution"
        add_baseline_and_extend_tail(wf.hole_contribution, 0, n_samples)
    end

    @series begin
        linecolor --> contact_id
        linewidth --> linewidth
        seriesalpha --> 0.5
        label --> "Sum"
        add_baseline_and_extend_tail(evt.waveforms[contact_id], 0, n_samples)
    end
end


text:
#**************************************************************************************
# make.jl
# =============== part of the GeoEfficiency.jl package.
#
# script for building documentation of the GeoEfficiency.jl package.
#
#  _args = ["clean", "doctest"]
#  include(raw"C:\Users\Mohamed\.julia\dev\GeoEfficiency\docs\mkdocs.jl")
#**************************************************************************************

using Documenter, DocumenterMarkdown
using GeoEfficiency

_args = @isdefined(_args) ? _args : ARGS
const PAGES = Any[
    "Home" => "index.md",
    "Manual" => [
		"manual/GeoEfficiency.md",
        #"manual/Error.md",
        #"manual/Input_Console.md",
        "manual/Physics_Model.md",
        "manual/Calculations.md",
        "manual/Output_Interface.md",
        "manual/Input_Batch.md",
     ],

    "Development" => "manual/Development.md",
    "Index" => "list.md",
    "LICENSE" => "LICENSE.md",
]

makedocs(
    format  = Markdown(),
    modules = [GeoEfficiency],
    clean   = "clean" in _args,
    doctest = "doctest" in _args,
    sitename= "GeoEfficiency.jl",
    authors = "Mohamed E. Krar",
    pages   = PAGES,
    #assets  = ["assets/custom.css"],
)

const REPO = "github.com/DrKrar/GeoEfficiency.jl.git" # "github.com/GeoEfficiency/GeoEfficiency.github.io.git"
const VERSIONS =["stable" => "v^", "v#.#", "dev" => "dev"]  # order of versions in drop down menu.
const BRANCH = "gh-pages" # "master"

# Deploy Markup pages
@info "Deploying MarkUp pages"
deploydocs(
    repo = REPO,
    branch = "gh-pages-Mark", #BRANCH, "gh-pages",
    target = "build",
	versions = VERSIONS,
)

@info "Deploying MkDocs pages"
deploydocs(
    repo = REPO,
    deps   = Deps.pip("mkdocs", "pygments", "python-markdown-math"),
    make   = () -> run(`mkdocs build`),
    branch = "gh-pages-MkDocs", #BRANCH, "gh-pages",
    target = "site"
)


text:
module BioRecordsProcessing

    using Glob, CodecZlib, BGZFStreams

    function process_directory(process_function, ReadType, input_directory, pattern, output_directory;
            prefix = "",
            max_records = Inf,
            paired = false,
        )

        files = glob(pattern, input_directory)
        println.(files)
        mkpath(output_directory)

        Threads.@threads for file in files
            process(process_function, ReadType, file, output_directory;
                prefix = prefix, max_records = max_records
            )
        end
    end

    function get_writer_reader(ReadType, file, extension, out_file)
        if extension == ".gz"
            reader = ReadType.Reader(GzipDecompressorStream(open(file)))
            writer = ReadType.Writer(GzipCompressorStream(open(out_file, "w")))

        elseif extension == ".bam"
            index_file = file * ".bai"
            if !isfile(index_file)
                @warn "Index file not found : $index_file"
                reader = ReadType.Reader(open(file))
            else
                reader = ReadType.Reader(open(file); index = index_file)
            end

            writer = ReadType.Writer(BGZFStream(out_file, "w"), reader.header)
        else
            reader = ReadType.Reader(open(file))
            writer = ReadType.Writer(open(out_file, "w"))
        end
        reader, writer
    end

    function process(process_function, ReadType, file, output_directory;
        prefix = "",
        max_records = Inf,
    )

        filename, extension = splitext(basename(file))
        out_file = joinpath(output_directory, string(filename, prefix, extension))
        @assert out_file != file

        reader, writer = get_writer_reader(ReadType, file, extension, out_file)

        record = ReadType.Record()
        k = 0
        while !eof(reader)

            read!(reader, record)
            out_record = process_function(record)
            !isnothing(out_record) && write(writer, out_record)

            k += 1
            k > max_records && break
            if mod(k, 100_000) == 0
                @info "$(Threads.threadid()), $(basename(file)) : Processed $(div(k, 1000))k records..."
            end
        end
        close(reader)
        close(writer)
    end

    # paired files

    function process_directory_paired(process_function, ReadType, input_directory, pattern, to_file2, output_directory;
        prefix = "",
        max_records = Inf,
        paired = false,
    )

        files1 = glob(pattern, input_directory)
        files2 = to_file2.(files1)
        println.(zip(files1, files2))
        mkpath(output_directory)

        Threads.@threads for i in eachindex(files1)
            process_paired(process_function, ReadType, files1[i], files2[i], output_directory;
                prefix = prefix, max_records = max_records
            )
        end
    end

    function process_paired(process_function, ReadType, file1, file2, output_directory;
        prefix = "",
        max_records = Inf,
    )

        filename1, extension1 = splitext(basename(file1))
        filename2, extension2 = splitext(basename(file2))
        out_file1 = joinpath(output_directory, string(filename1, prefix, extension1))
        out_file2 = joinpath(output_directory, string(filename2, prefix, extension2))
        @assert out_file1 != file1
        @assert out_file2 != file2

        reader1, writer1 = get_writer_reader(ReadType, file1, extension1, out_file1)
        reader2, writer2 = get_writer_reader(ReadType, file2, extension2, out_file2)

        record1 = ReadType.Record()
        record2 = ReadType.Record()
        k = 0
        while !(eof(reader1) || eof(reader2))

            read!(reader1, record1)
            read!(reader2, record2)
            out_record1, out_record2 = process_function(record1, record2)
            !isnothing(out_record1) && write(writer1, out_record1)
            !isnothing(out_record2) && write(writer2, out_record2)

            k += 2
            k > max_records && break
            if mod(k, 100_000) == 0
                @info "$(Threads.threadid()), $(basename(file1)) : Processed $(div(k, 1000))k records..."
            end
        end
        close(reader1)
        close(reader2)
        close(writer1)
        close(writer2)
    end

    # Grouped things

    struct RecordGrouper{R, G, F <: Function, FD <: Function}
        records::Vector{R}
        free_idx::Set{Int}
        groups::Dict{G, Vector{Int}}
        get_key::F
        isdone::FD

        RecordGrouper{R,G}(get_key::Function, isdone::Function) where {R, G} =
            new{R, G, typeof(get_key), typeof(isdone)}(
                [R() for i=1:1000],
                Set{Int}(1:1000),
                Dict{G, Vector{Int}}(),
                get_key,
                isdone
            )
    end

    function check_free_idx!(rg::RecordGrouper{R,G}) where {R,G}
        if isempty(rg.free_idx)
            N = length(rg.records)
            for i = N:2N
                push!(rg.records, R())
                push!(rg.free_idx, i)
            end
        end
    end

    function get_record(rg::RecordGrouper{R,G}) where {R,G}
        check_free_idx!(rg)
        idx = pop!(rg.free_idx)
        rg.records[idx], idx
    end

    function group_record!(rg::RecordGrouper, r, idx)
        key = rg.get_key(r)
        if haskey(rg.groups, key)
            push!(rg.groups[key], idx)
        else
            rg.groups[key] = [idx]
        end
        records = (rg.records[i] for i in rg.groups[key])
        key, rg.isdone(records)
    end

    function process(process_function, ReadType, rg::RecordGrouper, file, output_directory;
        prefix = "",
        max_records = Inf,
    )

        filename, extension = splitext(basename(file))
        out_file = joinpath(output_directory, string(filename,'.', prefix, extension))
        @assert out_file != file

        reader, writer = get_writer_reader(ReadType, file, extension, out_file)

        k = 0
        while !eof(reader)

            record, idx = get_record(rg)
            read!(reader, record)

            key, isdone = group_record!(rg, record, idx)

            if isdone
                records = (fg.records[i] for i in fg.groups[key])
                out_records = process_function(records...)
                for r in out_records
                    write(writer, r)
                end
            end

            k += 1
            k > max_records && break
            if mod(k, 100_000) == 0
                @info "$(Threads.threadid()), $(basename(f)) : Processed $(div(k, 1000))k records..."
            end
        end
        close(reader)
        close(writer)
    end

end

max_stars_repo_name: jonathanBieler/BioRecordsProcessing

max_issues_repo_name: jonathanBieler/BioRecordsProcessing

max_forks_repo_name: jonathanBieler/BioRecordsProcessing


text:
### A Pluto.jl notebook ###
# v0.14.0

using Markdown
using InteractiveUtils

# ╔═╡ b926bbf4-ab68-11eb-136c-952f90c757df
using DifferentialEquations, PyPlot

# ╔═╡ ab5e06ed-f1f3-41a3-99e2-2affda673af8
PyPlot.matplotlib.style.use("Solarize_Light2")

# ╔═╡ 1b8c79c1-d663-42dc-beda-29a881dcbd2a
begin
	# define matrices
	A = [0    1;
		 0  -1/5]
	B = [  0      0;
		 -3/10  -3/5]
	C = [0, 1/10]
	E = [3 6]
end

# ╔═╡ 3f443782-32d8-44f0-baf3-41645d949573
y_sp(t) = t > 0.0 ? 1.0 : 0.0

# ╔═╡ b9aa78a5-d7bb-4bd8-a38d-46ef023670ca
# right-hand side of ODE in x
f(x, h, params, t) = A * x + B * h(params, t - 0.3) + C * y_sp(t)

# ╔═╡ c8c53be9-39ed-4354-90ba-27ae3f9648a6
# initial condition
x₀ = [0.0, 0.0]

# ╔═╡ 874519d1-b28d-4749-89e8-fa32fa901a15
# history function.
#   stores history of solution to handle time delay
#   initiate as v, v′ = 0 for t < 0
h(params, t) = zeros(2)

# ╔═╡ 805a0a12-3780-49b7-8296-ca76da5ea247
begin
	final_time = 25.0
	prob = DDEProblem(f, x₀, h, (0.0, final_time))
	sol = solve(prob, d_discontinuities=[0.0, 0.3])
end

# ╔═╡ e2e938e4-6cf4-4cbb-95b4-f56beeea3048
begin
	# get solution at an array of times to plot it
	t = range(0.0, final_time, length=300)
	x = sol.(t) # internal state vectors
	# shift time
	t = t .+ 0.2
	# determine y at these times
	y = [(E * x_i)[1] for x_i in x] # y

	figure()
	plot(t, y_sp.(t), color="C1", label=L"$y_{sp}(t)$", linestyle="--")
	plot(t, y, color="C0", label=L"$y(t)$")
	plot([-1.0, 0.2], [0, 0], color="C0")
	plot([-1.0, 0.2], [0, 0], color="C1", linestyle="--")
	xlabel(L"$t$")
	ylabel(L"$y(t)$")
	legend()
	tight_layout()
	savefig("../images/time_delay_response.png", format="png", dpi=300)
	gcf()
end

# ╔═╡ 117c23bd-a09b-43f1-a0b1-1a031246c1ca
y

# ╔═╡ Cell order:
# ╠═b926bbf4-ab68-11eb-136c-952f90c757df
# ╠═ab5e06ed-f1f3-41a3-99e2-2affda673af8
# ╠═1b8c79c1-d663-42dc-beda-29a881dcbd2a
# ╠═3f443782-32d8-44f0-baf3-41645d949573
# ╠═b9aa78a5-d7bb-4bd8-a38d-46ef023670ca
# ╠═c8c53be9-39ed-4354-90ba-27ae3f9648a6
# ╠═874519d1-b28d-4749-89e8-fa32fa901a15
# ╠═805a0a12-3780-49b7-8296-ca76da5ea247
# ╠═e2e938e4-6cf4-4cbb-95b4-f56beeea3048
# ╠═117c23bd-a09b-43f1-a0b1-1a031246c1ca

max_stars_repo_name: SimonEnsemble/blah

max_issues_repo_name: SimonEnsemble/blah

max_forks_repo_name: SimonEnsemble/blah


text:
abstract type AbstractOperationsModel end

abstract type EconomicDispatch <: AbstractOperationsModel end

abstract type UnitCommitment <: AbstractOperationsModel end

abstract type CustomModel <: AbstractOperationsModel end

mutable struct PowerOperationModel{ M<:AbstractOperationsModel, T<:NetworkModel, S<:Union{Nothing,Array{NamedTuple{(:service, :formulation), Tuple{D,DataType}}}} where {D <: PowerSystems.Service}}
    psmodel::Type{M}
    generation::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}
    demand::Union{Nothing,Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}}
    storage::Union{Nothing,Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}}
    branches::Array{NamedTuple{(:device, :formulation), Tuple{DataType,DataType}}}
    transmission::Type{T}
    services::S
    system::PowerSystems.PowerSystem
    model::JuMP.AbstractModel
    dynamics::Bool
    ptdf::Union{Nothing,PTDFArray}
end

mutable struct PowerSimulationsModel{T<:AbstractOperationsModel, R<:Dates.Period}
    name::String
    model::PowerOperationModel{T}
    steps::Int64
    periods::Int64
    resolution::R
    date_from::DateTime
    date_to::DateTime
    lookahead_periods::Int64
    lookahead_resolution::R
    dynamic_analysis::Bool
    timeseries::Dict{Any,Any}
end

 mutable struct PowerResults
    ThermalGeneration::Union{Nothing,DataFrame}
    RenewableGEneration::Union{Nothing,DataFrame}
    HydroGeneration::Union{Nothing,DataFrame}
    Storage::Union{Nothing,DataFrame}
    Load::Union{Nothing,DataFrame}
    SolverOutput::Union{Nothing,Dict}
end


text:
import Pkg
Pkg.activate(".")

using DistancePlan
using StaticArrays
import Makie

scene = Makie.Scene()

B = 3.
O = 200
S = 200
origin = SVector{2}(0.5, 0.5)

bounds = DistancePlan.Box(SVector{2}([-B -B]), SVector{2}([B B]))
R = 2 * B
sample_bounds = DistancePlan.Box(SVector{2}([-R -R]), SVector{2}([R R]))
obstacles = [DistancePlan.sample(bounds) for _ in 1:O]

DistancePlan.evt(origin, bounds, obstacles, 2)

println("Gathering ", S, " samples")
tree = DistancePlan.evt(origin, bounds, obstacles, S)
println("Built tree")

DistancePlan.draw!(bounds)
DistancePlan.draw!.(obstacles)
DistancePlan.draw!(tree)

Makie.save("random_tree.png", scene, resolution=(1000, 1000))

max_stars_repo_name: krzentner/DistancePlan

max_issues_repo_name: krzentner/DistancePlan

max_forks_repo_name: krzentner/DistancePlan


text:
"""
    AbstractTabular{N}
    AbstractSeries === AbstractTabular{1}
    AbstractTable === AbstractTabular{2}

`AbstractTabular{N}`s are flexible, `N`-dimensional storage wrappers, designed to provide
a uniform interface for working with data in various (possibly nested) data structures.

Objects of this abstract type can be seen as an extension (or combination) of both
`AbstractArray` and `Associative`. An `AbstractTabular{N}` object behaves as a generic
mapping from `N` indices to it's elements..

Unlike `AbstractArray`, no assumptions about the indices or element type is made - excepting
that the indices for a given dimension must be unique, and amongst dimensions they are
"Cartesian" (e.g. a 2D tabular object must be rectangular in shape, not a generic nested
or ragged array).

The `AbstractTabular{N}` interface requires the following functions (WIP):

 - `indices(tabular)` returns a length-`N` tuple of the indices for each dimension.
 - `getindex(tabular, inds...)` fetches an element given the `N` indices.
 - `setindex!(tabular, value, inds...)` sets an element to `value` given the `N` indices.
"""
abstract type AbstractTabular{N}; end

const AbstractTable = AbstractTabular{2}
const AbstractSeries = AbstractTabular{1}

ndims(t::AbstractTabular{N}) where {N} = N
ndims(t::Type{<:AbstractTabular{N}}) where {N} = N

@inline size(t::AbstractTabular) = map(length, indices(t))
@inline size(t::AbstractTabular, i) = length(indices(t)[i])

function summary(t::AbstractSeries)
    string("s Tabular{$N}")
end

function summary(t::AbstractTable)
    string(join(size(t), "×"), " Table")
end

function summary(t::AbstractTabular{N}) where {N}
    string(join(size(t), "×"), " Tabular{$N}")
end

"""
    Series(data...)
    Table(data...)
    Tabular{N}(data...)

Wraps input `data` in a `Tabular` of the specified dimensionality.
Indexing behavior will depend on the data type, and assumes a nested
data access pattern. The immediate data container `data` is indexed by
the final index value, such that (like `Base.Array`) the data "closest"
in memory is represented by the first index, and "furthest" in memory
by the final index. By default, multidimensional arrays will preserve
their index order.

If you wish the indices applied in a different order, see `PermutedDimsTabular`
or consider using `transpose(table)` / `table'`.
"""
struct Tabular{N, D} <: AbstractTabular{N}
    data::D
end

const Table = Tabular{2}
const Series = Tabular{1}

@inline Tabular{N}(x::D) where {N,D} = Tabular{N,D}(x)
@inline Tabular{N}(x...) where {N} = Tabular{N}(x)
@inline Tabular{N}(x::Pair) where {N} = Tabular{N}((x,))

get(t::Tabular) = t.data


text:
"""
`module  JAC.ParityNonConservation`
    ... a submodel of JAC that contains all methods for computing P-odd and/or T-odd amplitudes between two
        bound-state levels.
"""
module ParityNonConservation

    using Printf,  ..Basics, ..Defaults, ..InteractionStrength, ..ManyElectron, ..Nuclear, ..Radial

    """
    `ParityNonConservation.schiffMomentAmplitude(finalLevel::Level, initialLevel::Level, nm::Nuclear.Model,
                                                     grid::Radial.Grid; display::Bool=false)`
        ... to compute the Schiff moment amplitude  <alpha_f J_f || H^(Schiff) || alpha_i J_i>  for the given final
            and initial level and for the nuclear density as given by the nuclear model. A value::ComplexF64 is returned.
    """
    function schiffMomentAmplitude(finalLevel::Level, initialLevel::Level, nm::Nuclear.Model, grid::Radial.Grid; display::Bool=false)
        nf = length(finalLevel.basis.csfs);    ni = length(initialLevel.basis.csfs)
        printstyled("Compute Schiff moment  matrix of dimension $nf x $ni in the final- and initial-state bases " *
                    "for the transition [$(initialLevel.index)- $(finalLevel.index)] ... ", color=:light_green)
        matrix = zeros(ComplexF64, nf, ni)
        #
        for  r = 1:nf
            for  s = 1:ni
                ##x wa = compute("angular coefficients: 1-p, Grasp92", 0, 1, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                # Calculate the spin-angular coefficients
                if  Defaults.saRatip()
                    waR = Basics.compute("angular coefficients: 1-p, Grasp92", 0, 1, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                    wa  = waR
                end
                if  Defaults.saGG()
                    subshellList = initialLevel.basis.subshells
                    opa = SpinAngular.OneParticleOperator(1, plus, true)
                    waG = SpinAngular.computeCoefficients(opa, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s], subshellList)
                    wa  = waG
                end
                if  Defaults.saRatip() && Defaults.saGG() && true
                    if  length(waR) != 0     println("\n>> Angular coeffients from GRASP/MCT   = $waR ")    end
                    if  length(waG) != 0     println(  ">> Angular coeffients from SpinAngular = $waG ")    end
                end
                #
                for  coeff in wa
                    tamp  = InteractionStrength.schiffMoment(finalLevel.basis.orbitals[coeff.a], initialLevel.basis.orbitals[coeff.b], nm, grid)
                    matrix[r,s] = matrix[r,s] + coeff.T * tamp
                end
            end
        end
        printstyled("done. \n", color=:light_green)
        amplitude = transpose(finalLevel.mc) * matrix * initialLevel.mc
        #
        if  display
            sa = @sprintf("%.8e", amplitude.re) * @sprintf("%.8e", amplitude.im)
            println("   Schiff moment amplitude:   "                                                                               *
                    "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || H^(Schiff) ($(nm.model)) ||"   *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            printSummary, iostream = Defaults.getDefaults("summary flag/stream")
            if  printSummary
                println(iostream, "   Schiff moment amplitude:   "                                                                               *
                                  "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || H^(Schiff) ($(nm.model)) ||"   *
                                  " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            end
        end

        return( amplitude )
    end

    """
    `ParityNonConservation.anapoleMomentAmplitude(finalLevel::Level, initialLevel::Level, grid::Radial.Grid; display::Bool=false)`
        ... to compute the anapole moment amplitude <(alpha_f J_f, kappa) J_i || H^(anapole) || alpha_i J_i>
            for the given final and initial level. A value::ComplexF64 is returned.
    """
    function anapoleMomentAmplitude(finalLevel::Level, initialLevel::Level, grid::Radial.Grid; display::Bool=false)
        nf = length(finalLevel.basis.csfs);    ni = length(initialLevel.basis.csfs)
        printstyled("Compute anapole moment matrix of dimension $nf x $ni in the final- and initial-state bases " *
                                       "for the transition [$(initialLevel.index)- $(finalLevel.index)] ... ", color=:light_green)
        matrix = zeros(ComplexF64, nf, ni)
        #
        #
        printstyled("done. \n", color=:light_green)
        #
        amplitude = 3.0 + 3.0*im
        #
        if  display
                sa = @sprintf("%.8e", amplitude.re) * @sprintf("%.8e", amplitude.im)
            println("   Anapole moment amplitude:  "                                                                  *
                    "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || H^(anapole) ||"   *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            printSummary, iostream = Defaults.getDefaults("summary flag/stream")
            if  printSummary
                println(iostream, "   Anapole moment amplitude:  "                                                                  *
                                  "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || H^(anapole) ||"   *
                                  " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            end
        end

        return( amplitude )
    end

    """
    `ParityNonConservation.weakChargeAmplitude(finalLevel::Level, initialLevel::Level, nm::Nuclear.Model,
                                                   grid::Radial.Grid; display::Bool=false)`
        ... to compute the weak-charge amplitude  <alpha_f J_f || H^(weak-charge) || alpha_i J_i>  for the given final and
            initial level and for the nuclear density as given by the nuclear model. A value::ComplexF64 is returned.
    """
    function weakChargeAmplitude(finalLevel::Level, initialLevel::Level, nm::Nuclear.Model, grid::Radial.Grid; display::Bool=false)
        nf = length(finalLevel.basis.csfs);    ni = length(initialLevel.basis.csfs)
        printstyled("Compute weak -charge matrix of dimension $nf x $ni in the final- and initial-state bases " *
                    "for the transition [$(initialLevel.index)- $(finalLevel.index)] ... ", color=:light_green)
        matrix = zeros(ComplexF64, nf, ni)
        #
        for  r = 1:nf
            for  s = 1:ni
                ##x wa = compute("angular coefficients: 1-p, Grasp92", 0, 1, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                # Calculate the spin-angular coefficients
                if  Defaults.saRatip()
                    waR = Basics.compute("angular coefficients: 1-p, Grasp92", 0, 1, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                    wa  = waR
                end
                if  Defaults.saGG()
                    subshellList = initialLevel.basis.subshells
                    opa = SpinAngular.OneParticleOperator(1, plus, true)
                    waG = SpinAngular.computeCoefficients(opa, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s], subshellList)
                    wa  = waG
                end
                if  Defaults.saRatip() && Defaults.saGG() && true
                    if  length(waR) != 0     println("\n>> Angular coeffients from GRASP/MCT   = $waR ")    end
                    if  length(waG) != 0     println(  ">> Angular coeffients from SpinAngular = $waG ")    end
                end
                #
                for  coeff in wa
                    tamp  = InteractionStrength.weakCharge(finalLevel.basis.orbitals[coeff.a], initialLevel.basis.orbitals[coeff.b], nm, grid)
                    matrix[r,s] = matrix[r,s] + coeff.T * tamp
                end
            end
        end
        printstyled("done. \n", color=:light_green)
        amplitude = transpose(finalLevel.mc) * matrix * initialLevel.mc
        #
        if  display
            sa = @sprintf("%.5e", amplitude.re) * "  " * @sprintf("%.5e", amplitude.im)
            println("   weak-charge amplitude:     "                                                         *
                    "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || H^(weak-charge) ($(nm.model)) ||"   *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            printSummary, iostream = Defaults.getDefaults("summary flag/stream")
            if  printSummary
                println(iostream, "   weak-charge amplitude:     "                                                         *
                                  "< level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] || H^(weak-charge) ($(nm.model)) ||"   *
                                  " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = " * sa)
            end
        end

        return( amplitude )
    end

end # module


text:
function test_primal()
    n    = 100
    dim  = 2
    X    = rand(n, dim)
    y    = rand(n) .>= 0.75
    data = Primal(X, y)

    @testset "primal dataset" begin
        @test typeof(data) <: AbstractData
        @test data.dim  == dim
        @test data.n    == n
        @test data.npos == sum(y)
        @test data.nneg == sum(.~y)
        @test data.n    == data.npos + data.nneg
        @test length(data.ind_pos) == data.npos
        @test length(data.ind_neg) == data.nneg
    end
end

function test_dual()
    models = [PatMat(0.5), TopPushK(5), TopPush()]

    @testset "$(nameof(typeof(model)))" for model in models
        test_dual(model)
    end
end

function test_dual(model::AbstractModel)
    n     = 100
    dim   = 2
    Xtrain = rand(n, dim)
    ytrain = rand(n) .>= 0.75

    Xvalid = rand(3*n, dim)
    yvalid = rand(3*n) .>= 0.75

    Xtest  = rand(2*n, dim)

    data_train = Dual(model, Xtrain, ytrain)
    data_valid = Dual(model, Xtrain, ytrain, Xvalid, yvalid)
    data_test  = Dual(model, Xtrain, ytrain, Xtest)

    @testset "dual train dataset" begin
        @test typeof(data_train) <: Dual{<:DTrain}
        @test issymmetric(data_train.K)
        test_dual(data_train)
    end

    @testset "dual validation dataset" begin
        @test typeof(data_valid) <: Dual{<:DValidation}
        @test size(data_valid.K, 2) == data_valid.n
        test_dual(data_valid)
    end

    @testset "dual test dataset" begin
        @test typeof(data_test) <: Dual{<:DTest}
        test_dual(data_test)
        @test size(data_test.K, 2) == data_test.n
    end
end

function test_dual(data::Dual{<:Union{DTrain, DValidation}})
    @test length(data.ind_α) == data.nα
    @test length(data.ind_β) == data.nβ
    @test size(data.K, 1)    == data.nα + data.nβ
    @test length(data.type.ind_pos) == data.type.npos
    @test length(data.type.ind_neg) == data.type.nneg
end

function test_dual(data::Dual{<:DTest})
    @test length(data.ind_α) == data.nα
    @test length(data.ind_β) == data.nβ
    @test size(data.K, 1)    == data.nα + data.nβ
end


text:
using Javis, Random, LinearAlgebra

function ground(args...)
    background("black")
    sethue("white")
end

function ptdraw(point, vertices)
    return (args...) -> begin
        maxnorm = norm(vertices[1] - O)
        colors = norm.(point .- vertices)
        colors = (maxnorm .- colors) .* (colors .== minimum(colors)) ./ maxnorm
        sethue(colors...)
        circle(point, 2, :fill)
    end
end

AL = 500
video = Video(500, 500)
Background(1:AL, ground)

# Main triangle
pts = ngon(Point(0, 40), 240, 3, deg2rad(30), vertices = true)
dpts = ngon(Point(0, 40), 245, 3, deg2rad(30), vertices = true) # draw a slightly bigger triangle to avoid drawing points on it.
Object(1:AL, (args...) -> poly(dpts, :stroke, close = true))

dots = []
points = [O]
for i in 1:(AL - 3)
    curr = rand(pts)
    nx = (curr.x + points[end].x) / 2
    ny = (curr.y + points[end].y) / 2
    push!(points, Point(nx, ny))
    push!(dots, Object(i:AL, ptdraw(points[end], pts)))
    Object(i:(i + 3), JLine(points[end], points[end - 1], color = "white", linewidth = 0.5))
end

render(video, pathname = "chaosgame.gif")


text:
export topology_id,restrain_geometry

# returns selected id of a given topology
function topology_id(t::AbstractTopology,k::Int64,i::Int64)
  if i<=nid(t) & i>0
    return t.attributes[k,i]
  else
    return void
  end
end

# returns all selected id of a given topology
function topology_id(t::AbstractTopology,i::Int64)
  if i<=nid(t) & i>0
    return t.attributes[:,i]
  else
    return void
  end
end

# reduces a vertex topology by restraining to the first N vertices
# not such a nice implementation, since I need to access the Geometry directly
function restrain_geometry(g::Geometry,N::Int64)
  xyz  = g.vertices[1:N,:]
  vvv  = VertexTopology(g.topology.elements[1:N,:],g.topology.attributes[1:N,:])
  return Geometry{dim(g)}(xyz,vvv)
end

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
# Create a function to convert a vector of vectors (e.g., list of list from YAML) to a matrix.
function Base.convert(::Type{Matrix{T}}, d::Vector{Vector{T}}) where {T}
    return hcat(d...)
end

setup(file_name::String) = setup(Scenario(), file_name);
setup(specs::Dict)       = setup(Scenario(), specs);

# Create a function to turn a JSON file name into a dictionary and then load that dictionary.
function setup(model, file_name::String, context::Module = @__MODULE__)
    println("Loading scenario specifications from '$file_name' in $context.")
    # The `do` pattern here take care of closing the file if anything goes wrong.
    if endswith(lowercase(file_name), ".json")
        specs = open(file_name, "r") do inputs
            JSON.parse(inputs)
        end
    elseif endswith(lowercase(file_name), ".yaml")
        specs = YAML.load_file(file_name)
    else
        error("I don't know how to load that file type.")
    end
    #display(specs)
    return setup(model, specs, context)
end

# setup fields matching keys with values from the dictionary.
function setup(model, specs::Dict, context::Module = @__MODULE__, spacing::String = "")

    # If this requires that we put in a new model, do so now. Afterwards, we'll skip the "model" field.
    if haskey(specs, "model")

        println(spacing, "Making new model (", specs["model"], ").")

        try

            # Turn the model text into an expression.
            println(spacing, "Turning ", specs["model"], " into an expression.")
            model_expr = parse(specs["model"])

            # Find the right place to evaluate this thing.
            if isa(model_expr, Symbol) && !isdefined(context, model_expr)
                mod = Overdot
            else
                mod = context
            end
            println(spacing, "Evaluating expression in ", mod)

            # We evaluate the expression to give a function, which ought to construct the model.
            model_constructor = Core.eval(mod, parse(specs["model"]))

            # Run the model constructor without arguments.
            model = model_constructor()

            println(spacing, "Created a default ", specs["model"], ".")

        catch err
            if isa(err, UndefVarError)
                println(spacing, "Model creation (", specs["model"], ") failed.???")
            else
                rethrow(err)
            end
            return nothing
        end

    end

    # Um, how did this happen?
    if model == nothing
        error("The model contained nothing.")
    end

    # Let's now try to map each key in the dictionary to a field of the current model.
    println(spacing, "Looking through fields referenced in the Dict...")
    for (key, value) in specs

        # Skip descriptions and model specs.
        # if contains(==, ["description", "model"], lowercase(key))
        if any((y->begin ==(y, lowercase(key)) end), ["description", "model"])
            continue
        end

        println(spacing, "Processing '$key'.")

        # Get the symbol for the field and make sure it exists.
        field = parse(key)
        if isdefined(model, field)

            # If the specs on the right are a dictionary, go deeper.
            if isa(value, Dict)

                println(spacing, "The '$key' field is a dictionary. Delving into it.")
                setfield!(model, field,
                setup(getfield(model, field), specs[key], context, spacing * "  "))

            # Otherwise, if it's a Vector, blow it away and replace it with what's in the specs.
            # A vector necessarily constructs everything from scratch, so they must
            # all specify which model to use.
            elseif isa(getfield(model, field), Vector) && !(eltype(getfield(model, field)) <: Number)

                println(spacing, "$field is a vector ($(typeof(getfield(model, field))))!")

                x = Vector{eltype(getfield(model, field))}(undef, length(value))
                for k = 1:length(value)
                    # try # We put the try here so that the whole assignment gets skipped.
                        x[k] = setup(nothing, value[k], context, spacing * "  ")
                    # catch err
                        # display(err)
                        # println(spacing, "Attempting to carry on.")
                    # end
                end
                setfield!(model, field, x)

                # Check for duplicated model names.
                model_names = [element.name for element in x]
                unique_model_names = unique(model_names)
                counts = [count(n -> n == model_name, model_names) for model_name in unique_model_names]
                if any(counts .> 1)
                    for model_name in unique_model_names[counts .> 1]
                        warn("Duplicate model name: " * model_name * ".")
                    end
                end

            # Otherwise, if this field is a function, make sure the right hand side evaluates to a function.
            elseif isa(getfield(model, field), Function) && isa(value, AbstractString)

                println("Trying to interpret the function: $value")
                converted = Core.eval(context, parse(value))
                setfield!(model, field, converted)

            # Otherwise, it's presumed to be a "shallow" type.
            else

                # Set the field of the struct.
                println(spacing, "Overwriting $key ($(typeof(getfield(model, field)))) with $(typeof(value))).")
                converted = convert_it(context, getfield(model, field), value, spacing)
                setfield!(model, field, converted)

            end

        # If the field isn't part of the struct, throw an error.
        else
            # error("Sorry, but the field '$key' isn't part of the $(typeof(model)) type.")
            println(spacing, "Sorry, but the field '$key' isn't part of the $(typeof(model)) type. (Skipping.)")
        end

    end # for key, value...

    println(spacing, "Done.")
    return model

end

function convert_it(context, target, value, spacing)

    # Try to convert the value.
    target_type = typeof(target)
    source_type = typeof(value)

    # If they're already the same type, assign.
    if source_type == target_type

        println(spacing, "Types match. Assigning directly.")
        # converted_value = value
        converted_value = value

    # Otherwise, if there's a conversion from the source type to the target
    # type, use it.
    elseif hasmethod(convert, Tuple{Type{target_type}, source_type})

        # Convert the value in the dictionary to the type contained in the
        # struct.
        println(spacing, "Converting $source_type to $target_type.")
        converted_value = convert(target_type, value)

    # If it's a string (and the target type isn't), see if we can evaluate it
    # into something useful.
    elseif isa(value, AbstractString)

        println(spacing, "Evaluating $value.")
        evaluated = Core.eval(context, parse(value))
        converted_value = convert_it(context, target, evaluated, spacing)

    # Otherwise, copy it and hope for the best.
    else
        println(spacing, "Copying")
        converted_value = deepcopy(value)
    end

    return converted_value

end


text:
"""
    LikelihoodRatioTest

Results of MixedModels.likelihoodratiotest

## Fields
* `formulas`: Vector of model formulae
* `models`: NamedTuple of the `dof` and `deviance` of the models
* `tests`: NamedTuple of the sequential `dofdiff`, `deviancediff`,
           and resulting `pvalues`

## Properties
* `deviance` : note that this is actually -2 log likelihood for linear models
               (i.e. without subtracting the constant for a saturated model)
* `pvalues`

"""
struct LikelihoodRatioTest
    formulas::AbstractVector{String}
    models::NamedTuple{(:dof, :deviance)}
    tests::NamedTuple{(:dofdiff, :deviancediff, :pvalues)}
    linear::Bool
end

function Base.propertynames(lrt::LikelihoodRatioTest, private::Bool=false)
    return (:deviance, :formulas, :models, :pvalues, :tests)
end

function Base.getproperty(lrt::LikelihoodRatioTest, s::Symbol)
    if s == :dof
        lrt.models.dof
    elseif s == :deviance
        lrt.models.deviance
    elseif s == :pvalues
        lrt.tests.pvalues
    elseif s == :formulae
        lrt.formulas
    else
        getfield(lrt, s)
    end
end

# backward syntactic but not type compatibility
Base.getindex(lrt::LikelihoodRatioTest, s::Symbol) = getfield(lrt, s)

"""
    likelihoodratiotest(m::MixedModel...)
    likelihoodratiotest(m0::LinearModel, m::MixedModel...)
    likelihoodratiotest(m0::GeneralizedLinearModel, m::MixedModel...)
    likelihoodratiotest(m0::TableRegressionModel{LinearModel}, m::MixedModel...)
    likelihoodratiotest(m0::TableRegressionModel{GeneralizedLinearModel}, m::MixedModel...)

Likeihood ratio test applied to a set of nested models.

!!! note
    The nesting of the models is not checked.  It is incumbent on the user
    to check this. This differs from `StatsModels.lrtest` as nesting in
    mixed models, especially in the random effects specification, may be non obvious.

!!! note
    For comparisons between mixed and non-mixed models, the deviance for the non-mixed
    model is taken to be -2 log likelihood, i.e. omitting the additive constant for the
    fully saturated model. This is in line with the computation of the deviance for mixed
    models.

This functionality may be deprecated in the future in favor of `StatsModels.lrtest`.
"""
function likelihoodratiotest(m::MixedModel...)
    _iscomparable(m...) ||
        throw(ArgumentError("""Models are not comparable: are the objectives, data
                               and, where appropriate, the link and family the same?
              """))

    m = collect(m)   # change the tuple to an array
    dofs = dof.(m)
    formulas = String.(Symbol.(getproperty.(m, :formula)))
    ord = sortperm(dofs)
    dofs = dofs[ord]
    formulas = formulas[ord]
    devs = objective.(m)[ord]
    dofdiffs = diff(dofs)
    devdiffs = .-(diff(devs))
    pvals = map(zip(dofdiffs, devdiffs)) do (dof, dev)
        if dev > 0
            ccdf(Chisq(dof), dev)
        else
            NaN
        end
    end

    return LikelihoodRatioTest(
        formulas,
        (dof=dofs, deviance=devs),
        (dofdiff=dofdiffs, deviancediff=devdiffs, pvalues=pvals),
        first(m) isa LinearMixedModel,
    )
end

_formula(::Union{LinearModel,GeneralizedLinearModel}) = "NA"
function _formula(x::TableRegressionModel{<:Union{LinearModel,GeneralizedLinearModel}})
    return String(Symbol(x.mf.f))
end

# for GLMMs we're actually looking at the deviance and additive constants are comparable
# (because GLM deviance is actually part of the GLMM deviance computation)
# for LMMs, we're always looking at the "deviance scale" but without the additive constant
# for the fully saturated model
function _criterion(
    x::Union{GeneralizedLinearModel,TableRegressionModel{<:GeneralizedLinearModel}}
)
    return deviance(x)
end
function _criterion(x::Union{LinearModel,TableRegressionModel{<:LinearModel}})
    return -2 * loglikelihood(x)
end

function likelihoodratiotest(
    m0::Union{
        TableRegressionModel{<:Union{LinearModel,GeneralizedLinearModel}},
        LinearModel,
        GeneralizedLinearModel,
    },
    m::MixedModel...,
)
    _iscomparable(m0, first(m)) ||
        throw(ArgumentError("""Models are not comparable: are the objectives, data
                               and, where appropriate, the link and family the same?
                            """))
    lrt = likelihoodratiotest(m...)
    devs = pushfirst!(lrt.deviance, _criterion(m0))
    formulas = pushfirst!(lrt.formulas, _formula(m0))
    dofs = pushfirst!(lrt.models.dof, dof(m0))
    devdiffs = pushfirst!(lrt.tests.deviancediff, devs[1] - devs[2])
    dofdiffs = pushfirst!(lrt.tests.dofdiff, dofs[2] - dofs[1])

    df, dev = first(dofdiffs), first(devdiffs)
    p = dev > 0 ? ccdf(Chisq(df), dev) : NaN
    pvals = pushfirst!(lrt.tests.pvalues, p)

    return LikelihoodRatioTest(
        formulas,
        (dof=dofs, deviance=devs),
        (dofdiff=dofdiffs, deviancediff=devdiffs, pvalues=pvals),
        lrt.linear,
    )
end

function Base.show(io::IO, ::MIME"text/plain", lrt::LikelihoodRatioTest)
    println(io, "Model Formulae")

    for (i, f) in enumerate(lrt.formulas)
        println(io, "$i: $f")
    end

    # the following was adapted from StatsModels#162
    # from nalimilan
    Δdf = lrt.tests.dofdiff
    Δdev = lrt.tests.deviancediff

    nc = 6
    nr = length(lrt.formulas)
    outrows = Matrix{String}(undef, nr + 1, nc)

    outrows[1, :] = [
        "", "model-dof", lrt.linear ? "-2 logLik" : "deviance", "χ²", "χ²-dof", "P(>χ²)"
    ] # colnms

    outrows[2, :] = [
        "[1]", string(lrt.dof[1]), Ryu.writefixed(lrt.deviance[1], 4), " ", " ", " "
    ]

    for i in 2:nr
        outrows[i + 1, :] = [
            "[$i]",
            string(lrt.dof[i]),
            Ryu.writefixed(lrt.deviance[i], 4),
            Ryu.writefixed(Δdev[i - 1], 4),
            string(Δdf[i - 1]),
            string(StatsBase.PValue(lrt.pvalues[i - 1])),
        ]
    end
    colwidths = length.(outrows)
    max_colwidths = [maximum(view(colwidths, :, i)) for i in 1:nc]
    totwidth = sum(max_colwidths) + 2 * 5

    println(io, '─'^totwidth)

    for r in 1:(nr + 1)
        for c in 1:nc
            cur_cell = outrows[r, c]
            cur_cell_len = length(cur_cell)

            padding = " "^(max_colwidths[c] - cur_cell_len)
            if c > 1
                padding = "  " * padding
            end

            print(io, padding)
            print(io, cur_cell)
        end
        print(io, "\n")
        r == 1 && println(io, '─'^totwidth)
    end
    print(io, '─'^totwidth)

    return nothing
end

Base.show(io::IO, lrt::LikelihoodRatioTest) = Base.show(io, MIME"text/plain"(), lrt)

function _iscomparable(m::LinearMixedModel...)
    isconstant(getproperty.(getproperty.(m, :optsum), :REML)) || throw(
        ArgumentError(
            "Models must all be fit with the same objective (i.e. all ML or all REML)"
        ),
    )

    if any(getproperty.(getproperty.(m, :optsum), :REML))
        isconstant(coefnames.(m)) || throw(
            ArgumentError(
                "Likelihood-ratio tests for REML-fitted models are only valid when the fixed-effects specifications are identical",
            ),
        )
    end

    isconstant(nobs.(m)) ||
        throw(ArgumentError("Models must have the same number of observations"))

    return true
end

function _samefamily(
    ::GeneralizedLinearMixedModel{<:AbstractFloat,S}...
) where {S<:Distribution}
    return true
end
_samefamily(::GeneralizedLinearMixedModel...) = false

function _iscomparable(m::GeneralizedLinearMixedModel...)
    # TODO: test that all models are fit with same fast/nAGQ option?
    _samefamily(m...) || throw(ArgumentError("Models must be fit to the same distribution"))

    isconstant(string.(Link.(m))) ||
        throw(ArgumentError("Models must have the same link function"))

    isconstant(nobs.(m)) ||
        throw(ArgumentError("Models must have the same number of observations"))

    return true
end

"""
    isnested(m1::MixedModel, m2::MixedModel; atol::Real=0.0)
Indicate whether model `m1` is nested in model `m2`, i.e. whether
`m1` can be obtained by constraining some parameters in `m2`.
Both models must have been fitted on the same data. This check
is conservative for `MixedModel`s and may reject nested models with different
parameterizations as being non nested.
"""
function StatsModels.isnested(m1::MixedModel, m2::MixedModel; atol::Real=0.0)
    try
        _iscomparable(m1, m2)
    catch e
        @error e.msg
        false
    end || return false

    # check that the nested fixef are a subset of the outer
    all(in.(coefnames(m1), Ref(coefnames(m2)))) || return false

    # check that the same grouping vars occur in the outer model
    grpng1 = fname.(m1.reterms)
    grpng2 = fname.(m2.reterms)

    all(in.(grpng1, Ref(grpng2))) || return false

    # check that every intercept/slope for a grouping var occurs in the
    # same grouping
    re1 = Dict(fname(re) => re.cnames for re in m1.reterms)
    re2 = Dict(fname(re) => re.cnames for re in m2.reterms)

    all(all(in.(val, Ref(re2[key]))) for (key, val) in re1) || return false

    return true
end

function _iscomparable(
    m1::TableRegressionModel{<:Union{LinearModel,GeneralizedLinearModel}}, m2::MixedModel
)
    _iscomparable(m1.model, m2) || return false

    # check that the nested fixef are a subset of the outer
    all(in.(coefnames(m1), Ref(coefnames(m2)))) || return false

    return true
end

# GLM isn't nested with in LMM and LM isn't nested within GLMM
_iscomparable(m1::Union{LinearModel,GeneralizedLinearModel}, m2::MixedModel) = false

function _iscomparable(m1::LinearModel, m2::LinearMixedModel)
    nobs(m1) == nobs(m2) || return false

    # XXX This reaches into the internal structure of GLM
    size(m1.pp.X, 2) <= size(m2.X, 2) || return false

    _isnested(m1.pp.X, m2.X) || return false

    !m2.optsum.REML ||
        throw(ArgumentError("REML-fitted models cannot be compared to linear models"))

    return true
end

function _iscomparable(m1::GeneralizedLinearModel, m2::GeneralizedLinearMixedModel)
    nobs(m1) == nobs(m2) || return false

    size(modelmatrix(m1), 2) <= size(modelmatrix(m2), 2) || return false

    _isnested(modelmatrix(m1), modelmatrix(m2)) || return false

    Distribution(m1) == Distribution(m2) ||
        throw(ArgumentError("Models must be fit to the same distribution"))

    Link(m1) == Link(m2) || throw(ArgumentError("Models must have the same link function"))

    return true
end

"""
    _isnested(x::AbstractMatrix, y::AbstractMatrix; atol::Real=0.0)

Test whether the column span of `x` is a subspace of (nested within)
the column span of y.

The nesting of the column span of the fixed-effects model matrices is a necessary,
but not sufficient condition for a linear model (whether mixed-effects or not)
to be nested within a linear mixed-effects model.

!!! note
    The `rtol` argument is an internal threshold and not currently
    compatible with the `atol` argument of `StatsModels.isnested`.
"""
function _isnested(x::AbstractMatrix, y::AbstractMatrix; rtol=1e-8, ranktol=1e-8)

    # technically this can return false positives if x or y
    # are rank deficient, but either they're rank deficient
    # in the same way (b/c same data) and we don't care OR
    # it's not the same data/fixef specification and we're
    # extra conservative
    size(x, 2) <= size(y, 2) || return false

    qy = qr(y).Q

    qrx = pivoted_qr(x)
    dvec = abs.(diag(qrx.R))
    fdv = first(dvec)
    cmp = fdv * ranktol
    r = searchsortedlast(dvec, cmp; rev=true)

    p = qy' * x

    nested = map(eachcol(p)) do col
        # if set Julia 1.6 as the minimum, we can use last(col, r)
        top = @view col[firstindex(col):(end - r - 1)]
        tail = @view col[(end - r):end]
        return norm(tail) / norm(top) < rtol
    end

    return all(nested)
end


text:
if basename(pwd()) == "aoc"
    cd("2016/20")
end

function loadranges(filename::AbstractString)
    a, ranges... = sort!(
        [range.(parse.(Int, split(line, '-'))...) for line in eachline(filename)],
        by = first
    )
    len = -1
    while len != length(ranges)
        seen = UnitRange[a]
        len = length(ranges)
        for b in ranges
            if last(a) + 1 >= first(b)
                a = first(a):max(last(a), last(b))
            else
                push!(seen, a)
                a = b
            end
        end
        push!(seen, a)
        a, ranges... = seen
    end
    ranges
end

part1(filename::AbstractString) = last(first(loadranges(filename))) + 1

@assert part1("example.txt") == 3
part1("input.txt")

function part2(filename::AbstractString, limit = Int(typemax(UInt32)) + 1)
    limit - sum(length.(loadranges(filename)))
end

@assert part2("example.txt", 10) == 2
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module CardThemes

# Notes for adding new card theme
#
# Steps to write new card theme:
#   1. specialize `cardtheme`, it should return two items:
#     * `templates`: a dict with following fields:
#        * "card": a Mustache template
#        * "section": a Mustache template
#     * `stylesheet_path`: absolute filepath to the stylesheet.
#   2. add the new theme name to `themelists`
#
# available placeholders for section template:
#   * `cards`: the card table/list/whatever contents
#   * `description`: section description contents
#
# available placeholders for card template:
#  * `description`: card description
#  * `covername`: the name of the card cover image
#  * `id`: reference id that reads by Documenter
#  * `title`: card title
#
using Mustache

"""
A list of DemoCards theme presets
"""
const themelists = ["grid", "list"]

# TODO: don't hardcode this
const max_coversize = (220, 200)

# add themes
include("grid/grid.jl")
include("list/list.jl")

"""
    cardtheme(theme = "grid";
              root = "<current-directory>",
              destination = "democards") -> templates, stylesheet_path

For given theme, return the templates and path to stylesheet.

`root` and `destination` should have the same value to that passed to `makedemos`.

Available themes are:
\n - "$(join(themelists, "\"\n  - \""))"
"""
function cardtheme(theme::AbstractString="grid";
                   root::AbstractString = Base.source_dir(),
                   destination::String = "democards")
    templates, src_stylesheet_path = cardtheme(Val(Symbol(theme)))

    # a copy is needed because Documenter only support relative path inside root
    absolute_root = joinpath(root, "src", destination)
    filename = "$(theme)theme.css"
    out_stylesheet_path = joinpath(absolute_root, filename)
    isdir(absolute_root) || mkpath(absolute_root)
    cp(src_stylesheet_path, out_stylesheet_path; force=true)

    return templates, joinpath(destination, filename)
end

# fallback if the theme name doesn't exist
function cardtheme(theme::Val{T}) where T
    themes_str = join(themelists, ", ")
    error("unrecognized card theme: $T\nAvaiable theme options are: $themes_str")
end

export cardtheme, max_coversize

end