text:
module NumericalFluxes

export NumericalFluxGradient,
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux,
    RoeNumericalFlux,
    HLLCNumericalFlux,
    RoeNumericalFluxMoist,
    CentralNumericalFluxGradient,
    CentralNumericalFluxFirstOrder,
    CentralNumericalFluxSecondOrder,
    CentralNumericalFluxDivergence,
    CentralNumericalFluxHigherOrder,
    LMARSNumericalFlux

using StaticArrays, LinearAlgebra
using ClimateMachine.VariableTemplates
using KernelAbstractions.Extras: @unroll
using ...BalanceLaws
import ...BalanceLaws:
    vars_state,
    boundary_state!,
    wavespeed,
    flux_first_order!,
    flux_second_order!,
    compute_gradient_flux!,
    compute_gradient_argument!,
    transform_post_gradient_laplacian!,
    boundary_conditions

"""
    NumericalFluxGradient

Any `P <: NumericalFluxGradient` should define methods for:

    numerical_flux_gradient!(
        gnf::P,
        balance_law::BalanceLaw,
        diffF, n⁻,
        Q⁻, Qstate_gradient_flux⁻, Qaux⁻,
        Q⁺, Qstate_gradient_flux⁺, Qaux⁺,
        t
    )

    numerical_boundary_flux_gradient!(
        gnf::P,
        balance_law::BalanceLaw,
        local_state_gradient_flux,
        n⁻,
        local_transform⁻, local_state_prognostic⁻, local_state_auxiliary⁻,
        local_transform⁺, local_state_prognostic⁺, local_state_auxiliary⁺,
        bctype,
        t
    )

"""
abstract type NumericalFluxGradient end

"""
    CentralNumericalFluxGradient <: NumericalFluxGradient

"""
struct CentralNumericalFluxGradient <: NumericalFluxGradient end

function numerical_flux_gradient!(
    ::CentralNumericalFluxGradient,
    balance_law::BalanceLaw,
    transform_gradient::MMatrix,
    normal_vector::AbstractArray,
    state_gradient⁻::AbstractArray,
    state_prognostic⁻::AbstractArray,
    state_auxiliary⁻::AbstractArray,
    state_gradient⁺::AbstractArray,
    state_prognostic⁺::AbstractArray,
    state_auxiliary⁺::AbstractArray,
    t,
)

    transform_gradient .=
        SVector(normal_vector) .* (state_gradient⁺ .+ state_gradient⁻)' ./ 2
end

function numerical_boundary_flux_gradient!(
    numerical_flux::CentralNumericalFluxGradient,
    bctype,
    balance_law::BalanceLaw,
    transform_gradient::MMatrix,
    normal_vector::SVector,
    state_gradient⁻::Vars{T},
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_gradient⁺::Vars{T},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    aux1⁻::Vars{A},
) where {D, T, S, A}
    boundary_state!(
        numerical_flux,
        bctype,
        balance_law,
        state_prognostic⁺,
        state_auxiliary⁺,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
        state1⁻,
        aux1⁻,
    )

    compute_gradient_argument!(
        balance_law,
        state_gradient⁺,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
    )
    transform_gradient .= normal_vector .* parent(state_gradient⁺)'
end

"""
    NumericalFluxFirstOrder

Any `N <: NumericalFluxFirstOrder` should define the a method for

    numerical_flux_first_order!(
        numerical_flux::N,
        balance_law::BalanceLaw,
        flux,
        normal_vector⁻,
        Q⁻, Qaux⁻,
        Q⁺, Qaux⁺,
        t
    )

where
- `flux` is the numerical flux array
- `normal_vector⁻` is the unit normal
- `Q⁻`/`Q⁺` are the minus/positive state arrays
- `t` is the time

An optional method can also be defined for

    numerical_boundary_flux_first_order!(
        numerical_flux::N,
        balance_law::BalanceLaw,
        flux,
        normal_vector⁻,
        Q⁻, Qaux⁻,
        Q⁺, Qaux⁺,
        bctype, t
    )

"""
abstract type NumericalFluxFirstOrder end

function numerical_flux_first_order! end

function numerical_boundary_flux_first_order!(
    numerical_flux::NumericalFluxFirstOrder,
    bctype,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
    state1⁻::Vars{S},
    aux1⁻::Vars{A},
) where {S, A}

    boundary_state!(
        numerical_flux,
        bctype,
        balance_law,
        state_prognostic⁺,
        state_auxiliary⁺,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
        state1⁻,
        aux1⁻,
    )

    numerical_flux_first_order!(
        numerical_flux,
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
end

"""
    RusanovNumericalFlux <: NumericalFluxFirstOrder

The RusanovNumericalFlux (aka local Lax-Friedrichs) numerical flux.

# Usage

    RusanovNumericalFlux()

Requires a `flux_first_order!` and `wavespeed` method for the balance law.
"""
struct RusanovNumericalFlux <: NumericalFluxFirstOrder end

update_penalty!(::RusanovNumericalFlux, ::BalanceLaw, _...) = nothing

function numerical_flux_first_order!(
    numerical_flux::RusanovNumericalFlux,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

    fluxᵀn = parent(fluxᵀn)
    wavespeed⁻ = wavespeed(
        balance_law,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
        direction,
    )
    wavespeed⁺ = wavespeed(
        balance_law,
        normal_vector,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    max_wavespeed = max.(wavespeed⁻, wavespeed⁺)
    penalty =
        max_wavespeed .* (parent(state_prognostic⁻) - parent(state_prognostic⁺))

    # TODO: should this operate on ΔQ or penalty?
    update_penalty!(
        numerical_flux,
        balance_law,
        normal_vector,
        max_wavespeed,
        Vars{S}(penalty),
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
    )

    fluxᵀn .+= penalty / 2
end

"""
    CentralNumericalFluxFirstOrder() <: NumericalFluxFirstOrder

The central numerical flux for nondiffusive terms

# Usage

    CentralNumericalFluxFirstOrder()

Requires a `flux_first_order!` method for the balance law.
"""
struct CentralNumericalFluxFirstOrder <: NumericalFluxFirstOrder end

function numerical_flux_first_order!(
    ::CentralNumericalFluxFirstOrder,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}

    FT = eltype(fluxᵀn)
    num_state_prognostic = number_states(balance_law, Prognostic())
    fluxᵀn = parent(fluxᵀn)

    flux⁻ = similar(fluxᵀn, Size(3, num_state_prognostic))
    fill!(flux⁻, -zero(FT))
    flux_first_order!(
        balance_law,
        Grad{S}(flux⁻),
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
        direction,
    )

    flux⁺ = similar(fluxᵀn, Size(3, num_state_prognostic))
    fill!(flux⁺, -zero(FT))
    flux_first_order!(
        balance_law,
        Grad{S}(flux⁺),
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

    fluxᵀn .+= (flux⁻ + flux⁺)' * (normal_vector / 2)
end

"""
    RoeNumericalFlux() <: NumericalFluxFirstOrder

A numerical flux based on the approximate Riemann solver of Roe

# Usage

    RoeNumericalFlux()

Requires a custom implementation for the balance law.
"""
struct RoeNumericalFlux <: NumericalFluxFirstOrder end

"""
    HLLCNumericalFlux() <: NumericalFluxFirstOrder

A numerical flux based on the approximate Riemann solver of the
HLLC method. The HLLC flux is a modification of the Harten, Lax, van-Leer
(HLL) flux, where an additional contact property is introduced in order
to restore missing rarefraction waves. The HLLC flux requires
model-specific information, hence it requires a custom implementation
based on the underlying balance law.

# Usage

    HLLCNumericalFlux()

Requires a custom implementation for the balance law.

 - [Toro2013](@cite)
"""
struct HLLCNumericalFlux <: NumericalFluxFirstOrder end

"""
    LMARSNumericalFlux <: NumericalFluxFirstOrder
Low Mach Number Approximate Riemann Solver. Upwind biased
first order flux function.

- [Chen2013](@cite)
"""
struct LMARSNumericalFlux <: NumericalFluxFirstOrder end

"""
    RoeNumericalFluxMoist <: NumericalFluxFirstOrder

A moist implementation of the numerical flux based on the approximate Riemann solver of Roe

Requires a custom implementation for the balance law.
"""
struct RoeNumericalFluxMoist <: NumericalFluxFirstOrder
    " set to true for low Mach number correction"
    LM::Bool
    " set to true for Hartman Hyman correction"
    HH::Bool
    " set to true for LeVeque correction"
    LV::Bool
    " set to true for Positivity preserving LeVeque Correction"
    LVPP::Bool
end

RoeNumericalFluxMoist(;
    LM::Bool = false,
    HH::Bool = false,
    LV::Bool = false,
    LVPP::Bool = false,
) = RoeNumericalFluxMoist(LM, HH, LV, LVPP)

"""
dispatch type for entropy conservative numerical fluxes (these are balance law
specific)
"""
struct EntropyConservative <: NumericalFluxFirstOrder end

"""
    numerical_volume_conservative_flux_first_order!(
        numflux::NumericalFluxFirstOrder,
        balancelaw::BalanceLaw,
        flux::Grad,
        state_1::Vars,
        aux_1::Vars,
        state_2::Vars,
        aux_2::Vars,
    )
Two point flux for use in the volume of entropy stable discretizations. Should
compute and store the conservative (symmetric) part of an entropy stable
flux-fluctuation and store the three vector components back to `flux`.
This should be symmetric in the sense that swapping the place of `state_1`,
`aux_1` with `state_2`, `aux_2` should result in the same flux.
This should be implemented with addition assignment.
Each balance law must implement a concrete implementation of this function.
"""
function numerical_volume_conservative_flux_first_order!(
    numflux::NumericalFluxFirstOrder,
    balancelaw::BalanceLaw,
    flux::AbstractArray{FT, 2},
    state_1::AbstractArray{FT, 1},
    aux_1::AbstractArray{FT, 1},
    state_2::AbstractArray{FT, 1},
    aux_2::AbstractArray{FT, 1},
) where {FT}
    numerical_volume_conservative_flux_first_order!(
        numflux,
        balancelaw,
        Grad{vars_state(balancelaw, Prognostic(), FT)}(flux),
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state_1),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux_1),
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state_2),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux_2),
    )
end
numerical_volume_conservative_flux_first_order!(::Nothing, _...) = nothing

"""
    numerical_volume_fluctuation_flux_first_order!(
        numflux::NumericalFluxFirstOrder,
        balancelaw::BalanceLaw,
        flux::Grad,
        state_1::Vars,
        aux_1::Vars,
        state_2::Vars,
        aux_2::Vars,
    )
Two point fluctuation flux for use in the volume of entropy stable
discretizations. Should compute and store the non-conservative (possibly
non-symmetric) part of an entropy stable flux-fluctuation and store the three
vector components back to `flux`.
This can be non-symmetric in the sense that swapping the place of `state_1`,
`aux_1` with `state_2`, `aux_2` can result in a different flux.
This should be implemented with addition assignment.
Each balance law must implement a concrete implementation of this function.
"""
function numerical_volume_fluctuation_flux_first_order!(
    numflux::NumericalFluxFirstOrder,
    balancelaw::BalanceLaw,
    fluctuation::AbstractArray{FT, 2},
    state_1::AbstractArray{FT, 1},
    aux_1::AbstractArray{FT, 1},
    state_2::AbstractArray{FT, 1},
    aux_2::AbstractArray{FT, 1},
) where {FT}
    numerical_volume_fluctuation_flux_first_order!(
        numflux,
        balancelaw,
        Grad{vars_state(balancelaw, Prognostic(), FT)}(fluctuation),
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state_1),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux_1),
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state_2),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux_2),
    )
end
numerical_volume_fluctuation_flux_first_order!(::Nothing, _...) = nothing

"""
    numerical_volume_flux_first_order!(
        numflux::NumericalFluxFirstOrder,
        balancelaw::BalanceLaw,
        flux::AbstractArray{FT, 2},
        state_1::AbstractArray{FT, 1},
        aux_1::AbstractArray{FT, 1},
        state_2::AbstractArray{FT, 1},
        aux_2::AbstractArray{FT, 1},
    )
Convenience function which calls
  `numerical_volume_conservative_flux_first_order!`
followed by
  `numerical_volume_fluctuation_flux_first_order!`
"""
function numerical_volume_flux_first_order!(
    numflux::NumericalFluxFirstOrder,
    balancelaw::BalanceLaw,
    flux::AbstractArray{FT, 2},
    state_1::AbstractArray{FT, 1},
    aux_1::AbstractArray{FT, 1},
    state_2::AbstractArray{FT, 1},
    aux_2::AbstractArray{FT, 1},
) where {FT}
    numerical_volume_conservative_flux_first_order!(
        numflux,
        balancelaw,
        flux,
        state_1,
        aux_1,
        state_2,
        aux_2,
    )
    numerical_volume_fluctuation_flux_first_order!(
        numflux,
        balancelaw,
        flux,
        state_1,
        aux_1,
        state_2,
        aux_2,
    )
end
numerical_volume_flux_first_order!(::Nothing, _...) = nothing

function numerical_flux_first_order!(
    numerical_flux::EntropyConservative,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}
    # Let's just work with straight arrays!
    fluxᵀn = parent(fluxᵀn)
    state_1 = parent(state_prognostic⁻)
    aux_1 = parent(state_auxiliary⁻)
    state_2 = parent(state_prognostic⁺)
    aux_2 = parent(state_auxiliary⁺)

    # create the storage for the volume flux
    num_state = length(fluxᵀn)
    FT = eltype(fluxᵀn)
    H = MArray{Tuple{3, num_state}, FT}(undef)
    fill!(H, -zero(FT))

    # Compute the volume flux
    numerical_volume_flux_first_order!(
        numerical_flux,
        balance_law,
        H,
        state_1,
        aux_1,
        state_2,
        aux_2,
    )

    # Multiply in the normal
    @inbounds fluxᵀn .+=
        normal_vector[1] * H[1, :] +
        normal_vector[2] * H[2, :] +
        normal_vector[3] * H[3, :]
end

# Convenience function for entropy stable discretizations
"""
    ave(a, b)
This computes the mean
    ave(a, b) = (a + b) / 2
"""
ave(a, b) = (a + b) / 2

# Convenience function for entropy stable discretizations
"""
    logave(a, b)
This computes the logarithmic mean
    logave(a, b) = (a - b) / (log(a) - log(b))
in a numerically stable way using the method in Appendix B. of Ishail and Roe
<doi:10.1016/j.jcp.2009.04.021>.
"""
function logave(a, b)
    ζ = a / b
    f = (ζ - 1) / (ζ + 1)
    u = f^2
    ϵ = eps(eltype(u))

    if u < ϵ
        F = @evalpoly(u, one(u), one(u) / 3, one(u) / 5, one(u) / 7, one(u) / 9)
    else
        F = log(ζ) / 2f
    end

    (a + b) / 2F
end

"""
    NumericalFluxSecondOrder

Any `N <: NumericalFluxSecondOrder` should define the a method for

    numerical_flux_second_order!(
        numerical_flux::N,
        balance_law::BalanceLaw,
        flux,
        normal_vector⁻,
        Q⁻, Qstate_gradient_flux⁻, Qaux⁻,
        Q⁺, Qstate_gradient_flux⁺, Qaux⁺,
        t
    )

where
- `flux` is the numerical flux array
- `normal_vector⁻` is the unit normal
- `Q⁻`/`Q⁺` are the minus/positive state arrays
- `Qstate_gradient_flux⁻`/`Qstate_gradient_flux⁺` are the minus/positive diffusive state arrays
- `Qstate_gradient_flux⁻`/`Qstate_gradient_flux⁺` are the minus/positive auxiliary state arrays
- `t` is the time

An optional method can also be defined for

    numerical_boundary_flux_second_order!(
        numerical_flux::N,
        balance_law::BalanceLaw,
        flux,
        normal_vector⁻,
        Q⁻, Qstate_gradient_flux⁻, Qaux⁻,
        Q⁺, Qstate_gradient_flux⁺, Qaux⁺,
        bctype,
        t
    )

"""
abstract type NumericalFluxSecondOrder end

function numerical_flux_second_order! end

function numerical_boundary_flux_second_order! end

"""
    CentralNumericalFluxSecondOrder <: NumericalFluxSecondOrder

The central numerical flux for diffusive terms

# Usage

    CentralNumericalFluxSecondOrder()

Requires a `flux_second_order!` for the balance law.
"""
struct CentralNumericalFluxSecondOrder <: NumericalFluxSecondOrder end

function numerical_flux_second_order!(
    ::CentralNumericalFluxSecondOrder,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector⁻::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
) where {S, D, HD, A}

    FT = eltype(fluxᵀn)
    num_state_prognostic = number_states(balance_law, Prognostic())
    fluxᵀn = parent(fluxᵀn)

    flux⁻ = similar(fluxᵀn, Size(3, num_state_prognostic))
    fill!(flux⁻, -zero(FT))
    flux_second_order!(
        balance_law,
        Grad{S}(flux⁻),
        state_prognostic⁻,
        state_gradient_flux⁻,
        state_hyperdiffusive⁻,
        state_auxiliary⁻,
        t,
    )

    flux⁺ = similar(fluxᵀn, Size(3, num_state_prognostic))
    fill!(flux⁺, -zero(FT))
    flux_second_order!(
        balance_law,
        Grad{S}(flux⁺),
        state_prognostic⁺,
        state_gradient_flux⁺,
        state_hyperdiffusive⁺,
        state_auxiliary⁺,
        t,
    )

    fluxᵀn .+= (flux⁻ + flux⁺)' * (normal_vector⁻ / 2)
end

abstract type DivNumericalPenalty end
struct CentralNumericalFluxDivergence <: DivNumericalPenalty end

function numerical_flux_divergence!(
    ::CentralNumericalFluxDivergence,
    balance_law::BalanceLaw,
    div_penalty::Vars{GL},
    normal_vector::SVector,
    grad⁻::Grad{GL},
    grad⁺::Grad{GL},
) where {GL}
    parent(div_penalty) .=
        (parent(grad⁺) .+ parent(grad⁻))' * (normal_vector / 2)
end

function numerical_boundary_flux_divergence!(
    numerical_flux::CentralNumericalFluxDivergence,
    bctype,
    balance_law::BalanceLaw,
    div_penalty::Vars{GL},
    normal_vector::SVector,
    grad⁻::Grad{GL},
    state_auxiliary⁻::Vars{A},
    grad⁺::Grad{GL},
    state_auxiliary⁺::Vars{A},
    t,
) where {GL, A}
    boundary_state!(
        numerical_flux,
        bctype,
        balance_law,
        grad⁺,
        state_auxiliary⁺,
        normal_vector,
        grad⁻,
        state_auxiliary⁻,
        t,
    )
    numerical_flux_divergence!(
        numerical_flux,
        balance_law,
        div_penalty,
        normal_vector,
        grad⁻,
        grad⁺,
    )
end

abstract type GradNumericalFlux end
struct CentralNumericalFluxHigherOrder <: GradNumericalFlux end

function numerical_flux_higher_order!(
    ::CentralNumericalFluxHigherOrder,
    balance_law::BalanceLaw,
    hyperdiff::Vars{HD},
    normal_vector::SVector,
    lap⁻::Vars{GL},
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    lap⁺::Vars{GL},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
) where {HD, GL, S, A}
    G = normal_vector .* (parent(lap⁺) .- parent(lap⁻))' ./ 2
    transform_post_gradient_laplacian!(
        balance_law,
        hyperdiff,
        Grad{GL}(G),
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
    )
end

function numerical_boundary_flux_higher_order!(
    numerical_flux::CentralNumericalFluxHigherOrder,
    bctype,
    balance_law::BalanceLaw,
    hyperdiff::Vars{HD},
    normal_vector::SVector,
    lap⁻::Vars{GL},
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    lap⁺::Vars{GL},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
) where {HD, GL, S, A}
    boundary_state!(
        numerical_flux,
        bctype,
        balance_law,
        state_prognostic⁺,
        state_auxiliary⁺,
        lap⁺,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        lap⁻,
        t,
    )
    numerical_flux_higher_order!(
        numerical_flux,
        balance_law,
        hyperdiff,
        normal_vector,
        lap⁻,
        state_prognostic⁻,
        state_auxiliary⁻,
        lap⁺,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
    )
end

numerical_boundary_flux_second_order!(
    numerical_flux::CentralNumericalFluxSecondOrder,
    bctype,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_gradient_flux⁻::Vars{D},
    state_hyperdiffusive⁻::Vars{HD},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_gradient_flux⁺::Vars{D},
    state_hyperdiffusive⁺::Vars{HD},
    state_auxiliary⁺::Vars{A},
    t,
    state1⁻::Vars{S},
    diff1⁻::Vars{D},
    aux1⁻::Vars{A},
) where {S, D, HD, A} = normal_boundary_flux_second_order!(
    numerical_flux,
    bctype,
    balance_law,
    fluxᵀn,
    normal_vector,
    state_prognostic⁻,
    state_gradient_flux⁻,
    state_hyperdiffusive⁻,
    state_auxiliary⁻,
    state_prognostic⁺,
    state_gradient_flux⁺,
    state_hyperdiffusive⁺,
    state_auxiliary⁺,
    t,
    state1⁻,
    diff1⁻,
    aux1⁻,
)

function normal_boundary_flux_second_order!(
    numerical_flux,
    bctype,
    balance_law::BalanceLaw,
    fluxᵀn::Vars{S},
    normal_vector,
    state_prognostic⁻,
    state_gradient_flux⁻,
    state_hyperdiffusive⁻,
    state_auxiliary⁻,
    state_prognostic⁺,
    state_gradient_flux⁺,
    state_hyperdiffusive⁺,
    state_auxiliary⁺,
    t,
    state1⁻,
    diff1⁻,
    aux1⁻,
) where {S}
    FT = eltype(fluxᵀn)
    num_state_prognostic = number_states(balance_law, Prognostic())
    fluxᵀn = parent(fluxᵀn)

    flux = similar(fluxᵀn, Size(3, num_state_prognostic))
    fill!(flux, -zero(FT))
    boundary_flux_second_order!(
        numerical_flux,
        bctype,
        balance_law,
        Grad{S}(flux),
        state_prognostic⁺,
        state_gradient_flux⁺,
        state_hyperdiffusive⁺,
        state_auxiliary⁺,
        normal_vector,
        state_prognostic⁻,
        state_gradient_flux⁻,
        state_hyperdiffusive⁻,
        state_auxiliary⁻,
        t,
        state1⁻,
        diff1⁻,
        aux1⁻,
    )

    fluxᵀn .+= flux' * normal_vector
end

# This is the function that my be overloaded for flux-based BCs
function boundary_flux_second_order!(
    numerical_flux::NumericalFluxSecondOrder,
    bctype,
    balance_law,
    flux,
    state_prognostic⁺,
    state_gradient_flux⁺,
    state_hyperdiffusive⁺,
    state_auxiliary⁺,
    normal_vector,
    state_prognostic⁻,
    state_gradient_flux⁻,
    state_hyperdiffusive⁻,
    state_auxiliary⁻,
    t,
    state1⁻,
    diff1⁻,
    aux1⁻,
)
    boundary_state!(
        numerical_flux,
        bctype,
        balance_law,
        state_prognostic⁺,
        state_gradient_flux⁺,
        state_hyperdiffusive⁺,
        state_auxiliary⁺,
        normal_vector,
        state_prognostic⁻,
        state_gradient_flux⁻,
        state_hyperdiffusive⁻,
        state_auxiliary⁻,
        t,
        state1⁻,
        diff1⁻,
        aux1⁻,
    )
    flux_second_order!(
        balance_law,
        flux,
        state_prognostic⁺,
        state_gradient_flux⁺,
        state_hyperdiffusive⁺,
        state_auxiliary⁺,
        t,
    )
end

# Array wrappers for numerical Fluxes
function numerical_flux_first_order!(
    numerical_flux_first_order,
    balance_law,
    flux::AbstractArray,
    normal_vector::AbstractArray,
    state_prognostic⁻::AbstractArray,
    state_auxiliary⁻::AbstractArray,
    state_prognostic⁺::AbstractArray,
    state_auxiliary⁺::AbstractArray,
    t,
    face_direction,
)
    FT = eltype(flux)
    numerical_flux_first_order!(
        numerical_flux_first_order,
        balance_law,
        Vars{vars_state(balance_law, Prognostic(), FT)}(flux),
        SVector(normal_vector),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁻),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁻),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁺),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁺),
        t,
        face_direction,
    )
end

function numerical_boundary_flux_first_order!(
    numerical_flux_first_order,
    bctag::Int,
    balance_law,
    flux::AbstractArray,
    normal_vector::AbstractArray,
    state_prognostic⁻::AbstractArray,
    state_auxiliary⁻::AbstractArray,
    state_prognostic⁺::AbstractArray,
    state_auxiliary⁺::AbstractArray,
    t,
    face_direction,
    state_prognostic_bottom1::AbstractArray,
    state_auxiliary_bottom1::AbstractArray,
)
    FT = eltype(flux)
    bcs = boundary_conditions(balance_law)
    # TODO: there is probably a better way to unroll this loop
    Base.Cartesian.@nif 7 d -> bctag == d <= length(bcs) d -> begin
        bc = bcs[d]
        numerical_boundary_flux_first_order!(
            numerical_flux_first_order,
            bc,
            balance_law,
            Vars{vars_state(balance_law, Prognostic(), FT)}(flux),
            SVector(normal_vector),
            Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁻),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁻),
            Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁺),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁺),
            t,
            face_direction,
            Vars{vars_state(balance_law, Prognostic(), FT)}(
                state_prognostic_bottom1,
            ),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(
                state_auxiliary_bottom1,
            ),
        )
    end d -> throw(BoundsError(bcs, bctag))
end

function numerical_flux_second_order!(
    numerical_flux_second_order,
    balance_law,
    flux::AbstractArray,
    normal_vector::AbstractArray,
    state_prognostic⁻::AbstractArray,
    state_gradient_flux⁻::AbstractArray,
    state_hyperdiffusive⁻::AbstractArray,
    state_auxiliary⁻::AbstractArray,
    state_prognostic⁺::AbstractArray,
    state_gradient_flux⁺::AbstractArray,
    state_hyperdiffusive⁺::AbstractArray,
    state_auxiliary⁺::AbstractArray,
    t,
)
    FT = eltype(flux)
    numerical_flux_second_order!(
        numerical_flux_second_order,
        balance_law,
        Vars{vars_state(balance_law, Prognostic(), FT)}(flux),
        SVector(normal_vector),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁻),
        Vars{vars_state(balance_law, GradientFlux(), FT)}(state_gradient_flux⁻),
        Vars{vars_state(balance_law, Hyperdiffusive(), FT)}(
            state_hyperdiffusive⁻,
        ),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁻),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁺),
        Vars{vars_state(balance_law, GradientFlux(), FT)}(state_gradient_flux⁺),
        Vars{vars_state(balance_law, Hyperdiffusive(), FT)}(
            state_hyperdiffusive⁺,
        ),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁺),
        t,
    )
end

function numerical_boundary_flux_second_order!(
    numerical_flux_second_order,
    bctag::Int,
    balance_law,
    flux::AbstractArray,
    normal_vector::AbstractArray,
    state_prognostic⁻::AbstractArray,
    state_gradient_flux⁻::AbstractArray,
    state_hyperdiffusive⁻::AbstractArray,
    state_auxiliary⁻::AbstractArray,
    state_prognostic⁺::AbstractArray,
    state_gradient_flux⁺::AbstractArray,
    state_hyperdiffusive⁺::AbstractArray,
    state_auxiliary⁺::AbstractArray,
    t,
    state_prognostic_bottom1::AbstractArray,
    state_auxiliary_bottom1::AbstractArray,
    state_gradient_flux_bottom1::AbstractArray,
)
    FT = eltype(flux)
    bcs = boundary_conditions(balance_law)
    # TODO: there is probably a better way to unroll this loop
    Base.Cartesian.@nif 7 d -> bctag == d <= length(bcs) d -> begin
        bc = bcs[d]
        numerical_boundary_flux_second_order!(
            numerical_flux_second_order,
            bc,
            balance_law,
            Vars{vars_state(balance_law, Prognostic(), FT)}(flux),
            SVector(normal_vector),
            Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁻),
            Vars{vars_state(balance_law, GradientFlux(), FT)}(
                state_gradient_flux⁻,
            ),
            Vars{vars_state(balance_law, Hyperdiffusive(), FT)}(
                state_hyperdiffusive⁻,
            ),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁻),
            Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁺),
            Vars{vars_state(balance_law, GradientFlux(), FT)}(
                state_gradient_flux⁺,
            ),
            Vars{vars_state(balance_law, Hyperdiffusive(), FT)}(
                state_hyperdiffusive⁺,
            ),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁺),
            t,
            Vars{vars_state(balance_law, Prognostic(), FT)}(
                state_prognostic_bottom1,
            ),
            Vars{vars_state(balance_law, GradientFlux(), FT)}(
                state_gradient_flux_bottom1,
            ),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(
                state_auxiliary_bottom1,
            ),
        )
    end d -> throw(BoundsError(bcs, bctag))
end

function numerical_boundary_flux_gradient!(
    numerical_flux_gradient,
    bctag::Int,
    balance_law,
    flux::AbstractArray,
    normal_vector::AbstractArray,
    state_gradient⁻::AbstractArray,
    state_prognostic⁻::AbstractArray,
    state_auxiliary⁻::AbstractArray,
    state_gradient⁺::AbstractArray,
    state_prognostic⁺::AbstractArray,
    state_auxiliary⁺::AbstractArray,
    t,
    state_prognostic_bottom1::AbstractArray,
    state_auxiliary_bottom1::AbstractArray,
)
    FT = eltype(flux)
    bcs = boundary_conditions(balance_law)
    # TODO: there is probably a better way to unroll this loop
    Base.Cartesian.@nif 7 d -> bctag == d <= length(bcs) d -> begin
        bc = bcs[d]
        # Computes G* incorporating boundary conditions
        numerical_boundary_flux_gradient!(
            numerical_flux_gradient,
            bc,
            balance_law,
            flux,
            SVector(normal_vector),
            Vars{vars_state(balance_law, Gradient(), FT)}(state_gradient⁻),
            Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁻),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁻),
            Vars{vars_state(balance_law, Gradient(), FT)}(state_gradient⁺),
            Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic⁺),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary⁺),
            t,
            Vars{vars_state(balance_law, Prognostic(), FT)}(
                state_prognostic_bottom1,
            ),
            Vars{vars_state(balance_law, Auxiliary(), FT)}(
                state_auxiliary_bottom1,
            ),
        )
    end d -> throw(BoundsError(bcs, bctag))
end

end # end of module


text:
const DOC_OBSERVATIONS =
    "on observations `ŷ`, "*
    "given ground truth values `y`. "
const DOC_WEIGHTS =
    "Optionally specify per-sample weights, `w`. "
const DOC_CLASS_WEIGHTS =
    "An optional `AbstractDict`, denoted `class_w` above, "*
    "keyed on `levels(y)`, specifies class weights. "

macro create_aliases(M_ex)
    esc(quote
        M = $M_ex
        for alias in Symbol.(instances(M))
        # isdefined(parentmodule(M), alias) || eval(:(const $alias = $M()))
        eval(:(const $alias = $M()))
        end
        end)
end

function detailed_doc_string(M; typename="", body="", footer="", scitype="")

    _instances = _decorate(instances(M))
    human_name = MLJBase.human_name(M)
    if isempty(scitype)
        scitype = target_scitype(M) |> string
    end

    if isempty(typename)
        ret = "    $M\n\n"
    else
        ret = "    MLJBase.$typename\n\n"
    end

    ret *= "A measure type for $(human_name)"
    isempty(_instances) ||
        (ret  *= ", which includes the instance(s): "*
         "$_instances")
    ret *= ".\n\n"
    ret *= "    $(name(M))()(ŷ, y)\n"
    supports_weights(M) &&
        (ret *= "    $(name(M))()(ŷ, y, w)\n")
    supports_class_weights(M) &&
        (ret *= "    $(name(M))()(ŷ, y, class_w)\n")
    ret *= "\n"
    if isempty(fieldnames(M))
            ret *= "Evaluate the $(human_name) "
    else
        ret *= "Evaluate the default instance of $(name(M)) "
    end
    ret *= "$DOC_OBSERVATIONS"
    supports_weights(M) &&
        (ret *= DOC_WEIGHTS)
    supports_class_weights(M) &&
        (ret *= DOC_CLASS_WEIGHTS)
    ret *= "\n\n"
    isempty(body) || (ret *= "$body\n\n")
    ret *= "Requires `scitype(y)` to be a subtype of $scitype; "
    ret *= "`ŷ` must be a $(prediction_type(M)) prediction. "
    isempty(footer) ||(ret *= "\n\n$footer")
    ret *= "\n\n"
    ret *= "For more information, run `info($(name(M)))`. "
    return ret
end

_err_create_docs() = error(
    "@create_docs syntax error. Usage: \n"*
    "@create_docs(MeasureType, typename=..., body=..., scitype=..., footer=...")
macro create_docs(M_ex, exs...)
    M_ex isa Symbol || _err_create_docs()
    t = ""
    b = ""
    s = ""
    f = ""
    for ex in exs
        ex.head == :(=) || _err_create_docs()
        ex.args[1] == :typename && (t = ex.args[2])
        ex.args[1] == :body &&     (b = ex.args[2])
        ex.args[1] == :scitype &&  (s = ex.args[2])
        ex.args[1] == :footer &&   (f = ex.args[2])
    end
    esc(quote
        "$(detailed_doc_string($M_ex, typename=$t, body=$b, scitype=$s, footer=$f))"
        function $M_ex end
        end)
end

# TODO: I wonder why this is not a macro?

"""
    metadata_measure(T; kw...)

Helper function to write the metadata for a single measure.
"""
function metadata_measure(T; name::String="",
                          human_name="",
                          instances::Vector{String}=String[],
                          target_scitype=Unknown,
                          prediction_type::Symbol=:unknown,
                          orientation::Symbol=:unknown,
                          reports_each_observation::Bool=true,
                          aggregation=Mean(),
                          is_feature_dependent::Bool=false,
                          supports_weights::Bool=false,
                          supports_class_weights::Bool=false,
                          docstring::String="",
                          distribution_type=missing)
    pred_str        = "$prediction_type"
    orientation_str = "$orientation"
#    dist = ifelse(ismissing(distribution_type), missing, "$distribution_type")
    ex = quote

        # traits common with models:
        if !isempty($name)
            MMI.name(::Type{<:$T}) = $name
        end
        if !isempty($docstring)
            MMI.docstring(::Type{<:$T}) = $docstring
        end
        MMI.target_scitype(::Type{<:$T}) = $target_scitype
        MMI.prediction_type(::Type{<:$T}) = Symbol($pred_str)
        MMI.supports_weights(::Type{<:$T}) = $supports_weights

        # traits specific to measures:
        if !isempty($instances)
            instances(::Type{<:$T}) = $instances
        end
        if !isempty($human_name)
            human_name(::Type{<:$T}) = $human_name
        end
        orientation(::Type{<:$T}) = Symbol($orientation_str)
        reports_each_observation(::Type{<:$T}) = $reports_each_observation
        aggregation(::Type{<:$T}) = $aggregation
        is_feature_dependent(::Type{<:$T}) = $is_feature_dependent
        supports_class_weights(::Type{<:$T}) = $supports_class_weights
        distribution_type(::Type{<:$T}) = $distribution_type

    end
    parentmodule(T).eval(ex)
end

ext: jl


text:
using Test
using GeoParams
@testset "EnergyParameters.jl" begin

# This tests the MaterialParameters structure
CharUnits_GEO   =   GEO_units(viscosity=1e19, length=10km);

# Heat capacity ---------

# Constant heat capacity
cp1      =   ConstantHeatCapacity()
@test cp1.cp.val == 1050J/kg/K

Nondimensionalize!(cp1,CharUnits_GEO)
@test cp1.cp.val ≈ 1.3368075000000002e22

@test ComputeHeatCapacity(100.0, cp1) ≈ 1.3368075000000002e22  # compute

# Temperature-dependent heat capacity
# dimensional
T        =   (250:100:1250)*K;
cp2      =   T_HeatCapacity_Whittacker()
Cp       =   ComputeHeatCapacity(T,cp2)
@test sum(Cp) ≈ 11667.035717418683J/kg/K

# nondimensional
cp2_nd   =   T_HeatCapacity_Whittacker()
Nondimensionalize!(cp2_nd,CharUnits_GEO)
T_nd     =   Float64.(T/CharUnits_GEO.Temperature)
Cp_nd    =   ComputeHeatCapacity(T_nd,cp2_nd)
@test sum(Cp_nd) ≈ 1.4853886523631602e23

# Dimensionalize again and double-check the results
@test ustrip(sum(abs.(Cp_nd*CharUnits_GEO.heatcapacity - Cp))) < 1e-11

# Test with arrays
T_array     =  ustrip.(T)*ones(10)'
Cp_array    =  zeros(size(T_array))
P_array     =  zeros(size(T_array))
ComputeHeatCapacity!(Cp_array, P_array, T_array,cp1)
@test Cp_array[1] ≈ 1.3368075000000002e22

T_array     =  ustrip.(T)*ones(10)'
Cp_array    =  zeros(size(T_array))
P_array     =  zeros(size(T_array))
ComputeHeatCapacity!(Cp_array, P_array, T_array,cp2)
@test sum(Cp_array[:,1]) ≈ 11667.035717418683

T_array     =  ustrip.(T)*ones(10)'
Cp_array    =  zeros(size(T_array))
P_array     =  zeros(size(T_array))
ComputeHeatCapacity!(Cp_array, T_array,cp2)
@test sum(Cp_array[:,1]) ≈ 11667.035717418683

# Check that it works if we give a phase array
MatParam    =   Array{MaterialParams, 1}(undef, 2);
MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=1,
                    HeatCapacity  = ConstantHeatCapacity());

MatParam[2] =   SetMaterialParams(Name="Crust", Phase=2,
                    HeatCapacity  = T_HeatCapacity_Whittacker());

# test computing material properties
n = 100;
Phases              = ones(Int64,n,n,n);
Phases[:,:,20:end] .= 2

Cp = zeros(size(Phases))
T  =  ones(size(Phases))*1500
P  =  zeros(size(Phases))

ComputeHeatCapacity!(Cp, Phases, P, T, MatParam)       # computation routine w/out P (not used in most heat capacity formulations)
@test sum(Cp[1,1,:]) ≈ 121399.0486067196

ComputeHeatCapacity!(Cp, Phases, T, MatParam)       # computation routine w/out P (not used in most heat capacity formulations)
@test sum(Cp[1,1,:]) ≈ 121399.0486067196

# test if we provide phase ratios
PhaseRatio  = zeros(n,n,n,3);
for i in CartesianIndices(Phases)
    iz = Phases[i]
    I = CartesianIndex(i,iz)
    PhaseRatio[I] = 1.0
end
ComputeHeatCapacity!(Cp, PhaseRatio, P, T, MatParam)
@test sum(Cp[1,1,:]) ≈ 121399.0486067196

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

# Conductivity ----------

# Constant

# Constant conductivity
cond      =   ConstantConductivity()
@test Value(cond.k) == 3.0Watt/K/m

Nondimensionalize!(cond,CharUnits_GEO)
@test  Value(cond.k) ≈ 3.8194500000000007

@test ComputeConductivity(100.0, cond) ≈ 3.8194500000000007 # compute

# Temperature-dependent conductivity
# dimensional
T        =   (250:100:1250)*K;
cond2    =   T_Conductivity_Whittacker()
k        =   ComputeConductivity(T,cond2)
@test sum(k) ≈ 27.503366436682285Watt/K/m

# nondimensional
cond2_nd =   T_Conductivity_Whittacker()
Nondimensionalize!(cond2_nd,CharUnits_GEO)
T_nd     =   Float64.(T/CharUnits_GEO.Temperature)
k_nd     =   ComputeConductivity(T_nd,cond2_nd)
@test sum(k_nd) ≈ 35.01591097886205

# Dimensionalize again and double-check the results
@test ustrip(sum(abs.(k_nd*CharUnits_GEO.conductivity - k))) < 1e-11

# Check if we use arrays
T_array     =  ustrip.(T)*ones(100)'
k_array     =  copy(T_array)
P_array     =  copy(T_array)

ComputeConductivity!(k_array,P_array,T_array, cond)
@test k_array[1] ≈ 3.8194500000000007

ComputeConductivity!(k_array,P_array,T_array, cond2)
@test sum(k_array) ≈ 2750.3366436682285

k_TP    =   Set_TP_Conductivity["LowerCrust"]
ComputeConductivity!(k_array, P_array, T_array, k_TP)
@test sum(k_array) ≈ 2055.7129327367625

# Check that it works if we give a phase array
MatParam    =   Array{MaterialParams, 1}(undef, 3);
MatParam[1] =   SetMaterialParams(Name="Mantle", Phase=1,
                    Conductivity  = ConstantConductivity());

MatParam[2] =   SetMaterialParams(Name="Crust", Phase=2,
                    Conductivity  = T_Conductivity_Whittacker());

MatParam[3] =   SetMaterialParams(Name="MantleLithosphere", Phase=3,
                    Conductivity  = Set_TP_Conductivity["Mantle"]);

# test computing material properties
n = 100;
Phases              = ones(Int64,n,n,n);
Phases[:,:,20:end] .= 2
Phases[:,:,60:end] .= 3

PhaseRatio  = zeros(n,n,n,3);
for i in CartesianIndices(Phases)
    iz = Phases[i]
    I = CartesianIndex(i,iz)
    PhaseRatio[I] = 1.0
end

k   =  zeros(size(Phases))
T   =  ones(size(Phases))*1500
P   =  zeros(size(Phases))

ComputeConductivity!(k, Phases, P, T, MatParam)
@test sum(k) ≈ 1.9216938849389635e6

ComputeConductivity!(k, PhaseRatio, P, T, MatParam)
@test sum(k) ≈ 1.9216938849389635e6

#

######

# TP-dependent conductivity for different predefines cases
T       =   (250:100:1250)*K;
P       = 1e6Pa*ones(size(T))
List    = ["LowerCrust"   "Mantle"        "OceanicCrust"  "UpperCrust"]
Sol_kT  = [20.55712932736763 28.700405819019323 20.55712932736763 19.940302462417037]*Watt/K/m
for i=1:length(List)
    k_TP    =   Set_TP_Conductivity[List[i]]
    k       =   ComputeConductivity(P,T,k_TP)
    @test sum(k) ≈ Sol_kT[i]

    k_TP_nd  =  deepcopy(k_TP)
    Nondimensionalize!(k_TP_nd,CharUnits_GEO)
    T_nd     =   Float64.(T/CharUnits_GEO.Temperature)
    P_nd     =   Float64.(P/CharUnits_GEO.stress)
    k_nd     =   ComputeConductivity(P_nd,T_nd,k_TP_nd)

    @test ustrip(sum(abs.(k_nd*CharUnits_GEO.conductivity - k))) < 1e-11

end

T = [200 300; 400 500]*K
k1        =   ComputeConductivity(T,cond2)

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

# Latent heat -----------
a = ConstantLatentHeat()
Q_L = ComputeLatentHeat(a)
@test Q_L == 400kJ/kg

Nondimensionalize!(a,CharUnits_GEO)
Q_L = ComputeLatentHeat(a)
@test Q_L ≈ 4e21
# -----------------------

# Radioactive heat ------
a = ConstantRadioactiveHeat()
H_r = ComputeRadioactiveHeat(a)
@test H_r ≈ 1.0e-6Watt/m^3

Nondimensionalize!(a,CharUnits_GEO)
H_r = ComputeRadioactiveHeat(a)
@test H_r == 0.1
# -----------------------

# Shear heating -------
Χ       = ConstantShearheating(1.0)

# Define parameters as vectors
τ       = [1 2 3 4]*1e6*Pa
ε       = [1 0.1 0.1 1]*1/s
ε_el    = [0.01 0.01 0.01 0.01]*1/s

τ_2D       = [1 2; 3 4]*1e6*Pa
ε_2D       = [1 0.1; 0.1 1]*1/s
ε_el_2D    = [0.01 0.01; 0.01 0.01]*1/s

# With elasticity
H_s1 = ComputeShearheating(τ,   ε,    ε_el,     Χ)
H_s2 = ComputeShearheating(τ_2D,ε_2D, ε_el_2D,  Χ)
@test H_s1 ≈ 5.4e6Pa/s
@test H_s2 ≈ 5.4e6Pa/s

# No elasticity
H_s3 = ComputeShearheating(τ,   ε,     Χ)
H_s4 = ComputeShearheating(τ_2D,ε_2D,  Χ)
@test H_s3 ≈ 5.5e6Pa/s
@test H_s4 ≈ 5.5e6Pa/s

# Now in non-dimensional units
τ       = [1 2 3 4]
ε       = [1 0.1 0.1 1]
ε_el    = [0.01 0.01 0.01 0.01]

τ_2D       = [1 2; 3 4]
ε_2D       = [1 0.1; 0.1 1]
ε_el_2D    = [0.01 0.01; 0.01 0.01]
Nondimensionalize!(Χ,CharUnits_GEO)

H_s1 = ComputeShearheating(τ,   ε,    ε_el,     Χ)
H_s2 = ComputeShearheating(τ_2D,ε_2D, ε_el_2D,  Χ)
H_s3 = ComputeShearheating(τ,   ε,     Χ)
H_s4 = ComputeShearheating(τ_2D,ε_2D,  Χ)
@test H_s1 ≈ 5.4
@test H_s2 ≈ 5.4
@test H_s3 ≈ 5.5
@test H_s4 ≈ 5.5
# -----------------------

end


text:
using Manifolds, Manopt, ManifoldsBase, Test

@testset "Plan Tests         " begin
    include("plans/test_options.jl")
    include("plans/test_debug.jl")
    include("plans/test_nelder_mead_plan.jl")
    include("plans/test_nonmutating.jl")
    include("plans/test_gradient_plan.jl")
    include("plans/test_primal_dual_plan.jl")
    include("plans/test_record.jl")
    include("plans/test_stopping_criteria.jl")
end
@testset "Function Tests     " begin
    include("functions/test_adjoint_differentials.jl")
    include("functions/test_bezier.jl")
    include("functions/test_differentials.jl")
    include("functions/test_costs.jl")
    include("functions/test_gradients.jl")
    include("functions/test_proximal_maps.jl")
    include("functions/test_manifold.jl")
end
@testset "Helper & Data Tests" begin
    include("helpers/test_error_measures.jl")
    include("helpers/test_data.jl")
end
@testset "Solver Tests       " begin
    include("solvers/test_alternating_gradient.jl")
    include("solvers/test_ChambollePock.jl")
    include("solvers/test_conjugate_gradient.jl")
    include("solvers/test_Douglas_Rachford.jl")
    include("solvers/test_cyclic_proximal_point.jl")
    include("solvers/test_gradient_descent.jl")
    include("solvers/test_Nelder_Mead.jl")
    include("solvers/test_quasi_Newton.jl")
    include("solvers/test_trust_regions.jl")
    include("solvers/test_particle_swarm.jl")
    include("solvers/test_stochastic_gradient_descent.jl")
    include("solvers/test_subgradient_method.jl")
end


text:
using BioformatsLoader
using JavaCall
using Test

@test BioformatsLoader.init() == nothing

oxr = OMEXMLReader()
@test !JavaCall.isnull(oxr.reader)
@test !JavaCall.isnull(oxr.meta)

nd2_file = bf_import(joinpath(@__DIR__,"test_format.nd2"))
czi_file = bf_import(joinpath(@__DIR__,"test_zen.CZI"))


text:
module FC

include("Format.jl")
using .Format

export fc

# Determines upper bound alpha with Floor-Ceiling Theorem
# Authors: Antara Hebbar and Jason Liu
function fc(m::Int64, s::Int64; output::Int64 = 1)
    output > 0 && printHeader(center("FLOOR-CEILING THEOREM"))

    if m <= 0 || s <= 0
        if output > 0
            printf("Floor-Ceiling Theorem does not apply", line = true)
            printEnd()
        end
        1 // 1
    elseif m < s
        alpha = m // s * fc(s, m, output = 0)
        if output > 0
            printfT(
                "Floor-Ceiling Theorem",
                "Upper bound α of muffins($m,$s) is $(formatFrac(alpha))",
            )
            printEnd()
        end
        alpha
    else
        alpha = m % s == 0 ? 1 // 1 :
            max(1 // 3, min(m // Int64(ceil(2m / s)s), 1 - m // Int64(floor(2m / s)s)))
        if output > 0
            printfT(
                "Floor-Ceiling Theorem",
                "Upper bound α of muffins($m,$s) is $(formatFrac(alpha))",
            )
            printEnd()
        end
        alpha
    end
end

end


text:
# ....................
# @license: %MIT License%:~ http://www.opensource.org/licenses/MIT
# @project: lr_with_outlier
# @file: /main.jl
# @created: Sunday, 5th April 2020
# @author: brentian (chuwzhang@126.com)
# @modified: brentian (chuwzhang@126.com>)
#    Sunday, 5th April 2020 4:14:00 pm
# ....................
# @description:
# real data
module MyApp

include("./linear_regression.jl")
using .lr_with_outlier
using Plots
# define size := m, n

Base.@ccallable function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    # Preparing an optimization model
    # === 2-d
    # n, m = 10, 2
    # x = [ ones(n) [i * 1.0 for i = 1:n]]
    # y = Array(2:n + 1)
    # y[3] = 100
    # === 1-d
    # n, m = 10, 1
    # x =  [i * 1.0 for i = 1:n]
    # y = Array(2:n + 1)
    # y[3] = 100
    y = [
        11.19
        10.52
        10.73
        10.26
        11.27
        11.54
        10.56
        10.48
        10.5
        10.34
        10.29
        10.43
        10.18
        10.43
        10.32
        10.14
        10.79
    ]

    x = [
        0.596
        0.604
        0.589
        0.589
        0.598
        0.595
        0.599
        0.601
        0.592
        0.600
        0.592
        0.587
        0.589
        0.588
        0.588
        0.586
        0.594
    ]

    # normalize
    x_min, x_max = minimum(x), maximum(x)
    x = (x .- x_min) ./ (x_max - x_min)
    y_min, y_max = minimum(y), maximum(y)
    y0 = (y .- y_min) ./ (y_max - y_min)
    n, m = size(x, 1), ndims(x)

    Plots.plotly()
    results = [lr_with_outlier.f(x, y0, n, m, i, x_min, x_max, y_min, y_max) for i in 1:2:n]
    plots = [i[1] for i in results]
    agg_p = Plots.plot(plots..., wsize = (1200, 900), show = false)
    # Plots.savefig(agg_p, "output_filename.html")
    # Plots.savefig(agg_p, "output_filename.png")
    Plots.savefig(agg_p, "figure.json")
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end
end

max_stars_repo_name: brentian/artoopt

max_issues_repo_name: brentian/artoopt

max_forks_repo_name: brentian/artoopt


text:
module Test_ldlt

include("../src/ldlt.jl")

using Test
using Random
using SparseArrays
using LinearAlgebra

@testset "ldlt" begin
    @testset "symbolic_ldlt" begin
        @testset "n = 1" begin
            Ap = [1,2]
            Ai = [1]
            Fp,Fi = symbolic_ldlt(Ap,Ai)
            @test Fp == [1,2]
            @test Fi == [1]
        end

        @testset "n = 2" begin
            @testset "A = [1 0; 0 1]" begin
                Ap = [1,2,3]
                Ai = [1,2]
                Fp,Fi = symbolic_ldlt(Ap,Ai)
                @test Fp == [1,2,3]
                @test Fi == [1,2]
            end
            @testset "A = [1 1; 1 0]" begin
                Ap = [1,3,3]
                Ai = [1,2]
                Fp,Fi = symbolic_ldlt(Ap,Ai)
                @test Fp == [1,3,4]
                @test Fi == [1,2,2]
            end
            @testset "A = [1 1; 1 1]" begin
                Ap = [1,3,4]
                Ai = [1,2,2]
                Fp,Fi = symbolic_ldlt(Ap,Ai)
                @test Fp == [1,3,4]
                @test Fi == [1,2,2]
            end
        end

        @testset "n in 3:100" begin
            @testset for seed = 1:100
                Random.seed!(seed)
                n = rand(3:100)
                A = sprand(n,n,rand())
                A = 2I + A + A'
                Ap,Ai = A.colptr,A.rowval
                Fref = lu(Matrix(A),Val(false)).L .!= 0
                Fp,Fi = symbolic_ldlt(Ap,Ai)
                F = SparseMatrixCSC(n,n,Fp,Fi,ones(Bool,length(Fi)))
                @test F == Fref
            end
        end
    end

    @testset "symbolic_ildlt_lof" begin
        @testset "c = 0" begin
            @testset for seed = 1:100
                Random.seed!(seed)
                n = rand(3:100)
                A = sprand(n,n,rand())
                A = tril(2I + A + A')
                Ap,Ai = A.colptr,A.rowval

                Fp,Fi,Fl = symbolic_ildlt_lof(Ap,Ai,0)
                @test (Fp,Fi) == (Ap,Ai)
            end
        end

        @testset "c = n" begin
            @testset for seed = 1:100
                Random.seed!(seed)
                n = rand(3:100)
                A = sprand(n,n,rand())
                A = 2I + A + A'
                Ap,Ai = A.colptr,A.rowval

                rFp,rFi = symbolic_ldlt(Ap,Ai)
                Fp,Fi,Fl = symbolic_ildlt_lof(Ap,Ai,n)
                @test (Fp,Fi) == (rFp,rFi)
            end
        end

        @testset "self-consistency" begin
            @testset for seed = 1:100
                Random.seed!(seed)
                n = rand(3:100)
                A = sprand(n,n,rand())
                A = 2I + A + A'
                Ap,Ai = A.colptr,A.rowval

                Fref = lu(Matrix(A),Val(false)).L .!= 0
                Fp,Fi,Fl = symbolic_ildlt_lof(Ap,Ai,n)
                F = SparseMatrixCSC(n,n,Fp,Fi,Fl)

                cmax = maximum(Fl)
                for c = 0:cmax
                    ref = SparseArrays.fkeep!(copy(F),(i,j,l)->(l<=c))
                    rFp,rFi,rFl = ref.colptr,ref.rowval,ref.nzval
                    iFp,iFi,iFl = symbolic_ildlt_lof(Ap,Ai,c)
                    @test (iFp,iFi,iFl) == (rFp,rFi,rFl)
                end
            end
        end
    end

    @testset "numeric_ldlt" begin
        function test(A,conj,adj)
            Ap,Ai,Av = A.colptr,A.rowval,A.nzval
            Fp,Fi = symbolic_ldlt(Ap,Ai)
            Fv = numeric_ldlt(Ap,Ai,Av,Fp,Fi; conj=conj)
            F = SparseMatrixCSC(size(A)...,Fp,Fi,Fv)

            L = tril(F,-1) + I; D = Diagonal(F);
            @test eltype(F) == eltype(A)
            @test L*D*adj(L) ≈ A
        end

        @testset for T in (Float32,Float64,ComplexF32,ComplexF64)
            @testset for (conj,adj) in ((conj,adjoint),(identity,transpose))
                @testset "n = 1" begin
                    Random.seed!(1)
                    Ap = [1,2]
                    Ai = [1]
                    Av = rand(T,1)
                    A = SparseMatrixCSC(1,1,Ap,Ai,Av)
                    A = 2I + A + adj(A)
                    test(A,conj,adj)
                end

                @testset "n = 2" begin
                    @testset "A = [1 0; 0 1]" begin
                        Random.seed!(1)
                        Ap = [1,2,3]
                        Ai = [1,2]
                        Av = rand(T,2)
                        A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                    @testset "A = [1 1; 1 0]" begin
                        Ap = [1,3,3]
                        Ai = [1,2]
                        Av = rand(T,2)
                        A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                    @testset "A = [1 1; 1 1]" begin
                        Ap = [1,3,4]
                        Ai = [1,2,2]
                        Av = rand(T,3)
                        A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                end

                @testset "n in 3:100" begin
                    @testset for seed = 1:100
                        Random.seed!(seed)
                        n = rand(3:100)
                        A = sprand(n,n,rand())
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                end
            end
        end
    end

    @testset "ildlt_tol" begin
        function test(A,conj,adj)
            A = tril(A)
            Ap,Ai,Av = A.colptr,A.rowval,A.nzval

            # Test that factorization is exact if tol = 0
            Fp,Fi,Fv = ildlt_tol(Ap,Ai,Av, 0.0; conj=conj)
            F = SparseMatrixCSC(size(A)...,Fp,Fi,Fv)
            L = tril(F,-1) + I; D = Diagonal(F);
            @test eltype(F) == eltype(A)
            @test tril(L*D*adj(L)) ≈ A

            # Test that factorization produces no fill-in for tol = Inf
            Fp,Fi,Fv = ildlt_tol(Ap,Ai,Av, Inf; conj=conj)
            @test (Fp,Fi) == (Ap,Ai)
        end

        @testset for T in (Float32,Float64,ComplexF32,ComplexF64)
            @testset for (conj,adj) in ((conj,adjoint),(identity,transpose))
                @testset "n = 1" begin
                    Random.seed!(1)
                    Ap = [1,2]
                    Ai = [1]
                    Av = rand(T,1)
                    A = SparseMatrixCSC(1,1,Ap,Ai,Av)
                    A = 2I + A + adj(A)
                    test(A,conj,adj)
                end

                @testset "n = 2" begin
                    @testset "A = [1 0; 0 1]" begin
                        Random.seed!(1)
                        Ap = [1,2,3]
                        Ai = [1,2]
                        Av = rand(T,2)
                        A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                    @testset "A = [1 1; 1 1]" begin
                        Ap = [1,3,4]
                        Ai = [1,2,2]
                        Av = rand(T,3)
                        A = SparseMatrixCSC(2,2,Ap,Ai,Av)
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                end

                @testset "n in 3:100" begin
                    @testset for seed = 1:100
                        Random.seed!(seed)
                        n = rand(3:100)
                        A = sprand(n,n,rand())
                        A = 2I + A + adj(A)
                        test(A,conj,adj)
                    end
                end
            end
        end
    end
end

end

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
@testset "orcv2" begin

@testset "ThreadSafeModule" begin
    ts_mod = ThreadSafeModule("jit")
    @test_throws LLVMException ts_mod() do mod
        error("Error")
    end
    run = Ref{Bool}(false)
    ts_mod() do mod
        run[] = true
    end
    @test run[]
end

@testset "JITDylib" begin
    LLJIT() do lljit
        es = ExecutionSession(lljit)

        @test LLVM.lookup_dylib(es, "my.so") === nothing

        jd = JITDylib(es, "my.so")
        jd_bare = JITDylib(es, "mybare.so", bare=true)

        @test LLVM.lookup_dylib(es, "my.so") === jd

        jd_main = JITDylib(lljit)

        prefix = LLVM.get_prefix(lljit)
        dg = LLVM.CreateDynamicLibrarySearchGeneratorForProcess(prefix)
        add!(jd_main, dg)

        addr = lookup(lljit, "jl_apply_generic")
        @test pointer(addr) != C_NULL
    end
end

@testset "Undefined Symbol" begin
    LLJIT() do lljit
        @test_throws LLVMException lookup(lljit, string(gensym()))
    end

    LLJIT(;tm=JITTargetMachine()) do lljit
        jd = JITDylib(lljit)

        ThreadSafeContext() do ts_ctx
            ctx = context(ts_ctx)
            mod = LLVM.Module("jit"; ctx)

            T_Int32 = LLVM.Int32Type(ctx)
            ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])
            fn = LLVM.Function(mod, "mysum", ft)
            linkage!(fn, LLVM.API.LLVMExternalLinkage)

            fname = "wrapper"
            wrapper = LLVM.Function(mod, fname, ft)
            # generate IR
            Builder(ctx) do builder
                entry = BasicBlock(wrapper, "entry"; ctx)
                position!(builder, entry)

                tmp = call!(builder, fn, [parameters(wrapper)...])
                ret!(builder, tmp)
            end

            # TODO: Get TM from lljit?
            tm = JITTargetMachine()
            triple!(mod, triple(lljit))
            ModulePassManager() do pm
                add_library_info!(pm, triple(mod))
                add_transform_info!(pm, tm)
                run!(pm, mod)
            end
            verify(mod)

            ts_mod = ThreadSafeModule(mod; ctx=ts_ctx)

            add!(lljit, jd, ts_mod)
            @test_throws LLVMException lookup(lljit, fname)
        end
    end
end

@testset "Loading ObjectFile" begin
    LLJIT(;tm=JITTargetMachine()) do lljit
        jd = JITDylib(lljit)

        sym = "SomeFunction"
        obj = ThreadSafeContext() do ts_ctx
            ctx = context(ts_ctx)

            mod = LLVM.Module("jit"; ctx)
            ft = LLVM.FunctionType(LLVM.VoidType(ctx))
            fn = LLVM.Function(mod, sym, ft)

            Builder(ctx) do builder
                entry = BasicBlock(fn, "entry"; ctx)
                position!(builder, entry)
                ret!(builder)
            end
            verify(mod)

            tm  = JITTargetMachine()
            return emit(tm, mod, LLVM.API.LLVMObjectFile)
        end
        add!(lljit, jd, MemoryBuffer(obj))

        addr = lookup(lljit, sym)

        @test pointer(addr) != C_NULL

        empty!(jd)
        @test_throws LLVMException lookup(lljit, sym)
    end

    LLJIT(;tm=JITTargetMachine()) do lljit
        jd = JITDylib(lljit)

        sym = "SomeFunction"
        obj = ThreadSafeContext() do ts_ctx
            ctx = context(ts_ctx)

            mod = LLVM.Module("jit"; ctx)
            ft = LLVM.FunctionType(LLVM.Int32Type(ctx))
            fn = LLVM.Function(mod, sym, ft)

            gv = LLVM.GlobalVariable(mod, LLVM.Int32Type(ctx), "gv")
            LLVM.extinit!(gv, true)

            Builder(ctx) do builder
                entry = BasicBlock(fn, "entry"; ctx)
                position!(builder, entry)
                val = load!(builder, gv)
                ret!(builder, val)
            end
            verify(mod)

            tm  = JITTargetMachine()
            return emit(tm, mod, LLVM.API.LLVMObjectFile)
        end

        data = Ref{Int32}(42)
        GC.@preserve data begin
            address = LLVM.API.LLVMOrcJITTargetAddress(
                reinterpret(UInt, Base.unsafe_convert(Ptr{Int32}, data)))
            flags = LLVM.API.LLVMJITSymbolFlags(
                LLVM.API.LLVMJITSymbolGenericFlagsExported, 0)
            symbol = LLVM.API.LLVMJITEvaluatedSymbol(
                address, flags)
            gv = LLVM.API.LLVMJITCSymbolMapPair(
                    mangle(lljit, "gv"), symbol)

            mu = LLVM.absolute_symbols(Ref(gv))
            LLVM.define(jd, mu)

            add!(lljit, jd, MemoryBuffer(obj))

            addr = lookup(lljit, sym)
            @test pointer(addr) != C_NULL

            @test ccall(pointer(addr), Int32, ()) == 42
            data[] = -1
            @test ccall(pointer(addr), Int32, ()) == -1
        end
        empty!(jd)
        @test_throws LLVMException lookup(lljit, sym)
    end
end

@testset "ObjectLinkingLayer" begin
    called_oll = Ref{Int}(0)

    ollc = LLVM.ObjectLinkingLayerCreator() do es, triple
        oll = ObjectLinkingLayer(es)
        register!(oll, GDBRegistrationListener())
        called_oll[] += 1
        return oll
    end

    GC.@preserve ollc begin
        builder = LLJITBuilder()
        linkinglayercreator!(builder, ollc)
        LLJIT(builder) do lljit
            jd = JITDylib(lljit)

            ThreadSafeContext() do ts_ctx
                ctx = context(ts_ctx)
                sym = "SomeFunctionOLL"

                mod = LLVM.Module("jit"; ctx)
                ft = LLVM.FunctionType(LLVM.VoidType(ctx))
                fn = LLVM.Function(mod, sym, ft)

                Builder(ctx) do builder
                    entry = BasicBlock(fn, "entry"; ctx)
                    position!(builder, entry)
                    ret!(builder)
                end
                verify(mod)

                ts_mod = ThreadSafeModule(mod; ctx=ts_ctx)
                add!(lljit, jd, ts_mod)

                addr = lookup(lljit, sym)

                @test pointer(addr) != C_NULL
            end
        end
    end
    @test called_oll[] >= 1
end

@testset "Lazy" begin
    LLJIT() do lljit
        jd = JITDylib(lljit)
        es = ExecutionSession(lljit)

        lctm = LLVM.LocalLazyCallThroughManager(triple(lljit), es)
        ism = LLVM.LocalIndirectStubsManager(triple(lljit))
        try
            # 1. define entry symbol
            entry_sym = "foo_entry"
            flags = LLVM.API.LLVMJITSymbolFlags(
                LLVM.API.LLVMJITSymbolGenericFlagsCallable |
                LLVM.API.LLVMJITSymbolGenericFlagsExported, 0)
            entry = LLVM.API.LLVMOrcCSymbolAliasMapPair(
                mangle(lljit, entry_sym),
                LLVM.API.LLVMOrcCSymbolAliasMapEntry(
                    mangle(lljit, "foo"), flags))

            mu = LLVM.reexports(lctm, ism, jd, Ref(entry))
            LLVM.define(jd, mu)

            # 2. Lookup address of entry symbol
            addr = lookup(lljit, entry_sym)
            @test pointer(addr) != C_NULL

            # 3. add MU that will call back into the compiler
            sym = LLVM.API.LLVMOrcCSymbolFlagsMapPair(mangle(lljit, "foo"), flags)

            function materialize(mr)
                syms = LLVM.get_requested_symbols(mr)
                @assert length(syms) == 1

                # syms contains mangled symbols
                # we need to emit an unmangled one

                TSM = ThreadSafeModule("jit")
                TSM() do mod
                    LLVM.apply_datalayout!(lljit, mod)

                    ctx = context(mod)
                    T_Int32 = LLVM.Int32Type(ctx)
                    ft = LLVM.FunctionType(T_Int32, [T_Int32, T_Int32])

                    fn = LLVM.Function(mod, "foo", ft)

                    # generate IR
                    Builder(ctx) do builder
                        entry = BasicBlock(fn, "entry"; ctx)
                        position!(builder, entry)

                        tmp = add!(builder, parameters(fn)...)
                        ret!(builder, tmp)
                    end
                end

                il = LLVM.IRTransformLayer(lljit)
                LLVM.emit(il, mr, TSM)

                return nothing
            end

            function discard(jd, sym)
            end

            mu = LLVM.CustomMaterializationUnit("fooMU", Ref(sym), materialize, discard)
            LLVM.define(jd, mu)

            @test ccall(pointer(addr), Int32, (Int32, Int32), 1, 2) == 3
        finally
            dispose(lctm)
            dispose(ism)
        end
    end
end

end


text:
export GitDiff, parse_git_diff_options, deltas, patches, diff_workdir, patch

type GitDiff
    ptr::Ptr{Void}

    function GitDiff(ptr::Ptr{Void})
        @assert ptr != C_NULL
        d = new(ptr)
        finalizer(d, free!)
        return d
    end
end

free!(d::GitDiff) = begin
    if d.ptr != C_NULL
        ccall((:git_diff_free, libgit2), Void, (Ptr{Void},), d.ptr)
        d.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, d::GitDiff) = d.ptr

type DiffStats
    files::Int
    adds::Int
    dels::Int
end
DiffStats() = DiffStats(0, 0, 0)

function cb_diff_file_stats(delta_ptr::Ptr{DiffDeltaStruct}, progress::Cfloat, payload::Ptr{Void})
    delta = unsafe_load(delta_ptr)::DiffDeltaStruct
    stats = unsafe_pointer_to_objref(payload)::DiffStats
    if delta.status == GitConst.DELTA_ADDED ||
       delta.status == GitConst.DELTA_DELETED ||
       delta.status == GitConst.DELTA_MODIFIED ||
       delta.status == GitConst.DELTA_RENAMED ||
       delta.status == GitConst.DELTA_COPIED ||
       delta.status == GitConst.DELTA_TYPECHANGE
        stats.files += 1
    end
    return GitErrorConst.GIT_OK
end

const c_cb_diff_file_stats = cfunction(cb_diff_file_stats, Cint,
                                       (Ptr{DiffDeltaStruct}, Cfloat, Ptr{Void}))

function cb_diff_line_stats(delta_ptr::Ptr{Void},
                            hunk_ptr::Ptr{Void},
                            line_ptr::Ptr{DiffLineStruct},
                            payload::Ptr{Void})
    line  = unsafe_load(line_ptr)
    stats = unsafe_pointer_to_objref(payload)::DiffStats
    if line.origin == GitConst.DIFF_LINE_ADDITION
        stats.adds += 1
    elseif line.origin == GitConst.DIFF_LINE_DELETION
        stats.dels += 1
    end
    return GitErrorConst.GIT_OK
end

const c_cb_diff_line_stats = cfunction(cb_diff_line_stats, Cint,
                                       (Ptr{Void}, Ptr{Void}, Ptr{DiffLineStruct}, Ptr{Void}))
Base.stat(d::GitDiff) = begin
    stats = DiffStats()
    ccall((:git_diff_foreach, libgit2), Void,
          (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Any),
          d, c_cb_diff_file_stats, C_NULL, c_cb_diff_line_stats, &stats)
    return stats
end

type DiffFile
    oid::Oid
    path::ByteString
    size::Int
    flags::Int
    mode::Int
end

DiffFile(struct::DiffFileStruct) = begin
    DiffFile(Oid(string(struct.id)),
             bytestring(struct.path),
             struct.size,
             struct.flags,
             struct.mode)
end

function delta_status_symbol(s::Integer)
    s == GitConst.DELTA_UNMODIFIED && return :unmodified
    s == GitConst.DELTA_ADDED      && return :added
    s == GitConst.DELTA_DELETED    && return :deleted
    s == GitConst.DELTA_MODIFIED   && return :modified
    s == GitConst.DELTA_RENAMED    && return :renamed
    s == GitConst.DELTA_COPIED     && return :copied
    s == GitConst.DELTA_IGNORED    && return :ignored
    s == GitConst.DELTA_UNTRACKED  && return :untracked
    s == GitConst.DELTA_TYPECHANGE && return :typechange
    throw(ArgumentError("Unknown status symbol :$s"))
end

type DiffDelta
    old_file::DiffFile
    new_file::DiffFile
    similarity::Int
    status::Symbol
    isbinary::Bool

    function DiffDelta(ptr::Ptr{DiffDeltaStruct})
        @assert ptr != C_NULL
        d = unsafe_load(ptr)::DiffDeltaStruct
        old_file_id = Oid(hex(d.old_file.id))
        fold = DiffFile(old_file_id,
                        d.old_file.path != C_NULL ? bytestring(d.old_file.path) : "",
                        int(d.old_file.size),
                        int(d.old_file.flags),
                        int(d.old_file.mode))
        new_file_id = Oid(hex(d.new_file.id))
        fnew = DiffFile(new_file_id,
                        d.new_file.path != C_NULL ? bytestring(d.new_file.path) : "",
                        int(d.new_file.size),
                        int(d.new_file.flags),
                        int(d.new_file.mode))
        return new(fold, fnew, int(d.similarity), delta_status_symbol(d.status),
                   (bool(d.flags & GitConst.DIFF_FLAG_NOT_BINARY) &&
                    bool(d.flags & GitConst.DIFF_FLAG_BINARY)))
    end
end

Base.length(d::GitDiff) = int(ccall((:git_diff_num_deltas, libgit2), Csize_t, (Ptr{Void},), d))

function deltas(d::GitDiff)
    ndelta = length(d)
    if ndelta == 0
        return nothing
    end
    ds = Array(DiffDelta, ndelta)
    for i in 1:ndelta
        delta_ptr = ccall((:git_diff_get_delta, libgit2), Ptr{DiffDeltaStruct},
                          (Ptr{Void}, Csize_t), d, i-1)
        ds[i] = DiffDelta(delta_ptr)
    end
    return ds
end

function patches(d::GitDiff)
    ndelta = length(d)
    if ndelta == 0
        return nothing
    end
    ps = GitPatch[]
    patch_ptr = Ptr{Void}[0]
    err = zero(Cint)
    for i in 1:ndelta
        err = ccall((:git_patch_from_diff, libgit2), Cint,
                    (Ptr{Ptr{Void}}, Ptr{Void}, Csize_t), patch_ptr, d, i-1)
        if err != GitErrorConst.GIT_OK
            break
        end
        push!(ps, GitPatch(patch_ptr[1]))
    end
    if err != GitErrorConst.GIT_OK
        throw(GitError(err))
    end
    return ps
end

#TODO: memory leaks?
function cb_diff_print(delta_ptr::Ptr{Void},
                       hunk_ptr::Ptr{Void},
                       line_ptr::Ptr{DiffLineStruct},
                       payload::Ptr{Void})
    l = unsafe_load(line_ptr)
    s = unsafe_pointer_to_objref(payload)::Array{Uint8,1}
    add_origin = false
    if l.origin == GitConst.DIFF_LINE_CONTEXT ||
       l.origin == GitConst.DIFF_LINE_ADDITION ||
       l.origin == GitConst.DIFF_LINE_DELETION
       add_origin = true
    end
    prev_len = length(s)
    if add_origin
        resize!(s, prev_len + l.content_len + 1)
        s[prev_len + 1] = l.origin
        for i in 1:l.content_len
            s[prev_len + i + 1] = unsafe_load(l.content, i)
        end
    else
        resize!(s, prev_len + l.content_len)
        for i in 1:l.content_len
            s[prev_len + i] = unsafe_load(l.content, i)
        end
    end
    return GitErrorConst.GIT_OK
end

const c_cb_diff_print = cfunction(cb_diff_print, Cint, (Ptr{Void}, Ptr{Void}, Ptr{DiffLineStruct}, Ptr{Void}))

function patch(d::GitDiff; format::Symbol=:patch)
    cformat = zero(Cint)
    if format === :name_status
        cformat = GitConst.DIFF_FORMAT_NAME_STATUS
    elseif format === :name_only
        cformat = GitConst.DIFF_FORMAT_NAME_ONLY
    elseif format === :raw
        cformat = GitConst.DIFF_FORMAT_RAW
    elseif format === :header
        cformat = GitConst.DIFF_FORMAT_PATCH_HEADER
    elseif format === :patch
        cformat = GitConst.DIFF_FORMAT_PATCH
    else
        throw(ArgumentError(("Unknown diff output format ($format)")))
    end
    s = Uint8[]
    ccall((:git_diff_print, libgit2), Cint,
           (Ptr{Void}, Cint, Ptr{Void}, Any), d, cformat, c_cb_diff_print, &s)
    return UTF8String(s)
end

# diffable GitTree, GitCommit, GitIndex, or Nothing
typealias Diffable Union(GitTree, GitCommit, GitIndex, Nothing)

Base.diff(repo::GitRepo, left::Nothing, right::Nothing, opts::MaybeDict=nothing) = nothing

Base.diff(repo::GitRepo, left::MaybeString, right::MaybeString, opts::MaybeDict=nothing) = begin
    l = left  != nothing ? revparse(repo, left)  : nothing
    r = right != nothing ? revparse(repo, right) : nothing
    if l != nothing
        return diff(repo, l, r, opts)
    elseif r != nothing
        opts = opts == nothing ? {} : opts
        return diff(repo, l, r,
                merge(opts, {:reverse => !get(opts, :reverse, true)}))
    end
    return nothing
end

Base.diff(repo::GitRepo, left::GitCommit, right::GitCommit, opts::MaybeDict=nothing) = begin
    return diff(repo, GitTree(left), GitTree(right), opts)
end

Base.diff(repo::GitRepo, left::GitCommit, right::Nothing, opts::MaybeDict=nothing) = begin
    return diff(repo, GitTree(left), nothing, opts)
end

Base.diff(repo::GitRepo, left::Nothing, right::GitCommit, opts::MaybeDict=nothing) = begin
    return diff(repo, nothing, GitTree(right), opts)
end

Base.diff(repo::GitRepo, c::GitCommit, opts::MaybeDict=nothing) = begin
    ps = parents(c)
    if length(ps) > 0
        return diff(repo, GitTree(c), GitTree(ps[1]), opts)
    else
        return diff(repo, GitTree(c), nothing, opts)
    end
end

Base.diff(repo::GitRepo, left::GitTree, right::String, opts::MaybeDict=nothing) = begin
    return diff(repo, left, revparse(repo, right), opts)
end

Base.diff(repo::GitRepo, left::GitTree, right::GitCommit, opts::MaybeDict=nothing) = begin
    return diff(repo, left, GitTree(right), opts)
end

typealias MaybeGitTree Union(Nothing, GitTree)

Base.diff(repo::GitRepo, left::MaybeGitTree, right::MaybeGitTree, opts::MaybeDict=nothing) = begin
    gopts = parse_git_diff_options(opts)
    diff_ptr = Ptr{Void}[0]
    @check ccall((:git_diff_tree_to_tree, libgit2), Cint,
                 (Ptr{Ptr{Void}},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{DiffOptionsStruct}),
                 diff_ptr,
                 repo,
                 left  != nothing ? left.ptr : C_NULL,
                 right != nothing ? right.ptr : C_NULL,
                 &gopts)
    free!(gopts.pathspec)
    return GitDiff(diff_ptr[1])
end

Base.diff(repo::GitRepo, left::Oid, right::Oid, opts::MaybeDict=nothing) = begin
    return diff(repo, lookup(repo, left), lookup(repo, right), opts)
end

Base.diff(repo::GitRepo, left::String, right::GitIndex, opts::MaybeDict=nothing) = begin
    return diff(repo, revparse(repo, left), right, opts)
end

Base.diff(repo::GitRepo, left::GitCommit, right::GitIndex, opts::MaybeDict=nothing) = begin
    return diff(repo, GitTree(left), right, opts)
end

Base.diff(repo::GitRepo, left::GitTree, right::GitIndex, opts::MaybeDict=nothing) = begin
    gopts = parse_git_diff_options(opts)
    diff_ptr = Ptr{Void}[0]
    @check ccall((:git_diff_tree_to_index, libgit2), Cint,
                 (Ptr{Ptr{Void}},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{DiffOptionsStruct}),
                 diff_ptr, repo, left, right, &gopts)
    free!(gopts.pathspec)
    return GitDiff(diff_ptr[1])
end

Base.diff(repo::GitRepo, idx::GitIndex, opts::MaybeDict=nothing) = begin
    return diff(repo, idx, nothing, opts)
end

Base.diff(repo::GitRepo, idx::GitIndex, other::Nothing, opts::MaybeDict=nothing) = begin
    gopts = parse_git_diff_options(opts)
    diff_ptr = Ptr{Void}[0]
    @check ccall((:git_diff_index_to_workdir, libgit2), Cint,
                  (Ptr{Ptr{Void}},
                   Ptr{Void},
                   Ptr{Void},
                   Ptr{DiffOptionsStruct}),
                  diff_ptr, repo, idx, &gopts)
    free!(gopts.pathspec)
    return GitDiff(diff_ptr[1])
end

Base.diff(repo::GitRepo, idx::GitIndex, other::GitCommit, opts::MaybeDict=nothing) = begin
    return diff(repo, idx, GitTree(other), opts)
end

Base.diff(repo::GitRepo, idx::GitIndex, other::GitTree, opts::MaybeDict=nothing) = begin
    gopts = parse_git_diff_options(opts)
    diff_ptr = Ptr{Void}[0]
    @check ccall((:git_diff_tree_to_index, libgit2), Cint,
                 (Ptr{Ptr{Void}},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{DiffOptionsStruct}),
                  diff_ptr, repo, other, idx, &gopts)
    free!(gopts.pathspec)
    return GitDiff(diff_ptr[1])
end

Base.merge!(d1::GitDiff, d2::GitDiff) = begin
    @check ccall((:git_diff_merge, libgit2), Cint, (Ptr{Void}, Ptr{Void}), d1, d2)
    return d1
end

function diff_workdir(repo::GitRepo, opts::MaybeDict=nothing)
    gopts = parse_git_diff_options(opts)
    diff_ptr = Ptr{Void}[0]
    @check ccall((:git_diff_index_to_workdir, libgit2), Cint,
                 (Ptr{Ptr{Void}},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{DiffOptionsStruct}),
                 diff_ptr, repo, C_NULL, &gopts)
    free!(gopts.pathspec)
    return GitDiff(diff_ptr[1])
end

function diff_workdir(repo::GitRepo, left::String, opts::MaybeDict=nothing)
    return diff_workdir(repo, revparse(repo, left), opts)
end

function diff_workdir(repo::GitRepo, left::GitCommit, opts::MaybeDict=nothing)
    return diff_workdir(repo, GitTree(left), opts)
end

function diff_workdir(repo::GitRepo, left::GitTree, opts::MaybeDict=nothing)
    gopts = parse_git_diff_options(opts)
    diff_ptr = Ptr{Void}[0]
    @check ccall((:git_diff_tree_to_workdir, libgit2), Cint,
                 (Ptr{Ptr{Void}},
                  Ptr{Void},
                  Ptr{Void},
                  Ptr{DiffOptionsStruct}),
                 diff_ptr, repo, left, &gopts)
    free!(gopts.pathspec)
    return GitDiff(diff_ptr[1])
end

parse_git_diff_options(o::Nothing) = DiffOptionsStruct()

function parse_git_diff_options(opts::Dict)
    max_size = zero(Coff_t)
    if haskey(opts, :max_size)
        max_size = convert(Coff_t, opts[:max_size])
    end
    context_lines = uint16(3)
    if haskey(opts, :context_lines)
        @assert opts[:context_lines] <= typemax(Uint16)
        context_lines = convert(Uint16, opts[:context_lines])
    end
    interhunk_lines = zero(Uint16)
    if haskey(opts, :interhunk_lines)
        @assert opts[:interhunk_lines] <= typemax(Uint16)
        interhunk_lines = convert(Uint16, opts[:interhunk_lines])
    end
    flags = zero(Uint32)
    if get(opts, :reverse, false)
        flags |= GitConst.DIFF_REVERSE
    end
    if get(opts, :force_text, false)
        flags |= GitConst.DIFF_FORCE_TEXT
    end
    if get(opts, :ignore_whitespace, false)
        flags |= GitConst.DIFF_IGNORE_WHITESPACE
    end
    if get(opts, :ignore_whitespace_change, false)
        flags |= GitConst.DIFF_IGNORE_WHITESPACE_CHANGE
    end
    if get(opts, :ignore_whitespace_eol, false)
        flags |= GitConst.DIFF_IGNORE_WHITESPACE_EOL
    end
    if get(opts, :ignore_submodules, false)
        flags |= GitConst.DIFF_IGNORE_SUBMODULES
    end
    if get(opts, :patience, false)
        flags |= GitConst.DIFF_PATIENCE
    end
    if get(opts, :minimal, false)
        flags |= GitConst.DIFF_MINIMAL
    end
    if get(opts, :include_ignored, false)
        flags |= GitConst.DIFF_INCLUDE_IGNORED
    end
    if get(opts, :include_untracked, false)
        flags |= GitConst.DIFF_INCLUDE_UNTRACKED
    end
    if get(opts, :include_unmodified, false)
        flags |= GitConst.DIFF_INCLUDE_UNMODIFIED
    end
    if get(opts, :recurse_untracked_dirs, false)
       flags |= GitConst.DIFF_RECURSE_UNTRACKED_DIRS
    end
    if get(opts, :disable_pathspec_match, false)
       flags |= GitConst.DIFF_DISABLE_PATHSPEC_MATCH
    end
    if get(opts, :show_untracked_content, false)
       flags |= GitConst.DIFF_SHOW_UNTRACKED_CONTENT
    end
    if get(opts, :skip_binary_check, false)
       flags |= GitConst.DIFF_SKIP_BINARY_CHECK
    end
    if get(opts, :include_typechange, false)
       flags |= GitConst.DIFF_INCLUDE_TYPECHANGE
    end
    if get(opts, :include_typechange_trees, false)
       flags |= GitConst.DIFF_INCLUDE_TYPECHANGE_TREES
    end
    if get(opts, :ignore_filemode, false)
       flags |= GitConst.DIFF_IGNORE_FILEMODE
    end
    if get(opts, :recurse_ignored_dirs, false)
       flags |= GitConst.DIFF_RECURSE_IGNORED_DIRS
    end
    if haskey(opts, :paths)
        paths = opts[:paths]
        if !(isa(paths, Vector{ASCIIString}) || isa(paths, Vector{UTF8String}))
            throw(ArgumentError("opts[:paths] must be of type Vector{ByteString}"))
        end
        pathspec = StrArrayStruct(paths)
    else
        pathspec = StrArrayStruct()
    end
    #TODO: rest of the options
    return DiffOptionsStruct(GitConst.DIFF_OPTIONS_VERSION,
                             flags,
                             GitConst.SUBMODULE_IGNORE_DEFAULT,
                             pathspec,
                             C_NULL,
                             C_NULL,
                             context_lines,  # defaults to 3
                             interhunk_lines,
                             0,
                             max_size,
                             C_NULL,
                             C_NULL)
end


text:
# Autogenerated wrapper script for MibS_jll for aarch64-apple-darwin
export libmibs, mibs

using CoinUtils_jll
using Osi_jll
using Clp_jll
using Cgl_jll
using Cbc_jll
using SYMPHONY_jll
using ALPS_jll
using BiCePS_jll
using CHiPPS_BLIS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MibS")
JLLWrappers.@declare_library_product(libmibs, "@rpath/libMibs.0.dylib")
JLLWrappers.@declare_executable_product(mibs)
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, Osi_jll, Clp_jll, Cgl_jll, Cbc_jll, SYMPHONY_jll, ALPS_jll, BiCePS_jll, CHiPPS_BLIS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libmibs,
        "lib/libMibs.0.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        mibs,
        "bin/mibs",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# refactored to JuliaReach/CarlemanLinearization.jl
# include("error_bounds.jl")

# extends functionality in JuliaReach/CarlemanLinearization.jl to work with LazySets.jl types
include("kronecker.jl")


text:
using TensorBasic
using Base.Test

# write your own tests here
@eval x = rand(2,3,4);
@eval y = tensor_permute(x, [1,2,3], [2,1,3]);
@test y == permutedims(x, [2,1,3]);


text:
######### Jags batch program example  ###########

using StatsPlots, Jags, Statistics

cd(ProjDir) do

  dyes = "
  model {
     for (i in 1:BATCHES) {
        for (j in 1:SAMPLES) {
           y[i,j] ~ dnorm(mu[i], tau.within);
        }
        mu[i] ~ dnorm(theta, tau.between);
     }

     theta ~ dnorm(0.0, 1.0E-10);
     tau.within ~ dgamma(0.001, 0.001);
     s2.within <- 1/tau.within;
     tau.between ~ dgamma(0.001, 0.001);
     s2.between <- 1/tau.between;
     s2.total <- s2.within + s2.between;
     f.within <- s2.within/s2.total;
     f.between <- s2.between/s2.total;
  }
  "

  data = Dict{String, Any}()
  data["y"] = reshape([
    1545, 1540, 1595, 1445, 1595,
    1520, 1440, 1555, 1550, 1440,
    1630, 1455, 1440, 1490, 1605,
    1595, 1515, 1450, 1520, 1560,
    1510, 1465, 1635, 1480, 1580,
    1495, 1560, 1545, 1625, 1445
  ], 6, 5)
  data["BATCHES"] = 6
  data["SAMPLES"] = 5

  inits = [
    Dict{String, Any}(
      "theta" => 1500,
      "tau.within" => 1,
      "tau.between" => 1
  	)
  ]

  monitors = Dict{String, Any}(
    "theta" => true,
    "s2.within" => true,
    "s2.between" => true
  	)

  jagsmodel = Jagsmodel(name="dyes", model=dyes,
    #ncommands=3, nchains=3,
    monitor=monitors)

  println("\nJagsmodel that will be used:")
  jagsmodel |> display
  println("Input observed data dictionary:")
  data |> display
  println("\nInput initial values dictionary:")
  inits |> display
  println()

  @time sim = jags(jagsmodel, data, inits, ProjDir)
  println()

  ## Plotting
  p = plot(sim)

end #cd


text:
module stoosc

    using Reexport

    include("FindPeaks/FindPeaks.jl")
    include("MiscFunctions/MiscFunctions.jl")

    @reexport begin

        # Unregistered packages
        using .FindPeaks
        using .MiscFunctions

        # Registered packages
        using DataFrames
        using DifferentialEquations
        using PyPlot
        using Statistics
        using StatsBase

        # import CurveFit
        import CSV  # https://github.com/JuliaData/CSV.jl/issues/981
        import ProgressMeter

    end

    export estimate_entrainment_area, plot_arnold_tongue
    include("arnold_tongue.jl")

    export save_data, load_data, save_figure
    include("input_output.jl")

    export load_model
    include("load_model.jl")

    export kfr, create_events, plot_events, estimate_period
    include("miscellaneous.jl")

    export set_initial_conditions, set_timespan, set_output, create_model,
    set_solver, get_parameter_index, get_parameter_value, set_parameter_value,
    set_input, simulate_population, scan, estimate_binary_arnold_tongue
    include("Model.jl")

    export plot_solution, select_time
    include("PopulationSolution.jl")

end  # module

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
# Global functions and data - https://mccdaq.github.io/daqhats/c.html

# mcc daqhats wrapper for c function calls
# for effective development, do the following
#    using Pkg
#    Pkg.add("Revise")
#    using Revise
#    includet("~/daqhatsJulia/daqhats.jl")
#
# then modify daqhats.jl save and check in the Julia REPL

# generic format of ccall is:
# ccall((:cFuncName, "library"), funcReturnType, (tupleInputTypes), argumentsPassedIn)
# ccall automatically converts argumentPassedIn to tupleInputTypes

# path to so library

# lib = "/usr/local/lib/libdaqhats.so" moved to src

struct HatInfoTemp
    address::UInt8 			# The board address.
    id::UInt16 				# The product ID, one of [HatIDs](@ref HatIDs)
    version::UInt16 			# The hardware version
    product_name::NTuple{256,UInt8}	# The product name (c - initialized to 256 characters)
end

struct HatInfo
    address::UInt8 			# The board address.
    id::String 				# The product ID, one of [HatIDs](@ref HatIDs)
    version::UInt16 		# The hardware version
    product_name::String	# The product name
end

"""
	function hat_list_dict()

This is a private function to give the dictionary for card type and card
number association.  It could be public in the future.
"""
function hat_list_dict()
	idDict = Dict{Symbol, UInt16}(         # match HAT symbol to UInt16
	:ANY                => 0,       # Match any DAQ HAT ID in hatlist().
	:MCC_118            => 0x0142,  # MCC 118 ID.
	:MCC_118_BOOTLOADER => 0x8142,  # MCC 118 in firmware update mode ID.
	:MCC_128            => 0x0146,  # MCC 128 ID.
	:MCC_134            => 0x0143,  # MCC 134 ID.
	:MCC_152            => 0x0144,  # MCC 152 ID.
	:MCC_172            => 0x0145)  # MCC 172 ID.

	ridDict = Dict(values(idDict) .=> keys(idDict))  # for reverse lookup
	return idDict, ridDict
end

"""
	hat_list_count(filter_id::Symbol)

Count detected DAQ HAT boards.

filter_id types: ``:ANY, :MCC_118, :MCC_118_BOOTLOADER, :MCC_128, :MCC_134, :MCC_152, :MCC_172``
"""
function hat_list_count(filter_id::Symbol)
		idDict, ridDict = hat_list_dict()
		# number of HATS
		count = ccall((:hat_list, "libdaqhats"),
		Cint, (UInt16, Ptr{Cvoid}), idDict[filter_id], C_NULL)
		return count
end

"""
	hat_list([filter_id::Symbol, [number]]; count = false)

Return a list of detected DAQ HAT boards.

`hat_list()` : Synonym for hat_list(:ANY)

`hat_list(filter_id::Symbol)` : All detected DAQ HAT boards of specified type.

`hat_list(filter_id::Symbol, number::Integer)` : List of detected DAQ HAT boards up to the specified number.

`filter_id types` : ``:ANY, :MCC_118, :MCC_118_BOOTLOADER, :MCC_128, :MCC_134, :MCC_152, :MCC_172``
"""
function hat_list(; count::Bool = false)
	hat_list(:ANY, count=count)	 # Match any DAQ HAT ID in hatlist()
end

function hat_list(filter_id::Symbol; count::Bool = false)
	# associate filter_id with number
	# idDict, ridDict = hat_list_dict()
	# get number of HATS of specified type
	number = hat_list_count(filter_id)
	# get the structure of information
	if count
		return number
	end
	list = hat_list(filter_id, number, count=count)
	return list
end

function hat_list(filter_id::Symbol, number::Integer; count::Bool = false)

	if number < 1
		error("Number of HATS must be >= 1")
	end
	idDict, ridDict = hat_list_dict()  # dictionary of hat types

	# number of HATS installed
	numbermax = ccall((:hat_list, "libdaqhats"),
	Cint, (UInt16, Ptr{Cvoid}), idDict[filter_id], C_NULL)

	if count
		return numbermax
	end

	if number <= numbermax
		# generate variables for ccall below
		listTemp = Vector{HatInfoTemp}(undef,number)
		list     = Vector{HatInfo}(undef,number)
		for i = 1:number
			listTemp[i] = HatInfoTemp(0, 0, 0, map(UInt8, (string(repeat(" ",256))...,)))
		end

		success = ccall((:hat_list, "libdaqhats"),
		Cint, (UInt16, Ptr{HatInfoTemp}), idDict[filter_id], listTemp)
		# get only the required text from .product_name
		for i = 1:number
			a = String([listTemp[i].product_name...])
			f = findfirst(Char(0x0), a)
			p = SubString(a, 1:f-1)
			keyvalue = string(ridDict[listTemp[i].id])
			list[i] = HatInfo(listTemp[i].address, keyvalue, listTemp[i].version, p)
		end
		return list
	else
		error("Requesting $number '$filter_id' HATS $numbermax HAT(S) installed")
	end
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
using CleanCMB
using Printf, CSV, DataFrames
using Mmap
using Statistics
using Tables
dir = "../data"
# %% Simulation parameters
nrz = 1500
ℓmin, ℓmax = 30, 260 # ℓ range for fitting
Alens = 1
# %% Specification of the experiments
ν = [27, 39, 93, 145, 225, 280, 350, 410, 850] # in GHz
#νused = [true, true, true, true, true, true, true, true, true] # Which frequencies to use for fitting?
νused = [true, true, true, true, true, true, false, false, false]
nν = length(ν)
kν = findall(x -> x == true, νused)
# %% Read in binned theory power spectra
clt_th = CSV.read(joinpath(dir, "tensor_eebb_binned.csv"), DataFrame)
cls_th = CSV.read(joinpath(dir, "scalar_eebb_binned.csv"), DataFrame)
ell_eff = clt_th.leff
nbands = length(ell_eff)
rclass = 0.01
# %% Read in covariance matrices and perform ILC
io = open(joinpath(dir, "cov_bb_foreg.dat"), "r")
cov3 = Mmap.mmap(io, Array{Float64,3}, (nν, nν, nbands))
close(io)
cl1, cl2, cl3 = zeros(nbands, nrz), zeros(nbands, nrz), zeros(nbands, nrz)
for irz = 1:nrz
    ioin = open(joinpath(dir, @sprintf("cov_bb_total_irz%04d.dat", irz)), "r")
    cov1 = Mmap.mmap(ioin, Array{Float64,3}, (nν, nν, nbands))
    close(ioin)
    ioin = open(joinpath(dir, @sprintf("cov_bb_noise_irz%04d.dat", irz)), "r")
    cov2 = Mmap.mmap(ioin, Array{Float64,3}, (nν, nν, nbands))
    close(ioin)
    w = ilc_weights(cov1[kν, kν, :])
    cl1[:, irz] = ilc_clean_cij(cov1[kν, kν, :], w)
    cl2[:, irz] = ilc_clean_cij(cov2[kν, kν, :], w)
    cl3[:, irz] = ilc_clean_cij(cov3[kν, kν, :], w)
end
m1, m2, m3 = mean(cl1, dims = 2), mean(cl2, dims = 2), mean(cl3, dims = 2)
v1 = var(cl1, dims = 2)

# %% Joint fit for the tensor-to-scalar ratio and the foreground amplitude
ii = findall(x -> x >= ℓmin && x <= ℓmax, ell_eff)
r, rw = zeros(nrz), zeros(nrz, 2)
y1 = clt_th.bb[ii] / rclass
y2 = m3[ii]
v = v1[ii]
Fij = [
    sum(y1 .* y1 ./ v) sum(y1 .* y2 ./ v)
    sum(y2 .* y1 ./ v) sum(y2 .* y2 ./ v)
] # 2x2 Fisher matrix for the tensor-to-scalar ratio and the foreground amplitude
Cij = inv(Fij) # Covariance matrix
for irz = 1:nrz
    x = cl1[ii, irz] .- m2[ii] .- Alens * cls_th.bb[ii]
    r[irz] = sum(x .* y1 ./ v) / sum(y1 .^ 2 ./ v)
    z = [sum(x .* y1 ./ v), sum(x .* y2 ./ v)]
    rw[irz, 1:2] = Fij \ z
end
println("Fitted ℓs: ", ell_eff[ii])
println("Used νs: ", ν[kν], " GHz")
println("Without foreground marginalisation:")
println("- r = ", mean(r), " ± ", std(r))
println("- Fisher error = ", 1 / √sum(y1 .^ 2 ./ v))
println("With foreground marginalisation:")
println("- r = ", mean(rw[:, 1]), " ± ", std(rw[:, 1]))
println("- Fisher error = ", √Cij[1, 1])
# println("- FG = ", mean(rw[:, 2]), " ± ", std(rw[:, 2]))
# println("- Fisher error = ", √Cij[2, 2])
# %% Save to the file
t = Tables.table([1:nrz r rw[:, 1]])
CSV.write("ilc_results.csv", t, header = ["irz", "r_wo_FGmarg", "r_w_FGmarg"])


text:
using IntervalArithmetic
using Test

let b

@testset "DecoratedInterval tests" begin
    a = DecoratedInterval(@interval(1, 2), com)
    @test decoration(a) == com

    b = sqrt(a)
    @test interval(b) == sqrt(interval(a))
    @test decoration(b) == com

    a = DecoratedInterval(@interval(-1, 1), com)
    b = sqrt(a)
    @test interval(b) == sqrt(Interval(0, 1))
    @test decoration(b) == trv

    d = DecoratedInterval(a, dac)
    @test decoration(d) == dac

    @test decoration(DecoratedInterval(1.1)) == com
    @test decoration(DecoratedInterval(1.1, dac)) == dac

    @test decoration(DecoratedInterval(2, 0.1, com)) == ill
    @test decoration(DecoratedInterval(2, 0.1)) == ill
    @test isnai(DecoratedInterval(2, 0.1))
    @test decoration(@decorated(2, 0.1)) == ill
    @test decoration(DecoratedInterval(big(2), big(1))) == ill
    @test isnai((DecoratedInterval(big(2), big(1))))
    @test isnai(@decorated(big(2), big(1)))

    # Disabling the following tests, because Julia 0.5 has some strange behaviour here
    # @test_throws ArgumentError DecoratedInterval(BigInt(1), 1//10)
    # @test_throws ArgumentError @decorated(BigInt(1), 1//10)

    # Tests related to powers of decorated Intervals
    @test @decorated(2,3) ^ 2 == DecoratedInterval(4, 9)
    @test @decorated(2,3) ^ -2 == DecoratedInterval(1/9,1/4)
    @test @decorated(-3,2) ^ 3 == DecoratedInterval(-27., 8.)
    @test @decorated(-3,-2) ^ -3 == DecoratedInterval(-1/8.,-1/27)
    @test @decorated(0,3) ^ 2 == DecoratedInterval(0, 9)
    @test @decorated(0,3) ^ -2 == DecoratedInterval(1/9, Inf, trv)
    @test @decorated(2,3)^Interval(0.0, 1.0) == DecoratedInterval(1.0,3.0)
    @test @decorated(2,3)^@decorated(0.0, 1.0) == DecoratedInterval(1.0,3.0)
    @test @decorated(0, 2)^Interval(0.0, 1.0) == DecoratedInterval(0.0,2.0, trv)
    @test @decorated(0, 2)^@decorated(0.0, 1.0) == DecoratedInterval(0.0,2.0, trv)
    @test @decorated(-3, 2)^Interval(0.0, 1.0) == DecoratedInterval(0.0,2.0, trv)
    @test @decorated(-3, 2)^@decorated(0.0, 1.0) == DecoratedInterval(0.0,2.0, trv)
    @test @decorated(-3, 2)^Interval(-1.0, 1.0) == DecoratedInterval(0.0,Inf, trv)
    @test @decorated(-3, 2)^@decorated(-1.0, 1.0) == DecoratedInterval(0.0, Inf, trv)
end

@testset "Convert string to DecoratedInterval" begin
    @test convert(DecoratedInterval{Float64}, "[1,2]") ==
                    DecoratedInterval(Interval(1, 2), com)

    @test convert(DecoratedInterval{Float64}, "[1,2]_dac") ==
                        DecoratedInterval(Interval(1, 2), dac)

end

end


text:
# RTM of the BP synthetic 2007 data set
# Author: Philipp Witte
# Date: March 2018
#

# TO DO
# Set up path where data will be saved
data_path = "/path/to/data/"

using JUDI, SegyIO, JLD, PyPlot, LinearAlgebra

# Load velocity model
if !isfile("bp_synthetic_2004_migration_velocity.jld")
    run(`wget ftp://slim.gatech.edu/data/SoftwareRelease/Imaging.jl/CompressiveLSRTM/bp_synthetic_2004_migration_velocity.jld`)
end
vp = load(join([pwd(), "/bp_synthetic_2004_migration_velocity.jld"]))["vp"] / 1f3

# Set up model structure
d = (6.25, 6.25)
o = (0., 0.)
m0 = (1f0 ./ vp).^2
n = size(m0)
model0 = Model(n, d, o, m0)

# Scan directory for segy files and create out-of-core data container
container = segy_scan(data_path, "bp_observed_data", ["GroupX", "GroupY", "RecGroupElevation", "SourceSurfaceElevation", "dt"])
d_obs = judiVector(container; segy_depth_key = "SourceDepth")

# Set up source
src_geometry = Geometry(container; key = "source")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.020)  # 27 Hz peak frequency
q = judiVector(src_geometry, wavelet)

# Info structure for linear operators
ntComp = get_computational_nt(src_geometry, d_obs.geometry, model0)    # no. of computational time steps
info = Info(prod(model0.n), d_obs.nsrc, ntComp)

###################################################################################################

# Enable optimal checkpointing
opt = Options(limit_m = true,
              buffer_size = 3000f0,
              isic = true,
			  dft_subsampling_factor=8)

# Setup operators
q_dist = generate_distribution(q)
F = judiModeling(info, model0, q.geometry, d_obs.geometry; options=opt)
J = judiJacobian(F, q)

# Set up random frequencies
nfreq = 10
J.options.frequencies = Array{Any}(undef, d_obs.nsrc)
J.options.dft_subsampling_factor = 8
for k=1:d_obs.nsrc
    J.options.frequencies[k] = select_frequencies(q_dist; fmin=0.003, fmax=0.04, nf=nfreq)
end

# Add random noise
d_lin = get_data(d_obs)
for j=1:d_lin.nsrc
    noise = randn(Float32, size(d_lin.data[j]))
    d_lin.data[j] += (noise/norm(vec(noise), Inf)*0.02f0)
end

# Topmute
Ml = judiMarineTopmute2D(35, d_lin.geometry)    # data topmute
d_lin = Ml*d_lin

# RTM
rtm = adjoint(J)*d_lin
save("bp_synethic_2004_rtm_frequency.jld", "rtm", rtm)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI