text:
@doc raw"""
`ParametersDGVIP1`: Parameters for right-hand side function of Discontinuous Galerkin Variational Integrator.

### Parameters

* `DT`: data type
* `TT`: parameter type
* `D`: dimension of the system
* `S`: number of basis nodes
* `R`: number of quadrature nodes

### Fields

* `Θ`:  function of the noncanonical one-form (∂L/∂v)
* `f`:  function of the force (∂L/∂q)
* `g`:  function of the projection ∇ϑ(q)⋅v
* `Δt`: time step
* `b`:  quadrature weights
* `c`:  quadrature nodes
* `m`:  mass matrix
* `a`:  derivative matrix
* `r⁻`: reconstruction coefficients, jump lhs value
* `r⁺`: reconstruction coefficients, jump rhs value
* `t`:  current time
* `q`:  current solution of qₙ
* `q⁻`: current solution of qₙ⁻
* `q⁺`: current solution of qₙ⁺
"""
mutable struct ParametersDGVIP1{DT, TT, D, S, R, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    Δt::TT

    b::Vector{TT}
    c::Vector{TT}
    m::Matrix{TT}
    a::Matrix{TT}
    r⁻::Vector{TT}
    r⁺::Vector{TT}

    t::TT

    q::Vector{DT}
    q⁻::Vector{DT}
    q⁺::Vector{DT}

    function ParametersDGVIP1{DT,D}(equs::ET, Δt::TT, b, c, m, a, r⁻, r⁺) where {DT, TT, D, ET <: NamedTuple}
        @assert length(b)  == length(c)
        @assert length(r⁻) == length(r⁺)
        new{DT,TT,D,length(r⁻),length(c),ET}(equs, Δt, b, c, m, a, r⁻, r⁺, zero(TT), zeros(DT,D), zeros(DT,D), zeros(DT,D))
    end

    function ParametersDGVIP1{DT,D}(equs::NamedTuple, Δt::TT,
                    basis::Basis{TT}, quadrature::QuadratureRule{TT}) where {DT,TT,D}

        # compute coefficients
        b = weights(quadrature)
        c = nodes(quadrature)
        m = zeros(TT, nnodes(quadrature), nbasis(basis))
        a = zeros(TT, nnodes(quadrature), nbasis(basis))
        r⁻= zeros(TT, nbasis(basis))
        r⁺= zeros(TT, nbasis(basis))

        for i in eachindex(basis)
            for j in eachindex(quadrature)
                m[j,i] = basis[nodes(quadrature)[j], i]
                a[j,i] = basis'[nodes(quadrature)[j], i]
            end
            r⁻[i] = basis[one(TT), i]
            r⁺[i] = basis[zero(TT), i]
        end

        if get_config(:verbosity) > 1
            println()
            println("  Discontinuous Galerkin Variational Integrator")
            println("  =============================================")
            println()
            println("    b = ", b)
            println("    c = ", c)
            println("    m = ", m)
            println("    a = ", a)
            println("    r⁻= ", r⁻)
            println("    r⁺= ", r⁺)
            println()
        end

        ParametersDGVIP1{DT,D}(equs, Δt, b, c, m, a, r⁻, r⁺)
    end
end

# function update_params!(params::ParametersDGVIP1, sol::AtomicSolutionPODE)
#     # set time for nonlinear solver and copy previous solution
#     params.t  = sol.t
#     params.q .= sol.q
#     # params.p .= sol.p
# end

function IntegratorCache{ST}(params::ParametersDGVIP1{DT,TT,D,S,R}; kwargs...) where {ST,DT,TT,D,S,R}
    IntegratorCacheDGVI{ST,D,S,R}(; kwargs...)
end

@inline CacheType(ST, params::ParametersDGVIP1{DT,TT,D,S,R}) where {DT,TT,D,S,R} = IntegratorCacheDGVI{ST,D,S,R}

@doc raw"""
`IntegratorDGVIP1`: Discontinuous Galerkin Variational Integrator. *EXPERIMENTAL*

### Parameters

### Fields

* `equation`: Implicit Ordinary Differential Equation
* `basis`: piecewise polynomial basis
* `quadrature`: numerical quadrature rule
* `Δt`: time step
* `params`: ParametersDGVIP1
* `solver`: nonlinear solver
* `iguess`: initial guess
* `q`: current solution vector for trajectory
* `p`: current solution vector for one-form
* `cache`: temporary variables for nonlinear solver
"""
struct IntegratorDGVIP1{DT, TT, D, S, R,
                BT <: Basis,
                PT <: ParametersDGVIP1{DT,TT,D,S},
                ST <: NonlinearSolver{DT},
                IT <: InitialGuessODE{TT}} <: IODEIntegrator{DT,TT}
    basis::BT
    quadrature::QuadratureRule{TT,R}

    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    q::Vector{DT}
    q⁻::Vector{DT}
    q⁺::Vector{DT}

    function IntegratorDGVIP1(basis::BT, quadrature::QuadratureRule{TT,R}, params::ParametersDGVIP1{DT,TT,D,S},
                    solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,BT,ST,IT}
        # create solution vectors
        q  = zeros(DT,D)
        q⁻ = zeros(DT,D)
        q⁺ = zeros(DT,D)

        new{DT, TT, D, S, R, BT, typeof(params), ST, IT}(basis, quadrature, params, solver, iguess, caches, q, q⁻, q⁺)
    end

    function IntegratorDGVIP1{DT,D}(equations::NamedTuple, basis::Basis{TT}, quadrature::QuadratureRule{TT,R}, Δt::TT;
                                  interpolation=HermiteExtrapolation{DT}) where {DT,TT,D,R}

        # get number of stages
        S = nbasis(basis)

        # create params
        params = ParametersDGVIP1{DT,D}(equations, Δt, basis, quadrature)

        # create cache dict
        caches = CacheDict(params)

        # create nonlinear solver
        solver = create_nonlinear_solver(DT, D*(S+2), params, caches)

        # create initial guess
        iguess = InitialGuessODE(get_config(:ig_extrapolation), equations[:v̄], Δt)

        # create integrator
        IntegratorDGVIP1(basis, quadrature, params, solver, iguess, caches)
    end

    function IntegratorDGVIP1(equation::Union{IODE{DT}, LODE{DT}}, basis::Basis, quadrature::QuadratureRule, Δt; kwargs...) where {DT}
        IntegratorDGVIP1{DT, ndims(equation)}(get_functions(equation), basis, quadrature, Δt; kwargs...)
    end
end

@inline equation(integrator::IntegratorDGVIP1, i::Symbol) = integrator.params.equs[i]
@inline equations(integrator::IntegratorDGVIP1) = integrator.params.equs
@inline timestep(integrator::IntegratorDGVIP1) = integrator.params.Δt
@inline Base.ndims(::IntegratorDGVIP1{DT,TT,D}) where {DT,TT,D} = D

function update_params!(params::ParametersDGVIP1, int::IntegratorDGVIP1)
    # set time for nonlinear solver and copy previous solution
    params.t  += int.params.Δt
    params.q  .= int.q
    params.q⁻ .= int.q⁻
    params.q⁺ .= int.q⁺
end

function initialize!(int::IntegratorDGVIP1, sol::AtomicSolutionPODE)
    # copy initial conditions from solution
    int.q  .= sol.q
    int.q⁻ .= int.q
    int.q⁺ .= int.q

    sol.t̄ = sol.t - timestep(int)

    equation(int, :v̄)(sol.t, sol.q, sol.v)

    # initialise initial guess
    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̄, sol.q̄, sol.v̄)
end

function initial_guess!(int::IntegratorDGVIP1{DT,TT, D, S, R}, sol::AtomicSolutionPODE{DT,TT},
                        cache::IntegratorCacheDGVI{DT}=int.caches[DT]) where {DT,TT,D,S,R}
    if nbasis(int.basis) > 0
        for i in eachindex(int.basis)
            evaluate!(int.iguess, sol.q̄, sol.v̄,
                                  sol.q, sol.v,
                                  cache.q̃,
                                  grid(int.basis)[i])

            for k in 1:D
                int.solver.x[D*(i-1)+k] = cache.q̃[k]
            end
        end
    else
        for i in 1:S
            for k in 1:D
                int.solver.x[D*(i-1)+k] = 0
            end
        end
    end

    evaluate!(int.iguess, sol.q̄, sol.v̄,
                          sol.q, sol.v,
                          cache.q̃,
                          one(TT))

    for k in 1:D
        int.solver.x[D*(S+0)+k] = cache.q̃[k]
        int.solver.x[D*(S+1)+k] = cache.q̃[k]
    end
end

"Compute stages of variational partitioned Runge-Kutta methods."
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersDGVIP1{DT,TT,D,S,R},
                caches::CacheDict) where {ST,DT,TT,D,S,R}
    @assert length(x) == length(b)

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache, params)

    # compute rhs b of nonlinear solver
    compute_rhs!(b, cache, params)
end

function compute_stages!(x, cache::IntegratorCacheDGVI{ST,D,S}, params::ParametersDGVIP1{DT,TT,D,S}) where {ST,DT,TT,D,S}
    # copy x to X
    for i in 1:S
        for k in 1:D
            cache.X[i][k] = x[D*(i-1)+k]
        end
    end

    # copy x to q̄=qₙ+₁
    for k in 1:D
        cache.q̄[k] = x[D*S+k]
    end

    # copy x to q̄⁺=qₙ+₁⁺
    for k in 1:D
        cache.q̄⁺[k] = x[D*(S+1)+k]
    end

    # compute Q, qₙ⁺ and qₙ₊₁⁻
    compute_stages_q!(cache, params)

    # compute V
    compute_stages_v!(cache, params)

    # compute P and F
    compute_stages_p!(cache, params)

    # compute jump
    compute_stages_λ!(cache, params)
end

"Compute solution at quadrature nodes and across jump."
function compute_stages_q!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local q::ST
    local q⁺::ST
    local q̄⁻::ST

    local X = cache.X
    local Q = cache.Q

    # copy q and q⁻
    cache.q  .= params.q
    cache.q⁻ .= params.q⁻

    # compute Q
    for i in 1:R
        for k in 1:D
            q = 0
            for j in 1:S
                q += params.m[i,j] * X[j][k]
            end
            Q[i][k] = q
        end
    end

    # compute qₙ⁺ and qₙ₊₁⁻
    for k in 1:D
        q⁺ = 0
        q̄⁻ = 0
        for i in 1:S
            q⁺ += params.r⁺[i] * X[i][k]
            q̄⁻ += params.r⁻[i] * X[i][k]
        end
        cache.q⁺[k] = q⁺
        cache.q̄⁻[k] = q̄⁻
    end
end

"Compute velocities at quadrature nodes."
function compute_stages_v!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local v::ST

    for i in 1:R
        for k in 1:D
            v = 0
            for j in 1:S
                v += params.a[i,j] * cache.X[j][k]
            end
            cache.V[i][k] = v / params.Δt
        end
    end
end

"Compute one-form and forces at quadrature nodes."
function compute_stages_p!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local tᵢ::TT

    # compute P=ϑ(Q) and F=f(Q)
    for i in 1:R
        tᵢ = params.t + params.Δt * params.c[i]
        params.equs[:ϑ](tᵢ, cache.Q[i], cache.V[i], cache.P[i])
        params.equs[:f](tᵢ, cache.Q[i], cache.V[i], cache.F[i])
    end
end

function compute_stages_λ!(cache::IntegratorCacheDGVI{ST,D,S,R},
                           params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local t₀::TT = params.t
    local t₁::TT = params.t + params.Δt

    # compute ϕ and ϕ̅
    cache.ϕ  .= 0.5 * (cache.q⁻ .+ cache.q⁺)
    cache.ϕ⁻ .= 0.5 * (cache.q⁻ .+ cache.q )
    cache.ϕ⁺ .= 0.5 * (cache.q  .+ cache.q⁺)

    cache.ϕ̅  .= 0.5 * (cache.q̄⁻ .+ cache.q̄⁺)
    cache.ϕ̅⁻ .= 0.5 * (cache.q̄⁻ .+ cache.q̄ )
    cache.ϕ̅⁺ .= 0.5 * (cache.q̄  .+ cache.q̄⁺)

    # compute λ and λ̄
    cache.λ  .= cache.q⁺ .- cache.q⁻
    cache.λ⁻ .= cache.q  .- cache.q⁻
    cache.λ⁺ .= cache.q⁺ .- cache.q

    cache.λ̄ .= cache.q̄⁺ .- cache.q̄⁻
    cache.λ̄⁻ .= cache.q̄  .- cache.q̄⁻
    cache.λ̄⁺ .= cache.q̄⁺ .- cache.q̄

    # compute ϑ
    params.equs[:ϑ](t₀, cache.q,  cache.q,  cache.θ)
    params.equs[:ϑ](t₀, cache.q⁻, cache.q⁻, cache.θ⁻)
    params.equs[:ϑ](t₀, cache.q⁺, cache.q⁺, cache.θ⁺)

    params.equs[:ϑ](t₁, cache.q̄,  cache.q̄,  cache.Θ̅)
    params.equs[:ϑ](t₁, cache.q̄⁻, cache.q̄⁻, cache.Θ̅⁻)
    params.equs[:ϑ](t₁, cache.q̄⁺, cache.q̄⁺, cache.Θ̅⁺)

    # compute projection
    params.equs[:g](t₀, cache.q,  cache.λ,  cache.g)
    params.equs[:g](t₀, cache.q⁻, cache.λ⁻, cache.g⁻)
    params.equs[:g](t₀, cache.q⁺, cache.λ⁺, cache.g⁺)

    params.equs[:g](t₁, cache.q̄,  cache.λ̄,  cache.ḡ)
    params.equs[:g](t₁, cache.q̄⁻, cache.λ̄⁻, cache.ḡ⁻)
    params.equs[:g](t₁, cache.q̄⁺, cache.λ̄⁺, cache.ḡ⁺)

    # # compute ϑ
    # params.equs[:ϑ](t₀, cache.ϕ⁻, cache.ϕ⁻, cache.θ⁻)
    # params.equs[:ϑ](t₀, cache.ϕ⁺, cache.ϕ⁺, cache.θ⁺)
    # params.equs[:ϑ](t₁, cache.ϕ̅⁻, cache.ϕ̅⁻, cache.Θ̅⁻)
    #
    # # compute projection
    # params.equs[:g](t₀, cache.ϕ⁻, cache.λ⁻, cache.g⁻)
    # params.equs[:g](t₀, cache.ϕ⁺, cache.λ⁺, cache.g⁺)
    # params.equs[:g](t₁, cache.ϕ̅⁻, cache.λ̄⁻, cache.ḡ⁻)
end

function compute_rhs!(b::Vector{ST}, cache::IntegratorCacheDGVI{ST,D,S,R},
                params::ParametersDGVIP1{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}

    local z::ST

    # compute b = - [(P-AF)]
    for i in 1:S
        for k in 1:D
            z = 0
            for j in 1:R
                z += params.b[j] * params.m[j,i] * cache.F[j][k] * params.Δt
                z += params.b[j] * params.a[j,i] * cache.P[j][k]
            end

            z += params.r⁺[i] * 0.5 * ( cache.θ[k] + cache.θ⁺[k] )
            z -= params.r⁻[i] * 0.5 * ( cache.Θ̅[k] + cache.Θ̅⁻[k] )

            z += params.r⁺[i] * 0.5 * cache.g⁺[k]
            z += params.r⁻[i] * 0.5 * cache.ḡ⁻[k]

            b[D*(i-1)+k] = z
        end
    end

    # compute b = qₙ⁺ - r⁺⋅X
    for k in 1:D
        b[D*S+k] = params.q⁺[k] - cache.q⁺[k]
    end

    for k in 1:D
        b[D*(S+1)+k] = cache.Θ̅⁺[k] - cache.Θ̅⁻[k] - cache.ḡ[k]
    end
end

function update_solution!(int::IntegratorDGVIP1{DT}, cache::IntegratorCacheDGVI{DT}) where {DT}
    int.q  .= cache.q̄
    int.q⁻ .= cache.q̄⁻
    int.q⁺ .= cache.q̄⁺
end

function integrate_step!(int::IntegratorDGVIP1{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCacheDGVI{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, int)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # copy solution from cache to integrator
    update_solution!(int, cache)
    sol.q = int.q

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
    # update_vector_fields!(int.iguess, sol.t, sol.q, int.q⁺)

    # take care of periodic solutions
    # cut_periodic_solution!(int.q,  int.equation.periodicity)
    # cut_periodic_solution!(int.q⁻, int.equation.periodicity)
    # cut_periodic_solution!(int.q⁺, int.equation.periodicity)

    # copy to solution
    # copy_solution!(sol, int.q, int.q⁺)
end


text:
"""
    read_walz(file)

Import a Walz GFS-3000 output file.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(PlantBiophysics))),"test","inputs","data","P1F20129.csv")
read_walz(file)
```
"""
function read_walz(file)
    df = CSV.read(file, DataFrame, header = 1, datarow = 3)

    if hasproperty(df, :Ttop)
        rename!(df, :Ttop => :Tmin)
    end

    df[!,:Comment] = locf(df[!,:Comment])
    dropmissing!(df, :VPD)

    # Renaming variables to fit the standard in the package:
    rename!(
        df,
        :GH2O => :gs, :ca => :Cₐ, :Tcuv => :T, :Pamb => :P, :rh => :Rh,
        :PARtop => :PPFD, :ci => :Cᵢ, :Comment => :curve
    )

    # Recomputing the variables to fit the units used in the package:
    df[!,:VPD] = round.(df[:,:VPD] .* df[:,:P] ./ 1000.0, digits = 3)
    df[!,:gs] = round.(gsw_to_gsc.(df[:,:gs]) ./ 1000.0, digits = 5)
    df[!,:AVPD] = df[:,:A] ./ (df[:,:Cₐ] .* sqrt.(df[:,:VPD]))
    df[!,:Rh] = df[!,:Rh] ./ 100.0

    return df
end

ext: jl


text:
"""
    PressureLevel <: ERA5Variable

Abstract supertype for Pressure-Level variables.  Contains the following fields:
- `varID` : The variable ID, that is also the identifier in the NetCDF files
- `lname` : The variable long-name, which is used to specify retrievals from CDS
- `vname` : The full-name of the variable
- `units` : The units of the variable
- `hPa`   : The pressure-level height of the variable in concern
"""
abstract type PressureLevel <: ERA5Variable end

"""
    PressureVariable <: PressureLevel

Subtype for Pressure-Level variables that can be directly retrieved from the CDS
"""
struct PressureVariable{ST<:AbstractString} <: PressureLevel
    varID :: ST
    lname :: ST
    vname :: ST
    units :: ST
    hPa   :: Int
end

"""
    PressureCustom <: PressureLevel

Subtype for custom user-defined Pressure-Level variables
"""
struct PressureCustom{ST<:AbstractString} <: PressureLevel
    varID :: ST
    lname :: ST
    vname :: ST
    units :: ST
    hPa   :: Int
end

"""
    PressureVariable(
        varID :: AbstractString,
        ST = String;
        hPa   :: Int = 0
        throw :: Bool = true
    ) -> evar :: SingleLevel

Retrieve the basic properties of the Pressure-Level variable defined by `varID` at pressure-height indicated by `hPa` and put them in the `evar` SingleLevel type structure.

Arguments
=========

- `varID` : variable ID (in string format) used in the NetCDF file

Keyword Arguments
=================

- `hPa` : Integer specifying pressure-level height in hPa
- `throw` : if `hPa` level does not exist and `throw` is true, throw error, otherwise find nearest pressure level
"""
function PressureVariable(
    varID :: AbstractString,
    ST = String;
    hPa   :: Int = 0,
    throw :: Bool = true
)

    isPressure(varID)

    @info "$(modulelog()) - Retrieving information for the PressureVariable defined by the ID \"$varID\""
    vlist,flist = listPressures();    ind = findall(varID.==vlist)[1]
    vtype = replace(flist[ind],".txt"=>"")
    fname = joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis",flist[ind])
    vlist = listera5variables(fname); ind = findall(varID.==vlist)[1]

    IDinfo = readdlm(fname,',',comments=true,comment_char='#')[ind,:]
    varID,lname,vname,units = IDinfo[[1,2,3,4]]

    if !iszero(hPa)
        prelist = era5Pressures()
        if sum(prelist.==hPa) != 1
            if throw
                error("$(modulelog()) - Pressure level specified in \"hPa\" argument is invalid, please check and see if you requested correctly")
            else
                @warn "$(modulelog()) - Pressure level specified in \"hPa\" argument does not exist, snapping to nearest pressure level"
                ipre = argmin(abs.(prelist.-hPa))
                hPa  = prelist[ipre]
            end
        end
    end

    if vtype == "pressurevariable"
        @info "$(modulelog()) - The ERA5Variable defined by \"$varID\" is of the PressureVariable type"
        return PressureVariable{ST}(varID,lname,vname,units,hPa)
    else
        @info "$(modulelog()) - The ERA5Variable defined by \"$varID\" is of the PressureCustom type"
        return PressureCustom{ST}(varID,lname,vname,units,hPa)
    end

end

"""
    PressureVariable(
        ST = String;
        varID :: AbstractString,
        lname :: AbstractString = "",
        vname :: AbstractString,
        units :: AbstractString,
        hPa   :: Int = 0,
        throw :: Bool = true
    ) -> evar :: PressureCustom

Create a custom Pressure-Level variable that is not in the default list exported by ERA5Reanalysis.jl.  These variables are not available in the CDS store, and so they must be separately calculated from other variables and analyzed.

Keyword Arguments
=================

- `varID` : variable ID (in string format) used in the NetCDF file
- `lname` : long-name for variable (used in specifying variable for CDS downloads)
- `vname` : user-defined variable name
- `units` : user-defined units of the variable
- `hPa`   : Pressure level specified in hPa. Default is 0, which indicates all levels.
- `throw` : if `hPa` level does not exist and `throw` is true, throw error, otherwise find nearest pressure level
"""
function PressureVariable(
    ST = String;
    varID :: AbstractString,
    lname :: AbstractString = "",
    vname :: AbstractString,
    units :: AbstractString,
    hPa   :: Int = 0,
    throw :: Bool = true
)

    if isPressure(varID,throw=false)
        error("$(modulelog()) - The PressureVariable \"$(varID)\" has already been defined, please use another identifier.")
    else
        @info "$(modulelog()) - Adding the PressureVariable \"$(varID)\" to the list."
    end

    open(joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis","pressurecustom.txt"),"a") do io
        write(io,"$varID,$lname,$vname,$units\n")
    end

    if !iszero(hPa)
        prelist = era5Pressures()
        if sum(prelist.==hPa) != 1
            if throw
                error("$(modulelog()) - Pressure level specified in \"hPa\" argument is invalid, please check and see if you requested correctly")
            else
                @warn "$(modulelog()) - Pressure level specified in \"hPa\" argument does not exist, snapping to nearest pressure level"
                ipre = argmin(abs.(prelist.-hPa))
                hPa  = prelist[ipre]
            end
        end
    end

    return PressureCustom{ST}(varID,lname,vname,units,hPa)

end

"""
    listPressures()

List all Pressure-Level Variables and the files the data are stored in.

Output
======

- `varlist` : List of all the Pressure-Level Variable IDs
- `fidlist` : List of the files that the Pressure-Level Variable information is stored in
"""
function listPressures()

    flist   = ["pressurevariable.txt","pressurecustom.txt"]
    varlist = []
    fidlist = []

    for fname in flist
        copyera5variables(fname)
        vvec = listera5variables(joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis",fname))
        varlist = vcat(varlist,vvec); fvec = fill(fname,length(vvec))
        fidlist = vcat(fidlist,fvec)
    end

    return varlist,fidlist

end

"""
    isPressure(
        varID :: AbstractString;
        throw :: Bool = true
    ) -> tf :: Bool

Extracts information of the Pressure-Level Variable with the ID `varID`.  If no Pressure-Level Variable with this ID exists, an error is thrown.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the Pressure-Level Variable.
        If the ID is not valid (i.e. not being used), then an error will be thrown.
- `throw` : If `true`, then throws an error if `RegID` is not a valid Pressure-Level Variable identifier instead of returning the Boolean `tf`

Returns
=======

- `tf` : True / False
"""
function isPressure(
    varID :: AbstractString;
    throw :: Bool = true
)

    @info "$(modulelog()) - Checking to see if the PressureVariable ID \"$varID\" is in use"
    vlist,_ = listPressures()

    if sum(vlist.==varID) == 0
        if throw
            error("$(modulelog()) - \"$(varID)\" is not a valid PressureVariable identifier, use the function PressureCustom() to add this ERA5Variable to the list.")
        else
            @warn "$(modulelog()) - \"$(varID)\" is not a valid PressureVariable identifier, use the function PressureCustom() to add this ERA5Variable to the list."
            return false
        end
    else
        @info "$(modulelog()) - The PressureVariable ID \"$varID\" is already in use"
        return true
    end

end

"""
    rmPressure( varID :: AbstractString ) -> nothing

Remove the Pressure-Level Variable with the ID `varID` from the lists.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the Pressure-Level Variable that is to be removed
"""
function rmPressure(varID::AbstractString)

    if isPressure(varID,throw=false)
        disable_logging(Logging.Warn)
        rmERA5Variable(PressureVariable(varID))
        disable_logging(Logging.Debug)
        @info "$(modulelog()) - Successfully removed the Pressure-Level variable defined by \"$(varID)\""
    else
        @warn "$(modulelog()) - No Pressure-Level variable defined by \"$(varID)\" exists, please make sure you specified the correct variable ID"
    end

    return nothing

end

"""
    resetPressures() -> nothing

Reset the list of Pressure-Level variables to the ERA5Reanalysis default.
"""
function resetPressures()

    @info "$(modulelog()) - Resetting the custom list of ERA5 PressureVariables back to the default"
    flist = ["pressurecustom.txt"]

    for fname in flist
        copyera5variables(fname,overwrite=true)
    end

    return nothing

end

"""
    era5Pressures() -> parray :: Vector{Int}

Returns the a vector containing the 37 pressure levels available in ERA5 in hPa units.

Returns
=======

- `parray` : vector containing list of pressures in Int format and hPa units
"""
era5Pressures() = [
    1,2,3,5,7,10,20,30,50,70,100,125,150,175,200,
    225,250,300,350,400,450,500,550,600,650,700,750,
    775,800,825,850,875,900,925,950,975,1000
]

function tablePressures()

    jfol = joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis"); mkpath(jfol);
    fvar = ["PressureVariable","PressureCustom"]
    fmat = []

    for fname in fvar
        fid  = joinpath(jfol,"$(lowercase(fname)).txt")
        try
            vmat = readdlm(fid,',',comments=true,comment_char='#')
            nvar = size(vmat,1); ff = fill(fname,nvar)
            vmat = cat(ff,vmat[:,[1,3,4,2]],dims=2)
            fmat = cat(fmat,vmat,dims=1)
        catch
        end
    end

    head = ["Variable Type","ID","Name","Units","ERA5 Long-Name"];

    pretty_table(
        fmat,head,
        alignment=[:c,:c,:l,:c,:l],
        crop = :none, tf = tf_compact
    );

end


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays

@testset "embeddedmap" begin
    m = 6; n = 5
    M = 10(1:m) .+ (1:n)'; L = LinearMap(M)
    offset = (3,4)

    BM = [zeros(offset...) zeros(offset[1], size(M,2));
        zeros(size(M,1), offset[2]) M]
    BL = @inferred LinearMap(L, size(BM); offset=offset)
    s1, s2 = size(BM)
    @test (@inferred Matrix(BL)) == BM == mul!(zero(BM), BL, 1)
    @test mul!(copy(BM), BL, 2, true, true) == 3BM
    @test (@inferred Matrix(BL')) == BM'
    @test (@inferred Matrix(transpose(BL))) == transpose(BM)

    @test_throws UndefKeywordError LinearMap(M, (10, 10))
    @test_throws ArgumentError LinearMap(M, (m, n), (0:m, 1:n))
    @test_throws ArgumentError LinearMap(M, (m, n), (0:m-1, 1:n))
    @test_throws ArgumentError LinearMap(M, (m, n), (1:m, 1:n+1))
    @test_throws ArgumentError LinearMap(M, (m, n), (1:m, 2:n+1))
    @test_throws ArgumentError LinearMap(M, (m, n), offset=(3,3))
    # @test_throws ArgumentError LinearMap(M, (m, n), (m:-1:1, 1:n))
    # @test_throws ArgumentError LinearMap(M, (m, n), (collect(m:-1:1), 1:n))
    @test size(@inferred LinearMap(M, (2m, 2n), (1:2:2m, 1:2:2n))) == (2m, 2n)
    @test @inferred !ishermitian(BL)
    @test @inferred !issymmetric(BL)
    @test @inferred LinearMap(L, size(BM), (offset[1] .+ (1:m), offset[2] .+ (1:n))) == BL
    Wc = @inferred LinearMap([2 im; -im 0]; ishermitian=true)
    Bc = @inferred LinearMap(Wc, (4,4); offset=(2,2))
    @test (@inferred ishermitian(Bc))

    x = randn(s2); X = rand(s2, 3)
    y = BM * x; Y = zeros(s1, 3)

    @test @inferred BL * x ≈ BM * x
    @test @inferred BL' * y ≈ BM' * y

    for α in (true, false, rand()),
            β in (true, false, rand()),
            t in (identity, adjoint, transpose)

        @test t(BL) * x ≈ mul!(copy(y), t(BL), x) ≈ t(BM) * x
        @test Matrix(t(BL) * X) ≈ mul!(copy(Y), t(BL), X) ≈ t(BM) * X
        y = randn(s1); Y = randn(s1, 3)
        @test (@inferred mul!(copy(y), t(BL), x, α, β)) ≈ mul!(copy(y), t(BM), x, α, β)
        @test (@inferred mul!(copy(Y), t(BL), X, α, β)) ≈ mul!(copy(Y), t(BM), X, α, β)
    end

    if VERSION >= v"1.8"
        M = rand(3,4)
        L = LinearMap(M)
        @test Matrix(reverse(L)) == reverse(M)
        for dims in (1, 2, (1,), (2,), (1, 2), (2, 1), :)
            @test Matrix(reverse(L, dims=dims)) == reverse(M, dims=dims)
        end
    end
end


text:
"""
Normalised mean square error on predictions
"""
function NMSE(test::Data, gp::GaussianProcess)
    yp, _ = posterior(test.X, gp)
    y = reduce(vcat, test.Y)

    var_y = sum((y .- (sum(y) / 200)).^2)
    N = size(y)[1]

    (100 / (var_y * N)) * sum((yp .- y).^2)
end


text:
using LookupTables
using Test

@testset "LookupTables.jl" begin
    nx = 100
    ny = 50

    y = range(13.7, step = 0.5, length = 50)
    x = range(0.0, step = 0.1, length = 100)

    vals = zeros(nx, ny)

    for i in eachindex(x)
        for j in eachindex(y)
            vals[i, j] = x[i]^2 - 3 * x[i] * y[j]
        end
    end

    ltable = LookupTable2D(x, y, vals)

    using BenchmarkTools

    @btime lookup(ltable, 4.4, 22.2)

    @btime [[lookup(ltable, xval, yval) for xval in x] for yval in y]

end


text:
__precompile__() # this module is safe to precompile
module SolarRotation
using PyCall

using ..SunAsAStar
using ..SunAsAStar: AstropyCoordinates, AstropyTime, Barycorrpy

using Dates, LinearAlgebra
using DataFrames

const SciPySpatialTransform = PyNULL()

function __init__()
	# Uncomment below and remove pyimport from CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic if start using it
	copy!(SciPySpatialTransform , pyimport("scipy.spatial.transform") )
	@assert SciPySpatialTransform != PyNULL()
end

"""
  `CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd; loc )`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
Inputs:
	jd: Julian Date (scalar)
	loc: Astropy Earth Location object. https://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html

Output:
	Delta:  F_obs**2 - F_sid**2 [(km/s)^2]

Based on Colier Cameron et al. (2019)
Uses ses JPL ephemeris rather than assuming a Keplerian orbit
Performs calculations in ecliptic reference frame.

Python version by Shubham Kanodia March 2021 (with input from Jason Wright and Eric Ford)
Translated to Julia by Eric Ford March 9, 2021.
"""
function CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd::Real; loc::PyObject )
#function CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd::Union{Real,Vector{<:Real}}; loc::PyObject )
#function CalculateFWHMDifference_SolarRotationGen(loc::PyObject, jd::AVR ) where { T<:Real, AVR<:AbstractVector{T} } # JDUTC::PyObject)
	# Functions and constants from AstroPy or SciPy that we'll use multiple times
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	get_body_barycentric = AstropyCoordinates.get_body_barycentric
	speed_of_light = Barycorrpy.PhysicalConstants.c
	R_sun = 695700  # km # should we use better, but non-standard value, e.g, 695660?

	SciPySpatialTransform = pyimport("scipy.spatial.transform")
	SciPyRotation = SciPySpatialTransform.Rotation

	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	################################
	######EARTH EPHEMERIS ############
	################################
	##### NUTATION, PRECESSION, ETC. #####

	# Observatory position wrt Geocenter
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
	r_eci = r_pint'  # [m]
	v_eci = v_pint'  # [m/s]

	##### EPHEMERIDES #####
	ephemeris = "de430"
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris) # [km]
	r_geo = earth_geo[1].xyz*1000.0 # [m]
	v_geo = earth_geo[2].xyz*1000.0/86400.0  # [m/s]

	PosVector_EarthSSB = r_eci + r_geo # [m]MethodError

	# Relativistic Addition of Velocities
	VelVector_EarthSSB = (v_eci + v_geo) / (1.0 + dot(v_eci,v_geo)/speed_of_light^2) # [m/s]

	################################
	######SOLAR EPHEMERIS ############
	################################
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)
	PosVector_SolSSB = solar_ephem[1].xyz*1000. #[m]
	VelVector_SolSSB = solar_ephem[2].xyz*1000.0/86400.0  # [m/s]

	################################
	####EQUATORIAL COORD VECTORS ####
	################################
	PosVector_EarthSol = PosVector_EarthSSB - PosVector_SolSSB
	PosMag_EarthSol = norm(PosVector_EarthSol)

	VelVector_EarthSol = (VelVector_EarthSSB - VelVector_SolSSB) / (1.0 + dot(VelVector_SolSSB,VelVector_EarthSSB)/speed_of_light^2)

	#OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	#if typeof(jd) <: Real
		OmegaVector = cross(vec(PosVector_EarthSol), vec(VelVector_EarthSol)) ./ dot(PosVector_EarthSol,PosVector_EarthSol)
	#else
	#	OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	#end
	################################
	################################

	#DeltaCentury = map(i->(Dates.year(JDUTC.datetime) .- 2000)./100, 1:length(jd) )
	DeltaCentury = (Dates.year(JDUTC.datetime) - 2000)/100

	# https://www2.mps.mpg.de/homes/fraenz/systems/systems3art.pdf

	#Eqn 14
	SolarInclination = 7.25
	SolarLongitude = 75.76 + 1.397*DeltaCentury
	EclipticEpsilon = 23.44

	# Need to perform Extrinsic Euler Rotations to rotate from equatorial to ecliptic
	#REcliptic = R.from_euler("X",  -EclipticEpsilon, degrees=true).as_matrix()
	REcliptic = SciPyRotation.from_euler("X",  -EclipticEpsilon, degrees=true).as_matrix()
	#REcliptic = SciPyRotation.from_euler("X",  -EclipticEpsilon, degrees=true).as_dcm()
	# may need to change to .as_matrix() with SciPy 1.4?

	# Intrinsic rotation to go from solar axis to ecliptic
	#RObliquity = R.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_matrix()
	RObliquity = SciPyRotation.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_matrix()
	#RObliquity = SciPyRotation.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_dcm()

	################################
	####ROTATED COORD VECTORS ####
	################################

	#RotatedPositionVector = map(i->REcliptic * PosVector_EarthSol[:,i], 1:size(PosVector_EarthSol,2) )
	#RotatedPositionHat = RotatedPositionVector ./ norm.(RotatedPositionVector)
	RotatedPositionVector = REcliptic * PosVector_EarthSol
	RotatedPositionHat = RotatedPositionVector / norm(RotatedPositionVector)

	#OmegaEarthVectorRotated = map(i->REcliptic * OmegaVector[i],1:size(OmegaVector,2))
	OmegaEarthVectorRotated = REcliptic * OmegaVector
	################################
	################################

	OmegaSolVector = [0.0, 0.0,  2.972 *1e-6]
	#OmegaSolHat = [0.0 , 0.0 , 1.0]

	# Rotated Solar rotation vector to ecliptic plane
	#OmegaSolVector = map(i->RObliquity * OmegaSolVector[:,i], 1:size(OmegaSolVector,2) )
	#OmegaSolHat = OmegaSolVector./norm.(OmegaSolVector)
	OmegaSolVector = RObliquity * OmegaSolVector
	OmegaSolHat = OmegaSolVector/norm(OmegaSolVector)

	sini = sqrt(1 - dot(OmegaSolHat,RotatedPositionHat)^2)
	Gamma = 1.04
	DeltaOmega = OmegaSolVector .- OmegaEarthVectorRotated

	Delta = ((Gamma* R_sun)^2) * (dot(DeltaOmega,DeltaOmega)*sini^2 - dot(OmegaSolVector,OmegaSolVector))
end

function CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd::Vector{<:Real}; loc::PyObject )
#function CalculateFWHMDifference_SolarRotationGen(loc::PyObject, jd::AVR ) where { T<:Real, AVR<:AbstractVector{T} } # JDUTC::PyObject)
	# Functions and constants from AstroPy or SciPy that we'll use multiple times
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	get_body_barycentric = AstropyCoordinates.get_body_barycentric
	speed_of_light = Barycorrpy.PhysicalConstants.c
	R_sun = 695700  # km # should we use better, but non-standard value, e.g, 695660?

	SciPySpatialTransform = pyimport("scipy.spatial.transform")
	SciPyRotation = SciPySpatialTransform.Rotation

	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	################################
	######EARTH EPHEMERIS ############
	################################
	##### NUTATION, PRECESSION, ETC. #####

	# Observatory position wrt Geocenter
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
	r_eci = r_pint'  # [m]
	v_eci = v_pint'  # [m/s]

	##### EPHEMERIDES #####
	ephemeris = "de430"
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris) # [km]
	r_geo = earth_geo[1].xyz*1000.0 # [m]
	v_geo = earth_geo[2].xyz*1000.0/86400.0  # [m/s]

	PosVector_EarthSSB = r_eci + r_geo # [m]MethodError

	# Relativistic Addition of Velocities
	VelVector_EarthSSB = (v_eci + v_geo) / (1.0 + dot(v_eci,v_geo)/speed_of_light^2) # [m/s]

	################################
	######SOLAR EPHEMERIS ############
	################################
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)
	PosVector_SolSSB = solar_ephem[1].xyz*1000. #[m]
	VelVector_SolSSB = solar_ephem[2].xyz*1000.0/86400.0  # [m/s]

	################################
	####EQUATORIAL COORD VECTORS ####
	################################
	PosVector_EarthSol = PosVector_EarthSSB - PosVector_SolSSB
	PosMag_EarthSol = norm(PosVector_EarthSol)

	VelVector_EarthSol = (VelVector_EarthSSB - VelVector_SolSSB) / (1.0 + dot(VelVector_SolSSB,VelVector_EarthSSB)/speed_of_light^2)

	OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	################################
	################################

	DeltaCentury = map(i->((JDUTC.datetime[i].year) .- 2000)./100, 1:length(jd) )
	#DeltaCentury = map(i->(Dates.year(JDUTC.datetime) .- 2000)./100, 1:length(jd) )
	#DeltaCentury = (Dates.year(JDUTC.datetime) - 2000)/100

	# https://www2.mps.mpg.de/homes/fraenz/systems/systems3art.pdf

	#Eqn 14
	SolarInclination = 7.25
	SolarLongitude = 75.76 .+ 1.397*DeltaCentury
	EclipticEpsilon = 23.44

	# Need to perform Extrinsic Euler Rotations to rotate from equatorial to ecliptic
	#REcliptic = R.from_euler("X",  -EclipticEpsilon, degrees=true).as_matrix()
	#REcliptic = SciPyRotation.from_euler("X",  -EclipticEpsilon, degrees=true).as_dcm()
	REcliptic = SciPyRotation.from_euler("X",  -EclipticEpsilon, degrees=true).as_dcm()
	# may need to change to .as_matrix() with SciPy 1.4?

	# Intrinsic rotation to go from solar axis to ecliptic
	#RObliquity = R.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_matrix()
	#RObliquity = SciPyRotation.from_euler("xz",  [SolarInclination, SolarLongitude], degrees=true).as_dcm()
	RObliquity = map(i->SciPyRotation.from_euler("xz",  [SolarInclination, SolarLongitude[i]], degrees=true).as_dcm(),1:length(SolarLongitude))

	################################
	####ROTATED COORD VECTORS ####
	################################

	RotatedPositionVector = map(i->REcliptic * PosVector_EarthSol[:,i], 1:size(PosVector_EarthSol,2) )
	RotatedPositionHat = RotatedPositionVector ./ norm.(RotatedPositionVector)
	#RotatedPositionVector = REcliptic * PosVector_EarthSol
	#RotatedPositionHat = RotatedPositionVector / norm(RotatedPositionVector)

	OmegaEarthVectorRotated = map(i->REcliptic * OmegaVector[i],1:size(OmegaVector,2))
	#OmegaEarthVectorRotated = REcliptic * OmegaVector
	################################
	################################

	OmegaSolVector = [0.0, 0.0,  2.972 *1e-6]
	#OmegaSolHat = [0.0 , 0.0 , 1.0]

	# Rotated Solar rotation vector to ecliptic plane
	OmegaSolVector = map(i->RObliquity[i] * OmegaSolVector[:,i], 1:size(OmegaSolVector,2) )
	OmegaSolHat = OmegaSolVector./norm.(OmegaSolVector)
	#OmegaSolVector = RObliquity * OmegaSolVector
	#OmegaSolHat = OmegaSolVector/norm(OmegaSolVector)

	sini = sqrt.(1 .- dot.(OmegaSolHat,RotatedPositionHat).^2)
	Gamma = 1.04
	DeltaOmega = OmegaSolVector .- OmegaEarthVectorRotated

	Delta = ((Gamma* R_sun)^2) * (dot.(DeltaOmega,DeltaOmega).*sini.^2 .- dot.(OmegaSolVector,OmegaSolVector))
end

"""
   `CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd; loc )`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
Inputs:
  jd: Julian Date (scalar)
  loc: Astropy Earth Location object. https://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html

Output:
  Delta:  F_obs**2 - F_sid**2 [(km/s)^2]

Based on Colier Cameron et al. (2019)
Uses ses JPL ephemeris rather than assuming a Keplerian orbit.
Performs calculations in equatorial reference frame.

Python version by Shubham Kanodia (with input from Jason Wright and Eric Ford) March 2021
Translated to Julia by Eric Ford March 9, 2021.
"""
function CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd::Real; loc::PyObject )
	# Functions and constants from AstroPy or SciPy that we'll use multiple times
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	get_body_barycentric = AstropyCoordinates.get_body_barycentric
	speed_of_light = Barycorrpy.PhysicalConstants.c
	R_sun = 695700  # km # should we use better, but non-standard value, e.g, 695660?

	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	#return JDUTC
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	################################
	######EARTH EPHEMERIS ############
	################################
	##### NUTATION, PRECESSION, ETC. #####

	# Observatory position wrt Geocenter
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
	r_eci = r_pint'  # [m]
	v_eci = v_pint'  # [m/s]

	##### EPHEMERIDES #####
	ephemeris = "de430"
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris) # [km]
	r_geo = earth_geo[1].xyz*1000.0 # [m]
	v_geo = earth_geo[2].xyz*1000.0/86400.0  # [m/s]

	PosVector_EarthSSB = r_eci + r_geo # [m]

	# Relativistic Addition of Velocities
	VelVector_EarthSSB = (v_eci + v_geo) / (1.0 + dot(v_eci,v_geo)/speed_of_light^2) # [m/s]

	################################
	######SOLAR EPHEMERIS ############
	################################
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)
	PosVector_SolSSB = solar_ephem[1].xyz*1000. #[m]
	VelVector_SolSSB = solar_ephem[2].xyz*1000.0/86400.0  # [m/s]

	################################
	####EQUATORIAL COORD VECTORS ####
	################################
	PosVector_EarthSol = PosVector_EarthSSB - PosVector_SolSSB
	PosMag_EarthSol = norm(PosVector_EarthSol)
	PosHat_EarthSol = PosVector_EarthSol ./ PosMag_EarthSol
	VelVector_EarthSol = (VelVector_EarthSSB - VelVector_SolSSB) / (1.0 + dot(VelVector_SolSSB,VelVector_EarthSSB)/speed_of_light^2)

	#OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	OmegaVector = cross(vec(PosVector_EarthSol), vec(VelVector_EarthSol)) ./ dot(PosVector_EarthSol,PosVector_EarthSol)

	################################
	################################
	# Rotation Axis of the Sun with respect to Equatorial System
	# https://www2.mps.mpg.de/homes/fraenz/systems/systems3art.pdf

	#Eqn 13
	alpha = 286.13 # Rotate along the z axis
	dec = 63.87 # Rotate along the x axis
	Theta = alpha * π/180
	Phi = dec * π/180
	EclipticEpsilon = 23.44
	SolarInclination = 7.25
	SolarLongitude = 75.76

	# Need to perform Extrinsic Euler Rotations
	# Transform to solar rotation axis
	REquatorial = [cos(Theta)*cos(Phi), sin(Theta)*cos(Phi), sin(Phi)]

	OmegaSolVector = [0.0, 0.0,  2.972 *1e-6]
	#OmegaSolHat = [0.0 , 0.0 , 1.0]

	OmegaSolHat = REquatorial
	OmegaSolVector = REquatorial * norm(OmegaSolVector)
	#OmegaSolHat = OmegaSolVector / norm(OmegaSolVector)

	sini = sqrt(1 - dot(OmegaSolHat,PosHat_EarthSol)^2)
	Gamma = 1.04
	DeltaOmega = OmegaSolVector .- OmegaVector

	Delta = ((Gamma* R_sun)^2) * (dot(DeltaOmega,DeltaOmega)*sini^2 - dot(OmegaSolVector,OmegaSolVector))

end

#=
# Still needs work due to pain of Python objects
function CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd::Vector{<:Real}; loc::PyObject )
	# Functions and constants from AstroPy or SciPy that we'll use multiple times
	get_body_barycentric_posvel = AstropyCoordinates.get_body_barycentric_posvel
	get_body_barycentric = AstropyCoordinates.get_body_barycentric
	speed_of_light = Barycorrpy.PhysicalConstants.c
	R_sun = 695700  # km # should we use better, but non-standard value, e.g, 695660?

	JDUTC = AstropyTime.Time(jd , format="jd", scale="utc")
	#return JDUTC
	# Convert times to obtain TDB and TT
	JDTDB = JDUTC.tdb
	JDTT = JDUTC.tt

	################################
	######EARTH EPHEMERIS ############
	################################
	##### NUTATION, PRECESSION, ETC. #####

	# Observatory position wrt Geocenter
	r_pint, v_pint = Barycorrpy.PINT_erfautils.gcrs_posvel_from_itrf(loc, JDUTC, JDTT)
	r_eci = r_pint'  # [m]
	v_eci = v_pint'  # [m/s]

	##### EPHEMERIDES #####
	ephemeris = "de430"
	earth_geo = get_body_barycentric_posvel("earth", JDTDB, ephemeris=ephemeris) # [km]
	r_geo = earth_geo[1].xyz*1000.0 # [m]
	v_geo = earth_geo[2].xyz*1000.0/86400.0  # [m/s]

	PosVector_EarthSSB = r_eci + r_geo # [m]

	# Relativistic Addition of Velocities
	VelVector_EarthSSB = (v_eci + v_geo) / (1.0 + dot(v_eci,v_geo)/speed_of_light^2) # [m/s]

	################################
	######SOLAR EPHEMERIS ############
	################################
	solar_ephem = get_body_barycentric_posvel("sun", JDTDB, ephemeris=ephemeris)
	PosVector_SolSSB = solar_ephem[1].xyz*1000. #[m]
	VelVector_SolSSB = solar_ephem[2].xyz*1000.0/86400.0  # [m/s]

	################################
	####EQUATORIAL COORD VECTORS ####
	################################
	PosVector_EarthSol = PosVector_EarthSSB - PosVector_SolSSB
	PosMag_EarthSol = norm(PosVector_EarthSol)
	PosHat_EarthSol = PosVector_EarthSol ./ PosMag_EarthSol
	VelVector_EarthSol = (VelVector_EarthSSB - VelVector_SolSSB) / (1.0 + dot(VelVector_SolSSB,VelVector_EarthSSB)/speed_of_light^2)

	OmegaVector = map(i->cross(PosVector_EarthSol[:,i], VelVector_EarthSol[:,i]) ./ dot(PosVector_EarthSol[:,i],PosVector_EarthSol[:,i]), 1:size(PosVector_EarthSol,2))
	#OmegaVector = cross(vec(PosVector_EarthSol), vec(VelVector_EarthSol)) ./ dot(PosVector_EarthSol,PosVector_EarthSol)

	################################
	################################
	# Rotation Axis of the Sun with respect to Equatorial System
	# https://www2.mps.mpg.de/homes/fraenz/systems/systems3art.pdf

	#Eqn 13
	alpha = 286.13 # Rotate along the z axis
	dec = 63.87 # Rotate along the x axis
	Theta = alpha * π/180
	Phi = dec * π/180
	EclipticEpsilon = 23.44
	SolarInclination = 7.25
	SolarLongitude = 75.76

	# Need to perform Extrinsic Euler Rotations
	# Transform to solar rotation axis
	REquatorial = [cos(Theta)*cos(Phi), sin(Theta)*cos(Phi), sin(Phi)]
	OmegaSolVector = [0.0, 0.0,  2.972 *1e-6]
	#OmegaSolHat = [0.0 , 0.0 , 1.0]

	OmegaSolHat = REquatorial
	OmegaSolVector = REquatorial * norm(OmegaSolVector)
	#OmegaSolHat = OmegaSolVector / norm(OmegaSolVector)

	println("size(PosHat_EarthSol) = ", size(PosHat_EarthSol))
	println("size(first(PosHat_EarthSol)) = ", size(first(PosHat_EarthSol)))
	sini = map(i->sqrt(1 - dot(OmegaSolHat,PosHat_EarthSol[i])^2), 1:length(PosHat_EarthSol))
	Gamma = 1.04
	DeltaOmega = map(i->OmegaSolVector .- OmegaVector[i],1:length(OmegaVector))
	println("# size(DeltaOmega) = ", size(DeltaOmega))
	println("# size(OmegaSolVector) = ", size(OmegaSolVector))
	println("# size(sini) = ", size(sini))

	Delta = ((Gamma* R_sun)^2) .* (dot.(DeltaOmega,DeltaOmega).*sini^.2 .- dot.(OmegaSolVector,OmegaSolVector))

end
=#

"""
	`CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd; long, lat, alt)`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
	Inputs:
		jd: Julian Date (scalar)
		long: longitude (degrees, west negative)
		lat: latitude (degrees)
		alt: altitude (m)
"""
function CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd::Union{Real,Vector{<:Real}}; long::Real, lat::Real, alt::Real )
#function CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd::Real; long::Real, lat::Real, alt::Real )
	loc = AstropyCoordinates.EarthLocation.from_geodetic(long, lat, height=alt)
	delta = CalculateFWHMDifference_SolarRotation_from_loc_Equatorial(jd, loc=loc)
	#delta = CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic(jd, loc=loc)
	return delta
end

"""
	`CalculateFWHMDifference_SolarRotation_from_obs(jd; obs)`
Calculate the difference between the Observed Solar FWHM and Sidereal Solar FWHM
Inputs:
	jd: Julian Date (scalar)
	obs: observatory name (symbol)
Currently has HARPS-N and NEID/WIYN coordinates.
"""
function CalculateFWHMDifference_SolarRotation_from_obs(jd::Union{Real,Vector{<:Real}}; obs::Symbol )
#function CalculateFWHMDifference_SolarRotation_from_obs(jd::Real; obs::Symbol )
	loc = get_obs_loc(obs)
	CalculateFWHMDifference_SolarRotation_from_long_lat_alt(jd, long=loc["lon"], lat=loc["lat"], alt=loc["elevation"])
end

calc_Δfwhm_solar_rotation(jd::Union{Real,Vector{<:Real}}; obs::Symbol ) = CalculateFWHMDifference_SolarRotation_from_obs(jd, obs=obs)
export calc_Δfwhm_solar_rotation

end # module SolarRotation


text:
using .ImageIO

"""
    function log_image(lg::CLogger, name::AbstractString, obj::AbstractArray{<:Colorant}; step=nothing, kwargs...)

Logs the image from an RGB array.
"""
function log_image(lg::CLogger, name::AbstractString, obj::AbstractArray{<:Colorant}; step=nothing, kwargs...)
    mktempdir(;prefix="jlcom_") do dirname
        fname = joinpath(dirname, name, "image.png")
        FileIO.save(fname, obj)
        log_image(lg, name, fname, step=step, kwargs...)
    end
end

process(lg::CLogger, name::AbstractString, obj::AbstractArray{<:Colorant}, step::Int) = log_image(lg, name, obj; step=step, copy_to_tmp=false)


text:
# LLVM IR generation

function irgen(@nospecialize(job::CompilerJob), method_instance::Core.MethodInstance)
    mod, compiled = @timeit_debug to "emission" compile_method_instance(job, method_instance)
    entry_fn = compiled[method_instance].specfunc
    ctx = context(mod)

    # clean up incompatibilities
    @timeit_debug to "clean-up" begin
        for llvmf in functions(mod)
            # only occurs in debug builds
            delete!(function_attributes(llvmf), EnumAttribute("sspstrong", 0; ctx))

            if Sys.iswindows()
                personality!(llvmf, nothing)
            end

            # remove the non-specialized jfptr functions
            if startswith(LLVM.name(llvmf), "jfptr_")
                unsafe_delete!(mod, llvmf)
            end
        end

        # remove the exception-handling personality function
        if Sys.iswindows() && "__julia_personality" in functions(mod)
            llvmf = functions(mod)["__julia_personality"]
            @compiler_assert isempty(uses(llvmf)) job
            unsafe_delete!(mod, llvmf)
        end
    end

    # target-specific processing
    process_module!(job, mod)
    entry = functions(mod)[entry_fn]

    # sanitize function names
    # FIXME: Julia should do this, but apparently fails (see maleadt/LLVM.jl#201)
    for f in functions(mod)
        LLVM.isintrinsic(f) && continue
        llvmfn = LLVM.name(f)
        startswith(llvmfn, "julia.") && continue # Julia intrinsics
        startswith(llvmfn, "llvm.") && continue # unofficial LLVM intrinsics
        llvmfn′ = safe_name(llvmfn)
        if llvmfn != llvmfn′
            @assert !haskey(functions(mod), llvmfn′)
            LLVM.name!(f, llvmfn′)
        end
    end

    # rename and process the entry point
    if job.source.name !== nothing
        LLVM.name!(entry, safe_name(string("julia_", job.source.name)))
    end
    if job.source.kernel
        LLVM.name!(entry, mangle_call(entry, job.source.tt))
    end
    entry = process_entry!(job, mod, entry)
    compiled[method_instance] =
        (; compiled[method_instance].ci, compiled[method_instance].func,
           specfunc=LLVM.name(entry))

    # minimal required optimization
    @timeit_debug to "rewrite" ModulePassManager() do pm
        global current_job
        current_job = job

        linkage!(entry, LLVM.API.LLVMExternalLinkage)

        # internalize all functions, but keep exported global variables
        exports = String[LLVM.name(entry)]
        for gvar in globals(mod)
            push!(exports, LLVM.name(gvar))
        end
        internalize!(pm, exports)

        # inline llvmcall bodies
        always_inliner!(pm)

        can_throw(job) || add!(pm, ModulePass("LowerThrow", lower_throw!))

        add_lowering_passes!(job, pm)

        run!(pm, mod)
    end

    return mod, compiled
end

## name mangling

# we generate function names that look like C++ functions, because many NVIDIA tools
# support them, e.g., grouping different instantiations of the same kernel together.

function mangle_param(t, substitutions)
    t == Nothing && return "v"

    if isa(t, DataType) || isa(t, Core.Function)
        tn = safe_name(t)

        # handle substitutions
        sub = findfirst(isequal(tn), substitutions)
        if sub === nothing
            str = "$(length(tn))$tn"
            push!(substitutions, tn)
        elseif sub == 1
            str = "S_"
        else
            str = "S$(sub-2)_"
        end

        # encode typevars as template parameters
        if !isempty(t.parameters)
            str *= "I"
            for t in t.parameters
                str *= mangle_param(t, substitutions)
            end
            str *= "E"
        end

        str
    elseif isa(t, Integer)
        "Li$(t)E"
    else
        tn = safe_name(t)
        "$(length(tn))$tn"
    end
end

function mangle_call(f, tt)
    fn = safe_name(f)
    str = "_Z$(length(fn))$fn"

    substitutions = String[]
    for t in tt.parameters
        str *= mangle_param(t, substitutions)
    end

    return str
end

# make names safe for ptxas
safe_name(fn::String) = replace(fn, r"[^A-Za-z0-9_]"=>"_")
safe_name(f::Union{Core.Function,DataType}) = safe_name(String(nameof(f)))
safe_name(f::LLVM.Function) = safe_name(LLVM.name(f))
safe_name(x) = safe_name(repr(x))

## exception handling

# this pass lowers `jl_throw` and friends to GPU-compatible exceptions.
# this isn't strictly necessary, but has a couple of advantages:
# - we can kill off unused exception arguments that otherwise would allocate or invoke
# - we can fake debug information (lacking a stack unwinder)
#
# once we have thorough inference (ie. discarding `@nospecialize` and thus supporting
# exception arguments) and proper debug info to unwind the stack, this pass can go.
function lower_throw!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "lower throw" begin

    throw_functions = [
        # unsupported runtime functions that are used to throw specific exceptions
        "jl_throw"                      => "exception",
        "jl_error"                      => "error",
        "jl_too_few_args"               => "too few arguments exception",
        "jl_too_many_args"              => "too many arguments exception",
        "jl_type_error"                 => "type error",
        "jl_type_error_rt"              => "type error",
        "jl_undefined_var_error"        => "undefined variable error",
        "jl_bounds_error"               => "bounds error",
        "jl_bounds_error_v"             => "bounds error",
        "jl_bounds_error_int"           => "bounds error",
        "jl_bounds_error_tuple_int"     => "bounds error",
        "jl_bounds_error_unboxed_int"   => "bounds error",
        "jl_bounds_error_ints"          => "bounds error",
        "jl_eof_error"                  => "EOF error",
    ]

    for f in functions(mod)
        fn = LLVM.name(f)
        for (throw_fn, name) in throw_functions
            occursin(throw_fn, fn) || continue

            for use in uses(f)
                call = user(use)::LLVM.CallInst

                # replace the throw with a PTX-compatible exception
                let builder = Builder(ctx)
                    position!(builder, call)
                    emit_exception!(builder, name, call)
                    dispose(builder)
                end

                # remove the call
                call_args = arguments(call)
                unsafe_delete!(LLVM.parent(call), call)

                # HACK: kill the exceptions' unused arguments
                #       this is needed for throwing objects with @nospecialize constructors.
                for arg in call_args
                    # peek through casts
                    if isa(arg, LLVM.AddrSpaceCastInst)
                        cast = arg
                        arg = first(operands(cast))
                        isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                    end

                    if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                        unsafe_delete!(LLVM.parent(arg), arg)
                    end
                end

                changed = true
            end

            @compiler_assert isempty(uses(f)) job
            break
         end
     end

    end
    return changed
end

# report an exception in a GPU-compatible manner
#
# the exact behavior depends on the debug level. in all cases, a `trap` will be emitted, On
# debug level 1, the exception name will be printed, and on debug level 2 the individual
# stack frames (as recovered from the LLVM debug information) will be printed as well.
function emit_exception!(builder, name, inst)
    job = current_job::CompilerJob
    bb = position(builder)
    fun = LLVM.parent(bb)
    mod = LLVM.parent(fun)
    ctx = context(mod)

    # report the exception
    if Base.JLOptions().debug_level >= 1
        name = globalstring_ptr!(builder, name, "exception")
        if Base.JLOptions().debug_level == 1
            call!(builder, Runtime.get(:report_exception), [name])
        else
            call!(builder, Runtime.get(:report_exception_name), [name])
        end
    end

    # report each frame
    if Base.JLOptions().debug_level >= 2
        rt = Runtime.get(:report_exception_frame)
        ft = convert(LLVM.FunctionType, rt; ctx)
        bt = backtrace(inst)
        for (i,frame) in enumerate(bt)
            idx = ConstantInt(parameters(ft)[1], i)
            func = globalstring_ptr!(builder, String(frame.func), "di_func")
            file = globalstring_ptr!(builder, String(frame.file), "di_file")
            line = ConstantInt(parameters(ft)[4], frame.line)
            call!(builder, rt, [idx, func, file, line])
        end
    end

    # signal the exception
    call!(builder, Runtime.get(:signal_exception))

    emit_trap!(job, builder, mod, inst)
end

function emit_trap!(@nospecialize(job::CompilerJob), builder, mod, inst)
    ctx = context(mod)
    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    call!(builder, trap)
end

## kernel promotion

@enum ArgumentCC begin
    BITS_VALUE  # bitstype, passed as value
    BITS_REF    # bitstype, passed as pointer
    MUT_REF     # jl_value_t*, or the anonymous equivalent
    GHOST       # not passed
end

function classify_arguments(@nospecialize(job::CompilerJob), codegen_ft::LLVM.FunctionType)
    source_sig = Base.signature_type(job.source.f, job.source.tt)::Type
    source_types = [source_sig.parameters...]

    codegen_types = parameters(codegen_ft)

    args = []
    codegen_i = 1
    for (source_i, source_typ) in enumerate(source_types)
        if isghosttype(source_typ) || Core.Compiler.isconstType(source_typ)
            push!(args, (cc=GHOST, typ=source_typ))
            continue
        end

        codegen_typ = codegen_types[codegen_i]
        if codegen_typ isa LLVM.PointerType && !issized(eltype(codegen_typ))
            push!(args, (cc=MUT_REF, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        elseif codegen_typ isa LLVM.PointerType && issized(eltype(codegen_typ)) &&
               !(source_typ <: Ptr) && !(source_typ <: Core.LLVMPtr)
            push!(args, (cc=BITS_REF, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        else
            push!(args, (cc=BITS_VALUE, typ=source_typ,
                         codegen=(typ=codegen_typ, i=codegen_i)))
        end
        codegen_i += 1
    end

    return args
end

if VERSION >= v"1.7.0-DEV.204"
    function is_immutable_datatype(T::Type)
        isa(T,DataType) && !Base.ismutabletype(T)
    end
else
    function is_immutable_datatype(T::Type)
        isa(T,DataType) && !T.mutable
    end
end

if VERSION >= v"1.7.0-DEV.204"
    function is_inlinealloc(T::Type)
        mayinlinealloc = (T.name.flags >> 2) & 1 == true
        # FIXME: To simple
        return mayinlinealloc
    end
else
    function is_inlinealloc(T::Type)
        return T.isinlinealloc
    end
end

function is_concrete_immutable(T::Type)
    is_immutable_datatype(T) && T.layout !== C_NULL
end

function is_pointerfree(T::Type)
    if !is_immutable_datatype(T)
        return false
    end
    return Base.datatype_pointerfree(T)
end

function deserves_stack(@nospecialize(T))
    if !is_concrete_immutable(T)
        return false
    end
    return is_inlinealloc(T)
end

deserves_argbox(T) = !deserves_stack(T)
deserves_retbox(T) = deserves_argbox(T)
function deserves_sret(T, llvmT)
    @assert isa(T,DataType)
    sizeof(T) > sizeof(Ptr{Cvoid}) && !isa(llvmT, LLVM.FloatingPointType) && !isa(llvmT, LLVM.VectorType)
end

# byval lowering
#
# some back-ends don't support byval, or support it badly, so lower it eagerly ourselves
# https://reviews.llvm.org/D79744
function lower_byval(@nospecialize(job::CompilerJob), mod::LLVM.Module, f::LLVM.Function)
    ctx = context(mod)
    ft = eltype(llvmtype(f))
    @compiler_assert return_type(ft) == LLVM.VoidType(ctx) job

    # find the byval parameters
    byval = BitVector(undef, length(parameters(ft)))
    if LLVM.version() >= v"12"
        for i in 1:length(byval)
            attrs = collect(parameter_attributes(f, i))
            byval[i] = any(attrs) do attr
                kind(attr) == kind(EnumAttribute("byval", 0; ctx))
            end
        end
    else
        # XXX: byval is not round-trippable on LLVM < 12 (see maleadt/LLVM.jl#186)
        #      so we need to re-classify the Julia arguments.
        #      remove this once we only support 1.7.
        has_kernel_state = kernel_state_type(job) !== Nothing
        orig_ft = if has_kernel_state
            # the kernel state has been added here already, so strip the first parameter
            LLVM.FunctionType(return_type(ft), parameters(ft)[2:end]; vararg=isvararg(ft))
        else
            ft
        end
        args = classify_arguments(job, orig_ft)
        filter!(args) do arg
            arg.cc != GHOST
        end
        for arg in args
            if arg.cc == BITS_REF
                # NOTE: +1 since this pass runs after introducing the kernel state
                byval[arg.codegen.i+has_kernel_state] = true
            end
        end
    end

    # fixup metadata
    #
    # Julia emits invariant.load and const TBAA metadta on loads from pointer args,
    # which is invalid now that we have materialized the byval.
    for (i, param) in enumerate(parameters(f))
        if byval[i]
            # collect all uses of the argument
            worklist = Vector{LLVM.Instruction}(user.(collect(uses(param))))
            while !isempty(worklist)
                value = popfirst!(worklist)

                # remove the invariant.load attribute
                md = metadata(value)
                if haskey(md, LLVM.MD_invariant_load)
                    delete!(md, LLVM.MD_invariant_load)
                end
                if haskey(md, LLVM.MD_tbaa)
                    delete!(md, LLVM.MD_tbaa)
                end

                # recurse on the output of some instructions
                if isa(value, LLVM.BitCastInst) ||
                   isa(value, LLVM.GetElementPtrInst) ||
                   isa(value, LLVM.AddrSpaceCastInst)
                    append!(worklist, user.(collect(uses(value))))
                end
            end
        end
    end

    # generate the new function type & definition
    new_types = LLVM.LLVMType[]
    for (i, param) in enumerate(parameters(ft))
        if byval[i]
            push!(new_types, eltype(param::LLVM.PointerType))
        else
            push!(new_types, param)
        end
    end
    new_ft = LLVM.FunctionType(return_type(ft), new_types)
    new_f = LLVM.Function(mod, "", new_ft)
    linkage!(new_f, linkage(f))
    for (arg, new_arg) in zip(parameters(f), parameters(new_f))
        LLVM.name!(new_arg, LLVM.name(arg))
    end

    # emit IR performing the "conversions"
    new_args = LLVM.Value[]
    Builder(ctx) do builder
        entry = BasicBlock(new_f, "entry"; ctx)
        position!(builder, entry)

        # perform argument conversions
        for (i, param) in enumerate(parameters(ft))
            if byval[i]
                # copy the argument value to a stack slot, and reference it.
                ptr = alloca!(builder, eltype(param))
                if LLVM.addrspace(param) != 0
                    ptr = addrspacecast!(builder, ptr, param)
                end
                store!(builder, parameters(new_f)[i], ptr)
                push!(new_args, ptr)
            else
                push!(new_args, parameters(new_f)[i])
                for attr in collect(parameter_attributes(f, i))
                    push!(parameter_attributes(new_f, i), attr)
                end
            end
        end

        # inline the old IR
        value_map = Dict{LLVM.Value, LLVM.Value}(
            param => new_args[i] for (i,param) in enumerate(parameters(f))
        )
        clone_into!(new_f, f; value_map,
                    changes=LLVM.API.LLVMCloneFunctionChangeTypeGlobalChanges)
        # NOTE: we need global changes because LLVM 12 wants to clone debug metadata

        # fall through
        br!(builder, blocks(new_f)[2])
    end

    # update uses of the kernel
    # NOTE: we assume kernel functions can't be called. on-device kernel launches,
    #       e.g. CUDA's dynamic parallelism, will pass the function to an API instead,
    #       and we update those constant expressions arguments here.
    for use in uses(f)
        val = user(use)
        if val isa LLVM.ConstantExpr && opcode(val) == LLVM.API.LLVMPtrToInt
            target = operands(val)[1]
            if target == f
                new_val = LLVM.const_ptrtoint(new_f, llvmtype(val))
                replace_uses!(val, new_val)

                # drop the old constant if it is unused
                # XXX: can we do this differently?
                if isempty(uses(val))
                    LLVM.unsafe_destroy!(val)
                end
            end
        end
    end

    # remove the old function
    # NOTE: if we ever have legitimate uses of the old function, create a shim instead
    fn = LLVM.name(f)
    @assert isempty(uses(f))
    # XXX: there may still be metadata using this function. RAUW updates those,
    #      but asserts on a debug build due to the updated function type.
    unsafe_delete!(mod, f)
    LLVM.name!(new_f, fn)

    return new_f
end

# kernel state arguments
#
# to facilitate passing stateful information to kernels without having to recompile, e.g.,
# the storage location for exception flags, or the location of a I/O buffer, we enable the
# back-end to specify a Julia object that will be passed to the kernel by-value, and to
# every called function by-reference. Access to this object is done using the
# `julia.gpu.state_getter` intrinsic. after optimization, these intrinsics will be lowered
# to refer to the state argument.
#
# note that we deviate from the typical Julia calling convention, by always passing the
# state objects by value instead of by reference, this to ensure that the state object
# is not copied to the stack (because LLVM doesn't see that all uses are read-only).
# in principle, `readonly byval` should be equivalent, but LLVM doesn't realize that.
# also see https://github.com/JuliaGPU/CUDA.jl/pull/1167 and the comments in that PR.
# once LLVM supports this pattern, consider going back to passing the state by reference,
# so that the julia.gpu.state_getter` can be simplified to return an opaque pointer.

# add a state argument to every function in the module, starting from the kernel entry point
function add_kernel_state!(@nospecialize(job::CompilerJob), mod::LLVM.Module,
                           entry::LLVM.Function)
    ctx = context(mod)
    entry_fn = LLVM.name(entry)

    # check if we even need a kernel state argument
    state = kernel_state_type(job)
    @assert job.source.kernel
    if state === Nothing
        return false
    end
    T_state = convert(LLVMType, state; ctx)

    # intrinsic returning an opaque pointer to the kernel state.
    # this is both for extern uses, and to make this transformation a two-step process.
    state_intr = kernel_state_intr(mod, T_state)

    # add a state argument to every function
    worklist = filter(!isdeclaration, collect(functions(mod)))
    workmap = Dict{LLVM.Function, LLVM.Function}()
    for f in worklist
        fn = LLVM.name(f)
        ft = eltype(llvmtype(f))
        LLVM.name!(f, fn * ".stateless")

        # create a new function
        new_param_types = [T_state, parameters(ft)...]
        new_ft = LLVM.FunctionType(return_type(ft), new_param_types)
        new_f = LLVM.Function(mod, fn, new_ft)
        LLVM.name!(parameters(new_f)[1], "state")
        linkage!(new_f, linkage(f))
        for (arg, new_arg) in zip(parameters(f), parameters(new_f)[2:end])
            LLVM.name!(new_arg, LLVM.name(arg))
        end

        workmap[f] = new_f
    end

    # clone and rewrite the function bodies
    for (f, new_f) in workmap
        # use a value mapper for rewriting function arguments
        value_map = Dict{LLVM.Value, LLVM.Value}()
        for (param, new_param) in zip(parameters(f), parameters(new_f)[2:end])
            LLVM.name!(new_param, LLVM.name(param))
            value_map[param] = new_param
        end

        # use a value materializer for replacing uses of the function in constants
        function materializer(val)
            if val isa LLVM.ConstantExpr && opcode(val) == LLVM.API.LLVMPtrToInt
                src = operands(val)[1]
                if haskey(workmap, src)
                    return LLVM.const_ptrtoint(workmap[src], llvmtype(val))
                end
            end
            return val
        end

        # we don't want module-level changes, because otherwise LLVM will clone metadata,
        # resulting in mismatching references between `!dbg` metadata and `dbg` instructions
        clone_into!(new_f, f; value_map, materializer,
                    changes=LLVM.API.LLVMCloneFunctionChangeTypeLocalChangesOnly)

        # we can't remove this function yet, as we might still need to rewrite any called,
        # but remove the IR already
        empty!(f)
    end

    # drop unused constants that may be referring to the old functions
    # XXX: can we do this differently?
    for f in worklist
        for use in uses(f)
            val = user(use)
            if val isa LLVM.ConstantExpr && isempty(uses(val))
                LLVM.unsafe_destroy!(val)
            end
        end
    end

    # update other uses of the old function, modifying call sites to pass the state argument
    function rewrite_uses!(f, new_f)
        # update uses
        Builder(ctx) do builder
            for use in uses(f)
                val = user(use)
                if val isa LLVM.CallInst || val isa LLVM.InvokeInst || val isa LLVM.CallBrInst
                    # NOTE: we unconditionally add the state argument, even if there's no uses,
                    #       assuming we'll perform dead arg elimination during optimization.

                    # forward the state argument
                    position!(builder, val)
                    state = call!(builder, state_intr, Value[], "state")
                    new_val = if val isa LLVM.CallInst
                        call!(builder, new_f, [state, arguments(val)...], operand_bundles(val))
                    else
                        # TODO: invoke and callbr
                        error("Rewrite of $(typeof(val))-based calls is not implemented: $val")
                    end
                    callconv!(new_val, callconv(val))

                    replace_uses!(val, new_val)
                    @assert isempty(uses(val))
                    unsafe_delete!(LLVM.parent(val), val)
                elseif val isa LLVM.ConstantExpr && opcode(val) == LLVM.API.LLVMBitCast
                    # XXX: why isn't this caught by the value materializer above?
                    target = operands(val)[1]
                    @assert target == f
                    new_val = LLVM.const_bitcast(new_f, llvmtype(val))
                    rewrite_uses!(val, new_val)
                    # we can't simply replace this constant expression, as it may be used
                    # as a call, taking arguments (so we need to rewrite it to pass the state)

                    # drop the old constant if it is unused
                    # XXX: can we do this differently?
                    if isempty(uses(val))
                        LLVM.unsafe_destroy!(val)
                    end
                else
                    error("Cannot rewrite unknown use of function: $val")
                end
            end
        end
    end
    for (f, new_f) in workmap
        rewrite_uses!(f, new_f)
        @assert isempty(uses(f))
        unsafe_delete!(mod, f)
    end

    return true
end

# lower calls to the state getter intrinsic. this is a two-step process, so that the state
# argument can be added before optimization, and that optimization can introduce new uses
# before the intrinsic getting lowered late during optimization.
function lower_kernel_state!(fun::LLVM.Function)
    job = current_job::CompilerJob
    mod = LLVM.parent(fun)
    ctx = context(fun)
    changed = false

    # check if we even need a kernel state argument
    state = kernel_state_type(job)
    if state === Nothing
        return false
    end

    # find the kernel state argument. this should be the first argument of the function.
    state_arg = parameters(fun)[1]
    T_state = convert(LLVMType, state; ctx)
    @assert llvmtype(state_arg) == T_state

    # fixup all uses of the state getter to use the newly introduced function state argument
    if haskey(functions(mod), "julia.gpu.state_getter")
        state_intr = functions(mod)["julia.gpu.state_getter"]

        Builder(ctx) do builder
            for use in uses(state_intr)
                inst = user(use)
                @assert inst isa LLVM.CallInst
                bb = LLVM.parent(inst)
                LLVM.parent(bb) == fun || continue

                position!(builder, inst)
                bb = LLVM.parent(inst)
                f = LLVM.parent(bb)

                replace_uses!(inst, state_arg)

                @assert isempty(uses(inst))
                unsafe_delete!(LLVM.parent(inst), inst)

                changed = true
            end
        end
    end

    return changed
end

function cleanup_kernel_state!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false

    # remove the getter intrinsic
    if haskey(functions(mod), "julia.gpu.state_getter")
        intr = functions(mod)["julia.gpu.state_getter"]
        if isempty(uses(intr))
            # if we're not emitting a kernel, we can't resolve the intrinsic to an argument.
            unsafe_delete!(mod, intr)
            changed = true
        end
    end

    return changed
end

function kernel_state_intr(mod::LLVM.Module, T_state)
    ctx = context(mod)

    state_intr = if haskey(functions(mod), "julia.gpu.state_getter")
        functions(mod)["julia.gpu.state_getter"]
    else
        LLVM.Function(mod, "julia.gpu.state_getter", LLVM.FunctionType(T_state))
    end
    push!(function_attributes(state_intr), EnumAttribute("readnone", 0; ctx))

    return state_intr
end

# run-time equivalent
function kernel_state_value(state)
    Context() do ctx
        T_state = convert(LLVMType, state; ctx)

        # create function
        llvm_f, _ = create_function(T_state)
        mod = LLVM.parent(llvm_f)

        # get intrinsic
        state_intr = kernel_state_intr(mod, T_state)

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

            val = call!(builder, state_intr, Value[], "state")

            ret!(builder, val)
        end

        call_function(llvm_f, state)
    end
end


text:
using Documenter, IntervalCensoredEventTimes

makedocs(
    modules = [IntervalCensoredEventTimes],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Hua Zhou",
    sitename = "IntervalCensoredEventTimes.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/Hua-Zhou/IntervalCensoredEventTimes.jl.git",
    push_preview = true
)


text:
@doc raw"""
`PODE`: Partitioned Ordinary Differential Equation

Defines a partitioned initial value problem
```math
\begin{align*}
\dot{q} (t) &= v(t, q(t), p(t)) , &
q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), p(t)) , &
p(t_{0}) &= p_{0} ,
\end{align*}
```
with vector fields ``v`` and ``f``, initial conditions ``(q_{0}, p_{0})`` and the solution
``(q,p)`` taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``.

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``v`` and ``f``
* `v`: function computing the vector field ``v``
* `f`: function computing the vector field ``f``
* `t₀`: initial time
* `q₀`: initial condition for `q`
* `p₀`: initial condition for `p`

The functions `v` and `f` must have the interface
```julia
    function v(t, q, p, v)
        v[1] = ...
        v[2] = ...
        ...
    end
```
and
```julia
    function f(t, q, p, f)
        f[1] = ...
        f[2] = ...
        ...
    end
```
where `t` is the current time, `q` and `p` are the current solution vectors
and `v` and `f` are the vectors which hold the result of evaluating the
vector fields ``v`` and ``f`` on `t`, `q` and `p`.
"""
struct PODE{dType <: Number, tType <: Number,
            vType <: Function, fType <: Function,
            pType <: Union{Tuple,Nothing}, N} <: Equation{dType, tType}

    d::Int
    n::Int
    v::vType
    f::fType
    t₀::tType
    q₀::Array{dType, N}
    p₀::Array{dType, N}
    parameters::pType
    periodicity::Vector{dType}

    function PODE(DT::DataType, N::Int, d::Int, n::Int, v::vType, f::fType,
                  t₀::tType, q₀::DenseArray{dType}, p₀::DenseArray{dType};
                  parameters=nothing, periodicity=zeros(DT,d)) where {
                        dType <: Number, tType <: Number, vType <: Function, fType <: Function}

        @assert d == size(q₀,1) == size(p₀,1)
        @assert n == size(q₀,2) == size(p₀,2)
        @assert ndims(q₀) == ndims(p₀) == N ∈ (1,2)

        new{DT, tType, vType, fType, typeof(parameters), N}(d, n, v, f, t₀,
                convert(Array{DT}, q₀), convert(Array{DT}, p₀),
                parameters, periodicity)
    end
end

function PODE(v, f, t₀, q₀::DenseArray{DT}, p₀::DenseArray{DT}; kwargs...) where {DT}
    PODE(DT, ndims(q₀), size(q₀,1), size(q₀,2), v, f, t₀, q₀, p₀; kwargs...)
end

function PODE(v, f, q₀, p₀; kwargs...)
    PODE(v, f, zero(eltype(q₀)), q₀, p₀; kwargs...)
end

Base.hash(ode::PODE, h::UInt) = hash(ode.d, hash(ode.n, hash(ode.v, hash(ode.f,
        hash(ode.t₀, hash(ode.q₀, hash(ode.p₀, hash(ode.periodicity, hash(ode.parameters, h)))))))))

Base.:(==)(ode1::PODE{DT1,TT1,VT1,FT1}, ode2::PODE{DT2,TT2,VT2,FT2}) where {DT1, DT2, TT1, TT2, VT1, VT2, FT1, FT2} = (
                                ode1.d == ode2.d
                             && ode1.n == ode2.n
                             && ode1.v == ode2.v
                             && ode1.f == ode2.f
                             && ode1.t₀ == ode2.t₀
                             && ode1.q₀ == ode2.q₀
                             && ode1.p₀ == ode2.p₀
                             && ode1.parameters == ode2.parameters
                             && ode1.periodicity == ode2.periodicity)

function Base.similar(ode::PODE, q₀, p₀; kwargs...)
    similar(ode, ode.t₀, q₀, p₀; kwargs...)
end

function Base.similar(ode::PODE, t₀::TT, q₀::DenseArray{DT}, p₀::DenseArray{DT};
                      parameters=ode.parameters, periodicity=ode.periodicity) where {DT  <: Number, TT <: Number}
    @assert ode.d == size(q₀,1) == size(p₀,1)
    PODE(ode.v, ode.f, t₀, q₀, p₀; parameters=parameters, periodicity=periodicity)
end

Base.ndims(ode::PODE) = ode.d


text:
s = """
function f(x)
    local a = 1
    local b = 2
    @lintpragma("Ignore unused a")
    return x+b
end
"""
msgs = lintstr(s)
@test_broken isempty(msgs)

s = """
function f(x)
    local a = 1
    local b = 2
    c = "a"
    @lintpragma("Ignore unused " * c)
    return x+b
end
"""
msgs = lintstr(s)
@test_broken length(msgs) == 2

# @lintpragma can also be used to generate messages
s = """
function f(x)
    local b = 2
    @lintpragma("Info type b")
    return x+b
end
"""
msgs = lintstr(s)
@test length(msgs) == 1
@test msgs[1].code == :I271
@test contains(msgs[1].message, "typeof(b) == Int")

s = """
function f(x)
    local b = 2
    @lintpragma("Warn type b")
    return x+b
end
"""
msgs = lintstr(s)
@test length(msgs) == 1
@test msgs[1].code == :W241
@test contains(msgs[1].message, "typeof(b) == Int")

s = """
function f(x)
    local b = 2
    @lintpragma("Error type b")
    return x+b
end
"""
msgs = lintstr(s)
@test length(msgs) == 1
@test msgs[1].code == :E221
@test contains(msgs[1].message, "typeof(b) == Int")

s = """
function f(x)
    local b = 2
    @lintpragma("Print type b")
    return x+b
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(x)
    local b = 2
    @lintpragma("Print type b[1")
    return x+b
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E138
@test contains(msgs[1].message, "incomplete pragma expression")

s = """
function f(x)
    local b = 2
    @lintpragma("Info me my own reminder")
    return x+b
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test contains(msgs[1].message, "my own reminder")

s = """
function f(x)
    local b = 2
    @lintpragma("Info me " * string(b))
    return x+b
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E137
@test contains(msgs[1].message, "lintpragma must be called using only string literals")

s = """
function f(x)
    local a = 1
    local b = 2
    @lintpragma("Ignore unused a")
    return x+a+b
end
"""
msgs = lintstr(s)
@test_broken msgs[1].code == :I381
@test_broken msgs[1].variable == "Ignore unused a"
@test_broken contains(msgs[1].message, "unused lintpragma")

s = """
@lintpragma()
"""
msgs = lintstr(s)
@test msgs[1].code == :E137
@test contains(msgs[1].message, "lintpragma must be called using only string literals")


text:
#!/usr/bin/env julia
using CSV, DataFrames, PGFPlotsX
it_covid19_pop = CSV.File("data/it_covid19_pop.csv") |> DataFrame
plotdicts = [
    Dict(:num => :CasesFemale, :den => :PopFemale, :lab => "incidens kvinnor"),
    Dict(:num => :CasesMale, :den => :PopMale, :lab => "incidens män"),
    Dict(:num => :DeathsFemale, :den => :PopFemale, :lab => "mortalitet kvinnor"),
    Dict(:num => :DeathsMale, :den => :PopMale, :lab => "mortalitet män")
]
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
p = @pgf Axis({ymode="log", xlabel="Ålder", title="COVID-19 Italien till 2020-03-19",
        legend_pos="outer north east", xmajorgrids, ymajorgrids, yminorgrids})
ages = it_covid19_pop[:AgeGrp10Start]
for pd in plotdicts
    rate = it_covid19_pop[pd[:num]] ./ it_covid19_pop[pd[:den]]
    @pgf push!(p, PlotInc(Table([ages, rate])), LegendEntry(pd[:lab]))
end
pgfsave("../../images/it_covid19.svg", p)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
#  Copyright 2016-2019 Eric S. Tellez <eric.tellez@infotec.mx>
#
#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at
#
#        http://www.apache.org/licenses/LICENSE-2.0
#
#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.

module Graph

using ..SimilaritySearch
using JSON
import ..SimilaritySearch:
    push!, search, fit, KnnResult, optimize!

export LocalSearchAlgorithm, NeighborhoodAlgorithm, SearchGraph,
    compute_aknn, find_neighborhood, push_neighborhood!

abstract type LocalSearchAlgorithm end
abstract type NeighborhoodAlgorithm end

mutable struct SearchGraph{T} <: Index
    db::Vector{T}
    recall::Float64
    k::Int
    links::Vector{Vector{Int32}}
    search_algo::LocalSearchAlgorithm
    neighborhood_algo::NeighborhoodAlgorithm
end

@enum VertexSearchState begin
    UNKNOWN = 0
    VISITED = 1
    EXPLORED = 2
end

function fit(::Type{SearchGraph}, dist::Function, dataset::AbstractVector{T}; recall=0.9, k=10, search_algo=BeamSearch(), neighborhood_algo=LogSatNeighborhood(1.1)) where T
    links = Vector{Int32}[]
    index = SearchGraph(T[], recall, k, links, search_algo, neighborhood_algo)
    for item in dataset
        push!(index, dist, item)
    end

    index
end

include("opt.jl")

## neighborhoods
include("neighborhood/fixedneighborhood.jl")
include("neighborhood/logneighborhood.jl")
include("neighborhood/logsatneighborhood.jl")
include("neighborhood/gallopingneighborhood.jl")
include("neighborhood/essencialneighborhood.jl")
include("neighborhood/satneighborhood.jl")
include("neighborhood/galsatneighborhood.jl")
include("neighborhood/vorneighborhood.jl")

## search algorithms
include("ihc.jl")
include("tihc.jl")
include("beamsearch.jl")

### Basic operations on the index
const OPTIMIZE_LOGBASE = 10

function find_neighborhood(index::SearchGraph{T}, dist::Function, item::T) where T
    n = length(index.db)
    n == 0 && return (NnResult(), Int32[])
    k = ceil(Int, log(OPTIMIZE_LOGBASE, 1+n))
    k1::Int = ceil(Int, log(OPTIMIZE_LOGBASE, 2+n))

    if n > 1 && k != k1
        optimize!(index, dist, recall=index.recall)
    end

    return neighborhood(index.neighborhood_algo, index, dist, item)
end

function push_neighborhood!(index::SearchGraph{T}, item::T, L::AbstractVector{Int32}, n::Int) where T
    for objID in L
        push!(index.links[objID], 1+n)
    end

    push!(index.links, L)
    push!(index.db, item)
end

function push!(index::SearchGraph{T}, dist::Function, item::T) where T
    knn, neighbors = find_neighborhood(index, dist, item)
    push_neighborhood!(index, item, neighbors, length(index.db))
    length(index.db) % 5000 == 0 && @debug "added n=$(length(index.db)), neighborhood=$(length(neighbors)), $(now())"
    knn, neighbors
end

const EMPTY_INT_VECTOR = Int[]

function search(index::SearchGraph{T}, dist::Function, q::T, res::KnnResult; hints::Vector{Int}=EMPTY_INT_VECTOR) where T
    length(index.db) == 0 && return res
    navigation_state = Dict{Int,VertexSearchState}()
    sizehint!(navigation_state, maxlength(res))
    # navigation_state = zeros(UInt8, length(index.db) + 1)
    search(index.search_algo, index, dist, q, res, navigation_state, hints)
end

function optimize!(index::SearchGraph{T}, dist::Function; recall=0.9, k=10, num_queries=128, perf=nothing) where T
    if perf == nothing
        perf = Performance(index.db, dist; expected_k=k, num_queries=num_queries)
    end
    optimize!(index.search_algo, index, dist, recall, perf)
end

## function compute_aknn(index::SearchGraph{T}, dist::Function, k::Int) where T
##     n = length(index.db)
##     aknn = [KnnResult(k) for i=1:n]
##     tabu = falses(length(index.db))
##
##     for i=1:n
##         if i > 1
##             fill!(tabu, false)
##         end
##
##         j = 1
##
##         q = index.db[i]
##         res = aknn[i]
##         for p in res
##             tabu[p.objID] = true
##         end
##
##         function oracle(q::T)
##             # this can be a very costly operation, or can be a very fast estimator, please do some research about it!!
##             if length(res) > 0
##                 a = Iterators.flatten(index.links[p.objID] for p in res)
##                 return Iterators.flatten((a, index.links[i]))
##             else
##                 return index.links[i]
##             end
##         end
##
##         beam_search(index.search_algo, index, dist, q, res, tabu, oracle)
##         for p in res
##             i < p.objID && push!(aknn[p.objID], i, p.dist)
##         end
##
##         if (i % 10000) == 1
##             @debug "algorithm=$(index.search_algo), neighborhood_factor=$(index.neighborhood_algo), k=$(k); advance $i of n=$n"
##         end
##     end
##
##     return aknn
## end

end


text:
@testset "maximumlikelihood.jl" begin
    @testset "fit(model, initialvalues)" begin
        n = 5000

        μ = 0.0
        σ = 1.0
        ξ = 0.1

        ϕ = log(σ)
        θ = [μ; ϕ; ξ]

        pd = GeneralizedExtremeValue(μ, σ, ξ)
        y = rand(pd, n)

        model = Extremes.BlockMaxima(Variable("y", y))

        # Initial value vector length != nparameter throws
        @test_throws AssertionError Extremes.fit(model, [0.0, 0.0, 0.0, 0.0])

        # Initial value vector invalid throws
        @test_throws AssertionError Extremes.fit(model, [Inf, Inf, Inf])

        # Initial value vector valid and length == nparameter does not throw
        @test_logs Extremes.fit(model, [0.0, 0.0, 0.0])

    end

    @testset "fit(model)" begin
        # No solution warn test
        n = 10

        μ = 0.0
        σ = 1.0
        ξ = .8

        ϕ = log(σ)
        θ = [μ; ϕ; ξ]

        pd = GeneralizedExtremeValue(μ, σ, ξ)
        y = rand(pd, n)

        model = Extremes.BlockMaxima(Variable("y", y))

        @test_logs (:warn,"The maximum likelihood algorithm did not find a solution. Maybe try with different initial values or with another method. The returned values are the initial values.") Extremes.fit(model)

        # stationary GEV fit by ML
        n = 5000

        μ = 0.0
        σ = 1.0
        ξ = 0.1

        ϕ = log(σ)
        θ = [μ; ϕ; ξ]

        pd = GeneralizedExtremeValue(μ, σ, ξ)
        y = rand(pd, n)

        model = Extremes.BlockMaxima(Variable("y", y))

        fm = Extremes.fit(model)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

        # non-stationary GEV fit by ML
        n = 5000

        x₁ = randn(n)
        x₂ = randn(n)
        x₃ = randn(n) / 10

        μ = 1.0 .+ x₁
        ϕ = -0.5 .+ x₂
        ξ = x₃

        ϕ = log(σ)
        θ = [1.0; 1.0; -0.5; 1.0; 0.0; 1.0]

        pd = GeneralizedExtremeValue.(μ, σ, ξ)
        y = rand.(pd)

        model = Extremes.BlockMaxima(Variable("y", y), locationcov = [Variable("x₁", x₁)], logscalecov = [Variable("x₂", x₂)], shapecov = [Variable("x₃", x₃)])

        fm = Extremes.fit(model)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

        # stationary GP fit by ML
        n = 5000

        σ = 1.0
        ξ = 0.1

        ϕ = log(σ)
        θ = [ϕ; ξ]

        pd = GeneralizedPareto(σ, ξ)
        y = rand(pd, n)

        model = Extremes.ThresholdExceedance(Variable("y", y))

        fm = Extremes.fit(model)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

        # non-stationary GP fit by ML
        n = 5000

        x₁ = randn(n) / 3
        x₂ = randn(n) / 10

        ϕ = -.5 .+ x₁
        ξ = x₂

        σ = exp.(ϕ)
        θ = [-.5; 1.0; 0.0; 1.0]

        pd = GeneralizedPareto.(σ, ξ)
        y = rand.(pd)

        model = Extremes.ThresholdExceedance(Variable("y", y), logscalecov = [Variable("x₁", x₁)], shapecov = [Variable("x₂", x₂)])

        fm = Extremes.fit(model)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end

    include(joinpath("maximumlikelihood", "maximumlikelihood_gev_test.jl"))
    include(joinpath("maximumlikelihood", "maximumlikelihood_gp_test.jl"))

end


text:
export search, _mk_query_args
"""
    search(q::AbstractQuery; client=get_search_client()) -> Result

Using a query objects, search RediSearch.

# Arguments
- `q::AbstractQuery`: Query object.

# KeyWords
- `client`: Client that will be searched.
"""
function search(q::AbstractQuery; client=get_search_client())
    isnothing(client) && return nothing
    args = _mk_query_args(q; client=client)
    start_time = time()
    res = Jedis.execute([SearchCommands.SEARCH_CMD, args...], client.client)

    return Results(
        res,
        !(q.no_content),
        time() - start_time;
        has_payload=q.with_payloads,
        with_scores=q.with_scores,
    )
end

"""
    _mk_query_args(query; client=get_search_client())

From a query object, generate the search arguments.
"""
function _mk_query_args(query; client=get_search_client())
    args = Vector{Union{Number,String}}([client.index_name])

    append!(args, get_args(query))
    return args
end


text:
# Autogenerated wrapper script for HarfBuzz_ICU_jll for x86_64-linux-gnu
export libharfbuzz_icu

using Cairo_jll
using Fontconfig_jll
using FreeType2_jll
using Glib_jll
using Graphite2_jll
using Libffi_jll
using HarfBuzz_jll
using ICU_jll
JLLWrappers.@generate_wrapper_header("HarfBuzz_ICU")
JLLWrappers.@declare_library_product(libharfbuzz_icu, "libharfbuzz-icu.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Cairo_jll, Fontconfig_jll, FreeType2_jll, Glib_jll, Graphite2_jll, Libffi_jll, HarfBuzz_jll, ICU_jll)
    JLLWrappers.@init_library_product(
        libharfbuzz_icu,
        "lib/libharfbuzz-icu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()