text:
"""Functions for computing CIB models."""

abstract type AbstractCIBModel{T<:Real} <: AbstractForegroundModel end

"""
    CIB_Planck2013{T}(; kwargs...)

Define CIB model parameters. Defaults are from Viero et al. 2013. All numbers
not typed are converted to type T. This model has the following parameters and default values:

* `nside::Int64 = 4096`
* `hod::String = "shang"`
* `Inu_norm = 0.3180384`
* `min_redshift = 0.0`
* `max_redshift = 5.0`
* `min_mass = 1e12`
* `box_size = 40000`
* `shang_zplat = 2.0`
* `shang_Td = 20.7`
* `shang_betan = 1.6`
* `shang_eta = 2.4`
* `shang_alpha = 0.2`
* `shang_Mpeak = 10^12.3`
* `shang_sigmaM = 0.3`
* `shang_Msmin = 1e11`
* `shang_Mmin = 1e10`
* `shang_I0 = 46`
* `jiang_gamma_1 = 0.13`
* `jiang_alpha_1 = -0.83`
* `jiang_gamma_2 = 1.33`
* `jiang_alpha_2 = -0.02`
* `jiang_beta_2 = 5.67`
* `jiang_zeta = 1.19`

"""
@with_kw struct CIB_Planck2013{T<:Real} <: AbstractCIBModel{T}
    @deftype T
    nside::Int64 = 4096
    hod::String = "shang"
    Inu_norm = 0.3180384
    min_redshift = 0.0
    max_redshift = 6.0
    min_mass = 2e10
    box_size = 40000

    # shang HOD
    shang_zplat = 2.0
    shang_Td = 20.7
    shang_beta = 1.6
    shang_eta = 2.4
    shang_alpha = 0.2
    shang_Mpeak = 10^12.3
    shang_sigmaM = 0.3
    shang_Msmin = 1e11
    shang_Mmin = 1e10
    shang_I0 = 46

    # jiang
    jiang_gamma_1 = 0.13
    jiang_alpha_1 = -0.83
    jiang_gamma_2 = 1.33
    jiang_alpha_2 = -0.02
    jiang_beta_2 = 5.67
    jiang_zeta = 1.19
end

function jiang_shmf(m, M_halo, model)
    dndm = (
        (
            (model.jiang_gamma_1 * ((m / M_halo)^model.jiang_alpha_1)) +
            (model.jiang_gamma_2 * ((m / M_halo)^model.jiang_alpha_2))
        ) * (exp(-(model.jiang_beta_2) * ((m / M_halo)^model.jiang_zeta)))
    )
    return dndm
end

"""
Build a linear interpolation function which maps log(M_h) to N_sat.
"""
function build_shang_interpolator(
    min_log_M::T,
    max_log_M::T,
    model::AbstractCIBModel;
    n_bin = 1000,
) where {T}

    x_m = LinRange(min_log_M, max_log_M, 1000)
    N_sat_i = zero(x_m)

    function integrand_m(lm, lM_halo)
        # lm is natural log of m
        m = exp(lm)
        M_halo = exp(lM_halo)
        dns_dm = jiang_shmf(m, M_halo, model)
        return dns_dm
    end

    for i = 1:size(x_m, 1)
        N_sat_i[i], err = quadgk(
            t -> integrand_m(t, x_m[i]),
            log(model.shang_Msmin),
            x_m[i],
            rtol = 1e-8,
        )
        N_sat_i[i] = convert(T, max(0.0, N_sat_i[i]))
    end

    return LinearInterpolation(x_m, N_sat_i)
end

function sigma_cen(m::T, model::AbstractCIBModel) where {T}
    return (
        exp(
            -(log10(m) - log10(model.shang_Mpeak))^2 /
            (T(2) * model.shang_sigmaM),
        ) * m
    ) / sqrt(T(2π) * model.shang_sigmaM)
end

function nu2theta(nu::T, z::T, model::AbstractCIBModel) where {T}
    phys_h = T(6.62606957e-27)   # erg.s
    phys_k = T(1.3806488e-16)    # erg/K
    Td = model.shang_Td * (one(T) + z)^model.shang_alpha
    xnu = phys_h * nu / phys_k / Td
    return xnu^(T(4) + model.shang_beta) / expm1(xnu) / nu / model.shang_I0
end

"""
<L_sat> interpolation values
"""
function integrand_L(lm, lM_halo, model::AbstractCIBModel)
    m = exp(lm)
    return sigma_cen(m, model) * jiang_shmf(m, exp(lM_halo), model)
end

"""
Build a linear interpolator that takes in ln(M_halo) and returns sigma.
"""
function build_sigma_sat_ln_interpolator(
    max_ln_M::T,
    model::AbstractCIBModel;
    n_bin = 1000,
) where {T}
    x = LinRange(log(model.shang_Mmin), max_ln_M, n_bin)
    L_mean = zero(x)

    for i = 1:n_bin
        L_mean[i], err = quadgk(
            t -> integrand_L(t, x[i], model),
            log(model.shang_Mmin),
            x[i],
            rtol = 1e-6,
        )
    end
    return LinearInterpolation(T.(x), T.(L_mean), extrapolation_bc = zero(T))
end

function build_muofn_interpolator(
    model;
    min_mu::T = -6.0f0,
    max_mu::T = 0.0f0,
    nbin = 1000,
) where {T}
    mu = T(10) .^ LinRange(min_mu, max_mu, nbin)
    n = zero(mu)
    for i = 1:nbin
        n[i], err = quadgk(
            lmu -> jiang_shmf(exp(lmu), one(T), model),
            log(mu[i]),
            0.0f0,
            rtol = 1.0f-6,
        )
    end
    return LinearInterpolation(T.(reverse(n)), T.(reverse(mu)))
end

"""
Compute redshift evolution factor for LF.
"""
function shang_z_evo(z::T, model::AbstractCIBModel) where {T}
    return (one(T) + min(z, model.shang_zplat))^model.shang_eta
end

"""
Construct the necessary interpolator set.
"""
function get_interpolators(
    model::AbstractCIBModel,
    cosmo::Cosmology.FlatLCDM{T},
    min_halo_mass::T,
    max_halo_mass::T,
) where {T}
    return (
        r2z = build_r2z_interpolator(
            model.min_redshift,
            model.max_redshift,
            cosmo,
        ),
        hod_shang = build_shang_interpolator(
            log(min_halo_mass),
            log(max_halo_mass),
            model,
        ),
        c_lnm2r = build_c_lnm2r_interpolator(),
        sigma_sat = build_sigma_sat_ln_interpolator(log(max_halo_mass), model),
        muofn = build_muofn_interpolator(model),
    )
end

# Fill up arrays with information related to CIB central sources.
function process_centrals!(
    model::AbstractCIBModel{T},
    cosmo::Cosmology.FlatLCDM{T},
    Healpix_res::Resolution;
    interp,
    hp_ind_cen,
    dist_cen,
    redshift_cen,
    lum_cen,
    n_sat_bar,
    n_sat_bar_result,
    halo_pos,
    halo_mass,
) where {T}

    N_halos = size(halo_mass, 1)

    Threads.@threads for i = 1:N_halos
        # location information for centrals
        hp_ind_cen[i] = Healpix.vec2pixRing(
            Healpix_res,
            halo_pos[1, i],
            halo_pos[2, i],
            halo_pos[3, i],
        )
        dist_cen[i] =
            sqrt(halo_pos[1, i]^2 + halo_pos[2, i]^2 + halo_pos[3, i]^2)
        redshift_cen[i] = interp.r2z(dist_cen[i])

        # compute HOD
        n_sat_bar[i] = interp.hod_shang(log(halo_mass[i]))
        n_sat_bar_result[i] =
            rand(Distributions.Poisson(Float64.(n_sat_bar[i])))

        # get central luminosity
        lum_cen[i] =
            sigma_cen(halo_mass[i], model) * shang_z_evo(redshift_cen[i], model)
    end
end

# Fill up arrays with information related to CIB satellites.
function process_sats!(
    model::AbstractCIBModel{T},
    cosmo::Cosmology.FlatLCDM{T},
    Healpix_res::Resolution;
    interp,
    hp_ind_sat,
    dist_sat,
    redshift_sat,
    lum_sat,
    cumsat,
    halo_mass,
    halo_pos,
    redshift_cen,
    n_sat_bar,
    n_sat_bar_result,
) where {T}

    N_halos = size(halo_mass, 1)
    Threads.@threads for i_halo = 1:N_halos
        r_cen = m2r(halo_mass[i_halo], cosmo)
        c_cen = mz2c(halo_mass[i_halo], redshift_cen[i_halo], cosmo)
        for j = 1:n_sat_bar_result[i_halo]
            # i is central halo index, j is index of satellite within each halo
            i_sat = cumsat[i_halo] + j # index of satellite in satellite arrays

            log_msat_inner = max(log(rand(T)), T(-7.0))
            r_sat =
                r_cen *
                interp.c_lnm2r(c_cen, log_msat_inner) *
                T(200.0^(-1.0 / 3.0))
            m_sat =
                (interp.muofn(rand(T) * n_sat_bar[i_halo]) * halo_mass[i_halo])

            phi = random_phi(T)
            theta = random_theta(T)
            x_sat = halo_pos[1, i_halo] + r_sat * sin(theta) * cos(phi)
            y_sat = halo_pos[2, i_halo] + r_sat * sin(theta) * sin(phi)
            z_sat = halo_pos[3, i_halo] + r_sat * cos(theta)
            dist_sat[i_sat] = sqrt(x_sat^2 + y_sat^2 + z_sat^2)
            redshift_sat[i_sat] = interp.r2z(dist_sat[i_sat])

            # lum_sat[i_sat] = interp.sigma_sat(log(m_sat)) * shang_z_evo(
            #     redshift_sat[i], model)
            lum_sat[i_sat] =
                sigma_cen(m_sat, model) *
                shang_z_evo(redshift_sat[i_sat], model)
            hp_ind_sat[i_sat] =
                Healpix.vec2pixRing(Healpix_res, x_sat, y_sat, z_sat)
        end
    end
end

"""
    generate_sources(model, cosmo, halo_pos_inp, halo_mass_inp; verbose=true)

Produce a source catalog from a model and halo catalog. This converts the
halo arrays into the type specified by `model`.

# Arguments:
- `model::AbstractCIBModel{T}`: source model parameters
- `cosmo::Cosmology.FlatLCDM{T}`: background cosmology
- `Healpix_res::Resolution`: Healpix map resolution
- `halo_pos_inp::AbstractArray{TH,2}`: halo positions with dims (3, nhalos)
- `halo_mass_inp::AbstractArray{TH,1}`: halo masses

# Keywords
- `verbose::Bool=true`: print out progress details
"""
function generate_sources(
    model::AbstractCIBModel{T},
    cosmo::Cosmology.FlatLCDM{T},
    halo_pos_inp::AbstractArray{TH,2},
    halo_mass_inp::AbstractArray{TH,1};
    verbose = true,
) where {T,TH}

    # make sure halo inputs are the CIB type
    halo_pos = convert(Array{T,2}, halo_pos_inp)
    halo_mass = convert(Array{T,1}, halo_mass_inp)

    # set up basics
    N_halos = size(halo_mass, 1)
    interp =
        get_interpolators(model, cosmo, minimum(halo_mass), maximum(halo_mass))
    res = Resolution(model.nside)

    verbose && println("Allocating for $(N_halos) centrals.")
    hp_ind_cen = Array{Int64}(undef, N_halos)  # healpix index of halo
    lum_cen = Array{T}(undef, N_halos)  # Lum of central w/o ν-dependence
    redshift_cen = Array{T}(undef, N_halos)
    dist_cen = Array{T}(undef, N_halos)
    n_sat_bar = Array{T}(undef, N_halos)
    n_sat_bar_result = Array{Int32}(undef, N_halos)

    # STEP 1: compute central properties -----------------------------------
    verbose && println("Processing centrals on $(Threads.nthreads()) threads.")
    process_centrals!(
        model,
        cosmo,
        res,
        interp = interp,
        hp_ind_cen = hp_ind_cen,
        dist_cen = dist_cen,
        redshift_cen = redshift_cen,
        lum_cen = lum_cen,
        n_sat_bar = n_sat_bar,
        n_sat_bar_result = n_sat_bar_result,
        halo_pos = halo_pos,
        halo_mass = halo_mass,
    )

    # STEP 2: Generate satellite arrays -----------------------------
    cumsat = generate_subhalo_offsets(n_sat_bar_result)
    total_n_sat = cumsat[end]
    hp_ind_sat = Array{Int64}(undef, total_n_sat)  # healpix index of halo
    lum_sat = Array{T}(undef, total_n_sat)  # Lum of central w/o ν-dependence
    redshift_sat = Array{T}(undef, total_n_sat)
    dist_sat = Array{T}(undef, total_n_sat)

    # STEP 3: compute satellite properties -----------------------------------
    verbose && println("Processing $(total_n_sat) satellites.")
    process_sats!(
        model,
        cosmo,
        res,
        interp = interp,
        hp_ind_sat = hp_ind_sat,
        dist_sat = dist_sat,
        redshift_sat = redshift_sat,
        lum_sat = lum_sat,
        cumsat = cumsat,
        halo_mass = halo_mass,
        halo_pos = halo_pos,
        redshift_cen = redshift_cen,
        n_sat_bar = n_sat_bar,
        n_sat_bar_result = n_sat_bar_result,
    )

    return (
        hp_ind_cen = hp_ind_cen,
        lum_cen = lum_cen,
        redshift_cen = redshift_cen,
        dist_cen = dist_cen,
        hp_ind_sat = hp_ind_sat,
        lum_sat = lum_sat,
        redshift_sat = redshift_sat,
        dist_sat = dist_sat,
        N_cen = N_halos,
        N_sat = total_n_sat,
    )
end

"""
    paint!(result_map, nu_obs, model, sources, fluxes_cen, fluxes_sat)

Paint a source catalog onto a map, recording the fluxes in
`fluxes_cen` and `fluxes_sat`.

# Arguments:
- `result_map::Map{T_map, RingOrder}`: Healpix map to paint
- `nu_obs`: frequency in Hz
- `model::AbstractCIBModel{T}`: source model parameters
- `sources`: NamedTuple containing source information from generate_sources
- `fluxes_cen::AbstractArray`: buffer for writing fluxes of centrals
- `fluxes_sat::AbstractArray`: buffer for writing fluxes of satellites
"""
function paint!(
    result_map::Map{T_map,RingOrder},
    nu_obs,
    model::AbstractCIBModel{T},
    sources,
    fluxes_cen::AbstractArray,
    fluxes_sat::AbstractArray,
) where {T_map,T}

    pixel_array = result_map.pixels
    fill!(pixel_array, zero(T))  # prepare the frequency map

    # process centrals for this frequency
    Threads.@threads for i = 1:sources.N_cen
        nu = (one(T) + sources.redshift_cen[i]) * nu_obs
        fluxes_cen[i] = l2f(
            sources.lum_cen[i] * nu2theta(nu, sources.redshift_cen[i], model),
            sources.dist_cen[i],
            sources.redshift_cen[i],
        )
        pixel_array[sources.hp_ind_cen[i]] += fluxes_cen[i]
    end

    # process satellites for this frequency
    Threads.@threads for i = 1:sources.N_sat
        nu = (one(T) + sources.redshift_sat[i]) * nu_obs
        fluxes_sat[i] = l2f(
            sources.lum_sat[i] * nu2theta(nu, sources.redshift_sat[i], model),
            sources.dist_sat[i],
            sources.redshift_sat[i],
        )
        pixel_array[sources.hp_ind_sat[i]] += fluxes_sat[i]
    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(result_map.resolution.nside)
    pixel_array .*= per_pixel_steradian
end

"""
Paint a source catalog onto a map.

This function creates the arrays for you.
"""
function paint!(
    result_map::Map{T,RingOrder},
    nu_obs::T,
    model::AbstractCIBModel,
    sources,
) where {T}

    fluxes_cen = Array{T,1}(undef, sources.N_cen)
    fluxes_sat = Array{T,1}(undef, sources.N_sat)

    paint!(result_map, nu_obs, model, sources, fluxes_cen, fluxes_sat)

    return fluxes_cen, fluxes_sat
end

"""
    estimate_fluxes!(  nu_obs, model, sources, fluxes_cen, fluxes_sat)

Paint a source catalog onto a map, recording the fluxes in
`fluxes_cen` and `fluxes_sat`.

# Arguments:
- `nu_obs`: frequency in Hz
- `model::AbstractCIBModel{T}`: source model parameters
- `sources`: NamedTuple containing source information from generate_sources
- `fluxes_cen::AbstractArray`: buffer for writing fluxes of centrals
- `fluxes_sat::AbstractArray`: buffer for writing fluxes of satellites
"""
function estimate_fluxes!(
    nu_obs,
    model::AbstractCIBModel{T},
    sources,
    fluxes_cen::AbstractArray,
    fluxes_sat::AbstractArray,
) where {T}

    # process centrals for this frequency
    Threads.@threads for i = 1:sources.N_cen
        nu = (one(T) + sources.redshift_cen[i]) * nu_obs
        fluxes_cen[i] = l2f(
            sources.lum_cen[i] * nu2theta(nu, sources.redshift_cen[i], model),
            sources.dist_cen[i],
            sources.redshift_cen[i],
        )
    end

    # process satellites for this frequency
    Threads.@threads for i = 1:sources.N_sat
        nu = (one(T) + sources.redshift_sat[i]) * nu_obs
        fluxes_sat[i] = l2f(
            sources.lum_sat[i] * nu2theta(nu, sources.redshift_sat[i], model),
            sources.dist_sat[i],
            sources.redshift_sat[i],
        )
    end
end

"""
    paint_with_mask!(result_map, nu_obs, model, sources, fluxes_cen, fluxes_sat,mask_cen, mask_sat)

Paint a source catalog onto a map, recording the fluxes in
`fluxes_cen` and `fluxes_sat`.

# Arguments:
- `result_map::Map{T_map, RingOrder}`: Healpix map to paint
- `nu_obs`: frequency in Hz
- `model::AbstractCIBModel{T}`: source model parameters
- `sources`: NamedTuple containing source information from generate_sources
- `fluxes_cen::AbstractArray`: buffer for writing fluxes of centrals
- `fluxes_sat::AbstractArray`: buffer for writing fluxes of satellites
- `mask_cen::AbstractArray`: buffer encoding  the index of the mask which satisfies the condition for  centrals
- `mask_sat::AbstractArray`: buffer encoding  the index of the mask which satisfies the condition for satellites
"""

function paint_with_mask!(
    result_map::Map{T_map,RingOrder},
    nu_obs,
    model::AbstractCIBModel{T},
    sources,
    fluxes_cen::AbstractArray,
    fluxes_sat::AbstractArray,
    mask_cen::AbstractArray,
    mask_sat::AbstractArray,
) where {T_map,T}

    pixel_array = result_map.pixels
    fill!(pixel_array, zero(T))  # prepare the frequency map

    # process centrals for this frequency
    Threads.@threads for i in mask_cen
        nu = (one(T) + sources.redshift_cen[i]) * nu_obs
        fluxes_cen[i] = l2f(
            sources.lum_cen[i] * nu2theta(nu, sources.redshift_cen[i], model),
            sources.dist_cen[i],
            sources.redshift_cen[i],
        )
        pixel_array[sources.hp_ind_cen[i]] += fluxes_cen[i]

    end

    # process satellites for this frequency
    Threads.@threads for i in mask_sat
        nu = (one(T) + sources.redshift_sat[i]) * nu_obs
        fluxes_sat[i] = l2f(
            sources.lum_sat[i] * nu2theta(nu, sources.redshift_sat[i], model),
            sources.dist_sat[i],
            sources.redshift_sat[i],
        )
        pixel_array[sources.hp_ind_sat[i]] += fluxes_sat[i]

    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(result_map.resolution.nside)
    pixel_array .*= per_pixel_steradian
end

"""
    estimate_normalization(
    result_map,
    ν_ref,
    unit_out)

estimate  normalization factor   so that both the fluxes and the map have consistent monopole term with
Fixen et al. 1998 model.

# Arguments:
- `result_map::Map{T_map, RingOrder}`: Healpix CIB map
- `ν_ref`: frequency in GHz (Default 857 GHz )
- `unit_out::Unitful.FreeUnits`: set the output units of the Fixsen monopole

"""

function estimate_normalization(
    result_map::Map{T,RingOrder},
    ν_ref::Unitful.Frequency,
    unit_out::Unitful.FreeUnits,
) where {T}

    fixsen_monopole = fixsen_model(ν_ref; unit_out = unit_out)
    maskpatch = findall(x -> (x != 0), result_map)

    fsky = size(maskpatch)[1] / size(result_map)[1]
    Smean = mean(result_map) / fsky * unit_out
    L0 = fixsen_monopole / Smean
    return L0
end

"""
    make_multifrequency_dataset(
    ν_obs,
    model,
    sources;
    ν_ref= 857u"GHz",
    output_μm= false,
    normalize = false,
    outputpath  = "./cib_data.hdf5",
    unit_out  = u"MJy/sr",
)

This function involves many subroutines aimed at generating at multiple frequency
a dataset of frequency maps at `nu_obs` frequencies, given the CIB model and the
luminosities painted from a halo catalog. As the painting procedure is _normalization free_,
the user can set  the fluxes and the maps  to match the Fixsen monopole at a reference frequency.

At a reference freq. `ν_ref`, we scale   both maps and fluxes generated with XGPaint so that :

`<I(ν_ref)> =  F( ν_ref, T)`

with `F( ν_ref, T) ` being the model of Fixsen. We then apply the normalization factor

`L0 = F( ν_ref, T) /<I(ν_ref)>`

to the other frequencies.

# Arguments:

- `ν_obs::AbstractArray`: array of frequencies to be simulated (in units of GHz)
- `ν_ref::Unitful.Frequency`: reference frequency in GHz (Default 857 GHz )
- `model::AbstractCIBModel{T}`: source model parameters
- `sources::NamedTuple`: containing source information from generate_sources
- `output_μm::Bool`: save the output catalog in units of μm (default outputs in GHz)
- `normalize::Bool`: apply the normalization if `true` (default `false`)
- `outputpath::String`: generic string encoding the path to the folder to save the
    data. It   looks for the string sequence `data` and replaces it with the frequency
    (or wavelength ) of the specific simulated dataset.
- `unit_out::Unitful.FreeUnits`: set the output units of the map and of the fluxes
    (default MJy/sr)
"""

function make_multifrequency_dataset(
    ν_obs::AbstractArray,
    model,
    sources;
    ν_ref::Unitful.Frequency = 857u"GHz",
    output_μm::Bool = false,
    normalize::Bool = false,
    outputpath::String = "./cib_data.hdf5",
    unit_out::Unitful.FreeUnits = u"MJy/sr",
)

    m = Map{Float64,RingOrder}(model.nside)  # create a Healpix map
    fluxes_cen = Array{Float32,1}(undef, sources.N_cen)
    fluxes_sat = Array{Float32,1}(undef, sources.N_sat)
    if normalize
        @time paint!(
            m,
            convert(Float32, ustrip(uconvert(u"Hz", ν_ref))),
            model,
            sources,
            fluxes_cen,
            fluxes_sat,
        )
        @time L0 = estimate_normalization(m, ν_ref, unit_out)
        print("estimated the normalization constant :$(ustrip(L0)).  \n")
    end

    c = uconvert(u"μm* GHz", float(Unitful.c0))
    lambdas = c ./ ν_obs

    for i in eachindex(ν_obs)
        ν = ustrip(uconvert(u"Hz", ν_obs[i]))

        if output_μm == true
            λ = round(typeof(1u"μm"), lambdas[i])
            fname = replace(outputpath, "data" => "$(ustrip(λ))$(unit(λ))")
            print("Saving map at $(λ)  in $fname \n")
        else
            fname = replace(
                outputpath,
                "data" => "$(ustrip(ν_obs[i]))$(unit(ν_obs[i]))",
            )
            print("Saving map at $(ν_obs[i])  in $fname \n")

        end
        @time paint!(
            m,
            convert(Float32, ν),
            model,
            sources,
            fluxes_cen,
            fluxes_sat,
        )
        if normalize
            fluxes_cen *= L0
            fluxes_sat *= L0
            m.pixels *= L0
        end
        write_to_hdf5(fname, fluxes_cen, fluxes_sat, sources, m.pixels)

    end

end


text:
using Adapt
using Test

# custom array type

struct CustomArray{T,N} <: AbstractArray{T,N}
    arr::Array
end

CustomArray(x::Array{T,N}) where {T,N} = CustomArray{T,N}(x)
Adapt.adapt_storage(::Type{<:CustomArray}, xs::Array) = CustomArray(xs)

Base.size(x::CustomArray, y...) = size(x.arr, y...)
Base.getindex(x::CustomArray, y...) = getindex(x.arr, y...)
Base.count(x::CustomArray) = count(x.arr)

const mat = CustomArray{Float64,2}(rand(2,2))
const vec = CustomArray{Float64,1}(rand(2))

const mat_bools = CustomArray{Bool,2}(rand(Bool,2,2))

macro test_adapt(to, src, dst, typ=nothing)
    quote
        @test adapt($to, $src) == $dst
        @test typeof(adapt($to, $src)) == typeof($dst)
        if $typ !== nothing
            @test typeof($dst) <: $typ
        end
    end
end

AnyCustomArray{T,N} = Union{CustomArray,WrappedArray{T,N,CustomArray,CustomArray{T,N}}}

# basic adaption
@test_adapt CustomArray mat.arr mat

# idempotency
@test_adapt CustomArray mat mat

# custom wrapper
struct Wrapper{T}
    arr::T
end
Wrapper(x::T) where T = Wrapper{T}(x)
Adapt.adapt_structure(to, xs::Wrapper) = Wrapper(adapt(to, xs.arr))
@test_adapt CustomArray Wrapper(mat.arr) Wrapper(mat)

## base wrappers

@test @inferred(adapt(nothing, NamedTuple())) == NamedTuple()
@test_adapt CustomArray (mat.arr,) (mat,)
@test @allocated(adapt(nothing, ())) == 0
@test @allocated(adapt(nothing, (1,))) == 0
@test @allocated(adapt(nothing, (1,2,3,4,5,6,7,8,9,10))) == 0

@test_adapt CustomArray (a=mat.arr,) (a=mat,)

@test_adapt CustomArray view(mat.arr,:,:) view(mat,:,:) AnyCustomArray
const inds = CustomArray{Int,1}([1,2])
@test_adapt CustomArray view(mat.arr,inds.arr,:) view(mat,inds,:) AnyCustomArray

# NOTE: manual creation of PermutedDimsArray because permutedims collects
@test_adapt CustomArray PermutedDimsArray(mat.arr,(2,1)) PermutedDimsArray(mat,(2,1)) AnyCustomArray

# NOTE: manual creation of ReshapedArray because Base.Array has an optimized `reshape`
@test_adapt CustomArray Base.ReshapedArray(mat.arr,(2,2),()) reshape(mat,(2,2)) AnyCustomArray

@test_adapt CustomArray Base.LogicalIndex(mat_bools.arr) Base.LogicalIndex(mat_bools) AnyCustomArray

@test_adapt CustomArray reinterpret(Int64,mat.arr) reinterpret(Int64,mat) AnyCustomArray

## doubly-wrapped

@test_adapt CustomArray reinterpret(Int64,view(mat.arr,:,:)) reinterpret(Int64,view(mat,:,:)) AnyCustomArray

@test_adapt CustomArray reshape(view(mat.arr,:,:), (2,2)) reshape(view(mat,:,:), (2,2)) AnyCustomArray
@test_adapt CustomArray reshape(reinterpret(Int64,mat.arr), (2,2)) reshape(reinterpret(Int64,mat), (2,2)) AnyCustomArray
@test_adapt CustomArray reshape(reinterpret(Int64,view(mat.arr,:,:)), (2,2)) reshape(reinterpret(Int64,view(mat,:,:)), (2,2)) AnyCustomArray

@test_adapt CustomArray view(reinterpret(Int64,mat.arr), :, :) view(reinterpret(Int64,mat), :, :) AnyCustomArray
@test_adapt CustomArray view(reinterpret(Int64,view(mat.arr,:,:)), :, :) view(reinterpret(Int64,view(mat,:,:)), :, :) AnyCustomArray
@test_adapt CustomArray view(Base.ReshapedArray(mat.arr,(2,2),()), :, :) view(reshape(mat, (2,2)), :, :) AnyCustomArray
@test_adapt CustomArray view(reshape(view(mat.arr,:,:), (2,2)), :, :) view(reshape(view(mat,:,:), (2,2)), :, :) AnyCustomArray
@test_adapt CustomArray view(reshape(reinterpret(Int64,mat.arr), (2,2)), :, :) view(reshape(reinterpret(Int64,mat), (2,2)), :, :) AnyCustomArray
@test_adapt CustomArray view(reshape(reinterpret(Int64,view(mat.arr,:,:)), (2,2)), :, :) view(reshape(reinterpret(Int64,view(mat,:,:)), (2,2)), :, :) AnyCustomArray

using LinearAlgebra

@test_adapt CustomArray mat.arr' mat' AnyCustomArray

@test_adapt CustomArray transpose(mat.arr) transpose(mat) AnyCustomArray

@test_adapt CustomArray LowerTriangular(mat.arr) LowerTriangular(mat) AnyCustomArray
@test_adapt CustomArray UnitLowerTriangular(mat.arr) UnitLowerTriangular(mat) AnyCustomArray
@test_adapt CustomArray UpperTriangular(mat.arr) UpperTriangular(mat) AnyCustomArray
@test_adapt CustomArray UnitUpperTriangular(mat.arr) UnitUpperTriangular(mat) AnyCustomArray

@test_adapt CustomArray Diagonal(vec.arr) Diagonal(vec) AnyCustomArray

const dl = CustomArray{Float64,1}(rand(2))
const du = CustomArray{Float64,1}(rand(2))
const d = CustomArray{Float64,1}(rand(3))
@test_adapt CustomArray Tridiagonal(dl.arr, d.arr, du.arr) Tridiagonal(dl, d, du) AnyCustomArray

@test ndims(LinearAlgebra.Transpose{Float64,Array{Float64,1}}) == 2

ext: jl


text:
function plot_LES_3D(frame, FILE_DIR, OUTPUT_PATH, axis_images; title="", colorscheme=:ice, rev=false)
    xy_file = jldopen(joinpath(FILE_DIR, "xy_slice.jld2"))
    xz_file = jldopen(joinpath(FILE_DIR, "xz_slice.jld2"))
    yz_file = jldopen(joinpath(FILE_DIR, "yz_slice.jld2"))
    instantaneous_statistics = jldopen(joinpath(FILE_DIR, "instantaneous_statistics.jld2"))

    frames = keys(xy_file["timeseries/T"])
    iteration = frames[frame]

    color_range_max = maximum(maximum.([xz_file["timeseries/T/$iteration"], yz_file["timeseries/T/$iteration"], xy_file["timeseries/T/$iteration"]]))
    color_range_min = minimum(minimum.([xz_file["timeseries/T/$iteration"], yz_file["timeseries/T/$iteration"], xy_file["timeseries/T/$iteration"]]))
    color_range = (color_range_min, color_range_max)

    xC = xz_file["grid/xC"][4:end-3]
    yC = xz_file["grid/yC"][4:end-3]
    zC = xz_file["grid/zC"][4:end-3]

    colormap = cgrad(colorscheme, scale=:log10, rev=rev)

    fig = Figure(resolution=(1920, 1080))
    ax = fig[1,1] = CairoMakie.Axis3(fig, aspect=(1, 1, 0.5), xlabel="x /m", ylabel="y /m", zlabel="z /m")

    x_xz = fill(xC[1], 128)
    y_xz = yC
    z_xz = zeros(length(x_xz), length(y_xz))

    for i in 1:size(z_xz, 1)
        z_xz[i,:] .= zC[i]
    end

    z_xz
    T_xz = transpose(hcat([xz_file["timeseries/T/$iteration"][:, :, i] for i in 1:length(zC)]... ))

    xz_surface = CairoMakie.surface!(ax, x_xz, y_xz, z_xz, color=T_xz, colormap=colormap)

    x_yz = xC
    y_yz = fill(yC[1], 128)
    z_yz = zeros(length(x_yz), length(y_yz))

    for i in 1:size(z_yz, 2)
        z_yz[:,i] .= zC[i]
    end

    z_yz

    T_yz = similar(z_yz)

    for i in 1:size(T_yz, 2)
        T_yz[:,i] = yz_file["timeseries/T/$iteration"][:, :, i]
    end
    T_yz

    yz_surface = CairoMakie.surface!(ax, x_yz, y_yz, z_yz, color=T_yz, colormap=colormap)

    x_xy = xC
    y_xy = yC
    z_xy = zeros(length(y_xy), length(x_xy))

    T_xy = xy_file["timeseries/T/$iteration"][:,:,1]

    xy_surface = CairoMakie.surface!(ax, x_xy, y_xy, z_xy, color=T_xy, colormap=colormap)

    xz_surface.colorrange = color_range
    yz_surface.colorrange = color_range
    xy_surface.colorrange = color_range

    ax_T_3D = fig[2,1] = CairoMakie.Axis(fig, aspect=DataAspect())

    rel_size = 40
    aspect = 1 / 4

    hidedecorations!(ax_T_3D)
    hidespines!(ax_T_3D)
    image!(ax_T_3D, axis_images.T_3D)
    rowsize!(fig.layout, 2, CairoMakie.Relative(1 / rel_size))

    plots_sublayout = fig[:,2] = GridLayout()

    colsize!(fig.layout, 2, CairoMakie.Relative(aspect))
    rowgap!(fig.layout, Relative(1 / 50 / 3))

    colorbar = CairoMakie.Colorbar(fig[3, 1], xy_surface, vertical=false)

    y_ax_u = plots_sublayout[1,1] = CairoMakie.Axis(fig, aspect=DataAspect())
    y_ax_v = plots_sublayout[3,1] = CairoMakie.Axis(fig, aspect=DataAspect())
    y_ax_T = plots_sublayout[5,1] = CairoMakie.Axis(fig, aspect=DataAspect())

    x_ax_u = plots_sublayout[2,2] = CairoMakie.Axis(fig, aspect=DataAspect())
    x_ax_v = plots_sublayout[4,2] = CairoMakie.Axis(fig, aspect=DataAspect())
    x_ax_T = plots_sublayout[6,2] = CairoMakie.Axis(fig, aspect=DataAspect())

    ax_u = plots_sublayout[1,2] = CairoMakie.Axis(fig)
    ax_v = plots_sublayout[3,2] = CairoMakie.Axis(fig)
    ax_T = plots_sublayout[5,2] = CairoMakie.Axis(fig)

    hidedecorations!(y_ax_u)
    hidedecorations!(y_ax_v)
    hidedecorations!(y_ax_T)

    hidedecorations!(x_ax_u)
    hidedecorations!(x_ax_v)
    hidedecorations!(x_ax_T)

    hidespines!(y_ax_u)
    hidespines!(y_ax_v)
    hidespines!(y_ax_T)

    hidespines!(x_ax_u)
    hidespines!(x_ax_v)
    hidespines!(x_ax_T)

    image!(x_ax_u, axis_images.u)
    image!(x_ax_v, axis_images.v)
    image!(x_ax_T, axis_images.T)

    image!(y_ax_u, axis_images.z)
    image!(y_ax_v, axis_images.z)
    image!(y_ax_T, axis_images.z)

    rowsize!(plots_sublayout, 2, CairoMakie.Relative(1 / rel_size))
    rowsize!(plots_sublayout, 4, CairoMakie.Relative(1 / rel_size))
    rowsize!(plots_sublayout, 6, CairoMakie.Relative(1 / rel_size))
    colsize!(plots_sublayout, 1, CairoMakie.Aspect(2, 1))

    colgap!(plots_sublayout, 1 / 50)
    rowgap!(plots_sublayout, Relative(1 / 50 / 3))

    u_line = CairoMakie.lines!(ax_u, instantaneous_statistics["timeseries/u/$iteration"][:], zC)
    v_line = CairoMakie.lines!(ax_v, instantaneous_statistics["timeseries/v/$iteration"][:], zC)
    T_line = CairoMakie.lines!(ax_T, instantaneous_statistics["timeseries/T/$iteration"][:], zC)

    if title !== ""
        supertitle = fig[0, :] = Label(fig, title, textsize=25)
    end

    trim!(fig.layout)

    close(xy_file)
    close(xz_file)
    close(yz_file)

    save(OUTPUT_PATH, fig, px_per_unit=4, pt_per_unit=2)
end


text:
Vector{UInt8}(
"GET /docs/index.html HTTP/1.1\r\n" *
"Host: www.nowhere123.com\r\n" *
"Accept: image/gif, image/jpeg, */*\r\n" *
"Accept-Language: en-us\r\n" *
"Accept-Encoding: gzip, deflate\r\n" *
"User-Agent: Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1)\r\n" *
"\r\n"
)


text:
function quickb(s::AbstractDiscreteSwitchedSystem, k::Integer=1, clb=true, cub=true)
    lb = 0.
    ub1 = 0.
    ub2 = 0.
    ub∞ = 0.
    for st in states(s)
        for sw in switchings(s, k, st)
            if clb && isperiodic(sw)
                psw = periodicswitching(sw)
                lb = max(lb, psw.growthrate)
                notifyperiodic!(s, psw)
                updatesmp!(s, psw)
            end
            if cub
                ub1 = max(ub1, opnorm(sw.A, 1)^(1/k))
                ub2 = max(ub2, opnorm(sw.A, 2)^(1/k))
                ub∞ = max(ub∞, opnorm(sw.A, Inf)^(1/k))
            end
        end
    end
    ub = min(ub1, ub2, ub∞)
    cub && updateub!(s, ub)
    lb, ub
end

function quickb(s::AbstractDiscreteSwitchedSystem, K::AbstractVector, args...)
    reduce((x, y) -> (max(x[1],y[1]), min(x[2], y[2])), quickb.(s, K, args...), init = (0., Inf))
end

ρA(s::DiscreteSwitchedLinearSystem) = ntransitions(s)

function quicklb(s::DiscreteSwitchedLinearSystem, k=1)
    quickb(s, k, true, false)[1]
end

function quickub(s::AbstractDiscreteSwitchedSystem, k=1)
    quickb(s, k, false, true)[2]
end


text:
show(io::IO, surr::T) where {T<:AbstractSurrogate} =
    print(io, "$(nameof(T))($(surr.ϑ))")

@with_kw_noshow struct Hinge{T<:Real} <: AbstractSurrogate
    ϑ::T = 1.0
    value::Function       = (s) -> max(0, 1 + ϑ*s)
    value_exact::Function = (s) -> Convex.max(0, 1 + ϑ*s)
    gradient::Function    = (s) -> 1 + ϑ*s <= 0 ? zero(ϑ) : ϑ
end

Hinge(ϑ) = Hinge(ϑ = ϑ)

@with_kw_noshow struct Quadratic{T<:Real} <: AbstractSurrogate
    ϑ::T = 1.0
    value::Function       = (s) -> max(0, 1 + ϑ*s)^2
    value_exact::Function = (s) -> Convex.square(Convex.max(0, 1 + ϑ*s))
    gradient::Function    = (s) -> 2*ϑ*max(0, 1 + ϑ*s)
end

Quadratic(ϑ) = Quadratic(ϑ = ϑ)


text:
"""
    RungeKutta3TimeStepper{FT, TG} <: AbstractTimeStepper

Holds parameters and tendency fields for a low storage, third-order Runge-Kutta-Wray
time-stepping scheme described by Le and Moin (1991).
"""
struct RungeKutta3TimeStepper{FT, TG} <: AbstractTimeStepper
    γ¹ :: FT
    γ² :: FT
    γ³ :: FT
    ζ² :: FT
    ζ³ :: FT
    Gⁿ :: TG
    G⁻ :: TG
end

"""
    RungeKutta3TimeStepper(float_type, arch, grid, tracers, χ=0.125;
                              Gⁿ = TendencyFields(arch, grid, tracers),
                              G⁻ = TendencyFields(arch, grid, tracers))

Return an `RungeKutta3TimeStepper` object with tendency fields on `arch` and
`grid`. The tendency fields can be specified via optional kwargs.
"""
function RungeKutta3TimeStepper(float_type, arch, grid, velocities, tracers;
                                Gⁿ = TendencyFields(arch, grid, tracers),
                                G⁻ = TendencyFields(arch, grid, tracers))

    γ¹ = 8 // 15
    γ² = 5 // 12
    γ³ = 3 // 4

    ζ² = -17 // 60
    ζ³ = -5 // 12

    return RungeKutta3TimeStepper{eltype(grid), typeof(Gⁿ)}(γ¹, γ², γ³, ζ², ζ³, Gⁿ, G⁻)
end

#####
##### Time steppping
#####

"""
    time_step!(model::IncompressibleModel{<:RungeKutta3TimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 3rd-order Runge-Kutta method.
The 3rd-order Runge-Kutta method takes three intermediate substep stages to
achieve a single timestep. A pressure correction step is applied at each intermediate
stage.
"""
function time_step!(model::IncompressibleModel{<:RungeKutta3TimeStepper}, Δt)

    γ¹ = model.timestepper.γ¹
    γ² = model.timestepper.γ²
    γ³ = model.timestepper.γ³

    ζ² = model.timestepper.ζ²
    ζ³ = model.timestepper.ζ³

    first_stage_Δt  = γ¹ * Δt
    second_stage_Δt = (γ² + ζ²) * Δt
    third_stage_Δt  = (γ³ + ζ³) * Δt

    arch = model.architecture
    grid = model.grid

    # Convert NamedTuples of Fields to NamedTuples of OffsetArrays
    velocities, tracers, pressures, diffusivities, Gⁿ, G⁻ =
        datatuples(model.velocities, model.tracers, model.pressures, model.diffusivities,
                   model.timestepper.Gⁿ, model.timestepper.G⁻)

    #
    # First stage
    #

    precomputations!(diffusivities, pressures, velocities, tracers, model)

    calculate_tendencies!(Gⁿ, velocities, tracers, pressures, diffusivities, model)

    rk3_substep!(velocities, tracers, arch, grid, Δt, γ¹, nothing, Gⁿ, nothing)

    calculate_pressure_correction!(pressures.pNHS, first_stage_Δt, velocities, model)
    pressure_correct_velocities!(velocities, arch, grid, first_stage_Δt, pressures.pNHS)

    tick!(model.clock, first_stage_Δt; stage=true)

    #
    # Second stage
    #

    precomputations!(diffusivities, pressures, velocities, tracers, model)

    store_tendencies!(G⁻, arch, grid, Gⁿ)
    calculate_tendencies!(Gⁿ, velocities, tracers, pressures, diffusivities, model)

    rk3_substep!(velocities, tracers, arch, grid, Δt, γ², ζ², Gⁿ, G⁻)

    calculate_pressure_correction!(pressures.pNHS, second_stage_Δt, velocities, model)
    pressure_correct_velocities!(velocities, arch, grid, second_stage_Δt, pressures.pNHS)

    tick!(model.clock, second_stage_Δt; stage=true)

    #
    # Third stage
    #

    precomputations!(diffusivities, pressures, velocities, tracers, model)

    store_tendencies!(G⁻, arch, grid, Gⁿ)
    calculate_tendencies!(Gⁿ, velocities, tracers, pressures, diffusivities, model)

    rk3_substep!(velocities, tracers, arch, grid, Δt, γ³, ζ³, Gⁿ, G⁻)

    calculate_pressure_correction!(pressures.pNHS, third_stage_Δt, velocities, model)
    pressure_correct_velocities!(velocities, arch, grid, third_stage_Δt, pressures.pNHS)

    tick!(model.clock, third_stage_Δt)

    return nothing
end

#####
##### Tracer time stepping and predictor velocity updating
#####

function rk3_substep!(U, C, arch, grid, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)

    workgroup, worksize = work_layout(grid, :xyz)

    barrier = Event(device(arch))

    substep_velocities_kernel! = rk3_substep_velocities!(device(arch), workgroup, worksize)
    substep_tracer_kernel! = rk3_substep_tracer!(device(arch), workgroup, worksize)

    velocities_event = substep_velocities_kernel!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻; dependencies=barrier)

    events = [velocities_event]

    for i in 1:length(C)
        @inbounds c = C[i]
        @inbounds Gcⁿ = Gⁿ[i+3]
        @inbounds Gc⁻ = isnothing(G⁻) ? nothing : G⁻[i+3] # so that Gc⁻===nothing for first substep.
        tracer_event = substep_tracer_kernel!(c, Δt, γⁿ, ζⁿ, Gcⁿ, Gc⁻, dependencies=barrier)
        push!(events, tracer_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

"""
Time step tracers via the 3rd-order Runge-Kutta method

    `c^{m+1} = c^m + Δt (γⁿ Gc^{m} + ζⁿ Gc^{m-1})`,

where `m` denotes the substage.
"""
@kernel function rk3_substep_tracer!(c, Δt, γⁿ, ζⁿ, Gcⁿ, Gc⁻)
    i, j, k = @index(Global, NTuple)

    @inbounds c[i, j, k] += Δt * (γⁿ * Gcⁿ[i, j, k] + ζⁿ * Gc⁻[i, j, k])
end

"""
Time step tracers from the first to the second stage via
the 3rd-order Runge-Kutta method

    `c^{2} = c^1 + Δt γ¹ Gc^{1}`.
"""
@kernel function rk3_substep_tracer!(c, Δt, γ¹, ::Nothing, Gc¹, ::Nothing)
    i, j, k = @index(Global, NTuple)

    @inbounds c[i, j, k] += Δt * γ¹ * Gc¹[i, j, k]
end

"""
Time step velocity fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_velocities!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U.u[i, j, k] += Δt * (γⁿ * Gⁿ.u[i, j, k] + ζⁿ * G⁻.u[i, j, k])
        U.v[i, j, k] += Δt * (γⁿ * Gⁿ.v[i, j, k] + ζⁿ * G⁻.v[i, j, k])
        U.w[i, j, k] += Δt * (γⁿ * Gⁿ.w[i, j, k] + ζⁿ * G⁻.w[i, j, k])
    end
end

"""
Time step velocity fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_velocities!(U, Δt, γ¹, ::Nothing, G¹, ::Nothing)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U.u[i, j, k] += Δt * γ¹ * G¹.u[i, j, k]
        U.v[i, j, k] += Δt * γ¹ * G¹.v[i, j, k]
        U.w[i, j, k] += Δt * γ¹ * G¹.w[i, j, k]
    end
end


text:
export qsfitversion, spectral_coverage

qsfitversion() = v"0.1.0"
qsfit_data() = artifact"qsfit_data"

gauss(x, μ, σ) = exp.(-0.5 .* ((x .- μ) ./ σ).^2) ./ sqrt(2pi) ./ σ

function planck(λ, T)
    h = 6.6260755  * 1e-27  # Planck's constant [erg s]
    c = 2.99792458 * 1e10   # Vacuum speed of light [cm / s]
    k = 1.380658   * 1e-16  # Boltzmann constant [erg / K]
    b = 2 * h * c^2. ./ (λ.^5.)
    d = h * c ./ (λ .* k .* T)
    return b ./ (exp.(d) .- 1)
end

function convol(v, _k)
    k = reverse(_k)
    nk = length(k)
    @assert nk < length(v)
    @assert mod(nk, 2) == 1
    r = div(nk, 2)
    out = fill(0., length(v))
    for i in r+1:length(v)-r
        out[i-r:i+r] .+= v[i] .* k
    end
    return out ./ sum(abs.(k))
end

function estimate_fwhm(λ, f)
    i = argmax(f)
    maxv = f[i]
    i = findall(f .>= (maxv/2))
    (i1, i2) = extrema(i)
    fwhm = λ[i2] - λ[i1]
    return fwhm
end

# Wavelength span is assumed to be equal to the initial FWHM plus the
# spectral resolution (in quadrature)
spectral_coverage(spec_λ, resolution, line::Union{SpecLineGauss, SpecLineLorentz, SpecLineVoigt}; kw...) =
    spectral_coverage(spec_λ, resolution, line.center.val, line.center.val * sqrt(line.fwhm.val^2 + resolution^2) / 3.e5; kw...)

function spectral_coverage(spec_λ, resolution, comp::ironuv; kw...)
    rr = extrema(ironuv_read()[1])
    return spectral_coverage(spec_λ, resolution, mean(rr), rr[2]-rr[1]; kw...)
end

function spectral_coverage(spec_λ, resolution, comp::ironopt; kw...)
    rr = extrema(QSFit.ironopt_read(comp.file)[:wavelength])
    return spectral_coverage(spec_λ, resolution, mean(rr), rr[2]-rr[1]; kw...)
end

function spectral_coverage(spec_λ::Vector{Float64}, resolution::Float64,
                           center_λ::Float64, span_λ::Float64; min_steps::Int=5)

    # Identify min/max wavelengths
    λmin = center_λ - span_λ / 2.
    λmax = center_λ + span_λ / 2.

    # Calculate step corresponding to the spectral resolution
    δ = resolution / 3e5 * center_λ

    # How many steps should be considered?
    steps = Int(ceil((λmax - λmin) / δ))
    (steps < min_steps)  &&  (steps = min_steps)
    bin_edges = range(λmin, λmax, length=steps+1)

    # How many intervals are sampled?
    good = 0
    for i in 1:steps
        good += sign(count(bin_edges[i] .<= spec_λ .< bin_edges[i+1]))
    end

    return (λmin, λmax, good / steps)
end


text:
"""

The exported names are
$(EXPORTS)
"""
module CellData

using Test
using DocStringExtensions
using FillArrays

using Gridap.Helpers
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials
using Gridap.Integration
using Gridap.ReferenceFEs
using Gridap.Algebra

import Gridap.Arrays: get_array
import Gridap.Arrays: lazy_append
import Gridap.Fields: evaluate
import Gridap.Fields: gradient
import Gridap.Fields: integrate
import Gridap.Helpers: operate
import Gridap.Integration: get_coordinates
import Gridap.Integration: get_weights
import Gridap.Fields: field_cache
import Gridap.Fields: evaluate_field!
import Gridap.Fields: evaluate_field_array
import Gridap.Arrays: array_cache
import Gridap.Arrays: getindex!
import Gridap.Arrays: reindex

export CellField
export get_cell_map
export get_cell_axes
export get_memo
export RefStyle
export is_in_ref_space
export is_in_physical_space
export MetaSizeStyle
export get_metasize
export is_basis
export is_test
export is_trial
export test_cell_field
export GenericCellField
export similar_object
export change_ref_style
export to_ref_space
export to_physical_space
export trialize_cell_basis
export convert_to_cell_field
export SkeletonCellField
export get_inward
export get_outward
export jump
export mean
export merge_cell_dofs_at_skeleton
export merge_cell_fields_at_skeleton

export CellQuadrature
export get_coordinates
export get_weights

export QPointCellField
export update_state_variables!

export CellDofBasis
export test_cell_dof_basis
export GenericCellDofBasis

export attach_dirichlet

export attach_constraints_rows
export attach_constraints_cols
export merge_cell_constraints_at_skeleton
export identity_constraints

export @law
export operate
export GridapType

include("CellFields.jl")

include("CellQuadratures.jl")

include("QPointCellFields.jl")

include("CellDofBases.jl")

include("AttachDirichlet.jl")

include("AttachConstraints.jl")

include("Law.jl")

end # module


text:
"""
    taccum(x::Incremental)

Returns an accumulated triangle out of an incremental triangle `x`.

# Example
```
julia> ti = Incremental([3 5 2; 4 2 0; 1 0 0])
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0

julia> taccum(ti)
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0
```
"""
function taccum(x::Incremental)
    acum = cumsum(x.claims, dims=2)
    for i in 1:size(acum, 2)
        for j in 1:size(acum, 1)
            if (i + j) > (size(acum, 1) + 1)
                acum[i, j] = 0
            end
        end
    end
    return Accumulated(acum)
end

"""
    tdeacum(x::Accumulated)

Returns an incremental triangle out of an accumulated triangle `x`.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0])
Accumulated Triangle {Int64}:
 3  8  10
 4  6   0
 1  0   0

julia> tdeacum(ta)
Incremental Triangle {Int64}:
 3  5  2
 4  2  0
 1  0  0
```
"""
function tdeacum(x::Accumulated)
    des = zeros(eltype(x.claims), size(x.claims, 1), size(x.claims, 1))
    des[:,1] = x.claims[:,1]
    for j in 2:size(x.claims, 1)
        for i in 1:size(x.claims, 1)
            des[i, j] = x.claims[i, j] - x.claims[i, j - 1]
            if (i + j) > (size(x.claims, 1) + 1)
                des[i, j] = 0
            end
        end
    end
    return Incremental(des)
end

"""
    incurred(x::Accumulated)

Returns the diagonal of an accumulated triangle `x`.

# Example
```
julia> ta = Accumulated([3 8 10; 4 6 0; 1 0 0]);

julia> incurred(ta)
3-element Vector{Any}:
 10
  6
  1
```
"""
function incurred(x::Accumulated)
    incurred = []
    for i in 1:size(x.claims, 1)
        push!(incurred, x.claims[i, size(x.claims,1)-i+1])
    end
    return incurred
end

"""
    randtri(n, d)

Returns a `n`x`n` random matrix, whose values are randomly
taken from the distribution `d`.

# Example
```
julia> d = IbnrClaims.Distributions.Gamma(5,5)
Distributions.Gamma{Float64}(α=5.0, θ=5.0)

julia> randtri(3, d)
3×3 Matrix{Float64}:
 18.162  19.348  29.123
 15.123  28.737   0.0
 18.066   0.0     0.0
```
"""
function randtri(n::Int, d)
    r = zeros(Float64, n, n)
    for i in 1:n
        for j in 1:(n - i + 1)
            r[i,j] = round(rand(d); digits=3)
        end
    end
    return r
end

# Auxiliary function. Helps checking no claims under the diagonal.
function siniestros_bajo_diagonal(t::Matrix)
    q = []
    for i in 1:size(t, 1)
        for j in (size(t, 1)-i+2):(size(t, 2))
            push!(q, t[i,j])
        end
    end
    sum(abs.(q)) != 0
end

# Auxiliary function. Helps calculating loss ratio vector.
function suma_movil_vector(x::Vector, per::Int)
    smv = []
    for i in 1:(length(x)-per+1)
        push!(smv, sum(getindex(x, i:(i+per-1))))
    end
    return smv
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 79168b5e-0805-4413-848a-8c501d81d91c
md"""
## 3.2
 $N$個の整数$a_0, a_1, \dots, a_{N-1}$のうち, 整数値$v$が何個含まれるかを求める$O(N)$のアルゴリズムを設計してください.
"""

# ╔═╡ e76f50bc-5cd6-481f-8fe8-b523534c1e11
function f₂(as, v)
	cnt = 0
	for a ∈ as
		if a == v
			cnt += 1
		end
	end
	return cnt
end

# ╔═╡ add92227-31f6-4408-a527-8f46f586c531
f₂([1 2 3 4 5 3], 3)

# ╔═╡ db1a3c61-5a3d-4431-b3e1-786415736801
f₂([1 2 3 4 5 3], 100)

# ╔═╡ 4349dc9c-cfcc-472e-b611-f87f22c9d4ee
md"""
## 3.3
 $N(\ge 2)$個の相異なる整数$a_0, a_1, \dots, a_{N-1}$が与えられます. このうち2番目に小さい値を求める$O(N)$のアルゴリズムを設計してください.
"""

# ╔═╡ 01f68ae8-e149-44dd-bf11-44e428f29500
function f₃(as)
	max₁, max₂ = maximum(as[1:2]), minimum(as[1:2])
	for a ∈ as
		if a > max₁
			max₂ = max₁
			max₁ = a
		elseif max₁ > a > max₂
			max₂ = a
		end
	end

	return max₂
end

# ╔═╡ b82ccb0e-c956-4e57-9470-7bd5ef4df4b3
f₃([1 2 3 4 5 6])

# ╔═╡ 8df5aae2-5299-40b0-a0db-ac83badd21b1
f₃([1 100 3 4 105 6])

# ╔═╡ f8afba4c-cf89-11eb-0c82-0fe26b9b6ff9
md"""
## 3.4
 $N$個の整数$a_0, a_1, \dots, a_{N-1}$が与えられます. この中から2つ選んで差をとります. その差の最大値を求める$O(N)$のアルゴリズムを設計してください.
"""

# ╔═╡ 6abeb865-69a0-43c9-bff3-03b17d2d6dfe
f₄(as) = maximum(as) - minimum(as)

# ╔═╡ 72bb5735-f769-4043-9305-41dd590f95e5
f₄([1 2 3 4 5 6])

# ╔═╡ f1b037f4-be79-4563-9e5f-97e2a3105d63
f₄([1 100 3 4 105 6])

# ╔═╡ Cell order:
# ╟─79168b5e-0805-4413-848a-8c501d81d91c
# ╠═e76f50bc-5cd6-481f-8fe8-b523534c1e11
# ╠═add92227-31f6-4408-a527-8f46f586c531
# ╠═db1a3c61-5a3d-4431-b3e1-786415736801
# ╟─4349dc9c-cfcc-472e-b611-f87f22c9d4ee
# ╠═01f68ae8-e149-44dd-bf11-44e428f29500
# ╠═b82ccb0e-c956-4e57-9470-7bd5ef4df4b3
# ╠═8df5aae2-5299-40b0-a0db-ac83badd21b1
# ╟─f8afba4c-cf89-11eb-0c82-0fe26b9b6ff9
# ╠═6abeb865-69a0-43c9-bff3-03b17d2d6dfe
# ╠═72bb5735-f769-4043-9305-41dd590f95e5
# ╠═f1b037f4-be79-4563-9e5f-97e2a3105d63


text:
include("../src/LLLplus.jl")
Pkg.add("PyPlot")
Pkg.add("CPUTime")
Pkg.add("ArrayViews")
Pkg.add("Docile")

using PyPlot
using CPUTime
using LLLplus

include("lrtest.jl")

lrtest(100,2.^[4;],[100;],[Int32,Int64,Int128,Float64,BigInt,BigFloat],"rand")
lrtest(200,2.^[4;],[100;],[Int32,Int64,Int128,Float64,BigInt,BigFloat],"rand")
savefig("benchmark/perfVsDataTypeN16.png") # run from root directory

lrtest(100,2.^[0:5;],[1;],[Float64],"randn")
lrtest(200,2.^[0:5;],[1;],[Float64],"randn")
savefig("benchmark/perfVsNfloat32.png") # run from root directory