text:
# Lower-level interface to the bzip2 library.

const bzip2 = "libbz2"

# Constants
# ---------

const BZ_RUN              = Cint(0)
const BZ_FLUSH            = Cint(1)
const BZ_FINISH           = Cint(2)

const BZ_OK               = Cint(0)
const BZ_RUN_OK           = Cint(1)
const BZ_FLUSH_OK         = Cint(2)
const BZ_FINISH_OK        = Cint(3)
const BZ_STREAM_END       = Cint(4)
const BZ_SEQUENCE_ERROR   = Cint(-1)
const BZ_PARAM_ERROR      = Cint(-2)
const BZ_MEM_ERROR        = Cint(-3)
const BZ_DATA_ERROR       = Cint(-4)
const BZ_DATA_ERROR_MAGIC = Cint(-5)
const BZ_IO_ERROR         = Cint(-6)
const BZ_UNEXPECTED_EOF   = Cint(-7)
const BZ_OUTBUFF_FULL     = Cint(-8)
const BZ_CONFIG_ERROR     = Cint(-9)

# BZStream
# --------

type BZStream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in_lo32::Cuint
    total_in_hi32::Cuint

    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out_lo32::Cuint
    total_out_hi32::Cuint

    state::Ptr{Void}

    bzalloc::Ptr{Void}
    bzfree::Ptr{Void}
    opaque::Ptr{Void}

    function BZStream()
        bzstream = new()
        bzstream.next_in        = C_NULL
        bzstream.avail_in       = 0
        bzstream.total_in_lo32  = 0
        bzstream.total_in_hi32  = 0
        bzstream.next_out       = C_NULL
        bzstream.avail_out      = 0
        bzstream.total_out_lo32 = 0
        bzstream.total_out_hi32 = 0
        bzstream.state          = C_NULL
        bzstream.bzalloc        = C_NULL
        bzstream.bzfree         = C_NULL
        bzstream.opaque         = C_NULL
        return bzstream
    end
end

# Functions
# ---------

function version()
    return unsafe_string(ccall((:BZ2_bzlibVersion, bzip2), Ptr{UInt8}, ()))
end

function init_compress!(
        bzstream::BZStream,
        blocksize100k::Integer,
        verbosity::Integer,
        workfactor::Integer)
    return ccall(
        (:BZ2_bzCompressInit, bzip2),
        Cint,
        (Ref{BZStream}, Cint, Cint, Cint),
        bzstream, blocksize100k, verbosity, workfactor)
end

function compress!(bzstream::BZStream, action::Integer)
    return ccall(
        (:BZ2_bzCompress, bzip2),
        Cint,
        (Ref{BZStream}, Cint),
        bzstream, action)
end

function end_compress!(bzstream::BZStream)
    return ccall((:BZ2_bzCompressEnd, bzip2), Cint, (Ref{BZStream},), bzstream)
end

function init_decompress!(bzstream::BZStream, verbosity::Integer, small::Integer)
    return ccall(
        (:BZ2_bzDecompressInit, bzip2),
        Cint,
        (Ref{BZStream}, Cint, Cint),
        bzstream, verbosity, small)
end

function decompress!(bzstream::BZStream)
    return ccall((:BZ2_bzDecompress, bzip2), Cint, (Ref{BZStream},), bzstream)
end

function end_decompress!(bzstream::BZStream)
    return ccall((:BZ2_bzDecompressEnd, bzip2), Cint, (Ref{BZStream},), bzstream)
end


text:
function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{CopperPlatePowerModel},
)
    buses = PSY.get_components(PSY.Bus, sys)
    bus_count = length(buses)

    if get_balance_slack_variables(psi_container.settings)
        add_slacks!(psi_container, CopperPlatePowerModel)
    end
    copper_plate(psi_container, :nodal_balance_active, bus_count)
    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{AreaBalancePowerModel},
)
    area_mapping = PSY.get_aggregation_topology_mapping(PSY.Area, sys)
    branches = get_available_components(PSY.Branch, sys)
    if get_balance_slack_variables(psi_container.settings)
        throw(IS.ConflictingInputsError("Slack Variables are not compatible with AreaBalancePowerModel"))
    end

    area_balance(psi_container, :nodal_balance_active, area_mapping, branches)
    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{StandardPTDFModel},
)
    buses = PSY.get_components(PSY.Bus, sys)
    ac_branches = get_available_components(PSY.ACBranch, sys)
    ptdf = get_PTDF(psi_container)

    if isnothing(ptdf)
        throw(ArgumentError("no PTDF matrix supplied"))
    end

    if get_balance_slack_variables(psi_container.settings)
        add_slacks!(psi_container, StandardPTDFModel)
    end

    ptdf_networkflow(psi_container, ac_branches, buses, :nodal_balance_active, ptdf)

    dc_branches = get_available_components(PSY.DCBranch, sys)
    dc_branch_types = typeof.(dc_branches)
    for btype in Set(dc_branch_types)
        typed_dc_branches = IS.FlattenIteratorWrapper(
            btype,
            Vector([[b for b in dc_branches if typeof(b) == btype]]),
        )
        flow_variables!(psi_container, StandardPTDFModel, typed_dc_branches)
    end
    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{T};
    instantiate_model = instantiate_nip_expr_model,
) where {T <: PM.AbstractPowerModel}
    if T in UNSUPPORTED_POWERMODELS
        throw(ArgumentError("$(T) formulation is not currently supported in PowerSimulations"))
    end

    if get_balance_slack_variables(psi_container.settings)
        add_slacks!(psi_container, T)
    end

    @debug "Building the $T network with $instantiate_model method"
    powermodels_network!(psi_container, T, sys, instantiate_model)
    add_pm_var_refs!(psi_container, T, sys)
    add_pm_con_refs!(psi_container, T, sys)
    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{T};
    instantiate_model = instantiate_bfp_expr_model,
) where {T <: PM.AbstractBFModel}
    if T in UNSUPPORTED_POWERMODELS
        throw(ArgumentError("$(T) formulation is not currently supported in PowerSimulations"))
    end

    get_balance_slack_variables(psi_container.settings) && add_slacks!(psi_container, T)

    @debug "Building the $T network with $instantiate_model method"
    powermodels_network!(psi_container, T, sys, instantiate_model)
    add_pm_var_refs!(psi_container, T, sys)
    add_pm_con_refs!(psi_container, T, sys)
    return
end

function construct_network!(
    psi_container::PSIContainer,
    sys::PSY.System,
    ::Type{T};
    instantiate_model = instantiate_vip_expr_model,
) where {T <: PM.AbstractIVRModel}
    if T in UNSUPPORTED_POWERMODELS
        throw(ArgumentError("$(T) formulation is not currently supported in PowerSimulations"))
    end

    if get_balance_slack_variables(psi_container.settings)
        add_slacks!(psi_container, T)
    end

    @debug "Building the $T network with $instantiate_model method"
    powermodels_network!(psi_container, T, sys, instantiate_model)
    add_pm_var_refs!(psi_container, T, sys)
    add_pm_con_refs!(psi_container, T, sys)
    return
end

ext: jl


text:
# ## Turbulence Closures
# In `turbulence.jl` we specify turbulence closures. Currently,
# pointwise models of the eddy viscosity/eddy diffusivity type are
# supported for turbulent shear and tracer diffusivity. Methods currently supported
# are:\
# [`ConstantViscosityWithDivergence`](@ref constant-viscosity)\
# [`SmagorinskyLilly`](@ref smagorinsky-lilly)\
# [`Vreman`](@ref vreman)\
# [`AnisoMinDiss`](@ref aniso-min-diss)\

#md # !!! note
#md #     Usage: This is a quick-ref guide to using turbulence models as a subcomponent
#md #     of `AtmosModel` \
#md #     $\nu$ is the kinematic viscosity, $C_smag$ is the Smagorinsky Model coefficient,
#md #     `turbulence=ConstantViscosityWithDivergence(ν)`\
#md #     `turbulence=SmagorinskyLilly(C_smag)`\
#md #     `turbulence=Vreman(C_smag)`\
#md #     `turbulence=AnisoMinDiss(C_poincare)`

using DocStringExtensions
using CLIMAParameters.Atmos.SubgridScale: inv_Pr_turb
export ConstantViscosityWithDivergence, SmagorinskyLilly, Vreman, AnisoMinDiss
export turbulence_tensors

# ### Abstract Type
# We define a `TurbulenceClosure` abstract type and
# default functions for the generic turbulence closure
# which will be overloaded with model specific functions.
abstract type TurbulenceClosure end

vars_state(::TurbulenceClosure, FT) = @vars()
vars_aux(::TurbulenceClosure, FT) = @vars()

function atmos_init_aux!(
    ::TurbulenceClosure,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
) end
function atmos_nodal_update_aux!(
    ::TurbulenceClosure,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function gradvariables!(
    ::TurbulenceClosure,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function diffusive!(
    ::TurbulenceClosure,
    ::Orientation,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end

"""
    ν, D_t, τ = turbulence_tensors(::TurbulenceClosure, orientation::Orientation, param_set::AbstractParameterSet, state::Vars, diffusive::Vars, aux::Vars, t::Real)

    Compute the kinematic viscosity (`ν`), the diffusivity (`D_t`) and SGS momentum flux tensor (`τ`).
"""
function turbulence_tensors end

turbulence_tensors(atmos::AtmosModel, args...) =
    turbulence_tensors(atmos.turbulence, atmos, args...)

turbulence_tensors(m::TurbulenceClosure, atmos::AtmosModel, args...) =
    turbulence_tensors(m, atmos.orientation, atmos.param_set, args...)

# We also provide generic math functions for use within the turbulence closures,
# commonly used quantities such as the [principal tensor invariants](@ref tensor-invariants), handling of
# [symmetric tensors](@ref symmetric-tensors) and [tensor norms](@ref tensor-norms)are addressed.

# ### [Pricipal Invariants](@id tensor-invariants)
"""
    principal_invariants(X)

Calculates principal invariants of a tensor `X`. Returns 3 element tuple containing the invariants.
"""
function principal_invariants(X)
    first = tr(X)
    second = (first^2 - tr(X .^ 2)) / 2
    third = det(X)
    return (first, second, third)
end

# ### [Symmetrize](@id symmetric-tensors)
# ```math
# \frac{\mathrm{X} + \mathrm{X}^{T}}{2}
# ```
"""
    symmetrize(X)

Compute `(X + X')/2`, returning a `SHermitianCompact` object.
"""
function symmetrize(X::StaticArray{Tuple{3, 3}})
    SHermitianCompact(SVector(
        X[1, 1],
        (X[2, 1] + X[1, 2]) / 2,
        (X[3, 1] + X[1, 3]) / 2,
        X[2, 2],
        (X[3, 2] + X[2, 3]) / 2,
        X[3, 3],
    ))
end

# ### [2-Norm](@id tensor-norms)
# ```math
# \sum_{i,j} S_{ij}^2
# ```
"""
    norm2(X)

Compute
```math
\\sum_{i,j} S_{ij}^2
```
"""
function norm2(X::SMatrix{3, 3, FT}) where {FT}
    abs2(X[1, 1]) +
    abs2(X[2, 1]) +
    abs2(X[3, 1]) +
    abs2(X[1, 2]) +
    abs2(X[2, 2]) +
    abs2(X[3, 2]) +
    abs2(X[1, 3]) +
    abs2(X[2, 3]) +
    abs2(X[3, 3])
end
function norm2(X::SHermitianCompact{3, FT, 6}) where {FT}
    abs2(X[1, 1]) +
    2 * abs2(X[2, 1]) +
    2 * abs2(X[3, 1]) +
    abs2(X[2, 2]) +
    2 * abs2(X[3, 2]) +
    abs2(X[3, 3])
end

# ### [Strain-rate Magnitude](@id strain-rate-magnitude)
# By definition, the strain-rate magnitude, as defined in
# standard turbulence modelling is computed such that
# ```math
# |\mathrm{S}| = \sqrt{2 \sum_{i,j} \mathrm{S}_{ij}^2}
# ```
"""
    strain_rate_magnitude(S)

Compute
```math
|S| = \\sqrt{2\\sum_{i,j} S_{ij}^2}
```
"""
function strain_rate_magnitude(S::SHermitianCompact{3, FT, 6}) where {FT}
    return sqrt(2 * norm2(S))
end

# ### [Constant Viscosity Model](@id constant-viscosity)
# `ConstantViscosityWithDivergence` requires a user to specify the constant viscosity (kinematic)
# and appropriately computes the turbulent stress tensor based on this term. Diffusivity can be
# computed using the turbulent Prandtl number for the appropriate problem regime.
# ```math
# \tau = - 2 \nu \mathrm{S}
# ```
"""
    ConstantViscosityWithDivergence <: TurbulenceClosure

Turbulence with constant dynamic viscosity (`ρν`).
Divergence terms are included in the momentum flux tensor.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ConstantViscosityWithDivergence{FT} <: TurbulenceClosure
    "Dynamic Viscosity [kg/m/s]"
    ρν::FT
end

vars_gradient(::ConstantViscosityWithDivergence, FT) = @vars()
vars_diffusive(::ConstantViscosityWithDivergence, FT) =
    @vars(S::SHermitianCompact{3, FT, 6})

function diffusive!(
    ::ConstantViscosityWithDivergence,
    ::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
end

function turbulence_tensors(
    m::ConstantViscosityWithDivergence,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    ν = m.ρν / state.ρ
    D_t = ν * _inv_Pr_turb
    τ = (-2 * ν) * S + (2 * ν / 3) * tr(S) * I
    return ν, D_t, τ
end

# ### [Smagorinsky-Lilly](@id smagorinsky-lilly)
# The Smagorinsky turbulence model, with Lilly's correction to
# stratified atmospheric flows, is included in CLIMA.
# The input parameter to this model is the Smagorinsky coefficient.
# For atmospheric flows, the coefficient `C_smag` typically takes values between
# 0.15 and 0.23. Flow dependent `C_smag` are currently not supported (e.g. Germano's
# extension). The Smagorinsky-Lilly model does not contain explicit filtered terms.
# ```math
# \nu = (C_{s} \mathrm{f}_{b} \Delta)^2 \sqrt{|\mathrm{S}|}
# ```
# with the stratification correction term
# ```math
# \mathrm{f}_{b}^{2} = \sqrt{1 - \frac{\mathrm{Ri}}{\mathrm{Pr}_{t}}}
# ```\
# $\mathrm{Ri}$ and $\mathrm{Pr}_{t}$ are the Richardson and
# turbulent Prandtl numbers respectively.  $\Delta$ is the mixing length in the
# relevant coordinate direction. We use the DG metric terms to determine the
# local effective resolution (see `src/Mesh/Geometry.jl`), and modify the vertical lengthscale by the
# stratification correction factor $\mathrm{f}_{b}$.
"""
    SmagorinskyLilly <: TurbulenceClosure

See § 1.3.2 in CliMA documentation

article{doi:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2,
  author = {Smagorinksy, J.},
  title = {General circulation experiments with the primitive equations},
  journal = {Monthly Weather Review},
  volume = {91},
  number = {3},
  pages = {99-164},
  year = {1963},
  doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  URL = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  eprint = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2}
  }

article{doi:10.1111/j.2153-3490.1962.tb00128.x,
  author = {LILLY, D. K.},
  title = {On the numerical simulation of buoyant convection},
  journal = {Tellus},
  volume = {14},
  number = {2},
  pages = {148-172},
  doi = {10.1111/j.2153-3490.1962.tb00128.x},
  url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2153-3490.1962.tb00128.x},
  eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2153-3490.1962.tb00128.x},
  year = {1962}
  }

Brunt-Vaisala frequency N² defined as in equation (1b) in
  Durran, D.R. and J.B. Klemp, 1982:
  On the Effects of Moisture on the Brunt-Väisälä Frequency.
  J. Atmos. Sci., 39, 2152–2158,
  https://doi.org/10.1175/1520-0469(1982)039<2152:OTEOMO>2.0.CO;2

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SmagorinskyLilly{FT} <: TurbulenceClosure
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end

vars_aux(::SmagorinskyLilly, FT) = @vars(Δ::FT)
vars_gradient(::SmagorinskyLilly, FT) = @vars(θ_v::FT)
vars_diffusive(::SmagorinskyLilly, FT) =
    @vars(S::SHermitianCompact{3, FT, 6}, N²::FT)

function atmos_init_aux!(
    ::SmagorinskyLilly,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end

function gradvariables!(
    m::SmagorinskyLilly,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end

function diffusive!(
    ::SmagorinskyLilly,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::SmagorinskyLilly,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    normS = strain_rate_magnitude(S)
    k̂ = vertical_unit_vector(orientation, param_set, aux)

    # squared buoyancy correction
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(FT(1) - Richardson * _inv_Pr_turb, FT(0), FT(1)))
    ν₀ = normS * (m.C_smag * aux.turbulence.Δ)^2 + FT(1e-5)
    ν = SVector{3, FT}(ν₀, ν₀, ν₀)
    ν_v = k̂ .* dot(ν, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

# ### [Vreman Model](@id vreman)
# Vreman's turbulence model for anisotropic flows, which provides a
# less dissipative solution (specifically in the near-wall and transitional regions)
# than the Smagorinsky-Lilly method. This model
# relies of first derivatives of the velocity vector (i.e., the gradient tensor).
# By design, the Vreman model handles transitional as well as fully turbulent flows adequately.
# The input parameter to this model is the Smagorinsky coefficient - the coefficient is modified
# within the model functions to account for differences in model construction.
# #### Equations
# ```math
# \nu = 2.5 \mathrm{C}_{smag} \sqrt{\frac{\mathrm{B}_{\beta}}{\alpha_{i}\alpha_{j}}}
# ```
# where $\mathrm{B}_{\beta}$ and $\alpha$ are functions of the velocity
# gradient tensor terms.
"""
    Vreman{FT} <: TurbulenceClosure

  §1.3.2 in CLIMA documentation
Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Smagorinsky coefficient
is specified and used to compute the equivalent Vreman coefficient.

1) ν_e = √(Bᵦ/(αᵢⱼαᵢⱼ)) where αᵢⱼ = ∂uⱼ∂uᵢ with uᵢ the resolved scale velocity component.
2) βij = Δ²αₘᵢαₘⱼ
3) Bᵦ = β₁₁β₂₂ + β₂₂β₃₃ + β₁₁β₃₃ - β₁₂² - β₁₃² - β₂₃²
βᵢⱼ is symmetric, positive-definite.
If Δᵢ = Δ, then β = Δ²αᵀα

@article{Vreman2004,
  title={An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications},
  author={Vreman, AW},
  journal={Physics of fluids},
  volume={16},
  number={10},
  pages={3670--3681},
  year={2004},
  publisher={AIP}
}

# Fields

$(DocStringExtensions.FIELDS)
"""
struct Vreman{FT} <: TurbulenceClosure
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end
vars_aux(::Vreman, FT) = @vars(Δ::FT)
vars_gradient(::Vreman, FT) = @vars(θ_v::FT)
vars_diffusive(::Vreman, FT) = @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)

function atmos_init_aux!(::Vreman, ::AtmosModel, aux::Vars, geom::LocalGeometry)
    aux.turbulence.Δ = lengthscale(geom)
end
function gradvariables!(
    m::Vreman,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function diffusive!(
    ::Vreman,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.turbulence.∇u = ∇transform.u
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::Vreman,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    α = diffusive.turbulence.∇u
    S = symmetrize(α)
    k̂ = vertical_unit_vector(orientation, param_set, aux)

    normS = strain_rate_magnitude(S)
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * _inv_Pr_turb, 0, 1))

    β = f_b² * (aux.turbulence.Δ)^2 * (α' * α)
    Bβ = principal_invariants(β)[2]

    ν₀ = m.C_smag^2 * FT(2.5) * sqrt(abs(Bβ / (norm2(α) + eps(FT))))

    ν = SVector{3, FT}(ν₀, ν₀, ν₀)
    ν_v = k̂ .* dot(ν, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

# ### [Anisotropic Minimum Dissipation](@id aniso-min-diss)
# This method is based Vreugdenhil and Taylor's minimum-dissipation eddy-viscosity model.
# The principles of the Rayleigh quotient minimizer are applied to the energy dissipation terms in the
# conservation equations, resulting in a maximum dissipation bound, and a model for
# eddy viscosity and eddy diffusivity.
# ```math
# \nu_e = (\mathrm{C}\delta)^2  \mathrm{max}\left[0, - \frac{\hat{\partial}_k \hat{u}_{i} \hat{\partial}_k \hat{u}_{j} \mathrm{\hat{S}}_{ij}}{\hat{\partial}_p \hat{u}_{q}} \right]
# ```
"""
    AnisoMinDiss{FT} <: TurbulenceClosure

  §1.3.2 in CLIMA documentation
Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Poincare coefficient
is specified and used to compute the equivalent AnisoMinDiss coefficient (computed as the solution to the
eigenvalue problem for the Laplacian operator).

@article{
doi:10.1063/1.5037039,
author = {Vreugdenhil,Catherine A.  and Taylor,John R. },
title = {Large-eddy simulations of stratified plane Couette flow using the anisotropic minimum-dissipation model},
journal = {Physics of Fluids},
volume = {30},
number = {8},
pages = {085104},
year = {2018},
doi = {10.1063/1.5037039},
URL = {
        https://doi.org/10.1063/1.5037039
},
}

Fields
$(DocStringExtensions.FIELDS)
"""
struct AnisoMinDiss{FT} <: TurbulenceClosure
    C_poincare::FT
end
vars_aux(::AnisoMinDiss, FT) = @vars(Δ::FT)
vars_gradient(::AnisoMinDiss, FT) = @vars(θ_v::FT)
vars_diffusive(::AnisoMinDiss, FT) = @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)
function atmos_init_aux!(
    ::AnisoMinDiss,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end
function gradvariables!(
    m::AnisoMinDiss,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function diffusive!(
    ::AnisoMinDiss,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.∇u = ∇transform.u
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end
function turbulence_tensors(
    m::AnisoMinDiss,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    k̂ = vertical_unit_vector(orientation, param_set, aux)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)

    ∇u = diffusive.turbulence.∇u
    S = symmetrize(∇u)
    normS = strain_rate_magnitude(S)

    δ = aux.turbulence.Δ
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * _inv_Pr_turb, 0, 1))

    δ_vec = SVector(δ, δ, δ)
    δ_m = δ_vec ./ transpose(δ_vec)
    ∇û = ∇u .* δ_m
    Ŝ = symmetrize(∇û)
    ν₀ =
        (m.C_poincare .* δ_vec) .^ 2 * max(
            FT(1e-5),
            -dot(transpose(∇û) * (∇û), Ŝ) / (dot(∇û, ∇û) .+ eps(normS)),
        )

    ν_v = k̂ .* dot(ν₀, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
module NFEdocs

# This package is to host NeuralFieldEq docs

end


text:
@testset "nodahilbert.jl" begin
    @test N(2) ≈ [0.0                0.3183098861837907;
                 -0.3183098861837907 0.0               ]
    @test N(3) ≈ [0.0                  0.3183098861837907 0.15915494309189535;
                 -0.3183098861837907   0.0                0.3183098861837907 ;
                 -0.15915494309189535 -0.3183098861837907 0.0                ]
    @test N(4) ≈ [0    1    1/2 1/3;
                 -1    0    1   1/2;
                 -1/2 -1    0   1  ;
                 -1/3 -1/2 -1   0  ] / π

    function f(i, j)
        if i == j
            return 0.0
        end
        return 1 / (π * (j - i))
    end
    for n in 5:10
        @test N(n) ≈ [f(i, j) for i in 1:n, j in 1:n]
    end

    @test N * [1, 2] ≈ [0.6366197723675814, -0.3183098861837907]
    @test N * [1, 2, 3] ≈ [1.1140846016432673, 0.6366197723675813, -0.7957747154594768]
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

import Logging: min_enabled_level, shouldlog, handle_message
using Base.CoreLogging, Logging, Printf

struct TimeLogger <: AbstractLogger
    initial_time::Float64
    file::Union{Nothing,IOStream}
    screen_log_level::Any
    io_log_level::Any
end

function TimeLogger(;
    initial_time::Float64,
    file::Union{Nothing,IOStream} = nothing,
    screen_log_level = CoreLogging.Info,
    io_log_level = CoreLogging.Info,
)::TimeLogger
    return TimeLogger(initial_time, file, screen_log_level, io_log_level)
end

min_enabled_level(logger::TimeLogger) = logger.io_log_level
shouldlog(logger::TimeLogger, level, _module, group, id) = true

function handle_message(
    logger::TimeLogger,
    level,
    message,
    _module,
    group,
    id,
    filepath,
    line;
    kwargs...,
)
    elapsed_time = time() - logger.initial_time
    time_string = @sprintf("[%12.3f] ", elapsed_time)

    if level >= Logging.Error
        color = :light_red
    elseif level >= Logging.Warn
        color = :light_yellow
    else
        color = :light_green
    end

    flush(stdout)
    flush(stderr)
    Base.Libc.flush_cstdio()
    if level >= logger.screen_log_level
        printstyled(time_string, color = color)
        println(message)
    end
    if logger.file !== nothing && level >= logger.io_log_level
        write(logger.file, time_string)
        write(logger.file, message)
        write(logger.file, "\n")
        flush(logger.file)
    end
    flush(stdout)
    flush(stderr)
    Base.Libc.flush_cstdio()
end

function setup_logger()
    initial_time = time()
    global_logger(TimeLogger(initial_time = initial_time))
    miplearn = pyimport("miplearn")
    miplearn.setup_logger(initial_time)
end

export TimeLogger


text:
using CoordinateTransformations: Translation, LinearMap
using StaticArrays
using Rotations: RotZ, RotX, MRP
# using Binary_BB_Trees
using LinearAlgebra  # TODO: is this used?
using RigidBodyDynamics
using MeshCat
using ColorTypes: RGBA, RGB
using MeshCatMechanisms
using PressureFieldContact
using FileIO

set_zero_subnormals(true)  # optional tweak to make simulation run faster
BLAS.set_num_threads(1)  # optional tweak to make simulation run faster (Disable if using MKL)
spoon_path = joinpath(dirname(pathof(PressureFieldContact)), "..", "test", "data", "spoon.obj")

### Create mechanism and temporary structure
my_mechanism = Mechanism(RigidBody{Float64}("world"); gravity=SVector{3,Float64}(0.0, 0.0, -9.8054))  # create empty mechanism
ts = TempContactStruct(my_mechanism)

### Properties + Names
name_lo = "box_lo"
name_up = "box_up"
name_spoon = "spoon"
rad_box = 0.02
i_prop_box = InertiaProperties(rho=500.0)
c_prop_box = ContactProperties(Ē=1.0e6, d=rad_box)

### Lower Box
eM_box_compliant = eMesh_box(rad_box, SVector{3,Float64}(0.0, 0.0, -rad_box))
add_contact!(ts, name_lo, eM_box_compliant, c_prop_box)

### Upper Box
joint_up = Prismatic(SVector{3,Float64}(0.0, 0.0, 1.0))
body_up, joint_up = add_body_contact!(ts, name_up, eM_box_compliant, c_prop_box, i_prop_box, joint_type=joint_up)

### Spoon
mesh_spoon = load(spoon_path)
scale_HomogenousMesh!(mesh_spoon, 0.01)  # make 0.01 of origional size or spoon will be HUGE
eMesh_spoon = eMesh(mesh_spoon, nothing, nothing)
# assume that the spoon in a shell of this density and thickness (d) for the purposes of inertia tensor calculation
i_prop_spoon = InertiaProperties(rho=900.0, d=0.0025)
_, spoon_joint = add_body_contact!(ts, "spoon", eMesh_spoon, nothing, i_prop_spoon)  # adds body and calculates inertia from the surface mesh

### Friction
add_friction_bristle!(ts, name_spoon, name_lo, μ=0.2, χ=0.2)
add_friction_bristle!(ts, name_spoon, name_up, μ=0.2, χ=0.2)
mech_scen = MechanismScenario(ts, calcXd!, n_quad_rule=1)  # there will be many intersections; first order is fine

### Set initial condition
# NOTE: All states are zero initialized including deformation state variables
set_state_spq!(mech_scen, spoon_joint, rot=RotX(pi/2), trans=SVector(0.0, 0.0, 0.000001))
set_configuration!(mech_scen.float.state, joint_up, [0.10])
x = get_state(mech_scen)

### Run forward dynamics
t_final = 1.6
rr = Radau_for_MechanismScenario(mech_scen)
data_time, data_state = integrate_scenario_radau(rr, mech_scen, x*1, t_final=t_final)
println("Finished compiling and running simulation beginning visualization")

### Add meshcat visualizer
if !@isdefined vis
    vis = Visualizer()
    open(vis)
end

### Add meshes to visualizer
mvis = MechanismVisualizer(my_mechanism, vis)
color_pad = RGBA{Float32}(1.0, 0.792, 0.588, 1.0)
set_mesh_visual!(mvis, mech_scen, name_lo, color_pad)
set_body_mesh_visual!(mvis, mech_scen, name_spoon, [0.5, 0.5, 0.5, 1.0])
set_body_mesh_visual!(mvis, mech_scen, name_up, color_pad)

### Move camera
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 27)
settransform!(vis["/Cameras/default"], Translation(0.0, 0.0, 0.04) ∘ LinearMap(RotZ(-π * 0.35)))

### Playback data
play_recorded_data(mvis, mech_scen, data_time, data_state, slowdown=1.0)


text:
module ParticleInCell

using Reexport
@reexport using GEMPIC
using LinearAlgebra
using Random
using SparseArrays
import Sobol

include("particle_1d1v.jl")
include("landau_damping.jl")
include("poisson_1d.jl")
include("poisson_2d.jl")
include("particle_mesh.jl")
include("particles.jl")
include("fdtd.jl")
include("pic.jl")

end


text:
using Minpack: hybrd
using Base.Test

@testset "analytic" begin

function residual(xvec, arg1, arg2)

    fvec = zeros(xvec)

    fvec[1] = xvec[1] - arg1
    fvec[2] = xvec[1] + xvec[2] - arg2

    return fvec
end

x0 = [1.0, 0.0]
args = (3.0, 1.0)
tol = 1e-8

x, f = hybrd(residual, x0, args, tol)

@test isapprox(x[1], 3.0; atol=1e-8)
@test isapprox(x[2], -2.0; atol=1e-8)
@test isapprox(f[1], 0.0; atol=1e-8)
@test isapprox(f[2], 0.0; atol=1e-8)

function residual2(xvec)

    fvec = zeros(xvec)

    fvec[1] = exp(-exp(-(xvec[1]+xvec[2]))) - xvec[2]*(1+xvec[1]^2)
    fvec[2] = xvec[1]*cos(xvec[2]) + xvec[2]*sin(xvec[1]) - 0.5

    return fvec
end

x0 = [0.0, 0.0]
args = ()
tol = 1e-8

x, f = hybrd(residual2, x0, args, tol)

@test isapprox(x[1], 0.3532; atol=1e-4)
@test isapprox(x[2], 0.6061; atol=1e-4)
@test isapprox(f[1], 0.0; atol=1e-8)
@test isapprox(f[2], 0.0; atol=1e-8)

end  # test set


text:
using Test
using Revise
using LinearAlgebra
using MFGnet
using Printf

function initializeWeights(d,m,ar=x->x)
    (K1,b1) = (ar(.01*randn(m,d+1)),ar(0.1*randn(m)))
    Θ1 = (K1,b1)
    (w0,A0,c0,z0) = (ar(ones(m)/sqrt(m)),ar(zeros(d+1,d+1)),ar(zeros(d+1)),ar([1.0]))
    return (w0,Θ1,A0,c0,z0)
end

Rs = [Float64,Float32]
for k = 1:length(Rs)
    R = Rs[k]

    d = 2
    m = 10
    nex = 20;

    N = SingleLayer()
    ar(x) = R.(x)
    (w,Θ,A,c,z0) = initializeWeights(d,m,ar)

    s = randn(R,d+1,nex)

    z = N(s,Θ)
    @test eltype(z) == R
    Φ = w'*z
    ds = randn(R,size(s))
    dΦ,d2Φ = getGradAndHessian(N,w,s,Θ)
    @test eltype(dΦ) == R
    @test eltype(d2Φ) == R
    dΦds = sum(ds.*dΦ,dims=1)
    curv = zeros(R,1,nex)
    for j=1:nex
        curv[1,j] = ds[:,j]'*d2Φ[:,:,j]*ds[:,j]
    end
    curvt = MFGnet.getHessmv(N,w,ds,s,Θ)
    trH = MFGnet.getTraceHess(N,w,s,Θ)
    trHt, = MFGnet.getTraceHessAndGrad(N,w,s,Θ)
    @test norm(trH-trHt)/norm(trH) < sqrt(eps(R))
    for j=1:nex
        @test abs(trH[j] - tr(d2Φ[:,:,j]))/abs(tr(d2Φ[:,:,j])) < sqrt(eps(R))
    end

    for j=1:15
        h = R(2.0^(-j))
        Φt = w'*N(s + h .* ds,Θ)

        E0 = norm(Φ-Φt)
        E1 = norm(Φ+h*dΦds - Φt)
        E2 = norm(Φ+h*dΦds+0.5*h^2*curv - Φt)

        @printf("h=%1.3e\t\tE0=%1.3e\tE1=%1.3e\tE1=%1.3e\n",h,E0,E1,E2)
    end
    @test true

    @testset "adjoint and mult. dispatch test  $(Rs[k])" begin
    d2 = 2
    V = randn(R,d+1,d2,nex)
    W = randn(R,m,d2,nex)
    S = randn(R,d+1,nex)
    JTW = MFGnet.getJSTmv(N,W,S,Θ)
    @test eltype(JTW) == R
    JV = MFGnet.getJSmv(N,V,S,Θ)
    @test eltype(JV) == R

    sumJS = R(0.0)
    sumJST = R(0.0)
    for j=1:d2
        # global sumJS, sumJST
        MFGnet.getJSmv(N,V[:,j,:],S,Θ)
        sumJS   = sumJS + norm(JV[:,j,:] - MFGnet.getJSmv(N,V[:,j,:],S,Θ))
        sumJST  = sumJST +  norm(JTW[:,j,:] - MFGnet.getJSTmv(N,W[:,j,:],S,Θ))
    end
    @test sumJS <= eps(R)
    @test sumJST <= eps(R)

    # adjoint test
    vJTw = sum(V.*JTW,dims=1)
    wJv  = sum(W.*JV, dims=1)
    @test sum(vJTw - wJv) <= sqrt(eps(R))

    # test JSJSTmv
    S = randn(R,d+1,nex)
    dZ = randn(R,m,1)
    d2Z = randn(R,m,m,1)

    H1 = MFGnet.getJSJSTmv(N,dZ,S,Θ)
    @test eltype(H1)==R
    H2 = MFGnet.getJSJSTmv(N,dZ,zeros(R,m,m),S,Θ)
    @test eltype(H2)==R
    @test norm(vec(H1) - vec(H2)) <= sqrt(eps(R))
    end

end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using Test
using SolidStateDetectors
using StaticArrays

import SolidStateDetectors.ConstructiveSolidGeometry as CSG
import SolidStateDetectors.ConstructiveSolidGeometry: Geometry, Dictionary

T = Float64

default_units = SolidStateDetectors.default_unit_tuple()
no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(CartesianVector{T}))

@testset "Test primitive parsing" begin
    @testset "Cone" begin
        for bot in (2.0, Dict("from" => 1.0, "to" => 2.0)),
            top in (2.0, 1.0, Dict("from" => 1.0, "to" => 2.0), Dict("from" => 2.0, "to" => 4.0)),
            φ in ((30,180),(0,360))

            dict = Dict("difference" => [
                Dict("cone"   => Dict(
                    "r"       => Dict("bottom" => bot, "top" => top),
                    "phi"     => Dict("from" => "$(φ[1])°", "to" => "$(φ[2])°"),
                    "h"       => 1.0))
                for i in 1:2
            ])

            c = Geometry(T, dict, default_units, no_translations)

            # Conversion from Geometry -> Dict and Dict -> Geometry should result in the same geometry
            output = Dictionary(c)
            name = collect(keys(output["difference"][1]))[1] # "tube" or "cone"
            if haskey(output["difference"][1][name], "phi")
                @test CSG._parse_value(T, output["difference"][1][name]["phi"]["to"], default_units.angle) == CSG._parse_value(T, "$(φ[2])°", default_units.angle)
            end
            # Internally, the primitive stores the offset of phi in a rotation, but it should not be part of the config file
            @test !haskey(output["difference"][1], "rotation")
            @test c.a.rotation ≈ Geometry(T, output, default_units, no_translations).a.rotation

            # No warnings or errors when decomposing the Cones into surfaces
            @test_nowarn CSG.surfaces(c.a)
            @test_nowarn CSG.surfaces(c.b)

            # Check if all Cones are saved the right way
            @test c.a isa CSG.Cone
            @test c.b isa CSG.Cone
        end
    end
    @testset "Torus" begin
        for r_tube in (2.0, Dict("from" => 1.0, "to" => 2.0)),
            φ in ((30,180),(0,360)),
            θmin in 0:45:90,
            θmax in 135:45:270

            dict = Dict("difference" => [
                Dict( "torus" => Dict(
                    "r_torus" => 10.0,
                    "r_tube"  => r_tube,
                    "phi"     => Dict("from" => "$(φ[1])°", "to" => "$(φ[2])°"),
                    "theta"   => Dict("from" => "$(θmin)°", "to" => "$(θmax)°")))
                for i in 1:2
            ])
            t = Geometry(T, dict, default_units, no_translations)

            # Conversion from Geometry -> Dict and Dict -> Geometry should result in the same geometry
            output = Dictionary(t)
            if haskey(output["difference"][1]["torus"], "phi")
                @test CSG._parse_value(T, output["difference"][1]["torus"]["phi"]["to"], default_units.angle) == CSG._parse_value(T, "$(φ[2])°", default_units.angle)
            end
            # Internally, the primitive stores the offset of phi in a rotation, but it should not be part of the config file
            @test !haskey(output["difference"][1], "rotation")
            @test t.a.rotation ≈ Geometry(T, output, default_units, no_translations).a.rotation

            # No warnings or errors when decomposing the Torus into surfaces
            @test_nowarn CSG.surfaces(t.a)
            @test_nowarn CSG.surfaces(t.b)

            # Check if all Torus are saved the right way
            @test t.a isa CSG.Torus
            @test t.b isa CSG.Torus
        end
    end
end


text:
function solve(data)
    volume(v) = prod(max(0, v[i+1] - v[i] + 1) for i=1:2:5)
    intersection(v1, v2) = [(i % 2 == 1 ? max : min)(v1[i], v2[i]) for i=1:6]
    cubes = []
    for (onoff, nums) in data
        function remove_volume(original, remove)
            # Function to remove one cube from another, by replacing it with up to 6 other cubes
            intersecting = intersection(original, remove)
            volume(intersecting) == 0 && return [original]
            original_divided = []
            add(v) = volume(v) != 0 && push!(original_divided, v)
            x, X, y, Y, z, Z = original
            xi, Xi, yi, Yi, zi, Zi = intersecting
            add([x, xi-1, y, Y, z, Z])
            add([Xi+1, X, y, Y, z, Z])
            add([xi, Xi, y, yi-1, z, Z])
            add([xi, Xi, Yi+1, Y, z, Z])
            add([xi, Xi, yi, Yi, z, zi-1])
            add([xi, Xi, yi, Yi, Zi+1, Z])
            return original_divided
        end
        # Remove the current cube from every previous cube (doesn't matter if on or off)
        !isempty(cubes) && (cubes = reduce(vcat, remove_volume(cube, nums) for cube in cubes))
        # Add the cube if it is turned on
        onoff && push!(cubes, nums)
    end
    # Sum all cubes, since none are intersecting after the previous operations
    return cubes .|> volume |> sum
end
function main()
    extract_ints(str) = parse.(Int, m.match for m in eachmatch(r"-?\d+", str))
    data = split.(readlines()) .|> sp -> (sp[1] == "on", extract_ints(sp[2]))
    data_part1 = filter(nums -> all(i ∈ -50:50 for i=nums[2]), data)

    println(solve(data_part1))
    println(solve(data))
end
main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
include("nodes.jl")
include("resolver.jl")

immutable ComposerError
    context::(@compat Union{AbstractString, Void})
    context_mark::(@compat Union{Mark, Void})
    problem::(@compat Union{AbstractString, Void})
    problem_mark::(@compat Union{Mark, Void})
    note::(@compat Union{AbstractString, Void})

    function ComposerError(context=nothing, context_mark=nothing,
                           problem=nothing, problem_mark=nothing,
                           note=nothing)
        new(context, context_mark, problem, problem_mark, note)
    end
end

type Composer
    input::EventStream
    anchors::Dict{AbstractString, Node}
    resolver::Resolver
end

function compose(events)
    composer = Composer(events, Dict{AbstractString, Node}(), Resolver())
    @assert typeof(forward!(composer.input)) == StreamStartEvent
    node = compose_document(composer)
    if typeof(peek(composer.input)) == StreamEndEvent
        forward!(composer.input)
    else
        @assert typeof(peek(composer.input)) == DocumentStartEvent
    end
    node
end

function compose_document(composer::Composer)
    @assert typeof(forward!(composer.input)) == DocumentStartEvent
    node = compose_node(composer, nothing, nothing)
    @assert typeof(forward!(composer.input)) == DocumentEndEvent
    empty!(composer.anchors)
    node
end

function compose_node(composer::Composer, parent::(@compat Union{Node, Void}),
                      index::(@compat Union{Int, Node, Void}))
    event = peek(composer.input)
    if typeof(event) == AliasEvent
        forward!(composer.input)
        anchor = event.anchor
        if !haskey(composer.anchors, anchor)
            throw(ComposerError(nothing, nothing, "found undefined alias $(anchor)",
                                event.start_mark))
        end
        return composer.anchors[anchor]
    end

    anchor = event.anchor
    if !is(anchor, nothing)
        if haskey(composer.anchors, anchor)
            throw(ComposerError(
                "found duplicate anchor $(anchor); first occurance",
                composer.anchors[anchor].start_mark, "second occurence",
                event.start_mark))
        end
    end

    node = nothing
    if typeof(event) == ScalarEvent
        node = compose_scalar_node(composer, anchor)
    elseif typeof(event) == SequenceStartEvent
        node = compose_sequence_node(composer, anchor)
    elseif typeof(event) == MappingStartEvent
        node = compose_mapping_node(composer, anchor)
    end

    node
end

function compose_scalar_node(composer::Composer, anchor::(@compat Union{AbstractString, Void}))
    event = forward!(composer.input)
    tag = event.tag
    if tag === nothing || tag == "!"
        tag = resolve(composer.resolver, ScalarNode,
                      event.value, event.implicit)
    end

    node = ScalarNode(tag, event.value, event.start_mark, event.end_mark,
                      event.style)
    if !is(anchor, nothing)
        composer.anchors[anchor] = node
    end

    node
end

function compose_sequence_node(composer::Composer, anchor::(@compat Union{AbstractString, Void}))
    start_event = forward!(composer.input)
    tag = start_event.tag
    if tag === nothing || tag == "!"
        tag = resolve(composer.resolver, SequenceNode,
                      nothing, start_event.implicit)
    end

    node = SequenceNode(tag, Any[], start_event.start_mark, nothing,
                        start_event.flow_style)
    if !is(anchor, nothing)
        composer.anchors[anchor] = node
    end

    index = 1
    while typeof(peek(composer.input)) != SequenceEndEvent
        push!(node.value, compose_node(composer, node, index))
        index += 1
    end

    end_event = forward!(composer.input)
    node.end_mark = end_event.end_mark

    node
end

function compose_mapping_node(composer::Composer, anchor::(@compat Union{AbstractString, Void}))
    start_event = forward!(composer.input)
    tag = start_event.tag
    if tag === nothing || tag == "!"
        tag = resolve(composer.resolver, MappingNode,
                      nothing, start_event.implicit)
    end

    node = MappingNode(tag, Any[], start_event.start_mark, nothing,
                       start_event.flow_style)
    if !is(anchor, nothing)
        composer.anchors[anchor] = node
    end

    while typeof(peek(composer.input)) != MappingEndEvent
        item_key = compose_node(composer, node, nothing)
        item_value = compose_node(composer, node, item_key)
        push!(node.value, (item_key, item_value))
    end

    end_event = forward!(composer.input)
    node.end_mark = end_event.end_mark

    node
end


text:
using Documenter, Mimi

makedocs(
    modules = [Mimi],
	sitename = "Mimi.jl",
	pages = [
		"Home" => "index.md",
		"Installation Guide" => "installation.md",
        "User Guide" => "userguide.md",
		"Tutorial" => "tutorial.md",
		"FAQ" => "faq.md",
		"Reference" => "reference.md",
		"Integration Guide: Port to v0.5.0" => "integrationguide.md"]
)

deploydocs(
    repo = "github.com/anthofflab/Mimi.jl.git",
)


text:
function MALA(likelihood_and_derivative_calculator,
              protein_at_observations,
              measurement_variance,
              number_of_samples,
              initial_position,
              step_size,
              proposal_covariance=I,
              thinning_rate=1,
              known_parameter_dict=Dict())

    mean_protein = mean(protein_at_observations[:,2])
    # initialise the covariance proposal matrix
    number_of_parameters = length(initial_position) - length(keys(known_parameter_dict))
    known_parameter_indices = [Int(known_parameter_dict[i][1]) for i in keys(known_parameter_dict)]
    known_parameter_values = [known_parameter_dict[i][2] for i in keys(known_parameter_dict)]
    unknown_parameter_indices = [i for i in 1:length(initial_position) if i ∉ known_parameter_indices]

    # check if default value is used, and set to q x q identity
    if proposal_covariance == I
        identity = true
    else
        identity = false
        proposal_cholesky = cholesky(proposal_covariance + 1e-8*I).L
    end

    proposal_covariance_inverse = inv(proposal_covariance)

    # initialise samples matrix and acceptance ratio counter
    accepted_moves = 0
    mcmc_samples = zeros(number_of_samples,number_of_parameters)
    mcmc_samples[1,:] = initial_position[unknown_parameter_indices]
    number_of_iterations = number_of_samples*thinning_rate

    # set LAP parameters
    k = 1
    c0 = 1.0
    c1 = log(10)/log(number_of_samples/5)

    # initial markov chain
    current_position = copy(initial_position)
    current_log_likelihood, current_log_likelihood_gradient = likelihood_and_derivative_calculator(protein_at_observations,
                                                                                current_position,
                                                                                mean_protein,
                                                                                measurement_variance)
    for iteration_index in 2:number_of_iterations
        # progress measure
        if iteration_index%(number_of_iterations//10)==0
            println(string("Progress: ",100*(iteration_index/number_of_iterations),'%'))
        end #if

        proposal = zeros(length(initial_position))
        if identity
            proposal[unknown_parameter_indices] .= ( current_position[unknown_parameter_indices] .+
                                                    step_size.*current_log_likelihood_gradient[unknown_parameter_indices]./2 +
                                                    sqrt(step_size)*randn(number_of_parameters) )
        else
            proposal[unknown_parameter_indices] .= ( current_position[unknown_parameter_indices] .+
                                                    step_size.*proposal_covariance*(current_log_likelihood_gradient[unknown_parameter_indices]./2) .+
                                                    sqrt(step_size).*proposal_cholesky*randn(number_of_parameters) )
        end

        # compute transition probabilities for acceptance step
        # fix known parameters
        # import pdb; pdb.set_trace()
        if length(known_parameter_dict) > 0
            proposal[known_parameter_indices] .= copy(known_parameter_values)
        end #if

        proposal_log_likelihood, proposal_log_likelihood_gradient = likelihood_and_derivative_calculator(protein_at_observations,
                                                                                proposal,
                                                                                mean_protein,
                                                                                measurement_variance)

        # if any of the parameters were negative we get -inf for the log likelihood
        if proposal_log_likelihood == -Inf
            if iteration_index%thinning_rate == 0
                mcmc_samples[Int64(iteration_index/thinning_rate),:] .= current_position[unknown_parameter_indices]
            end

            # LAP stuff also needed here
            acceptance_probability = 0
            if iteration_index%k == 0 && iteration_index > 1
                gamma_1 = 1/(iteration_index^c1)
                gamma_2 = c0*gamma_1
                log_step_size_squared = log(step_size^2) + gamma_2*(acceptance_probability - 0.574)
                step_size = sqrt(exp(log_step_size_squared))
            end
            continue
        end

        forward_helper_variable = ( proposal[unknown_parameter_indices] .- current_position[unknown_parameter_indices] .-
                                    step_size.*proposal_covariance*(current_log_likelihood_gradient[unknown_parameter_indices]./2) )
        backward_helper_variable = ( current_position[unknown_parameter_indices] .- proposal[unknown_parameter_indices] .-
                                     step_size.*proposal_covariance*(proposal_log_likelihood_gradient[unknown_parameter_indices]./2) )
        transition_kernel_pdf_forward = (-transpose(forward_helper_variable)*proposal_covariance_inverse*forward_helper_variable) /(2*step_size)
        transition_kernel_pdf_backward = (-transpose(backward_helper_variable)*proposal_covariance_inverse*backward_helper_variable)/(2*step_size)

        # accept-reject step
        acceptance_probability = min(1,exp(proposal_log_likelihood - transition_kernel_pdf_forward - current_log_likelihood + transition_kernel_pdf_backward))
        # print(acceptance_probability)
        if rand() < acceptance_probability
            current_position .= proposal
            current_log_likelihood = proposal_log_likelihood
            current_log_likelihood_gradient .= proposal_log_likelihood_gradient
            accepted_moves += 1
        end

        if iteration_index%thinning_rate == 0
            mcmc_samples[Int64(iteration_index/thinning_rate),:] .= current_position[unknown_parameter_indices]
        end

        # LAP stuff
        if iteration_index%k == 0 && iteration_index > 1
            gamma_1 = 1/(iteration_index^c1)
            gamma_2 = c0*gamma_1
            log_step_size_squared = log(step_size^2) + gamma_2*(acceptance_probability - 0.574)
            step_size = sqrt(exp(log_step_size_squared))
        end
    end # for
    println(string("Acceptance ratio: ",accepted_moves/number_of_iterations))

    return mcmc_samples
end

"""
A function which gives a (hopefully) decent MALA output for a given dataset with known or
unknown parameters. If a previous output already exists, this will be used to create a
proposal covariance matrix, otherwise one will be constructed with a two step warm-up
process.
"""
function run_mala_for_dataset(data_filename,
                              protein_at_observations,
                              measurement_variance,
                              number_of_parameters,
                              known_parameter_dict,
                              step_size = 1,
                              number_of_chains = 8,
                              number_of_samples = 80000)
    # make sure all data starts from time "zero"
    for i in 1:size(protein_at_observations,1)
        protein_at_observations[i,1] -= protein_at_observations[1,1]
    end

    mean_protein = mean(protein_at_observations[:,2])

    loading_path = joinpath(dirname(pathof(DelayedKalmanFilter)),"../test/data")
    saving_path = joinpath(dirname(pathof(DelayedKalmanFilter)),"../test/output")

    # if we already have mcmc samples, we can use them to construct a covariance matrix to directly sample
    if ispath(joinpath(saving_path,string("final_parallel_mala_output_",data_filename)))
        println("Posterior samples already exist, sampling directly without warm up...")

        mala_output = load(joinpath(saving_path,string("final_parallel_mala_output_",data_filename)),"mcmc_samples")
        previous_number_of_chains = size(mala_output,1)
        previous_number_of_samples = size(mala_output,2)
        previous_number_of_parameters = size(mala_output,3)

        # construct proposal covariance matrix
        new_number_of_samples = previous_number_of_samples - Int64(floor(previous_number_of_samples/2))
        burn_in = Int64(floor(previous_number_of_samples/2))+1

        samples_with_burn_in = reshape(mala_output[:,burn_in:end,:],
                                       (new_number_of_samples*previous_number_of_chains,previous_number_of_parameters))
        proposal_covariance = cov(samples_with_burn_in)

        # start from mean
        states = zeros((number_of_chains,7))
        states[:,[3,4]] = [log(log(2)/30),log(log(2)/90)]
        states[:,[1,2,5,6,7]] = mean(samples_with_burn_in,dims=1)

        # turn into array of arrays
        initial_states = [states[i,:] for i in 1:size(states,1)]

        # pool_of_processes = mp_pool.ThreadPool(processes = number_of_chains)
        # process_results = [ pool_of_processes.apply_async(MALA,
        #                                                   args=(protein_at_observations,
        #                                                         measurement_variance,
        #                                                         number_of_samples,
        #                                                         initial_state,
        #                                                         step_size,
        #                                                         proposal_covariance,
        #                                                         1,
        #                                                         known_parameters))
        #                     for initial_state in initial_states ]
        # ## Let the pool know that these are all so that the pool will exit afterwards
        # # this is necessary to prevent memory overflows.
        # pool_of_processes.close()
        #
        # array_of_chains = zeros(number_of_chains,number_of_samples,number_of_parameters)
        # for chain_index, process_result in enumerate(process_results):
        #     this_chain = process_result.get()
        #     array_of_chains[chain_index,:,:] = this_chain
        # pool_of_processes.join()
        array_of_chains = reshape(MALA(log_likelihood_and_derivative_with_prior_and_transformation,
                                       protein_at_observations,
                                       measurement_variance,
                                       number_of_samples,
                                       initial_states[1],
                                       step_size,
                                       proposal_covariance,
                                       1,
                                       known_parameter_dict),
                                  (1,number_of_samples,number_of_parameters))

        save(joinpath(saving_path,string("final_parallel_mala_output_",data_filename)),"mcmc_samples",array_of_chains)

    else
        # warm up chain
        println(string("New data set, initial warm up with ",string(Int64(floor(number_of_samples*0.3)))," samples..."))
        # Initialise by Latin Hypercube sampling
        println("Latin Hypercube sampling initial positions...")
        plan, _ = LHCoptim(20,number_of_parameters,1000)
        scaled_plan = scaleLHC(plan,[(minimum(protein_at_observations[:,2])/2,2*mean(protein_at_observations[:,2])),
                                     (2.5,5.5),
                                     (log(0.1),log(100.)),
                                     (log(0.1),log(35.)),
                                     (5.,35.)])
        # sample subset - one for each chain
        sampling_indices = sample(1:20,number_of_chains;replace=false)
        scaled_plan = scaled_plan[sampling_indices,:]

        states = zeros(number_of_chains,7)
        if number_of_chains > 1
            states[:,[3,4]] .= [log(log(2)/30),log(log(2)/90)]
        else
            states[:,[3,4]] = [log(log(2)/30),log(log(2)/90)]
        end
        states[:,[1,2,5,6,7]] .= scaled_plan
        # turn into array of arrays
        initial_states = [states[i,:] for i in 1:size(states,1)]

        println(string("Warming up with ",string(Int64(0.3*number_of_samples))," samples..."))
        initial_burnin_number_of_samples = Int64(0.3*number_of_samples)

        # pool_of_processes = mp_pool.ThreadPool(processes = number_of_chains)
        # process_results = [ pool_of_processes.apply_async(hes_inference.kalman_mala,
        #                                                   args=(protein_at_observations,
        #                                                         measurement_variance,
        #                                                         initial_burnin_number_of_samples,
        #                                                         initial_state,
        #                                                         step_size,
        #                                                         np.power(np.diag([2*mean_protein,4,9,8,39]),2),# initial variances are width of prior squared
        #                                                         1, # thinning rate
        #                                                         known_parameters))
        #                     for initial_state in initial_states ]
        # ## Let the pool know that these are all so that the pool will exit afterwards
        # # this is necessary to prevent memory overflows.
        # pool_of_processes.close()
        #
        # array_of_chains = zeros(number_of_chains,initial_burnin_number_of_samples,number_of_parameters)
        # for chain_index, process_result in enumerate(process_results):
        #     this_chain = process_result.get()
        #     array_of_chains[chain_index,:,:] = this_chain
        # pool_of_processes.join()

        proposal_covariance = diagm([2*mean_protein,4,9,8,39]).^2

        array_of_chains = reshape(MALA(log_likelihood_and_derivative_with_prior_and_transformation,
                                       protein_at_observations,
                                       measurement_variance,
                                       initial_burnin_number_of_samples,
                                       initial_states[1],
                                       step_size,
                                       proposal_covariance,
                                       1,
                                       known_parameter_dict),
                                  (1,initial_burnin_number_of_samples,number_of_parameters))

        save(joinpath(saving_path,string("first_parallel_mala_output_",data_filename)),"mcmc_samples",array_of_chains)

        println(string("Second warm up with ",string(Int64(number_of_samples*0.7))," samples..."))
        second_burnin_number_of_samples = Int64(0.7*number_of_samples)

        # construct proposal covariance matrix
        new_number_of_samples = initial_burnin_number_of_samples - Int64(floor(initial_burnin_number_of_samples/2))
        burn_in = Int64(floor(initial_burnin_number_of_samples/2))+1
        samples_with_burn_in = reshape(array_of_chains[:,burn_in:end,:],
                                       (new_number_of_samples*number_of_chains,number_of_parameters))
        proposal_covariance = cov(samples_with_burn_in)

        # make new initial states
        println("Latin Hypercube sampling initial positions...")
        plan, _ = LHCoptim(20,number_of_parameters,1000)
        scaled_plan = scaleLHC(plan,[(minimum(protein_at_observations[:,2])/2,2*mean(protein_at_observations[:,2])),
                                     (2.5,5.5),
                                     (log(0.1),log(100.)),
                                     (log(0.1),log(35.)),
                                     (5.,35.)])
        # sample subset - one for each chain
        sampling_indices = sample(1:20,number_of_chains;replace=false)
        scaled_plan = scaled_plan[sampling_indices,:]

        states = zeros(number_of_chains,7)
        if number_of_chains > 1
            states[:,[3,4]] .= [log(log(2)/30),log(log(2)/90)]
        else
            states[:,[3,4]] = [log(log(2)/30),log(log(2)/90)]
        end
        states[:,[1,2,5,6,7]] .= scaled_plan
        # turn into array of arrays
        initial_states = [states[i,:] for i in 1:size(states,1)]

        # pool_of_processes = mp_pool.ThreadPool(processes = number_of_chains)
        # process_results = [ pool_of_processes.apply_async(hes_inference.kalman_mala,
        #                                                   args=(protein_at_observations,
        #                                                         measurement_variance,
        #                                                         second_burnin_number_of_samples,
        #                                                         initial_state,
        #                                                         step_size,
        #                                                         proposal_covariance,
        #                                                         1, # thinning rate
        #                                                         known_parameters))
        #                     for initial_state in initial_states ]
        # ## Let the pool know that these are all finished so that the pool will exit afterwards
        # # this is necessary to prevent memory overflows.
        # pool_of_processes.close()
        #
        # array_of_chains = zeros(number_of_chains,second_burnin_number_of_samples,number_of_parameters)
        # for chain_index, process_result in enumerate(process_results):
        #     this_chain = process_result.get()
        #     array_of_chains[chain_index,:,:] = this_chain
        # pool_of_processes.join()

        array_of_chains = reshape(MALA(log_likelihood_and_derivative_with_prior_and_transformation,
                                       protein_at_observations,
                                       measurement_variance,
                                       second_burnin_number_of_samples,
                                       initial_states[1],
                                       step_size,
                                       proposal_covariance,
                                       1,
                                       known_parameter_dict),
                                  (1,second_burnin_number_of_samples,number_of_parameters))

        save(joinpath(saving_path,string("second_parallel_mala_output_",data_filename)),"mcmc_samples",array_of_chains)

        # sample directly
        println("Now sampling directly...")
        new_number_of_samples = second_burnin_number_of_samples - Int64(floor(second_burnin_number_of_samples/2))
        burn_in = Int64(floor(second_burnin_number_of_samples/2))+1
        samples_with_burn_in = reshape(array_of_chains[:,burn_in:end,:],
                                       (new_number_of_samples*number_of_chains,number_of_parameters))
        proposal_covariance = cov(samples_with_burn_in)

        # make new initial states
        println("Latin Hypercube sampling initial positions...")
        plan, _ = LHCoptim(20,number_of_parameters,1000)
        scaled_plan = scaleLHC(plan,[(minimum(protein_at_observations[:,2])/2,2*mean(protein_at_observations[:,2])),
                                     (2.5,5.5),
                                     (log(0.1),log(100.)),
                                     (log(0.1),log(35.)),
                                     (5.,35.)])
        # sample subset - one for each chain
        sampling_indices = sample(1:20,number_of_chains;replace=false)
        scaled_plan = scaled_plan[sampling_indices,:]

        states = zeros(number_of_chains,7)
        if number_of_chains > 1
            states[:,[3,4]] .= [log(log(2)/30),log(log(2)/90)]
        else
            states[:,[3,4]] = [log(log(2)/30),log(log(2)/90)]
        end
        states[:,[1,2,5,6,7]] .= scaled_plan
        # turn into array of arrays
        initial_states = [states[i,:] for i in 1:size(states,1)]

        # pool_of_processes = mp_pool.ThreadPool(processes = number_of_chains)
        # process_results = [ pool_of_processes.apply_async(hes_inference.kalman_mala,
        #                                                   args=(protein_at_observations,
        #                                                         measurement_variance,
        #                                                         number_of_samples,
        #                                                         initial_state,
        #                                                         step_size,
        #                                                         proposal_covariance,
        #                                                         1, # thinning rate
        #                                                         known_parameters))
        #                     for initial_state in initial_states ]
        # ## Let the pool know that these are all finished so that the pool will exit afterwards
        # # this is necessary to prevent memory overflows.
        # pool_of_processes.close()
        #
        # array_of_chains = np.zeros((number_of_chains,number_of_samples,number_of_parameters))
        # for chain_index, process_result in enumerate(process_results):
        #     this_chain = process_result.get()
        #     array_of_chains[chain_index,:,:] = this_chain
        # pool_of_processes.join()

        array_of_chains = reshape(MALA(log_likelihood_and_derivative_with_prior_and_transformation,
                                       protein_at_observations,
                                       measurement_variance,
                                       number_of_samples,
                                       initial_states[1],
                                       step_size,
                                       proposal_covariance,
                                       1,
                                       known_parameter_dict),
                                  (1,number_of_samples,number_of_parameters))

        save(joinpath(saving_path,string("final_parallel_mala_output_",data_filename)),"mcmc_samples",array_of_chains)
    end #if-else
end


text:
mutable struct Program
    nodes::Vector{Any}
    args::Vector{Vector{Int}}
    return_node::Int
    graph::DirectedGraph
end

struct ArgumentNode end

struct CallNode
    f::Any
    args::Vector{Any}
end

struct TupleNode
    args::Vector{Any}
end

struct ZeroNode
    node_id::Int
end

struct TrackedValue
    id::Int
    program::Program
end

Base.show(io::IO, tv::TrackedValue) = print(io, "TrackedValue($(tv.id))")

function add_call!(program::Program, f, args::Vector{Any})
    push!(program.nodes, CallNode(f, args))

    node_id = create_node!(program.graph)

    for arg in args
        if arg isa TrackedValue
            add_edge!(program.graph, arg.id => node_id)
        end
    end

    return TrackedValue(node_id, program)
end

function add_call!(f, args::Vector{Any})
    program = nothing

    for arg in args
        if arg isa TrackedValue
            program = arg.program
            break
        end
    end

    program::Program

    return add_call!(program, f, args)
end

function add_zero!(program::Program, node_id::Int)
    push!(program.nodes, ZeroNode(node_id))

    zero_node_id = create_node!(program.graph)
    add_edge!(program.graph, node_id => zero_node_id)

    # TODO: this feels a bit weird
    return TrackedValue(zero_node_id, program)
end

function add_tuple!(program::Program, args::Vector{Any})
    push!(program.nodes, TupleNode(args))

    node_id = create_node!(program.graph)

    for arg in args
        if arg isa TrackedValue
            add_edge!(program.graph, arg.id => node_id)
        end
    end

    return TrackedValue(node_id, program)
end

function record_function(@nospecialize(f), nargs::Int)
    program = Program(Any[], Vector{Int}[], 0, DirectedGraph())

    args = TrackedValue[]
    arg_ids = Int[]
    push!(program.args, arg_ids)

    for i in 1:nargs
        push!(args, TrackedValue(i, program))

        push!(program.nodes, ArgumentNode())
        node_id = create_node!(program.graph)
        push!(arg_ids, node_id)
    end

    return_value = f(args...)

    if return_value isa NTuple{N,TrackedValue} where N
        return_value = add_tuple!(program, Any[return_value...])
    end

    if !(return_value isa TrackedValue)
        error("unsupported return value") # TODO: better error message
    end

    program.return_node = return_value.id

    return program
end

function compact!(program)
    required_nodes = Set(reduce(vcat, program.args))
    union!(required_nodes, ancestors(program.graph, Set(program.return_node)))

    remove_nodes!(program.graph, setdiff(nodes(program.graph), required_nodes))

    return program
end

function create_subprograms(program::Program)
    subprograms = Vector{Int}[]

    compact!(program)

    executable_order = topological_sort(program.graph)

    unassigned_nodes = copy(nodes(program.graph))

    # TODO: this seems confusing
    for args in reverse(program.args)
        subprogram = copy(args)
        setdiff!(unassigned_nodes, subprogram)

        nodes = descendants(program.graph, Set(args))
        intersect!(nodes, unassigned_nodes)
        setdiff!(unassigned_nodes, nodes)

        append!(subprogram, intersect(executable_order, nodes))
        push!(subprograms, subprogram)
    end

    @assert isempty(unassigned_nodes)

    return reverse(subprograms)
end


text:
# Interval
# ========
#
# Base interval types and utilities.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Note, just to be clear: this shadows IntervalTrees.Interval
"""
    struct Interval{T} <: IntervalTrees.AbstractInterval{Int64}

The first three fields (`seqname`, `first`, and `last`) are mandatory arguments when constructing the [`Interval`](@ref Interval) object.

# Fields
- `seqname::String`: the sequence name associated with the interval.
- `first::Int64`: the leftmost position.
- `last::Int64`: the rightmost position.
- `strand::Strand`: the [`strand`](@ref Strand).
- `metadata::T`
"""
struct Interval{T} <: IntervalTrees.AbstractInterval{Int64}
    seqname::String
    first::Int64
    last::Int64
    strand::Strand
    metadata::T
end

function Interval(seqname::AbstractString, first::Integer, last::Integer, strand::Union{Strand,Char}=STRAND_BOTH, metadata::T=nothing) where T
    return Interval{T}(seqname, first, last, strand, metadata)
end

function Interval(seqname::AbstractString, range::UnitRange{R}, strand::Union{Strand,Char}=STRAND_BOTH, metadata::T=nothing) where {T,R<:Integer}
    return Interval{T}(seqname, first(range), last(range), strand, metadata)
end

"""
    Interval{T}(data)

The returned data is converted to Interval{T} if there is an implemented [`Base.convert`](https://docs.julialang.org/en/v1/base/base/#Base.convert) function for the type of data.
This method provides a useful hook for converting custom types to Interval{T}.
"""
function Interval{T}(data) :: Interval{T} where T
    return data #Note: the returned data is converted to Interval{T}.
end

function BioGenerics.seqname(i::Interval)
    return i.seqname
end

function BioGenerics.metadata(i::Interval)
    return i.metadata
end

function strand(i::Interval)
    return i.strand
end

"""
    leftposition(i::Interval)

Return the leftmost position of `i`.
"""
function BioGenerics.leftposition(i::Interval)
    return i.first
end

"""
    rightposition(i::Interval)

Return the rightmost position of `i`.
"""
function BioGenerics.rightposition(i::Interval)
    return i.last
end

IntervalTrees.first(i::Interval) = leftposition(i)
IntervalTrees.last(i::Interval) = rightposition(i)

function Base.isless(a::Interval, b::Interval, seqname_isless::Function=isless)
    if seqname(a) != seqname(b)
        return seqname_isless(seqname(a), seqname(b))::Bool
    end

    if leftposition(a) != leftposition(b)
        return leftposition(a) < leftposition(b)
    end

    if rightposition(a) != rightposition(b)
        return rightposition(a) < rightposition(b)
    end

    if strand(a) != strand(b)
        return strand(a) < strand(b)
    end

    return false
end

"""
Check if two intervals are well ordered.

Intervals are considered well ordered if seqname(a) <= seqname(b) and leftposition(a) <= leftposition(b).
"""
function isordered(a::Interval, b::Interval, seqname_isless::Function=isless)
    if seqname(a) != seqname(b)
        return seqname_isless(seqname(a), seqname(b))::Bool
    end

    if leftposition(a) != leftposition(b)
        return leftposition(a) < leftposition(b)
    end

    return true
end

"""
Return true if interval `a` entirely precedes `b`.
"""
function precedes(a::Interval, b::Interval, seqname_isless::Function=isless)
    return (rightposition(a) < leftposition(b) && seqname(a) == seqname(b)) || seqname_isless(seqname(a), seqname(b))::Bool
end

function Base.:(==)(a::Interval, b::Interval)
    return seqname(a)       == seqname(b) &&
           leftposition(a)  == leftposition(b) &&
           rightposition(a) == rightposition(b) &&
           strand(a)        == strand(b) &&
           metadata(a)      == metadata(b)
end

"Return true if interval `a` overlaps interval `b`, with no consideration to strand"
function BioGenerics.isoverlapping(a::Interval, b::Interval)
    return leftposition(a) <= rightposition(b) &&
           leftposition(b) <= rightposition(a) &&
           seqname(a)      == seqname(b)
end

function Base.show(io::IO, i::Interval)
    if get(io, :compact, false)
        print(io, seqname(i), ":", leftposition(i), "-", rightposition(i), "  ", strand(i), "  ", metadata(i) === nothing ? "nothing" : metadata(i))
    else
        println(io, summary(i), ':')
        println(io, "  sequence name: ", seqname(i))
        println(io, "  leftmost position: ", leftposition(i))
        println(io, "  rightmost position: ", rightposition(i))
        println(io, "  strand: ", strand(i))
          print(io, "  metadata: ", metadata(i) === nothing ? "nothing" : metadata(i))
    end
end

function metadatatype(::Type{T}) where T
    return _metadatatype(eltype(T))
end

function metadatatype(x::Any)
    return metadatatype(typeof(x))
end

function _metadatatype(::Type{Interval{T}}) where T
    return T
end

function _metadatatype(::Type{T}) where T
    return T
end


text: include("Notebooks.jl")

max_stars_repo_name: MattAmos/Example-Pluto-CI

max_issues_repo_name: MattAmos/Example-Pluto-CI

max_forks_repo_name: MattAmos/Example-Pluto-CI


text:
# Extract the bow shock location from 2D equatorial run outputs and save into file.
#
# To run in multi-threading mode
#   julia -t nthreads demo_bowshock_series.jl
# or alternatively
#   JULIA_NUM_THREADS=$nthreads ./julia demo_bowshock_series.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Glob
using JLD2: jldsave

# Upstream solar wind temperature
const Tsw = 0.5e6 #[K]

function extract_bowshock_position(files; verbose=true)
   nfiles = length(files)

   verbose && println("Number of files: $nfiles")

   meta = load(files[1])

   x = LinRange{Float32}(meta.coordmin[1], meta.coordmax[1], meta.ncells[1])
   y = LinRange{Float32}(meta.coordmin[2], meta.coordmax[2], meta.ncells[2])

   close(meta.fid)

   # Only extract bow shock location near the front between y = ±20Re
   Re = Vlasiator.Re
   ymin_ = findlast(<(-20Re), y)
   ymax_ = findfirst(>(20Re), y)

   x_crossing = zeros(Float32, ymax_-ymin_+1, nfiles)
   y_crossing = y[ymin_:ymax_]

   Threads.@threads for ifile = 1:nfiles
      verbose && println("$(files[ifile]) on thread $(Threads.threadid())")
      f = load(files[ifile])
      # Obtain thermal temperature
      T = f["T"]
      close(f.fid)
      T = reshape(T, f.ncells[1], f.ncells[2])
      # Extract bow shock location from the 1st point which fulfills
      # the threshold: T > 4 * Tsw
      for (i,j) in enumerate(ymin_:ymax_) # scan in y direction
         ind_ = findlast(>(4*Tsw), @view T[:,j]) # count from upstream
         x_crossing[i,ifile] = x[ind_]
      end
   end

   return x_crossing, y_crossing
end

#####
files = glob("bulk*.vlsv", "run_rho2_bz-5_timevarying_startfrom300s")

@time x_crossing, y_crossing = extract_bowshock_position(files)

jldsave("example.jld2"; x_crossing, y_crossing)


text:
using NNlib: PoolDims

# descriptor

mutable struct PoolDesc
    ptr::cudnnPoolingDescriptor_t
end

unsafe_free!(pd::PoolDesc)=cudnnDestroyPoolingDescriptor(pd.ptr)

Base.unsafe_convert(::Type{cudnnPoolingDescriptor_t}, pd::PoolDesc)=pd.ptr

function PoolDesc(nd, window, padding, stride, mode, maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN)
    pd = Ref{cudnnPoolingDescriptor_t}()
    cudnnCreatePoolingDescriptor(pd)
    cudnnSetPoolingNdDescriptor(pd[],cudnnPoolingMode_t(mode),maxpoolingNanOpt,nd,pdsize(window,nd),pdsize(padding,nd),pdsize(stride,nd))
    this = PoolDesc(pd[])
    finalizer(unsafe_free!, this)
    return this
end

function PoolDesc(pdims::PoolDims, mode, maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN)
    pd = NNlib.padding(pdims)
    if !all(pd[1:2:end] .== pd[2:2:end])
        @warn("CuDNN does not support asymmetric padding; defaulting to symmetric choice")
    end
    return PoolDesc(NNlib.spatial_dims(pdims), NNlib.kernel_size(pdims), pd[1:2:end],
                    NNlib.stride(pdims), mode, maxpoolingNanOpt)
end

# wrappers

function cudnnPoolingForward(y::CuArray{T,N}, x::CuArray{T,N}, pdims::PoolDims;
                             alpha=1, mode=0) where {T,N}
    beta = 0
    cudnnPoolingForward(handle(), PoolDesc(pdims, mode), Ref(T(alpha)), TensorDesc(x), x, Ref(T(beta)), TensorDesc(y), y)
    return y
end

function cudnnPoolingBackward(dx::CuArray{T,N}, dy::CuArray{T,N}, x::CuArray{T,N}, y::CuArray{T,N},
                              pdims::PoolDims; alpha=1, mode=0) where {T,N}
    if alpha!=1 && mode==0; error("Gradient of pool(alpha!=1,mode=0) broken in CUDNN"); end
    beta = 0
    cudnnPoolingBackward(handle(),
        PoolDesc(pdims, mode), Ref(T(alpha)), TensorDesc(y), y,
        TensorDesc(dy), dy, TensorDesc(x), x, Ref(T(beta)), TensorDesc(dx), dx
    )
    return dx
end


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector
using LinearAlgebra: norm
using AdmittanceModels: lossless_modes_dense, apply_transform, get_Y

# Modes of a rectangular box with constant μ and ϵ.
box_modes(m, n, l, a, b, c, μ, ϵ) = sqrt((m * π/a)^2 + (n * π/b)^2 +
    (l * π/c)^2)/(2π * sqrt(μ * ϵ))
# The lowest mode of our box is the 011 mode.
a, b, c = 10, 12, 14
correct_A(x, y, z) = sin((π/b) * y) * sin((π/c) * z) * [1, 0, 0]

file_name = joinpath(@__DIR__, "box.geo")
geo_write!(file_name, characteristic_length_factor=1,
    footer="""
    Box(1) = {0, 0, 0, $a, $b, $c};
    """)

initialize!()
gmsh_open!(file_name)
N, K = 3, 4
mesh!(K)
node_tags, points, tcomp = get_triangulated_complex(N, K)
group_dict = get_physical_groups(node_tags, points)
@test typeof(tcomp) <: TriangulatedComplex{N, K}
comp = tcomp.complex
orient!(comp)
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

μ⁻, ϵ = 2, 3
μ⁻_form = get_material(comp, μ⁻, 3)
Λ⁻_form = get_material(comp, 0, 2)
σ_form = get_material(comp, 0, 2)
ϵ_form = get_material(comp, ϵ, 2)

pso, null_basis = electrodynamics_pso(m, mesh, Vector{Cell{N}}[], boundary(comp),
    μ⁻_form, Λ⁻_form, σ_form, ϵ_form)
constrained_pso = apply_transform(pso, null_basis)

density(mat) = count(!iszero, mat)/(size(mat, 1) * size(mat, 2))
@test density(null_basis) < .002
@test all(map(density, get_Y(constrained_pso)) .< .05)

λs, vs = lossless_modes_dense(constrained_pso)
freq = imag(λs[1])/(2π)
correct_freq = box_modes(0,1,1,a,b,c,1/μ⁻,ϵ)

comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
correct_vec_A = [correct_A(p.coords...) for p in comp_points]
correct_vec_A /= maximum(norm.(correct_vec_A))
vec_A = sharp(m, comp, null_basis * vs[:,1])
vec_A /= maximum(norm.(vec_A))
i = argmax(norm.(correct_vec_A))
vec_A *= sign(transpose(vec_A[i]) * correct_vec_A[i])

@testset "box modes" begin
    @test isapprox(freq, correct_freq, rtol=5e-3)
    @test maximum([norm(x - y) for (x, y) in zip(vec_A, correct_vec_A)]) < .25
end

if false # set to true to display plots with Gmsh
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    add_field!("Vector potential", node_tags, vec_A[ordering])
    add_field!("Correct vector potential", node_tags, correct_vec_A[ordering])
    gui!()
end
rm(file_name)