text:
const dummy_function = ()->nothing

num_linkconstraints(modeledge::OptiEdge) = length(modeledge.linkconstraints)
function _caching_optimizer(modelnode::OptiNode)
    if isa(modelnode.model.moi_backend,MOIU.CachingOptimizer)
        return modelnode.model.moi_backend
    else
        return modelnode.model.moi_backend.optimizer
    end
end
moi_optimizer(modelnode::OptiNode) = _caching_optimizer(modelnode).optimizer.model
_caching_optimizer(model::Any) = model.moi_backend
function set_g_link!(linkedge::OptiEdge,l,gl,gu)
    cnt = 1
    for (ind,linkcon) in linkedge.linkconstraints
        l[cnt] = 0. # need to implement dual start later
        if linkcon.set isa MOI.EqualTo
            gl[cnt] = linkcon.set.value
            gu[cnt] = linkcon.set.value
        elseif linkcon.set isa MOI.GreaterThan
            gl[cnt] = linkcon.set.lower
            gu[cnt] = Inf
        elseif linkcon.set isa MOI.LessThan
            gl[cnt] = -Inf
            gu[cnt] = linkcon.set.upper
        else
            gl[cnt] = linkcon.set.lower
            gu[cnt] = linkcon.set.upper
        end
        cnt += 1
    end
end

function hessian_lagrangian_structure(graph::OptiGraph,I,J,ninds,nnzs_hess_inds,modelnodes)
    @blas_safe_threads for k=1:length(modelnodes)
        isempty(nnzs_hess_inds[k]) && continue
        offset = ninds[k][1]-1
        II = view(I,nnzs_hess_inds[k])
        JJ = view(J,nnzs_hess_inds[k])
        hessian_lagrangian_structure(
            moi_optimizer(modelnodes[k]),II,JJ)
        II.+= offset
        JJ.+= offset
    end
end

function jacobian_structure(linkedge::OptiEdge,I,J,ninds,x_index_map,g_index_map)
    offset=1
    for linkcon in link_constraints(linkedge)
        offset += jacobian_structure(linkcon,I,J,ninds,x_index_map,g_index_map,offset)
    end
end
function jacobian_structure(linkcon,I,J,ninds,x_index_map,g_index_map,offset)
    cnt = 0
    for var in get_vars(linkcon)
        I[offset+cnt] = g_index_map[linkcon]
        J[offset+cnt] = x_index_map[var]
        cnt += 1
    end
    return cnt
end

function jacobian_structure(
    graph::OptiGraph,I,J,ninds,minds,pinds,nnzs_jac_inds,nnzs_link_jac_inds,
    x_index_map,g_index_map,modelnodes,linkedges)

    @blas_safe_threads for k=1:length(modelnodes)
        isempty(nnzs_jac_inds[k]) && continue
        offset_i = minds[k][1]-1
        offset_j = ninds[k][1]-1
        II = view(I,nnzs_jac_inds[k])
        JJ = view(J,nnzs_jac_inds[k])
        jacobian_structure(
            moi_optimizer(modelnodes[k]),II,JJ)
        II.+= offset_i
        JJ.+= offset_j
    end

    @blas_safe_threads for q=1:length(linkedges)
        isempty(nnzs_link_jac_inds[q]) && continue
        II = view(I,nnzs_link_jac_inds[q])
        JJ = view(J,nnzs_link_jac_inds[q])
        jacobian_structure(
            linkedges[q],II,JJ,ninds,x_index_map,g_index_map)
    end
end
function eval_objective(graph::OptiGraph,x,ninds,x_index_map,modelnodes)
    obj = Threads.Atomic{Float64}(0.)
    @blas_safe_threads for k=1:length(modelnodes)
         Threads.atomic_add!(obj,eval_objective(
             moi_optimizer(modelnodes[k]),view(x,ninds[k])))
    end
    return obj.value + eval_function(graph.objective_function,x,ninds,x_index_map)
end
function eval_objective_gradient(graph::OptiGraph,f,x,ninds,modelnodes)
    @blas_safe_threads for k=1:length(modelnodes)
        eval_objective_gradient(moi_optimizer(modelnodes[k]),view(f,ninds[k]),view(x,ninds[k]))
    end
end

function eval_function(aff::GenericAffExpr,x,ninds,x_index_map)
    function_value = aff.constant
    for (var,coef) in aff.terms
        function_value += coef*x[x_index_map[var]]
    end
    return function_value
end
function eval_constraint(linkedge::OptiEdge,c,x,ninds,x_index_map)
    cnt = 1
    for linkcon in link_constraints(linkedge)
        c[cnt] = eval_function(get_func(linkcon),x,ninds,x_index_map)
        cnt += 1
    end
end
get_func(linkcon) = linkcon.func
function eval_constraint(graph::OptiGraph,c,x,ninds,minds,pinds,x_index_map,modelnodes,linkedges)
    @blas_safe_threads for k=1:length(modelnodes)
        eval_constraint(moi_optimizer(modelnodes[k]),view(c,minds[k]),view(x,ninds[k]))
    end
    @blas_safe_threads for q=1:length(linkedges)
        eval_constraint(linkedges[q],view(c,pinds[q]),x,ninds,x_index_map)
    end
end

function eval_hessian_lagrangian(graph::OptiGraph,hess,x,sig,l,
                                 ninds,minds,nnzs_hess_inds,modelnodes)
    @blas_safe_threads for k=1:length(modelnodes)
        isempty(nnzs_hess_inds) && continue
        eval_hessian_lagrangian(moi_optimizer(modelnodes[k]),
                                view(hess,nnzs_hess_inds[k]),view(x,ninds[k]),sig,
                                view(l,minds[k]))
    end
end

function eval_constraint_jacobian(linkedge::OptiEdge,jac,x)
    offset=0
    for linkcon in link_constraints(linkedge)
        offset+=eval_constraint_jacobian(linkcon,jac,offset)
    end
end
function eval_constraint_jacobian(linkcon,jac,offset)
    cnt = 0
    for coef in get_coeffs(linkcon)
        cnt += 1
        jac[offset+cnt] = coef
    end
    return cnt
end
get_vars(linkcon) = keys(linkcon.func.terms)
get_coeffs(linkcon) = values(linkcon.func.terms)

function eval_constraint_jacobian(graph::OptiGraph,jac,x,
                                  ninds,minds,nnzs_jac_inds,nnzs_link_jac_inds,modelnodes,linkedges)
    @blas_safe_threads for k=1:length(modelnodes)
        eval_constraint_jacobian(
            moi_optimizer(modelnodes[k]),view(jac,nnzs_jac_inds[k]),view(x,ninds[k]))
    end
    @blas_safe_threads for q=1:length(linkedges)
        eval_constraint_jacobian(linkedges[q],view(jac,nnzs_link_jac_inds[q]),x)
    end
end
get_nnz_link_jac(linkedge::OptiEdge) = sum(
    length(linkcon.func.terms) for (ind,linkcon) in linkedge.linkconstraints)

struct GraphModel{T} <: AbstractNLPModel{T,Vector{T}}
    ninds::Vector{UnitRange{Int}}
    minds::Vector{UnitRange{Int}}
    pinds::Vector{UnitRange{Int}}
    nnzs_jac_inds::Vector{UnitRange{Int}}
    nnzs_hess_inds::Vector{UnitRange{Int}}
    nnzs_link_jac_inds::Vector{UnitRange{Int}}

    x_index_map::Dict #
    g_index_map::Dict #

    modelnodes::Vector{OptiNode}
    linkedges::Vector{OptiEdge}

    meta::NLPModelMeta{T, Vector{T}}
    graph::OptiGraph
    counters::NLPModelsCounters
    ext::Dict{Symbol,Any}
end

obj(nlp::GraphModel, x::AbstractVector) =  eval_objective(nlp.graph,x,nlp.ninds,nlp.x_index_map,nlp.modelnodes)
grad!(nlp::GraphModel, x::AbstractVector, f::AbstractVector) =eval_objective_gradient(nlp.graph,f,x,nlp.ninds,nlp.modelnodes)
cons!(nlp::GraphModel, x::AbstractVector, c::AbstractVector) = eval_constraint(
    nlp.graph,c,x,nlp.ninds,nlp.minds,nlp.pinds,nlp.x_index_map,nlp.modelnodes,nlp.linkedges)
function hess_coord!(nlp::GraphModel,x::AbstractVector,l::AbstractVector,hess::AbstractVector;obj_weight=1.)
    eval_hessian_lagrangian(
        nlp.graph,hess,x,obj_weight,l,nlp.ninds,nlp.minds,
        nlp.nnzs_hess_inds,nlp.modelnodes,
    )
end
function jac_coord!(nlp::GraphModel,x::AbstractVector,jac::AbstractVector)
    eval_constraint_jacobian(
        nlp.graph,jac,x,nlp.ninds,nlp.minds,nlp.nnzs_jac_inds,nlp.nnzs_link_jac_inds,
        nlp.modelnodes,nlp.linkedges,
    )
end
function hess_structure!(nlp::GraphModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
    hessian_lagrangian_structure(
        nlp.graph,I,J,nlp.ninds,nlp.nnzs_hess_inds,nlp.modelnodes,
    )
end
function jac_structure!(nlp::GraphModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
    jacobian_structure(
        nlp.graph,I,J,nlp.ninds,nlp.minds,nlp.pinds,nlp.nnzs_jac_inds,nlp.nnzs_link_jac_inds,
        nlp.x_index_map,nlp.g_index_map,nlp.modelnodes,nlp.linkedges,
    )
end

function GraphModel(graph::OptiGraph)
    modelnodes = all_nodes(graph)
    linkedges = all_edges(graph)

    for modelnode in modelnodes
        num_variables(modelnode) == 0 && error("Empty node exist! Delete the empty nodes.")
    end

    @blas_safe_threads for k=1:length(modelnodes)
        set_optimizer(modelnodes[k],Optimizer)
        if modelnodes[k].model.nlp_data !== nothing
            MOI.set(modelnodes[k].model, MOI.NLPBlock(),
                    _create_nlp_block_data(modelnodes[k].model))
            empty!(modelnodes[k].model.nlp_data.nlconstr_duals)
        end
        MOIU.attach_optimizer(modelnodes[k].model)
        MOI.initialize(moi_optimizer(modelnodes[k]).nlp_data.evaluator,[:Grad,:Hess,:Jac])
    end

    K = length(modelnodes)
    ns= [num_variables(modelnode) for modelnode in modelnodes]
    n = sum(ns)
    ns_cumsum = cumsum(ns)
    ms= [num_constraints(modelnode) for modelnode in modelnodes]
    ms_cumsum = cumsum(ms)
    m = sum(ms)

    nnzs_hess = [get_nnz_hess(moi_optimizer(modelnode)) for modelnode in modelnodes]
    nnzs_hess_cumsum = cumsum(nnzs_hess)
    nnz_hess = sum(nnzs_hess)

    nnzs_jac = [get_nnz_jac(moi_optimizer(modelnode)) for modelnode in modelnodes]
    nnzs_jac_cumsum = cumsum(nnzs_jac)
    nnz_jac = sum(nnzs_jac)

    nnzs_link_jac = [get_nnz_link_jac(linkedge) for linkedge in linkedges]
    nnzs_link_jac_cumsum = cumsum(nnzs_link_jac)
    nnz_link_jac = isempty(nnzs_link_jac) ? 0 : sum(nnzs_link_jac)

    ninds = [(i==1 ? 0 : ns_cumsum[i-1])+1:ns_cumsum[i] for i=1:K]
    minds = [(i==1 ? 0 : ms_cumsum[i-1])+1:ms_cumsum[i] for i=1:K]
    nnzs_hess_inds = [(i==1 ? 0 : nnzs_hess_cumsum[i-1])+1:nnzs_hess_cumsum[i] for i=1:K]
    nnzs_jac_inds = [(i==1 ? 0 : nnzs_jac_cumsum[i-1])+1:nnzs_jac_cumsum[i] for i=1:K]

    Q = length(linkedges)
    ps= [num_linkconstraints(modeledge) for modeledge in linkedges]
    ps_cumsum =  cumsum(ps)
    p = sum(ps)
    pinds = [(i==1 ? m : m+ps_cumsum[i-1])+1:m+ps_cumsum[i] for i=1:Q]
    nnzs_link_jac_inds =
        [(i==1 ? nnz_jac : nnz_jac+nnzs_link_jac_cumsum[i-1])+1: nnz_jac + nnzs_link_jac_cumsum[i] for i=1:Q]

    x =Vector{Float64}(undef,n)
    xl=Vector{Float64}(undef,n)
    xu=Vector{Float64}(undef,n)

    l =Vector{Float64}(undef,m+p)
    gl=Vector{Float64}(undef,m+p)
    gu=Vector{Float64}(undef,m+p)

    @blas_safe_threads for k=1:K
        set_x!(moi_optimizer(modelnodes[k]),view(x,ninds[k]),view(xl,ninds[k]),view(xu,ninds[k]))
        set_g!(moi_optimizer(modelnodes[k]),view(l,minds[k]),view(gl,minds[k]),view(gu,minds[k]))
    end

    @blas_safe_threads for q=1:Q
        set_g_link!(linkedges[q],view(l,pinds[q]),view(gl,pinds[q]),view(gu,pinds[q]))
    end

    modelmap=Dict(modelnodes[k].model=> k for k=1:K)
    x_index_map = Dict(
        # var=>ninds[modelmap[var.model]][backend(var.model).optimizer.model_to_optimizer_map[var.index].value]
        var=>ninds[modelmap[var.model]][_caching_optimizer(getnode(var)).model_to_optimizer_map[var.index].value]
        for modelnode in modelnodes for var in all_variables(modelnode))
    cnt = 0
    g_index_map = Dict(con=> m + (cnt+=1) for linkedge in linkedges for (ind,con) in linkedge.linkconstraints)

    jac_constant = true
    hess_constant = true
    for node in modelnodes
        j,h = is_jac_hess_constant(moi_optimizer(node))
        jac_constant = jac_constant & j
        hess_constant = hess_constant & h
    end

    ext = Dict{Symbol,Any}(:n=>n,:m=>m,:p=>p,:ninds=>ninds,:minds=>minds,:pinds=>pinds,
                           :linkedges=>linkedges,:jac_constant=>jac_constant,:hess_constant=>hess_constant)

    return GraphModel(
        ninds,minds,pinds,nnzs_jac_inds,nnzs_hess_inds,nnzs_link_jac_inds,
        x_index_map,g_index_map,modelnodes,linkedges,
        NLPModelMeta(
            n,
            ncon = m+p,
            x0 = x,
            lvar = xl,
            uvar = xu,
            y0 = l,
            lcon = gl,
            ucon = gu,
            nnzj = nnz_jac + nnz_link_jac,
            nnzh = nnz_hess,
            minimize = true # graph.objective_sense == MOI.MIN_SENSE
        ),
        graph,NLPModelsCounters(),
        ext
    )

end

function get_part(graph::OptiGraph,nlp::GraphModel)
    n = nlp.ext[:n]
    m = nlp.ext[:m]
    p = nlp.ext[:p]

    ninds = nlp.ninds
    minds = nlp.minds
    pinds = nlp.pinds

    ind_ineq = findall(get_lcon(nlp).!=get_ucon(nlp))
    l = length(ind_ineq)

    part = Vector{Int}(undef,n+m+l+p)

    for k=1:length(ninds)
        part[ninds[k]].=k
    end
    for k=1:length(minds)
        part[minds[k].+n.+l].=k
    end
    cnt = 0

    for linkedge in nlp.ext[:linkedges]
        for (ind,con) in linkedge.linkconstraints
            cnt+=1
            attached_node_idx = graph.node_idx_map[con.attached_node]
            part[n+l+m+cnt] = attached_node_idx != nothing ? attached_node_idx : error("All the link constraints need to be attached to a node")
        end
    end

    cnt = 0
    for q in ind_ineq
        cnt+=1
        part[n+cnt] = part[n+l+q]
    end

    return part
end

function optimize!(graph::OptiGraph; option_dict = Dict{Symbol,Any}(), kwargs...)
    graph.objective_function.constant = 0.
    nlp = GraphModel(graph)

    K = num_all_nodes(graph)
    if (haskey(kwargs,:schwarz_custom_partition) && kwargs[:schwarz_custom_partition]) ||
        (haskey(option_dict,:schwarz_custom_partition) && option_dict[:schwarz_custom_partition])

        part= get_part(graph,nlp)
        option_dict[:schwarz_part] = part
        option_dict[:schwarz_num_parts] = num_all_nodes(graph)

    elseif (haskey(kwargs,:schur_custom_partition) && kwargs[:schur_custom_partition]) ||
        (haskey(option_dict,:schur_custom_partition) && option_dict[:schur_custom_partition])

        part= get_part(graph,nlp)
        part[part.>K].=0
        option_dict[:schur_part] = part
        option_dict[:schur_num_parts] = K
    end

    option_dict[:jacobian_constant] = nlp.ext[:jac_constant]
    option_dict[:hessian_constant] = nlp.ext[:hess_constant]
    ips = InteriorPointSolver(nlp;option_dict=option_dict,kwargs...)
    result = optimize!(ips)

    graph.optimizer = ips
    graph.objective_function.constant = graph.optimizer.obj_val

    @blas_safe_threads for k=1:K
        moi_optimizer(nlp.modelnodes[k]).result = MadNLPExecutionStats(
            ips.status,view(result.solution,nlp.ninds[k]),
            ips.obj_val,
            view(result.constraints,nlp.minds[k]),
            ips.inf_du, ips.inf_pr,
            view(result.multipliers,nlp.minds[k]),
            view(result.multipliers_L,nlp.ninds[k]),
            view(result.multipliers_U,nlp.ninds[k]),
            ips.cnt.k, ips.nlp.counters,ips.cnt.total_time)
            # TODO: quick hack to specify to JuMP that the
            # model is not dirty (so we do not run in `OptimizeNotCalled`
            # exception).
            nlp.modelnodes[k].model.is_model_dirty = false
    end
end


text:
using Pkg

Pkg.activate("coveragetempenv", shared = true)

Pkg.add(PackageSpec(name = "CoverageTools"))

using CoverageTools

function main()
    help = """
    code coverage processing.

    codecov.jl <package>
    """

    root_directory = dirname(@__DIR__)

    if length(ARGS) == 0
        package_names = readdir(joinpath(root_directory, "lib"))
        push!(package_names, "Yao")
    else
        package_names = ARGS
    end
    dirs = map(package_names) do name
        return if name == "Yao"
            joinpath(root_directory, "src")
        else
            joinpath(root_directory, "lib", name, "src")
        end
    end
    pfs = mapreduce(process_folder, vcat, dirs)
    LCOV.writefile(joinpath(root_directory, "lcov.info"), pfs)
    return 0
end

main()

ext: jl


text:
"""
`module  JAC.Atomic`
    ... a submodel of JAC that contains all methods to set-up and process (simple) atomic and SCF computations as well as
        atomic representations.
"""
module Atomic

    ## using Interact
    using ..Basics, ..Radial, ..ManyElectron, ..Nuclear
    using ..Einstein, ..Hfs, ..IsotopeShift, ..PlasmaShift, ..LandeZeeman, ..AlphaVariation, ..FormFactor, ..DecayYield,
          ..MultipolePolarizibility, ..PhotoEmission, ..PhotoIonization, ..PhotoExcitation, ..PhotoRecombination,
          ..AutoIonization, ..DoubleAutoIonization, ..Dielectronic, ..ImpactExcitation, ..CoulombExcitation, ..CoulombIonization,
          ..PhotoDoubleIonization, ..PhotoExcitationFluores, ..PhotoExcitationAutoion, ..RayleighCompton, ..MultiPhotonDeExcitation,
          ..PhotoIonizationFluores, ..PhotoIonizationAutoion, ..ImpactExcitationAutoion, ..RadiativeAuger,
          ..MultiPhotonIonization, ..MultiPhotonDoubleIon, ..InternalConversion

    """
    `struct  Computation`
        ... defines a type for defining  (the model of simple) atomic computation of a single multiplet,
            including the SCF and CI as well as level properties and transition property calculations.

        + name                           ::String                          ... A name associated to the computation.
        + nuclearModel                   ::Nuclear.Model                   ... Model, charge and parameters of the nucleus.
        + grid                           ::Radial.Grid                     ... The radial grid to be used for the computation.
        + properties                     ::Array{AbstractLevelProperty,1}  ... List of atomic properties to be calculated.
        + configs                        ::Array{Configuration,1}          ... A list of non-relativistic configurations.
        + asfSettings                    ::AsfSettings
            ... Provides the settings for the SCF process and for the CI and QED calculations.
        + initialConfigs                 ::Array{Configuration,1}
            ... A list of initial-state configurations for some transition property calculation, such as radiative transition, Auger, etc.
        + initialAsfSettings             ::AsfSettings                     ... Provides the SCF and CI settings for the initial-state multiplet.
        + intermediateConfigs            ::Array{Configuration,1}          ... A list of initial-state configurations.
        + intermediateAsfSettings        ::AsfSettings                     ... Provides the SCF settings for the intermediate-state multiplet.
        + finalConfigs                   ::Array{Configuration,1}          ... A list of final-state configurations.
        + finalAsfSettings               ::AsfSettings                     ... Provides the SCF and CI settings for the final-state multiplet.
        + alphaSettings                  ::AlphaVariation.Settings         ... Settings for alpha-variation parameter calculations.
        + einsteinSettings               ::Einstein.Settings               ... Settings for Einstein coefficient calculations.
        + formSettings                   ::FormFactor.Settings             ... Settings for atomic form factor calculations.
        + hfsSettings                    ::Hfs.Settings                    ... Settings for hyperfine parameter calculations.
        + isotopeSettings                ::IsotopeShift.Settings           ... Settings for isotope shift parameter calculations.
        + plasmaSettings                 ::PlasmaShift.Settings            ... Settings for plasma-shift calculations.
        + polaritySettings               ::MultipolePolarizibility.Settings .. Settings for polarizibility calculations.
        + yieldSettings                  ::DecayYield.Settings             ... Settings for fluoresence and Auger yield calculations.
        + zeemanSettings                 ::LandeZeeman.Settings            ... Settings for Lande-Zeeman coefficient calculations.
        + process                        ::Basic.AbstractProcess
            ... An (additional) process for which the properties are to be evaluated for the given initial- and final-state configurations.
        + processSettings                ::Union{PhotoEmission.Settings, AutoIonization.Settings, PlasmaShift.AugerSettings,
                                                 PhotoIonization.Settings, PlasmaShift.PhotoSettings,
                                                 PhotoExcitation.Settings, PhotoExcitationAutoion.Settings, PhotoRecombination.Settings,
                                                 DoubleAutoIonization.Settings, PhotoDoubleIonization.Settings,
                                                 ImpactExcitation.Settings, Dielectronic.Settings, RadiativeAuger.Settings,
                                                 PairAnnihilation1Photon.Settings, ImpactExcitationAutoion.Settings,
                                                 MultiPhotonDeExcitation.Settings, CoulombExcitation.Settings,
                                                 CoulombIonization.Settings}
            ... Provides the settings for the selected process.
    """
    struct  Computation
        name                           ::String
        nuclearModel                   ::Nuclear.Model
        grid                           ::Radial.Grid
        properties                     ::Array{AbstractLevelProperty,1}
        configs                        ::Array{Configuration,1}
        asfSettings                    ::AsfSettings
        initialConfigs                 ::Array{Configuration,1}
        initialAsfSettings             ::AsfSettings
        intermediateConfigs            ::Array{Configuration,1}
        intermediateAsfSettings        ::AsfSettings
        finalConfigs                   ::Array{Configuration,1}
        finalAsfSettings               ::AsfSettings
        alphaSettings                  ::AlphaVariation.Settings
        einsteinSettings               ::Einstein.Settings
        formSettings                   ::FormFactor.Settings
        hfsSettings                    ::Hfs.Settings
        isotopeSettings                ::IsotopeShift.Settings
        plasmaSettings                 ::PlasmaShift.Settings
        polaritySettings               ::MultipolePolarizibility.Settings
        yieldSettings                  ::DecayYield.Settings
        zeemanSettings                 ::LandeZeeman.Settings
        process                        ::AbstractProcess
        processSettings                ::Union{PhotoEmission.Settings, AutoIonization.Settings, PlasmaShift.AugerSettings,
                                               PhotoIonization.Settings, PlasmaShift.PhotoSettings,
                                               PhotoRecombination.Settings, Dielectronic.Settings, ImpactExcitation.Settings,
                                               CoulombExcitation.Settings, CoulombIonization.Settings,
                                               PhotoExcitation.Settings, PhotoExcitationFluores.Settings,
                                               PhotoExcitationAutoion.Settings, RayleighCompton.Settings,
                                               DoubleAutoIonization.Settings, PhotoDoubleIonization.Settings,
                                               MultiPhotonDeExcitation.Settings, PhotoIonizationFluores.Settings,
                                               PhotoIonizationAutoion.Settings, ImpactExcitationAutoion.Settings,
                                               RadiativeAuger.Settings, MultiPhotonIonization.Settings,
                                               MultiPhotonDoubleIon.Settings, InternalConversion.Settings} #= ,
                                               #
                                               PairAnnihilation1Photon.Settings } =#
    end

    """
    `Atomic.Computation()`  ... constructor for an 'empty' instance::Atomic.Computation.
    """
    function Computation()
        Computation("", Nuclear.Model(1.), Radial.Grid(), AbstractLevelProperty[],
                    Configuration[], AsfSettings(),
                    Configuration[], AsfSettings(),
                    Configuration[], AsfSettings(),
                    Configuration[], AsfSettings(),
                    AlphaVariation.Settings(), Einstein.Settings(), FormFactor.Settings(), Hfs.Settings(), IsotopeShift.Settings(),
                    PlasmaShift.Settings(), MultipolePolarizibility.Settings(), DecayYield.Settings(),  LandeZeeman.Settings(),
                    Basics.NoProcess(), PhotoEmission.Settings() )
    end

    """
    `Atomic.Computation(comp::Atomic.Computation;`

        name=..,                nuclearModel=..,            grid=..,                    configs=..,                   asfSettings=..,
        initialConfigs=..,      initialAsfSettings=..,      intermediateConfigs=..,     intermediateAsfSettings=..,
        finalConfigs=..,        finalAsfSettings=..,        alphaSettings=..,           einsteinSettings=..,
        formSettings=..,        hfsSettings=..,             isotopeSettings=..,         plasmaSettings=..,
        polaritySettings=..,    yieldSettings::=..,         zeemanSettings=..,
        process=..,             processSettings=..,         printout::Bool=false)

        ... constructor for modifying the given Atomic.Computation by 'overwriting' the previously selected parameters.
    """
    function Computation(comp::Atomic.Computation;
        name::Union{Nothing,String}=nothing,                                        nuclearModel::Union{Nothing,Nuclear.Model}=nothing,
        grid::Union{Nothing,Radial.Grid}=nothing,                                   properties::Union{Nothing,Array{AbstractLevelProperty,1},Any}=nothing,
        configs::Union{Nothing,Array{Configuration,1}}=nothing,                     asfSettings::Union{Nothing,AsfSettings}=nothing,
        initialConfigs::Union{Nothing,Array{Configuration,1}}=nothing,              initialAsfSettings::Union{Nothing,AsfSettings}=nothing,
        intermediateConfigs::Union{Nothing,Array{Configuration,1}}=nothing,         intermediateAsfSettings::Union{Nothing,AsfSettings}=nothing,
        finalConfigs::Union{Nothing,Array{Configuration,1}}=nothing,                finalAsfSettings::Union{Nothing,AsfSettings}=nothing,
        alphaSettings::Union{Nothing,AlphaVariation.Settings}=nothing,              einsteinSettings::Union{Nothing,Einstein.Settings}=nothing,
        formSettings::Union{Nothing,FormFactor.Settings}=nothing,                   hfsSettings::Union{Nothing,Hfs.Settings}=nothing,
        isotopeSettings::Union{Nothing,IsotopeShift.Settings}=nothing,              plasmaSettings::Union{Nothing,PlasmaShift.Settings}=nothing,
        polaritySettings::Union{Nothing,MultipolePolarizibility.Settings}=nothing,  yieldSettings::Union{Nothing,DecayYield.Settings}=nothing,
        zeemanSettings::Union{Nothing,LandeZeeman.Settings}=nothing,
        process::Union{Nothing,Basics.AbstractProcess}=nothing,                     processSettings::Union{Nothing,Any}=nothing,
        printout::Bool=false)

        if  name                    == nothing  namex                    = comp.name                    else  namex                    = name                     end
        if  nuclearModel            == nothing  nuclearModelx            = comp.nuclearModel            else  nuclearModelx            = nuclearModel             end
        if  grid                    == nothing  gridx                    = comp.grid                    else  gridx                    = grid                     end
        if  properties              == nothing  propertiesx              = comp.properties              else  propertiesx              = properties               end
        if  configs                 == nothing  configsx                 = comp.configs                 else  configsx                 = configs                  end
        if  asfSettings             == nothing  asfSettingsx             = comp.asfSettings             else  asfSettingsx             = asfSettings              end
        if  initialConfigs          == nothing  initialConfigsx          = comp.initialConfigs          else  initialConfigsx          = initialConfigs           end
        if  initialAsfSettings      == nothing  initialAsfSettingsx      = comp.initialAsfSettings      else  initialAsfSettingsx      = initialAsfSettings       end
        if  intermediateConfigs     == nothing  intermediateConfigsx     = comp.intermediateConfigs     else  intermediateConfigsx     = intermediateConfigs      end
        if  intermediateAsfSettings == nothing  intermediateAsfSettingsx = comp.intermediateAsfSettings else  intermediateAsfSettingsx = intermediateAsfSettings  end
        if  finalConfigs            == nothing  finalConfigsx            = comp.finalConfigs            else  finalConfigsx            = finalConfigs             end
        if  finalAsfSettings        == nothing  finalAsfSettingsx        = comp.finalAsfSettings        else  finalAsfSettingsx        = finalAsfSettings         end
        if  alphaSettings           == nothing  alphaSettingsx           = comp.alphaSettings           else  alphaSettingsx           = alphaSettings            end
        if  einsteinSettings        == nothing  einsteinSettingsx        = comp.einsteinSettings        else  einsteinSettingsx        = einsteinSettings         end
        if  formSettings            == nothing  formSettingsx            = comp.formSettings            else  formSettingsx            = formSettings             end
        if  hfsSettings             == nothing  hfsSettingsx             = comp.hfsSettings             else  hfsSettingsx             = hfsSettings              end
        if  isotopeSettings         == nothing  isotopeSettingsx         = comp.isotopeSettings         else  isotopeSettingsx         = isotopeSettings          end
        if  plasmaSettings          == nothing  plasmaSettingsx          = comp.plasmaSettings          else  plasmaSettingsx          = plasmaSettings           end
        if  polaritySettings        == nothing  polaritySettingsx        = comp.polaritySettings        else  polaritySettingsx        = polaritySettings         end
        if  yieldSettings           == nothing  yieldSettingsx           = comp.yieldSettings           else  yieldSettingsx           = yieldSettings            end
        if  zeemanSettings          == nothing  zeemanSettingsx          = comp.zeemanSettings          else  zeemanSettingsx          = zeemanSettings           end
        if  process                 == nothing  processx                 = comp.process                 else  processx                 = process                  end
        if  processSettings         == nothing  prsx                     = comp.processSettings         else  prsx                     = processSettings          end

        if      processx == Basics.NoProcess                                                              prsx = PhotoEmission.Settings()
        elseif  processx == Basics.Auger            && typeof(prsx) != AutoIonization.Settings            prsx = AutoIonization.Settings()
        elseif  processx == Basics.AugerInPlasma    && typeof(prsx) != PlasmaShift.AugerSettings          prsx = PlasmaShift.AugerSettings()
        elseif  processx == Basics.Radiative        && typeof(prsx) != PhotoEmission.Settings             prsx = PhotoEmission.Settings()
        elseif  processx == Basics.PhotoExc         && typeof(prsx) != PhotoExcitation.Settings           prsx = PhotoExcitation.Settings()
        elseif  process == Basics.Photo             && typeof(prsx) != PhotoIonization.Settings           prsx = PhotoIonization.Settings()
        elseif  process == Basics.PhotoInPlasma     && typeof(prsx) != PlasmaShift.PhotoSettings          prsx = PlasmaShift.PhotoSettings()
        elseif  process == Basics.Rec               && typeof(prsx) != PhotoRecombination.Settings        prsx = PhotoRecombination.Settings()
        elseif  process == Basics.Dierec            && typeof(prsx) != Dielectronic.Settings              prsx = Dielectronic.Settings()
        elseif  process == Basics.PhotoExcFluor     && typeof(prsx) != PhotoExcitationFluores.Settings    prsx = PhotoExcitationFluores.Settings()
        elseif  process == Basics.PhotoExcAuto      && typeof(prsx) != PhotoExcitationAutoion.Settings    prsx = PhotoExcitationAutoion.Settings()
        elseif  process == Basics.PhotoIonFluor     && typeof(prsx) != PhotoIonizationFluores.Settings    prsx = PhotoIonizationFluores.Settings()
        elseif  process == Basics.PhotoIonAuto      && typeof(prsx) != PhotoIonizationAutoion.Settings    prsx = PhotoIonizationAutoion.Settings()
        elseif  process == Basics.MultiPhotonDE     && typeof(prsx) != MultiPhotonDeExcitation.Settings   prsx = MultiPhotonDeExcitation.Settings()
        elseif  process == Basics.Compton           && typeof(prsx) != RayleighCompton.Settings           prsx = RayleighCompton.Settings()
        elseif  process == Basics.Eimex             && typeof(prsx) != ImpactExcitation.Settings          prsx = ImpactExcitation.Settings()
        elseif  process == Basics.ImpactExcAuto     && typeof(prsx) != ImpactExcitationAutoion.Settings   prsx = ImpactExcitationAutoion.Settings()
        elseif  process == Basics.RAuger            && typeof(prsx) != RadiativeAuger.Settings            prsx = RadiativeAuger.Settings()
        elseif  process == Basics.MultiPI           && typeof(prsx) != MultiPhotonIonization.Settings     prsx = MultiPhotonIonization.Settings()
        elseif  process == Basics.MultiPDI          && typeof(prsx) != MultiPhotonDoubleIon.Settings      prsx = MultiPhotonDoubleIon.Settings()
        elseif  process == Basics.InternalConv      && typeof(prsx) != InternalConversion.Settings        prsx = InternalConversion.Settings()
        elseif  process == Basics.Coulex            && typeof(prsx) != CoulombExcitation.Settings         prsx = CoulombExcitation.Settings()
        elseif  process == Basics.Coulion           && typeof(prsx) != CoulombIonization.Settings         prsx = CoulombIonization.Settings()
        end

        cp = Computation(namex, nuclearModelx, gridx, propertiesx, configsx, asfSettingsx, initialConfigsx, initialAsfSettingsx,
                         intermediateConfigsx,  intermediateAsfSettingsx, finalConfigsx, finalAsfSettingsx,
                         alphaSettingsx, einsteinSettingsx , formSettingsx, hfsSettingsx, isotopeSettingsx, plasmaSettingsx,
                         polaritySettingsx, yieldSettingsx, zeemanSettingsx,  processx, prsx)

        if printout  Base.show(cp)      end
        return( cp )
    end

    """
    `Atomic.Computation( ... example for SCF computations)`

            grid     = Radial.Grid(true)
            nuclearM = Nuclear.Model(18., "Fermi")
            settings = AsfSettings(AsfSettings(), selectLevelsCI = true, selectedLevelsCI = [1,2, 4,5, 7,8], jjLS = LSjjSettings(false) )
            configs  = [Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6")]
            Atomic.Computation(Atomic.Computation(), name="Example", grid=grid, nuclearModel=nuclearM, configs=configs, asfSettings=settings )

    `Atomic.Computation( ... example for the computation of atomic properties)`

            grid        = Radial.Grid(true)
            nuclearM    = Nuclear.Model(26., "Fermi", 58., 3.81, AngularJ64(5//2), 1.0, 1.0)
            hfsSettings = Hfs.Settings(true, true, false, false, false, false, false, Int64[] )
            configs     = [Configuration("[Ne] 3s"), Configuration("[Ne] 3p"), Configuration("[Ne] 3d")]
            Atomic.Computation(Atomic.Computation(), name="Example", grid=grid, nuclearModel=nuclearM, configs=configs, properties=[HFS()],
                               hfsSettings=hfsSettings )

    `Atomic.Computation( ... example for the computation of one atomic process)`

            grid           = Radial.Grid(true)
            initialConfigs = [Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")]
            finalConfigs   = [Configuration("[Ne] 3s^2 3p^5")]
            photoSettings  = PhotoEmission.Settings(PhotoEmission.Settings(), multipoles=[E1, M1], gauges=[UseCoulomb], printBefore=true)
            Atomic.Computation(Atomic.Computation(), name="Example", grid=grid, nuclearModel=nuclearM;
                              initialConfigs=initialConfigs, finalConfigs=finalConfigs,
                              process = Radiative(), processSettings=photoSettings );

        ... These simple examples can be further improved by overwriting the corresponding parameters.
    """
    function Computation(wa::Bool)
        Atomic.Computation()
    end

    # `Base.string(comp::Atomic.Computation)`  ... provides a String notation for the variable comp::Atomic.Computation.
    function Base.string(comp::Atomic.Computation)
        sa = "Atomic computation:    $(comp.name) for Z = $(comp.nuclearModel.Z), "
        if  comp.process != NoProcess   sa = sa * "for the process $(comp.process) and with the \ninitial configurations:    "
            for  config  in  comp.initialConfigs   sa = sa * string(config) * ",  "         end
            if  length(comp.intermediateConfigs) > 0
                sa = sa * "\nintermediate configurations:"
                for  config  in  comp.finalConfigs     sa = sa * string(config) * ",  "     end
            end
            sa = sa * "\nfinal configurations:      "
            for  config  in  comp.finalConfigs     sa = sa * string(config) * ",  "         end
        else                                sa = sa * "for the properties $(comp.properties) and with the \nconfigurations:        "
            for  config  in  comp.configs   sa = sa * string(config) * ",  "                end
        end
        return( sa )
    end

    # `Base.show(io::IO, comp::Atomic.Computation)`  ... prepares a printout of comp::Atomic.Computation.
    function Base.show(io::IO, comp::Atomic.Computation)
        sa = Base.string(comp);             print(io, sa, "\n")
        println(io, "nuclearModel:          $(comp.nuclearModel)  ")
        println(io, "grid:                  $(comp.grid)  ")
        #
        # For the computation of some given atomic process
        if  comp.process != NoProcess
        println(io, "processSettings:              \n$(comp.processSettings)  ")
        if  comp.initialAsfSettings != AsfSettings()
        println(io, "initialAsfSettings:           \n$(comp.initialAsfSettings)  ")         end
        if  comp.intermediateAsfSettings != AsfSettings()
        println(io, "intermediateAsfSettings:      \n$(comp.intermediateAsfSettings)  ")    end
        if  comp.finalAsfSettings != AsfSettings()
        println(io, "finalAsfSettings:             \n$(comp.finalAsfSettings)  ")           end
        #
        else
        # For the computation of no or several atomic properties
        if  comp.asfSettings != AsfSettings()
        println(io, "asfSettings:                  \n$(comp.asfSettings)  ")                end
        if  AlphaX in comp.properties     &&  comp.alphaSettings    != AlphaVariation.Settings()
        println(io, "alphaSettings:                \n$(comp.alphaSettings)  ")              end
        if  EinsteinX in comp.properties  &&  comp.einsteinSettings != Einstein.Settings()
        println(io, "einsteinSettings:             \n$(comp.einsteinSettings)  ")           end
        if  FormF in comp.properties      &&  comp.formSettings     != FormFactor.Settings()
        println(io, "formSettings:                 \n$(comp.formSettings)  ")               end
        if  HFS in comp.properties        &&  comp.hfsSettings      != Hfs.Settings()
        println(io, "hfsSettings:                  \n$(comp.hfsSettings)  ")                end
        if  Isotope in comp.properties    &&  comp.isotopeSettings  != IsotopeShift.Settings()
        println(io, "isotopeSettings:              \n$(comp.isotopeSettings)  ")            end
        if  Plasma in comp.properties     &&  comp.plasmaSettings   != PlasmaShift.Settings()
        println(io, "plasmaSettings:               \n$(comp.plasmaSettings)  ")             end
        if  Polarity in comp.properties   &&  comp.polaritySettings != MultipolePolarizibility.Settings()
        println(io, "polaritySettings:             \n$(comp.polaritySettings)  ")           end
        if  Yields in comp.properties     &&  comp.yieldSettings    != DecayYield.Settings()
        println(io, "yieldSettings:                \n$(comp.yieldSettings)  ")              end
        if  Zeeman in comp.properties     &&  comp.zeemanSettings   != LandeZeeman.Settings()
        println(io, "zeemanSettings:               \n$(comp.zeemanSettings)  ")             end
        end
    end

end # module


text:
partial(x::TangentBundle, i) = x.partials[i]
partial(x::TaylorBundle{1}, i) = x.coeffs[i]
partial(x::UniformBundle, i) = x.partial
partial(x::CompositeBundle{N, B}, i) where {N, B} = Tangent{B}(map(x->partial(x, i), x.tup)...)
partial(x::ZeroTangent, i) = ZeroTangent()
primal(x::AbstractTangentBundle) = x.primal
primal(z::ZeroTangent) = ZeroTangent()

first_partial(x::TangentBundle{1}) = getfield(getfield(x, :partials), 1)
first_partial(x::TaylorBundle{1}) = getfield(getfield(x, :coeffs), 1)
first_partial(x::UniformBundle) = getfield(x, :partial)
first_partial(x::CompositeBundle) = map(first_partial, getfield(x, :tup))

# TODO: Which version do we want in ChainRules?
function my_frule(args::ATB{1}...)
    frule(map(first_partial, args), map(primal, args)...)
end

# Fast path for some hot cases
my_frule(::ZeroBundle{1, typeof(frule)}, args::ATB{1}...) = nothing
my_frule(::ZeroBundle{1, typeof(my_frule)}, args::ATB{1}...) = nothing

(::∂☆{N})(::ZeroBundle{N, typeof(my_frule)}, ::ZeroBundle{N, ZeroBundle{1, typeof(frule)}}, args::ATB{N}...) where {N} = ZeroBundle{N}(nothing)
(::∂☆{N})(::ZeroBundle{N, typeof(my_frule)}, ::ZeroBundle{N, ZeroBundle{1, typeof(my_frule)}}, args::ATB{N}...) where {N} = ZeroBundle{N}(nothing)

shuffle_down(b::UniformBundle{N, B, U}) where {N, B, U} =
    UniformBundle{minus1(N), <:Any, U}(UniformBundle{1, B, U}(b.primal, b.partial), b.partial)

function shuffle_down(b::TangentBundle{N, B}) where {N, B}
    # N.B: This depends on the special properties of the canonical tangent index order
    TangentBundle{N-1}(
        TangentBundle{1}(b.primal, (partial(b, 1),)),
        ntuple(2^(N-1)-1) do i
            TangentBundle{1}(partial(b, 2*i), (partial(b, 2*i+1),))
        end)
end

function shuffle_down(b::TaylorBundle{N, B}) where {N, B}
    TaylorBundle{N-1}(
        TangentBundle{1}(b.primal, (b.coeffs[1],)),
        ntuple(N-1) do i
            TangentBundle{1}(b.coeffs[i], (b.coeffs[i+1],))
        end)
end

function shuffle_down(b::CompositeBundle{N, B}) where {N, B}
    z = CompositeBundle{N-1, CompositeBundle{1, B}}(
        (CompositeBundle{N-1, Tuple}(
            map(shuffle_down, b.tup)
        ),)
    )
    z
end

function shuffle_up(r::CompositeBundle{1})
    z₀ = primal(r.tup[1])
    z₁ = partial(r.tup[1], 1)
    z₂ = primal(r.tup[2])
    z₁₂ = partial(r.tup[2], 1)
    if z₁ == z₂
        return TaylorBundle{2}(z₀, (z₁, z₁₂))
    else
        return TangentBundle{2}(z₀, (z₁, z₂, z₁₂))
    end
end

function taylor_compatible(a::ATB{N}, b::ATB{N}) where {N}
    primal(b) === a[TaylorTangentIndex(1)] || return false
    return all(1:(N-1)) do i
        b[TaylorTangentIndex(i)] === a[TaylorTangentIndex(i+1)]
    end
end

# Check whether the tangent bundle element is taylor-like
isswifty(::TaylorBundle) = true
isswifty(::UniformBundle) = true
isswifty(b::CompositeBundle) = all(isswifty, b.tup)
isswifty(::Any) = false

function shuffle_up(r::CompositeBundle{N}) where {N}
    a, b = r.tup
    if isswifty(a) && isswifty(b) && taylor_compatible(a, b)
        return TaylorBundle{N+1}(primal(a),
            ntuple(i->i == N+1 ?
                b[TaylorTangentIndex(i-1)] : a[TaylorTangentIndex(i)],
            N+1))
    else
        return TangentBundle{N+1}(r.tup[1].primal,
            (r.tup[1].partials..., primal(b),
            ntuple(i->partial(b,i), 2^(N+1)-1)...))
    end
end

function shuffle_up(r::UniformBundle{N, B, U}) where {N, B, U}
    (a, b) = primal(r)
    if r.partial === b
        u = b
    elseif b == NoTangent() && U === ZeroTangent
        u = b
    else
        error()
    end
    UniformBundle{N+1}(a, u)
end
@ChainRulesCore.non_differentiable shuffle_up(r::UniformBundle)

struct ∂☆internal{N}; end
struct ∂☆shuffle{N}; end

shuffle_base(r) = TangentBundle{1}(r[1], (r[2],))

function (::∂☆internal{1})(args::AbstractTangentBundle{1}...)
    r = my_frule(args...)
    if r === nothing
        return ∂☆recurse{1}()(args...)
    else
        return shuffle_base(r)
    end
end

function (::∂☆shuffle{N})(args::AbstractTangentBundle{N}...) where {N}
    ∂☆p = ∂☆{minus1(N)}()
    ∂☆p(ZeroBundle{minus1(N)}(my_frule), map(shuffle_down, args)...)
end

function (::∂☆internal{N})(args::AbstractTangentBundle{N}...) where {N}
    r = ∂☆shuffle{N}()(args...)
    if primal(r) === nothing
        return ∂☆recurse{N}()(args...)
    else
        return shuffle_up(r)
    end
end
(::∂☆{N})(args::AbstractTangentBundle{N}...) where {N} = ∂☆internal{N}()(args...)

# Special case rules for performance
@Base.aggressive_constprop function (::∂☆{N})(f::ATB{N, typeof(getfield)}, x::TangentBundle{N}, s::AbstractTangentBundle{N}) where {N}
    s = primal(s)
    TangentBundle{N}(getfield(primal(x), s),
        map(x->lifted_getfield(x, s), x.partials))
end

@Base.aggressive_constprop function (::∂☆{N})(f::ATB{N, typeof(getfield)}, x::TaylorBundle{N}, s::AbstractTangentBundle{N}) where {N}
    s = primal(s)
    TaylorBundle{N}(getfield(primal(x), s),
        map(y->lifted_getfield(y, s), x.coeffs))
end

@Base.aggressive_constprop function (::∂☆{N})(::ATB{N, typeof(getfield)}, x::CompositeBundle{N}, s::AbstractTangentBundle{N, Int}) where {N}
    x.tup[primal(s)]
end

@Base.aggressive_constprop function (::∂☆{N})(::ATB{N, typeof(getfield)}, x::CompositeBundle{N, B}, s::AbstractTangentBundle{N, Symbol}) where {N, B}
    x.tup[Base.fieldindex(B, primal(s))]
end

@Base.aggressive_constprop function (::∂☆{N})(f::ATB{N, typeof(getfield)}, x::ATB{N}, s::ATB{N}, inbounds::ATB{N}) where {N}
    s = primal(s)
    TangentBundle{N}(getfield(primal(x), s, primal(inbounds)),
        map(x->lifted_getfield(x, s), x.partials))
end

@Base.aggressive_constprop function (::∂☆{N})(f::ATB{N, typeof(getfield)}, x::UniformBundle{N, <:Any, U}, s::AbstractTangentBundle{N}) where {N, U}
    UniformBundle{N,<:Any,U}(getfield(primal(x), primal(s)), x.partial)
end

@Base.aggressive_constprop function (::∂☆{N})(f::ATB{N, typeof(getfield)}, x::UniformBundle{N, <:Any, U}, s::AbstractTangentBundle{N}, inbounds::AbstractTangentBundle{N}) where {N, U}
    UniformBundle{N,<:Any,U}(getfield(primal(x), primal(s), primal(inbounds)), x.partial)
end

function (::∂☆{N})(f::ATB{N, typeof(tuple)}, args::AbstractTangentBundle{N}...) where {N}
    ∂vararg{N}()(args...)
end

struct FwdMap{N, T<:AbstractTangentBundle{N}}
    f::T
end
(f::FwdMap{N})(args::AbstractTangentBundle{N}...) where {N} = ∂☆{N}()(f.f, args...)

function (::∂☆{N})(::ZeroBundle{N, typeof(map)}, f::ATB{N}, tup::CompositeBundle{N, <:Tuple}) where {N}
    ∂vararg{N}()(map(FwdMap(f), tup.tup)...)
end

function (::∂☆{N})(::ZeroBundle{N, typeof(map)}, f::ATB{N}, args::ATB{N, <:AbstractArray}...) where {N}
    # TODO: This could do an inplace map! to avoid the extra rebundling
    rebundle(map(FwdMap(f), map(unbundle, args)...))
end

function (::∂☆{N})(::ZeroBundle{N, typeof(map)}, f::ATB{N}, args::ATB{N}...) where {N}
    ∂☆recurse{N}()(ZeroBundle{N, typeof(map)}(map), f, args...)
end

function (::∂☆{N})(f::ZeroBundle{N, typeof(ifelse)}, arg::ATB{N, Bool}, args::ATB{N}...) where {N}
    ifelse(arg.primal, args...)
end

struct FwdIterate{N, T<:AbstractTangentBundle{N}}
    f::T
end
function (f::FwdIterate)(arg::ATB{N}) where {N}
    r = ∂☆{N}()(f.f, arg)
    primal(r) === nothing && return nothing
    (∂☆{N}()(ZeroBundle{N}(getindex), r, ZeroBundle{N}(1)),
     primal(∂☆{N}()(ZeroBundle{N}(getindex), r, ZeroBundle{N}(2))))
end
function (f::FwdIterate)(arg::ATB{N}, st) where {N}
    r = ∂☆{N}()(f.f, arg, ZeroBundle{N}(st))
    primal(r) === nothing && return nothing
    (∂☆{N}()(ZeroBundle{N}(getindex), r, ZeroBundle{N}(1)),
     primal(∂☆{N}()(ZeroBundle{N}(getindex), r, ZeroBundle{N}(2))))
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(Core._apply_iterate)}, iterate::ATB{N}, f::ATB{N}, args::ATB{N}...) where {N}
    Core._apply_iterate(FwdIterate(iterate), this, (f,), args...)
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(iterate)}, t::CompositeBundle{N, <:Tuple}) where {N}
    r = iterate(t.tup)
    r === nothing && return ZeroBundle{N}(nothing)
    ∂vararg{N}()(r[1], ZeroBundle{N}(r[2]))
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(iterate)}, t::CompositeBundle{N, <:Tuple}, a::ATB{N}, args::ATB{N}...) where {N}
    r = iterate(t.tup, primal(a), map(primal, args)...)
    r === nothing && return ZeroBundle{N}(nothing)
    ∂vararg{N}()(r[1], ZeroBundle{N}(r[2]))
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(Base.indexed_iterate)}, t::CompositeBundle{N, <:Tuple}, i::ATB{N}) where {N}
    r = Base.indexed_iterate(t.tup, primal(i))
    ∂vararg{N}()(r[1], ZeroBundle{N}(r[2]))
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(Base.indexed_iterate)}, t::CompositeBundle{N, <:Tuple}, i::ATB{N}, st1::ATB{N}, st::ATB{N}...) where {N}
    r = Base.indexed_iterate(t.tup, primal(i), primal(st1), map(primal, st)...)
    ∂vararg{N}()(r[1], ZeroBundle{N}(r[2]))
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(Base.indexed_iterate)}, t::TangentBundle{N, <:Tuple}, i::ATB{N}, st::ATB{N}...) where {N}
    ∂vararg{N}()(this(ZeroBundle{N}(getfield), t, i), ZeroBundle{N}(primal(i) + 1))
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(getindex)}, t::CompositeBundle{N, <:Tuple}, i::ZeroBundle) where {N}
    t.tup[primal(i)]
end

function (this::∂☆{N})(::ZeroBundle{N, typeof(typeof)}, x::ATB{N}) where {N}
    DNEBundle{N}(typeof(primal(x)))
end

function (this::∂☆{N})(f::ZeroBundle{N, Core.IntrinsicFunction}, args::ATB{N}...) where {N}
    ff = primal(f)
    if ff === Base.not_int
        DNEBundle{N}(ff(map(primal, args)...))
    else
        error("Missing rule for intrinsic function $ff")
    end
end


text:
module MieScattering

    using SpecialFunctions

    export scattering, absorption, extinction, backscattering

    export Scatterer, mx
    include("interface.jl")

    export Mie, mie_scattering, scattering_function
    include("mie.jl")

    export Rayleigh, rayleigh_scattering
    include("rayleigh.jl")

end

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# The methods below are specialized on the mortar type
# and called from the basic `create_cache` method at the top.
function create_cache(mesh::P4estMesh{3}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype)
  # TODO: Taal compare performance of different types
  fstar_threaded = [Array{uEltype, 4}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2), 4)
                    for _ in 1:Threads.nthreads()]

  fstar_tmp_threaded = [Array{uEltype, 3}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                         for _ in 1:Threads.nthreads()]
  u_threaded          = [Array{uEltype, 3}(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2))
                         for _ in 1:Threads.nthreads()]

  (; fstar_threaded, fstar_tmp_threaded, u_threaded)
end

#     index_to_start_step_3d(index::Symbol, index_range)
#
# Given a symbolic `index` and an `indexrange` (usually `eachnode(dg)`),
# return `index_start, index_step_i, index_step_j`, i.e., a tuple containing
# - `index_start`,   an index value to begin a loop
# - `index_step_i`,  an index step to update during an `i` loop
# - `index_step_j`,  an index step to update during a `j` loop
# The resulting indices translate surface indices to volume indices.
#
# !!! warning
#     This assumes that loops using the return values are written as
#
#     i_volume_start, i_volume_step_i, i_volume_step_j = index_to_start_step_3d(symbolic_index_i, index_range)
#     j_volume_start, j_volume_step_i, j_volume_step_j = index_to_start_step_3d(symbolic_index_j, index_range)
#     k_volume_start, k_volume_step_i, k_volume_step_j = index_to_start_step_3d(symbolic_index_k, index_range)
#
#     i_volume, j_volume, k_volume = i_volume_start, j_volume_start, k_volume_start
#     for j_surface in index_range
#       for i_surface in index_range
#         # do stuff with `(i_surface, j_surface)` and `(i_volume, j_volume, k_volume)`
#
#         i_volume += i_volume_step_i
#         j_volume += j_volume_step_i
#         k_volume += k_volume_step_i
#       end
#       i_volume += i_volume_step_j
#       j_volume += j_volume_step_j
#       k_volume += k_volume_step_j
#     end
@inline function index_to_start_step_3d(index::Symbol, index_range)
  index_begin = first(index_range)
  index_end   = last(index_range)

  if index === :begin
    return index_begin, 0, 0
  elseif index === :end
    return index_end, 0, 0
  elseif index === :i_forward
    return index_begin, 1, index_begin - index_end - 1
  elseif index === :i_backward
    return index_end, -1, index_end + 1 - index_begin
  elseif index === :j_forward
    return index_begin, 0, 1
  else # if index === :j_backward
    return index_end, 0, -1
  end
end

# Extract the two varying indices from a symbolic index tuple.
# For example, `surface_indices((:i_forward, :end, :j_forward)) == (:i_forward, :j_forward)`.
@inline function surface_indices(indices::NTuple{3, Symbol})
  i1, i2, i3 = indices
  index = i1
  (index === :begin || index === :end) && return (i2, i3)

  index = i2
  (index === :begin || index === :end) && return (i1, i3)

  # i3 in (:begin, :end)
  return (i1, i2)
end

function prolong2interfaces!(cache, u,
                             mesh::P4estMesh{3},
                             equations, surface_integral, dg::DG)
  @unpack interfaces = cache
  index_range = eachnode(dg)

  @threaded for interface in eachinterface(dg, cache)
    # Copy solution data from the primary element using "delayed indexing" with
    # a start value and two step sizes to get the correct face and orientation.
    # Note that in the current implementation, the interface will be
    # "aligned at the primary element", i.e., the indices of the primary side
    # will always run forwards.
    primary_element = interfaces.neighbor_ids[1, interface]
    primary_indices = interfaces.node_indices[1, interface]

    i_primary_start, i_primary_step_i, i_primary_step_j = index_to_start_step_3d(primary_indices[1], index_range)
    j_primary_start, j_primary_step_i, j_primary_step_j = index_to_start_step_3d(primary_indices[2], index_range)
    k_primary_start, k_primary_step_i, k_primary_step_j = index_to_start_step_3d(primary_indices[3], index_range)

    i_primary = i_primary_start
    j_primary = j_primary_start
    k_primary = k_primary_start
    for j in eachnode(dg)
      for i in eachnode(dg)
        for v in eachvariable(equations)
          interfaces.u[1, v, i, j, interface] = u[v, i_primary, j_primary, k_primary, primary_element]
        end
        i_primary += i_primary_step_i
        j_primary += j_primary_step_i
        k_primary += k_primary_step_i
      end
      i_primary += i_primary_step_j
      j_primary += j_primary_step_j
      k_primary += k_primary_step_j
    end

    # Copy solution data from the secondary element using "delayed indexing" with
    # a start value and two step sizes to get the correct face and orientation.
    secondary_element = interfaces.neighbor_ids[2, interface]
    secondary_indices = interfaces.node_indices[2, interface]

    i_secondary_start, i_secondary_step_i, i_secondary_step_j = index_to_start_step_3d(secondary_indices[1], index_range)
    j_secondary_start, j_secondary_step_i, j_secondary_step_j = index_to_start_step_3d(secondary_indices[2], index_range)
    k_secondary_start, k_secondary_step_i, k_secondary_step_j = index_to_start_step_3d(secondary_indices[3], index_range)

    i_secondary = i_secondary_start
    j_secondary = j_secondary_start
    k_secondary = k_secondary_start
    for j in eachnode(dg)
      for i in eachnode(dg)
        for v in eachvariable(equations)
          interfaces.u[2, v, i, j, interface] = u[v, i_secondary, j_secondary, k_secondary, secondary_element]
        end
        i_secondary += i_secondary_step_i
        j_secondary += j_secondary_step_i
        k_secondary += k_secondary_step_i
      end
      i_secondary += i_secondary_step_j
      j_secondary += j_secondary_step_j
      k_secondary += k_secondary_step_j
    end
  end

  return nothing
end

function calc_interface_flux!(surface_flux_values,
                              mesh::P4estMesh{3},
                              nonconservative_terms,
                              equations, surface_integral, dg::DG, cache)
  @unpack neighbor_ids, node_indices = cache.interfaces
  @unpack contravariant_vectors = cache.elements
  index_range = eachnode(dg)

  @threaded for interface in eachinterface(dg, cache)
    # Get element and side information on the primary element
    primary_element  = neighbor_ids[1, interface]
    primary_indices  = node_indices[1, interface]
    primary_direction = indices2direction(primary_indices)

    i_primary_start, i_primary_step_i, i_primary_step_j = index_to_start_step_3d(primary_indices[1], index_range)
    j_primary_start, j_primary_step_i, j_primary_step_j = index_to_start_step_3d(primary_indices[2], index_range)
    k_primary_start, k_primary_step_i, k_primary_step_j = index_to_start_step_3d(primary_indices[3], index_range)

    i_primary = i_primary_start
    j_primary = j_primary_start
    k_primary = k_primary_start

    # Get element and side information on the secondary element
    secondary_element = neighbor_ids[2, interface]
    secondary_indices = node_indices[2, interface]
    secondary_direction = indices2direction(secondary_indices)
    secondary_surface_indices = surface_indices(secondary_indices)

    # Get the surface indexing on the secondary element.
    # Note that the indices of the primary side will always run forward but
    # the secondary indices might need to run backwards for flipped sides.
    i_secondary_start, i_secondary_step_i, i_secondary_step_j = index_to_start_step_3d(secondary_surface_indices[1], index_range)
    j_secondary_start, j_secondary_step_i, j_secondary_step_j = index_to_start_step_3d(secondary_surface_indices[2], index_range)
    i_secondary = i_secondary_start
    j_secondary = j_secondary_start

    for j in eachnode(dg)
      for i in eachnode(dg)
        # Get the normal direction from the primary element.
        # Note, contravariant vectors at interfaces in negative coordinate direction
        # are pointing inwards. This is handled by `get_normal_direction`.
        normal_direction = get_normal_direction(primary_direction, contravariant_vectors,
                                                i_primary, j_primary, k_primary,
                                                primary_element)

        calc_interface_flux!(surface_flux_values, mesh, nonconservative_terms, equations,
                             surface_integral, dg, cache,
                             interface, normal_direction,
                             i, j, primary_direction, primary_element,
                             i_secondary, j_secondary, secondary_direction, secondary_element)

        # Increment the primary element indices
        i_primary += i_primary_step_i
        j_primary += j_primary_step_i
        k_primary += k_primary_step_i
        # Increment the secondary element surface indices
        i_secondary += i_secondary_step_i
        j_secondary += j_secondary_step_i
      end
      # Increment the primary element indices
      i_primary += i_primary_step_j
      j_primary += j_primary_step_j
      k_primary += k_primary_step_j
      # Increment the secondary element surface indices
      i_secondary += i_secondary_step_j
      j_secondary += j_secondary_step_j
    end
  end

  return nothing
end

# Inlined function for interface flux computation for conservative flux terms
@inline function calc_interface_flux!(surface_flux_values,
                                      mesh::P4estMesh{3},
                                      nonconservative_terms::Val{false}, equations,
                                      surface_integral, dg::DG, cache,
                                      interface_index, normal_direction,
                                      primary_i_node_index, primary_j_node_index,
                                      primary_direction_index, primary_element_index,
                                      secondary_i_node_index, secondary_j_node_index,
                                      secondary_direction_index, secondary_element_index)
  @unpack u = cache.interfaces
  @unpack surface_flux = surface_integral

  u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_i_node_index, primary_j_node_index, interface_index)

  flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)

  for v in eachvariable(equations)
    surface_flux_values[v, primary_i_node_index, primary_j_node_index,
                        primary_direction_index, primary_element_index] = flux_[v]
    surface_flux_values[v, secondary_i_node_index, secondary_j_node_index,
                        secondary_direction_index, secondary_element_index] = -flux_[v]
  end
end

# Inlined function for interface flux computation for flux + nonconservative terms
@inline function calc_interface_flux!(surface_flux_values,
                                      mesh::P4estMesh{3},
                                      nonconservative_terms::Val{true}, equations,
                                      surface_integral, dg::DG, cache,
                                      interface_index, normal_direction,
                                      primary_i_node_index, primary_j_node_index,
                                      primary_direction_index, primary_element_index,
                                      secondary_i_node_index, secondary_j_node_index,
                                      secondary_direction_index, secondary_element_index)
  @unpack u = cache.interfaces
  surface_flux, nonconservative_flux = surface_integral.surface_flux

  u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_i_node_index, primary_j_node_index, interface_index)

  flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)

  # Compute both nonconservative fluxes
  # In general, nonconservative fluxes can depend on both the contravariant
  # vectors (normal direction) at the current node and the averaged ones.
  # However, both are the same at watertight interfaces, so we pass the
  # `normal_direction` twice.
  noncons_primary   = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations)
  noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, normal_direction, equations)

  # Store the flux with nonconservative terms on the primary and secondary elements
  for v in eachvariable(equations)
    # Note the factor 0.5 necessary for the nonconservative fluxes based on
    # the interpretation of global SBP operators coupled discontinuously via
    # central fluxes/SATs
    surface_flux_values[v, primary_i_node_index, primary_j_node_index,
                        primary_direction_index, primary_element_index] = flux_[v] + 0.5 * noncons_primary[v]
    surface_flux_values[v, secondary_i_node_index, secondary_j_node_index,
                        secondary_direction_index, secondary_element_index] = -(flux_[v] + 0.5 * noncons_secondary[v])
  end
end

function prolong2boundaries!(cache, u,
                             mesh::P4estMesh{3},
                             equations, surface_integral, dg::DG)
  @unpack boundaries = cache
  index_range = eachnode(dg)

  @threaded for boundary in eachboundary(dg, cache)
    # Copy solution data from the element using "delayed indexing" with
    # a start value and two step sizes to get the correct face and orientation.
    element      = boundaries.neighbor_ids[boundary]
    node_indices = boundaries.node_indices[boundary]

    i_node_start, i_node_step_i, i_node_step_j = index_to_start_step_3d(node_indices[1], index_range)
    j_node_start, j_node_step_i, j_node_step_j = index_to_start_step_3d(node_indices[2], index_range)
    k_node_start, k_node_step_i, k_node_step_j = index_to_start_step_3d(node_indices[3], index_range)

    i_node = i_node_start
    j_node = j_node_start
    k_node = k_node_start
    for j in eachnode(dg)
      for i in eachnode(dg)
        for v in eachvariable(equations)
          boundaries.u[v, i, j, boundary] = u[v, i_node, j_node, k_node, element]
        end
        i_node += i_node_step_i
        j_node += j_node_step_i
        k_node += k_node_step_i
      end
      i_node += i_node_step_j
      j_node += j_node_step_j
      k_node += k_node_step_j
    end
  end

  return nothing
end

function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing,
                             mesh::P4estMesh{3},
                             equations, surface_integral, dg::DG)
  @unpack boundaries = cache
  @unpack surface_flux_values, node_coordinates, contravariant_vectors = cache.elements
  @unpack surface_flux = surface_integral
  index_range = eachnode(dg)

  @threaded for local_index in eachindex(boundary_indexing)
    # Use the local index to get the global boundary index from the
    # pre-sorted list
    boundary = boundary_indexing[local_index]

    # Get information on the adjacent element, compute the surface fluxes,
    # and store them
    element       = boundaries.neighbor_ids[boundary]
    node_indices  = boundaries.node_indices[boundary]
    direction     = indices2direction(node_indices)

    i_node_start, i_node_step_i, i_node_step_j = index_to_start_step_3d(node_indices[1], index_range)
    j_node_start, j_node_step_i, j_node_step_j = index_to_start_step_3d(node_indices[2], index_range)
    k_node_start, k_node_step_i, k_node_step_j = index_to_start_step_3d(node_indices[3], index_range)

    i_node = i_node_start
    j_node = j_node_start
    k_node = k_node_start
    for j in eachnode(dg)
      for i in eachnode(dg)
        # Extract solution data from boundary container
        u_inner = get_node_vars(boundaries.u, equations, dg, i, j, boundary)

        # Outward-pointing normal direction (not normalized)
        normal_direction = get_normal_direction(direction, contravariant_vectors,
                                                i_node, j_node, k_node, element)

        # Coordinates at boundary node
        x = get_node_coords(node_coordinates, equations, dg,
                            i_node, j_node, k_node, element)

        flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations)

        # Copy flux to element storage in the correct orientation
        for v in eachvariable(equations)
          surface_flux_values[v, i, j, direction, element] = flux_[v]
        end

        i_node += i_node_step_i
        j_node += j_node_step_i
        k_node += k_node_step_i
      end
      i_node += i_node_step_j
      j_node += j_node_step_j
      k_node += k_node_step_j
    end
  end
end

function prolong2mortars!(cache, u,
                          mesh::P4estMesh{3}, equations,
                          mortar_l2::LobattoLegendreMortarL2,
                          surface_integral, dg::DGSEM)
  @unpack fstar_tmp_threaded = cache
  @unpack neighbor_ids, node_indices = cache.mortars
  index_range = eachnode(dg)

  @threaded for mortar in eachmortar(dg, cache)
    # Copy solution data from the small elements using "delayed indexing" with
    # a start value and two step sizes to get the correct face and orientation.
    small_indices = node_indices[1, mortar]

    i_small_start, i_small_step_i, i_small_step_j = index_to_start_step_3d(small_indices[1], index_range)
    j_small_start, j_small_step_i, j_small_step_j = index_to_start_step_3d(small_indices[2], index_range)
    k_small_start, k_small_step_i, k_small_step_j = index_to_start_step_3d(small_indices[3], index_range)

    for position in 1:4
      i_small = i_small_start
      j_small = j_small_start
      k_small = k_small_start
      element = neighbor_ids[position, mortar]
      for j in eachnode(dg)
        for i in eachnode(dg)
          for v in eachvariable(equations)
            cache.mortars.u[1, v, position, i, j, mortar] = u[v, i_small, j_small, k_small, element]
          end
          i_small += i_small_step_i
          j_small += j_small_step_i
          k_small += k_small_step_i
        end
        i_small += i_small_step_j
        j_small += j_small_step_j
        k_small += k_small_step_j
      end
    end

    # Buffer to copy solution values of the large element in the correct orientation
    # before interpolating
    u_buffer = cache.u_threaded[Threads.threadid()]
    # temporary buffer for projections
    fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

    # Copy solution of large element face to buffer in the
    # correct orientation
    large_indices = node_indices[2, mortar]

    i_large_start, i_large_step_i, i_large_step_j = index_to_start_step_3d(large_indices[1], index_range)
    j_large_start, j_large_step_i, j_large_step_j = index_to_start_step_3d(large_indices[2], index_range)
    k_large_start, k_large_step_i, k_large_step_j = index_to_start_step_3d(large_indices[3], index_range)

    i_large = i_large_start
    j_large = j_large_start
    k_large = k_large_start
    element = neighbor_ids[5, mortar]
    for j in eachnode(dg)
      for i in eachnode(dg)
        for v in eachvariable(equations)
          u_buffer[v, i, j] = u[v, i_large, j_large, k_large, element]
        end
        i_large += i_large_step_i
        j_large += j_large_step_i
        k_large += k_large_step_i
      end
      i_large += i_large_step_j
      j_large += j_large_step_j
      k_large += k_large_step_j
    end

    # Interpolate large element face data from buffer to small face locations
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, :, mortar),
                            mortar_l2.forward_lower,
                            mortar_l2.forward_lower,
                            u_buffer,
                            fstar_tmp)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, :, mortar),
                            mortar_l2.forward_upper,
                            mortar_l2.forward_lower,
                            u_buffer,
                            fstar_tmp)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 3, :, :, mortar),
                            mortar_l2.forward_lower,
                            mortar_l2.forward_upper,
                            u_buffer,
                            fstar_tmp)
    multiply_dimensionwise!(view(cache.mortars.u, 2, :, 4, :, :, mortar),
                            mortar_l2.forward_upper,
                            mortar_l2.forward_upper,
                            u_buffer,
                            fstar_tmp)
  end

  return nothing
end

function calc_mortar_flux!(surface_flux_values,
                           mesh::P4estMesh{3},
                           nonconservative_terms, equations,
                           mortar_l2::LobattoLegendreMortarL2,
                           surface_integral, dg::DG, cache)
  @unpack neighbor_ids, node_indices = cache.mortars
  @unpack contravariant_vectors = cache.elements
  @unpack fstar_threaded, fstar_tmp_threaded = cache
  index_range = eachnode(dg)

  @threaded for mortar in eachmortar(dg, cache)
    # Choose thread-specific pre-allocated container
    fstar = fstar_threaded[Threads.threadid()]
    fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

    # Get index information on the small elements
    small_indices = node_indices[1, mortar]
    small_direction = indices2direction(small_indices)

    i_small_start, i_small_step_i, i_small_step_j = index_to_start_step_3d(small_indices[1], index_range)
    j_small_start, j_small_step_i, j_small_step_j = index_to_start_step_3d(small_indices[2], index_range)
    k_small_start, k_small_step_i, k_small_step_j = index_to_start_step_3d(small_indices[3], index_range)

    for position in 1:4
      i_small = i_small_start
      j_small = j_small_start
      k_small = k_small_start
      element = neighbor_ids[position, mortar]
      for j in eachnode(dg)
        for i in eachnode(dg)
          # Get the normal direction on the small element.
          # Note, contravariant vectors at interfaces in negative coordinate direction
          # are pointing inwards. This is handled by `get_normal_direction`.
          normal_direction = get_normal_direction(small_direction, contravariant_vectors,
                                                  i_small, j_small, k_small, element)

          calc_mortar_flux!(fstar, mesh, nonconservative_terms, equations,
                            surface_integral, dg, cache,
                            mortar, position, normal_direction,
                            i, j)

          i_small += i_small_step_i
          j_small += j_small_step_i
          k_small += k_small_step_i
        end
        i_small += i_small_step_j
        j_small += j_small_step_j
        k_small += k_small_step_j
      end
    end

    # Buffer to interpolate flux values of the large element to before
    # copying in the correct orientation
    u_buffer = cache.u_threaded[Threads.threadid()]

    mortar_fluxes_to_elements!(surface_flux_values,
                               mesh, equations, mortar_l2, dg, cache,
                               mortar, fstar, u_buffer, fstar_tmp)
  end

  return nothing
end

# Inlined version of the mortar flux computation on small elements for conservation fluxes
@inline function calc_mortar_flux!(fstar,
                                   mesh::P4estMesh{3},
                                   nonconservative_terms::Val{false}, equations,
                                   surface_integral, dg::DG, cache,
                                   mortar_index, position_index, normal_direction,
                                   i_node_index, j_node_index)
  @unpack u = cache.mortars
  @unpack surface_flux = surface_integral

  u_ll, u_rr = get_surface_node_vars(u, equations, dg, position_index, i_node_index, j_node_index, mortar_index)

  flux = surface_flux(u_ll, u_rr, normal_direction, equations)

  # Copy flux to buffer
  set_node_vars!(fstar, flux, equations, dg, i_node_index, j_node_index, position_index)
end

# Inlined version of the mortar flux computation on small elements for conservation fluxes
# with nonconservative terms
@inline function calc_mortar_flux!(fstar,
                                   mesh::P4estMesh{3},
                                   nonconservative_terms::Val{true}, equations,
                                   surface_integral, dg::DG, cache,
                                   mortar_index, position_index, normal_direction,
                                   i_node_index, j_node_index)
  @unpack u = cache.mortars
  surface_flux, nonconservative_flux = surface_integral.surface_flux

  u_ll, u_rr = get_surface_node_vars(u, equations, dg, position_index, i_node_index, j_node_index, mortar_index)

  # Compute conservative flux
  flux = surface_flux(u_ll, u_rr, normal_direction, equations)

  # Compute nonconservative flux and add it to the flux scaled by a factor of 0.5 based on
  # the interpretation of global SBP operators coupled discontinuously via
  # central fluxes/SATs
  noncons = nonconservative_flux(u_ll, u_rr, normal_direction, normal_direction, equations)
  flux_plus_noncons = flux + 0.5 * noncons

  # Copy to buffer
  set_node_vars!(fstar, flux_plus_noncons, equations, dg, i_node_index, j_node_index, position_index)
end

@inline function mortar_fluxes_to_elements!(surface_flux_values,
                                            mesh::P4estMesh{3}, equations,
                                            mortar_l2::LobattoLegendreMortarL2,
                                            dg::DGSEM, cache, mortar, fstar, u_buffer, fstar_tmp)
  @unpack neighbor_ids, node_indices = cache.mortars
  index_range = eachnode(dg)

  # Copy solution small to small
  small_indices   = node_indices[1, mortar]
  small_direction = indices2direction(small_indices)

  for position in 1:4
    element = neighbor_ids[position, mortar]
    for j in eachnode(dg), i in eachnode(dg)
      for v in eachvariable(equations)
        surface_flux_values[v, i, j, small_direction, element] = fstar[v, i, j, position]
      end
    end
  end

  # Project small fluxes to large element.
  multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_lower, mortar_l2.reverse_lower,
    view(fstar, .., 1),
    fstar_tmp)
  add_multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_upper, mortar_l2.reverse_lower,
    view(fstar, .., 2),
    fstar_tmp)
  add_multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_lower, mortar_l2.reverse_upper,
    view(fstar, .., 3),
    fstar_tmp)
  add_multiply_dimensionwise!(
    u_buffer,
    mortar_l2.reverse_upper, mortar_l2.reverse_upper,
    view(fstar, .., 4),
    fstar_tmp)

  # The flux is calculated in the outward direction of the small elements,
  # so the sign must be switched to get the flux in outward direction
  # of the large element.
  # The contravariant vectors of the large element (and therefore the normal
  # vectors of the large element as well) are four times as large as the
  # contravariant vectors of the small elements. Therefore, the flux needs
  # to be scaled by a factor of 4 to obtain the flux of the large element.
  u_buffer .*= -4

  # Copy interpolated flux values from buffer to large element face in the
  # correct orientation.
  # Note that the index of the small sides will always run forward but
  # the index of the large side might need to run backwards for flipped sides.
  large_element = neighbor_ids[5, mortar]
  large_indices  = node_indices[2, mortar]
  large_direction = indices2direction(large_indices)
  large_surface_indices = surface_indices(large_indices)

  i_large_start, i_large_step_i, i_large_step_j = index_to_start_step_3d(large_surface_indices[1], index_range)
  j_large_start, j_large_step_i, j_large_step_j = index_to_start_step_3d(large_surface_indices[2], index_range)

  # Note that the indices of the small sides will always run forward but
  # the large indices might need to run backwards for flipped sides.
  i_large = i_large_start
  j_large = j_large_start
  for j in eachnode(dg)
    for i in eachnode(dg)
      for v in eachvariable(equations)
        surface_flux_values[v, i_large, j_large, large_direction, large_element] = u_buffer[v, i, j]
      end
      i_large += i_large_step_i
      j_large += j_large_step_i
    end
    i_large += i_large_step_j
    j_large += j_large_step_j
  end

  return nothing
end

function calc_surface_integral!(du, u,
                                mesh::P4estMesh{3},
                                equations,
                                surface_integral::SurfaceIntegralWeakForm,
                                dg::DGSEM, cache)
  @unpack boundary_interpolation = dg.basis
  @unpack surface_flux_values = cache.elements

  # Note that all fluxes have been computed with outward-pointing normal vectors.
  # Access the factors only once before beginning the loop to increase performance.
  # We also use explicit assignments instead of `+=` to let `@muladd` turn these
  # into FMAs (see comment at the top of the file).
  factor_1 = boundary_interpolation[1,          1]
  factor_2 = boundary_interpolation[nnodes(dg), 2]
  @threaded for element in eachelement(dg, cache)
    for m in eachnode(dg), l in eachnode(dg)
      for v in eachvariable(equations)
        # surface at -x
        du[v, 1,          l, m, element] = (
          du[v, 1,          l, m, element] + surface_flux_values[v, l, m, 1, element] * factor_1)

        # surface at +x
        du[v, nnodes(dg), l, m, element] = (
          du[v, nnodes(dg), l, m, element] + surface_flux_values[v, l, m, 2, element] * factor_2)

        # surface at -y
        du[v, l, 1,          m, element] = (
          du[v, l, 1,          m, element] + surface_flux_values[v, l, m, 3, element] * factor_1)

        # surface at +y
        du[v, l, nnodes(dg), m, element] = (
          du[v, l, nnodes(dg), m, element] + surface_flux_values[v, l, m, 4, element] * factor_2)

        # surface at -z
        du[v, l, m, 1,          element] = (
          du[v, l, m, 1,          element] + surface_flux_values[v, l, m, 5, element] * factor_1)

        # surface at +z
        du[v, l, m, nnodes(dg), element] = (
          du[v, l, m, nnodes(dg), element] + surface_flux_values[v, l, m, 6, element] * factor_2)
      end
    end
  end

  return nothing
end

end # @muladd


text:
using DemoPackageTEH
using Test

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


text:
module  TestSibship

using PopGenCore
using PopGenSims
using Test

cats = @nancycats ;

@testset "Sibship Simulations" begin
    @testset "all sibship" begin
        sims_all = simulatekin(cats, fullsib = 50, halfsib = 50, unrelated = 50, parentoffspring = 50)
        @test sims_all isa PopData
        @test sims_all.metadata.samples == 400
        @test length(sims_all.genodata.name) == 3600
    end

    @testset "unrelated" begin
        sims_un = simulatekin(cats, unrelated = 50)
        @test sims_un isa PopData
        @test sims_un.metadata.samples == 100
        @test length(sims_un.genodata.name) == 900
    end

    @testset "halfsib" begin
        sims_hs = simulatekin(cats, halfsib = 50)
        @test sims_hs isa PopData
        @test sims_hs.metadata.samples == 100
        @test length(sims_hs.genodata.name) == 900
    end

    @testset "fullsib" begin
        sims_fs = simulatekin(cats, fullsib = 50)
        @test sims_fs isa PopData
        @test sims_fs.metadata.samples == 100
        @test length(sims_fs.genodata.name) == 900
    end

    @testset "parent offspring" begin
        sims_po = simulatekin(cats, parentoffspring = 50)
        @test sims_po isa PopData
        @test sims_po.metadata.samples == 100
        @test length(sims_po.genodata.name) == 900
    end
end

end # module


text:
module OrihimeTests

include("test_samples.jl")
include("test_readablestring.jl")

end  # module OrihimeTests


text:
module UrlDownload

using ProgressMeter
using HTTP

export urldownload

const ext2sym = Dict(
    r".feather" => :FEATHER,
    r".[ct]sv" => :CSV,
    r".jpe?g" => :PIC,
    r".png" => :PIC,
    r".json" => :JSON,
    r".tga" => :PIC,
    r".gif" => :PIC,
    r".bmp" => :PIC,
    r".pcx" => :PIC
)

const sym2func = Dict(
    :FEATHER => (x; kw...) -> load_feather(x; kw...),
    :PIC => (x; kw...) -> load_pic(x; kw...),
    :CSV => (x; kw...) -> load_csv(x; kw...),
    :TSV => (x; kw...) -> load_csv(x; kw...),
    :JSON => (x; kw...) -> load_json(x; kw...)
)

const Compressor = Dict(
    :gzip => (lib = :CodecZlib, stream = :GzipDecompressorStream, transcode = :GzipDecompressor),
    :zstd => (lib = :CodecZstd, stream = :ZstdDecompressorStream, transcode = :ZstdDecompressor),
    :xz => (lib = :CodecXz, stream = :XzDecompressorStream, transcode = :XzDecompressor),
    :lz4 => (lib = :CodecLz4, stream = :LZ4FrameDecompressorStream, transcode = :LZ4FrameDecompressor),
    :bzip2 => (lib = :CodecBzip2, stream = :Bzip2DecompressorStream, transcode = :Bzip2Decompressor)
)

abstract type Resource end

struct File <: Resource
    name::String
end

struct URL <: Resource
    name::String
end

"""
    f_str(name)

Use this macro to explicitly show that downloaded url is actually local file resource. It is useful if resource type autodetection fails.

# Example
```
using UrlDownload: @f_str
url = f"/tmp/data"
res = urldownload(url)

# Alternatively
using UrlDownload: File
url = File("/tmp/data")
res = urldownload(url)
```
"""
macro f_str(p)
    return quote
        File($(esc(p)))
    end
end

"""
    u_str(name)

Use this macro to explicitly show that downloaded url is remote http resource. It is useful if resource type autodetection fails.

# Example
```
using UrlDownload: @u_str
url = u"https://example.com/data.csv"
res = urldownload(url)

# Alternatively
using UrlDownload: URL
url = URL("https://example.com/data.csv")
res = urldownload(url)
```
"""
macro u_str(p)
    return quote
        URL($(esc(p)))
    end
end

function load_feather(buf; kw...)
    lib = checked_import(:Feather)
    return Base.invokelatest(lib.read, buf)
end

function load_csv(buf; kw...)
    lib = checked_import(:CSV)
    return Base.invokelatest(lib.File, buf; kw...)
end

function load_pic(buf; kw...)
    lib = checked_import(:ImageMagick)
    return Base.invokelatest(lib.load_, _getdata(buf))
end

function load_json(buf; kw...)
    lib = checked_import(:JSON3)
    return Base.invokelatest(lib.read, _getdata(buf))
end

# Borrowed directly from FileIO

# This one is outdated.
# is_installed(pkg::Symbol) = get(Pkg.installed(), string(pkg), nothing) != nothing

function _findmod(f::Symbol)
    for (u,v) in Base.loaded_modules
        (Symbol(v) == f) && return u, v
    end
    nothing, nothing
end

function loadmod(pkg::Symbol)
    @eval Base.__toplevel__  import $pkg
    return _findmod(pkg)
end

function checked_import(pkg::Symbol)
    if isdefined(Main, pkg)
        m1 = getfield(Main, pkg)
        isa(m1, Module) && return m1
    end
    mid, m = _findmod(pkg)
    mid == nothing || return m
    _, m = loadmod(pkg)
    return m
end

getext(url) = splitext(url)
getext(url::IO) = splitext(url.name[2:end-1])
getext(url::Resource) = getext(url.name)

function datatype(url)
    _, ext = getext(url)
    ext = lowercase(ext)
    for (k, v) in ext2sym
        occursin(k, ext) && return v
    end
end

# this set of functions is needed for distinguish between raw data and zipstreams,
# we do not want to wrap/unwrap zipstreams.
createbuffer(data::Vector{UInt8}) = IOBuffer(data)
createbuffer(data) = data

_getdata(buf::IOBuffer) = buf.data
_getdata(buf) = buf

function wrapdata(url, data, format, parser, error_on_undetected_format = true; kw...)
    if parser === nothing
        buf = createbuffer(data)
        dtype = format === nothing ? datatype(url) : format

        if dtype in keys(sym2func)
            return sym2func[dtype](buf; kw...)
        else
            dtype = dtype === nothing ? :unknown : dtype
            if error_on_undetected_format
                @error "Data format $dtype is not supported."
            else
                @warn "Data format $dtype is not supported."
                return data
            end
        end
    else
        return parser(data; kw...)
    end
end

maybezip(url) = occursin(".zip", url)
maybezip(url::Resource) = occursin(".zip", url.name)
maybezip(url::IO) = occursin(".zip", url.name)

# Check the file format of a stream.
# Taken from TableReader.jl
function checkformat(data, url)
    if length(data) >= 6
        magic = data[1:6]
    else
        error("Incomplete data received")
    end
    if magic[1:6] == b"\xFD\x37\x7A\x58\x5A\x00"
        return :xz
    elseif magic[1:2] == b"\x1f\x8b"
        return :gzip
    elseif magic[1:4] == b"\x28\xb5\x2f\xfd"
        return :zstd
    elseif magic[1:3] == b"\x42\x5A\x68"
        return :bzip2
    elseif magic[1:4] == b"\x04\x22\x4D\x18"
        return :lz4
    end

    # there is a possibility of false positive, but
    # better to be sorry
    if maybezip(url) &
        ((magic[1:4] == b"\x50\x4b\x03\x04") | (magic[1:4] == b"\x50\x4b\x05\x06") | (magic[1:4] == b"\x50\x4b\x07\x08"))
            return :zip
    end

    # we are giving up
    return :none
end

wrap_httpkw(httpkw::Pair) = [httpkw]
wrap_httpkw(httpkw) = httpkw

saveraw(io::Nothing, data) = nothing
saveraw(io, data) = write(io, data)
function saveraw(io::AbstractString, data)
    iostream = open(io, "w")
    saveraw(iostream, data)
    close(iostream)
end

get_data(url::IO, headers, httpkw, progress, update_period) = read(url)
get_data(url::AbstractString, headers, httpkw, progress, update_period) =
    get_data(autodetect_uri_type(url), headers, httpkw, progress, update_period)
get_data(url::File, headers, httpkw, progress, update_period) = read(url.name)
function get_data(url::URL, headers, httpkw, progress, update_period)
    local body
    HTTP.open("GET", url.name, headers; httpkw...) do stream
        resp = startread(stream)
        eof(stream) && return # nothing to process yet

        body = UInt8[]
        if progress
            total_bytes = floor(parse(Float64, HTTP.header(resp, "Content-Length", "NaN")))
            total_bytes = isnan(total_bytes) ? typemax(Int) : Int(total_bytes)
            p = Progress(total_bytes, update_period)
        end

        while !eof(stream)
            append!(body, readavailable(stream))
            if progress
                update!(p, length(body))
            end
        end
    end

    return body
end

autodetect_uri_type(url) = url
function autodetect_uri_type(url::AbstractString)
    m = match(r"^https?://.*", url)
    if m === nothing
        return File(url)
    else
        return URL(url)
    end
end

"""
    urldownload(url, progress = false;
                parser = nothing, format = nothing, save_raw = nothing,
                compress = :auto, multifiles = false, headers = HTTP.Header[],
                httpkw = Pair[], update_period = 1, kw...)

Download file from the corresponding `url` in memory and process it to the necessary data structure.

# Arguments
* `url`: url of download
* `progress`: show `ProgressMeter`, by default it is not shown
* `parser`: custom parser, function that should accept one positional argument of the type `Vector{UInt8}` and optional keyword arguments and return necessary data structure. If parser is set than it overrides all other settings, such as `format`. If parser is not set, than internal parsers are used for data process.
* `format`: one of the fixed formats (:CSV, :PIC, :FEATHER, :JSON), if set overrides autodetection mechanism.
* `save_raw`: if set to `String` or `IO` then downloaded raw data is stored in corresponding stream.
* `compress`: :auto by default, can be one of :none, :xz, :gzip, :bzip2, :lz4, :zstd, :zip. Determines whether file is compressed and compression type. Decompressed data is processed either by custom `parser` or by internal parser. By default for any compression type except of `:zip` internal parser is `CSV.File`, for `:zip` usual rules applies. If `compress` is `:none` than custom parser should decompress data on its own.
* `multifiles`: `false` by default, for `:zip` compressed data defines, whether process only first file inside archive or return an array of decompressed and processed objects.
* `headers`: `HTTP.jl` arguments that set http header of the request.
* `httpkw`: `HTTP.jl` additional keyword arguments that is passed to the `GET` function. Should be supplied as a vector of pairs.
* `update_period`: period of `ProgressMeter` update, by default 1 sec
* `kw...`: any keyword arguments that should be passed to the data parser.
"""
function urldownload(url, progress = false;
        parser = nothing, format = nothing,
        compress = :auto, multifiles = false,
        save_raw = nothing,
        headers = HTTP.Header[],
        update_period = 1, httpkw = Pair[], kw...)

    httpkw = wrap_httpkw(httpkw)
    body = get_data(url, headers, httpkw, progress, update_period)
    saveraw(save_raw, body)

    compress = compress == :auto ? checkformat(body, url) : compress

    if compress == :none
        # skip unzipping entirely, it's parser responisbility to process the data
        wrapdata(url, body, format, parser; kw...)
    elseif compress == :zip
        zlib = checked_import(:ZipFile)
        zread = Base.invokelatest(getfield(zlib, :Reader), IOBuffer(body)).files
        if multifiles
            return [wrapdata(z.name, z, format, parser, false; kw...) for z in zread]
        else
            if length(zread) > 1
                @warn "More than one file in zip archive, returning first."
            elseif length(zread) == 0
                @error "Zip archive is empty."
            end
            zread = zread[1]
            # This one can easily fail for non csv files, nothing I can do about it
            wrapdata(zread.name, zread, format, parser; kw...)
        end
    elseif compress in keys(Compressor)
        # it's one of the TranscodingStreams.jl streams, not much to do here,
        # defaults to CSV/custom parser
        if parser === nothing
            lib = checked_import(Compressor[compress].lib)
            stream = getfield(lib, Compressor[compress].stream)
            csvlib = checked_import(:CSV)
            csvlibfile = getfield(csvlib, :File)
            return Base.invokelatest(csvlibfile, Base.invokelatest(stream, IOBuffer(body)); kw...)
        else
            lib = checked_import(Compressor[compress].lib)
            transcoder = getfield(lib, Compressor[compress].transcode)
            data = Base.invokelatest(transcode, transcoder, body)
            return parser(data; kw...)
        end
    else
        error("Unknown compress format: $compress")
    end
end

end # module


text:
#=

  "ClaDS: Cladogenetic XX" submodule package

=#

module ClaDS

using Random: randexp
using Distributions: Bernoulli, LogNormal, Weights, sample, Normal, InverseGamma
using RCall: @rput, @rget, reval
using Statistics: mean, var
using StringDistances: Levenshtein
# other submodules dependencies
using ..Utils

# files
include("clads/tree_class.jl")
include("clads/clads_output.jl")
include("clads/enhance_tree_utils.jl")
include("clads/sim_ClaDS_utils.jl")
include("clads/sim_ClaDS.jl")
include("clads/load_tree.jl")
include("clads/plot_ClaDS.jl")
include("clads/infer_ClaDS2.jl")
include("clads/mcmc_ClaDS2_utils.jl")
include("clads/sample_fractions.jl")
include("clads/LTT.jl")
include("clads/n_tips.jl")
include("clads/Tree_utils.jl")
include("clads/branch_lengths.jl")
include("clads/RTT.jl")
include("clads/enhance_tree.jl")
include("clads/rates.jl")
include("clads/gibbs.jl")
include("clads/export_in_R.jl")
include("clads/prune.jl")

end # module ClaDS

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA