text:
module Materials

using StaticArrays: SVector

abstract type AbstractMaterial end

abstract type AbstractTransparentMaterial <: AbstractMaterial end
# add absorption etc. methods, they'll be useful

# ============================================================================ #
# ====                               Air                                  ==== #
# ============================================================================ #

struct Air <: AbstractTransparentMaterial end

name(::Air) = :Air

function index(
    ::Air,
    ::T;
    temperature::T = T(20.0),
    pressure::T = T(1.0),
)::T where {T <: Real}
    return one(T)
end

function absorption(
    ::Air,
    ::T;
    temperature::T = T(20.0),
    pressure::T = T(1.0),
)::T where {T <: Real}
    return zero(T)
end

"""λ in m."""
function absairindex(
    λ::T;
    temperature::T = T(20.0),
    pressure::T = T(1.0),
)::T where {T <: Real}
    # convert to required units
    λ *= 1e-6 # convert to um
    n_ref =
        one(T) + (
            (
                6432.8 +
                ((2949810.0 * λ^2) / (146.0 * λ^2 - one(T))) +
                ((25540.0 * λ^2) / (41.0 * λ^2 - one(T)))
            ) * 1e-8
        )
    n_rel =
        one(T) +
        ((n_ref - one(T)) / (one(T) + (temperature - 15.0) * 0.0034785)) *
        pressure
    return n_rel
end

# assume wavelengths are in m
function polyfit_indices(
    wavelengths::Union{AbstractRange{<:Real}, AbstractArray{<:Real, 1}},
    indices::AbstractArray{<:Number, 1};
    degree::Int = 5,
)
    w = 1e-6 .* wavelengths # convert to um
    okay = (indices .> 0.0)
    if !any(okay)
        return (ones(Float64, size(w)) .* NaN, nothing)
    end
    xs = range(-1.0, stop = 1.0, length = length(w[okay]))
    poly = fit(xs, indices[okay], degree)
    interp_indices = poly.(xs)
    # ensure output has all entries
    out = ones(Float64, size(w)) .* NaN
    out[okay] = interp_indices
    return (out, poly)
end

# ============================================================================ #
# ====                          Perfect Abs.                              ==== #
# ============================================================================ #

struct PerfectAbsorber <: AbstractMaterial end

name(::PerfectAbsorber) = :PerfectAbsorber

function absorption(
    ::M,
    ::T;
    temperature::T = T(20.0),
    pressure::T = T(1.0),
)::T where {M <: AbstractMaterial, T <: Real}
    return one(T)
end

# ============================================================================ #
# ====                             Glass                                  ==== #
# ============================================================================ #

struct Glass <: AbstractTransparentMaterial
    name::Symbol
    dispform::Int
    C1::Float64
    C2::Float64
    C3::Float64
    C4::Float64
    C5::Float64
    C6::Float64
    C7::Float64
    C8::Float64
    C9::Float64
    C10::Float64
    λmin::Float64
    λmax::Float64
    D₀::Float64
    D₁::Float64
    D₂::Float64
    E₀::Float64
    E₁::Float64
    λₜₖ::Float64
    temp::Float64
    ΔPgF::Float64
    PR::Float64
    relcost::Float64
    TCE::Float64
    CR::Float64
    status::Int
    SR::Float64
    transmission::Union{Nothing, SVector{100, SVector{3, Float64}}}
    transmissionN::Int
    Nd::Float64
    AR::Float64
    FR::Float64
    exclude_sub::Int
    Vd::Float64
    ignore_thermal_exp::Int
    p::Float64
    meltfreq::Int
    function Glass(
        name::Symbol,
        dispform,
        C1,
        C2,
        C3,
        C4,
        C5,
        C6,
        C7,
        C8,
        C9,
        C10,
        λmin,
        λmax,
        D₀,
        D₁,
        D₂,
        E₀,
        E₁,
        λₜₖ,
        temp,
        ΔPgF,
        PR,
        relcost,
        TCE,
        CR,
        status,
        SR,
        transmission,
        Nd,
        AR,
        FR,
        exclude_sub,
        Vd,
        ignore_thermal_exp,
        p,
        meltfreq,
    )
        # need a constructor to massage the transmission data
        if transmission === nothing
            transmission_s = nothing
            transmissionN = -1
        else
            fill =
                [SVector(0.0, 0.0, 0.0) for _ in 1:(100 - length(transmission))]
            transmission_s =
                SVector{100, SVector{3, Float64}}(transmission..., fill...)
            transmissionN = length(transmission)
        end
        return new(
            name,
            dispform,
            C1,
            C2,
            C3,
            C4,
            C5,
            C6,
            C7,
            C8,
            C9,
            C10,
            λmin,
            λmax,
            D₀,
            D₁,
            D₂,
            E₀,
            E₁,
            λₜₖ,
            temp,
            ΔPgF,
            PR,
            relcost,
            TCE,
            CR,
            status,
            SR,
            transmission_s,
            transmissionN,
            Nd,
            AR,
            FR,
            exclude_sub,
            Vd,
            ignore_thermal_exp,
            p,
            meltfreq,
        )
    end
end

name(g::Glass) = g.name

"""
    info([io::IO], glass::AbstractGlass)

Print out all data associated with `glass` in an easily readable format.

# Examples
```julia-repl
julia> info(GlassCat.RPO.IG4)
name:                                              :AGF_52
Dispersion formula:                                Schott (1)
Dispersion formula coefficients:
     a₀:                                           6.91189161
     a₁:                                           -0.000787956404
     a₂:                                           -4.22296071
     a₃:                                           142.900646
     a₄:                                           -1812.32748
     a₅:                                           7766.33028
Valid wavelengths:                                 3.0μm to 12.0μm
Reference temperature:                              20.0°C
Thermal ΔRI coefficients:
     D₀:                                           3.24e-5
     D₁:                                           0.0
     D₂:                                           0.0
     E₀:                                           0.0
     E₁:                                           0.0
     λₜₖ:                                          0.0
TCE (÷1e-6):                                       20.4
Ignore thermal expansion:                          false
Density (p):                                       4.47g/m³
ΔPgF:                                              0.0
RI at sodium D-Line (587nm):                       1.0
Abbe Number:                                       0.0
Cost relative to N_BK7:                              ?
Status:                                            Standard (0)
Melt frequency:                                    0
Exclude substitution:                              false
```
"""
function info(io::IO, glass::Glass)
    D = glass.dispform

    println(io, "$(rpad("name:", 48)) $(glass.name)")

    if (D == -2)
        println(io, "$(rpad("Dispersion formula:", 50)) Cauchy (-2)")
    elseif (D == -1)
        println(
            io,
            "$(rpad("Dispersion formula:", 50)) Fitted for model/MIL glass",
        )
    else
        println(
            io,
            "$(rpad("Dispersion formula:", 50)) $(DISPFORM_NAMES[D]) ($D)",
        )
    end
    println(io, "Dispersion formula coefficients:")
    if (D == -2) # Cauchy
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
        println(io, "     $(rpad("E:", 45)) $(glass.C5)")
        println(io, "     $(rpad("F:", 45)) $(glass.C6)")
    elseif (D == -1)
        println(io, "     $(rpad("C₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("C₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("C₃:", 45)) $(glass.C4)")
    elseif (D == 1) # Schott
        println(io, "     $(rpad("a₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("a₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("a₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("a₃:", 45)) $(glass.C4)")
        println(io, "     $(rpad("a₄:", 45)) $(glass.C5)")
        println(io, "     $(rpad("a₅:", 45)) $(glass.C6)")
    elseif (D == 2)  # Sellmeier1
        println(io, "     $(rpad("K₁:", 45)) $(glass.C1)")
        println(io, "     $(rpad("L₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("K₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("L₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("K₃:", 45)) $(glass.C5)")
        println(io, "     $(rpad("L₃:", 45)) $(glass.C6)")
    elseif (D == 3)  # Herzberger
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
        println(io, "     $(rpad("E:", 45)) $(glass.C5)")
        println(io, "     $(rpad("F:", 45)) $(glass.C6)")
    elseif (D == 4)  # Sellmeier2
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("λ₁:", 45)) $(glass.C3)")
        println(io, "     $(rpad("B₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("λ₂:", 45)) $(glass.C5)")
    elseif (D == 5)  # Conrady
        println(io, "     $(rpad("n₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("A:", 45)) $(glass.C2)")
        println(io, "     $(rpad("B:", 45)) $(glass.C3)")
    elseif (D == 6)  # Sellmeier3
        println(io, "     $(rpad("K₁:", 45)) $(glass.C1)")
        println(io, "     $(rpad("L₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("K₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("L₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("K₃:", 45)) $(glass.C5)")
        println(io, "     $(rpad("L₃:", 45)) $(glass.C6)")
        println(io, "     $(rpad("K₄:", 45)) $(glass.C7)")
        println(io, "     $(rpad("L₄:", 45)) $(glass.C8)")
    elseif (D == 7) || (D == 8)  # HandbookOfOptics1/2
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
    elseif (D == 9)  # Sellmeier4
        println(io, "     $(rpad("A:", 45)) $(glass.C1)")
        println(io, "     $(rpad("B:", 45)) $(glass.C2)")
        println(io, "     $(rpad("C:", 45)) $(glass.C3)")
        println(io, "     $(rpad("D:", 45)) $(glass.C4)")
        println(io, "     $(rpad("E:", 45)) $(glass.C5)")
    elseif (D == 10) || (D == 12) # Extended1/2
        println(io, "     $(rpad("a₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("a₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("a₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("a₃:", 45)) $(glass.C4)")
        println(io, "     $(rpad("a₄:", 45)) $(glass.C5)")
        println(io, "     $(rpad("a₅:", 45)) $(glass.C6)")
        println(io, "     $(rpad("a₆:", 45)) $(glass.C7)")
        println(io, "     $(rpad("a₇:", 45)) $(glass.C8)")
    elseif (D == 11)  # Sellmeier5
        println(io, "     $(rpad("K₁:", 45)) $(glass.C1)")
        println(io, "     $(rpad("L₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("K₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("L₂:", 45)) $(glass.C4)")
        println(io, "     $(rpad("K₃:", 45)) $(glass.C5)")
        println(io, "     $(rpad("L₃:", 45)) $(glass.C6)")
        println(io, "     $(rpad("K₄:", 45)) $(glass.C7)")
        println(io, "     $(rpad("L₄:", 45)) $(glass.C8)")
        println(io, "     $(rpad("K₅:", 45)) $(glass.C9)")
        println(io, "     $(rpad("L₅:", 45)) $(glass.C10)")
    elseif (D == 13)  # Extended3
        println(io, "     $(rpad("a₀:", 45)) $(glass.C1)")
        println(io, "     $(rpad("a₁:", 45)) $(glass.C2)")
        println(io, "     $(rpad("a₂:", 45)) $(glass.C3)")
        println(io, "     $(rpad("a₃:", 45)) $(glass.C4)")
        println(io, "     $(rpad("a₄:", 45)) $(glass.C5)")
        println(io, "     $(rpad("a₅:", 45)) $(glass.C6)")
        println(io, "     $(rpad("a₆:", 45)) $(glass.C7)")
        println(io, "     $(rpad("a₇:", 45)) $(glass.C8)")
        println(io, "     $(rpad("a₈:", 45)) $(glass.C9)")
    else
        println(io, "     INVALID DISPERSION FORMULA!!")
    end
    println(
        io,
        "$(rpad("Valid wavelengths:", 50)) $(glass.λmin)μm to $(glass.λmax)μm",
    )
    println(io, "$(rpad("Reference temperature:", 50)) $(glass.temp)°C")

    if !isnan(glass.D₀) && (
        glass.D₀ != 0 ||
        glass.D₁ != 0 ||
        glass.D₂ != 0 ||
        glass.E₀ != 0 ||
        glass.E₁ != 0
    )
        println(io, "Thermal ΔRI coefficients:")
        println(io, "     $(rpad("D₀:", 45)) $(glass.D₀)")
        println(io, "     $(rpad("D₁:", 45)) $(glass.D₁)")
        println(io, "     $(rpad("D₂:", 45)) $(glass.D₂)")
        println(io, "     $(rpad("E₀:", 45)) $(glass.E₀)")
        println(io, "     $(rpad("E₁:", 45)) $(glass.E₁)")
        println(io, "     $(rpad("λₜₖ:", 45)) $(glass.λₜₖ)")
    end

    println(io, "$(rpad("TCE (÷1e-6):", 50)) $(glass.TCE)")
    println(
        io,
        "$(rpad("Ignore thermal expansion:", 50)) $(glass.ignore_thermal_exp == 1)",
    )

    println(io, "$(rpad("Density (p):", 50)) $(glass.p)g/m³")
    println(io, "$(rpad("ΔPgF:", 50)) $(glass.ΔPgF)")

    println(io, "$(rpad("RI at sodium D-Line (587nm):", 50)) $(glass.Nd)")
    println(io, "$(rpad("Abbe Number:", 50)) $(glass.Vd)")

    println(
        io,
        "$(rpad("Cost relative to N_BK7:", 50)) $(glass.relcost == -1 ? "?" : glass.relcost)",
    )

    if glass.CR != -1 ||
       glass.FR != -1 ||
       glass.SR != -1 ||
       glass.AR != -1 ||
       glass.PR != -1
        println(io, "Environmental resistance:")
        println(
            io,
            "     $(rpad("Climate (CR):", 45)) $(glass.CR == -1 ? "?" : glass.CR)",
        )
        println(
            io,
            "     $(rpad("Stain (FR):", 45)) $(glass.FR == -1 ? "?" : glass.FR)",
        )
        println(
            io,
            "     $(rpad("Acid (SR):", 45)) $(glass.SR == -1 ? "?" : glass.SR)",
        )
        println(
            io,
            "     $(rpad("Alkaline (AR):", 45)) $(glass.AR == -1 ? "?" : glass.AR)",
        )
        println(
            io,
            "     $(rpad("Phosphate (PR):", 45)) $(glass.PR == -1 ? "?" : glass.PR)",
        )
    end

    println(
        io,
        "$(rpad("Status:", 50)) $(STATUS[glass.status + 1]) ($(glass.status))",
    )
    println(
        io,
        "$(rpad("Melt frequency:", 50)) $(glass.meltfreq == -1 ? "?" : glass.meltfreq)",
    )
    println(
        io,
        "$(rpad("Exclude substitution:", 50)) $(glass.exclude_sub == 1)",
    )

    if glass.transmission !== nothing
        println(io, "Transmission data:")
        println(
            io,
            "$(lpad("Wavelength", 15))$(lpad("Transmission", 15))$(lpad("Thickness", 15))",
        )
        for i in 1:(glass.transmissionN)
            λ, t, τ = glass.transmission[i]
            println(
                io,
                "$(lpad("$(λ)μm", 15))$(lpad(t, 15))$(lpad("$(τ)mm", 15))",
            )
        end
    end
end

NBK7 = Glass(
    :NBK7,
    2,
    1.03961212,
    0.00600069867,
    0.231792344,
    0.0200179144,
    1.01046945,
    103.560653,
    0.0,
    0.0,
    NaN,
    NaN,
    0.3,
    2.5,
    1.86e-6,
    1.31e-8,
    -1.37e-11,
    4.34e-7,
    6.27e-10,
    0.17,
    20.0,
    -0.0009,
    2.3,
    1.0,
    7.1,
    1.0,
    1,
    1.0,
    [
        (0.3, 0.05, 25.0),
        (0.31, 0.25, 25.0),
        (0.32, 0.52, 25.0),
        (0.334, 0.78, 25.0),
        (0.35, 0.92, 25.0),
        (0.365, 0.971, 25.0),
        (0.37, 0.977, 25.0),
        (0.38, 0.983, 25.0),
        (0.39, 0.989, 25.0),
        (0.4, 0.992, 25.0),
        (0.405, 0.993, 25.0),
        (0.42, 0.993, 25.0),
        (0.436, 0.992, 25.0),
        (0.46, 0.993, 25.0),
        (0.5, 0.994, 25.0),
        (0.546, 0.996, 25.0),
        (0.58, 0.995, 25.0),
        (0.62, 0.994, 25.0),
        (0.66, 0.994, 25.0),
        (0.7, 0.996, 25.0),
        (1.06, 0.997, 25.0),
        (1.53, 0.98, 25.0),
        (1.97, 0.84, 25.0),
        (2.325, 0.56, 25.0),
        (2.5, 0.36, 25.0),
    ],
    1.5168,
    2.3,
    0.0,
    0,
    64.17,
    0,
    2.51,
    0.0,
)

"""
    modelglass(Nd::Float64, Vd::Float64, ΔPgF::Float64, name::Symbol) -> Glass

Generates a glass object for the given refractive index at d-light (587.5618nm),
 `Nd`, the Abbe number also at d-light, `Vd`, and partial dispersion, `ΔPgF`.
The mean error to measured data for these models is typically small -
usually < 0.0001.
Behavior may differ from other optical simulation tools when using model glasses.

The approximate dispersion calculation used for these glasses is generally only
valid for visible wavelengths, in this case a limit of 360nm to 750nm is imposed.

# Examples
```julia-repl
julia> index(modelglass(1.5168, 64.17, 0.0), 0.5875618)
1.5168003970108495

julia> index(modelglass(1.2344, 61.57, 0.003), 0.678)
1.2329425902693352
```
"""
function modelglass(Nd::Float64, Vd::Float64, ΔPgF::Float64, name::Symbol)
    # from Schott "TIE-29: Refractive Index and Dispersion"
    a = ΔPgF + 0.6438 - 0.001682 * Vd
    # Using fitting results from https://www.gnu.org/software/goptical/manual/Material_Abbe_class_reference.html
    C1 = a * -6.11873891971188577088 + 1.17752614766485175224
    C2 = a * 18.27315722388047447566 + -8.93204522498095698779
    C3 = a * -14.55275321129051135927 + 7.91015964461522003148
    C4 = a * 3.48385106908642905310 + -1.80321117937358499361
    return Glass(
        name,
        -1,
        C1,
        C2,
        C3,
        C4,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.36,
        0.75,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        TEMP_REF,
        ΔPgF,
        -1.0,
        -1.0,
        0.0,
        -1.0,
        0,
        -1.0,
        nothing,
        Nd,
        -1.0,
        -1.0,
        0,
        Vd,
        1,
        0.0,
        0.0,
    )
end

"""
    glassfromMIL(glasscode::Union{Float64,Int}, name::Symbol) -> Glass

Generates a glass object for the given glass code based on U.S. military
standard MIL-G-174, see
[the MIL specification](http://www.newportglass.com/GeneCd.htm)
for further details.

The glass code is a six-digit number specifying the glass according to its
refractive index `Nd` at d-light (587.5618nm), and its Abbe number `Vd` also
taken at d-light.
The resulting glass code is the value of `Nd - 1` rounded to three digits,
followed by `Vd` rounded to three digits, with all decimal points ignored.
For example, N_BK7 has `Nd = 1.5168` and `Vd = 64.17`, giving a six-digit glass
code of `517642`.

For `Nd > 1.999` the format `1.123642` can be used representing `Nd = 2.123` and
`Vd = 64.2`.

**Accuracy is poor given the low precision of the input parameters**, the mean
error to measured data may be significant.
Behavior may differ from other optical simulation tools when using MIL glasses.
The approximate dispersion calculation used these glasses is generally only
valid  for visible wavelengths, in this case a limit of 360nm to 750nm is
imposed.

# Examples
```julia-repl
julia> index(glassfromMIL(517642), 0.5875618)
1.5170003960064509

julia> index(glassfromMIL(1.134642), 0.5875618)
2.1340008686098946
```
"""
function glassfromMIL(glasscode::Int, name::Symbol)
    Nd = floor(Int, glasscode / 1000) / 1000 + 1
    Vd = (glasscode - floor(Int, glasscode / 1000) * 1000) / 10
    g = modelglass(Nd, Vd, 0.0, name)
    return g
end

function glassfromMIL(glasscode::Float64, name::Symbol)
    @assert glasscode > 1.0
    glasscodeid = round(Int, glasscode * 1000000)
    Nd = floor(Int, glasscode * 1000) / 1000 + 1
    Vd = round(
        (glasscode * 1000 - floor(Int, glasscode * 1000)) * 100,
        digits = 1,
    )
    g = modelglass(Nd, Vd, 0.0, name)
    return g
end

info(g::Glass) = info(stdout, g)

"""
    absorption(glass::Glass, wavelength; temperature=20.0, pressure=1.0)

Compute the intensity absorption per mm of `glass` at `wavelength`,
optionally at specified `temperature` and `pressure`.
Transmission values are linearly interpolated from the adjacent values in the
data table of `glass`, if `wavelength` is below the minimum or above the
maximum in the table then the nearest value is taken.

Absorption is defined as ``\\frac{-\\log(t)}{\\tau}`` where ``t`` is the
transmission value and ``\\tau`` is the thickness, both of which are provided
in the data table.

Wavelength, temp and pressure in m, °C and Atm.
"""
function absorption(
    glass::Glass,
    λ::T;
    temperature::T = T(20.0),
    pressure::T = T(1.0),
)::T where {T <: Real}
    # if the glass has no transmission data then assume no absorption
    if glass.transmission === nothing
        return zero(T)
    end

    reference_temp = T(glass.temp)

    # to work out the wavelength at the reference temperature we need the RIs of air at system temp and at reference temp
    n_air_at_sys =
        absairindex(λ, temperature = temperature, pressure = pressure)
    n_air_at_ref = absairindex(λ, temperature = reference_temp)

    # scale the wavelength to air at the reference temperature/pressure (and um)
    λ = λ * 1e-6 * (n_air_at_sys / n_air_at_ref)

    tdata = glass.transmission
    N = glass.transmissionN
    if λ < tdata[1][1]
        t = tdata[1][2]
        τ = tdata[1][3]
        return T(-log1p(t - 1.0) / τ)
    elseif λ > tdata[N][1]
        t = tdata[N][2]
        τ = tdata[N][3]
        return T(-log1p(t - 1.0) / τ)
    else
        let λlow = 0.0,
            tlow = 0.0,
            τlow = 0.0,
            λhigh = 0.0,
            thigh = 0.0,
            τhigh = 0.0

            for i in 2:N
                if λ <= tdata[i][1]
                    λlow, tlow, τlow = tdata[i - 1]
                    λhigh, thigh, τhigh = tdata[i]
                    break
                end
            end
            λhigh = T(λhigh)
            λlow = T(λlow)
            δλ = λhigh - λlow
            @assert τlow == τhigh
            t = (tlow * (λhigh - λ) / δλ) + (thigh * (λ - λlow) / δλ)
            return -log1p(t - 1.0) / τhigh
        end
    end
end

"""
    index(glass::AbstractGlass, wavelength; temperature=20°C, pressure=1Atm)

Compute the refractive index of `glass` at `wavelength`, optionally at specified `temperature` and `pressure`.
Result is relative to the refractive index of air at given temperature and pressure.

Arguments are interpretted as m, °C and Atm respectively.

**This is defined to always equal 1.0 for Air at any temperature and pressure**,
 use [`absairindex`](@ref) for the absolute refractive index of air at a given
 temperature and pressure.
"""
function index(
    glass::Glass,
    λ::T;
    temperature::T = T(20.0),
    pressure::T = T(1.0),
)::T where {T <: Real}
    # all calculations for the material must be done at the refernce temperature
    reference_temp = T(glass.temp)

    # to work out the wavelength at the reference temperature we need the RIs of air at system temp and at reference temp
    n_air_at_sys =
        absairindex(λ, temperature = temperature, pressure = pressure)
    n_air_at_ref = absairindex(λ, temperature = reference_temp)

    # scale the wavelength to air at the reference temperature/pressure
    λabs = λ * n_air_at_sys * 1e-6 # (and um)
    λ = λabs / n_air_at_ref

    if (λ < glass.λmin) || (λ > glass.λmax)
        error(
            "Cannot calculate an index for the specified wavelength: $λ, valid range: [$(glass.λmin), $(glass.λmax)].\n",
        )
    end

    if glass.dispform == -2
        # Cauchy
        n_rel =
            T(glass.C1) +
            (glass.C2 * λ^(-2)) +
            (glass.C3 * λ^(-4)) +
            (glass.C4 * λ^(-6)) +
            (glass.C5 * λ^(-8)) +
            (glass.C6 * λ^(-10))
    elseif glass.dispform == -1
        # use fitted result from GOptical:
        n_rel = T(
            glass.Nd +
            (glass.Nd - one(T)) / glass.Vd *
            (glass.C1 + glass.C2 / λ + glass.C3 / λ^2 + glass.C4 / λ^3),
        )
    elseif glass.dispform == 1
        # Schott
        formula_rhs =
            T(glass.C1) +
            (glass.C2 * λ^2) +
            (glass.C3 * λ^(-2)) +
            (glass.C4 * λ^(-4)) +
            (glass.C5 * λ^(-6)) +
            (glass.C6 * λ^(-8))
        n_rel = sqrt(formula_rhs)
    elseif glass.dispform == 2
        # Sellmeier1
        formula_rhs =
            (glass.C1 * λ^2 / (λ^2 - glass.C2)) +
            (glass.C3 * λ^2 / (λ^2 - glass.C4)) +
            (glass.C5 * λ^2 / (λ^2 - glass.C6))
        n_rel = sqrt(formula_rhs + one(T))
    elseif glass.dispform == 3
        # Herzberger
        L = one(T) / (λ^2 - T(0.028))
        n_rel =
            T(glass.C1) +
            (glass.C2 * L) +
            (glass.C3 * L^2) +
            (glass.C4 * λ^2) +
            (glass.C5 * λ^4) +
            (glass.C6 * λ^6)
    elseif glass.dispform == 4
        # Sellmeier2
        formula_rhs =
            T(glass.C1) +
            (glass.C2 * λ^2 / (λ^2 - (glass.C3)^2)) +
            (glass.C4 * λ^2 / (λ^2 - (glass.C5)^2))
        n_rel = sqrt(formula_rhs + one(T))
    elseif glass.dispform == 5
        # Conrady
        n_rel = T(glass.C1) + (glass.C2 / λ) + (glass.C3 / λ^3.5)
    elseif glass.dispform == 6
        # Sellmeier3
        formula_rhs =
            (glass.C1 * λ^2 / (λ^2 - glass.C2)) +
            (glass.C3 * λ^2 / (λ^2 - glass.C4)) +
            (glass.C5 * λ^2 / (λ^2 - glass.C6)) +
            (glass.C7 * λ^2 / (λ^2 - glass.C8))
        n_rel = sqrt(formula_rhs + one(T))
    elseif glass.dispform == 7
        # HandbookOfOptics1
        formula_rhs =
            T(glass.C1) + (glass.C2 / (λ^2 - glass.C3)) - (glass.C4 * λ^2)
        n_rel = sqrt(formula_rhs)
    elseif glass.dispform == 8
        # HandbookOfOptics2
        formula_rhs =
            T(glass.C1) + (glass.C2 * λ^2 / (λ^2 - glass.C3)) - (glass.C4 * λ^2)
        n_rel = sqrt(formula_rhs)
    elseif glass.dispform == 9
        # Sellmeier4
        formula_rhs =
            T(glass.C1) +
            (glass.C2 * λ^2 / (λ^2 - glass.C3)) +
            (glass.C4 * λ^2 / (λ^2 - glass.C5))
        n_rel = sqrt(formula_rhs)
    elseif glass.dispform == 10
        # Extended1
        formula_rhs =
            T(glass.C1) +
            (glass.C2 * λ^2) +
            (glass.C3 * λ^(-2)) +
            (glass.C4 * λ^(-4)) +
            (glass.C5 * λ^(-6)) +
            (glass.C6 * λ^(-8)) +
            (glass.C7 * λ^(-10)) +
            (glass.C8 * λ^(-12))
        n_rel = sqrt(formula_rhs)
    elseif glass.dispform == 11
        # Sellmeier5
        formula_rhs =
            (glass.C1 * λ^2 / (λ^2 - glass.C2)) +
            (glass.C3 * λ^2 / (λ^2 - glass.C4)) +
            (glass.C5 * λ^2 / (λ^2 - glass.C6)) +
            (glass.C7 * λ^2 / (λ^2 - glass.C8)) +
            (glass.C9 * λ^2 / (λ^2 - glass.C10))
        n_rel = sqrt(formula_rhs + one(T))
    elseif glass.dispform == 12
        # Extended2
        formula_rhs =
            T(glass.C1) +
            (glass.C2 * λ^2) +
            (glass.C3 * λ^(-2)) +
            (glass.C4 * λ^(-4)) +
            (glass.C5 * λ^(-6)) +
            (glass.C6 * λ^(-8)) +
            (glass.C7 * λ^4) +
            (glass.C8 * λ^6)
        n_rel = sqrt(formula_rhs)
    elseif glass.dispform == 13
        # Extended3
        formula_rhs =
            T(glass.C1) +
            (glass.C2 * λ^2) +
            (glass.C3 * λ^(4)) +
            (glass.C4 * λ^(-2)) +
            (glass.C5 * λ^(-4)) +
            (glass.C6 * λ^(-6)) +
            (glass.C7 * λ^(-8)) +
            (glass.C8 * λ^(-10)) +
            (glass.C9 * λ^(-12))
        n_rel = sqrt(formula_rhs)
    else
        @error "Invalid glass dispersion formula"
    end

    # get the absolute index of the material
    n_abs = n_rel * n_air_at_ref

    # If "TD" is included in the glass data, then include pressure and temperature dependence of the lens
    # environment. From Schott"s technical report "TIE-19: Temperature Coefficient of the Refractive Index".
    # The above "n_rel" data are assumed to be from the reference temperature T_ref. Now we add a small change
    # delta_n to it due to a change in temperature.
    ΔT = temperature - reference_temp
    if !isnan(glass.D₀) &&
       abs(ΔT) > 0.0 &&
       (
           glass.D₀ != 0 ||
           glass.D₁ != 0 ||
           glass.D₂ != 0 ||
           glass.E₀ != 0 ||
           glass.E₁ != 0
       )
        Sₜₖ = glass.λₜₖ < 0.0 ? -one(T) : one(T)
        Δn_abs =
            ((n_rel^2 - one(T)) / (2.0 * n_rel)) * (
                glass.D₀ * ΔT +
                glass.D₁ * ΔT^2 +
                glass.D₂ * ΔT^3 +
                ((glass.E₀ * ΔT + glass.E₁ * ΔT^2) / (λ^2 - Sₜₖ * glass.λₜₖ^2))
            )
        n_abs = n_abs + Δn_abs
    end

    # make the index relative to the RI of the air at the system temperature/pressure again
    n_rel = n_abs / n_air_at_sys
    return n_rel
end

# ============================================================================ #
# ====                        Material Table                              ==== #
# ============================================================================ #

MaterialTable = Dict{Symbol, AbstractMaterial}(
    :Air => Air(),
    :PerfectAbsorber => PerfectAbsorber(),
    :NBK7 => NBK7,
)

function registerMaterial(name::Symbol, mat::T) where {T <: AbstractMaterial}
    if name in MaterialTable
        error("name '$name' already registered.")
    end
    return MaterialTable[name] = mat
end

end # module


text:
center(w) = flex_row(w)
center(w::Widget) = w
center(w::Widget{:toggle}) = flex_row(w)

manipulatelayout(::WidgetTheme) = t -> node(:div, map(center, values(components(t)))..., map(center, t.output))

function widget(::WidgetTheme, x::Observable; label = nothing)
    if label === nothing
        x
    else
        Widget{:observable}(["label" => label], output = x, layout = t -> flex_row(t["label"], t.output))
    end
end

ext: jl


text:
# Independent Component Analysis

#### FastICA type

mutable struct ICA{T<:Real}
    mean::Vector{T}   # mean vector, of length m (or empty to indicate zero mean)
    W::Matrix{T}      # component coefficient matrix, of size (m, k)
end

indim(M::ICA) = size(M.W, 1)
outdim(M::ICA) = size(M.W, 2)
Base.mean(M::ICA) = fullmean(indim(M), M.mean)

transform(M::ICA, x::AbstractVecOrMat) = At_mul_B(M.W, centralize(x, M.mean))

#### core algorithm

# the abstract type for all g functions:
#
# Let f be an instance of such type, then
#
#   evaluate(f, x) --> (v, d)
#
# It returns a function value v, and derivative d
#
abstract type ICAGDeriv{T<:Real} end

struct Tanh{T} <: ICAGDeriv{T}
    a::T
end

evaluate(f::Tanh{T}, x::T) where {T<:Real} = (a = f.a; t = tanh(a * x); (t, a * (1 - t * t)))

struct Gaus{T} <: ICAGDeriv{T} end
evaluate(f::Gaus{T}, x::T) where {T<:Real} = (x2 = x * x; e = exp(-x2/2); (x * e, (1 - x2) * e))

## a function to get a g-fun

icagfun(fname::Symbol, ::Type{T} = Float64) where {T<:Real}=
    fname == :tanh ? Tanh{T}(1.0) :
    fname == :gaus ? Gaus{T}() :
    error("Unknown gfun $(fname)")

icagfun(fname::Symbol, a::T) where {T<:Real} =
    fname == :tanh ? Tanh(a) :
    fname == :gaus ? error("The gfun $(fname) has no parameters") :
    error("Unknown gfun $(fname)")

# Fast ICA
#
# Reference:
#
#   Aapo Hyvarinen and Erkki Oja
#   Independent Component Analysis: Algorithms and Applications.
#   Neural Network 13(4-5), 2000.
#
function fastica!(W::DenseMatrix{T},      # initialized component matrix, size (m, k)
                  X::DenseMatrix{T},      # (whitened) observation sample matrix, size(m, n)
                  fun::ICAGDeriv{T},      # approximate neg-entropy functor
                  maxiter::Int,           # maximum number of iterations
                  tol::Real,              # convergence tolerance
                  verbose::Bool) where {T<:Real}          # whether to show iterative info

    # argument checking
    m = size(W, 1)
    k = size(W, 2)
    size(X, 1) == m || throw(DimensionMismatch("Sizes of W and X mismatch."))
    n = size(X, 2)
    k <= min(m, n) || throw(DimensionMismatch("k must not exceed min(m, n)."))

    if verbose
        @printf("FastICA Algorithm (m = %d, n = %d, k = %d)\n", m, n, k)
        println("============================================")
    end

    # pre-allocated storage
    Wp = Matrix{T}(m, k)    # to store previous version of W
    U  = Matrix{T}(n, k)    # to store w'x & g(w'x)
    Y  = Matrix{T}(m, k)    # to store E{x g(w'x)} for components
    E1 = Vector{T}(k)       # store E{g'(w'x)} for components

    # normalize each column
    for j = 1:k
        w = view(W,:,j)
        scale!(w, 1.0 / sqrt(sum(abs2, w)))
    end

    # main loop
    chg = NaN
    t = 0
    converged = false
    while !converged && t < maxiter
        t += 1
        copy!(Wp, W)

        # apply W of previous step
        At_mul_B!(U, X, W)  # u <- w'x

        # compute g(w'x) --> U and E{g'(w'x)} --> E1
        _s = 0.0
        for j = 1:k
            for i = 1:n
                u, v = evaluate(fun, U[i,j])
                U[i,j] = u
                _s += v
            end
            E1[j] = _s / n
        end

        # compute E{x g(w'x)} --> Y
        scale!(A_mul_B!(Y, X, U), 1.0 / n)

        # update w: E{x g(w'x)} - E{g'(w'x)} w := y - e1 * w
        for j = 1:k
            w = view(W,:,j)
            y = view(Y,:,j)
            e1 = E1[j]
            for i = 1:m
                w[i] = y[i] - e1 * w[i]
            end
        end

        # symmetric decorrelation: W <- W * (W'W)^{-1/2}
        copy!(W, W * _invsqrtm!(W'W))

        # compare with Wp
        chg = 0.0
        for j = 1:k
            s = 0.0
            w = view(W,:,j)
            wp = view(Wp,:,j)
            for i = 1:m
                s += abs(w[i] - wp[i])
            end
            if s > chg
                chg = s
            end
        end
        converged = (chg < tol)

        if verbose
            @printf("Iter %4d:  change = %.6e\n", t, chg)
        end
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))
    return W
end

#### interface function

function fit(::Type{ICA}, X::DenseMatrix{T},                # sample matrix, size (m, n)
                          k::Int;                           # number of independent components
                          alg::Symbol=:fastica,             # choice of algorithm
                          fun::ICAGDeriv=icagfun(:tanh, T), # approx neg-entropy functor
                          do_whiten::Bool=true,             # whether to perform pre-whitening
                          maxiter::Integer=100,             # maximum number of iterations
                          tol::Real=1.0e-6,                 # convergence tolerance
                          mean=nothing,                     # pre-computed mean
                          winit::Matrix{T}=zeros(T,0,0),    # init guess of W, size (m, k)
                          verbose::Bool=false) where {T<:Real}             # whether to display iterations

    # check input arguments
    m, n = size(X)
    n > 1 || error("There must be more than one samples, i.e. n > 1.")
    k <= min(m, n) || error("k must not exceed min(m, n).")

    alg == :fastica || error("alg must be :fastica")
    maxiter > 1 || error("maxiter must be greater than 1.")
    tol > 0 || error("tol must be positive.")

    # preprocess data
    mv = preprocess_mean(X, mean)
    Z::Matrix{T} = centralize(X, mv)

    W0= zeros(T, 0,0)    # whitening matrix
    if do_whiten
        C = scale!(A_mul_Bt(Z, Z), 1.0 / (n - 1))
        Efac = eigfact(C)
        ord = sortperm(Efac.values; rev=true)
        (v, P) = extract_kv(Efac, ord, k)
        W0 = scale!(P, 1 ./ sqrt.(v))
        # println(W0' * C * W0)
        Z = W0'Z
    end

    # initialize
    W = (isempty(winit) ? randn(T, size(Z,1), k) : copy(winit))

    # invoke core algorithm
    fastica!(W, Z, fun, maxiter, tol, verbose)

    # construct model
    if do_whiten
        W = W0 * W
    end
    return ICA(mv, W)
end


text:
using PSSFSS
using PSSFSS.Elements: s₁s₂2β₁β₂
using PSSFSS.Zint: filljk!, vtxcrd
using PSSFSS.RWG: setup_rwg
#using StaticArrays
using Cubature
using LinearAlgebra: ⋅, norm
using Test

"""
    duffy(f,reltol=1e-10)

Perform the integral  ``\\int_0^1 \\int_0^{1-y} f(x,y) \\, dx \\, dy`` via a call to `pcubature` of the
`Cubature` package using a Duffy transform.
"""
function duffy(f,reltol=1e-10)
    a = [0.,0.]
    b = [1., 1.]
    hcubature(x -> (1-x[2]) * f([(1-x[2])*x[1],x[2]]), a, b; reltol)
end

"""
    duffy2(f,reltol=1e-10)

Perform the integral  ``\\int_0^1 \\int_0^{1-y} f(x,y) \\, dx \\, dy`` via a call to `hcubature` of the
`Cubature` package using a Duffy transform.  Here `f(ξ,η)` returns a 2-vector.
"""
function duffy2(f,reltol=1e-10)
    a = [0.,0.]
    b = [1., 1.]
    hcubature(2, (x,v) -> v[:] = (1-x[2]) * f([(1-x[2])*x[1],x[2]]), a, b; reltol)
end

function fρ_r(ξ,η,rm,rs)
    r = rs[1] + ξ*(rs[2]-rs[1]) + η*(rs[3]-rs[1]) - rm
    return r/norm(r)
end

function frinv(ξ,η,rm,rs)
    r = rs[1] + ξ*(rs[2]-rs[1]) + η*(rs[3]-rs[1]) - rm
    return 1/norm(r)
end

function twice_area(r1,r2,r3)
    r21 = r2-r1
    r31 = r3-r1
    return abs(r21[1]*r31[2] - r21[2]*r31[1])
end
twice_area(rvec) = twice_area(rvec[1],rvec[2],rvec[3])

sheet = polyring(ntri=100, sides=32, a=[0.4957], b=[0.5557],
                 s1=[1.1414, 0.0], s2=[0.5707, 0.9885], units=inch)
rwgdata = setup_rwg(sheet)
nface = size(sheet.fv,2)
nufp = rwgdata.nufp
i2s = CartesianIndices((nface,nface))
ψ₁, ψ₂ = 0., 0.  # Incremental phase shifts (rad)
s₁, s₂ = sheet.s₁, sheet.s₂ # Lattice vectors (mil)
β₁, β₂ = s₁s₂2β₁β₂(s₁,s₂) # Reciprocal lattice vectors (1/mil)
u = 0.4 * max(norm(β₁), norm(β₂)) #  smoothing parameter (1/mil)
sheet.u = u
sheet.ψ₁ = ψ₁
sheet.ψ₂ = ψ₂
extract = true
filljk!(sheet, rwgdata, extract)
conv = 1e-7

@testset "filljk! Test" begin
    for iufp = [1,2,3, nufp-3, nufp-2, nufp-1, nufp]
        ifmifs = rwgdata.ufp2fp[iufp][1]  # Obtain index into face/face matrix
        obs_src = i2s[ifmifs]
        ifm, ifs = obs_src[1], obs_src[2] # face indices: match and source
        rs = vtxcrd(ifs, sheet) # Obtain coordinates of source tri's vertices
        area2 = twice_area(rs)
        rm = sum(vtxcrd(ifm, sheet))/3 # Obtain coordinates of match tri's centroid
        if ifm == ifs
            # Singularity must occur at integration limit for hcubature
            rsnew = [rm, rs[1], rs[2]]
            a1 = twice_area(rsnew)
            rinv1, err_rinv1 = duffy(ξη -> frinv(ξη[1],ξη[2],rm,rsnew))
            ρ_r1, err_ρ1 = duffy2(ξη -> fρ_r(ξη[1],ξη[2],rm,rsnew))
            rsnew = [rm, rs[1], rs[3]]
            a2 = twice_area(rsnew)
            rinv2, err_rinv2 = duffy(ξη -> frinv(ξη[1],ξη[2],rm,rsnew))
            ρ_r2, err_ρ2 = duffy2(ξη -> fρ_r(ξη[1],ξη[2],rm,[rm, rs[1], rs[3]]))
            rsnew = [rm, rs[2], rs[3]]
            a3 = twice_area(rsnew)
            rinv3, err_rinv3 = duffy(ξη -> frinv(ξη[1],ξη[2],rm,rsnew))
            ρ_r3, err_ρ3 = duffy2(ξη -> fρ_r(ξη[1],ξη[2],rm,rsnew))
            rinv_duffy = (a1*rinv1 + a2*rinv2 + a3*rinv3) / area2
            ρ_r_duffy = (a1*ρ_r1 + a2*ρ_r2 + a3*ρ_r3) / area2
        else
            rinv_duffy, err_rinv = duffy(ξη -> frinv(ξη[1],ξη[2],rm,rs))
            ρ_r_duffy, err_ρ_r = duffy2(ξη -> fρ_r(ξη[1],ξη[2],rm,rs))
        end
        @test abs((u*sheet.rinv[iufp] - rinv_duffy) / rinv_duffy) < conv
        @test norm((sheet.ρ_r[iufp] - ρ_r_duffy) / ρ_r_duffy) < conv
    end
end


text:
# some views are in colorchannels.jl
using Colors, FixedPointNumbers, ImageCore, OffsetArrays, Test
using OffsetArrays: IdentityUnitRange
using PaddedViews: filltype

@testset "rawview" begin
    a = map(N0f8, rand(3,5))
    a[2,2] = N0f8(0.5)
    v = rawview(a)
    @test v[2,2] === a[2,2].i
    v[1,3] = 0xff
    @test a[1,3] === N0f8(1)
    v[1,3] = 0x01
    @test a[1,3] === N0f8(1/255)
    s = view(a, 1:2, 1:2)
    v = rawview(s)
    @test v[2,2] === a[2,2].i
    v[2,2] = 0x0f
    @test a[2,2].i == 0x0f
    @test rawview(v) === v
end

@testset "normedview" begin
    a = rand(UInt8, 3, 5)
    a[2,2] = 0x80
    v = normedview(a)
    @test v[2,2] === N0f8(0.5)
    v[1,3] = 1
    @test a[1,3] === 0xff
    v[1,3] = 1/255
    @test a[1,3] === 0x01
    s = view(a, 1:2, 1:2)
    v = normedview(s)
    @test v[2,2] === N0f8(0.5)
    v[2,2] = 15/255
    @test a[2,2] == 0x0f
    @test normedview(v) === v
    @test normedview(N0f8, v) === v
end

@testset "StackedView" begin
    for (A, B, T) = (([1 3;2 4], [-1 -5; -2 -3], Int),
                     ([1 3;2 4], [-1.0 -5.0; -2.0 -3.0], Float64))
        V = @inferred(StackedView(A, B))
        @test eltype(V) == T
        @test size(V) == (2, 2, 2)
        @test axes(V) === (Base.OneTo(2), Base.OneTo(2), Base.OneTo(2))
        @test @inferred(V[1,1,1]) === T(1)
        @test @inferred(V[2,1,1]) === T(-1)
        @test V[1,:,:] == A
        @test V[2,:,:] == B
        @test_throws BoundsError V[0,1,1]
        @test_throws BoundsError V[3,1,1]
        @test_throws BoundsError V[2,0,1]
        @test_throws BoundsError V[2,3,1]
        @test_throws BoundsError V[1,1,0]
        @test_throws BoundsError V[1,1,3]
        V32 = @inferred(StackedView{Float32}(A, B))
        @test eltype(V32) == Float32
        @test V32[1,1,2] == Float32(3)
        V[1,2,2] = 0
        @test A[2,2] == 0
        V[2,1,2] = 11
        @test B[1,2] == 11

        V = @inferred(StackedView(A, zeroarray, B))
        @test eltype(V) == T
        @test size(V) == (3, 2, 2)
        @test axes(V) === (Base.OneTo(3), Base.OneTo(2), Base.OneTo(2))
        @test V[1,:,:] == A
        @test all(iszero, V[2,:,:])
        @test V[3,:,:] == B
        @test_throws ErrorException V[2,1,1] = 7
        V32 = @inferred(StackedView{Float32}(A, zeroarray, B))
        @test eltype(V32) == Float32
        @test V32[1,1,2] == Float32(3)
    end

    # With mixed grayscale/real arrays
    a, b = Gray{N0f8}[0.1 0.2; 0.3 0.4], [0.5 0.6; 0.7 0.8]
    V = @inferred(StackedView(a, b))
    @test eltype(V) == Float64
    @test V[1,1,1] === Float64(N0f8(0.1))
    @test V[2,1,1] === 0.5
    V = @inferred(StackedView{N0f8}(a, b))
    @test eltype(V) == N0f8
    @test V[1,1,1] === N0f8(0.1)
    @test V[2,1,1] === N0f8(0.5)
    @test b[1,1] === 0.5
    V = @inferred(StackedView(a, zeroarray, b))
    @test eltype(V) == Float64
    V = @inferred(StackedView{N0f8}(a, zeroarray, b))
    @test eltype(V) == N0f8

    # With colorview
    a = [0.1 0.2; 0.3 0.4]
    b = [0.5 0.6; 0.7 0.8]
    z = zeros(2,2)  # because setindex! won't work with zeroarray
    # GrayA
    v = @inferred(colorview(GrayA, a, b))
    @test eltype(v) == GrayA{Float64}
    @test @inferred(v[2,1]) === GrayA(0.3,0.7)
    v = @inferred(colorview(GrayA{N0f8}, a, b))
    @test @inferred(v[2,1]) === GrayA{N0f8}(0.3,0.7)
    v[1,2] = GrayA(0.25, 0.25)
    @test @inferred(v[1,2]) === GrayA{N0f8}(0.25, 0.25)
    v = @inferred(colorview(GrayA{N0f8}, a, zeroarray))
    @test @inferred(v[2,1]) === GrayA{N0f8}(0.3,0)
    @test_throws ErrorException (v[1,2] = GrayA(0.25, 0.25))
    # RGB
    v = @inferred(colorview(RGB{N0f8}, a, zeroarray, b))
    @test @inferred(v[2,1]) === RGB{N0f8}(0.3,0,0.7)
    v = colorview(RGB{N0f8}, a, z, b)
    @test @inferred(v[2,1]) === RGB{N0f8}(0.3,0,0.7)
    v[2,1] = RGB(0,0.9,0)
    @test @inferred(v[2,1]) === RGB{N0f8}(0,0.9,0)
    @test a[2,1] == b[2,1] == 0
    @test z[2,1] == N0f8(0.9)
    # RGBA
    v = @inferred(colorview(RGBA{N0f8}, a, zeroarray, b, b))
    @test @inferred(v[2,2]) === RGBA{N0f8}(0.4,0,0.8,0.8)
    v = colorview(RGBA{N0f8}, a, copy(z), b, copy(z))
    v[2,2] = RGBA(0.75,0.8,0.75,0.8)
    @test @inferred(v[2,2]) === RGBA{N0f8}(0.75,0.8,0.75,0.8)
    @test a[2,2] == b[2,2] == N0f8(0.75)

    @test eltype(zeroarray) == Union{}
    @test_throws ErrorException colorview(RGB{N0f8}, zeroarray, zeroarray, zeroarray)
    @test_throws DimensionMismatch StackedView(rand(2,3), rand(2,5))

    # With padding
    r = reshape([0.1,0.2], 2, 1)
    g = [false true]
    cv = colorview(RGB, paddedviews(0, r, g, zeroarray)...)
    @test cv[1,1] === RGB(0.1,0,0)
    @test cv[2,1] === RGB(0.2,0,0)
    @test cv[1,2] === RGB(0,1.0,0)
    @test cv[2,2] === RGB(0,0.0,0)
    @test axes(cv) === (Base.OneTo(2), Base.OneTo(2))

    a = [0.1 0.2; 0.3 0.4]
    b = OffsetArray([0.1 0.2; 0.3 0.4], 0:1, 2:3)
    @test_throws DimensionMismatch colorview(RGB, a, b, zeroarray)
    cv = colorview(RGB, paddedviews(0, a, b, zeroarray)...)
    @test axes(cv) == (0:2, 1:3)
    @test red.(cv[axes(a)...]) == a
    @test green.(cv[axes(b)...]) == b
    @test parent(copy(cv)) == [RGB(0,0,0)   RGB(0,0.1,0)   RGB(0,0.2,0);
                               RGB(0.1,0,0) RGB(0.2,0.3,0) RGB(0,0.4,0);
                               RGB(0.3,0,0) RGB(0.4,0,0)   RGB(0,0,0)]
    chanv = channelview(cv)
    @test @inferred(axes(chanv)) == (IdentityUnitRange(1:3), IdentityUnitRange(0:2), IdentityUnitRange(1:3))
    @test chanv[1,1,1] == 0.1
    @test chanv[2,1,2] == 0.3

    @test_throws ErrorException paddedviews(0, zeroarray, zeroarray)

    # Just to boost coverage. These methods are necessary to
    # make inference happy.
    @test_throws ErrorException ImageCore._unsafe_getindex(1, (1,1))
    @test_throws ErrorException ImageCore._unsafe_setindex!(1, (1,1), 0)
end

@testset "Multi-component colorview" begin
    r, g, b = rand(Gray{N0f16},5,5), rand(Gray{N0f16},5,5), rand(Gray{N0f16},5,5)
    A = colorview(RGB, r, g, b)
    @test A[1,1] == RGB(r[1,1], g[1,1], b[1,1])
    o = ones(2, 2)
    A = colorview(RGB, o, o, zeroarray)
    @test A[1,1] == RGB(1,1,0)
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2),Base.OneTo(2))
    A = colorview(RGB, o, zeroarray, o)
    @test A[1,1] == RGB(1,0,1)
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2),Base.OneTo(2))
    A = colorview(RGB, zeroarray, o, o)
    @test A[1,1] == RGB(0,1,1)
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2),Base.OneTo(2))
end

@testset "mosaicviews" begin
    # only test image cases
    @testset "2D inputs" begin
        A1 = fill(Gray(1.), 2, 2)
        A2 = fill(RGB(1., 0., 0.), 3, 3)
        A3 = fill(RGB(0., 1., 0.), 3, 3)
        out = mosaic(A1, A2, A3) |> collect
        @test_reference "references/mosaicviews/2d_opaque_1.png" out by=isequal
        out = mosaic(A1, A2, A3; npad=2, fillvalue=Gray(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/2d_opaque_2.png" out by=isequal
        out = mosaic(A1, A2, A3; npad=2, fillvalue=Gray(0.), nrow=2, rowmajor=true) |> collect
        @test_reference "references/mosaicviews/2d_opaque_3.png" out by=isequal

        A1 = fill(GrayA(1.), 2, 2)
        A2 = fill(RGBA(1., 0., 0.), 3, 3)
        A3 = fill(RGBA(0., 1., 0.), 3, 3)
        out = mosaic(A1, A2, A3) |> collect
        @test_reference "references/mosaicviews/2d_transparent_1.png" out by=isequal
        out = mosaic(A1, A2, A3; npad=2, fillvalue=GrayA(0., 0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/2d_transparent_2.png" out by=isequal
        out = mosaic(A1, A2, A3; npad=2, fillvalue=GrayA(0., 0.), nrow=2, rowmajor=true) |> collect
        @test_reference "references/mosaicviews/2d_transparent_3.png" out by=isequal
    end

    @testset "3D inputs" begin
        A = fill(RGB(0., 0., 0.), 2, 2, 3)
        A[:, :, 1] .= RGB(1., 0., 0.)
        A[:, :, 2] .= RGB(0., 1., 0.)
        A[:, :, 3] .= RGB(0., 0., 1.)
        out = mosaic(A) |> collect
        @test_reference "references/mosaicviews/3d_opaque_1.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=Gray(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/3d_opaque_2.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=Gray(0.), nrow=2, rowmajor=true) |> collect
        @test_reference "references/mosaicviews/3d_opaque_3.png" out by=isequal
        out = mosaic(A, A; npad=2, fillvalue=Gray(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/3d_opaque_4.png" out by=isequal

        A = fill(RGBA(0., 0., 0.), 2, 2, 3)
        A[:, :, 1] .= RGBA(1., 0., 0.)
        A[:, :, 2] .= RGBA(0., 1., 0.)
        A[:, :, 3] .= RGBA(0., 0., 1.)
        out = mosaic(A) |> collect
        @test_reference "references/mosaicviews/3d_transparent_1.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=GrayA(0., 0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/3d_transparent_2.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=GrayA(0.), nrow=2, rowmajor=true) |> collect
        @test_reference "references/mosaicviews/3d_transparent_3.png" out by=isequal
        out = mosaic(A, A; npad=2, fillvalue=GrayA(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/3d_transparent_4.png" out by=isequal
    end

    @testset "4D inputs" begin
        A = fill(RGB(0., 0., 0.), 2, 2, 2, 2)
        A[1, :, 1, 1] .= RGB(1., 0., 0.)
        A[:, :, 1, 2] .= RGB(0., 1., 0.)
        A[:, :, 2, 1] .= RGB(0., 0., 1.)
        out = mosaic(A) |> collect
        @test_reference "references/mosaicviews/4d_opaque_1.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=Gray(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/4d_opaque_2.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=Gray(0.), nrow=2, rowmajor=true) |> collect
        @test_reference "references/mosaicviews/4d_opaque_3.png" out by=isequal
        out = mosaic(A, A; npad=2, fillvalue=Gray(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/4d_opaque_4.png" out by=isequal

        A = fill(RGBA(0., 0., 0.), 2, 2, 2, 2)
        A[1, :, 1, 1] .= RGBA(1., 0., 0.)
        A[:, :, 1, 2] .= RGBA(0., 1., 0.)
        A[:, :, 2, 1] .= RGBA(0., 0., 1.)
        out = mosaic(A) |> collect
        @test_reference "references/mosaicviews/4d_transparent_1.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=GrayA(0., 0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/4d_transparent_2.png" out by=isequal
        out = mosaic(A; npad=2, fillvalue=GrayA(0., 0.), nrow=2, rowmajor=true) |> collect
        @test_reference "references/mosaicviews/4d_transparent_3.png" out by=isequal
        out = mosaic(A, A; npad=2, fillvalue=GrayA(0.), nrow=2) |> collect
        @test_reference "references/mosaicviews/4d_transparent_4.png" out by=isequal
    end
end

@testset "PaddedViews" begin
    # cases that don't promote array eltype:
    #   * (Number, Colorant)
    #   * (Gray, Gray)
    #   * (Gray, Color3)
    #   * (Color3, Color3)
    #   * (Color3, TransparentColor)
    #   * (TransparentColor, TransparentColor)
    @test @inferred(filltype(Float32, Gray{N0f8})) == Gray{N0f8}
    @test @inferred(filltype(Float32, RGB{N0f8})) == RGB{N0f8}
    @test @inferred(filltype(Gray{Float32}, Gray{N0f8})) === Gray{N0f8}
    @test @inferred(filltype(Gray{N0f8}, Gray{Float32})) === Gray{Float32}
    @test @inferred(filltype(Gray{Float32}, RGB{N0f8})) === RGB{N0f8}
    @test @inferred(filltype(Gray{N0f8}, RGB{Float32})) === RGB{Float32}
    @test @inferred(filltype(RGB{Float32}, RGB{N0f8})) === RGB{N0f8}
    @test @inferred(filltype(RGB{N0f8}, RGB{Float32})) === RGB{Float32}
    @test @inferred(filltype(Lab{Float32}, RGB{N0f8})) === RGB{N0f8}
    @test @inferred(filltype(RGB{N0f8}, Lab{Float32})) === Lab{Float32}

    @test @inferred(filltype(Gray{N0f8}, AGray{Float32})) === AGray{Float32}
    @test @inferred(filltype(Gray{Float32}, AGray{N0f8})) === AGray{N0f8}
    @test @inferred(filltype(Gray{N0f8}, ARGB{Float32})) === ARGB{Float32}
    @test @inferred(filltype(Gray{Float32}, ARGB{N0f8})) === ARGB{N0f8}
    @test @inferred(filltype(BGR{N0f8}, ARGB{Float32})) === ARGB{Float32}
    @test @inferred(filltype(BGR{Float32}, ARGB{N0f8})) === ARGB{N0f8}
    @test @inferred(filltype(AGray{N0f8}, ARGB{Float32})) === ARGB{Float32}
    @test @inferred(filltype(AGray{Float32}, ARGB{N0f8})) === ARGB{N0f8}

    # cases that promote both colorant type and storage type
    #   * (Colorant, Number) -- numeric arrays are treated like gray images
    #   * (Color3, Gray)
    #   * (TransparentColor, Colorant)
    @test @inferred(filltype(Gray{N0f8}, Float32)) === Gray{Float32}
    @test @inferred(filltype(Gray{Float64}, Float32)) === Gray{Float32}
    @test @inferred(filltype(RGB{N0f8}, Float32)) === RGB{Float32}
    @test @inferred(filltype(RGB{Float64}, Float32)) === RGB{Float32}
    for T in (Gray24, AGray32, RGB24, ARGB32)
        @test @inferred(filltype(T, Float32)) === @inferred(filltype(Float32, T)) === T
    end

    @test @inferred(filltype(RGB{N0f8}, Gray{Float32})) === RGB{Float32}
    @test @inferred(filltype(RGB{Float32}, Gray{N0f8})) === RGB{Float32}
    @test @inferred(filltype(Lab{Float32}, Gray{N0f8})) === Lab{Float32}
    @test @inferred(filltype(Lab{Float32}, Gray{Float64})) === Lab{Float64}

    @test @inferred(filltype(AGray{N0f8}, Gray{Float32})) === AGray{Float32}
    @test @inferred(filltype(AGray{Float32}, Gray{N0f8})) === AGray{Float32}
    @test @inferred(filltype(AGray{N0f8}, RGB{Float32})) === ARGB{Float32}
    @test @inferred(filltype(AGray{Float32}, RGB{N0f8})) === ARGB{Float32}
    @test @inferred(filltype(AGray{N0f8}, Lab{Float32})) === ALab{Float32}
    @test @inferred(filltype(AGray{Float64}, Lab{Float32})) === ALab{Float64}

    @test @inferred(filltype(ARGB{N0f8}, Gray{Float32})) === ARGB{Float32}
    @test @inferred(filltype(ARGB{Float32}, Gray{N0f8})) === ARGB{Float32}
    @test @inferred(filltype(ARGB{N0f8}, BGR{Float32})) === ABGR{Float32}
    @test @inferred(filltype(ARGB{Float32}, BGR{N0f8})) === ABGR{Float32}
    @test @inferred(filltype(ARGB{N0f8}, Lab{Float32})) === ALab{Float32}
    @test @inferred(filltype(ARGB{Float64}, Lab{Float32})) === ALab{Float64}

    A = Gray{N0f8}[Gray(0.) Gray(0.3); Gray(0.6) Gray(1.)]
    fv = Gray{N0f8}(0)
    Ap = PaddedView(fv, A, (-1:4, -1:4))
    @test eltype(Ap) == Gray{N0f8}
    @test @inferred(getindex(Ap, -1, -1)) === fv
    @test @inferred(getindex(Ap, 2, 2)) === Gray{N0f8}(1.)
    @test Ap[axes(A)...] == A

    fv = RGB{Float32}(1., 0., 0.)
    Ap = PaddedView(fv, A, (-1:4, -1:4))
    @test eltype(Ap) == RGB{Float32}
    @test @inferred(getindex(Ap, -1, -1)) === fv
    @test @inferred(getindex(Ap, 2, 2)) === RGB{Float32}(1., 1., 1.)
    @test Ap[axes(A)...] == RGB{Float32}.(A)

    fv = AGray{Float32}(0.5, 0)
    Ap = PaddedView(fv, A, (-1:4, -1:4))
    @test eltype(Ap) == AGray{Float32}
    @test @inferred(getindex(Ap, -1, -1)) === fv
    @test @inferred(getindex(Ap, 2, 2)) === AGray{Float32}(1., 1.)
    @test Ap[axes(A)...] == AGray{Float32}.(A)

    fv = ARGB{Float32}(0.5, 0., 0., 0.)
    Ap = PaddedView(fv, A, (-1:4, -1:4))
    @test eltype(Ap) == ARGB{Float32}
    @test @inferred(getindex(Ap, -1, -1)) === fv
    @test @inferred(getindex(Ap, 2, 2)) === ARGB{Float32}(1., 1., 1., 1.)
    @test Ap[axes(A)...] == ARGB{Float32}.(A)

    # order irrelevant
    A = mosaic(rand(Float32, 4, 4), rand(RGB{N0f8}, 4, 4))
    @test eltype(A) == RGB{Float32}
    A = mosaic(rand(RGB{N0f8}, 4, 4), rand(Float32, 4, 4))
    @test eltype(A) == RGB{Float32}

    A = mosaic(rand(Float32, 4, 4), rand(Gray{Float64}, 4, 4))
    @test eltype(A) == Gray{Float64}
    A = mosaic(rand(Gray{Float64}, 4, 4), rand(Float32, 4, 4))
    @test eltype(A) == Gray{Float64}

    A = mosaic(rand(Float32, 4, 4), rand(Gray{Float64}, 4, 4), rand(RGB{N0f8}, 4, 4))
    @test eltype(A) == RGB{Float64}
    A = mosaic(rand(Float32, 4, 4), rand(RGB{N0f8}, 4, 4), rand(Gray{Float64}, 4, 4))
    @test eltype(A) == RGB{Float64}

    A = mosaic(rand(Gray, 4, 4), rand(RGB, 4, 4))
    @test eltype(A) == RGB{Float64} # the filltype is always a concrete type
end

nothing


text:
using Neurodevo
using Test

include("test/controller.jl")
include("test/step.jl")
include("test/model.jl")
include("test/data.jl")


text:
# Compare all four implementations of demand-oblivious routing; they should
# give (approximately) the same result.

include("topologies.jl")

using PrettyTables, LightGraphs, DelimitedFiles
using JuMP, CPLEX, Gurobi, Mosek, MosekTools

opt = CPLEX.Optimizer
# opt = Gurobi.Optimizer
# opt = optimizer_with_attributes(Gurobi.Optimizer, "NonConvex" => 2)
# opt = Mosek.Optimizer

algos = [(CuttingPlane(), FlowFormulation(), false),
         (CuttingPlane(), PathFormulation(), false),
         (CuttingPlane(), PathFormulation(), true),
         (DualReformulation(), FlowFormulation(), false),
         (DualReformulation(), PathFormulation(), false),
         (DualReformulation(), PathFormulation(), true)]
headers = ["Cutting plane, flow formulation", "Cutting plane, path formulation", "Cutting plane, column generation",
           "Dual reformulation, flow formulation", "Dual reformulation, path formulation", "Dual reformulation, column generation"]

ftopologies = [t for t in topologies if ! (t.name in ["att", "abilene"])]
table_ratio = Matrix{Any}(undef, length(ftopologies), length(algos))
table_time = Matrix{Any}(undef, length(ftopologies), length(algos))

for (i, t) in enumerate(ftopologies)
    println("=> $(t.name)")
    g, k, dm = topology_to_graphs(t)

    for (j, (impl_algo, impl_form, impl_cg)) in enumerate(algos)
        println("  -> $(impl_algo), $(impl_form), column generation: $(impl_cg)")

        mt = ModelType(Load(), MinimumMaximum(), impl_form, impl_cg, impl_algo, ObliviousUncertainty(), UncertainDemand())
        rd = RoutingData(g, k, opt, mt, name=t.name, logmessage=(x)->nothing)
        sol = compute_routing(rd)

        table_ratio[i, j] = sol.objectives[sol.n_iter]
        table_time[i, j] = sol.total_time_ms
    end
end

@ptconfclean
@ptconf tf = unicode_rounded
println("\n\n  Ratios\n==========")
@pt :header = headers table_ratio
println("\n\n  Total times [ms]\n====================")
@pt :header = headers table_time


text:
using CauseMap

# specify parameter ranges to test over
E_vals     = 2:10  # range to test of system dimensionality
tau_s_vals = 1:1 # range for lag length for manifold reconstruction
tau_p_vals = 0:15  # range to test for time lag of causal effect
#

# read in data
para, didi = load_example_data("ParaDidi")

# run analysis
makeoptimizationplots(para, didi,
                                        E_vals, tau_s_vals, tau_p_vals,
                                        "Para.", "Didi.";
                                        nreps=10, left_E=2, left_tau_p=0,   # optional
                                        right_E=7, right_tau_p=12, lagunit=.5,
                                        unit="days", show_tau_s=false # optional
                                    )


text:
using Base.Meta: parse

#################
# Overload of ~ #
#################

"""
    struct Model{pvars, dvars, F, TD}
        f::F
        data::TD
    end

A `Model` struct with parameter variables `pvars`, data variables `dvars`, inner
function `f` and `data::NamedTuple`.
"""
struct Model{pvars, dvars, F, TD}
    f::F
    data::TD
end
function Model{pvars, dvars}(f::F, data::TD) where {pvars, dvars, F, TD}
    return Model{pvars, dvars, F, TD}(f, data)
end
pvars(m::Model{params}) where {params} = Tuple(params.types)
dvars(m::Model{params, data}) where {params, data} = Tuple(data.types)
@generated function inpvars(::Val{sym}, ::Model{params}) where {sym, params}
    return sym in params.types ? :(true) : :(false)
end
@generated function indvars(::Val{sym}, ::Model{params, data}) where {sym, params, data}
    return sym in data.types ? :(true) : :(false)
end

(model::Model)(args...; kwargs...) = model.f(args..., model; kwargs...)

# TODO: Replace this macro, see issue #514
"""
Usage: @VarName x[1,2][1+5][45][3]
  return: (:x,[1,2],6,45,3)
"""
macro VarName(expr::Union{Expr, Symbol})
    ex = deepcopy(expr)
    isa(ex, Symbol) && return var_tuple(ex)
    (ex.head == :ref) || throw("VarName: Mis-formed variable name $(expr)!")
    inds = :(())
    while ex.head == :ref
        if length(ex.args) >= 2
            pushfirst!(inds.args, Expr(:vect, ex.args[2:end]...))
            end
        ex = ex.args[1]
        isa(ex, Symbol) && return var_tuple(ex, inds)
    end
    throw("VarName: Mis-formed variable name $(expr)!")
end
function var_tuple(sym::Symbol, inds::Expr=:(()))
    return esc(:($(QuoteNode(sym)), $inds, $(QuoteNode(gensym()))))
end

function wrong_dist_errormsg(l)
    return "Right-hand side of a ~ must be subtype of Distribution or a vector of" *
        "Distributions on line $(l)."
end

"""
    generate_observe(observation, dist, model_info)

Generate an observe expression for observation `observation` drawn from
a distribution or a vector of distributions (`dist`).
"""
function generate_observe(observation, dist, model_info)
    main_body_names = model_info[:main_body_names]
    vi = main_body_names[:vi]
    sampler = main_body_names[:sampler]
    return quote
        isdist = if isa($dist, AbstractVector)
            # Check if the right-hand side is a vector of distributions.
            all(d -> isa(d, Distribution), $dist)
        else
            # Check if the right-hand side is a distribution.
            isa($dist, Distribution)
        end
        @assert isdist @error($(wrong_dist_errormsg(@__LINE__)))
        $vi.logp += Turing.observe($sampler, $dist, $observation, $vi)
    end
end

"""
    generate_assume(var, dist, model_info)

Generate an assume expression for parameters `var` drawn from
a distribution or a vector of distributions (`dist`).
"""
function generate_assume(var::Union{Symbol, Expr}, dist, model_info)
    main_body_names = model_info[:main_body_names]
    vi = main_body_names[:vi]
    sampler = main_body_names[:sampler]

    varname = gensym(:varname)
    sym, idcs, csym = gensym(:sym), gensym(:idcs), gensym(:csym)
    csym_str, indexing, syms = gensym(:csym_str), gensym(:indexing), gensym(:syms)

    if var isa Symbol
        varname_expr = quote
            $sym, $idcs, $csym = @VarName $var
            $csym = Symbol($(model_info[:name]), $csym)
            $syms = Symbol[$csym, $(QuoteNode(var))]
            $varname = Turing.VarName($vi, $syms, "")
        end
    else
        varname_expr = quote
            $sym, $idcs, $csym = @VarName $var
            $csym_str = string($(model_info[:name]))*string($csym)
            $indexing = mapfoldl(string, *, $idcs, init = "")
            $varname = Turing.VarName($vi, Symbol($csym_str), $sym, $indexing)
        end
    end

    lp = gensym(:lp)
    return quote
        $varname_expr
        isdist = if isa($dist, AbstractVector)
            # Check if the right-hand side is a vector of distributions.
            all(d -> isa(d, Distribution), $dist)
        else
            # Check if the right-hand side is a distribution.
            isa($dist, Distribution)
        end
        @assert isdist @error($(wrong_dist_errormsg(@__LINE__)))

        ($var, $lp) = if isa($dist, AbstractVector)
            Turing.assume($sampler, $dist, $varname, $var, $vi)
        else
            Turing.assume($sampler, $dist, $varname, $vi)
        end
        $vi.logp += $lp
    end
end

"""
    tilde(left, right, model_info)

The `tilde` function generates observation expression for data variables and assumption expressions for parameter variables, updating `model_info` in the process.
"""
function tilde(left, right, model_info)
    return generate_observe(left, right, model_info)
end
function tilde(left::Union{Symbol, Expr}, right, model_info)
    return _tilde(getvsym(left), left, right, model_info)
end

function _tilde(vsym, left, dist, model_info)
    main_body_names = model_info[:main_body_names]
    model_name = main_body_names[:model]

    if vsym in model_info[:arg_syms]
        if !(vsym in model_info[:tent_dvars_list])
            @debug " Observe - `$(vsym)` is an observation"
            push!(model_info[:tent_dvars_list], vsym)
        end

        return quote
            if Turing.indvars($(Val(vsym)), $model_name)
                $(generate_observe(left, dist, model_info))
            else
                $(generate_assume(left, dist, model_info))
            end
        end
    else
        # Assume it is a parameter.
        if !(vsym in model_info[:tent_pvars_list])
            msg = " Assume - `$(vsym)` is a parameter"
            if isdefined(Main, vsym)
                msg  *= " (ignoring `$(vsym)` found in global scope)"
            end

            @debug msg
            push!(model_info[:tent_pvars_list], vsym)
        end

        return generate_assume(left, dist, model_info)
    end
end

#################
# Main Compiler #
#################

"""
    @model(body)

Macro to specify a probabilistic model.

Example:

Model definition:

```julia
@model model_generator(x = default_x, y) = begin
    ...
end
```

Expanded model definition

```julia
# Allows passing arguments as kwargs
model_generator(; x = nothing, y = nothing)) = model_generator(x, y)
function model_generator(x = nothing, y = nothing)
    pvars, dvars = Turing.get_vars(Tuple{:x, :y}, (x = x, y = y))
    data = Turing.get_data(dvars, (x = x, y = y))

    inner_function(sampler::Turing.AnySampler, model) = inner_function(model)
    function inner_function(model)
        return inner_function(Turing.VarInfo(), Turing.SampleFromPrior(), model)
    end
    function inner_function(vi::Turing.VarInfo, model)
        return inner_function(vi, Turing.SampleFromPrior(), model)
    end
    # Define the main inner function
    function inner_function(vi::Turing.VarInfo, sampler::Turing.AnySampler, model)
        local x
        if isdefined(model.data, :x)
            x = model.data.x
        else
            x = default_x
        end
        local y
        if isdefined(model.data, :y)
            y = model.data.y
        else
            y = nothing
        end

        vi.logp = zero(Real)
        ...
    end
    model = Turing.Model{pvars, dvars}(inner_function, data)
    return model
end
```

Generating a model: `model_generator(x_value)::Model`.
"""
macro model(input_expr)
    build_model_info(input_expr) |> translate_tilde! |> update_args! |> build_output
end

"""
    build_model_info(input_expr)

Builds the `model_info` dictionary from the model's expression.
"""
function build_model_info(input_expr)
    # Extract model name (:name), arguments (:args), (:kwargs) and definition (:body)
    modeldef = MacroTools.splitdef(input_expr)
    # Function body of the model is empty
    warn_empty(modeldef[:body])
    # Construct model_info dictionary

    arg_syms = [(arg isa Symbol) ? arg : arg.args[1] for arg in modeldef[:args]]
    model_info = Dict(
        :name => modeldef[:name],
        :input_expr => input_expr,
        :main_body => modeldef[:body],
        :arg_syms => arg_syms,
        :args => modeldef[:args],
        :kwargs => modeldef[:kwargs],
        :tent_dvars_list => Symbol[],
        :tent_pvars_list => Symbol[],
        :main_body_names => Dict(
            :vi => gensym(:vi),
            :sampler => gensym(:sampler),
            :model => gensym(:model),
            :pvars => gensym(:pvars),
            :dvars => gensym(:dvars),
            :data => gensym(:data),
            :inner_function => gensym(:inner_function)
        )
    )

    return model_info
end

"""
    translate_tilde!(model_info)

Translates ~ expressions to observation or assumption expressions, updating `model_info`.
"""
function translate_tilde!(model_info)
    ex = model_info[:main_body]
    ex = MacroTools.postwalk(x -> @capture(x, L_ ~ R_) ? tilde(L, R, model_info) : x, ex)
    model_info[:main_body] = ex
    return model_info
end

"""
    update_args!(model_info)

Extracts default argument values and replaces them with `nothing`.
"""
function update_args!(model_info)
    fargs = model_info[:args]
    fargs_default_values = Dict()
    for i in 1:length(fargs)
        if isa(fargs[i], Symbol)
            fargs_default_values[fargs[i]] = :nothing
            fargs[i] = Expr(:kw, fargs[i], :nothing)
        elseif isa(fargs[i], Expr) && fargs[i].head == :kw
            fargs_default_values[fargs[i].args[1]] = fargs[i].args[2]
            fargs[i] = Expr(:kw, fargs[i].args[1], :nothing)
        else
            throw("Unsupported argument type $(fargs[i]).")
        end
    end
    model_info[:args] = fargs
    model_info[:arg_defaults] = fargs_default_values

    return model_info
end

"""
    build_output(model_info)

Builds the output expression.
"""
function build_output(model_info)
    # Construct user-facing function
    main_body_names = model_info[:main_body_names]
    vi_name = main_body_names[:vi]
    model_name = main_body_names[:model]
    sampler_name = main_body_names[:sampler]
    data_name = main_body_names[:data]
    pvars_name = main_body_names[:pvars]
    dvars_name = main_body_names[:dvars]
    inner_function_name = main_body_names[:inner_function]

    args = model_info[:args]
    arg_syms = model_info[:arg_syms]
    outer_function_name = model_info[:name]
    tent_pvars_list = model_info[:tent_pvars_list]
    tent_dvars_list = model_info[:tent_dvars_list]
    main_body = model_info[:main_body]
    arg_defaults = model_info[:arg_defaults]

    if length(tent_dvars_list) == 0
        tent_dvars_nt = :(NamedTuple())
    else
        tent_dvars_nt = :($([:($var = $var) for var in tent_dvars_list]...),)
    end

    #= Does the following for each of the tentative dvars
        local x
        if isdefined(model.data, :x)
            x = model.data.x
        else
            x = default_x
        end
    =#
    unwrap_data_expr = Expr(:block)
    for var in tent_dvars_list
        push!(unwrap_data_expr.args, quote
            local $var
            if isdefined($model_name.data, $(QuoteNode(var)))
                $var = $model_name.data.$var
            else
                $var = $(arg_defaults[var])
            end
        end)
    end

    return esc(quote
        # Allows passing arguments as kwargs
        $outer_function_name(;$(args...)) = $outer_function_name($(arg_syms...))
        # Outer function with `nothing` as default values
        function $outer_function_name($(args...))
            # Adds variables equal to `nothing` to pvars and the rest to dvars
            # `tent_pvars_list` is the tentative list of pvars
            # `tent_dvars_nt` is the tentative named tuple of dvars
            $pvars_name, $dvars_name = Turing.get_vars($(Tuple{tent_pvars_list...}), $(tent_dvars_nt))
            # Filter out the dvars equal to `nothing`
            $data_name = Turing.get_data($dvars_name, $tent_dvars_nt)

            # Define fallback inner functions
            function $inner_function_name($sampler_name::Turing.AnySampler, $model_name)
                return $inner_function_name($model_name)
            end
            function $inner_function_name($model_name)
                return $inner_function_name(Turing.VarInfo(), Turing.SampleFromPrior(), $model_name)
            end
            function $inner_function_name($vi_name::Turing.VarInfo, $model_name)
                return $inner_function_name($vi_name, Turing.SampleFromPrior(), $model_name)
            end

            # Define the main inner function
            function $inner_function_name(
                $vi_name::Turing.VarInfo,
                $sampler_name::Turing.AnySampler,
                $model_name
                )

                $unwrap_data_expr
                $vi_name.logp = zero(Real)
                $main_body
            end
            $model_name = Turing.Model{$pvars_name, $dvars_name}($inner_function_name, $data_name)
            return $model_name
        end
    end)
end

@generated function get_vars(tent_pvars::Type{Tpvars}, tent_dvars_nt::NamedTuple) where {Tpvars <: Tuple}
    tent_pvar_syms = [Tpvars.types...]
    tent_dvar_syms = [tent_dvars_nt.names...]
    dvar_types = [tent_dvars_nt.types...]
    append!(tent_pvar_syms, [tent_dvar_syms[i] for i in 1:length(tent_dvar_syms) if dvar_types[i] == Nothing])
    setdiff!(tent_dvar_syms, tent_pvar_syms)
    pvars_tuple = Tuple{tent_pvar_syms...}
    dvars_tuple = Tuple{tent_dvar_syms...}

    return :($pvars_tuple, $dvars_tuple)
end

@generated function get_data(::Type{Tdvars}, nt) where Tdvars
    dvars = Tdvars.types
    args = []
    for var in dvars
        push!(args, :($var = nt.$var))
    end
    if length(args) == 0
        return :(NamedTuple())
    else
        return :($(args...),)
    end
end

function warn_empty(body)
    if all(l -> isa(l, LineNumberNode), body.args)
        @warn("Model definition seems empty, still continue.")
    end
    return
end

####################
# Helper functions #
####################

getvsym(s::Symbol) = s
function getvsym(expr::Expr)
    @assert expr.head == :ref "expr needs to be an indexing expression, e.g. :(x[1])"
    return getvsym(expr.args[1])
end


text:
@testset "Lennard Jones 150 Cluster Test" begin
    ϵ = 1.0u"eV"
    σ = 1.0u"Å"
    rcutoff = 100.0u"Å"
    species = [:Ar, :H]
    p = LennardJones(ϵ, σ, rcutoff, species)

    box = [[8.0, 0.0, 0.0], [0.0, 8.0, 0.0], [0.0, 0.0, 8.0]] * 1u"Å"
    bcs = [DirichletZero(), DirichletZero(), DirichletZero()]
    l = readlines("integration/lj_clusters/150.xyz")
    r = [ parse.(Float64, split(li)) for li in l ]
    atoms = [Atom(:Ar, ri * u"Å") for ri in r]
    system = FlexibleSystem(atoms, box, bcs)

    e_a_tol = auconvert(1e-5 * 1u"eV")
    es, fs = energy_and_force(system, p)
    @test es ≈ auconvert(-893.310258 * 1u"eV") atol = e_a_tol
    fi = auconvert(0.0 * 1u"eV/Å")
    f_a_tol = auconvert(1e-10 * 1u"eV/Å")
    @test sum(fs)[1] ≈ fi  atol=f_a_tol
    @test sum(fs)[2] ≈ fi  atol=f_a_tol
    @test sum(fs)[3] ≈ fi  atol=f_a_tol

    for ff in fs
        for ffi in ff
            @test ffi ≈ fi atol = f_a_tol
        end
    end
    @test virial(system, p) ≈ auconvert(0.0 * 1u"eV") atol=e_a_tol

end


text:
"""
    poormans_tokenize

Tokenizes by removing punctuation and splitting on spaces
"""
function poormans_tokenize(source::AbstractString)
    cleaned = filter(s->!ispunct(s), source)
    split(cleaned)
end

"""
    punctuation_space_tokenize

Tokenizes by removing punctuation, unless it occurs inside of a word.
"""
function punctuation_space_tokenize(source::AbstractString)
    preprocced = source
    pass1=replace(preprocced,r"[[:punct:]]*[[:space:]][[:punct:]]*" => " ")
    pass2=replace(pass1,r"[[:punct:]]*$|^[[:punct:]]*" => "")
    split(pass2)
end


text:
@testset "bias correction for time series" begin
    for i in 1:10
        rad = Array{Union{Float64, Missing}}(rand(1:100.0, 100))
        for j in 1:10
            rad[rand(1:100)] = missing
        end
        clouds = rand(1:31, 100)
        @test length(bias_correction(rad, clouds)) == length(rad)
    end
end

@testset "bias correction for datacube" begin
    for i in 1:10
        x = rand(10:30)
        y = rand(10:30)
        z = rand(10:30)
        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x, y, z))
        for j in 1:rand(1:20)
            rad[rand(1:x), rand(1:y), rand(1:z)] = missing
        end
        clouds = rand(1:30, x, y, z)
        mask = rand(0:1, x, y)
        @test size(bias_correction(rad, clouds, mask)) == (x, y, z)
    end
end


text:
include("..\\libs\\FunctionsRobot.jl")

function perimetrAroundObject!(r::Robot, side)
    check = false

    while !check
        if ismarker(r)
            break
        end
        while isborder(r, nextSideConterclockwise(side))
            putmarker!(r)
            move!(r, side)
        end
        putmarker!(r)
        side = nextSideConterclockwise(side)
        move!(r, side)
    end
end

function perimetr_of_object!(r::Robot)
    sides = moveAndReturnDirections!(r)
    perimetrAroundObject!(r, searchObject!(r))
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end


text:
include("hodgkinhuxley.jl")
include("prinzneuron.jl")
include("simplesynapse.jl")


text:
# This test is derived from test cmcompute.c in Calceph version 2.3.2
# the test data files are copied from calceph-2.3.2.tar.gz
function testFunction1(testFile,ephFiles,pflag)

    eph = Ephem(ephFiles)
    if pflag
        prefetch(eph)
    end

    con = constants(eph)
    AU = con[:AU]

    f = open(testFile);
    for ln in eachline(f)
        elts=split(ln)
        jd0=parse(Float64,elts[1])
        target=parse(Int,elts[2])
        center=parse(Int,elts[3])
        dt = jd0 - trunc(Int,jd0);
        jd0 = trunc(Int,jd0) + 2.4515450000000000000E+06
        ref = [parse(Float64, x) for x in elts[4:end]]
        val = compute(eph,jd0,dt,target,center)

        ϵ = 1.0e-8
        val0 = val[:]
        if (target==15)
            ϵ = 1.0e-7
            while val[3]>2π
                val[3]-=2π
            end
            while val[3]<=0
                val[3]+=2π
            end
        end
        for i in 1:6
            @test abs(ref[i]-val[i]) < ϵ
        end
        ref = val0
        ϵ = 3.0e-15
        if target ∉ [15,16,17]
            val = compute(eph,jd0,dt,target,center,unitAU+unitDay)
            for i in 1:6
                @test abs(ref[i]-val[i]) < ϵ
            end
            val = compute(eph,jd0,dt,target,center,unitAU+unitSec)
            for i in 1:6
                if i>3
                    val[i]*=86400
                end
                @test abs(ref[i]-val[i]) < ϵ
            end
            ϵ = 3.0e-14;
            val = compute(eph,jd0,dt,target,center,unitKM+unitDay)
            for i in 1:6
                @test abs(ref[i]-val[i]/AU) < ϵ
            end

            val = compute(eph,jd0,dt,target,center,unitKM+unitSec)
            for i in 1:6
                if i>3
                    val[i]*=86400
                end
                @test abs(ref[i]-val[i]/AU) < ϵ
            end

            ϵ = 3.0e-15
            val = compute(eph,jd0,dt,target,center,unitDay+unitAU,3)
            @test length(val)==12
            for i in 1:6
                @test abs(ref[i]-val[i]) < ϵ
            end

            ref = val
            val = compute(eph,jd0,dt,target,center,unitDay+unitAU,2)
            @test length(val)==9
            for i in 1:9
                @test abs(ref[i]-val[i]) < ϵ
            end
            val = compute(eph,jd0,dt,target,center,unitDay+unitAU,1)
            @test length(val)==6
            for i in 1:6
                @test abs(ref[i]-val[i]) < ϵ
            end
            val = compute(eph,jd0,dt,target,center,unitDay+unitAU,0)
            @test length(val)==3
            for i in 1:3
                @test abs(ref[i]-val[i]) < ϵ
            end

            ϵ = 3.0e-14
            val = compute(eph,jd0,dt,target,center,unitSec+unitKM,3)
            @test length(val)==12
            for i in 1:12
                if i>3
                    val[i]*=86400
                end
                if i>6
                    val[i]*=86400
                end
                if i>9
                    val[i]*=86400
                end
                @test abs(ref[i]-val[i]/AU) < ϵ
            end

        elseif target == 15
            val = compute(eph,jd0,dt,target,center,unitRad+unitDay)
            for i in 1:6
                @test abs(ref[i]-val[i]) < ϵ
            end
            val = compute(eph,jd0,dt,target,center,unitRad+unitSec)
            for i in 1:6
                if i>3
                    val[i]*=86400
                end
                @test abs(ref[i]-val[i]) < ϵ
            end
        elseif target ∈ [16,17]
            ϵ = 1e-18
            val = compute(eph,jd0,dt,target,center,unitSec)
            @test abs(ref[1]-val[1]) < ϵ
            val = compute(eph,jd0,dt,target,center,unitDay)
            @test abs(ref[1]-val[1]*86400) < ϵ*86400
        end

    end
    close(f)

end


text:
module CUDA

using ..CuArrays

if CuArrays.libcudnn !== nothing  # TODO: use CuArrays.has_cudnn()
  using CuArrays: CUDNN
  include("curnn.jl")
  include("cudnn.jl")
else
  @warn "CUDNN is not installed, some functionality will not be available."
end

end


text:
function calDiffAdv_QUICKEST!(
    Xflx_U  :: AbstractArray{Float64},
    Xflx_V  :: AbstractArray{Float64},
    Xflx_W  :: AbstractArray{Float64},

    X       :: AbstractArray{Float64},
    u_U     :: AbstractArray{Float64},
    v_V     :: AbstractArray{Float64},
    w_W     :: AbstractArray{Float64},
    amo     :: AdvancedMatrixOperators,
    Kh_U    :: AbstractArray{Float64},
    Kh_V    :: AbstractArray{Float64},
    Kv      :: Float64,
    Δt      :: Float64,
    wksp    :: Workspace,
)

    # x direction
    calFluxDensity_abstract!(;
        Xflx_bnd            = Xflx_U,

        X_T                 = X,
        u_bnd               = u_U,
        K_bnd               = Kh_U,
        Δt                  = Δt,
        Δx_bnd              = amo.gd.Δx_U,

        bnd_∂x_T         = amo.U_∂x_T,
        T_∂x_bnd         = amo.T_∂x_U,
        bnd_pos_dir_T    = amo.bmo.U_E_T,
        bnd_neg_dir_T    = amo.bmo.U_W_T,
        bnd_interp_T     = amo.U_interp_T,
        bnd_mask_bnd     = amo.U_mask_U,

        grid_bnd = :U,
        wksp = wksp,
    )

    calFluxDensity_abstract!(;
        Xflx_bnd            = Xflx_V,

        X_T                 = X,
        u_bnd               = v_V,
        K_bnd               = Kh_V,
        Δt                  = Δt,
        Δx_bnd              = amo.gd.Δy_V,

        bnd_∂x_T         = amo.V_∂y_T,
        T_∂x_bnd         = amo.T_∂y_V,
        bnd_pos_dir_T    = amo.bmo.V_N_T,
        bnd_neg_dir_T    = amo.bmo.V_S_T,
        bnd_interp_T     = amo.V_interp_T,
        bnd_mask_bnd     = amo.V_mask_V,

        grid_bnd = :V,
        wksp = wksp,
    )

    calFluxDensity_abstract!(;
        Xflx_bnd            = Xflx_W,

        X_T                 = X,
        u_bnd               = w_W,
        K_bnd               = Kv,
        Δt                  = Δt,
        Δx_bnd              = amo.gd.Δz_W,

        bnd_∂x_T         = amo.W_∂z_T,
        T_∂x_bnd         = amo.T_∂z_W,
        bnd_pos_dir_T    = amo.bmo.W_UP_T,
        bnd_neg_dir_T    = amo.bmo.W_DN_T,
        bnd_interp_T     = amo.W_interp_T,
        bnd_mask_bnd   = amo.W_mask_W,

        grid_bnd = :W,
        wksp = wksp,
    )

end

function calFluxDensity_abstract!(;

    Xflx_bnd              :: AbstractArray{Float64},

    X_T                       :: AbstractArray{Float64},
    u_bnd                     :: AbstractArray{Float64},
    K_bnd                     :: Union{Float64, Array{Float64}},
    Δt                        :: Float64,
    Δx_bnd                    :: AbstractArray{Float64},

    bnd_∂x_T                  :: AbstractArray{Float64, 2},
    T_∂x_bnd                  :: AbstractArray{Float64, 2},
    bnd_pos_dir_T             :: AbstractArray{Float64, 2},
    bnd_neg_dir_T             :: AbstractArray{Float64, 2},
    bnd_interp_T              :: AbstractArray{Float64, 2},
    bnd_mask_bnd              :: AbstractArray{Float64, 2},

    grid_bnd :: Symbol,
    wksp     :: Workspace,
)

    # I keep the original memory intensive code for reference.
    # The original code does not utilize Workspace. When running
    # with 4 cores, 35 vertical layers, gx1v6 grid and 8 substeps,
    # it takes about 12.5 secs while new code using Workspace
    # takes about 7 secs. This is a huge improvement.

    #=
    _X_T    = view(X_T, :)
    _u_bnd  = view(u_bnd, :)
    _Δx_bnd = view(Δx_bnd, :)

    GRAD_bnd = bnd_∂x_T * _X_T
    CURV_T   = T_∂x_bnd * GRAD_bnd

    CURV_on_pos_dir_bnd = bnd_neg_dir_T * CURV_T
    CURV_on_neg_dir_bnd = bnd_pos_dir_T * CURV_T

    pos_u_mask_bnd = _u_bnd .>= 0

    CURV_r_bnd = CURV_on_neg_dir_bnd .* pos_u_mask_bnd + CURV_on_pos_dir_bnd .* (1 .- pos_u_mask_bnd)

    uΔt_bnd    = u_bnd * Δt
    X_star_bnd = bnd_interp_T * _X_T

    X_star_bnd = bnd_interp_T * _X_T - uΔt_bnd / 2.0 .* GRAD_bnd + ( K * Δt / 2.0 .- (_Δx_bnd.^2.0)/6.0 + (uΔt_bnd.^2.0) / 6.0 ) .* CURV_r_bnd

    Xflx_bnd[:] = bnd_mask_bnd * ( u_bnd .* X_star_bnd - K * ( GRAD_bnd - uΔt_bnd / 2.0 .* CURV_r_bnd ) )

    =#
    _X_T    = view(X_T, :)
    _u_bnd  = view(u_bnd, :)
    _Δx_bnd = view(Δx_bnd, :)
    _Xflx_bnd = view(Xflx_bnd, :)

    GRAD_bnd            = getSpace!(wksp, grid_bnd, true)
    CURV_T              = getSpace!(wksp, :T, true)
    CURV_on_pos_dir_bnd = getSpace!(wksp, grid_bnd, true)
    CURV_on_neg_dir_bnd = getSpace!(wksp, grid_bnd, true)
    pos_u_mask_bnd      = getSpace!(wksp, grid_bnd, true)
    CURV_r_bnd          = getSpace!(wksp, grid_bnd, true)
    uΔt_bnd             = getSpace!(wksp, grid_bnd, true)
    X_star_bnd          = getSpace!(wksp, grid_bnd, true)
    tmp_bnd             = getSpace!(wksp, grid_bnd, true)

    mul!(GRAD_bnd, bnd_∂x_T, _X_T)
    mul!(CURV_T,   T_∂x_bnd, GRAD_bnd)

    mul!(CURV_on_pos_dir_bnd, bnd_neg_dir_T, CURV_T)
    mul!(CURV_on_neg_dir_bnd, bnd_pos_dir_T, CURV_T)

    @. pos_u_mask_bnd = _u_bnd >= 0

    @. CURV_r_bnd = CURV_on_neg_dir_bnd * pos_u_mask_bnd + CURV_on_pos_dir_bnd * (1 - pos_u_mask_bnd)

    @. uΔt_bnd = u_bnd * Δt

    mul!(X_star_bnd, bnd_interp_T, _X_T)
    @. X_star_bnd += - uΔt_bnd / 2.0 * GRAD_bnd + ( K_bnd * Δt / 2.0 - (_Δx_bnd^2.0)/6.0 + (uΔt_bnd^2.0) / 6.0 ) * CURV_r_bnd

    @. tmp_bnd = u_bnd .* X_star_bnd - K_bnd * ( GRAD_bnd - uΔt_bnd / 2.0 .* CURV_r_bnd )
    mul!(_Xflx_bnd, bnd_mask_bnd, tmp_bnd)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module LickometerAnalysis
using DelimitedFiles
using PyPlot
using Peaks
using Images
using Statistics

export x_ms2min, x_ms2sec, setthresh, roll1d, baselinecorrection, findlicks, filterlicks, findlickratePerBout, min2time, time2min, find_bout_boundary, lickcount

# Write your package code here.
"""convert milliseconds to miniutes"""
function x_ms2min(cycles, waittime)
  xmin = 0:(waittime/1000/60):(cycles-1)*(waittime/1000/60) #Convert to minutes
  return(xmin)
end

"""convert milliseconds to seconds"""
function x_ms2sec(cycles, waittime)
  xsec = 0:(waittime/1000):(cycles-1)*(waittime/1000) #Convert to minutes
  return(xsec)
end

"""threshold is the `frac` of the `maxidx`-th highest value in the `dat` vector"""
function setthresh(dat, maxidx, frac)
  thresh = sort(dat, rev = true)[maxidx]*frac
  return(thresh)
end

function roll1d(f, v, win)
  output = zeros(length(v)-win)
  for (i, val) in enumerate(Int(1-win/2):Int(length(v)-win-win/2))
    output[i] = f(v[val:val+win-1])
  end
  output
end

function baselinecorrection(dat, win, filt)
  tmp = padarray(dat, Pad(:symmetric, win))
  val  = roll1d(filt, tmp, win*2)
  datn = dat .- val
  return(datn)
end

"""Finding local maxima that are higher than the threshold"""
function findlicks(v, thresh, win)
  vout = Vector{Int}()
  v_jitter = v .+ rand(length(v)).*1e-10
  localmax_idx, localmax_val = findmaxima(v_jitter, win) ## Critial factor = window size
  for (i, val) in enumerate(localmax_val)
    if val > thresh
      push!(vout, localmax_idx[i])
    end
  end
  vout
end

"""Finding local maxima that are higher than the threshold with `filterlicks`"""
function findlicks(vec, thresh, win, filterwinsize)
  vout = findlicks(vec, thresh, win)
  lickidx = filterlicks(vec, thresh, vout, filterwinsize)
  return(lickidx)
end

""" Find local maxima whose neighbors are lower than the threshold"""
function filterlicks(vec::AbstractVector, thresh::Number, idx::Number, winsize::Number)
  vec_thresh = vec .> thresh
  return(!all(vec_thresh[idx-winsize:idx+1]))
end

function filterlicks(vec::AbstractVector, thresh::Number, lickidx::AbstractVector, winsize::Number)
  keep = falses(length(lickidx))
  for (n, i) in enumerate(lickidx)
    keep[n] = filterlicks(vec, thresh, i, winsize)
  end
  return(lickidx[keep])
end

function findlickratePerBout(lickidx, waittime, boutinterval_thresh)
  bout_boundary = findall(diff(lickidx) .> boutinterval_thresh)
  dts = [(lickidx[bout_boundary[i+1]] - lickidx[bout_boundary[i]+1])*waittime for i in 1:length(bout_boundary)-1]
  nlicksperbout = [bout_boundary[i+1]-bout_boundary[i]+1+1 for i in 1:length(bout_boundary)-1]
  lickratePerBout = nlicksperbout./dts*1000
  lickratePerBout
end

""" hour:min:sec to minutes. `offset` in second"""
function min2time(minute; offset = 0, scale = 1)
  sec = minute*60
  sec = sec-offset
  minute = sec/60
  minute = minute * scale
  hour = minute/60
  h = floor(Int, hour)
  m = hour%1*60
  s = m%1*60
  m = floor(Int, m)
  return(h, m, s)
end

""" minute to hour:min:sec `offset` in second"""
function time2min(h, m, s; offset = 0)
  s = s-offset
  m = m + s/60
  return(h*60+m)
end

"""find lick boundary: A bout is consist of multiple licks"""
function find_bout_boundary(lickidx, boutinterval_thresh)
  d_lickidx = diff(lickidx)
  bout_boundary = findall(d_lickidx .> boutinterval_thresh)
  v = Vector{UnitRange}(undef, length(bout_boundary));
  for i in eachindex(bout_boundary)
    if i == 1
      v[i] = lickidx[i]:lickidx[bout_boundary[i]]
    else
      v[i] = d_lickidx[bout_boundary[i-1]] + lickidx[bout_boundary[i-1]]: lickidx[bout_boundary[i]]
    end
  end
  v
end

"""
Count licks within a give time window.
Time window can be set by `timefirst` and `timelast` in ms, sec, or min, or else depending on the unit of `x_axis`.
"""
function lickcount(x_axis, lickidx, timefirst, timelast)
  keep = x_axis .>= timefirst .&& x_axis .<= timelast #Time window from 12am - 8am.
  firstidx = findfirst(keep .== 1)
  lastidx = findlast(keep .== 1)
  valout = sum(lickidx .>= firstidx .&& lickidx .<= lastidx)
  return(valout)
end

end


text:
using UMAP

using Pkg
Pkg.activate(".")
using HypergraphModularity

using Arpack
using BSON
using LinearAlgebra
using SparseArrays
using Random

function regularized_spectral_embedding(A; tol=1e-10, maxiter=300, nev=3)
    dC = vec(sum(A,dims=1))
    dR = vec(sum(A,dims=2))
    τC = sum(dC) / length(dC)
    τR = sum(dR) / length(dR)
    DC = Diagonal(1.0 ./ sqrt.(dC .+ τC))
    DR = Diagonal(1.0 ./ sqrt.(dR .+ τR))
    N = DR * A * DC
    s = svds(N, nsv=nev, tol=tol, maxiter=maxiter)
    U, V = s[1].U, s[1].V
    return U, V
end

function read_dataset(dataset)
    H, labels = read_hypergraph_data(dataset, 25, 2)
    names = read_hypergraph_label_names(dataset)
    II = Int64[]
    JJ = Int64[]
    edge_ind = 1
    for E in values(H.E), (edge, _) in E
        for v in edge
            push!(II, v)
            push!(JJ, edge_ind)
	end
        edge_ind += 1
    end
    B = sparse(II, JJ, 1, maximum(II), maximum(JJ));
    size(B)
    return H, labels, names, B
end

function gen_embeddings(dataset)
    H, labels, names, B = read_dataset(dataset)
    U, V = regularized_spectral_embedding(B, tol=1e-10, maxiter=500, nev=25)
    xy_coords = umap([U' V'], 2)
    bson("$dataset-xy.bson", Dict("xy" => xy_coords, "U" => U, "V" => V, "B" => B, "labels" => labels, "names" => names))
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
using RestrictProlong, Test

@testset "restrict" begin
    oddset =  (([1,0,0,0,0], [0.5,0,0]),
               ([0,1,0,0,0], [0.25,0.25,0]),
               ([0,0,1,0,0], [0,0.5,0]),
               ([0,0,0,1,0], [0,0.25,0.25]),
               ([0,0,0,0,1], [0,0,0.5]))
    evenset = (([1,0,0,0,0,0], [0.5,0,0]),
               ([0,1,0,0,0,0], [3/8,1/8,0]),
               ([0,0,1,0,0,0], [1/8,3/8,0]),
               ([0,0,0,1,0,0], [0,3/8,1/8]),
               ([0,0,0,0,1,0], [0,1/8,3/8]),
               ([0,0,0,0,0,1], [0,0,0.5]))

    # 1D
    for (y, Y) in oddset
        @test restrict(y) ≈ Y
        @test restrict(y, 1) ≈ Y
        @test restrict(y, 2) == y
    end
    for (y, Y) in evenset
        @test restrict(y) ≈ Y
        @test restrict(y, 1) ≈ Y
        @test restrict(y, 2) == y
    end

    # 2D
    for yset in (oddset, evenset)
        for xset in (oddset, evenset)
            for (y, Y) in yset
                for (x, X) in xset
                    @test restrict(x.*y') ≈ X.*Y'
                    @test restrict(x.*y', 1) ≈ X.*y'
                    @test restrict(x.*y', 2) ≈ x.*Y'
                    @test restrict(x.*y', [1,3]) ≈ X.*y'
                    @test restrict(x.*y', ()) == x.*y'
                end
            end
        end
    end
end

@testset "prolong" begin
    oddset =  (([1,0,0], [0.5,0.25,0,0,0]),
               ([0,1,0], [0,0.25,0.5,0.25,0]),
               ([0,0,1], [0,0,0,0.25,0.5]))
    evenset = (([1,0,0], [1/2,3/8,1/8,0,0,0]),
               ([0,1,0], [0,1/8,3/8,3/8,1/8,0]),
               ([0,0,1], [0,0,0,1/8,3/8,1/2]))

    # 1D
    for (Y, y) in oddset
        @test prolong(Y) ≈ y
        @test prolong(Y, length(y)) ≈ y
        @test prolong(Y, length(y), 1) ≈ y
        @test prolong(Y, 1, 2) == Y
    end
    for (Y, y) in evenset
        @test prolong(Y, length(y)) ≈ y
        @test prolong(Y, length(y), 1) ≈ y
        @test prolong(Y, 1, 2) == Y
    end

    # 2D
    for yset in (oddset, evenset)
        for xset in (oddset, evenset)
            for (Y, y) in yset
                for (X, x) in xset
                    @test prolong(X.*Y', (length(x),length(y))) ≈ x.*y'
                    @test prolong(X.*Y', length(x), 1) ≈ x.*Y'
                    @test prolong(X.*Y', length(y), 2) ≈ X.*y'
                end
            end
        end
    end
end

nothing


text:
"""
    SeisSort(in, out;<keyword arguments>)

Sort a seis file using its header words

# Arguments
* `in`: input filename >> a text file with information about data extent, data and header file names; a binary file containing data and a binary file containing headers.
* `out`: output filename

# Keyword arguments
* `key=["imx","imy"]`
* `rev=false` : sort headers in decreasing order
* `ntrace=1000` : number of traces to read at a time

# Output
file `out` is created with data sorted.

*Credits: AS, 2015*
"""
function SeisSort(in, out;key=["imx","imy"],rev=false,ntrace=100000)
    filename_h = ParseHeaderName(in)
    stream_h = open(filename_h)
    seek(stream_h, header_count["n1"])
    nt = read(stream_h,Int32)
    nx = convert(Int64,filesize(stream_h)/(4*length(fieldnames(Header))))
    h = Header[]
    # find min and max for each key
    h1 = GrabHeader(stream_h,1)
    minval = Array(Float32,length(key))

    for ikey=1:length(key)
	minval[ikey] = getfield(h1,Symbol(key[ikey]))
    end
    for j=2:nx
	h1 = GrabHeader(stream_h,j)
	for ikey=1:length(key)
	    key_val = abs(getfield(h1,Symbol(key[ikey])))
	    if (key_val < minval[ikey])
		minval[ikey] = key_val
	    end
	end
    end
    mykey = vec(zeros(Float64,nx))
    seekstart(stream_h)
    for j=1:nx
	h1 = GrabHeader(stream_h,j)
	for ikey=1:length(key)
	    mykey[j] += ((getfield(h1,Symbol(key[ikey])) + minval[ikey] + 1)*
                         (10^(6*(length(key)-ikey))))
	end
    end
    close(stream_h)
    p = convert(Array{Int32,1},sortperm(mykey,rev=rev))
    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
    filename_d_out = join([DATAPATH out "@data@"])
    filename_h_out = join([DATAPATH out "@headers@"])
    nhead = length(fieldnames(Header))
    stream_h = open(filename_h_out)
    nx = Int(floor(filesize(stream_h)/(nhead*4)))
    close(stream_h)
    extent = ReadTextHeader(in)
    extent.n2 = nx
    extent.n3 = 1
    extent.n4 = 1
    extent.n5 = 1
    WriteTextHeader(out,extent,"native_float",4,filename_d_out,filename_h_out)
    FetchHeaders(filename_h,out,p,nx)
    Seismic.FetchTraces(in,out)
    tmp = join(["tmp_SeisSort_",string(Int(floor(rand()*100000)))])
    cp(out,tmp,remove_destination=true);
    @compat SeisProcessHeaders(out, tmp, [UpdateHeader],
                               [Dict(:itmin=>1,:itmax=>nt)])
    filename_h_tmp = join([DATAPATH tmp "@headers@"])
    filename_h_out = join([DATAPATH out "@headers@"])
    cp(filename_h_tmp,filename_h_out,remove_destination=true);
    rm(filename_h_tmp);
    rm(tmp);
end

function FetchHeaders(filename_h_in::AbstractString, filename_out::AbstractString,
                      p::Array{Int32,1}, nx)
    stream_h = open(filename_h_in)
    h = Header[]
    for j = 1:nx
	append!(h,[GrabHeader(stream_h,p[j])])
    end
    SeisWriteHeaders(filename_out,h,update_tracenum=false)
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
strptr(kstr::KString) = strptr(kstr.s,1,kstr.l)
strptr(p::Ptr{Int8},st::Integer,ed::Integer) = stringfrompointer(convert(Ptr{UInt8},p),Int64(st),Int64(ed))
strptr(p::Ptr{UInt8},st::Integer,ed::Integer) = stringfrompointer(p,Int64(st),Int64(ed))
function stringfrompointer(p::Ptr{UInt8},st::Int64,ed::Int64)
    len = ed-st+1
    @assert len>0
    data = Array{UInt8,1}(len)
    for i in st:ed
        data[i-st+1] = unsafe_load(p,i)
    end
    ASCIIString(data)
end


text:
#### Bors

"""
    Bors

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Bors
    "bors checks"
    checks::Vector{String}
    "repo ID"
    id::String = ""
end

function toml(bors::Bors)
    checks = bors.checks
    status = join(["\"$(x)\"" for x in checks], ",\n")
    return "status = [
$status
]
delete_merged_branches = true
timeout_sec = 86400
block_labels = [ \"do-not-merge-yet\" ]
cut_body_after = \"<!--\"
"
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates