text:
"""
From Arb docs:

The mag_t type holds an unsigned floating-point number with a fixed-precision mantissa (30 bits)
and an arbitrary-precision exponent (represented as an fmpz_t), suited for representing magnitude bounds.
The special values zero and positive infinity are supported, but not NaN.

Operations that involve rounding will always produce a valid upper bound, or a lower bound if
the function name has the suffix lower. For performance reasons, no attempt is made to compute
the best possible bounds: in general, a bound may be several ulps larger/smaller than the optimal bound.
Some functions such as mag_set() and mag_mul_2exp_si() are always exact and therefore do not require
separate lower versions.
""" Mag

# using the same field names that Nemo.jl uses

mutable struct Mag
    exp::Int        # fmpz       radius, exponent of 2 (2^exp)
    man::UInt       # mp_limb_t  radius, unsigned by definition

    function Mag()
        z = new()
        ccall(@libarb(mag_init), Cvoid, (Ref{Mag},), z)
        finalizer(mag_clear, z)
        return z
    end
end

mag_clear(x::Mag) = ccall(@libarb(mag_clear), Cvoid, (Ref{Mag},), x)

Mag(x::Mag) = x

swap(x::Mag, y::Mag) where {P} = ccall(@libarb(mag_swap), Cvoid, (Ref{Mag}, Ref{Mag}), x, y)

function copy(x::Mag)
    z = Mag()
    ccall(@libarb(mag_set), Cvoid, (Ref{Mag}, Ref{Mag}), z, x)
    return z
end

function Mag(x::Float64)
    z = Mag()
    ccall(@libarb(mag_set_d), Cvoid, (Ref{Mag}, Cdouble), z, x)
    return z
end

# Mag(-1.0, RoundFromZero)
function Mag(x::Base.IEEEFloat, ::Type{RoundingMode{:FromZero}})
    z = Mag()
    ccall(@libarb(mag_set_d), Cvoid, (Ref{Mag}, Cdouble), z, x)
    return z
end
Mag(x::Base.IEEEFloat, rounding::RoundingMode{:FromZero}) = Mag(x, RoundingMode{:FromZero})
# Mag(-1.0, RoundToZero)
function Mag(x::Base.IEEEFloat, ::Type{RoundingMode{:ToZero}})
    z = Mag()
    ccall(@libarb(mag_set_d_lower), Cvoid, (Ref{Mag}, Cdouble), z, x)
    return z
end
Mag(x::Base.IEEEFloat, rounding::RoundingMode{:ToZero}) = Mag(x, RoundingMode{:ToZero})
# Mag(-1.0)
Mag(x::Base.IEEEFloat) = Mag(x, RoundingMode{:FromZero})

function Mag(x::Int64)
    unsd = UInt64(x)
    z = Mag()
    ccall(@libarb(mag_set_si), Cvoid, (Ref{Mag}, Clong), z, x)
    return z
end

function Mag(x::UInt64)
    z = Mag()
    ccall(@libarb(mag_set_ui), Cvoid, (Ref{Mag}, Culong), z, x)
    return z
end

# Mag(2, RoundFromZero)
function Mag(x::Unsigned, ::Type{RoundingMode{:FromZero}})
    z = Mag()
    ccall(@libarb(mag_set_ui), Cvoid, (Ref{Mag}, Culong), z, x)
    return z
end
Mag(x::Unsigned, rounding::RoundingMode{:FromZero}) = Mag(x, RoundingMode{:FromZero})
# Mag(2, RoundToZero)
function Mag(x::Unsigned, ::Type{RoundingMode{:ToZero}})
    z = Mag()
    ccall(@libarb(mag_set_ui_lower), Cvoid, (Ref{Mag}, Culong), z, x)
    return z
end
Mag(x::Unsigned, rounding::RoundingMode{:ToZero}) = Mag(x, RoundingMode{:ToZero})
# Mag(2)
Mag(x::Unsigned) = Mag(x, RoundingMode{:ToZero})

Mag(x::Signed, ::Type{RoundingMode{:FromZero}}) = Mag(Unsigned(x), RoundingMode{:FromZero})
Mag(x::Signed, rounding::RoundingMode{:FromZero}) = Mag(Unsigned(x), RoundingMode{:FromZero})
Mag(x::Signed, ::Type{RoundingMode{:ToZero}}) = Mag(Unsigned(x), RoundingMode{:ToZero})
Mag(x::Signed, rounding::RoundingMode{:ToZero}) = Mag(Unsigned(x), RoundingMode{:ToZero})
Mag(x::Signed) = Mag(Unsigned(x), RoundingMode{:FromZero})

Float32(x::Mag) = Float32(Float64(x))
Float16(x::Mag) = Float16(Float64(x))


text:
"""
    advance_or_update_ode_integrator!(integrator::DDEIntegrator[, always_calc_begin = false])

Advance or update the ODE integrator of `integrator` to the next time interval by updating its
values and interpolation data with the current values and a full set of interpolation data of
`integrator`.
"""
function advance_or_update_ode_integrator!(integrator, always_calc_begin = false)
  ode_integrator = integrator.integrator

  if ode_integrator.t != integrator.t + integrator.dt
    advance_ode_integrator!(integrator, always_calc_begin)
  else
    update_ode_integrator!(integrator, always_calc_begin)
  end
end

"""
    advance_ode_integrator!(integrator::DDEIntegrator[, always_calc_begin = false])

Advance the ODE integrator of `integrator` to the next time interval by updating its values
and interpolation data with the current values and a full set of interpolation data of
`integrator`.
"""
function advance_ode_integrator!(integrator::DDEIntegrator, always_calc_begin = false)
  @unpack f, u, t, p, k, dt, uprev, alg, cache = integrator
  ode_integrator = integrator.integrator

  # algorithm only works if current time of DDE integrator equals final time point
  # of solution
  t != ode_integrator.sol.t[end] && error("cannot advance ODE integrator")

  # complete interpolation data of DDE integrator for time interval [t, t+dt]
  # and copy it to ODE integrator
  # has to be done before updates to ODE integrator, otherwise history function
  # is incorrect
  if iscomposite(alg)
    DiffEqBase.addsteps!(k, t, uprev, u, dt, f, p, cache.caches[cache.current],
                         always_calc_begin, true, true)
  else
    DiffEqBase.addsteps!(k, t, uprev, u, dt, f, p, cache, always_calc_begin, true, true)
  end
  @inbounds for i in 1:length(k)
    copyat_or_push!(ode_integrator.k, i, k[i])
  end

  # move ODE integrator to interval [t, t+dt]
  ode_integrator.t = t + dt
  ode_integrator.tprev = t
  ode_integrator.dt = dt
  if iscomposite(alg)
    ode_integrator.cache.current = cache.current
  end
  if isinplace(integrator.sol.prob)
    recursivecopy!(ode_integrator.u, u)
  else
    ode_integrator.u = u
  end

  # u(t) is not modified hence we do not have to copy it
  ode_integrator.uprev = ode_integrator.sol.u[end]

  # update prev_idx to index of t and u(t) in solution
  integrator.prev_idx = length(ode_integrator.sol.t)

  nothing
end

"""
    update_ode_integrator!(integrator::DDEIntegrator[, always_calc_begin = false])

Update the ODE integrator of `integrator` by updating its values and interpolation data
with the current values and a full set of interpolation data of `integrator`.
"""
function update_ode_integrator!(integrator::DDEIntegrator, always_calc_begin = false)
  @unpack f, u, t, p, k, dt, uprev, alg, cache = integrator
  ode_integrator = integrator.integrator

  # algorithm only works if the ODE integrator is already moved to the current integration
  # interval
  ode_integrator.t != t + dt && error("cannot update ODE integrator")

  if iscomposite(alg)
    addsteps!(k, t, uprev, u, dt, f, p, cache.caches[cache.current],
              always_calc_begin, true, true)
  else
    addsteps!(k, t, uprev, u, dt, f, p, cache,
              always_calc_begin, true, true)
  end
  @inbounds for i in 1:length(k)
    copyat_or_push!(ode_integrator.k, i, k[i])
  end

  # update state of the dummy ODE solver
  if isinplace(integrator.sol.prob)
    recursivecopy!(ode_integrator.u, u)
  else
    ode_integrator.u = integrator.u
  end

  nothing
end

"""
    move_back_ode_integrator!(integrator::DDEIntegrator)

Move the ODE integrator of `integrator` one integration step back by reverting its values
and interpolation data to the values saved in the dense history.
"""
function move_back_ode_integrator!(integrator::DDEIntegrator)
  ode_integrator = integrator.integrator
  @unpack sol = ode_integrator

  # set values of the ODE integrator back to the values in the solution
  if isinplace(sol.prob)
    recursivecopy!(ode_integrator.u, sol.u[end])
  else
    ode_integrator.u = sol.u[end]
  end
  ode_integrator.t = sol.t[end]
  ode_integrator.tprev = sol.t[integrator.prev2_idx]

  # u(tprev) is not modified hence we do not have to copy it
  ode_integrator.uprev = sol.u[integrator.prev2_idx]

  # revert to the previous time step
  ode_integrator.dt = ode_integrator.dtcache

  # we do not have to reset the interpolation data in the initial time step since always a
  # constant extrapolation is used (and interpolation data of solution at initial
  # time point is not complete!)
  if length(sol.t) > 1
    recursivecopy!(ode_integrator.k, sol.k[end])
  end

  nothing
end

#=
Dealing with discontinuities

If we hit a discontinuity (this is checked in `apply_step!`), then we remove the
discontinuity, additional discontinuities at the current time point (if present), and
maybe also discontinuities and time stops coming shortly after the current time point
in `handle_discontinuities!`. The order of the discontinuity at the current time point is
defined as the lowest order of all these discontinuities.

If the problem is not neutral, we will only add additional discontinuities if
this order is less or equal to the order of the algorithm in
`add_next_discontinuities!`. If we add discontinuities, we add discontinuities
of the next order caused by constant lags (these we can calculate explicitly and
just add them to `d_discontinuities` and `tstops`) and we add the current
discontinuity to `tracked_discontinuities` which is the array of old
discontinuities that are checked by a `DiscontinuityCallback` (if existent).
=#

# handle discontinuities at the current time point of the `integrator`
function OrdinaryDiffEq.handle_discontinuities!(integrator::DDEIntegrator)
    # remove all discontinuities at current time point and calculate minimal order
    # of these discontinuities
    d = OrdinaryDiffEq.pop_discontinuity!(integrator)
    order = d.order
    tdir_t = integrator.tdir * integrator.t
    while OrdinaryDiffEq.has_discontinuity(integrator) &&
        OrdinaryDiffEq.first_discontinuity(integrator) == tdir_t

        d2 = OrdinaryDiffEq.pop_discontinuity!(integrator)
        order = min(order, d2.order)
    end

    # remove all discontinuities close to the current time point as well and
    # calculate minimal order of these discontinuities
    # integrator.EEst has unitless type of integrator.t
    if typeof(integrator.EEst) <: AbstractFloat
        maxΔt = 10eps(integrator.t)

        while OrdinaryDiffEq.has_discontinuity(integrator) &&
            abs(OrdinaryDiffEq.first_discontinuity(integrator).t - tdir_t) < maxΔt

            d2 = OrdinaryDiffEq.pop_discontinuity!(integrator)
            order = min(order, d2.order)
        end

        # also remove all corresponding time stops
        while OrdinaryDiffEq.has_tstop(integrator) &&
            abs(OrdinaryDiffEq.first_tstop(integrator) - tdir_t) < maxΔt

            OrdinaryDiffEq.pop_tstop!(integrator)
        end
    end

    # add discontinuities of next order to integrator
    add_next_discontinuities!(integrator, order)

  nothing
end

"""
    add_next_discontinuities!(integrator::DDEIntegrator, order[, t=integrator.t])

Add discontinuities of next order that are propagated from discontinuity of
order `order` at time `t` in `integrator`, but only if `order` is less or equal
than the order of the applied method or the problem is neutral.

Discontinuities caused by constant delays are immediately calculated, and
discontinuities caused by dependent delays are tracked by a callback.
"""
function add_next_discontinuities!(integrator, order, t=integrator.t)
  neutral = integrator.sol.prob.neutral
  next_order = neutral ? order : order + 1

  # only track discontinuities up to order of the applied method
  alg_maximum_order = OrdinaryDiffEq.alg_maximum_order(integrator.alg)
  next_order <= alg_maximum_order + 1 || return

  # discontinuities caused by constant lags
  if has_constant_lags(integrator)
    constant_lags = integrator.sol.prob.constant_lags
    maxlag = integrator.tdir * (integrator.sol.prob.tspan[end] - t)

    for lag in constant_lags
      if integrator.tdir * lag < maxlag
        # calculate discontinuity and add it to heap of discontinuities and time stops
        d = Discontinuity(integrator.tdir * (t + lag), next_order)
        push!(integrator.d_discontinuities_propagated, d)
        push!(integrator.tstops_propagated, d.t)
      end
    end
  end

  # track propagated discontinuities with callback
  push!(integrator.tracked_discontinuities, Discontinuity(integrator.tdir * t, order))

  nothing
end

# Interface for accessing and removing next time stops and discontinuities
function OrdinaryDiffEq.has_tstop(integrator::DDEIntegrator)
    return _has(integrator.opts.tstops, integrator.tstops_propagated)
end
function OrdinaryDiffEq.first_tstop(integrator::DDEIntegrator)
    return _first(integrator.opts.tstops, integrator.tstops_propagated)
end
function OrdinaryDiffEq.pop_tstop!(integrator::DDEIntegrator)
    return _pop!(integrator.opts.tstops, integrator.tstops_propagated)
end

function OrdinaryDiffEq.has_discontinuity(integrator::DDEIntegrator)
    return _has(integrator.opts.d_discontinuities,
                integrator.d_discontinuities_propagated)
end
function OrdinaryDiffEq.first_discontinuity(integrator::DDEIntegrator)
    return _first(integrator.opts.d_discontinuities,
                  integrator.d_discontinuities_propagated)
end
function OrdinaryDiffEq.pop_discontinuity!(integrator::DDEIntegrator)
    return _pop!(integrator.opts.d_discontinuities,
                 integrator.d_discontinuities_propagated)
end

_has(x, y) = !isempty(x) || !isempty(y)
function _first(x, y)
    if isempty(x)
        return first(y)
    elseif isempty(y)
        return first(x)
    else
        return min(first(x), first(y))
    end
end
function _pop!(x, y)
    if isempty(x)
        pop!(y)
    elseif isempty(y)
        pop!(x)
    else
        if first(x) < first(y)
            pop!(x)
        else
            pop!(y)
        end
    end
end

ext: jl

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
export CartesianTrajectory3D, cartesian3dNodes, cartesian3dDensity

"""
    CartesianTrajectory3D(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , numSlices=1
                  , kargs...)

returns a 3d cartesian trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`numSlices=1`)                 - number of slices
"""
function CartesianTrajectory3D(::Type{T}, numProfiles, numSamplingPerProfile
                  ; TE=0.0
                  , AQ=1.e-3
                  , numSlices=1
                  , kargs...) where T
  nodes = cartesian3dNodes(T, numProfiles, numSamplingPerProfile; numSlices=numSlices)
  times = readoutTimes(T, numProfiles, numSamplingPerProfile, numSlices; TE=TE, AQ=AQ)
  return  Trajectory("Cartesian3D", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, numSlices, true, false)
end

function cartesian3dNodes(::Type{T}, numProfiles, numSamplingPerProfile
                          ; numSlices=1, kargs...) where T
  nodes = zeros(T, 3, numSamplingPerProfile, numProfiles, numSlices)
  posX = collect( -ceil(Int64, (numSamplingPerProfile-1)/2.):floor(Int64, (numSamplingPerProfile-1)/2.) ) / numSamplingPerProfile
  posY = collect( -ceil(Int64, (numProfiles-1)/2.):floor(Int64, (numProfiles-1)/2.) ) / numProfiles
  posZ = collect( -ceil(Int64, (numSlices-1)/2.):floor(Int64, (numSlices-1)/2.) ) / numSlices

  for j = 1:numSlices
    for l = 1:numProfiles
      for k = 1:numSamplingPerProfile
        nodes[1,k,l,j] = posX[k]
        nodes[2,k,l,j] = posY[l]
        nodes[3,k,l,j] = posZ[j]
      end
    end
  end
  return reshape(nodes, 3, numSamplingPerProfile*numProfiles*numSlices)
end

function cartesian3dDensity(::Type{T}, numSamplingPerProfile::Int64, numProfiles::Int64, numSlices::Int64) where T
  density = zeros(T, numSamplingPerProfile, numProfiles, numSlices)
  for j = 1:numSlices
    for l = 1:numProfiles
      for k = 1:numSamplingPerProfile
        density[k,l,j] = 1
      end
    end
  end
  return vec(density)
end


text:
@reexport using JLD
import JLD: save, load

export WritePlan,StorePlan, initialize_storage, store_data!

"""
    WritePlan(file,write_Δt,varlist)

Create a plan for writing data to file. The filename `file` is specified,
to be written to every `write_Δt` time units. The variable names to be written are
specified as a set of dictionary pairs in `varlist`.
"""
struct WritePlan
    filen::String
    write_Δt::Float64
    vardict::Dict
    WritePlan(filen,write_Δt,varlist...) = new(filen,write_Δt,_get_typedict(varlist))
end

"""
    save(t,R::WritePlan,v)

Check if time `t` is appropriate for writing to file, according to the WritePlan `R`,
and if so, write the specified variables `v` (which may be separated by commas)
in `R` to the file in `R`.
"""
function save(t,R::WritePlan,v...)
    tol = 1e-8
    if (isapprox(mod(t,R.write_Δt),0,atol=tol) ||
        isapprox(mod(t,R.write_Δt),R.write_Δt,atol=tol))
        save(R.filen,_writelist(R,v...)...)
    end
end

"""
    load(R::WritePlan[,var::String])

Load the data stored in the file specified by the WritePlan `R`. If
variable `var` is specified, outputs the specified variable.
"""
load(R::WritePlan) = load(R.filen)
load(R::WritePlan,v...) = load(R.filen,v...)

function _writelist(R::WritePlan,varlist...)
    @assert length(R.vardict)==length(varlist)
    list = ()
    for v in varlist
        @assert haskey(R.vardict,v.first) "No such variable exists in this list"
        @assert typeof(v.second) == R.vardict[v.first] "Invalid type of data for this entry"
        list = (list...,v.first,v.second)
    end
    return list
end

"""
    StorePlan(sample_time,varlist[;htype=RegularHistory,min_time=-Inf,max_time=Inf])

Create a plan for storing history data. The storage of data is specied to occur
at a sample interval given by `sample_time`. Optionally, one can also specify the
initial time `min_time` and final time `max_time` (included) for a sampling window.
These default to all times. The list of variables to be stored is
specified as a list of variables `varlist`, provided as comma-separated named pairs,
e.g. `"varname1"-> v1,"varname2"-> v2`. The optional argument `htype` can be used to set the
history data to `PeriodicHistory` or `RegularHistory` (the default). In the
case of `PeriodicHistory`, the data is assumed to repeat with a period equal
to length(history)+1.
"""
struct StorePlan{H<:HistoryType}
    min_t::Float64
    max_t::Float64
    store_Δt::Float64
    vardict::Dict
    StorePlan(store_Δt,varlist...;htype=RegularHistory,min_time=-Inf,max_time=Inf) =
              new{htype}(min_time,max_time,store_Δt,_get_typedict(varlist))
end

"""
    initialize_storage(S::StorePlan) -> Dict

Initialize a storage data stack for the storage plan `S`. The output is
an empty vector of `History` vectors.
"""
function initialize_storage(S::StorePlan{H}) where {H}
    data = Dict{String,History}()
    for v in S.vardict
        get!(data,v.first,History(v.second,htype=H))
    end
    return data
end

"""
    store_data!(data,t,S::StorePlan,v)

Check whether time `t` is a time for saving for storage as described by plan `S`,
and if so, push the variables specified in `v` (a list of pairs) onto the `data` stack.
Note that a deepcopy of the values of each variable are stored automatically.
"""
function store_data!(data,t,S::StorePlan,v...)
  tol = 1e-8
  if t >= (S.min_t-tol) && t <= (S.max_t + tol) &&
      ((isapprox(mod(t,S.store_Δt),0,atol=tol) ||
        isapprox(mod(t,S.store_Δt),S.store_Δt,atol=tol)))
        _store_data!(data,S::StorePlan,v...)
  end
  return data
end

function _store_data!(data,S::StorePlan,varlist...)
    for v in varlist
        @assert haskey(data,v.first) "No such variable exists in this list"
        @assert typeof(v.second) == S.vardict[v.first] "Invalid type of data for this entry"
        push!(data[v.first],deepcopy(v.second))
    end
    return data
end

function _get_typedict(varlist)
    newdict = Dict{String,Type}()
    for l in varlist
        get!(newdict,l.first,typeof(l.second))
    end
    return newdict
end


text:
## Main function to parse a template This works in several steps: each
## character parsed into a token, the tokens are squashed, then
## nested, then rendered.

"""
   Mustache.parse(template, tags = ("{{", "}}"))

Parse a template into tokens.

* `template`: a string containing a template
* `tags`: the tags used to indicate a variable. Adding interior braces (`{`,`}`) around the
variable will prevent HTML escaping. (That is for the default tags, `{{{varname}}}` is used; for
tags like `("<<",">>")` then `<<{varname}>>` is used.)

"""
function parse(template, tags = ("{{", "}}"))
    tokens = make_tokens(template, tags)
    out = nestTokens(tokens)
    MustacheTokens(out)
end

"""
    mt"string"

Macro to parse tokens from a string.

"""
macro mt_str(s)
    parse(s)
end

"""
    jmt"string"

Macro that interpolates values escaped by dollar signs, then parses strings.

Note: modified from [HypertextLiteral](https://github.com/clarkevans/HypertextLiteral.jl/blob/master/src/macros.jl)

Example:

```
x = 1
toks = jmt"\$(2x) by {{:a}}"
toks(a=2) # "2 by 2"
```
"""
macro jmt_str(expr::String)
    # Essentially this is an ad-hoc scanner of the string, splitting
    # it by `$` to find interpolated parts and delegating the hard work
    # to `Meta.parse`, treating everything else as a literal string.

    args = Any[]
    start = idx = 1
    strlen = length(expr)
    while true
        idx = findnext(isequal('$'), expr, start)
        if idx == nothing
           push!(args, expr[start:strlen])
           break
        end
        push!(args, expr[start:idx-1])
        start = idx + 1
        (nest, tail) = Meta.parse(expr, start; greedy=false)
        if nest == nothing
            throw("missing interpolation expression")
        end
        if !(expr[start] == '(' || nest isa Symbol)
            throw(DomainError(nest,
             "interpolations must be symbols or parenthesized"))
        end
        if Meta.isexpr(nest, :(=))
            throw(DomainError(nest,
             "assignments are not permitted in an interpolation"))
        end
        if nest isa String
            # this is an interpolated string literal
            nest = Expr(:string, nest)
        end
        push!(args, nest)
        start = tail
    end

    quote
        Mustache.parse(string($(map(esc, args)...)))
    end

end

## Dict for storing parsed templates
TEMPLATES = Dict{AbstractString, MustacheTokens}()

"""
    Mustache.load(filepath, args...)

Load file specified through  `filepath` and return the compiled tokens.
Tokens are memoized for efficiency,

Additional arguments are passed to `Mustache.parse` (for adjusting the tags).
"""
function load(filepath, args...)

    isfile(filepath) || throw(ArgumentError("File $filepath not found"))

    key = string(mtime(filepath)) * filepath * string(hash(args))
    haskey(TEMPLATES,key) && return  TEMPLATES[key]

    open(filepath) do s
        global tpl = parse(read(s, String), args...)
    end

    TEMPLATES[key] = tpl
    tpl

end

# old names. To be deprecated?
const template_from_file = load


text:
const FloatOrComplex = Union{Float32, Float64, ComplexF32, ComplexF64}

"""
   define function as linear operator
"""
struct LinearOperator{T<:DataType, Ti<:Int64, F1, F2}

    data_format       :: T

    # forward operator
    forward_operator  :: F1
    forward_params    :: Union{Nothing, NamedTuple}
    is_forward_mutate :: Bool

    # adjoint operator
    adjoint_operator  :: F2
    adjoint_params    :: Union{Nothing, NamedTuple}
    is_adjoint_mutate :: Bool

    # size of the linear operator
    m                 :: Ti
    n                 :: Ti

    # some properties
    is_symmetrical    :: Bool
    is_hermitian      :: Bool
end

"""
   Constructor for linear operator which consist of forward and its adjoint operators
"""
function LinearOperator(data_format::DataType, forward_op::F1, adjoint_op::F2, m::Ti, n::Ti;
                        forward_params=nothing, adjoint_params=nothing)  where {Ti<:Int64, F1, F2}

    # check data fromat
    data_format <: FloatOrComplex || error("only support float or complex")

    # forward operator mutating or not
    nargs = first(methods(forward_op)).nargs
    (nargs == 2 || nargs == 3) || error("wrong number of input to forward operator")
    is_forward_mutate =  (nargs == 3)

    # adjoint operator mutating or not
    nargs = first(methods(adjoint_op)).nargs
    (nargs == 2 || nargs == 3) || error("wrong number of input to adjoint operator")
    is_adjoint_mutate =  (nargs == 3)

    # proerties
    is_symmetrical = false
    is_hermitian   = false

    return LinearOperator(data_format,
                          forward_op, forward_params, is_forward_mutate,
                          adjoint_op, adjoint_params, is_adjoint_mutate,
                          m, n,
                          is_symmetrical, is_hermitian)
end

"""
   Constructor for symmetrical or hermitian linear operator
"""
function LinearOperator(data_format::DataType, forward_op::F1, m::Ti; forward_params=nothing) where {Ti<:Int64, F1, F2}

    # check data fromat
    data_format <: FloatOrComplex || error("only support float or complex")

    # operator support in-place operation
    is_forward_mutate = first(methods(forward_op)).nargs == 3

    # adjoint is same as forward operator
    adjoint_op     = forward_op
    adjoint_params = forward_params
    is_adjoint_mutate = is_forward_mutate
    n = m

    # proerties
    if data_format <: AbstractFloat
       is_symmetrical = true
       is_hermitian   = false
    else
       is_symmetrical = false
       is_hermitian   = true
    end

    return LinearOperator(data_format,
                          forward_op, forward_params, is_forward_mutate,
                          adjoint_op, adjoint_params, is_adjoint_mutate,
                          m, n,
                          is_symmetrical, is_hermitian)
end

"""
   Customerize printing function for linear operator
"""
function Base.show(io::IO, L::LinearOperator)
    @printf("data_format = %10s\n", L.data_format)
    @printf("forward operator = %20s, ismutating = %6s\n", L.forward_operator, L.is_forward_mutate)
    @printf("adjoint operator = %20s, ismutating = %6s\n", L.adjoint_operator, L.is_adjoint_mutate)
    @printf("size of linear operator = (%d, %d)\n", L.m, L.n)
    @printf("issymmetrical = %6s, ishermitian = %6s\n", L.is_symmetrical, L.is_hermitian)
    return nothing
end

# enquire the properties of linear operator
Base.eltype(L::LinearOperator) = L.data_format
Base.ndims(L::LinearOperator)  = 2
Base.size(L::LinearOperator)   = (L.m, L.n)
Base.size(L::LinearOperator, n::Integer) = (n==1 || n==2) ? size(L)[n] : error("Linear operator only have two dimensions")
Base.isreal(L::LinearOperator) = L.data_format <: AbstractFloat
LinearAlgebra.issymmetric(L::LinearOperator) = L.is_symmetrical
LinearAlgebra.ishermitian(L::LinearOperator) = L.is_hermitian
ismutating(L::LinearOperator) = (L.is_forward_mutate, L.is_adjoint_mutate)

"""
   Apply the linear operator to an Abstract vector
"""
function A_mul_b!(y::AbstractVector, L::LinearOperator, x::AbstractVector)

    # check the length of input and output
    (length(x) == L.n && length(y) == L.m) || throw(DimensionMismatch("check the length of x or y"))

    # apply forward operator
    if ismutating(L)[1]

       L.forward_params == nothing ? L.forward_operator(y, x) : L.forward_operator(y, x; L.forward_params...)

    else

       L.forward_params == nothing ? copyto!(y,L.forward_operator(x)) : copyto!(y,L.forward_operator(x; L.forward_params...))
    end

    return nothing
end

"""
   define the multiplication operator
"""
function Base.:(*)(L::LinearOperator, x::AbstractVector)

    # check the length of input
    length(x) == L.n || throw(DimensionMismatch())

    if ismutating(L)[1]

       # allocate memory
       y = similar(x, promote_type(eltype(L), eltype(x)), L.m)
       L.forward_params == nothing ? L.forward_operator(y, x) : L.forward_operator(y, x; L.forward_params...)

    else

       L.forward_params == nothing ? y = L.forward_operator(x) : y = L.forward_operator(x; L.forward_params...)

    end

    return y
end

"""
   apply the adjoint operator
"""
function Ac_mul_b!(y::AbstractVector, L::LinearOperator, x::AbstractVector)

    # check the length of input and output
    (length(x) == L.m && length(y) == L.n) || throw(DimensionMismatch("check the length of x or y"))

    # apply adjoint operator
    if ismutating(L)[2]

       L.adjoint_params == nothing ? L.adjoint_operator(y, x) : L.adjoint_operator(y, x; L.adjoint_params...)

    else

       L.adjoint_params == nothing ? copyto!(y,L.adjoint_operator(x)) : copyto!(y,L.adjoint_operator(x; L.adjoint_params...))
    end

    return nothing
end

"""
   apply the adjoint operator, At_mulb! is equal to Ac_mul_b! when the linear operator is real
"""
function At_mul_b!(y::AbstractVector, L::LinearOperator, x::AbstractVector)

    # check the length of input and output
    (length(x) == L.m && length(y) == L.n) || throw(DimensionMismatch("check the length of x or y"))

    if !isreal(L)
       x = conj(x)
    end

    # apply adjoint operator
    if ismutating(L)[2]

       L.adjoint_params == nothing ? L.adjoint_operator(y, x) : L.adjoint_operator(y, x; L.adjoint_params...)

    else

       L.adjoint_params == nothing ? copyto!(y,L.adjoint_operator(x)) : copyto!(y,L.adjoint_operator(x; L.adjoint_params...))
    end

    if !isreal(L)
       conj!(y)
    end

    return nothing
end

"""
   get the corresponding matrix of a linear operator
"""
function Base.Matrix(L::LinearOperator)

    # size of linear operator
    (m, n) = size(L)

    # allocate space
    A = Matrix{L.data_format}(undef, L.m, L.n)

    # input vector
    x = zeros(L.data_format, n)

    # probing column by column
    @inbounds for i = 1 : n
         x[i] = one(L.data_format)
         A_mul_b!(view(A,:,i), L, x)

         # prepare for next column
         x[i] = zero(L.data_format)
    end

    return A
end

"""
   get the corresponding sparse matrix of a linear operator
"""
function SparseArrays.sparse(L::LinearOperator)

    # size of linear operator
    (m, n) = size(L)

    # data format
    Tv = eltype(L)

    # the variable for sparse matrix
    rowval = Int64[]
    colptr = zeros(Int64, n+1)
    nzval  = Tv[]

    # vector input
    x = zeros(Tv, n)
    y = Vector{Tv}(undef, m)

    # loop over all columns
    @inbounds for i = 1 : n
        x[i] = one(Tv)
        A_mul_b!(y, L, x)
        js = findall(!iszero, y)
        colptr[i] = length(nzval) + 1

        if length(js) > 0
           append!(rowval, js)
           append!(nzval, y[js])
        end

        x[i] = zero(Tv)
    end
    colptr[n+1] = length(nzval) + 1

    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end


text:
using CoulombIntegrals
using CompactBases
using FillArrays
using LazyArrays
using LinearAlgebra
using BandedMatrices
using Parameters
using Test

include("exact_poisson.jl")
include("hydrogenic_orbitals.jl")
include("integrals.jl")
include("test_convergence_rates.jl")

@testset "Basic test" begin
    function get_grid(grid_type, rmax, Z; ρ=0.1, ρmax=0.6, α=0.002)
        ρ /= Z

        if grid_type == :fd_uniform
            N = ceil(Int, rmax/ρ)
            FiniteDifferences(N, ρ)
        elseif grid_type == :fd_loglin
            StaggeredFiniteDifferences(0.1ρ, ρmax, α, rmax)
        elseif grid_type == :implicit_fd
            N = ceil(Int, rmax/ρ)
            ImplicitFiniteDifferences(N, ρ)
        elseif grid_type == :fedvr
            t = range(0, stop=rmax, length=41)
            FEDVR(t, Vcat(10,Fill(7,length(t)-2)))[:,2:end-1]
        elseif grid_type == :bsplines
            t = ExpKnotSet(5, -2.0, log10(rmax), 200)
            BSpline(t)[:,2:end-1]
        end
    end

    n,ℓ = 3,0
    n′,ℓ′ = 2,0
    k = abs(ℓ′-ℓ) # Multipole order

    Z = 1.0
    rmax = 40(max(n,n′)^1.05)/Z

    grid_types = [(:fd_uniform, 2e-4, 4e-8),
                  (:fd_loglin, 1e-6, 5e-10),
                  (:implicit_fd, 2e-5, 4e-8),
                  (:fedvr, 3e-7, 8e-11),
                  (:bsplines, 5e-9, 2e-9)]

    # All in the name of type-stability
    Yk_scaled(Yk_f) = (r,Z) -> Yk_f(Z*r)
    Yk_ref = exact_Yᵏs[((n,ℓ),(n′,ℓ′),k)]
    fYk = Base.Fix2(Yk_scaled(Yk_ref), Z)

    ϕf = hydredwfn(n, ℓ, Z)
    ϕ′f = hydredwfn(n′, ℓ′, Z)

    @testset "Grid = $(grid_type)" for (grid_type,tol,asymtol) in grid_types
        R = get_grid(grid_type, rmax, Z, ρ=0.1, ρmax=1.0)
        r = axes(R,1)

        Ỹᵏ = R \ fYk.(r)

        ϕ = R \ ϕf.(r)
        ϕ′ = R \ ϕ′f.(r)

        o = R \ one.(r)
        out1 = similar(o)
        out2 = similar(o)

        ρ = Density(applied(*,R,ϕ), applied(*,R,ϕ′))
        potential = CoulombRepulsionPotential(R, k)
        copyto!(potential, ρ)
        Y = potential.poisson.Y
        @test Y ≈ Ỹᵏ atol=tol

        mpot = Matrix(potential)
        @test mpot isa (R isa CompactBases.BSplineOrRestricted ? BandedMatrix : Diagonal)

        mul!(out1, potential, o)
        mul!(out2, LinearOperator(mpot, R), o)
        if R isa CompactBases.BSplineOrRestricted
            # Matricization does not yield the same behaviour for the
            # first few elements since the Dirichlet0 boundary
            # condition at r=0 does not allow the accurate
            # representation of Y/r.
            @test out1[55:end] ≈ out2[55:end] atol=tol
        else
            @test out1 ≈ out2 atol=tol
        end

        @testset "AsymptoticPoissonProblem" begin
            R̃ = R[:,1:end-5]
            apoisson = AsymptoticPoissonProblem(R, k, R̃)
            CoulombIntegrals.solve!(apoisson, ρ)
            @test Y ≈ apoisson.Y atol=asymtol
        end
    end
end

include("poisson_problem_hermiticity.jl")

@testset "Accuracy tests" begin
    rₘₐₓ = 100.0

    orbital_modes = [:symbolic]
    # # Arnoldi is not stable with very small grid spacings;
    # # shift-and-invert has to be implemented.
    # orbital_modes = [:symbolic, :arnoldi]

    coulomb_modes = [:poisson]
    # # Direct integration is simply too slow.
    # coulomb_modes = [:poisson, :direct]
    Z = 1.0

    fedvr_order = 10

    fedvr_inputs(exponents,errors) =
        ("FEDVR",
         (rₘₐₓ,N) -> begin
         t = range(0,stop=rₘₐₓ,length=N)
         amended_order = vcat(fedvr_order+2, Fill(fedvr_order,length(t)-2))
         R = FEDVR(t, amended_order)[:,2:end-1]
         R,(t[2]-t[1])/fedvr_order
         end,
         unique(ceil.(Int, 2 .^ exponents)),
         errors)

    fd_inputs(exponents,errors) =
        ("Finite-differences",
         (rₘₐₓ,N) -> begin
         R = StaggeredFiniteDifferences(rₘₐₓ,N)
         R,step(R)
         end,
         unique(ceil.(Int, 2 .^ exponents)),
         errors)

    # TODO: Add B-spline convergence tests

    @testset "Yᵏ convergence rates" begin
        @testset "$Basis" for (Basis,Rf,Ns,errors) in [
            fedvr_inputs(range(1,stop=5,length=20),
                         [(6, "local", 2e-5, 7.0)]),
            fd_inputs(range(6.5,stop=10,length=15),
                      [(6, "local", 1e-3, 2.0)])
        ]
            @test test_convergence_rates(
                Rf, rₘₐₓ, Ns, 2,
                errors,
                ["N", "ρ", "a", "b", "k", "Local δ", "Global δ", "Time"],
                orbital_modes, coulomb_modes,
                sort(collect(keys(exact_Yᵏs)))) do cache, orbital_mode, coulomb_mode, ((n,ℓ), (n′,ℓ′), k)
                    (n′ > n || n′ == n && ℓ′ > ℓ) && return nothing
                    Yᵏ_error(cache,
                             k, n, ℓ, n′, ℓ′, Z,
                             orbital_mode, coulomb_mode)
                end
        end
    end

    @testset "Fᵏ convergence rates" begin
        @testset "$Basis" for (Basis,Rf,Ns,errors) in [
            fedvr_inputs(range(2,stop=6,length=20),
                         [(11,"δF",2e-6,6.0),
                          (13,"δF′",2e-6,6.0)]),
            fd_inputs(range(7,stop=11,length=15),
                      [(11,"δF",2e-4,2.0),
                       (13,"δF′",2e-4,2.0)])
        ]
            @test test_convergence_rates(
                Rf, rₘₐₓ, Ns, 2,
                errors,
                ["N", "ρ", "a", "b", "c", "d", "k", "Exact", "Exact",
                 "F", "δ", "F′", "δ′", "F-F′", "Time"],
                orbital_modes, coulomb_modes, collect(keys(exact_Fᵏs))) do cache, orbital_mode, coulomb_mode, ((n,ℓ),(n′,ℓ′),k)
                    Fᵏ_error(cache,
                             k, n, ℓ, n′, ℓ′,
                             Z, orbital_mode, coulomb_mode)
                end
        end
    end

    @testset "Gᵏ convergence rates" begin
        @testset "$Basis" for (Basis,Rf,Ns,errors) in [
            fedvr_inputs(range(2,stop=6,length=20),
                         [(11,"δG",2e-6,6.0)]),
            fd_inputs(range(7,stop=12,length=15),
                      [(11,"δG",2e-4,2.0)])
        ]
            @test test_convergence_rates(
                Rf, rₘₐₓ, Ns, 2,
                errors,
                ["N", "ρ", "a", "b", "c", "d", "k", "Exact", "Exact",
                 "G", "δ", "Time"],
                orbital_modes, coulomb_modes, collect(keys(exact_Gᵏs))) do cache, orbital_mode, coulomb_mode, ((n,ℓ),(n′,ℓ′),k)
                    Gᵏ_error(cache,
                             k, n, ℓ, n′, ℓ′,
                             Z, orbital_mode, coulomb_mode)
                end
        end
    end
end


text:
using KissABC
using Documenter

makedocs(;
    modules = [KissABC],
    authors = "Francesco Alemanno <francescoalemanno710@gmail.com> and contributors",
    repo = "https://github.com/francescoalemanno/KissABC.jl/blob/{commit}{path}#L{line}",
    sitename = "KissABC.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://francescoalemanno.github.io/KissABC.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
#        "Example: Gaussian Mixture" => "example_1.md",
        #"Reference" => "reference.md",
    ],
)

deploydocs(; repo = "github.com/francescoalemanno/KissABC.jl", push_preview = true)


text:
# Script

using Distributions, PyPlot, BayesianNonparametricStatistics

β=0.5
θ = sumoffunctions(vcat([faberschauderone],[faberschauder(j,k) for j in 0:4 for k in 1:2^j]),vcat([1.0],[(-1)^(j*k)*2^(-β*j) for j in 0:4 for k in 1:2^j]))

x = 0.0:0.001:1.0

y = θ.(x)

# Uncomment the following lines to plot θ.
# clf()
# plot(x,y)

sde = SDEWithConstantVariance(θ, 1.0, 0.0, 1000.0, 0.01)

X = rand(sde)

# Uncomment the following lines to plot a sample from sde.
# clf()
# plot(X)

M = SDEModel(1.0, 0.0)

Π = FaberSchauderExpansionWithGaussianCoefficients([2^(β*j) for j in 0:4])
postΠ = calculateposterior(Π, X, M )

for k in 1:100
    f = rand(postΠ)
    y = f.(x)
    plot(x,y)
end


text:
# module Pointwise

export pw, l, ₚ, PwVar, liftapply, PointwiseStyle
import AndTraits

"""
Pointwise application.

Pointwise function application gives meaning to expressions such as `x + y`  when `x` and `y` are functions.
That is `x + y` is the function `ω -> x(ω) + y(ω)`.

An argument can be either __lifted__ or __not lifted__.
For example in `x = 1 ~  Normal(0, 1); y = pw(+, x, 3)`, `x` will be lifted but `3` will not be in the sense that
`y` will resolve to `ω -> x(ω) + 3` and not `ω -> x(ω) + 3(ω)`.

`pw` uses some reasonable defaults for what to lift and what not to lift, but to have more explicit control use `l` and `dl` to
lift and dont lift respectively.

Example:
```
using OmegaCore
x = 1 ~ StdUniform{Float64}()
y = pw(+, x, 4)

flip(ω::Ω) = x(ω) > 0.5

f(ω::Ω) = flip(ω) ? sqrt : sin
randsample(pw(map, f, [0, 1, 2]))
randsample(pw(map, sqrt, [0, 1, 2])) # Will error!
randsample(pw(map, dl(sqrt), [0, 1, 2]))
randsample(pw(f, 3))

g(ϵ) = ω::Ω -> x(ω) + ϵ
u = 1 ~ StdNormal{Float64}()
g.(u) := ω::Ω -> g(u(ω))(ω)

g(ϵ) = (i, ω::Ω) -> i ~ Normal(ω, 0, 1) + ϵ
u = 1 ~ StdNormal{Float64}()
g.(u) := (i, ω::Ω) -> g(u(ω))(i, ω) # This is what I'd want, and the result should be a class, but i cant tell that based on
# Types of u or types of g (well, in principle I could for type of g but not in Julia)
```
"""
# 1. f(::AbstractVariableOrClass, ω) = ...
# 2. Different types, so we decide when we construct them
struct Pw{ARGS, D} <: AbstractVariableOrClass
  f::D
  args::ARGS
  Pw(f::F, args::A) where {F, A} = new{A, F}(f, args)
  Pw(f::Type{T}, args::A) where {T, A} = new{A, Type{T}}(f, args)
end

Var.traitvartype(::AbstractVariableOrClass) = TraitIsVariableOrClass
Var.traitvartype(::Type{<:AbstractVariableOrClass}) = TraitIsVariableOrClass

@inline pw(f, args...) = Pw(f, args)

@inline function (pwv::Pw)(i, ω::AbstractΩ)
  f_ = liftapply(pwv.f, i, ω)
  args_ = map(a -> liftapply(a, i, ω), pwv.args)
  ret = f_(args_...)
  liftapply(ret, i, ω)
end

@inline function Var.recurse(pwv::Pw, ω::AbstractΩ)
  f_ = liftapply(pwv.f, ω)
  args_ = map(a -> liftapply(a, ω), pwv.args)
  ret = f_(args_...)
  liftapply(ret, ω)
end
# struct PwClass{ARGS, D} <: AbstractClass
#   f::D
#   args::ARGS
#   PwClass(f::F, args::A) where {F, A} = new{A, F}(f, args)
#   PwClass(f::Type{T}, args::A) where {T, A} = new{A, Type{T}}(f, args)
# end

# (p::PwClass{Tuple{T1}})(i, ω) where {T1} =
#   lift_output(p.f(liftapply(p.args[1], i, ω)), i, ω)

# (p::PwClass{Tuple{T1, T2}})(i, ω) where {T1, T2} =
#   lift_output(p.f(liftapply(p.args[1], i, ω), liftapply(p.args[2], i, ω)), i, ω)

# struct PwVar{ARGS, D} <: AbstractVariable
#     f::D
#   args::ARGS
#   PwVar(f::F, args::A) where {F, A} = new{A, F}(f, args)
#   PwVar(f::Type{T}, args::A) where {T, A} = new{A, Type{T}}(f, args)
# end

# # FIXME: What about more than two arguments
# pw(f::F, arg1::A1) where {F, A1} =
#   inferpwtype(f, arg1, AndTraits.conjointraits(traitvartype(F), traitvartype(A1)))
# pw(f::F, arg1::A1, arg2::A2) where {F, A1, A2} =
#   inferpwtype(f, arg1, arg2, AndTraits.conjointraits(traitvartype(F), traitvartype(A1), traitvartype(A2)))

# # FIXME: Feel like this wil likely break type inference
# pw(f, args...) =
#   inferpwtype(f, args..., AndTraits.conjointraits(map(typeof, args)...))

# # inferpwtype(AndTraits.traitmatch(TraitIsVariable, TraitIsClass)) = PwClass(f, args)

# inferpwtype(f, arg1, ::AndTraits.traitmatch(TraitIsVariable, TraitIsClass)) = PwClass(f, (arg1,))
# inferpwtype(f, arg1, ::AndTraits.traitmatch(TraitIsClass)) = PwClass(f, (arg1,))
# inferpwtype(f, arg1, ::AndTraits.traitmatch(TraitIsVariable)) = PwVar(f, (arg1,))
# inferpwtype(f, arg1, ::AndTraits.traitmatch(TraitUnknownVariableType, TraitIsVariable, TraitIsClass)) = PwClass(f, (arg1,))

# inferpwtype(f, arg1, arg2, ::AndTraits.traitmatch(TraitIsVariable, TraitIsClass)) = PwClass(f, (arg1, arg2))
# inferpwtype(f, arg1, arg2, ::AndTraits.traitmatch(TraitIsClass)) = PwClass(f, (arg1, arg2))
# inferpwtype(f, arg1, arg2, ::AndTraits.traitmatch(TraitIsVariable)) = PwVar(f, (arg1, arg2))
# inferpwtype(f, arg1, arg2,  ::AndTraits.traitmatch(TraitUnknownVariableType, TraitIsVariable, TraitIsClass)) = PwClass(f, (arg1, arg2))

Base.show(io::IO, p::Pw) = print(io, p.f, "ₚ", p.args)

# Lifting
struct LiftBox{T}# <: ABox
  val::T
end
"`l(x)` constructs object that indicates that `x` should be applied pointwise.  See `pw`"
l(x) = LiftBox(x)

@inline unbox(x::LiftBox) = x.val
@inline unbox(x) = x

@inline liftapply(f::T, ω) where T = liftapplyt(traitvartype(T), f, ω)
@inline liftapply(f::T, i, ω) where T = liftapplyt(traitvartype(T), f, i, ω)

@inline liftapply(f::Ref, ω) = f[]
@inline liftapply(f::Ref, i, ω) = f[]

@inline liftapply(f::LiftBox, i, ω) = liftapply(unbox(f), i, ω)

# Random Variable
@inline liftapplyt(::AndTraits.traitmatch(TraitIsVariable), f, ω) = f(ω)
@inline liftapplyt(::AndTraits.traitmatch(TraitUnknownVariableType), f, ω) = f
@inline liftapplyt(::AndTraits.traitmatch(TraitIsVariableOrClass), f, ω) = f(ω)

# Class
@inline liftapplyt(::AndTraits.traitmatch(TraitIsVariable), f, i, ω) = f(ω)
@inline liftapplyt(::AndTraits.traitmatch(TraitIsClass), f, i, ω) = f(i, ω)
@inline liftapplyt(::AndTraits.traitmatch(TraitUnknownVariableType), f, i, ω) = f
@inline liftapplyt(::AndTraits.traitmatch(TraitIsVariableOrClass), f, i, ω) = f(i, ω)

## Broadcasting
struct PointwiseStyle <: Broadcast.BroadcastStyle end
Base.BroadcastStyle(::Type{<:AbstractVariableOrClass}) = PointwiseStyle()
Base.broadcastable(x::AbstractVariableOrClass) = x

Base.broadcasted(::PointwiseStyle, f, args...)  = pw(f, args...)
Base.BroadcastStyle(::PointwiseStyle, ::Base.Broadcast.DefaultArrayStyle{0}) = PointwiseStyle()

# Handle `f` is random variable over functions
Base.broadcast(f::AbstractVariableOrClass, args...) = pw(f, args...)
Base.broadcast(f::AbstractVariableOrClass, args::Vararg{Number}) = pw(f, args...)

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
import Base.push!
using LinearAlgebra
using Plots
using Random
using ProgressMeter

mutable struct Body
    r::Vector{Float64}
    v::Vector{Float64}
    m::Float64
end

Base.rand(rng::AbstractRNG, ::Random.SamplerType{Body}) = Body(
    2rand(Float32, 2) .- 1.,
    2rand(Float32, 2) .- 1.,
    1.)

struct Quad
    xb::Float64
    yb::Float64
    xt::Float64
    yt::Float64
end

mutable struct BHtree
    r::Vector{Float64}
    m::Float64

    quad::Quad
    nw::Union{BHtree, Body, Nothing}
    ne::Union{BHtree, Body, Nothing}
    sw::Union{BHtree, Body, Nothing}
    se::Union{BHtree, Body, Nothing}

    function BHtree(quad::Quad)
        new([0, 0],
            0,
            quad,
            nothing,
            nothing,
            nothing,
            nothing)
    end
end

BHnode = Union{BHtree, Body, Nothing}

function init_tree(quad=Quad(-10, -10, 10, 10))
    BHtree(quad)
end

function build_tree(ensemble)
    root = init_tree()
    for particle in ensemble
        try
            push!(root, particle)
        catch exc
            @debug "skipping body" exc
        end
    end
    return root
end

function balance_center_of_mass!(t::BHnode)
    if t isa BHtree
        t.m = 0
        leafs = [t.nw, t.ne, t.sw, t.se]
        for l in leafs
            if l != nothing
                t.r = (t.r * t.m + l.r * l.m) / (t.m + l.m)
                t.m = t.m + l.m
            end
        end
    end
end

function push!(t::BHtree, b::Body)

    quad_symbol = fit(t.quad, b.r)
    leaf = getproperty(t, quad_symbol)
    if leaf isa Nothing
        setproperty!(t, quad_symbol, b)
    elseif leaf isa BHtree
        push!(leaf, b)
    elseif leaf isa Body
        b2 = leaf
        child = BHtree(
            quadrant_from_symbol(t.quad, quad_symbol)
        )
        setproperty!(t, quad_symbol, child)
        push!(child, b2)
        push!(child, b)
    end

    balance_center_of_mass!(t)
end

function quadrant_from_symbol(parent::Quad, symbol::Symbol)
    mid_x = (parent.xt + parent.xb) / 2
    mid_y = (parent.yt + parent.yb) / 2
    if symbol == :ne
        Quad(mid_x, mid_y, parent.xt, parent.yt)
    elseif symbol == :se
        Quad(mid_x, parent.yb, parent.xt, mid_y)
    elseif symbol == :nw
        Quad(parent.xb, mid_y, mid_x, parent.yt)
    elseif symbol == :sw
        Quad(parent.xb, parent.yb, mid_x, mid_y)
    end
end

function fit(quad::Quad, point)
    x = point[1]
    y = point[2]
    mid_x = (quad.xt + quad.xb) / 2
    mid_y = (quad.yt + quad.yb) / 2
    if !(quad.xb <= x < quad.xt
        && quad.yb <= y < quad.yt)
        throw(DomainError(point, "outside quadrant"))
    elseif mid_x > x && mid_y <= y
        :nw
    elseif mid_x <= x && mid_y <= y
        :ne
    elseif mid_x > x && mid_y > y
        :sw
    elseif mid_x <= x && mid_y > y
        :se
    else
        error("no quad")
    end
end

function rectangle_from_coords(xb,yb,xt,yt)
    [
        xb yb
        xt yb
        xt yt
        xb yt
        xb yb
        NaN NaN
    ]
end

rect(q::Quad) = Shape([q.xb, q.xt, q.xt, q.xb], [q.yb, q.yb, q.yt, q.yt])
@recipe f(::Type{Quad}, q::Quad) = rect(q)
@recipe f(::Type{Vector{Body}}, e::Vector{Body}) = vcat(collect(map(b -> b.r, e))...)
# @recipe f(::Type{Body}, b::Body) = |> vcat()

@recipe function f(::Type{Val{:body_distribution}}, x, y, z)
    seriestype := :scatter
    @show x, y, z
    x := y[:, 1]
    y := y[:, 2]

    ()
end

function width(q::Quad)
    abs(q.xt - q.xb)
end

function plot_bhnode!(my_plt, n::BHnode)
    if n isa BHtree
        r = rect(n.quad)
        v = collect(zip(r.x, r.y))
        for (x, y) in v
            push!(my_plt, x, y)
        end
        push!(my_plt, v[1][1], v[1][2])
        push!(my_plt, NaN, NaN)
        plot_bhnode!(my_plt, n.nw)
        plot_bhnode!(my_plt, n.ne)
        plot_bhnode!(my_plt, n.sw)
        plot_bhnode!(my_plt, n.se)
    end
    my_plt
end

function plot_bodies!(my_plt, ensemble::Vector{Body})

    for (i, p) in enumerate(ensemble)
        push!(my_plt, p.r[1], p.r[2])
    end
    my_plt
end

ensemble = [
    Body([0.563299, -0.135254, 0.115279], [-0.547044, 0.256697, -0.648452], 1.0),
    Body([-0.062891, -0.501245, -0.568135], [0.816119, 0.143996, 0.753906], 1.0),
    Body([0.473176, -0.5876, 0.56664], [0.793547, 0.510543, 0.0547936], 1.0),
    Body([-0.734561, 0.354197, 0.980844], [0.674463, -0.435813, 0.622977], 1.0),
    Body([-0.401376, 0.382082, 0.86666], [-0.929288, -0.407342, 0.513163], 1.0),
    Body([0.843067, -0.417499, -0.334527], [-0.924203, -0.216044, 0.111438], 1.0),
    Body([0.445615, -0.936267, 0.417503], [0.753769, 0.551013, -0.629933], 1.0),
    Body([-0.249496, -0.753273, 0.935374], [-0.43562, -0.366953, -0.764842], 1.0),
    Body([0.150283, 0.119875, -0.659424], [0.239267, 0.316603, -0.244475], 1.0),
    Body([0.802724, -0.429913, 0.527407], [0.216259, 0.724374, 0.795849], 1.0)]

θ = 0.5
ε = 0.0001
dt = 0.0001

function rk4!(p::Body, a, h)
    k1_v = a
    k1_r = p.v

    k2_v = a .* (k1_r*h/2)
    k2_r = p.v .* k1_v*h/2

    k3_v = a .* (k2_r*h/2)
    k3_r = p.v .* k2_v*h/2

    k4_v = a .* (k3_r*h)
    k4_r = p.v .* k3_v*h

    p.r += h/6 * (k1_r + 2*k2_r + 2*k3_r + k4_r)
    p.v += h/6 * (k1_v + 2*k2_v + 2*k3_v + k4_v)
end

function leapfrog!(p::Body, a, h)
    p.r += 0.5h * p.v
    p.v += h * a
    p.r += 0.5h * p.v
end

function calculate_force(body, node::BHnode, θ=0.5)
    f = [0., 0.]
    for child in [node.nw, node.ne, node.sw, node.se]
        if child isa Body && body != child
            dr = body.r - child.r
            r2 = dot(dr, dr)
            r3 = r2 * √r2 + ε
            f -= (dr / r3) * node.m
        elseif child isa BHtree
            s = width(child.quad)
            dr = body.r - child.r
            r2 = dot(dr, dr)
            d = abs2(r2)

            if s/d < θ
                r3 = r2 * √r2 + ε
                f -= (dr / r3) * node.m
            else
                f -= calculate_force(body, child)
            end
        end
    end
    return f
end

function step!(ensemble, tree, h)
    for p1 in ensemble
        f = calculate_force(p1, tree)
        rk4!(p1, f, h)
    end
end

ensemble = rand(Body, 100)
lim = (-2, 2)
marker = (:cross, 3)

function perform!(n, ensemble::Vector{Body}, h=0.0001)
    hidden_n = 100
    output_n = Integer(n / hidden_n)

    println(stderr, "Performing simulation")
    prog = Progress(output_n, 1, "Computing...")
    tree = undef
    plt_quads = undef
    anim = Animation()
    for i=1:output_n
        for j=1:hidden_n
            tree = build_tree(ensemble)
            step!(ensemble, tree, h)
        end
        plt_quads = plot(1, framestyle=:zerolines, fill=false, c=false, linecolor=:blue)
        plot_bhnode!(plt_quads, tree)
        ehh = collect(transpose(hcat(collect(map(b -> b.r, ensemble))...)))

        plt_scatter = scatter!(ehh[:,1], ehh[:,2], m = marker,
            xlim=lim, ylim=lim,
            title="step $(i*hidden_n)")

        frame(anim, plt_scatter)
        plt_scatter
        next!(prog)
    end
    display(plt_quads)
    gif(anim)
end

# perform!(100000, ensemble)


text:
part1(filename::AbstractString) = sum(area(line) for line in eachline(filename))

function area(line::AbstractString)
    x, y, z = parse.(Int, split(line, 'x'))
    sides = [x * y, y * z, x * z]
    2 * sum(sides) + minimum(sides)
end

@assert area("2x3x4") == 58
@assert area("1x1x10") == 43

part2(filename::AbstractString) = sum(ribbonlength(line) for line in eachline(filename))

function ribbonlength(line::AbstractString)
    dims = parse.(Int, split(line, 'x'))
    2 * sum(dims) - 2 * maximum(dims) + prod(dims)
end

@assert ribbonlength("2x3x4") == 34
@assert ribbonlength("1x1x10") == 14

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
"""
    Wagner{TF} <: AbstractModel

Aerodynamic model based on Wagner's function with state variables ``\\lambda_1,
\\lambda_2``, inputs ``u, v, \\omega``, and parameters ``a, b, a_0, \\alpha_0,
c_{d0}, c_{m0}``
"""
struct Wagner{TF} <: AbstractModel
    C1::TF
    C2::TF
    eps1::TF
    eps2::TF
end

# --- Constructors --- #

"""
    Wagner(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3)

Initialize an object of type [`Wagner`](@ref)
"""
function Wagner(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3)
    return Wagner(C1, C2, eps1, eps2)
end

# --- Traits --- #

number_of_states(::Type{<:Wagner}) = 2
number_of_inputs(::Type{<:Wagner}) = 3
number_of_parameters(::Type{<:Wagner}) = 6

inplaceness(::Type{<:Wagner}) = OutOfPlace()

rate_jacobian_type(::Type{<:Wagner}) = Identity()
state_jacobian_type(::Type{<:Wagner}) = Linear()
input_jacobian_type(::Type{<:Wagner}) = Nonlinear()
parameter_jacobian_type(::Type{<:Wagner}) = Nonlinear()
time_gradient_type(::Type{<:Wagner}) = Zeros()

# --- Methods --- #

function get_residual(model::Wagner, dx, x, y, p, t)
    # extract rates
    dλ1, dλ2 = dx
    # extract states
    λ1, λ2 = x
    # extract inputs
    u, v, ω = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    C1 = model.C1
    C2 = model.C2
    ε1 = model.eps1
    ε2 = model.eps2
    # calculate rates
    return wagner_residual(dλ1, dλ2, λ1, λ2, u, v, ω, a, b, α0, C1, C2, ε1, ε2)
end

# --- Performance Overloads --- #

function get_state_jacobian(model::Wagner, dx, x, y, p, t)
    # extract inputs
    u, v, ω = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    ε1 = model.eps1
    ε2 = model.eps2
    # jacobian with respect to aerodynamic states
    return wagner_state_jacobian(u, a, b, ε1, ε2)
end

function get_input_jacobian(model::Wagner, dx, x, y, p, t)
    # extract states
    λ1, λ2 = x
    # extract inputs
    u, v, ω = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    C1 = model.C1
    C2 = model.C2
    ε1 = model.eps1
    ε2 = model.eps2
    # return jacobian
    return wagner_input_jacobian(λ1, λ2, u, v, ω, a, b, α0, C1, C2, ε1, ε2)
end

function get_parameter_jacobian(model::Wagner, dx, x, y, p, t)
    # extract states
    λ1, λ2 = x
    # extract inputs
    u, v, ω = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    C1 = model.C1
    C2 = model.C2
    ε1 = model.eps1
    ε2 = model.eps2
    # return jacobian
    return wagner_parameter_jacobian(λ1, λ2, u, v, ω, b, α0, C1, C2, ε1, ε2)
end

# --- Convenience Methods --- #

set_states!(y, model::Wagner; lambda) = y .= lambda

function set_inputs!(y, model::Wagner; u, v, omega)

    y[1] = u
    y[2] = v
    y[3] = omega

    return y
end

function set_parameters!(p, model::Wagner; a, b, a0, alpha0, cd0, cm0)

    p[1] = a
    p[2] = b
    p[3] = a0
    p[4] = alpha0
    p[5] = cd0
    p[6] = cm0

    return p
end

function separate_states(model::Wagner, x)

    return (lambda = x,)
end

function separate_inputs(model::Wagner, y)

    return (u=y[1], v=y[2], omega=y[3])
end

function separate_parameters(model::Wagner, p)

    return (a=p[1], b=p[2], a0=p[3], alpha0=p[4], cd0=p[5], cm0=p[6])
end

# --- Internal Methods --- #

function wagner_residual(dλ1, dλ2, λ1, λ2, u, v, ω, a, b, α0, C1, C2, ε1, ε2)
    r1 = dλ1 + ε1*u/b*λ1 - C1*ε1*u/b*(v + (b/2-a*b)*ω - u*α0)
    r2 = dλ2 + ε2*u/b*λ2 - C2*ε2*u/b*(v + (b/2-a*b)*ω - u*α0)
    return SVector(r1, r2)
end

wagner_state_jacobian(u, a, b, ε1, ε2) = @SMatrix [ε1*u/b 0; 0 ε2*u/b]

function wagner_input_jacobian(λ1, λ2, u, v, ω, a, b, α0, C1, C2, ε1, ε2)

    tmp1 = v/b + (1/2-a)*ω - 2*α0*u/b
    λ1dot_u = ε1/b*λ1 - C1*ε1*tmp1
    λ2dot_u = ε2/b*λ2 - C2*ε2*tmp1

    tmp2 = u/b
    λ1dot_v = -C1*ε1*tmp2
    λ2dot_v = -C2*ε2*tmp2

    tmp3 = u*(1/2-a)
    λ1dot_ω = -C1*ε1*tmp3
    λ2dot_ω = -C2*ε2*tmp3

    return @SMatrix [λ1dot_u λ1dot_v λ1dot_ω; λ2dot_u λ2dot_v λ2dot_ω]
end

function wagner_parameter_jacobian(λ1, λ2, u, v, ω, b, α0, C1, C2, ε1, ε2)

    λ1dot_a = C1*ε1*u*ω
    λ1dot_b = -ε1*u/b^2*λ1 + C1*ε1/b^2*(u*v - u^2*α0)
    λ1dot_a0 = 0
    λ1dot_α0 = C1*ε1*u^2/b
    λ1dot_cd0 = 0
    λ1dot_cm0 = 0

    λ2dot_a = C2*ε2*u*ω
    λ2dot_b = -ε2*u/b^2*λ2 + C2*ε2/b^2*(u*v - u^2*α0)
    λ2dot_a0 = 0
    λ2dot_α0 = C2*ε2*u^2/b
    λ2dot_cd0 = 0
    λ2dot_cm0 = 0

    return @SMatrix [
        λ1dot_a λ1dot_b λ1dot_a0 λ1dot_α0 λ1dot_cd0 λ1dot_cm0;
        λ2dot_a λ2dot_b λ2dot_a0 λ2dot_α0 λ2dot_cd0 λ2dot_cm0;
        ]
end

function wagner_loads(a, b, ρ, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # Wagner's function at t = 0.0
    ϕ0 = 1 - C1 - C2
    # normal force at reference point
    N = tmp1*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    # axial force at reference point
    A = -a0*ρ*b*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)^2 + ρ*b*u^2*cd0
    # moment at reference point
    M = -tmp2*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) + 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N

    return SVector(N, A, M)
end

function wagner_loads_λ(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)
    d = b/2 - a*b
    ϕ0 = 1 - C1 - C2
    N_λ = a0*ρ*u*b
    A_λ = -2*a0*ρ*b*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)
    M_λ = (b/2 + a*b)*N_λ
    return @SMatrix [N_λ N_λ; A_λ A_λ; M_λ M_λ]
end
wagner_loads_λdot() = @SMatrix [0 0; 0 0; 0 0]

wagner_loads_udot() = SVector(0, 0, 0)

function wagner_loads_vdot(a, b, ρ)
    # non-circulatory load factor
    tmp = pi*ρ*b^3
    # normal force at reference point
    N_vdot = tmp/b
    # axial force at reference point
    A_vdot = 0
    # moment at reference point
    M_vdot = -tmp/2 + (b/2 + a*b)*N_vdot

    return SVector(N_vdot, A_vdot, M_vdot)
end

function wagner_loads_ωdot(a, b, ρ)
    # non-circulatory load factor
    tmp = pi*ρ*b^3
    # normal force at reference point
    N_ωdot = -a*tmp
    # axial force at reference point
    A_ωdot = 0.0
    # moment at reference point
    M_ωdot = -tmp*(b/8 - a*b/2) + (b/2 + a*b)*N_ωdot

    return SVector(N_ωdot, A_ωdot, M_ωdot)
end

function wagner_loads_a(a, b, ρ, a0, α0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)

    tmp1 = a0*ρ*u*b
    tmp2 = pi*ρ*b^3

    ϕ0 = 1 - C1 - C2

    N = tmp1*((v + (b/2 - a*b)*ω - u*α0)*ϕ0 + λ1 + λ2) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    N_a = -tmp1*b*ω*ϕ0 - tmp2*ωdot

    A_a = 2*a0*ρ*b^2*ω*ϕ0*((v + (b/2 - a*b)*ω - u*α0)*ϕ0 + λ1 + λ2)

    M_a = tmp2*b/2*ωdot + b*N + (b/2 + a*b)*N_a

    return SVector(N_a, A_a, M_a)
end

function wagner_loads_b(a, b, ρ, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)

    tmp1 = a0*ρ*u*b
    tmp1_b = a0*ρ*u

    tmp2 = pi*ρ*b^3
    tmp2_b = 3*pi*ρ*b^2

    d = b/2 - a*b
    d_b = 1/2 - a

    ϕ0 = 1 - C1 - C2

    N = tmp1*((v + (b/2 - a*b)*ω - u*α0)*ϕ0 + λ1 + λ2) + tmp2*(vdot/b + u/b*ω - a*ωdot)

    N_b = tmp1_b*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) + tmp1*d_b*ω*ϕ0 +
        tmp2_b*(vdot/b + u/b*ω - a*ωdot) + tmp2*(-vdot/b^2 - u/b^2*ω)
    A_b = -a0*ρ*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)^2 -
        2*a0*ρ*b*d_b*ω*ϕ0*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) + ρ*u^2*cd0
    M_b = -tmp2_b*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) - tmp2*(1/8 - a/2)*ωdot +
        4*ρ*b*u^2*cm0 + (1/2 + a)*N + (b/2 + a*b)*N_b

    return SVector(N_b, A_b, M_b)
end

function wagner_loads_ρ(a, b, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)

    tmp1_ρ = a0*u*b

    tmp2_ρ = pi*b^3

    d = b/2 - a*b

    ϕ0 = 1 - C1 - C2

    N_ρ = tmp1_ρ*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) + tmp2_ρ*(vdot/b + u/b*ω - a*ωdot)
    A_ρ = -a0*b*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)^2 + b*u^2*cd0
    M_ρ = -tmp2_ρ*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) + 2*b^2*u^2*cm0 + (b/2 + a*b)*N_ρ

    return SVector(N_ρ, A_ρ, M_ρ)
end

function wagner_loads_a0(a, b, ρ, α0, C1, C2, u, v, ω, λ1, λ2)

    tmp1_a0 = ρ*u*b

    tmp2 = pi*ρ*b^3

    d = b/2 - a*b

    ϕ0 = 1 - C1 - C2

    N_a0 = tmp1_a0*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)
    A_a0 = -ρ*b*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)^2
    M_a0 = (b/2 + a*b)*N_a0

    return SVector(N_a0, A_a0, M_a0)
end

function wagner_loads_α0(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)

    tmp1 = a0*ρ*u*b

    d = b/2 - a*b

    ϕ0 = 1 - C1 - C2

    N_α0 = -tmp1*u*ϕ0
    A_α0 = 2*a0*ρ*b*u*ϕ0*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)
    M_α0 = (b/2 + a*b)*N_α0

    return SVector(N_α0, A_α0, M_α0)
end

function wagner_loads_cd0(b, ρ, u)
    N_cd0 = 0
    A_cd0 = ρ*b*u^2
    M_cd0 = 0
    return SVector(N_cd0, A_cd0, M_cd0)
end

function wagner_loads_cm0(b, ρ, u)
    N_cm0 = 0
    A_cm0 = 0
    M_cm0 = 2*ρ*b^2*u^2
    return SVector(N_cm0, A_cm0, M_cm0)
end

function wagner_loads_u(a, b, ρ, a0, α0, cd0, cm0, C1, C2, u, v, ω, λ1, λ2)
    tmp1 = a0*ρ*u*b
    tmp1_u = a0*ρ*b
    tmp2 = pi*ρ*b^3
    d = b/2 - a*b
    ϕ0 = 1 - C1 - C2
    N_u = tmp1_u*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) - tmp1*α0*ϕ0 + tmp2/b*ω
    A_u = 2*a0*ρ*b*α0*ϕ0*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2) + 2*ρ*b*u*cd0
    M_u = -tmp2*ω + 4*ρ*b^2*u*cm0 + (b/2 + a*b)*N_u
    return SVector(N_u, A_u, M_u)
end

function wagner_loads_v(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)
    ϕ0 = 1 - C1 - C2
    d = b/2 - a*b
    N_v = a0*ρ*u*b*ϕ0
    A_v = -2*a0*ρ*b*ϕ0*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)
    M_v = (b/2 + a*b)*N_v

    return SVector(N_v, A_v, M_v)
end

function wagner_loads_ω(a, b, ρ, a0, α0, C1, C2, u, v, ω, λ1, λ2)
    tmp1 = a0*ρ*u*b
    tmp2 = pi*ρ*b^3
    d = b/2 - a*b
    ϕ0 = 1 - C1 - C2
    N_ω = tmp1*d*ϕ0 + tmp2*u/b
    A_ω = -2*a0*ρ*b*d*ϕ0*((v + d*ω - u*α0)*ϕ0 + λ1 + λ2)
    M_ω = -tmp2*u + (b/2 + a*b)*N_ω

    return SVector(N_ω, A_ω, M_ω)
end


text:
macro incrmark!(x)
    quote
        $(esc(x)) = incrmark($(esc(x)))
    end
end

macro decrmark!(x)
    quote
        $(esc(x)) = decrmark($(esc(x)))
    end
end

macro setmark!(x,m)
    quote
        $(esc(x)) = setmark($(esc(x)),$(esc(m)))
    end
end

macro unmark!(x)
    quote
        $(esc(x)) = unmark($(esc(x)))
    end
end


text:
@testset "Testing Vec Constructors" begin
	@test Vec{2,ComplexF64}(1+im,2) === SVector{2,ComplexF64}(1+im,2)
	@test Vec3(1,2,3) === SVector{3,Float}(1,2,3)
	@test Vec2(1,2) === SVector{2,Float}(1,2)
	@test Vec1(1) === SVector{1,Float}(1)
end


text:
operators = [:(+), :(-), :(*), :(/)]

for operator in operators
    f = Meta.parse("Base.:$(operator)")
    elementwise_operator = Meta.parse(".$(operator)")

    funcs = quote
        """
            $($f)(a::CertainValue, b::AbstractUncertainValue; n::Int = 30000)

        Compute `a $($operator) b`. Treats the certain value as a scalar and performs the
        operation element-wise on a default of `n = 30000` realizations of `b`.

        To tune the number of draws to `n`, use the `$($operator)(a, b, n::Int)` syntax.
        """
        function $(f)(a::CertainValue, b::AbstractUncertainValue; n::Int = 30000)
            $(elementwise_operator)(a.value, b, n)
        end

        """
            $($f)(a::AbstractUncertainValue, b::CertainValue; n::Int = 30000)

        Compute `a $($operator) b`. Treats the certain value as a scalar and performs the
        operation element-wise on a default of `n = 30000` realizations of `a`.

        To tune the number of draws to `n`, use the `$($operator)(a, b, n::Int)` syntax.
        """
        function $(f)(a::AbstractUncertainValue, b::CertainValue; n::Int = 30000)
            $(elementwise_operator)(a, b.value, n)
        end

        """
        $($f)(a::AbstractUncertainValue, b::CertainValue; n::Int)

        Compute `a $($operator) b`. Treats the certain value as a scalar and performs the
        operation element-wise on `n` realizations of `a`.

        This function is called with the `$($operator)(a, b, n::Int)` syntax.
            """
        function $(f)(a::AbstractUncertainValue, b::CertainValue, n::Int)
            $(elementwise_operator)(a, b.value, n)
        end

        """
            $($f)(a::CertainValue, b::AbstractUncertainValue, n::Int)

        Compute `a $($operator) b`. Treats the certain value as a scalar and performs the
        operation element-wise on `n` realizations of `b`.

        This function is called with the `$($operator)(a, b, n::Int)` syntax.
        """
        function $(f)(a::CertainValue, b::AbstractUncertainValue, n::Int)
            $(elementwise_operator)(a.value, b, n)
        end
    end

    eval(funcs)
end


text:
# from computations/differentiation.jl
export differentiation,
    antidifferentiation,
    derivative_dict,
    antiderivative_dict

"""
The differentation function returns an operator that can be used to differentiate
a function in the dictionary, with the result as an expansion in a second dictionary.

The differentiation operator is efficient for some combinations of source and
destination dictionary. If a dictionary supports at least one differentiation
operator, than `hasderivative(dict)` is true and `derivative_dict(dict)` returns
the destination dictionary.

The order of the differentation can be passed as an optional second argument.

Examples are:
```
differentiation(Φ::Dictionary)
```
This will use the default destination dictionary, order 1 and dimension 1.
Another example is:
```
differentiation(Φ, 1)
```
or
```
differentiation(Φ, (0,1))
```
"""
differentiation

"""
The antidifferentiation function returns an operator that can be used to find the antiderivative
of a function in the dictionary, with the result an expansion in a second dictionary.

See also: `differentiation`.
"""
antidifferentiation

orderiszero(order::Int) = order==0
orderiszero(order) = sum(order) == 0

"Determine the order of differentiation from the arguments given by the user."
difforder(Φ::Dictionary; order = nothing, dim = nothing, options...) =
    _difforder(Φ, order, dim)

# Integer order is given: we use that
_difforder(Φ::Dictionary, order::Int, dim::Nothing) = order
# Integer order and dimension are given: we make a tuple of the form (0,0,order,0)
_difforder(Φ::Dictionary, order::Int, dim::Int) = order .* dimension_tuple(dimension(Φ), dim)
# Nothing is given: assume integer order 1 and try again
_difforder(Φ::Dictionary, order::Nothing, dim::Nothing) =
    dimension(Φ) == 1 ? _difforder(Φ, 1, dim) : _difforder(Φ, 1, 1)
# Only dimension is given: use order 1 in that dimension.
_difforder(Φ::Dictionary, order::Nothing, dim::Int) = _difforder(Φ, 1, dim)
# None of the above: we don't know what to do.
_difforder(Φ::Dictionary, order, dim) = error("Supplied derivative order and dimension not understood.")

# Assign a default order if none is given
"A dictionary in the same span as the differentiated dictionary"
derivative_dict(Φ::Dictionary; options...) =
    derivative_dict(Φ, difforder(Φ; options...); options...)
antiderivative_dict(Φ::Dictionary; options...) =
    antiderivative_dict(Φ, difforder(Φ; options...); options...)
"""
    A dictionary that evaluates in the same way as the diffentiated dictionary,
    but constructed using the `derivative_dict` and a transformation matrix that maps the coefficients.
"""
function differentiation_dict(dict::Dictionary, args...; options...)
	@assert hasderivative(dict)
	differentiation(dict, args...; options...) * dict
end

# If a set has a differentiation operator, then we can represent the set of derivatives
# by an OperatedDict.
# @deprecate derivative(dict::Dictionary, args...; options...) diff(dict, args...; options...)
"""
The diffentiated dictionary
"""
diff(dict::Dictionary; options...) =
    diff(dict, difforder(dict; options...); options...)
function diff(dict::Dictionary, args...; options...)
	if hasderivative(dict)
        differentiation_dict(dict, args...; options...)
    else
        differentiation_dict(dict, args...; options...)
        # error("Default not yet implemented")
    end
end

# Insert the operator eltype if it is not given
differentiation(dict::Dictionary, args...; options...) =
    differentiation(operatoreltype(dict), dict, args...; options...)
differentiation(src::Dictionary, dest::Dictionary, args...; options...) =
    differentiation(operatoreltype(src, dest), src, dest, args...; options...)

antidifferentiation(dict::Dictionary, args...; options...) =
    antidifferentiation(operatoreltype(dict), dict, args...; options...)
antidifferentiation(src::Dictionary, dest::Dictionary, args...; options...) =
    antidifferentiation(operatoreltype(src, dest), src, dest, args...; options...)

# Assign a default order if none is given
differentiation(::Type{T}, src::Dictionary; options...) where {T} =
    differentiation(T, src, difforder(src; options...); options...)
differentiation(::Type{T}, src::Dictionary, dest::Dictionary; options...) where {T} =
    differentiation(T, src, dest, difforder(src; options...); options...)

antidifferentiation(::Type{T}, src::Dictionary; options...) where {T} =
    antidifferentiation(T, src, difforder(src; options...); options...)
antidifferentiation(::Type{T}, src::Dictionary, dest::Dictionary; options...) where {T} =
    antidifferentiation(T, src, dest, difforder(src; options...); options...)

function differentiation(::Type{T}, src::Dictionary, order::Int; options...) where {T}
    if orderiszero(order)
        IdentityOperator{T}(src)
    else
        if hasderivative(src, order)
            dest = derivative_dict(src, order; options...)
            differentiation(T, src, dest, order; options...)
        else
            # If order==1 and the test above failed, there is an issue with the dictionary.
            order == 1 && error("Dictionary does not support differentiation.")
            # From here on, we assume that order > 1
            @assert order > 1
            dest1 = derivative_dict(src, 1; options...)
            D1 = differentiation(T, src, dest1, 1; options...)
            D2 = differentiation(T, dest1, order-1; options...)
            D2*D1
        end
    end
end

function differentiation(::Type{T}, src::Dictionary, order; options...) where {T}
    if orderiszero(order)
        IdentityOperator{T}(src)
    else
        @assert hasderivative(src, order)
        dest = derivative_dict(src, order; options...)
        differentiation(T, src, dest, order; options...)
    end
end

function antidifferentiation(::Type{T}, src::Dictionary, order::Int; options...) where {T}
    if orderiszero(order)
        IdentityOperator{T}(src)
    else
        if hasantiderivative(src, order)
            dest = antiderivative_dict(src, order; options...)
            antidifferentiation(T, src, dest, order; options...)
        else
            # If order==1 and the test above failed, there is an issue with the dictionary.
            order == 1 && error("Dictionary does not support antidifferentiation.")
            # From here on, we assume that order > 1
            @assert order > 1
            dest1 = antiderivative_dict(src, 1; options...)
            D1 = antidifferentiation(T, src, dest1, 1; options...)
            D2 = antidifferentiation(T, dest1, order-1; options...)
            D2*D1
        end
    end
end

function antidifferentiation(::Type{T}, src::Dictionary, order; options...) where {T}
    if orderiszero(order)
        IdentityOperator{T}(src)
    else
        @assert hasantiderivative(src, order)
        dest = antiderivative_dict(src, order; options...)
        antidifferentiation(T, src, dest, order; options...)
    end
end


text:
using GreatCircle
using Geodesy: LatLon
using Base.Test

# One test decimal degree is 111000m
latitude = 40.0
longitude = -76.0
position = LatLon(lat=latitude, lon=longitude)

"""
    Great Circle Tests
"""
azimuth = 90.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone to the right
@test new_gc["longitude"] > longitude + 0.9

a = great_circle(111320., azimuth, 0., 0.)
@test abs(a["latitude"]-0.)<=eps(typeof(a["latitude"]))
for i = 1:2
    a = great_circle(111320., azimuth, a["latitude"], a["longitude"])
    @test abs(a["latitude"]-0.)<=eps(typeof(a["latitude"]))
end
aposition = LatLon(0., 0.)
a = great_circle(111320., azimuth, aposition)
@test abs(a.lat-0.)<=eps(typeof(a.lat))
for i = 1:2
    a = great_circle(111320., azimuth, a)
    @test abs(a.lat-0.)<=eps(typeof(a.lat))
end

azimuth = 270.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone to the left
@test new_gc["longitude"] < longitude - 0.9

azimuth = 180.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone down
@test new_gc["latitude"] < latitude - 0.9

azimuth = 0.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up
@test new_gc["latitude"] > latitude + 0.9

azimuth = 315.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up and to the left
@test new_gc["latitude"] > latitude + 0.45
@test new_gc["longitude"] < longitude - 0.45

new_gc = great_circle(111000., azimuth, position)
@test new_gc.lat > latitude + 0.45
@test new_gc.lon < longitude - 0.45

# One decimal degree is 111000m
latitude = collect(40.0:5.:60.0)
longitude = collect(-40.0:-5.:-60.0)

azimuth = ones(longitude) * 90.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone to the right
@test all(new_gc["longitude"] .> longitude + 0.9)

azimuth = 270.
distance = 111000. * ones(longitude)
new_gc = great_circle(distance, azimuth, latitude, longitude)
# We should have gone to the left
@test all(new_gc["longitude"] .< longitude - 0.9)

azimuth = 180.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone down
@test all(new_gc["latitude"] .< latitude - 0.9)

azimuth = 0.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up
@test all(new_gc["latitude"] .> latitude + 0.9)

azimuth = 315.
new_gc = great_circle(111000., azimuth, latitude, longitude)
# We should have gone up and to the left
@test all(new_gc["latitude"] .> latitude + 0.45)
@test all(new_gc["longitude"] .< longitude - 0.45)

"""
    Great Distance Tests
"""
# One decimal degree at the equator is about 111.32km
latitude_start  = 0.
latitude_end    = 0.
longitude_start = 50.
longitude_end   = 52.
gd = great_distance(latitude_start, longitude_start, latitude_end, longitude_end)
@test round(gd["distance"] ./ 1000, 2) .== 111.32 .* 2

a = LatLon(lat=latitude_start, lon=longitude_start)
b = LatLon(lat=latitude_end, lon=longitude_end)
gd = great_distance(a, b)
@test round(gd["distance"] ./ 1000, 2) .== 111.32 .* 2

# One decimal degree is 111000m
latitude_start  = 0.
latitude_end    = 0.
longitude_start = [49., 75.]
longitude_end   = [50., 76.]
gd = great_distance(latitude_start, longitude_start, latitude_end, longitude_end)
@test all( round(gd["distance"] ./ 1000, 2) .== 111.32 )


text:
function Reagents.channel(::Type{A} = Any, ::Type{B} = A) where {A,B}
    bag1 = Bag{Message{A}}()
    bag2 = Bag{Message{B}}()
    return Swap(bag1, bag2), Swap(bag2, bag1)
end
# TODO: is it better to return `Endpoint` which then lifted to `Swap`?

# mutable struct Message{A,B}
mutable struct Message{A}
    payload::A
    reaction::Reaction
    continuation::Reactable
    offer::Offer
end

function isdeleted(msg::Message)
    offer = msg.offer
    # TODO: don't wait?
    return !(offer.state[] isa Union{Pending,Waiting})
end

struct Swap{A,B} <: Reagent
    msgs::Bag{Message{A}}
    dual::Bag{Message{B}}
end

hascas(::Swap) = true
maysync(::Swap) = true

function tryreact!(
    actr::Reactor{<:Swap{A,B}},
    a,
    rx::Reaction,
    offer::Union{Offer,Nothing},
) where {A,B}
    a = convert(A, a)
    (; msgs, dual) = actr.reagent
    if offer !== nothing  # && !maysync(actr.continuation)
        let m = Message{A}(a, rx, actr.continuation, offer)
            push!(msgs, m)
            @trace(
                label = :pushed_offer,
                offerid = offerid(offer),
                taskid = objectid(current_task()),
                still_in_msg = m in msgs,
                offer,
                msgs,
            )
            # Add backedges for cleaning up catalysts. Block the reaction if the
            # catalyst is already removed.
            backedge!(offer, msgs) && return Block()
        end
    else
        @trace(
            label = :swap_without_offer,
            offerid = offerid(offer),
            taskid = objectid(current_task()),
            offer,
            msgs,
        )
    end
    retry = false
    for msg in dual
        if msg.offer === offer || has(rx.offers, msg.offer)
            continue
        end
        ans = tryreact_together!(msg, actr.continuation, a, rx, offer)
        @trace(
            label = :done_tryreact_together,
            offerid = offerid(offer),
            taskid = objectid(current_task()),
            msg,
            ans,
            a,
            rx,
            offer,
        )
        if ans isa Retry
            maybe_redissolve!(msg)
            retry = true
        elseif ans isa SomehowBlocked
            # continue
        else
            return ans
        end
        let dual_offer = msg.offer
            if dual_offer isa Waiter
                if msg.reaction.restart_on_failure
                    tryput!(dual_offer, Rescinded())
                end
            end
        end
    end
    for _ in msgs  # TODO: better cleanup
    end
    return retry ? Retry() : Block()
end

function tryreact_together!(msg::Message, k::Reactable, a, rx, offer)
    dual_offer = msg.offer
    function commit_dual_offer(b)  # `b` is the output of `msg.continuation`
        # TODO: don't wait for `CASing` here?
        old = dual_offer.state[]
        @trace(
            label = :commit_dual_offer,
            offer_state = old,
            taskid = objectid(current_task()),
            dual_taskid = objectid(dual_offer.task),
            offerid = offerid(dual_offer),
        )
        old isa Union{Waiting,Pending} || return Retry()
        function fulfil_dual_offer(old)
            function wake_dual(_)
                @trace(
                    label = :fulfil_wake_dual,
                    taskid = objectid(current_task()),
                    dual_taskid = objectid(dual_offer.task),
                    offerid = offerid(dual_offer),
                )
                schedule(dual_offer.task)
            end
            if old isa Pending
                return CAS(dual_offer.state, old, b)
            elseif old isa Waiting
                return CAS(dual_offer.state, old, b) ⨟ PostCommit(wake_dual)
            else
                return Retry()
            end
        end
        return Read(dual_offer.state) ⨟ Computed(fulfil_dual_offer)
    end
    actr = then(
        Reagent(msg.continuation) ⨟  # execute the dual continuation with `a` as the input
        Computed(commit_dual_offer) ⨟  # use `CAS` reagent to set/commit the dual offer
        Return(msg.payload),  # use the dual's value as the input for this continuation `k`
        k,
    )
    rx2 = withoffer(combine(rx, msg.reaction), msg.offer)
    rx2 = @set rx2.restart_on_failure = rx.restart_on_failure
    return tryreact!(actr, a, rx2, offer)
end

function Base.show(io::IO, ::MIME"text/plain", @nospecialize(swap::Swap{A,B})) where {A,B}
    nsenders = nitems(swap.msgs)
    nreceivers = nitems(swap.dual)
    print(io, "<Swap ", A, " → ", B, ": ")
    print(io, nsenders, " sender(s) ")
    print(io, nreceivers, " receiver(s)")
    print(io, '>')
end


text:
isbanded(A::AbstractTriangular) = isbanded(parent(A))
bandwidths(A::Union{UpperTriangular,UnitUpperTriangular}) =
    (min(0,bandwidth(parent(A),1)), bandwidth(parent(A),2))
bandwidths(A::Union{LowerTriangular,UnitLowerTriangular}) =
    (bandwidth(parent(A),1), min(0,bandwidth(parent(A),2)))

triangularlayout(::Type{Tri}, ML::BandedColumns) where {Tri} = Tri(ML)
triangularlayout(::Type{Tri}, ML::BandedRows) where {Tri} = Tri(ML)
triangularlayout(::Type{Tri}, ML::ConjLayout{<:BandedRows}) where {Tri} = Tri(ML)

function tribandeddata(::TriangularLayout{'U'}, A)
    B = triangulardata(A)
    u = bandwidth(B,2)
    D = bandeddata(B)
    view(D, 1:u+1, :)
end

function tribandeddata(::TriangularLayout{'L'}, A)
    B = triangulardata(A)
    l,u = bandwidths(B)
    D = bandeddata(B)
    view(D, u+1:l+u+1, :)
end

tribandeddata(A) = tribandeddata(MemoryLayout(A), A)

Base.replace_in_print_matrix(A::Union{UpperTriangular{<:Any,<:AbstractBandedMatrix},
                                      UnitUpperTriangular{<:Any,<:AbstractBandedMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

Base.replace_in_print_matrix(A::Union{LowerTriangular{<:Any,<:AbstractBandedMatrix},
                                      UnitLowerTriangular{<:Any,<:AbstractBandedMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    -bandwidth(A,1) ≤ j-i ≤ bandwidth(A,2) ? s : Base.replace_with_centered_mark(s)

# Mul
@lazylmul UpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul UnitUpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul LowerTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazylmul UnitLowerTriangular{T, <:AbstractBandedMatrix{T}} where T

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{'U',UNIT,<:BandedColumnMajor},
                                   <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbmv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{'L',UNIT,<:BandedColumnMajor},
                                   <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbmv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                   <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbmv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatMulVec{<:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                   <:AbstractStridedLayout, T, T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbmv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
end

# Ldiv
@lazyldiv UpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv UnitUpperTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv LowerTriangular{T, <:AbstractBandedMatrix{T}} where T
@lazyldiv UnitLowerTriangular{T, <:AbstractBandedMatrix{T}} where T

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector{T},
         M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:BandedColumnMajor},
                                   <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbsv!('U', 'N', UNIT, size(A,1), bandwidth(A,2), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:BandedColumnMajor},
                                   <:AbstractStridedLayout,T,T}) where {UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbsv!('L', 'N', UNIT, size(A,1), bandwidth(A,1), tribandeddata(A), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,BandedRowMajor},
                                   <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbsv!(UPLO, 'T', UNIT, transpose(tribandeddata(A)), dest)
end

@inline function _copyto!(::AbstractStridedLayout, dest::AbstractVector,
         M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,ConjLayout{BandedRowMajor}},
                                   <:AbstractStridedLayout,T,T}) where {UPLO,UNIT,T <: BlasFloat}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    tbsv!(UPLO, 'C', UNIT, tribandeddata(A)', dest)
end

## generic fallback

#Generic solver using naive substitution
# manually hoisting x[j] significantly improves performance as of Dec 2015
# manually eliding bounds checking significantly improves performance as of Dec 2015
# directly indexing A.data rather than A significantly improves performance as of Dec 2015
# replacing repeated references to A.data with [Adata = A.data and references to Adata]
# does not significantly impact performance as of Dec 2015
# replacing repeated references to A.data[j,j] with [Ajj = A.data[j,j] and references to Ajj]
# does not significantly impact performance as of Dec 2015
function banded_naivesub!(::TriangularLayout{'U','N'}, A, b::AbstractVector, x::AbstractVector = b)
    require_one_based_indexing(A, b, x)
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    D = triangulardata(A)
    u = bandwidth(D,2)
    @inbounds for j in n:-1:1
        iszero(D[j,j]) && throw(SingularException(j))
        xj = x[j] = D[j,j] \ b[j]
        for i in j-1:-1:max(1,j-u) # counterintuitively 1:j-1 performs slightly better
            b[i] -= D[i,j] * xj
        end
    end
    x
end
function banded_naivesub!(::TriangularLayout{'U','U'}, A, b::AbstractVector, x::AbstractVector = b)
    require_one_based_indexing(A, b, x)
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    D = triangulardata(A)
    u = bandwidth(D,2)
    @inbounds for j in n:-1:1
        xj = x[j] = b[j]
        for i in j-1:-1:max(1,j-u) # counterintuitively 1:j-1 performs slightly better
            b[i] -= D[i,j] * xj
        end
    end
    x
end
function banded_naivesub!(::TriangularLayout{'L','N'}, A, b::AbstractVector, x::AbstractVector = b)
    require_one_based_indexing(A, b, x)
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    D = triangulardata(A)
    l = bandwidth(D,1)
    @inbounds for j in 1:n
        iszero(D[j,j]) && throw(SingularException(j))
        xj = x[j] = D[j,j] \ b[j]
        for i in j+1:min(j+l,n)
            b[i] -= D[i,j] * xj
        end
    end
    x
end
function banded_naivesub!(::TriangularLayout{'L','U'}, A, b::AbstractVector, x::AbstractVector = b)
    require_one_based_indexing(A, b, x)
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    D = triangulardata(A)
    l = bandwidth(D,1)
    @inbounds for j in 1:n
        xj = x[j] = b[j]
        for i in j+1:min(j+l,n)
            b[i] -= D[i,j] * xj
        end
    end
    x
end

function _copyto!(_, dest::AbstractVector,
                M::MatLdivVec{<:TriangularLayout{UPLO,UNIT,<:AbstractBandedLayout}}) where {UPLO,UNIT}
    A,x = M.args
    x ≡ dest || copyto!(dest, x)
    banded_naivesub!(MemoryLayout(A), A, x)
end


text:
#=

Filename: Model.jl.
Author: Shunsuke Hori

Defining an instance "Model" storing parameters, functional forms, and grid.

=#

struct Model{TF <: AbstractFloat, TR <: Real, TI <: Integer}
    alpha::TR    # Productivity parameter
    beta::TF     # Discount factor
    gamma::TR    # risk aversion
    mu::TR       # First parameter in lognorm(mu, sigma)
    s::TR        # Second parameter in lognorm(mu, sigma)
    grid_min::TR # Smallest grid point
    grid_max::TR # Largest grid point
    grid_size::TI # Number of grid points
    u::Function       # utility function
    u_prime::Function # derivative of utility function
    f::Function       # production function
    f_prime::Function # derivative of production function
    grid::Vector{TR}   # grid
end

"""
construct Model instance using the information of parameters and functional form

arguments: see above

return: Model type instance
"""
function Model(; alpha::Real=0.65,   # Productivity parameter
                 beta::AbstractFloat=0.95,    # Discount factor
                 gamma::Real=1.0,    # risk aversion
                 mu::Real=0.0,       # First parameter in lognorm(mu, sigma)
                 s::Real=0.1,        # Second parameter in lognorm(mu, sigma)
                 grid_min::Real=1e-6,# Smallest grid point
                 grid_max::Real=4.0, # Largest grid point
                 grid_size::Integer=200, # Number of grid points
                 u::Function= c->(c^(1-gamma)-1)/(1-gamma), # utility function
                 u_prime::Function = c-> c^(-gamma),        # u'
                 f::Function = k-> k^alpha,                 # production function
                 f_prime::Function = k -> alpha*k^(alpha-1) # f'
                 )

    grid=collect(linspace(grid_min, grid_max, grid_size))

    if gamma == 1   # when gamma==1, log utility is assigned
        u_log(c) = log(c)
        m = Model(alpha, beta, gamma, mu, s, grid_min, grid_max,
                grid_size, u_log, u_prime, f, f_prime, grid)
    else
        m = Model(alpha, beta, gamma, mu, s, grid_min, grid_max,
                  grid_size, u, u_prime, f, f_prime, grid)
    end
    return m
end


text:
using SparseArrays, LinearAlgebra
using JLD, CSV, DataFrames, DataFramesMeta
#using DataDeps
using ArtifactUtils, Pkg.Artifacts

artifacts_path = joinpath(pkgdir(VarianceComponentsHDFE), "benchmark\\Artifacts.toml")

# add in code to generate .jld, .csv, etc. in the benchmark/data directory
# This should only generate if the file doesn't exist, or if this `force_generate = true`
force_generate = false

pkg_dir = pkgdir(VarianceComponentsHDFE)

# Obtain environment variables
run_large_benchmark = get(ENV, "VCHDFE_LARGE_BENCHMARK", "false")  == "true" ? true : false
run_huge_benchmark = get(ENV, "VCHDFE_HUGE_BENCHMARK", "false")  == "true" ? true : false

## Medium-Sized Network Generator
function rademacher!(R; demean = false)
    R .= R - (R .== 0)

    if demean == true
        R .= R .- mean(R, dims = 2)
    end
    return nothing
end

function compute_X_No_Controls(data)
    id = data.id
    firmid = data.firmid
    y = data.y

    NT = size(y,1);
    J = maximum(firmid);
    N = maximum(id);
    K = 0
    nparameters = N + J + K

    #Worker Dummies
    D = sparse(collect(1:NT),id,1);

    #Firm Dummies
    F = sparse(collect(1:NT),firmid,1);

    # N+J x N+J-1 restriction matrix
    S= sparse(1.0I, J-1, J-1);
    S=vcat(S,sparse(-zeros(1,J-1)));
    X_Laplacian = hcat(D, -F)
    X_GroundedLaplacian = hcat(D, -F*S)
    S_xx = Symmetric(X_GroundedLaplacian'*X_GroundedLaplacian)
    S_xx_lap = Symmetric(SparseMatrixCSC{Float64,Int64}(X_Laplacian'*X_Laplacian))

    A, diag = adj(SparseMatrixCSC{Float64,Int64}(X_Laplacian'*X_Laplacian))

    return X_Laplacian, X_GroundedLaplacian, S_xx, S_xx_lap, A
end

function compute_X_Controls(data)
    id = data.id
    firmid = data.firmid
    y = data.y

    NT=size(y,1);
    J=maximum(firmid);
    N=maximum(id);
    K = 2
    nparameters = N + J + K

    #Worker Dummies
    D=sparse(collect(1:NT),id,1);

    #Firm Dummies
    F=sparse(collect(1:NT),firmid,1);

    # N+J x N+J-1 restriction matrix
    S= sparse(1.0I, J-1, J-1);
    S=vcat(S,sparse(-zeros(1,J-1)));

    #Assuming columns 5 and 6 are the controls
    controls = hcat(data.control1[id], data.control2[id])

    Xcontrols = hcat(D,F*S,controls)
    S_xx = Symmetric(Xcontrols'*Xcontrols)

    return Xcontrols, S_xx
end

# We don't use the symmetric flag here, it is making the cluster kill the process
# Also flagged X_laplacian as a sparse matrix of Floats
function compute_X_No_Controls_Huge(data)
    id = data.id
    firmid = data.firmid
    y = data.y

    NT = size(y,1);
    J = maximum(firmid);
    N = maximum(id);
    K = 0
    nparameters = N + J + K

    #Worker Dummies
    D = sparse(collect(1:NT),id,1);

    #Firm Dummies
    F = sparse(collect(1:NT),firmid,1);

    # N+J x N+J-1 restriction matrix
    S= sparse(1.0I, J-1, J-1);
    S=vcat(S,sparse(-zeros(1,J-1)));
    X_Laplacian = SparseMatrixCSC{Float64,Int64}(hcat(D, -F))
    X_GroundedLaplacian = hcat(D, -F*S)
    S_xx = X_GroundedLaplacian'*X_GroundedLaplacian
    S_xx_lap = X_Laplacian'*X_Laplacian
    A, diag = adj(S_xx_lap)

    return X_Laplacian, X_GroundedLaplacian, S_xx, S_xx_lap, A
end

if ~isfile(pkg_dir*"/benchmark/data/medium_main.jld") || force_generate
    #data = CSV.read(datadep"VarianceComponentsHDFE/medium_nocontrols_pruned.csv"; header=true)
    data_hash = artifact_hash("medium_nocontrols_pruned",artifacts_path)
    data = CSV.read(joinpath(artifact_path(data_hash), "medium_nocontrols_pruned.csv"))

    X_Laplacian, X_GroundedLaplacian, S_xx , S_xx_lap, A = compute_X_No_Controls(data)
    save(pkg_dir*"/benchmark/data/medium_main.jld", "X_Laplacian", X_Laplacian, "X_GroundedLaplacian", X_GroundedLaplacian, "S_xx", S_xx, "S_xx_lap", S_xx_lap, "A", A)
end

if ~isfile(pkg_dir*"/benchmark/data/medium_controls_main.jld") || force_generate
    #data = CSV.read(datadep"VarianceComponentsHDFE/medium_controls_pruned.csv"; header=true)
    data_hash = artifact_hash("medium_controls_pruned",artifacts_path)
    data = CSV.read(joinpath(artifact_path(data_hash), "medium_controls_pruned.csv"))

    Xcontrols, S_xx = compute_X_Controls(data)
    save(pkg_dir*"/benchmark/data/medium_controls_main.jld", "Xcontrols", Xcontrols, "S_xx", S_xx)
end

if run_large_benchmark && (~isfile(pkg_dir*"/benchmark/data/large_main.jld") || force_generate)
    #data = CSV.read(datadep"VarianceComponentsHDFE/large_nocontrols_pruned.csv"; header=true)
    data_hash = artifact_hash("large_nocontrols_pruned",artifacts_path)
    data = CSV.read(joinpath(artifact_path(data_hash), "large_nocontrols_pruned.csv"))

    X_Laplacian, X_GroundedLaplacian, S_xx , S_xx_lap, A = compute_X_No_Controls(data)
    save(pkg_dir*"/benchmark/data/large_main.jld", "X_Laplacian", X_Laplacian, "X_GroundedLaplacian", X_GroundedLaplacian, "S_xx", S_xx, "S_xx_lap", S_xx_lap, "A", A)
 end

 if run_large_benchmark && (~isfile(pkg_dir*"/benchmark/data/large_controls_main.jld") || force_generate)
    #lar_csv = CSV.read(joinpath(artifact"large_nocontrols_pruned", "large_nocontrols_pruned.csv"))
    data_hash = artifact_hash("large_controls_pruned",artifacts_path)
    data = CSV.read(joinpath(artifact_path(data_hash), "large_controls_pruned.csv"))

    Xcontrols, S_xx = compute_X_Controls(data)
    save(pkg_dir*"/benchmark/data/large_controls_main.jld", "Xcontrols", Xcontrols, "S_xx", S_xx)
 end

 if run_huge_benchmark && (~isfile(pkg_dir*"/benchmark/data/huge_main.jld") || force_generate)
     #data = CSV.read(datadep"VarianceComponentsHDFE/huge_pruned_main.csv"; header=true)
     data_hash = artifact_hash("huge_nocontrols_pruned",artifacts_path)
     data = CSV.read(joinpath(artifact_path(data_hash), "huge_nocontrols_pruned.csv"))

     X_Laplacian, X_GroundedLaplacian, S_xx , S_xx_lap, A = compute_X_No_Controls_Huge(data)
     save(pkg_dir*"/benchmark/data/huge_main.jld", "X_Laplacian", X_Laplacian, "X_GroundedLaplacian", X_GroundedLaplacian, "S_xx", S_xx, "S_xx_lap", S_xx_lap, "A", A)
 end


text:
alg_maximum_order(alg::AbstractMethodOfStepsAlgorithm) = alg_maximum_order(alg.alg)

isconstrained(alg::AbstractMethodOfStepsAlgorithm{constrained}) where constrained =
    constrained