text:
module Acor

using FFTW
using Statistics

""" Returns an estimate of the autocorrelation function of the given
series.  """
function acf(xs::Array{Float64, 1})
    n = size(xs, 1)
    N = 1
    while N < 2*n
        N = N << 1
    end

    mu = mean(xs)

    ys = zeros(N)
    for i in 1:n
        ys[i] = xs[i] - mu
    end

    ys_tilde = rfft(ys)
    ac = irfft(abs2.(ys_tilde), N)

    ac[1:n]/ac[1]
end

""" Returns an estimate of the autocorrelation length of the given
series, or `Inf` if this is not possible to accurately estimate.  """
function acl(xs::Array{Float64, 1})
    ac = acf(xs)

    n = size(ac, 1)

    cumac = 2.0.*cumsum(ac) .- 1.0
    for i in 1:div(n,2)
        if cumac[i] < i/5
            return cumac[i]
        end
    end

    return Inf
end

""" Returns an array giving the estimated ACL for each parameter,
treating `xs` as the output of an ensemble MCMC (shape `(ndim,
nwalkers, nsteps)`).  """
function acl(xs::Array{Float64, 3})
    ndim = size(xs, 1)
    nwalk = size(xs, 2)
    nstep = size(xs, 3)

    acls = zeros(ndim)
    means = zeros(nstep)

    for i in 1:ndim
        for k in 1:nstep
            means[k] = 0.0
            for j in 1:nwalk
                means[k] += xs[i,j,k]
            end
            means[k] /= nwalk
        end
        acls[i] = acl(means)
    end

    acls
end

""" Returns an array of shape `(ndim, nwalkers*nsteps)` where suitable
for thinning to generate independent samples.

This involves both a flattening and a dimension-permutation, so that
successive samples from the same ensemble don't sit adjacent to each
other in memory."""
function flatten(pts)
    p = permutedims(pts, (1, 3, 2))
    reshape(p, (size(p,1), size(p,2)*size(p,3)))
end

""" Thins an ensemble by the maximum autocorrelation length of the
    parameters in the ensemble. """
function thin(pts)
    l = maximum(acl(pts))

    if l == Inf
        warn("Cannot thin because ACL is infinite.")
        zeros(size(pts, 1), size(pts, 2), 0)
    else
        n = convert(Int, ceil(l))
        pts[:,:,1:n:end]
    end
end

""" Returns an estimate of the Gelman-Rubin R statistic for each
parameter, treating each walker as an independent simulation.  The R
statistic estimates the factor by which the variance of the estimated
distribution can be reduced by running the simulation longer."""
function gelman_rubin_rs(ps)
    ndim, nwalkers, nsteps = size(ps)

    rs = zeros(ndim)
    means = zeros(ndim, nwalkers)
    vars = zeros(ndim, nwalkers)

    for i in 1:ndim
        for j in 1:nwalkers
            for k in 1:nsteps
                means[i,j] += ps[i,j,k]
            end
            means[i,j] /= nsteps

            for k in 1:nsteps
                x = ps[i,j,k] - means[i,j]
                vars[i,j] += x*x
            end
            vars[i,j] /= nsteps
        end

        Bon = var(means[i,:])
        W = mean(vars[i,:])

        s2 = (nsteps-1.0)/nsteps*W + Bon

        rs[i] = s2/W
    end

    rs
end

""" An estimate of the AIC (Akaike Information Criterion) from the
MCMC samples of the log-likelihood.  The AIC is `-2.0*max(lnps) +
nparams`.  The estimate assumes Gaussianity of the likelihood
function; under this assumption the mean of the log-likelihood is
`nparams/2.0` below the maximum and the variance of the log-likelihood
is `nparams/2.0`.  So the estimate of the AIC is

    waic = -2.0*mean(lnps) + 2.0*var(lnps)

This estimate is related to (or equal to?) the WAIC described by XXXX
(Gelman ref.).
"""
waic(lnps) = -2.0*mean(lnps) + 2.0*var(lnps)

"""Returns an array of acceptance rates for each walker from the given
chain. """
function acceptance_rate(ps)
    nw = size(ps, 2)
    nt = size(ps, 3)

    ar = zeros(nw)

    for j in 1:nw
        acc = 0
        for k in 1:nt-1
            if ps[:,j,k] == ps[:,j,k+1]
                # Pass
            else
                acc += 1
            end
        end
        ar[j] = float(acc)/float(nt)
    end

    ar
end

""" Returns the acceptance rate given that the chain `ps` has been
thinned by `thin`."""
function acceptance_rate(ps, thin)
    ar = acceptance_rate(ps)

    1 - exp(log1p(-ar)/thin)
end

end


text:
using .PlotAxes
export logrange

function PlotAxes.asplotable(x::AuditorySpectrogram,args...;quantize=(100,128),
    kwds...)
  args = replace(collect(args),:freq => (:freq => logrange))

  asplotable(getcontents(x),args...;quantize=quantize,kwds...)
end

default_quantize(x) = (100,)
default_quantize(x,y) = (100,128,)
default_quantize(x,y,args...) = (100,128,fill(10,length(args))...)

function PlotAxes.asplotable(x::MetaAxisArray;kwds...)
  asplotable(x,axisnames(x)...;kwds...)
end
function PlotAxes.asplotable(x::MetaAxisArray,ax1,axes...;
  quantize=default_quantize(ax1,axes...),kwds...)
  # TODO: handle time or freq not being present
  args = (ax1,axes...)
  args = replace(collect(args),:freq => (:freq => logrange))
  infront = filter(x -> x !== nothing,indexin([:time,:freq => logrange],args))
  others = setdiff(1:length(args),infront)
  args = args[vcat(infront,others)]

  asplotable(getcontents(x),args...;quantize=quantize,kwds...)
end

ext: jl


text:
module NamedVectors

using Printf: @printf
using ..Factories: Inference, TypeFactory, FunctionFactory, addargs!, extendbody!
using ..Factories: MixEscaped, Escaped, UnEscaped
using ..TypeTraits: efficientoperations

export NamedVector, @namedvector
export HomoNamedVector, @homonamedvector

"""
    NamedVector

Abstract type for all named vectors.
"""
abstract type NamedVector end

"""
    getindex(nv::NamedVector, index::Int)

Get the value by the `[]` syntax.
"""
Base.getindex(nv::NamedVector, index::Int) = getfield(nv, index)

"""
    setindex!(nv::NamedVector, value, index::Int)

Set the value by the `[]` syntax if mutable.
"""
Base.setindex!(nv::NamedVector, value, index::Int) = setfield!(nv, index, value)

"""
    ==(nv1::NamedVector, nv2::NamedVector) -> Bool

Overloaded equivalent operator. Two named vector are equal to each other if and only if their keys as well as their values are equal to each other.
!!! note
    It is not necessary for two named vectors to be of the same concrete type to be equal to each other.
"""
Base.:(==)(nv1::NamedVector, nv2::NamedVector) = (keys(nv1) == keys(nv2)) && (values(nv1 ) == values(nv2))

"""
    isequal(nv1::NamedVector, nv2::NamedVector) -> Bool

Overloaded equivalent operator. Two named vector are equal to each other if and only if their keys as well as their values are equal to each other.
!!! note
    It is not necessary for two named vectors to be of the same concrete type to be equal to each other.
"""
Base.isequal(nv1::NamedVector, nv2::NamedVector) = isequal(keys(nv1), keys(nv2)) && isequal(values(nv1), values(nv2))

"""
    <(nv1::NamedVector, nv2::NamedVector) -> Bool

Compare two named vectors and judge whether the first is less than the second.
"""
Base.:<(nv1::NamedVector, nv2::NamedVector) = <(efficientoperations, nv1, nv2)

"""
    isless(nv1::NamedVector, nv2::NamedVector) -> Bool

Compare two named vectors and judge whether the first is less than the second.
"""
Base.isless(nv1::NamedVector, nv2::NamedVector) = isless(efficientoperations, nv1, nv2)

"""
    show(io::IO, nv::NamedVector)

Show a concrete `NamedVector`.
"""
Base.show(io::IO, nv::NamedVector) = @printf io "%s(%s)" nameof(typeof(nv)) join(repr.(values(nv)), ", ")

"""
    hash(nv::NamedVector, h::UInt)

Hash a concrete `NamedVector`.
"""
Base.hash(nv::NamedVector, h::UInt) = hash(values(nv), h)

"""
    convert(::Type{Tuple}, nv::NamedVector) -> Tuple
    convert(::Type{NV}, nv::Tuple) where {NV<:NamedVector} -> NV

Convert a named vector to tuple and vice versa.
"""
@generated Base.convert(::Type{Tuple}, nv::NamedVector) = Expr(:tuple, (:(getfield(nv, $i)) for i = 1:fieldcount(nv))...)
function Base.convert(::Type{NV}, nv::Tuple) where NV<:NamedVector
    @assert fieldcount(NV) == length(nv) "convert error: dismatched length between $NV($(fieldcount(NV))) and input tuple($(length(nv)))."
    return NV(nv...)
end

"""
    length(::Type{NV}) where NV<:NamedVector -> Int
    length(nv::NamedVector) -> Int

Get the length of a concrete `NamedVector`.
"""
Base.length(::Type{NV}) where {NV<:NamedVector} = fieldcount(NV)
Base.length(nv::NamedVector) = length(typeof(nv))

"""
    zero(::Type{NV}) where NV<:NamedVector -> NV
    zero(nv::NamedVector) -> typeof(nv)

Get a concrete `NamedVector` with all values being zero.
"""
@generated function Base.zero(::Type{NV}) where NV<:NamedVector
    zeros = (zero(fieldtype(NV, i)) for i = 1:fieldcount(NV))
    return :(NV($(zeros...)))
end
Base.zero(nv::NamedVector) = zero(typeof(nv))

"""
    iterate(nv::NamedVector, state=1)
    iterate(rv::Iterators.Reverse{<:NamedVector}, state=length(rv.itr))

Iterate or reversely iterate over the values of a concrete `NamedVector`.
"""
Base.iterate(nv::NamedVector, state=1) = (state > length(nv)) ? nothing : (getfield(nv, state), state+1)
Base.iterate(rv::Iterators.Reverse{<:NamedVector}, state=length(rv.itr)) = (state < 1) ? nothing : (getfield(rv.itr, state), state-1)

"""
    keys(nv::NamedVector) -> NTuple{length(nv), Symbol}

Iterate over the names.
"""
Base.keys(nv::NamedVector) = fieldnames(typeof(nv))

"""
    values(nv::NamedVector) -> NTuple{length(nv)}

Iterate over the values.
"""
Base.values(nv::NamedVector) = convert(Tuple, nv)

"""
    pairs(nv::NamedVector) -> Base.Generator

Iterate over the name-value pairs.
"""
Base.pairs(nv::NamedVector) = Base.Generator(=>, keys(nv), values(nv))

"""
    replace(nv::NamedVector; kwargs...) -> typeof(nv)

Return a copy of a concrete `NamedVector` with some of the field values replaced by the keyword arguments.
"""
Base.replace(nv::NamedVector; kwargs...) = replace(efficientoperations, nv; kwargs...)

"""
    map(f, nvs::NV...) where NV<:NamedVector -> NV

Apply function `f` elementwise on the input named vectors.
"""
@generated function Base.map(f, nvs::NV...) where NV<:NamedVector
    exprs = Vector{Expr}(undef, fieldcount(NV))
    for i = 1:length(exprs)
        tmp = [:(nvs[$j][$i]) for j = 1:length(nvs)]
        exprs[i] = :(f($(tmp...)))
    end
    return :(($NV)($(exprs...)))
end

"""
    @namedvector structdef::Expr

Decorate a "raw" struct to be a subtype of `NamedVector`. Here, "raw" means that the input struct has no explicit supertype and no inner constructors.
"""
macro namedvector(structdef::Expr)
    tf = TypeFactory(structdef)
    @assert tf.supertype == Inference(:Any) "@namedvector error: no explicit supertype except `Any` is allowed."
    @assert length(tf.constructors) == 0 "@namedvector error: no inner constructor is allowed."
    tf.supertype = Inference(:NamedVector)
    paramnames = tuple((param.name for param in tf.params)...)
    fieldnames = tuple((field.name for field in tf.fields)...)
    fldnm = FunctionFactory(name = :(Base.fieldnames))
    addargs!(fldnm, :(::Type{<:$(tf.name)}))
    extendbody!(fldnm, Expr(:tuple, QuoteNode.(fieldnames)...))
    structdef = tf(MixEscaped(Escaped(tf.name), UnEscaped(paramnames..., :NamedVector)))
    fldnmdef = fldnm(MixEscaped(Escaped(:tuple)))
    return Expr(:block, :(Base.@__doc__($structdef)), fldnmdef)
end

"""
    HomoNamedVector{T}

Abstract type for all homogeneous named vectors.
"""
abstract type HomoNamedVector{T} <: NamedVector end

"""
    eltype(::Type{<:HomoNamedVector{T}}) where T
    eltype(nv::HomoNamedVector)

Get the type parameter of a concrete `HomoNamedVector`.
"""
Base.eltype(::Type{<:HomoNamedVector{T}}) where T = T
Base.eltype(nv::HomoNamedVector) = eltype(typeof(nv))

"""
    @homonamedvector typename fieldnames dtype::Union{Expr, Symbol}=:nothing mutable::Union{Expr, Bool}=false

Construct a concrete homogeneous named vector with the type name being `typename` and the fieldnames specified by `fieldnames`, and optionally, the type parameters specified by `dtype`.`mutable` can be used as a keyword argument to determine whether the concrete type is mutable.
"""
macro homonamedvector(typename, fieldnames, dtype::Union{Expr, Symbol}=:nothing, mutable::Union{Expr, Bool}=false)
    typename = Symbol(typename)
    fieldnames = tuple(eval(fieldnames)...)
    @assert all(isa(name, Symbol) for name in fieldnames) "homonamedvector error: every field name should be a `Symbol`."
    isa(dtype, Expr) && (dtype.head == :(<:)) && (@assert dtype.args |> length == 1 "homonamedvector error: wrong `dtype`.")
    dname, dscope = (dtype == :nothing) ? (:T, :Any) : (isa(dtype, Expr) && (dtype.head == :(<:))) ? (:T, dtype.args[1]) : (dtype, :concrete)
    if isa(mutable, Expr)
        @assert (mutable.head == :(=)) && (mutable.args[1] == :mutable) && isa(mutable.args[2], Bool) "homonamedvector error: wrong `mutable`."
        mutable = mutable.args[2]
    end
    if dscope == :concrete
        new = :($(esc(typename)))
        super = :(HomoNamedVector{$(esc(dname))})
        body = (:($field::$(esc(dname))) for field in fieldnames)
    else
        new = :($(esc(typename)){$dname<:$(esc(dscope))})
        super = :(HomoNamedVector{$dname})
        body = (:($field::$dname) for field in fieldnames)
    end
    structdef = Expr(:struct, mutable, Expr(:<:, new, super), Expr(:block, body...))
    functions = :(Base.fieldnames(::Type{<:$(esc(typename))}) = $fieldnames)
    return Expr(:block, :(Base.@__doc__($structdef)), functions)
end

end #module


text:
# Default functions

## 1. Variables creation
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String) where T
  return
end

function create_vars!(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String) where T
  return
end

## 2. Power balance
function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T
  return ["", Polynomial(), 0, 0]
end

function Sorig(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T
  return Polynomial()
end

function Sdest(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T
  return Polynomial()
end

## 3. Constraints creation
function constraint(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T
  return Dict{String, Constraint}()
end

function constraint(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}, scenario::String, OPFpbs::OPFProblems) where T
  return Dict{String, Constraint}()
end

## 4. Element cost
function cost(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, Snode::Polynomial, lb, ub) where T
  return 0
end

function cost(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{Link, Variable}, Snode::Polynomial, lb, ub) where T
  return 0
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP, MathOptInterface, DataStructures
import JuMP: value, fix, set_name

"""
    function build_model(;
        instance::UnitCommitmentInstance,
        optimizer = nothing,
        variable_names::Bool = false,
    )::JuMP.Model

Build the JuMP model corresponding to the given unit commitment instance.

Arguments
=========
- `instance`:
    the instance.
- `optimizer`:
    the optimizer factory that should be attached to this model (e.g. Cbc.Optimizer).
    If not provided, no optimizer will be attached.
- `variable_names`:
    If true, set variable and constraint names. Important if the model is going
    to be exported to an MPS file. For large models, this can take significant
    time, so it's disabled by default.
"""
function build_model(;
    instance::UnitCommitmentInstance,
    optimizer = nothing,
    formulation = Formulation(),
    variable_names::Bool = false,
)::JuMP.Model
    @info "Building model..."
    time_model = @elapsed begin
        model = Model()
        if optimizer !== nothing
            set_optimizer(model, optimizer)
        end
        model[:obj] = AffExpr()
        model[:instance] = instance
        _setup_transmission(model, formulation.transmission)
        for l in instance.lines
            _add_transmission_line!(model, l, formulation.transmission)
        end
        for b in instance.buses
            _add_bus!(model, b)
        end
        for g in instance.units
            _add_unit!(model, g, formulation)
        end
        for ps in instance.price_sensitive_loads
            _add_price_sensitive_load!(model, ps)
        end
        _add_system_wide_eqs!(model)
        @objective(model, Min, model[:obj])
    end
    @info @sprintf("Built model in %.2f seconds", time_model)
    if variable_names
        _set_names!(model)
    end
    return model
end


text:
# # Eady instability
#
# This script simulates the instability of a sheared, free-surface
# flow using `ClimateMachine.Ocean.HydrostaticBoussinesqSuperModel`.

using Printf

ENV["GKSwstype"] = "nul"
using Plots
using JLD2
using Plots
using ClimateMachine
using CUDA

ClimateMachine.init()

using ClimateMachine.Ocean
using ClimateMachine.Ocean.Domains
using ClimateMachine.Ocean.Fields

using ClimateMachine.GenericCallbacks: EveryXSimulationTime
using ClimateMachine.GenericCallbacks: EveryXSimulationSteps
using ClimateMachine.Ocean: current_step, Δt, current_time
using CLIMAParameters: AbstractEarthParameterSet, Planet

# # The Eady instability problem
#
# Our objective is to model the instability
# of a sinusoidal geostrophic jet with the streamfunction
#
# ```math
# ψ(y, z) = α λ cos(y / λ) (z + H)
# ```
#
# where ``α`` is geostrophic shear, ``λ`` is the width
# of the sinusoid, and ``H`` is depth. We set ``λ = L / π``,
# where ``L`` is the width of the domain, which fits a
# half wavelength of the sinusoid.
#
# The jet is vertically sheared with the total velocity
# field ``u = - ∂_y ψ``, such that
#
# ```math
# u = α sin(y / λ) (z + H) \, .
# ```
#
# The total buoyancy field is ``b = f ∂_z ψ + N² z``, where
# ``f ∂_z ψ`` is the geostrophic component of buoyancy and
# ``N² z`` reflects a stable stratification with buoyancy
# frequency ``N``. We thus have
#
# ```math
# b = α λ cos(y / λ) + N² z
# ```
#
# The baroclinic component of ``u`` is the "thermal wind"
# associated with ``b``. The barotropic component of ``u``,
# ``⟨u⟩ = \frac{1}{H} ∫ u \rm{d} z``, however, is
# geostrophically-balanced with the sea surface displacement
# via the north-south (``y``) momentum balance,
#
# ```math
# f ⟨u⟩ = - g ∂_y η \, ,
# ```
#
# where ``g`` is gravitational acceleration. Using
# ``⟨u⟩ = α H sin(y / λ) / 2``, we obtain
#
# ```math
# η = α f λ H / (2 g) cos(y / λ) \, .
# ```
#
# We add a periodic barotropic perturbation to
# this unstable state to simulate the growth of
# instability. The periodic perturbation has the
# streamfunction
#
# ```math
# ψ̃(x, y) = \cos(k x) \cos(ℓ y) \exp((y - y_0)^2 / 2L̃^2) \, ,
# ```
#
# where ``k = 12π / L``, ``ℓ = 6π / L``, ``y_0 = L/2``, and
# ``L̃ = 5 \times 10^4`` m. The velocity field and height
# perturbation associated with this streamfunction is
#
# ```math
# \begin{gather}
# ũ = - ∂_y ψ̃ = (ℓ \tan(ℓ y) + (y - y₀) / L̃^2 ) ψ̃
# ṽ = + ∂_x ψ̃ = k \tan(k x) ψ̃
# η = + f ψ̃  / g
# \end{gather}
#```
#
# We set the amplitude of the perturbation to 1% of the
# basic state.
#
# # Parameters
#
# We now choose parameters for the Eady instability problem.

Nh = 32
Nz = 4
Np = 4

L = 1e6 # Domain width (m)
H = 1e3 # Domain height (m)
f = 1e-4 # Coriolis parameter (s⁻¹)
α = 5f # Geostrophic shear (s⁻¹)
N² = 1e-5 # Initial buoyancy gradient (s⁻²)
νh = κh = 1e3 # Horizontal viscosity and diffusivity (m² s⁻¹)
νz = κh = 1e-2 # Vertical viscosity and diffusivity (m² s⁻¹)

minutes = 60.0
hours = 60minutes
days = 24hours
years = 365days

stop_time = 60days # Simulation stop time

struct EarthParameters <: AbstractEarthParameterSet end
Planet.grav(::EarthParameters) = 0.1
g = Planet.grav(EarthParameters())

# # The domain

domain = RectangularDomain(
    Ne = (Int(Nh / Np), Int(Nh / Np), Int(Nz / Np)),
    Np = Np,
    x = (0, L),
    y = (0, L),
    z = (-H, 0),
    periodicity = (true, false, false),
)

# # Boundary conditions

no_slip = OceanBC(Impenetrable(NoSlip()), Insulating())
free_surface = OceanBC(Penetrable(FreeSlip()), Insulating())

# # Initial conditions
#
# We idealize θ as buoyancy. We initialize the flow with a
# vertically-sheared, horizontally-sinusoidal jet in geostrophic
# balance. We add a small amount of noise to the buoyancy field.

## Geostrophic jet:
λ = L / π # Sinusoidal jet width (m)

U(y, z) = +α * sin(y / λ) * (z + H)
B(y, z) = +α * f * λ * cos(y / λ) + N² * z
h(y, z) = α * f * λ * H / 2g * cos(y / λ)

# Geostrophic mode-3 perturbation:
width = 5e4
kx = 12π / L
ky = 6π / L
y₀ = L / 2

ψ̃(x, y) = cos(kx * x) * cos(ky * y) * exp(-(y - y₀)^2 / 2width^2)

# ũ = - ∂y ψ̃
ũ(x, y) = (ky * tan(ky * y) + (y - y₀)^2 / width^2) * ψ̃(x, y)

# ṽ = ∂x ψ̃
ṽ(x, y) = -kx * tan(kx * x) * ψ̃(x, y)

a = 0.01
uᵢ(x, y, z) = U(y, z) + a * ũ(x, y)
vᵢ(x, y, z) = a * ṽ(x, y)
ηᵢ(x, y, z) = h(y, z) + a * f / g * ψ̃(x, y)
θᵢ(x, y, z) = B(y, z)

initial_conditions = InitialConditions(u = uᵢ, v = vᵢ, θ = θᵢ, η = ηᵢ)

minutes = 60.0

model = Ocean.HydrostaticBoussinesqSuperModel(
    domain = domain,
    time_step = 5minutes,
    initial_conditions = initial_conditions,
    parameters = EarthParameters(),
    # array_type = CuArray,
    buoyancy = (αᵀ = 1 / g,),
    coriolis = (f₀ = f, β = 0),
    turbulence_closure = (νʰ = νh, κʰ = κh, νᶻ = νz, κᶻ = νz),
    rusanov_wave_speeds = (cʰ = sqrt(g * H), cᶻ = 1e-2),
    boundary_tags = ((0, 0), (1, 1), (1, 2)),
    boundary_conditions = (no_slip, free_surface),
)

# We prepare a callback that periodically fetches the horizontal velocity and
# tracer concentration for later animation,

fetched_states = []

realdata = Array(model.state.realdata)
u = SpectralElementField(domain, model.grid, view(realdata, :, 1, :))

volume = assemble(u)
x = volume.x[:, 1, 1]
y = volume.y[1, :, 1]

start_time = time_ns()

data_fetcher = EveryXSimulationTime(days) do

    realdata = Array(model.state.realdata)
    u = SpectralElementField(domain, model.grid, view(realdata, :, 1, :))
    v = SpectralElementField(domain, model.grid, view(realdata, :, 2, :))
    η = SpectralElementField(domain, model.grid, view(realdata, :, 3, :))
    θ = SpectralElementField(domain, model.grid, view(realdata, :, 4, :))

    umax = maximum(abs, u)
    elapsed = (time_ns() - start_time) * 1e-9

    step = @sprintf("Step: %d", current_step(model))
    sim_time = @sprintf("time: %.2f days", current_time(model) / days)
    wall_time = @sprintf("wall time: %.2f min", elapsed / 60)

    @info "$step, $sim_time, $wall_time"

    isnan(umax) && error("NaN'd out.")

    u_assembly = assemble(data.(u.elements))
    v_assembly = assemble(data.(v.elements))
    θ_assembly = assemble(data.(θ.elements))

    push!(
        fetched_states,
        (
            u = u_assembly,
            v = v_assembly,
            θ = θ_assembly,
            time = current_time(model),
        ),
    )
end

# and then run the simulation.

model.solver_configuration.timeend = stop_time

try
    result = ClimateMachine.invoke!(
        model.solver_configuration;
        user_callbacks = [data_fetcher],
    )
catch err
    @warn "Simulation stopped early because " print(showerror, err)
end

# Finally, we make an animation of the evolving shear instability.

animation = @animate for (i, state) in enumerate(fetched_states)
    local u
    local v
    local θ

    @info "Plotting frame $i of $(length(fetched_states))..."

    kwargs = (
        xlim = domain.x,
        ylim = domain.y,
        linewidth = 0,
        aspectratio = 1,
        xlabel = "x (m)",
        ylabel = "y (m)",
    )

    u = state.u[:, :, end]
    v = state.v[:, :, end]
    θ = state.θ[:, :, end]

    ulim = 1.2 * α * H
    θmin = minimum(θ)
    θmax = maximum(θ)

    ulevels = range(-ulim, ulim, length = 31)
    θlevels = range(θmin, θmax, length = 31)

    u_plot = heatmap(
        x,
        y,
        clamp.(u, -ulim, ulim)';
        #levels = ulevels,
        color = :balance,
        clim = (-ulim, ulim),
        kwargs...,
    )

    v_plot = contourf(
        x,
        y,
        clamp.(v, -ulim, ulim)';
        #levels = ulevels,
        color = :balance,
        clim = (-ulim, ulim),
        kwargs...,
    )

    θ_plot = heatmap(
        x,
        y,
        clamp.(θ, θmin, θmax)';
        #levels = θlevels,
        color = :thermal,
        clim = (θmin, θmax),
        kwargs...,
    )

    u_title =
        @sprintf("surface u (m s⁻¹) at t = %.2f days", state.time / days)
    v_title =
        @sprintf("surface v (m s⁻¹) at t = %.2f days", state.time / days)
    θ_title =
        @sprintf("surface b (m s⁻²) at t = %.2f days", state.time / days)

    plot(
        u_plot,
        v_plot,
        θ_plot,
        title = [u_title v_title θ_title],
        layout = (1, 3),
        size = (1500, 400),
    )
end

name =
    @sprintf("eady_dg_Np%d_Nh%d_Nz%d_αf%.2e_ν%.2e.gif", Np, Nh, Nz, α / f, νh)

gif(animation, name, fps = 8)


text:
const USE_GPU = false  # Use GPU? If this is set false, then no GPU needs to be available
const GPU_ID  = 0
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    CUDA.device!(GPU_ID) # select GPU
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using MAT, Random, Plots, Printf, Statistics
##################################################
@parallel_indices (ix,iy) function compute_1!(Yf::Data.Array, v1::Data.Number, v2::Data.Number, a::Data.Number, b::Data.Number, dx::Data.Number, dy::Data.Number)
    if (ix<=size(Yf,1) && iy<=size(Yf,2))  Yf[ix,iy] = Yf[ix,iy] + a*sin( dx*(ix-0.5)*v1 + dy*(iy-0.5)*v2 ) + b*cos( dx*(ix-0.5)*v1 + dy*(iy-0.5)*v2 )  end
    return
end

@parallel function compute_2!(Yf::Data.Array, c::Data.Number)
    @all(Yf) = @all(Yf)*c
    return
end
##################################################
@views function RandomField2De()
    do_viz  = false
    do_save = false
    Random.seed!(1234)     # Resetting the random seed
    # Physics
    lx, ly = 100.0, 100.0  # domain size
    sf     = 1.0           # standard deviation
    cl     = (10.0, 8.0)   # correlation lengths in [x,y,z]
    # Numerics
    nh     = 10000         # inner parameter, number of harmonics
    nx, ny = 128, 128      # numerical grid resolution
    # Derived numerics
    dx, dy = lx/nx, ly/ny  # numerical grid step size
    c      = sf/sqrt(nh)
    # Array allocation
    Yf     = @zeros(nx, ny)
    # Scalar allocations
    ϕ=0.0; k=0.0; d=0.0; θ=0.0; v1=0.0; v2=0.0; a=0.0; b=0.0
    # Preparation of visualisation
    if do_viz
        X, Y = -lx/2:dx:lx/2, -ly/2:dy:ly/2
    end
    println("Starting 2D RandomField generation (anisotropic exponential covariance function)...")
    # Loop over nh harmonics
    for ih = 1:nh
        if (ih==501)  global wtime0 = Base.time()  end
        ϕ    = 2.0*pi*rand()
        # Gaussian spectrum
        flag = true
        while flag
            k  = tan(pi*0.5*rand())
            d  = (k*k)/(1.0+(k*k))
            if (rand()<d)  flag = false  end
        end
        θ    = acos(1.0-2.0*rand())
        v1   = k*sin(ϕ)*sin(θ)/cl[1]
        v2   = k*cos(ϕ)*sin(θ)/cl[2]
        a, b = randn(), randn()
        @parallel compute_1!(Yf, v1, v2, a, b, dx, dy)
    end
    @parallel compute_2!(Yf, c)
    # Performance
    wtime    = Base.time() - wtime0
    A_eff    = 2/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    wtime_it = wtime/(nh-500)                   # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                   # Effective memory throughput [GB/s]
    @printf("Total harmonic iters=%d, time=%1.3e sec (@ T_eff = %1.2f GB/s) \n", nh, wtime, round(T_eff, sigdigits=2))
    # Visualisation
    if do_viz
        display(heatmap(X, Y, Array(Yf)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:hot, title="2D RandomField"))
    end
    if do_save  file = matopen("grf2D_expon.mat", "w"); write(file, "grf2D", Array(Yf)); close(file)  end
    return
end

RandomField2De()


text:
# Common utilities

##### common types

"""
    ClusteringResult

Base type for the output of clustering algorithm.
"""
abstract type ClusteringResult end

# generic functions

"""
    nclusters(R::ClusteringResult) -> Int

Get the number of clusters.
"""
nclusters(R::ClusteringResult) = length(R.counts)

"""
    counts(R::ClusteringResult) -> Vector{Int}

Get the vector of cluster sizes.

`counts(R)[k]` is the number of points assigned to the ``k``-th cluster.
"""
counts(R::ClusteringResult) = R.counts

"""
    wcounts(R::ClusteringResult) -> Vector{Float64}
    wcounts(R::FuzzyCMeansResult) -> Vector{Float64}

Get the weighted cluster sizes as the sum of weights of points assigned to each
cluster.

For non-weighted clusterings assumes the weight of every data point is 1.0,
so the result is equivalent to `convert(Vector{Float64}, counts(R))`.
"""
wcounts(R::ClusteringResult) = convert(Vector{Float64}, counts(R))

"""
    assignments(R::ClusteringResult) -> Vector{Int}

Get the vector of cluster indices for each point.

`assignments(R)[i]` is the index of the cluster to which the ``i``-th point
is assigned.
"""
assignments(R::ClusteringResult) = R.assignments

##### convert display symbol to disp level
const DisplayLevels = Dict(:none => 0, :final => 1, :iter => 2)

display_level(s::Symbol) = get(DisplayLevels, s) do
    throw(ArgumentError("Invalid value for the 'display' option: $s."))
end

##### update minimum value

function updatemin!(r::AbstractArray, x::AbstractArray)
    n = length(r)
    length(x) == n || throw(DimensionMismatch("Inconsistent array lengths."))
    @inbounds for i = 1:n
        xi = x[i]
        if xi < r[i]
            r[i] = xi
        end
    end
    return r
end


text:
# Types that are taken from conll file and used in parser
typealias Word String           # stands for words in sentence
typealias PosTag  UInt8         # 17 universal part-of-speech tags
typealias DepRel UInt8          # 37 universal depedency relations
typealias Position Int16        # sentence position
typealias Cost Position         # [0:nword]
typealias Move Int              # [1:nmove]
typealias WordId Int32          # [1:nvocab]
typealias Pvec Vector{Position} # used for stack, head ? don't sure what it is
typealias Dvec Vector{DepRel}   # used for deprel

typealias XPosTag UInt8		# language specific pos-tags

# Universal POS tags (17)
const UPOSTAG = Dict{String,PosTag}(
"ADJ"   => 1, # adjective
"ADP"   => 2, # adposition
"ADV"   => 3, # adverb
"AUX"   => 4, # auxiliary
"CCONJ" => 5, # coordinating conjunction
"DET"   => 6, # determiner
"INTJ"  => 7, # interjection
"NOUN"  => 8, # noun
"NUM"   => 9, # numeral
"PART"  => 10, # particle
"PRON"  => 11, # pronoun
"PROPN" => 12, # proper noun
"PUNCT" => 13, # punctuation
"SCONJ" => 14, # subordinating conjunction
"SYM"   => 15, # symbol
"VERB"  => 16, # verb
"X"     => 17, # other
)

# Universal Dependency Relations (37)
const UDEPREL = Dict{String,DepRel}(
"root"       => 1,  # root
"acl"        => 2,  # clausal modifier of noun (adjectival clause)
"advcl"      => 3,  # adverbial clause modifier
"advmod"     => 4,  # adverbial modifier
"amod"       => 5,  # adjectival modifier
"appos"      => 6,  # appositional modifier
"aux"        => 7,  # auxiliary
"case"       => 8,  # case marking
"cc"         => 9,  # coordinating conjunction
"ccomp"      => 10, # clausal complement
"clf"        => 11, # classifier
"compound"   => 12, # compound
"conj"       => 13, # conjunct
"cop"        => 14, # copula
"csubj"      => 15, # clausal subject
"dep"        => 16, # unspecified dependency
"det"        => 17, # determiner
"discourse"  => 18, # discourse element
"dislocated" => 19, # dislocated elements
"expl"       => 20, # expletive
"fixed"      => 21, # fixed multiword expression
"flat"       => 22, # flat multiword expression
"goeswith"   => 23, # goes with
"iobj"       => 24, # indirect object
"list"       => 25, # list
"mark"       => 26, # marker
"nmod"       => 27, # nominal modifier
"nsubj"      => 28, # nominal subject
"nummod"     => 29, # numeric modifier
"obj"        => 30, # object
"obl"        => 31, # oblique nominal
"orphan"     => 32, # orphan
"parataxis"  => 33, # parataxis
"punct"      => 34, # punctuation
"reparandum" => 35, # overridden disfluency
"vocative"   => 36, # vocative
"xcomp"      => 37, # open clausal complement
)

const XPOSTAG = Dict{String, XPosTag}()	# dictionary for xpostags

immutable Vocab
    cdict::Dict{Char, Int}         # character vocabulary
    idict::Dict{String, Int}       # word dictionary (input) obtained from conll file
    odict::Dict{String, Int}       # word dictionary (output) obtained from lm training
    sosword::String                # start of sentence word (<s>)
    eosword::String                # end of sentence word (</s>)
    unkword::String                # unknown word (<unk>)
    sowchar::Char                  # sow character chosen from hardware character
    eowchar::Char                  # eow character chosen from hardware character
    unkchar::Char
    postags::Dict{String, PosTag}
    deprels::Dict{String, DepRel}

    xpostags::Dict{String, XPosTag}
end

                                # CONLLU FORMAT
type Sentence                   # 1. ID: Word index, integer starting at 1 for each new sentence
    word::Vector{Word}          # 2. FORM: Word form or punctuation symbol
    # stem::Vector{Stem}        # 3. LEMMA: Lemma for m or punctiation symbol
    postag::Vector{PosTag}      # 4. UPOSTAG: Universal part-of-speech tag
    xpostag::Vector{XPosTag}  # 5. Language specific pos-tag
    # feats::Vector{Feats}      # 6. FEATS: List of morphological features from the universal feature inventory or from a defined language-specific extension; underscore if not available
    head::Vector{Position}      # 7. HEAD: Head of the current word, which is either a vvalue of ID or zero
    deprel::Vector{DepRel}       # 8. DEPREL: Universal dependency relation to the HEAD(root iff HEAD=0) or a defined language-specific subtype
    # deps::Vector{Deps}        # 9. DEPS: Enhanced dependency graph in the form of a list of head-deprel pairs.
    # misc::Vector{Misc}        # 10. MISC: Any other annotation.

    # language model dependent features
    wvec::Vector                # word vectors
    fvec::Vector                # forw context vectors
    bvec::Vector                # backw context vectors
    vocab::Vocab                # go to 13
    parse

    Sentence(v::Vocab) = new([],[],[], [], [], [], [], [], v, nothing)
end

Base.length(s::Sentence) = length(s.word)

# add-hoc solution for parser.jl not sure whether it is need?
typealias Corpus AbstractVector{Sentence}


text:
tests = [
         "univariate",
         "bivariate",
         "interp",
         ]

println("Running tests:")

for t in tests
    test_fn = "$t.jl"
    println(" * $test_fn")
    include(test_fn)
end


text:
using SparseArrays, CUDA, LinearAlgebra

CUDA.allowscalar(false)

N = 400
uA2 = cu(zeros(N,N))
A2u = cu(zeros(N,N))
A2 = cu(sprand(N, N, 0.01))
u = cu(rand(N, N))

mul!(A2u, A2, u)
mul!(uA2, u, A2')


text:
abstract type K80 <: NASM end

function P(mod::K80, t::Array{Float64})
  return [P(mod, i) for i in t]
end


text:
# Exports
export get_abs_time, get_rel_time

"""
Computes the absolute time
"""
function get_abs_time(problem::SatPlanningProblem, t_since::Real)
    return problem.t_start + t_since
end

"""
Computes the relative time
"""
function get_rel_time(problem::SatPlanningProblem, epc::Epoch)
    return epc - problem.t_start
end


text:
"""
Implements a pub-sub model, where you can subscribe to a feed with an optional set of filters
and get only those callbacks where all your filterset overlaps with the matching filterset
of the publisher.  Once registered, filters are fixed.  You must unregister and then register
again to update filters.

# this will set up and register a publisher which matches subscribers listening for UID(1) and EDGX
# and also subscribers that don't filter on those filterkeys
pub = Publisher(somefunction, Filters([:uid, UID(1)], [:exch, EDGX]))

# this subscribes the receiving object on all callbacks to somefunction.  the callbacks will occur
# whenever there is a publish command on somefunction with uid filter that includes either UID(1) or UID(2)
sub = subscribe(somefunction, recvObject, Filters([:uid, UID(1), UID(2)]))

# this generates all the callbacks for connected listeners, effectively calling:
#     somefunction(recvObject, args...)
# for each connected subscriber
publish(pub, args...)

# done listening...
unregister(sub)

# done publishing
unregister(pub)

"""

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

typealias FilterSet Set{Any}
typealias FilterMap Dict{Symbol,FilterSet}

immutable Filters
  d::FilterMap
end
Filters() = Filters(FilterMap())
Filters(filterlists::Vector...) = Filters(FilterMap(map(x -> (Symbol(x[1]), FilterSet(x[2:end])), filterlists)))

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

# a Publisher is set up from the broadcast side
type Publisher{F<:Function}
  f::F
  listeners::Vector{Any}
  filters::Filters
end
Publisher(f::Function, filters::Filters = Filters()) = register(Publisher(f, Any[], filters))

Base.isempty(publisher::Publisher) = isempty(publisher.listeners)

# call this to actually trigger the callbacks
function publish(publisher::Publisher, args...)
    for listener in publisher.listeners
        publisher.f(listener, args...)
    end
end

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

immutable Subscriber{F<:Function,L}
  f::F
  listener::L
  filters::Filters
end

# call this to start listening
function subscribe(f::Function, listener, filters::Filters = Filters())
  register(Subscriber(f, listener, filters))
end

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

"""
A wrapper to schedule method calls for the future.  For example, a minute timer:

```
    cb = Callback(calculate, my_obj)
    schedule(TimeOfDay("9:31") : TimeOfDay("00:01") : stopTime, cb)
```

This example will call `calculate(my_obj)` once per minute in the simulation.
"""
immutable Callback
    pub::Publisher
    subs::Vector{Subscriber}
end

function Callback(f::Function, objs...)
    filters = Filters(vcat(:obj, objs))
    subs = [subscribe(f, obj, filters) for obj in objs]
    pub = Publisher(f, filters)
    Callback(pub, subs)
end

# function Callback(f::Function, args...)
#     obj = args[1]
#     filters = Filters(vcat(:obj, obj))
#     subs = [subscribe(f, obj, filters) for obj in args]
#     pub = Publisher(f, filters)
#     Callback(pub, subs)
# end

publish(cb::Callback, args...) = publish(cb.pub, args...)

# # subscribe ourselves
# filter = Filters([:obj, top])
# subscribe(calculate, top, filter)

# # set up a timer with the same callback function and filter (I'll wrap this functionality eventually)
# timer = Publisher(calculate, filter)
# schedule(TimeOfDay("9:31") : TimeOfDay("00:01") : stopTime, timer)

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

immutable Hub
  subscribers::Set{Subscriber}
  publishers::Set{Publisher}
end

const HUB = Hub(Set{Subscriber}(), Set{Publisher}())

function reset_hub()
  for publisher in HUB.publishers
    empty!(publisher.listeners)
  end
  empty!(HUB.publishers)
  empty!(HUB.subscribers)
end

function register(subscriber::Subscriber)
  # add to subscribers list
  push!(HUB.subscribers, subscriber)

  # match subscriber filters to publisher's filters... add this subscriber's listener to matching publishers
  for publisher in HUB.publishers
    if matches(publisher, subscriber)
      push!(publisher.listeners, subscriber.listener)
    end
  end

  subscriber
end

function unregister(subscriber::Subscriber)
  delete!(HUB.subscribers, subscriber)

  for publisher in HUB.publishers

    # TODO: this should be a simple "delete!" call, but doesn't work for vectors
    delidx = 0
    for (i,listener) in enumerate(publisher.listeners)
      if listener === subscriber.listener
        delidx = i
        break
      end
    end
    if delidx > 0
      deleteat!(publisher.listeners, delidx)
    end
  end
end

function register(publisher::Publisher)
  # add to publishers list
  push!(HUB.publishers, publisher)

  # match publisher's filters to subscriber's filters... add matching subscribers' anonfuns to publisher.
  for subscriber in HUB.subscribers
    if matches(publisher, subscriber)
      push!(publisher.listeners, subscriber.listener)
    end
  end

  publisher
end

function unregister(publisher::Publisher)
  delete!(HUB.publishers, publisher)
end

function matches(publisher::Publisher, subscriber::Subscriber)
  # return instantly if it's not the right function
  publisher.f == subscriber.f || return false

  # now lets evaluate the publisher's filters and, for each filter symbol, ensure the subscriber:
  #   1) doesn't have a filter for that symbol, or
  #   2) has a non-empty intersection between the filter sets for that symbol
  # if both are false, return false
  for (filterkey, pubset) in publisher.filters.d
    if haskey(subscriber.filters.d, filterkey)
      # both sub and pub have this filter key... if the values don't overlap then there's no match
      subscriberset = subscriber.filters.d[filterkey]
      isempty(intersect(pubset, subscriberset)) && return false
    end
  end

  # if we got this far, return true
  return true
end


text:
type TokenData
  id::Int32
  freq::Int32
end

TokenData() = new(0, 0)

type TfidfModel
  voc::Dict{String,TokenData}
  size::Int64
  filter_low::Int
  filter_high::Float64
  config::TextConfig
end

function save(ostream, model::TfidfModel)
    M = Dict{String,Any}(
        "length" => length(model.voc),
        "size" => model.size,
        "filter_low" => model.filter_low,
        "filter_high" => model.filter_high,
    )
    write(ostream, JSON.json(M), "\n")
    save(ostream, model.config)
    for (k, v) in model.voc
        write(ostream, JSON.json((k, v.id, v.freq)), "\n")
    end
end

function load(istream, ::Type{TfidfModel})
    m = JSON.parse(readline(istream))
    config = load(istream, TextConfig)
    voc = Dict{String,TokenData}()
    for i in 1:m["length"]
        k, id, freq = JSON.parse(readline(istream))
        voc[k] = TokenData(id, freq)
    end
    TfidfModel(voc, m["size"], m["filter_low"], m["filter_high"], config)
end

TfidfModel() = TfidfModel(Dict{String,TokenData}(), 0, 1, 0.70, TextConfig())

function TfidfModel(corpus, config::TextConfig)
  model = TfidfModel()
  model.config = config
  V = Dict{String,Int}()
  for text in corpus
    voc, maxfreq = compute_vocabulary(text, model.config)
    for (token, freq) in voc
      V[token] = get(V, token, 0) + freq
    end
    model.size += 1
  end

  for (token, freq) in V
    if freq < model.filter_low || freq > model.filter_high * model.size
      continue
    end
    model.voc[token] = TokenData(length(model.voc)+1, freq)
  end

  model
end

function compute_vocabulary(text::String, config::TextConfig)
    voc = Dict{String,Int}()
    maxfreq = 0
    for token in expand_tokens(text, config)
        freq = get(voc, token, 0) + 1
        voc[token] = freq
        maxfreq = max(maxfreq, freq)
    end

    voc, maxfreq
end

function vectorize(text::String, model::TfidfModel; maxlength=typemax(Int))
  localvoc = Dict{String,Int32}()
  vec = Term[]
  localvoc, maxfreq = compute_vocabulary(text, model.config)
  for (token, freq) in localvoc
    tokendata = try
      model.voc[token]
    catch err
      continue
    end

    tf = freq / maxfreq
    idf = log(model.size / tokendata.freq)
    push!(vec, Term(tokendata.id, tf * idf))
    # push!(vec, (token, tfidf))
  end

  if length(vec) > maxlength
     sort!(vec, by=(x) -> -x.weight)
     vec = vec[1:maxlength]
  end
  sort!(vec, by=(x) -> x.id)

  DocumentType(vec)
end

max_stars_repo_name: INGEOTEC/ranix

max_issues_repo_name: INGEOTEC/ranix

max_forks_repo_name: INGEOTEC/ranix


text:
using Markdown
using InteractiveUtils
using CSV, DataFrames, Query
using Dates, MD5
using ArgParse

  function parse_commandline_verify_downloads()
     s = ArgParseSettings( description = "Verify FITS file downloads match contents of meta.csv.")
     @add_arg_table! s begin
         "inputdir"
            help = "Directory with downloaded FITS files."
            arg_type = String
            required = true
#=
         "--root"
             help = "Path to root of data directories"
             arg_type = String
             default = "/gpfs/group/ebf11/default/RISE_NEID_ebf/data"
         "--input"
             help = "Path to inputs FITS files (fits)"
             arg_type = String
             default = "solar_L2/v1.1"
         "--output"
             help = "Path to daily CCFs (jld2)"
             arg_type = String
             default = "output"
=#
         "--checksums"
            help = "Compute md5 checksums for input FITS files."
            action = :store_true
         "--max_num_spectra"
            help = "Maximum number of spectra to process."
            arg_type = Int
            default = 300
         "--crawl"
            help = "Run on each date's subdirectory."
            action = :store_true
         "--quiet"
            help = "Don't write output to terminal"
            action = :store_true
         "--overrwrite"
            help = "Specify it's ok to overwrite files."
            action = :store_true
      end

     return parse_args(s)
 end
 args = parse_commandline_verify_downloads()

compute_checksums = args["checksums"]
input_path = abspath(args["inputdir"])

begin
	"""
	make_meta_redownload( input_path)
	Compare MD5 checksums of data in input_path in meta.csv from NEID archive to actual (or cached values in md5.csv)
	"""
	function make_meta_redownload end

	function make_meta_redownload(input_path::String ; write_md5s::Bool = true, write_meta_redownload::Bool = true, max_files::Integer = 300, checksums::Bool = true )
		meta_filename = joinpath(input_path,"meta.csv")
		manifest_filename = joinpath(input_path,"md5.csv")
		meta_out_filename = joinpath(input_path,"meta_redownload.csv")
		verified_filename = joinpath(input_path,"0_download_verified")
                mixed_versions_filename = joinpath(input_path,"0_mixed_versions")

		if isfile(verified_filename)
                        time_verified = mtime(verified_filename)
                else
                        time_verified = 0
		end
		if isfile(meta_filename)
                        time_meta = mtime(meta_filename)
                else
			@warn("Can't open $meta_filename.")
			return DataFrame()
		end
	        manifest = DataFrame(:Filename=>readdir(input_path,join=true)) |> @filter( occursin(r"\.fits$",_.Filename) ) |> DataFrame
                manifest.mtime = mtime.(manifest.Filename)
                max_fits_mtime = maximum(manifest.mtime)

		if isfile(manifest_filename) && (mtime(manifest_filename) > max_fits_mtime)
                        args["quiet"] || println("# Reusing $manifest_filename.")
                        time_manifest = mtime(manifest_filename)
                        manifest = CSV.read(manifest_filename,DataFrame)
                else
                        # In principle, could only update md5 sums for files modified since creating last md5.csv
                        time_manifest = datetime2unix(now())
                        if checksums
			   args["quiet"] || println("# Computing md5sums's for ", size(manifest,1), " files.")
                           manifest = add_md5s!(manifest, max_files=300)
                           if write_md5s
			      args["quiet"] || println("# Writing $manifest_filename.")
                              CSV.write(manifest_filename,manifest)
                           end
                        end
		end

                if (time_verified >= time_meta) && (time_verified >= time_manifest) && (time_verified >= max_fits_mtime)
			return DataFrame()
                end

                meta = CSV.read(meta_filename,DataFrame)
	        #if args["target"] == "Sun"  # TODO: Eventually need to generalize for non-solar observations
                   meta = meta |> @filter( occursin(r"^Sun$",_.object) ) |> DataFrame
                #end
                m = match(r"(v\d+\.\d+)\.\d+\/", input_path)
                if m != nothing
                   drp_latest_minor = Base.thisminor(VersionNumber(m[1]))
                else
                   drp_latest_minor = Base.thisminor(maximum(VersionNumber.(meta.swversion)))
                end
                if !args["quiet"]
                   println("# Checking minor version matches: ", drp_latest_minor)
                end
                if ! all(Base.thisminor.(VersionNumber.(meta.swversion)) .== drp_latest_minor)
                   msg = "meta.csv contains multiple DRP versions: " * string(unique(meta.swversion))
                   @warn(msg)
                   #touch(mixed_versions_filename)
                   open(mixed_versions_filename, "w") do io
                       write(io, msg)
                   end;
                   # TODO: Should we automatically download new metadata here?
                   # Hard since we won't know all the right query flags (perhaps save them to the directory?)
                   # For now let a script find this file and trigger a reprocess
                end
                manifest.filename = basename.(manifest.Filename)
		need_to_redownload_df = make_meta_redownload(meta,manifest, input_path, checksums=checksums,max_files=max_files)

		if checksums && size(need_to_redownload_df,1) == 0
			touch(verified_filename)
			rm(meta_out_filename,force=true)
		end
		if write_meta_redownload && (size(need_to_redownload_df,1)>=1)
			CSV.write(meta_out_filename,need_to_redownload_df)
		end
		need_to_redownload_df
	end

	function add_md5s!(manifest_df::DataFrame; max_files::Integer = 300, checksums::Bool = true)
		download_stats_df = manifest_df |> @take(max_files) |> @map({_.Filename, mtime_download=mtime(_.Filename), md5_download=string(bytes2hex(open(md5,_.Filename))) })  |> DataFrame
        end

	function make_meta_redownload(meta_df::DataFrame, manifest_df::DataFrame, path::String; max_files::Integer = 300, checksums::Bool = true, suspect_dirname = "suspect" )
                if size(manifest_df,1) == 0
                   return meta_df
                   @warn("Empty dataframe meta: " * string(size(meta_df,1)) * ", manifest: " * string(size(manifest_df,1)) )
                end
                if size(meta_df,1) == 0 || size(manifest_df,1) == 0
                   @warn("Empty dataframe meta: " * string(size(meta_df,1)) * ", manifest: " * string(size(manifest_df,1)) )
                   return DataFrame()
                end
		meta_missing_files = meta_df |> @filter( !(_.l2filename ∈ manifest_df.filename) ) |> DataFrame
		need_to_redownload_df = copy(meta_missing_files)
		if checksums
                        download_success_df = manifest_df |> @join(meta_df, _.filename, _.l2filename, {_.filename, success= _.md5_download == __.l2checksum}) |> DataFrame
			meta_bad_checksum_df = DataFrame(download_success_df |> @filter(!_.success) |> @join(meta_df, _.filename, _.l2filename, {_.filename, meta=__}) |> DataFrame).meta
                        if size(meta_bad_checksum_df,1) >= 1
			   append!(need_to_redownload_df,meta_bad_checksum_df)
                           suspect_files_dir = joinpath(path,suspect_dirname)
                           isdir(suspect_files_dir) || mkdir(suspect_files_dir)
                           for file in meta_bad_checksum_df.filename
                               mv(file, joinpath(suspect_files_dir,file), force=true)
                           end
                        end
                 end
                 begin
                        drp_latest_minor = Base.thisminor(maximum(VersionNumber.(meta_df.swversion)))
                        meta_old_drp_df = meta_df |> @filter( Base.thisminor(VersionNumber(_.swversion)) < drp_latest_minor ) |> DataFrame
                        if size(meta_old_drp_df,1) >= 1
			   append!(need_to_redownload_df,meta_old_drp_df)
                           suspect_files_dir = joinpath(path,suspect_dirname)
                           isdir(suspect_files_dir) || mkdir(suspect_files_dir)
                           for file in meta_old_drp_df.filename
                               mv(file, joinpath(suspect_files_dir,file), force=true)
                           end
                        end
		end
		need_to_redownload_df
	end
end

need_to_redownload_df = DataFrame()

if !args["crawl"]

   if !args["quiet"]
      println("# calling make_meta_redownload(",input_path,", ", "checksums=",compute_checksums,")" )
   end
   need_to_redownload_df = make_meta_redownload(input_path,checksums=compute_checksums)

   if !args["quiet"]
      println("# verify_downloads.jl $input_path identified ", size(need_to_redownload_df,1), " files to download.")
   end

else
  subdirs = DataFrame(:Filename=>readdir(input_path,join=true))
  #println(subdirs.Filename)

  subdirs = subdirs |> @filter( occursin(r"\d{4}\-\d{2}\-\d{2}",_.Filename) ) |> DataFrame
  #println("# Filtered")
  #println(subdirs)

  for dir in subdirs.Filename

      if !args["quiet"]
         println("# calling make_meta_redownload(",dir,", ", "checksums=",compute_checksums,")" )
      end
      global need_to_redownload_df = make_meta_redownload(dir,checksums=compute_checksums)
      if !args["quiet"]
         println("# verify_downloads.jl $dir identified ", size(need_to_redownload_df,1), " files to download.")
      end
  #break
  end

end


text:
@testset "BLAS" begin
    @testset "gemm" begin
        rng = MersenneTwister(1)
        dims = 3:5
        for m in dims, n in dims, p in dims, tA in ('N', 'T'), tB in ('N', 'T')
            α = randn(rng)
            A = randn(rng, tA === 'N' ? (m, n) : (n, m))
            B = randn(rng, tB === 'N' ? (n, p) : (p, n))
            C = gemm(tA, tB, α, A, B)
            ȳ = randn(rng, size(C)...)
            rrule_test(gemm, ȳ, (tA, nothing), (tB, nothing), (α, randn(rng)),
                       (A, randn(rng, size(A))), (B, randn(rng, size(B))))
        end
    end
    @testset "gemv" begin
        rng = MersenneTwister(2)
        for n in 3:5, m in 3:5, t in ('N', 'T')
            α = randn(rng)
            A = randn(rng, m, n)
            x = randn(rng, t === 'N' ? n : m)
            y = α * (t === 'N' ? A : A') * x
            ȳ = randn(rng, size(y)...)
            rrule_test(gemv, ȳ, (t, nothing), (α, randn(rng)), (A, randn(rng, size(A))),
                       (x, randn(rng, size(x))))
        end
    end
end


text:
include("guiding_center_4d_settings.jl")

set_config(:quasi_newton_refactorize, 1)
set_config(:nls_atol, 4eps())
set_config(:nls_rtol, 2eps())
set_config(:nls_stol, 1E-15)
set_config(:nls_nmin, 2)
set_config(:nls_nmax, 20)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GLEW"
version = v"2.1.0"

# Collection of sources required to build GLEW
sources = [
    ArchiveSource("https://github.com/nigels-com/glew/releases/download/glew-$(version)/glew-$(version).tgz",
                  "04de91e7e6763039bc11940095cd9c7f880baba82196a7765f727ac05a993c95")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/glew-*
EXTRA_VARS=()
if [[ "${target}" == *-linux-* ]] || [[ "${target}" == *-freebsd* ]]; then
    # On Linux and FreeBSD this variable by default does `-L/usr/lib`
    EXTRA_VARS+=(LDFLAGS.EXTRA="")
fi
make INCLUDE="-Iinclude -I${includedir}" \
    GLEW_DEST="${prefix}" \
    "${EXTRA_VARS[@]}" \
    install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libGLEW", "glew32"], :libGLEW),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Libglvnd_jll"),
    Dependency("Xorg_libXi_jll"),
    BuildDependency("Xorg_xorgproto_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Test
using JLD2
using Symbolics
using LinearAlgebra
using Random
using BenchmarkTools
using RoboDojo

# Utils
include("utils.jl")

# Solver
include("solver/lu.jl")
include("solver/random_qp.jl")
include("solver/soc.jl")

# Robots
include("robots/hopper.jl")
include("robots/biped.jl")
include("robots/quadruped.jl")

# Simulator
include("simulator/utils.jl")
include("simulator/simulate.jl")
include("simulator/codegen.jl")


text:
# read.jl: Silo variable reading routines

export siloread

function siloread{T1<:AbstractString}(file::silofile, name::T1)
  const curve_int = unsafe_load(DBGetCurve(file.dbfile, name))
#   dump(curve_int)
  const npts = curve_int.npts
#   println("Importing $npts points from ", bytestring(curve_int.title))
#   dump(curve_int.x)
  const y_values_eltype = curve_int.datatype
#   println(y_values_eltype)
#   datatype = Float32
  if haskey(DBdatatype.SiloJuliatypemap, y_values_eltype)
    datatype = DBdatatype.SiloJuliatypemap[y_values_eltype]
  else
    error("Datatype $y_values_eltype currently unsupported.")
  end
  xv=pointer_to_array(reinterpret(Ptr{datatype}, curve_int.x), npts)
  yv=pointer_to_array(reinterpret(Ptr{datatype}, curve_int.y), npts)
  return xv, yv
end


text:
using Gadfly, Distributions

set_default_plot_size(6.6inch, 3.3inch)

x = [0.025, 0.028, 0.035, 0.036, 0.037, 0.039, 0.039, 0.041, 0.041, 0.041, 0.042, 0.043, 0.045, 0.047, 0.049, 0.051, 0.053, 0.053, 0.055, 0.056, 0.057, 0.057, 0.059, 0.063, 0.063, 0.067, 0.067, 0.069, 0.07, 0.072, 0.073, 0.075, 0.076, 0.076, 0.078, 0.082, 0.083, 0.083, 0.084, 0.085, 0.085, 0.086, 0.09, 0.092, 0.092, 0.093, 0.095, 0.095, 0.096, 0.103, 0.103,
0.104, 0.104, 0.105, 0.105, 0.106, 0.107, 0.107, 0.108, 0.109, 0.11, 0.117, 0.119, 0.119, 0.12, 0.126, 0.133, 0.133, 0.135, 0.135, 0.135, 0.137, 0.138, 0.141, 0.144, 0.144, 0.145, 0.145, 0.146, 0.146, 0.149, 0.151, 0.155, 0.155, 0.159, 0.162, 0.162, 0.163, 0.164, 0.167, 0.17, 0.171, 0.172, 0.179, 0.181, 0.187, 0.187, 0.188, 0.188, 0.201, 0.204, 0.204, 0.205, 0.208, 0.208, 0.213, 0.214, 0.214, 0.218, 0.221, 0.222, 0.23, 0.231, 0.238, 0.239, 0.245, 0.245, 0.25, 0.256, 0.256, 0.258, 0.262, 0.262, 0.263, 0.264, 0.267, 0.269, 0.273, 0.273, 0.276, 0.279, 0.288, 0.292, 0.296, 0.296, 0.297, 0.298, 0.301, 0.322, 0.324, 0.328, 0.332, 0.336, 0.341, 0.347, 0.354, 0.357, 0.365, 0.376]

y = [0.018, 4.0, 0.288, 0.157, 0.037, 0.163, 0.06499999, 1.0, 0.09199999, 0.203, 0.26, 0.028, 0.243, 4.0, 2.0, 0.852, 58.0, 4.0, 0.023, 2.48, 4.0, 3.0, 14.0, 37.3, 0.121, 2.0, 12.0,
4.0, 19.0, 4.0, 0.038, 1.0, 2.0, 0.146, 0.061, 0.127, 2.0, 203.0, 2.0, 2.0, 0.278, 2.0, 0.229, 2.0, 0.3, 4.0, 4.0, 0.146, 289.0001, 4.0, 2.0, 2.0, 0.163, 1.0, 0.266, 4.0, 3.0, 1.0, 9.0, 3.0, 0.171, 1.0, 0.264, 10.0, 0.374, 1.0, 2.0, 2.0, 25.0, 2.0, 0.582, 2.0, 15.0, 2.0, 2.0, 3.0, 2.0, 2.0, 2.0, 1.93, 9.0, 2.0, 19.0, 0.684, 4.0, 2.0, 2.0, 2.0, 19.0, 2.0, 6.0, 2.0, 5.0, 2.0, 4.0, 1351.0, 3.0, 2.0, 4.0, 4.0, 12.0, 9.0, 6.0, 6.0, 6.0, 1291.0, 5.0, 79.00001, 28.0, 7.0, 5.0, 60.0, 16.0, 15.0, 14.0, 49.0, 9.0, 19.0, 26.0, 35.0, 79.00001, 26.0, 4.0, 192.0, 79.00001, 26.0, 35.0, 148.0, 31.0, 46.99999, 79.00001, 456.0, 88.0, 21.0, 149.0, 31.0, 465.0001, 19.0, 81.00001, 48.0, 139.0, 938.9999, 303.0, 214.0, 440.9999, 168.0, 367.0, 196.0, 437.9999]

p1 = plot(x=rand(Rayleigh(2),1000), y=rand(Rayleigh(2),1000),
    Geom.density2d(levels = x->maximum(x)*0.5.^collect(1:2:8)), Geom.point,
    Theme(key_position=:none),
    Scale.color_continuous(colormap=x->colorant"red"))

p2 = plot(x=x, y=y, Geom.point,
    layer(Geom.density2d(levels=7), order=2),
    Scale.y_log10, Scale.color_continuous)

hstack(p1, p2)


text:
__precompile__()
module FINUFFT

## Export
export nufft1d1, nufft1d2, nufft1d3
export nufft2d1, nufft2d2, nufft2d3
export nufft3d1, nufft3d2, nufft3d3

export nufft1d1!, nufft1d2!, nufft1d3!
export nufft2d1!, nufft2d2!, nufft2d3!
export nufft3d1!, nufft3d2!, nufft3d3!

export nufftf1d1!, nufftf1d2!, nufftf1d3!
export nufftf2d1!, nufftf2d2!, nufftf2d3!
export nufftf3d1!, nufftf3d2!, nufftf3d3!

export finufft_default_opts
export nufft_opts
export nufft_c_opts # backward-compability

## External dependencies
using finufft_jll

const libfinufft = finufft_jll.libfinufft

const BIGINT = Int64 # defined in src/finufft.h

## FINUFFT opts struct from src/finufft.h
"""
    mutable struct nufft_opts
        modeord             :: Cint
        chkbnds             :: Cint
        debug               :: Cint
        spread_debug        :: Cint
        showwarn            :: Cint
        nthreads            :: Cint
        fftw                :: Cint
        spread_sort         :: Cint
        spread_kerevalmeth  :: Cint
        spread_kerpad       :: Cint
        upsampfac           :: Cdouble
        spread_thread       :: Cint
        maxbatchsize        :: Cint
        spread_nthr_atomic  :: Cint
        spread_max_sp_size  :: Cint
    end

Options struct passed to the FINUFFT library.

# Fields

## Data handling opts

    modeord :: Cint

(type 1,2 only):    0: CMCL-style increasing mode order,
                    1: FFT-style mode order

    chkbnds :: Cint

0: don't check NU pts in [-3pi,3pi),
1: do (<few % slower)

## Diagnostic opts

    debug :: Cint

0: silent,
1: some timing/debug,
2: more

    spread_debug :: Cint

0: silent,
1: some timing/debug,
2: tonnes

    showwarn :: Cint

0: don't print warnings to stderr,
1: do

## Algorithm performance opts

    nthreads :: Cint

number of threads to use, or 0 uses all available

    fftw :: Cint

plan flags to FFTW (`FFTW_ESTIMATE`=64, `FFTW_MEASURE`=0,...)

    spread_sort :: Cint

0: don't sort,
1: do,
2: heuristic choice

    spread_kerevalmeth :: Cint

0: exp(sqrt()) spreading kernel,
1: Horner piecewise poly (faster)

    spread_kerpad :: Cint

option only for exp(sqrt()).
0: don't pad kernel to 4n,
1: do

    upsampfac :: Cdouble

upsampling ratio sigma: 2.0 std, 1.25 small FFT, 0.0 auto

    spread_thread :: Cint

(vectorized ntr>1 only):    0: auto, 1: seq multithreaded,
                            2: parallel single-thread spread

    maxbatchsize :: Cint

option for vectorized ntr>1 only:
max transform batch, 0 auto

    spread_nthr_atomic :: Cint

if >=0, threads above which spreader OMP critical goes atomic

    spread_max_sp_size :: Cint

if >0, overrides spreader (dir=1) max subproblem size

"""
mutable struct nufft_opts
    modeord             :: Cint
    chkbnds             :: Cint
    #
    debug               :: Cint
    spread_debug        :: Cint
    showwarn            :: Cint
    #
    nthreads            :: Cint
    fftw                :: Cint
    spread_sort         :: Cint
    spread_kerevalmeth  :: Cint
    spread_kerpad       :: Cint
    upsampfac           :: Cdouble
    spread_thread       :: Cint
    maxbatchsize        :: Cint
    spread_nthr_atomic  :: Cint
    spread_max_sp_size  :: Cint
end

const nufft_c_opts = nufft_opts # backward compability

"""
    finufft_default_opts()

Return a [`nufft_opts`](@ref) struct with the default FINUFFT settings.\\
See: <https://finufft.readthedocs.io/en/latest/usage.html#options>
"""
function finufft_default_opts()
    opts = nufft_opts(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
    ccall( (:finufft_default_opts, libfinufft),
           Nothing,
           (Ref{nufft_opts},),
           opts
           )
    # default to number of julia threads
    opts.nthreads = Threads.nthreads()
    return opts
end

### Error handling
const WARN_EPS_TOO_SMALL            = 1
const ERR_MAXNALLOC                 = 2
const ERR_SPREAD_BOX_SMALL          = 3
const ERR_SPREAD_PTS_OUT_RANGE      = 4
const ERR_SPREAD_ALLOC              = 5
const ERR_SPREAD_DIR                = 6
const ERR_UPSAMPFAC_TOO_SMALL       = 7
const HORNER_WRONG_BETA             = 8
const ERR_NDATA_NOTVALID            = 9
const ERR_TYPE_NOTVALID             = 10
# some generic internal allocation failure...
const ERR_ALLOC                     = 11
const ERR_DIM_NOTVALID              = 12
const ERR_SPREAD_THREAD_NOTVALID    = 13

struct FINUFFTError <: Exception
    errno::Cint
    msg::String
end
Base.showerror(io::IO, e::FINUFFTError) = print(io, "FINUFFT Error ($(e.errno)): ", e.msg)

function check_ret(ret)
    # Check return value and output error messages
    if ret==0
        return
    elseif ret==WARN_EPS_TOO_SMALL
        @warn "requested tolerance epsilon too small to achieve"
        return
    elseif ret==ERR_MAXNALLOC
        msg = "attemped to allocate internal array larger than MAX_NF (defined in defs.h)"
    elseif ret==ERR_SPREAD_BOX_SMALL
        msg = "spreader: fine grid too small compared to spread (kernel) width"
    elseif ret==ERR_SPREAD_PTS_OUT_RANGE
        msg = "spreader: if chkbnds=1, a nonuniform point coordinate is out of input range [-3pi,3pi]^d"
    elseif ret==ERR_SPREAD_ALLOC
        msg = "spreader: array allocation error"
    elseif ret==ERR_SPREAD_DIR
        msg = "spreader: illegal direction (should be 1 or 2)"
    elseif ret==ERR_UPSAMPFAC_TOO_SMALL
        msg = "upsampfac too small (should be >1.0)"
    elseif ret==HORNER_WRONG_BETA
        msg = "upsampfac not a value with known Horner poly eval rule (currently 2.0 or 1.25 only)"
    elseif ret==ERR_NDATA_NOTVALID
        msg = "ntrans not valid in many (vectorized) or guru interface (should be >= 1)"
    elseif ret==ERR_TYPE_NOTVALID
        msg = "transform type invalid"
    elseif ret==ERR_ALLOC
        msg = "general allocation failure"
    elseif ret==ERR_DIM_NOTVALID
        msg = "dimension invalid"
    elseif ret==ERR_SPREAD_THREAD_NOTVALID
        msg = "spread_thread option invalid"
    else
        msg = "unknown error"
    end
    throw(FINUFFTError(ret, msg))
end

### Simple Interfaces (allocate output)

## Type-1

"""
    nufft1d1(xj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             ms      :: Integer
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-1 1D complex nonuniform FFT.
"""
function nufft1d1(xj::Array{T},
                  cj::Array{Complex{T}},
                  iflag::Integer,
                  eps::T,
                  ms::Integer,
                  opts::nufft_opts=finufft_default_opts()) where T <: Union{Float32,Float64}
    fk = Array{Complex{T}}(undef, ms)
    nufft1d1!(xj, cj, iflag, eps, fk, opts)
    return fk
end

"""
    nufft2d1(xj      :: Array{Float64},
             yj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             ms      :: Integer,
             mt      :: Integer,
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-1 2D complex nonuniform FFT.
"""
function nufft2d1(xj      :: Array{T},
                  yj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: T,
                  ms      :: Integer,
                  mt      :: Integer,
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    fk = Array{Complex{T}}(undef, ms, mt)
    nufft2d1!(xj, yj, cj, iflag, eps, fk, opts)
    return fk
end

"""
    nufft3d1(xj      :: Array{Float64},
             yj      :: Array{Float64},
             zj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             ms      :: Integer,
             mt      :: Integer,
             mu      :: Integer,
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-1 3D complex nonuniform FFT.
"""
function nufft3d1(xj      :: Array{T},
                  yj      :: Array{T},
                  zj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: T,
                  ms      :: Integer,
                  mt      :: Integer,
                  mu      :: Integer,
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    fk = Array{Complex{T}}(undef, ms, mt, mu)
    nufft3d1!(xj, yj, zj, cj, iflag, eps, fk, opts)
    return fk
end

## Type-2

"""
    nufft1d2(xj      :: Array{Float64},
             iflag   :: Integer,
             eps     :: Float64,
             fk      :: Array{ComplexF64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-2 1D complex nonuniform FFT.
"""
function nufft1d2(xj      :: Array{T},
                  iflag   :: Integer,
                  eps     :: T,
                  fk      :: Array{Complex{T}},
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    nj = length(xj)
    cj = Array{Complex{T}}(undef, nj)
    nufft1d2!(xj, cj, iflag, eps, fk, opts)
    return cj
end

"""
    nufft2d2(xj      :: Array{Float64},
             yj      :: Array{Float64},
             iflag   :: Integer,
             eps     :: Float64,
             fk      :: Array{ComplexF64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-2 2D complex nonuniform FFT.
"""
function nufft2d2(xj      :: Array{T},
                  yj      :: Array{T},
                  iflag   :: Integer,
                  eps     :: T,
                  fk      :: Array{Complex{T}},
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    nj = length(xj)
    cj = Array{Complex{T}}(undef, nj)
    nufft2d2!(xj, yj, cj, iflag, eps, fk, opts)
    return cj
end

"""
    nufft3d2(xj      :: Array{Float64},
             yj      :: Array{Float64},
             zj      :: Array{Float64},
             iflag   :: Integer,
             eps     :: Float64,
             fk      :: Array{ComplexF64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-2 3D complex nonuniform FFT.
"""
function nufft3d2(xj      :: Array{T},
                  yj      :: Array{T},
                  zj      :: Array{T},
                  iflag   :: Integer,
                  eps     :: T,
                  fk      :: Array{Complex{T}},
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    nj = length(xj)
    cj = Array{Complex{T}}(undef, nj)
    nufft3d2!(xj, yj, zj, cj, iflag, eps, fk, opts)
    return cj
end

## Type-3

"""
    nufft1d3(xj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             sk      :: Array{Float64},
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-3 1D complex nonuniform FFT.
"""
function nufft1d3(xj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: T,
                  sk      :: Array{T},
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    fk = Array{Complex{T}}(undef, nk)
    nufft1d3!(xj, cj, iflag, eps, sk, fk, opts)
    return fk
end

"""
    nufft2d3(xj      :: Array{Float64},
             yj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             sk      :: Array{Float64},
             tk      :: Array{Float64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-3 2D complex nonuniform FFT.
"""
function nufft2d3(xj      :: Array{T},
                  yj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: T,
                  sk      :: Array{T},
                  tk      :: Array{T},
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    fk = Array{Complex{T}}(undef, nk)
    nufft2d3!(xj, yj, cj, iflag, eps, sk, tk, fk, opts)
    return fk
end

"""
    nufft3d3(xj      :: Array{Float64},
             yj      :: Array{Float64},
             zj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             sk      :: Array{Float64},
             tk      :: Array{Float64}
             uk      :: Array{Float64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-3 3D complex nonuniform FFT.
"""
function nufft3d3(xj      :: Array{T},
                  yj      :: Array{T},
                  zj      :: Array{T},
                  cj      :: Array{Complex{T}},
                  iflag   :: Integer,
                  eps     :: T,
                  sk      :: Array{T},
                  tk      :: Array{T},
                  uk      :: Array{T},
                  opts    :: nufft_opts = finufft_default_opts()) where T <: Union{Float32,Float64}
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    fk = Array{Complex{T}}(undef, nk)
    nufft3d3!(xj, yj, zj, cj, iflag, eps, sk, tk, uk, fk, opts)
    return fk
end

### Direct interfaces (No allocation)

## 1D

"""
    nufft1d1!(xj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-1 1D complex nonuniform FFT. Output stored in fk.
"""
function nufft1d1!(xj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    ms = length(fk)
    # Calling interface
    # int finufft1d1(BIGINT nj,FLT* xj,CPX* cj,int iflag,FLT eps,BIGINT ms,
    # 	       CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft1d1, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, cj, iflag, eps, ms, fk, opts
                )
    check_ret(ret)
end

"""
    nufft1d2!(xj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-2 1D complex nonuniform FFT. Output stored in cj.
"""
function nufft1d2!(xj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    ms = length(fk)
    # Calling interface
    # int finufft1d2(BIGINT nj,FLT* xj,CPX* cj,int iflag,FLT eps,BIGINT ms,
    #                CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft1d2, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, cj, iflag, eps, ms, fk, opts
                )
    check_ret(ret)
end

"""
    nufft1d3!(xj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              sk      :: Array{Float64},
              fk      :: Array{ComplexF64},
              [, opts :: nufft_opts]
             )

Compute type-3 1D complex nonuniform FFT. Output stored in fk.
"""
function nufft1d3!(xj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   sk      :: Array{Float64},
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(fk)==nk
    # Calling interface
    # int finufft1d3(BIGINT nj,FLT* x,CPX* c,int iflag,FLT eps,BIGINT nk, FLT* s, CPX* f, nufft_opts opts);
    ret = ccall( (:finufft1d3, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                Ref{Cdouble},
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, cj, iflag, eps, nk, sk, fk, opts
                )
    check_ret(ret)
end

## 2D

"""
    nufft2d1!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-1 2D complex nonuniform FFT. Output stored in fk.
"""
function nufft2d1!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    ms, mt = size(fk)
    # Calling interface
    # int finufft2d1(BIGINT nj,FLT* xj,FLT *yj,CPX* cj,int iflag,FLT eps,
    #                BIGINT ms, BIGINT mt, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft2d1, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                BIGINT,
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                )
    check_ret(ret)
end

"""
    nufft2d2!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-2 2D complex nonuniform FFT. Output stored in cj.
"""
function nufft2d2!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    ms, mt = size(fk)
    # Calling interface
    # int finufft2d2(BIGINT nj,FLT* xj,FLT *yj,CPX* cj,int iflag,FLT eps,
    #                BIGINT ms, BIGINT mt, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft2d2, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                BIGINT,
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                )
    check_ret(ret)
end

"""
    nufft2d3!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              sk      :: Array{Float64},
              tk      :: Array{Float64},
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
             )

Compute type-3 2D complex nonuniform FFT. Output stored in fk.
"""
function nufft2d3!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   sk      :: Array{Float64},
                   tk      :: Array{Float64},
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(tk)==nk
    @assert length(fk)==nk
    # Calling interface
    # int finufft2d3(BIGINT nj,FLT* x,FLT *y,CPX* cj,int iflag,FLT eps,BIGINT nk, FLT* s, FLT* t, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft2d3, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, yj, cj, iflag, eps, nk, sk, tk, fk, opts
                )
    check_ret(ret)
end

## 3D

"""
    nufft3d1!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              zj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-1 3D complex nonuniform FFT. Output stored in fk.
"""
function nufft3d1!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   zj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    ms, mt, mu = size(fk)
    # Calling interface
    # int finufft3d1(BIGINT nj,FLT* xj,FLT *yj,FLT *zj,CPX* cj,int iflag,FLT eps,
    # 	       BIGINT ms, BIGINT mt, BIGINT mu, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft3d1, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                BIGINT,
                BIGINT,
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                )
    check_ret(ret)
end

"""
    nufft3d2!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              zj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-2 3D complex nonuniform FFT. Output stored in cj.
"""
function nufft3d2!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   zj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    ms, mt, mu = size(fk)
    # Calling interface
    # int finufft3d2(BIGINT nj,FLT* xj,FLT *yj,FLT *zj,CPX* cj,int iflag,FLT eps,
    #                BIGINT ms, BIGINT mt, BIGINT mu, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft3d2, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                BIGINT,
                BIGINT,
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                )
    check_ret(ret)
end

"""
    nufft3d3!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              zj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              sk      :: Array{Float64},
              tk      :: Array{Float64},
              uk      :: Array{Float64},
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
             )

Compute type-3 3D complex nonuniform FFT. Output stored in fk.
"""
function nufft3d3!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   zj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   sk      :: Array{Float64},
                   tk      :: Array{Float64},
                   uk      :: Array{Float64},
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(tk)==nk
    @assert length(uk)==nk
    @assert length(fk)==nk
    # Calling interface
    # int finufft3d3(BIGINT nj,FLT* x,FLT *y,FLT *z, CPX* cj,int iflag,
    #                FLT eps,BIGINT nk,FLT* s, FLT* t, FLT *u,
    #                CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft3d3, libfinufft),
                Cint,
                (BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Cint,
                Cdouble,
                BIGINT,
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{Cdouble},
                Ref{ComplexF64},
                Ref{nufft_opts}),
                nj, xj, yj, zj, cj, iflag, eps, nk, sk, tk, uk, fk, opts
                )
    check_ret(ret)
end

# Load single precision interfaces
include("single.jl")
include("guru.jl")

end # module


text:
using BenchmarkTools
using FinancialMonteCarlo

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2;
mc = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

@btime FwdPrice = pricer(Model, rfCurve, mc, FwdData);
@btime EuPrice = pricer(Model, rfCurve, mc, EUData);
@btime AmPrice = pricer(Model, rfCurve, mc, AMData);
@btime BarrierPrice = pricer(Model, rfCurve, mc, BarrierData);
@btime AsianPrice1 = pricer(Model, rfCurve, mc, AsianFloatingStrikeData);
@btime AsianPrice1 = pricer(Model, rfCurve, mc, AsianFixedStrikeData);

optionDatas = [FwdData, EUData, AMData, BarrierData, AsianFloatingStrikeData, AsianFixedStrikeData]

@btime (FwdPrice, EuPrice, AMPrice, BarrierPrice, AsianPrice1, AsianPrice2) = pricer(Model, rfCurve, mc, optionDatas)


text:
using SavitzkyGolay
using Plots

# Plain call
t = 0:20
y = collect(0:20)
y_ = savitzky_golay(y, 11, 2)
plot(t, [y y_.y], label=["Original signal" "Filtered signal"], ylabel="", xlabel="t", legend=:topleft)

t = LinRange(-4, 4, 500)
y = exp.(-t.^2) .+ 0.05 .* (1.0 .+ randn(length(t)))
y2_ = savitzky_golay(y, 21, 4)
plot(t, [y y2_.y], label=["Original signal" "Filtered signal"], ylabel="", xlabel="t", legend=:topleft)

# With constructor
t = 0:20
y = collect(0:20)
f1 = SGolay(11, 2)
y1 = f1(y)
plot(t, [y y1.y], label=["Original signal" "Filtered signal"], ylabel="", xlabel="t", legend=:topleft)

t = LinRange(-4, 4, 500)
y = exp.(-t.^2) .+ 0.05 .* (1.0 .+ randn(length(t)))
f2 = SGolay(21, 4)
y2 = f2(y)
plot(t, [y y2.y], label=["Original signal" "Filtered signal"], ylabel="", xlabel="t", legend=:topleft)


text:
struct List
    block::Block
    items::Vector{Word}
    scroll::Int
end

function draw(list::List, area::Rect, buf::Buffer)
    draw(list.block, area, buf);
    list_area = inner(list.block, area)

    if width(list_area) < 1 || height(list_area) < 1
        return;
    end

    x, y = left(list_area), top(list_area)

    for (i, item) in enumerate(list.items)
        if i < list.scroll || i - list.scroll > height(list_area)
            continue
        end
        set(buf, x + 1, y, '*')
        set(buf, x + 3, y, item.text, item.style)
        y += 1
    end
end

struct SelectableList
    block::Block
    items::Vector{Word}
    scroll::Int
    selected::Int
end

function draw(list::SelectableList, area::Rect, buf::Buffer)
    list_area = inner(list.block, area)

    draw(list.block, area, buf)
    list_area = inner(list.block, area)

    if width(list_area) < 1 || height(list_area) < 1
        return;
    end

    x, y = left(list_area), top(list_area)

    for (i, item) in enumerate(list.items)
        if i < list.scroll || i - list.scroll > height(list_area)
            continue
        end
        if i == list.selected
            set(buf, x + 1, y, '>')
            set(buf, x + 3, y, item.text, item.style * Crayon(bold = true, italics = true))
        else
            set(buf, x + 1, y, '*')
            set(buf, x + 3, y, item.text, item.style)
        end
        y += 1
    end
end


text:
module PyBase

include("docs.jl")
include("macro.jl")
include("core.jl")
include("julia_api.jl")
include("default.jl")
include("plain.jl")
include("ufcs.jl")
include("supported_methods.jl")

end # module


text:
print_spec(io::IO,t) = print(io,t) #fallback

print_spec(io::IO,x::Type{MOLAR}) = print(io,"Molar ")
print_spec(io::IO,x::Type{MASS}) = print(io,"Mass ")
print_spec(io::IO,x::Type{TOTAL}) = print(io,"Total ")
print_spec(io::IO,x::Type{OneMol}) = print(io,"One mol")
print_spec(io::IO,x::Type{SingleComponent}) = print(io,"Single component")

print_spec(io::IO,x::Type{HumidityDewPoint}) = print(io,"Humidity dew point")
print_spec(io::IO,x::Type{HumidityRatio}) = print(io,"Humidity ratio")
print_spec(io::IO,x::Type{WetBulbTemperature}) = print(io,"Web bulb temperature")
print_spec(io::IO,x::Type{RelativeHumidity}) = print(io,"Relative humidity")
print_spec(io::IO,x::Type{MolarHumidity}) = print(io,"Humidity wet molar fraction")
print_spec(io::IO,x::Type{MassHumidity}) = print(io,"Humidity wet molar fraction")

print_spec(io::IO,::Type{Pressure}) = print(io,"Pressure")
print_spec(io::IO,::Type{Temperature}) = print(io,"Temperature")

function print_spec(io::IO,::Type{Gibbs{T2}}) where T2
    print_spec(io,T2)
    print(io,"Gibbs energy")
end

function print_spec(io::IO,::Type{VaporFraction{T2}}) where T2
    print_spec(io,T2)
    print(io,"vapor fraction")
end

function print_spec(io::IO,::Type{HumiditySpec{T2}}) where T2
    print_spec(io,T2)
end

function print_spec(io::IO,::Type{Helmholtz{T2}}) where T2
    print_spec(io,T2)
    print(io,"Helmholtz energy")
end

function print_spec(io::IO,::Type{Enthalpy{T2}}) where T2
    print_spec(io,T2)
    print(io,"enthalpy")
end

function print_spec(io::IO,::Type{Entropy{T2}}) where T2
    print_spec(io,T2)
    print(io,"entropy")
end

function print_spec(io::IO,::Type{InternalEnergy{T2}}) where T2
    print_spec(io,T2)
    print(io,"internal energy")
end

function print_spec(io::IO,::Type{MaterialAmount{MOLAR}})
    print(io,"Moles")
end

function print_spec(io::IO,::Type{MaterialAmount{MASS}})
    print(io,"Mass")
end

function print_spec(io::IO,::Type{MaterialCompounds{MASS,FRACTION}})
    print(io,"Mass fraction")
end

function print_spec(io::IO,::Type{PhaseTag})
    print(io,"Phase specification")
end

function print_spec(io::IO,::Type{MaterialCompounds{MOLAR,FRACTION}})
    print(io,"Molar fraction")
end

function print_spec(io::IO,::Type{MaterialCompounds{MASS,TOTAL_AMOUNT}})
    print(io,"Mass amounts")
end

function print_spec(io::IO,::Type{MaterialCompounds{MOLAR,TOTAL_AMOUNT}})
    print(io,"Molar amounts")
end

function print_spec(io::IO,::Type{VolumeAmount{T,VOLUME}}) where T
    print_spec(io,T)
    print(io,"volume")
end

function print_spec(io::IO,::Type{VolumeAmount{T,DENSITY}}) where T
    print_spec(io,T)
    print(io,"density")
end

function print_spec(io::IO,::Type{TwoPhaseEquilibrium})
    print(io,"Two phase equilibrium")
end

function print_spec(io::IO,::Type{Options})
    print(io,"Options")
end

function Base.show(io::IO, sp::Spec{T}) where T
    print(io,"spec(")

    a = value(sp)
    units = default_units(T)

    if a isa VariableSpec
        print(io,"x₀)(")
        a = "x₀"
    end
    print(io,string(SPEC_TO_KW[specification(sp)])," = ")

    if is_unitful(a)
        units = unit(a)
        a = ustrip(a)
    end
    print(io,a)
    if is_real(a)
        if !(units == Unitful.NoUnits)
            printstyled(io,'[',color = :light_black)
            print(IOContext(io,:compact => true),units)
            printstyled(io,']',color = :light_black)
        end
    end
    print(io,")")
end

function print_long_spec(io::IO, sp::Spec{T}) where T
    print_spec(io,T)
    print(io," : ")
    a = value(sp)
    units = default_units(T)

    if is_unitful(a)
        units = unit(a)
        a = ustrip(a)
    end
    print(io,a)
    if is_real(a)
        if !(units == Unitful.NoUnits)
            printstyled(io,'[',color = :light_black)
            print(IOContext(io,:compact => true),units)
            printstyled(io,']',color = :light_black)
        end
    end
end

#==function Base.show(::IO, sp::T) where T <:AbstractSpec
    print(io,"::" )
    print_spec(io,T)
end
==#

function Base.show(io::IO, ::MIME"text/plain", sp::ThermodynamicState)
    @nospecialize sp
    lcall = length(sp.callables)
    lsp = length(sp.specs)
    subscripts = '₁':'₉'

    print(io,"ThermodynamicState")
    if !iszero(lcall)
        print(io,"(")
        for (i,spec) in enumerate(sp.callables)
            (i !=1) && print(io,",")
            print(io,"x",subscripts[i])
        end
    print(io,")")
    end

    if isone(lcall+lsp)
        println(io," with 1 property:")
    else
        println(io," with ",lcall+lsp, " properties:")
    end

    if !iszero(lcall)
        for (i,spec) in enumerate(sp.callables)
            (i !=1) && println(io)
            print("  ")
            print_spec(io,typeof(spec))
            print(io," : ")
            print(io,"x",subscripts[i])
        end
        !iszero(lsp) && println(io)
    end

    if !iszero(lsp)
        for (i,spec) in enumerate(sp.specs)
            (i !=1) && println(io)
            print("  ")
            print_long_spec(io,spec)
        end
    end
end

function Base.show(io::IO, sp::ThermodynamicState)
    lcall = length(sp.callables)
    lsp = length(sp.specs)
    subscripts = '₁':'₉'

    print(io,"state(")

    if !iszero(lcall)
        for (i,spec) in enumerate(sp.callables)
            (i !=1) && print(io,",")
            print(io,"x",subscripts[i])
        end
        print(io,")(")

        for (i,spec) in enumerate(sp.callables)
            (i !=1) && print(io,", ")
            print(io,string(SPEC_TO_KW[spec])," = ")
            print(io,"x",subscripts[i])
        end
        !iszero(lsp) && print(io,", ")
    end

    if !iszero(lsp)
        for (i,spec) in enumerate(sp.specs)
            (i !=1) && print(io,", ")
            print(io,string(SPEC_TO_KW[specification(spec)])," = ")
            a = value(spec)
            units = default_units(typeof(specification(spec)))
            if is_unitful(a)
                units = unit(a)
                a = ustrip(a)
            end
            if is_real(a)
                print(IOContext(io,:compact => true),a)
                if !(units == Unitful.NoUnits)
                    printstyled(io,'[',color = :light_black)
                    print(IOContext(io,:compact => true),units)
                    printstyled(io,']',color = :light_black)
                end
            else
                print(io,a)
            end
        end
    end

    print(io,")")
end