text:
export NormalNormalTransformer

@doc raw"""
    NormalNormalTransformer(normal_sampler, new_params, old_params=(0.0, 1.0))

Generate normal random variables from another normal sampler.

Transform samples from a normal distribution sampler into samples from a normal
distribution with different parameters. This uses the result that if ``Z`` is a
standard normal random variable (``Z \sim \text{Normal}(0, 1)``) then
``X =  \sigma Z + \mu`` has distribution ``\text{Normal}(\mu, \sigma^2)``.

# Arguments
* `normal_sampler::sampler`: a sampler producing normally distributed samples.
* `new_params::Tuple{Real, Real}`: a tuple (mu, sigma2) giving the mean and
    variance of the target distribution
* `old_params::Tuple{Real, Real}=(0.0, 1.0)`: a tuple (mu, sigma2) giving the
    mean and variance of the normal sampler. Defaults to a standard normal
    distribution.

# Examples
```julia
# Transforming ten standard normal random variables into N(1, 2)
s = BoxMullerTransformSampler()
t = NormalNormalTransformer(s, (1, 2))
transformed_samples = sample(t, 10)
```
"""
struct NormalNormalTransformer <: Sampler
    normal_sampler::Sampler
    new_params::NamedTuple{(:mu, :sigma2), Tuple{Real, Real}}
    old_params::NamedTuple{(:mu, :sigma2), Tuple{Real, Real}}
    function NormalNormalTransformer(
        normal_sampler::Sampler,
        new_params::Tuple{Real, Real},
        old_params::Tuple{Real, Real}=(0.0, 1.0)
    )
        # Name parameter tuples
        new_params = (mu = new_params[1], sigma2 = new_params[2])
        old_params = (mu = old_params[1], sigma2 = old_params[2])
        # Validation
        if !(old_params[:sigma2] > 0 && new_params[:sigma2] > 0)
            error("variance must be positive")
        end
        new(normal_sampler, new_params, old_params)
    end
end

@doc raw"""
    sample(s::NormalNormalTransformer, N)

Transform normal samples into normal samples with new parameters.

# Arguments
* `s::NormalNormalTransformer`: a normal-normal transformer.
* `N::Integer`: the number of samples to generate.

# Returns
* `samples::Array{Float64}`: a 2-D array in which each column is a sample.
"""
function sample(s::NormalNormalTransformer, N::Integer)
    old_samples = sample(s.normal_sampler, N)
    std_samples = (old_samples .- s.old_params[:mu]) ./
        sqrt(s.old_params[:sigma2])
    new_samples = sqrt(s.new_params[:sigma2]) .* std_samples .+
        s.new_params[:mu]

    return new_samples
end


text:
# Ket Constructors
# qubit constructors
export bloch_qubit_ket
export mirror_symmetric_qubit_kets, planar_symmetric_qubit_kets
export trine_qubit_kets, bb84_qubit_kets, sic_qubit_kets

# 2-qubit constructors
export bell_kets

# general constructors
export computational_basis_kets
export generalized_bell_kets

"""
    computational_basis_kets( dim :: Int64 ) :: Vector{Ket{Int64}}

The computational basis vectors for the Hilbert space of dimension, `dim`.

```math
\\Psi = \\{ |j\\rangle \\}_{j=0}^{(d-1)}
```
"""
function computational_basis_kets(dim::Int64) :: Vector{Ket{Int64}}
    Ket.(computational_basis_vectors(dim))
end

"""
    bloch_qubit_ket( θ :: Real, ϕ :: Real ) :: Ket{Complex{Float64}}

Returns the qubit ket for the specified spherical coordinate on the surface of
bloch sphere, `(r=1, θ, ϕ)`:

```math
|\\psi\\rangle = \\cos(\\theta/2)|0\\rangle + e^{i\\phi}\\sin(\\theta/2)|1\\rangle
```

If the ket does not have a phase then a real-valued `Ket` is constructed as:

```julia
bloch_qubit_ket( θ :: Real ) :: Ket{Float64}
```

Inputs:

* `θ ∈ [0, 2π]`: polar angle (w.r.t z-axis)
* `ϕ ∈ [0, 2π]`: azimuthal angle (x-y plane)

A `DomainError` is thrown if inputs `θ` and/or `ϕ` do are not within the valid range.
"""
function bloch_qubit_ket(θ::Real, ϕ::Real) :: Ket{Complex{Float64}}
    if !(0 <= θ <= 2π)
        throw(DomainError(θ, "polar bloch-angle (θ) must be in domain [0,2π]"))
    elseif !(0 <= ϕ <= 2π)
        throw(DomainError(ϕ, "azimuthal angle (ϕ) must be in domain [0,2π]"))
    end

    Ket([cos(θ/2),exp(im*ϕ)*sin(θ/2)])
end

function bloch_qubit_ket(θ::Real) :: Ket{Float64}
    if !(-2π <= θ <= 2π)
        throw(DomainError(θ, "polar bloch-angle (θ) must be in domain [0,2π]"))
    end

    Ket([cos(θ/2),sin(θ/2)])
end

"""
    bell_kets() :: Vector{Ket{Float64}}

The Bell basis kets for maximally entangled bipartite qubit systems.
These kets are ordered as ``\\{|\\Phi^+\\rangle, |\\Phi^-\\rangle, |\\Psi^+\\rangle, |\\Psi^-\\rangle \\}``, where

```math
\\begin{matrix}
    |\\Phi^+\\rangle = \\frac{1}{\\sqrt{2}}(|00\\rangle + |11\\rangle), &
|\\Phi^-\\rangle = \\frac{1}{\\sqrt{2}}(|00\\rangle - |11\\rangle), \\\\
    |\\Psi^+\\rangle = \\frac{1}{\\sqrt{2}}(|01\\rangle + |10\\rangle), &
    |\\Psi^-\\rangle = \\frac{1}{\\sqrt{2}}(|01\\rangle - |10\\rangle). \\\\
\\end{matrix}
```
"""
bell_kets() :: Vector{Ket{Float64}} = Ket.([
    1/sqrt(2)*[1,0,0,1],
    1/sqrt(2)*[1,0,0,-1],
    1/sqrt(2)*[0,1,1,0],
    1/sqrt(2)*[0,1,-1,0]
])

"""
    generalized_bell_kets( dim :: Int64 ) :: Vector{Ket{Float64}}

The Bell basis for entangled bipartite quantum states each of dimension `dim`. Each state is constructed by

```math
|\\Psi^p_c\\rangle = \\frac{1}{\\sqrt{d}}\\sum_{j=0}^{d-1} e^{i2\\pi pj/d}|j\\rangle |\\mod(j+c,d)\\rangle
```

where ``p,c\\in \\{0,\\cdots, (d-1)\\}`` and ``d`` is `dim`. When iterated, ``c`` is the major
index and ``p`` is the minor index.

A `DomainError` is thrown if `dim ≥ 2` is not satisfied.
"""
function generalized_bell_kets(dim :: Int64) :: Vector{Ket{Complex{Float64}}}
    if !(dim ≥ 2)
        throw(DomainError(dim, "Hilbert space dimension must satisfy `dim ≥ 2`"))
    end

    basis = computational_basis_vectors(dim)

    kets = map(
        c -> map(
            p -> Ket(sum(
                j -> exp(im*2*π*p*j/dim)*kron(basis[j+1],basis[mod(j + c, dim) + 1]),
                0:dim-1
            )/sqrt(dim)),
        0:dim-1),
    0:dim-1)

    collect(flatten(kets))
end

"""
    mirror_symmetric_qubit_kets( θ :: Real ) :: Vector{Ket{Float64}}

Returns the triplet of qubit kets in the x-z plane of bloch sphere. The first ket
is  ``|0\\rangle`` and the other two are symmetric across the z-axis,
``|\\pm\\rangle = \\cos(\\theta)|0 \\rangle \\pm \\sin(\\theta)|1\\rangle``.

Input:

* `θ ∈ [0,π/2]`: the hilbert space angle between ``|0\\rangle`` and symmetric kets.
"""
function mirror_symmetric_qubit_kets(θ::Real) :: Vector{Ket{Float64}}
    if !(0 <= θ <= π/2)
        throw(DomainError(θ, "Hilbert space angle θ must satisfy `0 <= θ <= π/2`"))
    end

    ψ1 = Ket([1.,0.])
    ψ2 = bloch_qubit_ket(2*θ)
    ψ3 = bloch_qubit_ket(-2*θ)

    [ψ1,ψ2,ψ3]
end

"""
    planar_symmetric_qubit_kets( n :: Int64 ) :: Vector{Ket{Float64}}

Constructs a set of `n ``Ket`s oriented symmetrically in the x-z-plane.
Each ket is separated by a bloch angle of `2π/n`. The `Ket`s are constructed
with the form:

```math
|\\psi_j \\rangle = \\cos(j \\pi/n) | 0\\rangle + \\sin(j \\pi/n)|1\\rangle
```

where ``j \\in \\{0,\\cdots, (n-1)\\}``.

A `DomainError` is thrown if `n < 2`.
"""
function planar_symmetric_qubit_kets(n :: Int64) :: Vector{Ket{Float64}}
    if !(n ≥ 2)
        throw(DomainError(n, "n must satisfy `n ≥ 2`"))
    end

    map(θ -> (θ ≤ π) ? bloch_qubit_ket(θ) : bloch_qubit_ket(θ - 2π), 0:2π/n:2π*(n-1)/n+π/(2*n))
end

"""
    trine_qubit_kets() :: Vector{Ket{Float64}}

The triplet of `Ket`s separated by equal angles in the x-z plane of bloch sphere.

```math
    |\\psi_1\\rangle = |0\\rangle, \\quad
    |\\psi_2\\rangle = \\frac{1}{2}|0\\rangle + \\frac{\\sqrt{3}}{2}|1\\rangle, \\quad
    |\\psi_3\\rangle = \\frac{1}{2}|0\\rangle - \\frac{\\sqrt{3}}{2}|1\\rangle
```
"""
trine_qubit_kets() :: Vector{Ket{Float64}} = Ket.([
    [1.0, 0], [0.5, sqrt(3)/2], [0.5, -sqrt(3)/2]
])

"""
    sic_qubit_kets() :: Vector{Ket{Complex{Float64}}}

The quadruplet of symmetric informationally complete (SIC) qubits. This set of qubits
correspond to the vertices of a tetrahedron inscribed on bloch sphere.

```math
\\begin{matrix}
    |\\psi_1\\rangle = |0\\rangle, & \\quad |\\psi_2\\rangle = \\frac{1}{\\sqrt{3}}|0\\rangle + \\sqrt{\\frac{2}{3}}|1\\rangle, \\\\
    |\\psi_3\\rangle = \\frac{1}{\\sqrt{3}}|0\\rangle + \\sqrt{\\frac{2}{3}} e^{i 2\\pi/3}|1\\rangle, & \\quad |\\psi_4\\rangle = \\frac{1}{\\sqrt{3}}|0\\rangle + \\sqrt{\\frac{2}{3}} e^{i 4\\pi/3}|1\\rangle
\\end{matrix}
```
"""
sic_qubit_kets() :: Vector{Ket{Complex{Float64}}} = Ket.([
    [1., 0im], [1/sqrt(3), sqrt(2/3)+ 0im],
    [1/sqrt(3), sqrt(2/3)*exp(im*2π/3)], [1/sqrt(3), sqrt(2/3)*exp(im*4π/3)]
])

"""
    bb84_qubit_kets :: Vector{Ket{Float64}}

The quadruplet of qubit kets used in the BB84 Quantum Key Distribution protocol. The
states are ``|0\\rangle``, ``|+\\rangle``, ``|1\\rangle``, and ``|- \\rangle``.
"""
bb84_qubit_kets() :: Vector{Ket{Float64}} = Ket.([
    [1,0], [1,1]/sqrt(2), [0,1], [1,-1]/sqrt(2)
])

ext: jl


text:
"""
$(TYPEDEF)

Partial prior belief (absolute data) on any variable, given `<:SamplableBelief` and which dimensions of the intended variable.

Notes
- If using [`AMP.ManifoldKernelDensity`](@ref), don't double partial.  Only define the partial in this `PartialPrior` container.
  - Future TBD, consider using `AMP.getManifoldPartial` for more general abstraction.
"""
struct PartialPrior{T <: SamplableBelief,P <: Tuple} <: AbstractPrior
  Z::T
  partial::P
end

getManifold(pp::PartialPrior{<:PackedManifoldKernelDensity}) = pp.Z.manifold

getSample(cf::CalcFactor{<:PartialPrior}) = samplePoint(cf.factor.Z)

"""
$(TYPEDEF)

Serialization type for `PartialPrior`.
"""
mutable struct PackedPartialPrior <: PackedInferenceType
  Z::String
  partials::Vector{Int}
end

function convert(::Type{PackedPartialPrior}, d::PartialPrior)
  PackedPartialPrior(convert(PackedSamplableBelief, d.Z), [d.partial...;])
end
function convert(::Type{PartialPrior}, d::PackedPartialPrior)
  PartialPrior(convert(SamplableBelief, d.Z),(d.partials...,))
end


text:
using FactNSF
using Test

@testset "FactNSF.jl" begin
    # Write your tests here.
end

max_stars_repo_name: xinshuoyang/FactNSF

max_issues_repo_name: xinshuoyang/FactNSF

max_forks_repo_name: xinshuoyang/FactNSF


text:
"""
    LiuScore{TI<:Integer,TR<:Real}

    LiuScore(
      samplespace::SampleSpace,
      p0::Real, pmcrv::Real, prior,
      alpha::Real, beta::Real,
      fs::Real, fp::Real;
      npriorpivots::Integer = 50,
      npivots::Integer = 15,
      minconditionalpower::Real = 0.0,
      MONOTONECONDITIONALPOWER::Bool = false,
      label::String = ""
    )

This type represents a set of parameters for finding optimal two-stage designs
minimizing a variant of the criterion proposed by Liu et al. in [1].

> [1] Liu G, Frank G, Zhu C, and Lu C. Evaluating the adaptive performance of flexible sample size designs with treatment difference in an interval. `Statistics in Medicine` 2008; 27(4): 584-596.

# Parameters

| Parameter    | Description |
| -----------: | :---------- |
| samplespace  | a sample space object |
| p0           | upper boundary of the null hypothesis |
| pmcrv        | point alternative to power on |
| prior        | prior distribution on true response rate, expected score under the prior conditional on a relevant effect is minimized |
| alpha        | maximal tolerable type one error rate |
| beta         | maximal tolerable type two error rate on p1 |
| fs           | scaling parameter as described in [1] |
| fp           | scaling parameter as described in [1] |
| npriorpivots | number of pivot points for integration of prior |
| npivots      | number of pivot points piecewise linear approximation of local score |
| minconditionalpower | minimal conditional power upon continuation to stage two |
| MONOTONECONDITIONALPOWER | if true, the conditional power must be monotonously increasing, this constraint is only relevant if nmax is set very restrictively |
| label        | descriptive string for parameter set |
"""
mutable struct LiuScore{TI<:Integer,TR<:Real} <: VagueAlternative

  label::String
  samplespace::SampleSpace
  p0::TR
  pmcrv::TR
  prior
  mtoer::TR
  beta::TR
  mincondpower::TR
  fs::TR
  fp::TR
  MONOTONECONDITIONALPOWER::Bool
  npriorpivots::TI # number of pivots for prior evaluation
  npivots::TI # number of pivots for approximation of g

  function LiuScore{TI,TR}(
    label::String,
    samplespace::SampleSpace,
    p0::TR, pmcrv::TR, prior,
    alpha::TR, beta::TR,
    mincondpower::TR,
    fs::TR, fp::TR,
    MONOTONECONDITIONALPOWER::Bool,
    npriorpivots::TI, # number of pivots for prior evaluation
    npivots::TI, # number of pivots for approximation of g
  ) where {TI<:Integer,TR<:Real}

    @checkprob alpha beta p0 pmcrv mincondpower
    z = quadgk(prior, 0, 1, abstol = .0001)[1]
    abs(z - 1) > .001 ?
      error(@sprintf("prior must integrate to one, is %.6f", z)) : nothing
    # compute conditional prior as weight function
    zz = quadgk(prior, pmcrv, 1, abstol = 0.0001)[1]
    conditionalprior = p -> (p > pmcrv ? prior(p) / zz : 0)
    fp < 0 ?
      error(@sprintf("fp must be positive, is %.4f", fp)) : nothing
    fp >= 1 ?
      error(@sprintf("fp must be smaller than 1, is %.4f", fp)) : nothing
    fs <= 1 ?
      error(@sprintf("df must be greater than 1, is %.4f", fs)) : nothing
    new(
      label,
      samplespace,
      p0, pmcrv, conditionalprior,
      alpha, beta,
      mincondpower,
      fs, fp,
      MONOTONECONDITIONALPOWER,
      npriorpivots, npivots
    )

  end # inner constructor

end # LiuScore

function LiuScore(
  samplespace::SampleSpace,
  p0::Real, pmcrv::Real, prior,
  alpha::Real, beta::Real,
  fs::Real, fp::Real;
  npriorpivots::Integer = 50,
  npivots::Integer = 15,
  minconditionalpower::Real = 0.0,
  MONOTONECONDITIONALPOWER::Bool = false,
  label::String = ""
)

  # unify types
  p0, pmcrv, fs, fp, alpha, beta, minconditionalpower =
    promote(p0, pmcrv, fs, fp, alpha, beta, minconditionalpower)
  npriorpivots, npivots = promote(npriorpivots, npivots)
  LiuScore{typeof(npivots),typeof(p0)}(
    label, samplespace, p0, pmcrv, prior, alpha, beta, minconditionalpower, fs, fp,
    MONOTONECONDITIONALPOWER, npriorpivots, npivots
  )

end # LiuScore

maxsamplesize(params::LiuScore) = maxsamplesize(params.samplespace)

hasmonotoneconditionalpower(params::LiuScore) = params.MONOTONECONDITIONALPOWER

minconditionalpower(params::LiuScore) = params.mincondpower

targetpower(params::LiuScore) = 1 - params.beta

Base.show(io::IO, params::LiuScore) = print("LiuScore")

function ros(design::Design, params::LiuScore, p1::T) where {T<:Real}

    @checkprob p1
    powerreq = targetpower(params)
    n1 = interimsamplesize(design)
    nreq__(p, powerreq) = nreq_(p, powerreq, params.p0, mtoer(params))
    denom = nreq__(p1, powerreq)
    if denom == 0
        return 100
    end
    ros = 0.0
    for x1 in 0:n1
        ros += dbinom(x1, n1, p1) * max(0.0, samplesize(design, x1)/denom - 1)
    end
    return min(100, ros / (params.fs - 1))

end # ros

function rup(design::Design, params::LiuScore, p1::T) where {T<:Real}

    @checkprob p1
    if p1 < params.pmcrv
        return 0.0
    end
    powerreq = targetpower(params) # we can minimize this if we set power to zero!
    n1 = interimsamplesize(design)
    nreq__(p, powerreq) = nreq_(p, powerreq, params.p0, mtoer(params))
    nom   = max(0, nreq__(p1, powerreq) - nreq__(p1, power(design, p1)))
    denom = nreq__(p1, powerreq) - nreq__(p1, (1 - params.fp)*powerreq)
    if (denom == 0.0) & (nom > 0)
        return 100
    end
    if (denom == 0.0) & (nom == 0)
        return 0.0
    end
    res = nom / denom
    return min(100, res)

end # rup

function score(
  design::Design, params::LiuScore, p1::T
) where {T<:Real}

  return ros(design, params, p1) + rup(design, params, p1)

end # score

function score(design::Design, params::LiuScore)

    function f(p)
        res = params.prior(p)*score(design, params, p)
        return min(100, res)
    end

    return quadgk(f, params.pmcrv, 1, reltol = .001, maxevals = 1e5)[1]

end # score

function completemodel(ipm::IPModel, params::LiuScore, n1::T) where {T<:Integer}

    ss = samplespace(params)
    !possible(n1, ss) ? error("n1 and sample space incompatible") : nothing
    # extract ip model
    m             = ipm.m
    y             = ipm.y
    nvals         = ipm.nvals
    cvals         = ipm.cvals
    cvalsfinite   = ipm.cvalsfinite
    cvalsinfinite = [-Inf; Inf]
    # extract other parameters
    nmax          = maxsamplesize(ss, n1)
    p0            = null(params)
    prior         = params.prior
    npriorpivots  = params.npriorpivots
    npivots      = params.npivots
    # get grid for prior and conditional prior
    priorpivots, dcdf   = findgrid(prior, 0, 1, npriorpivots)
    cpriorpivots, dccdf = findgrid(prior, params.pmcrv, 1, 1000) # not performance critical!
    # add type one error rate constraint
    JuMP.@constraint(m,
        sum(dbinom(x1, n1, p0)*_cpr(x1, n1, n, c, p0)*y[x1, n, c] for
            x1 in 0:n1, n in nvals, c in cvals
        ) <= mtoer(params)
    )
    # add conditional type two error rate constraint (power)
    for x1 in 0:n1
        posterior1 = dbinom.(x1, n1, cpriorpivots) .* dccdf
        z1 = sum(posterior1)
        JuMP.@constraint(m,
            sum(sum(posterior1 .* _cpr.(x1, n1, n, c, cpriorpivots))/z1*y[x1, n, c] for
                n in nvals, c in cvalsfinite
            ) + sum(y[x1, n, c] for
                n in nvals, c in cvalsinfinite
            ) >= minconditionalpower(params)
        )
        # ensure monotonicity if required
        if (x1 >= 1) & hasmonotoneconditionalpower(params)
            posterior2 = dbinom.(x1 - 1, n1, cpriorpivots) .* dccdf
            z2 = sum(posterior2)
            JuMP.@constraint(m,
                sum(sum(posterior1 .* _cpr.(x1, n1, n, c, cpriorpivots))/z1*y[x1, n, c]
                        - sum(posterior2 .* _cpr.(x1 - 1, n1, n, c, cpriorpivots))/z2*y[x1 - 1, n, c] for
                    n in nvals, c in cvals
                ) >= 0
            )
        end
    end
    # add optimality criterion
    nreq__(p, powerreq) = nreq_(p, powerreq, p0, mtoer(params))
    # construct expressions for ROS, upper bound necessary to guarantee finteness!
    JuMP.@expression(m, ros[p in priorpivots],
        1/(params.fs - 1) * sum(
            dbinom(x1, n1, p) * max(0, min(100.0, (n / nreq__(p, 1 - params.beta) - 1))) * y[x1, n, c] for
            x1 in 0:n1, n in nvals, c in cvals
        )
    )
    # construct expressions for power
    JuMP.@expression(m, designpower[p in priorpivots],
        sum(dbinom(x1, n1, p) * _cpr(x1, n1, n, c, p) * y[x1, n, c] for
            x1 in 0:n1, n in nvals, c in cvals
        )
    )
    pivots = [collect(linspace(0, 1 - params.beta, npivots - 1)); 1] # lambda formulation requires edges!, rup is 0 from powerreq to 1
    JuMP.@variable(m, 0 <= lambda[priorpivots, pivots] <= 1)
    function frup(power, p; rupmax = 100.0)
        if p < params.pmcrv
            return 0.0
        end
        nom   = max(0, nreq__(p, 1 - params.beta) - nreq__(p, power))
        denom = nreq__(p, 1 - params.beta) - nreq__(p, (1 - params.fp)*(1 - params.beta))
        if (denom == 0.0) & (nom > 0)
            return rupmax
        end
        if (denom == 0.0) & (nom == 0)
            return 0.0
        end
        res = max(0, min(rupmax, nom / denom))
        return res
    end
    JuMP.@variable(m, rup[priorpivots])
    for p in priorpivots
        JuMP.addSOS2(m, [lambda[p, piv] for piv in pivots])
        JuMP.@constraint(m, sum(lambda[p, piv] for piv in pivots) == 1)
        JuMP.@constraint(m, sum(lambda[p, piv]*piv for piv in pivots) == designpower[p]) # defines lambdas!
        JuMP.@constraint(m, sum(lambda[p, piv]*frup(piv, p) for piv in pivots) == rup[p])
    end
    JuMP.@objective(m, Min,
        sum((rup[priorpivots[i]] + ros[priorpivots[i]])*dcdf[i] for i in 1:npriorpivots) # normalized version
    )
    return true

end # completemodel

# ToDo: implement
_isfeasible(design::Design, params::LiuScore) = true

# utility
function nreq_(p, power, p0, alpha; nmax = 10000.0)

    if alpha >= power
        return 0.0 # qnorm(1 - alpha) + qnorm(power) < 0 ....
    end
    if p <= p0
        return 0.0 # no effect
    end
    if power == 0
        return 0.0
    end
    if abs(power - 1.0) < 5 * eps(Float64)
        return nmax
    end
    res = (1 - p) * p * ( (qnorm(1 - alpha) + qnorm(power)) / (p0 - p) )^2
    res = min(nmax, max(0, res))
    return res

end

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
plan_flows_fft(a::Array, effort) = plan_fft(a; flags=effort)
plan_flows_rfft(a::Array, effort) = plan_rfft(a; flags=effort)

"""
    ZeroDGrid()

Constructs a placeholder grid object for "0D" problems (in other words, systems of ODEs).
"""
struct ZeroDGrid{T, Ta} <: AbstractGrid{T, Ta} end

function getaliasedwavenumbers(nk, nkr, aliasfraction)
  # Index endpoints for aliased i, j wavenumbers
  # 1/3 aliasfraction => upper 1/6 of +/- wavenumbers (1/3 total) are set to 0 after performing fft.
  # 1/2 aliasfraction => upper 1/4 of +/- wavenumbers (1/2 total) are set to 0 after performing fft.
  L = (1 - aliasfraction)/2 # (1 - 1/3) / 2 + 1 = 1/3.
  R = (1 + aliasfraction)/2 # (1 + 1/3) / 2 - 1 = 2/3.
  iL = floor(Int, L*nk) + 1
  iR =  ceil(Int, R*nk)

  aliasfraction < 1 || error("`aliasfraction` must be less than 1") # aliasfraction=1 is not sensible.
   kalias = (aliasfraction > 0) ? (iL:iR) : Int(nx/2+1)
  kralias = (aliasfraction > 0) ? (iL:nkr) : nkr

  kalias, kralias
end

"""
    OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE)

Constructs a OneDGrid object with size `Lx`, resolution `nx`, and leftmost
position `x0`. FFT plans are generated for `nthreads` CPUs using
FFTW flag `effort`.
"""
struct OneDGrid{T<:AbstractFloat, Ta<:AbstractArray, Tfft, Trfft} <: AbstractGrid{T, Ta}
        nx :: Int
        nk :: Int
       nkr :: Int

        dx :: T
        Lx :: T

         x :: Ta
         k :: Ta
        kr :: Ta
    invksq :: Ta
   invkrsq :: Ta

   fftplan :: Tfft
  rfftplan :: Trfft

  # Range objects that access the aliased part of the wavenumber range
    kalias :: UnitRange{Int}
   kralias :: UnitRange{Int}
end

function OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE, T=Float64, dealias=1/3,
                  ArrayType=Array)

  dx = Lx/nx
  x = ArrayType{T}(range(x0, step=dx, length=nx))

  nk = nx
  nkr = Int(nx/2+1)

  i₁ = 0:Int(nx/2)
  i₂ = Int(-nx/2+1):-1
   k = ArrayType{T}(2π/Lx*cat(i₁, i₂; dims=1))
  kr = ArrayType{T}(2π/Lx*cat(i₁; dims=1))

   invksq = @. 1/k^2
  invkrsq = @. 1/kr^2
   invksq[1] = 0
  invkrsq[1] = 0

  FFTW.set_num_threads(nthreads)
   fftplan = plan_flows_fft(ArrayType{Complex{T}, 1}(undef, nx), effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 1}(undef, nx), effort)

  kalias, kralias = getaliasedwavenumbers(nk, nkr, dealias)

  Ta = typeof(x)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)

  OneDGrid{T, Ta, Tfft, Trfft}(nx, nk, nkr, dx, Lx, x, k, kr,
                               invksq, invkrsq, fftplan, rfftplan, kalias, kralias)
end

"""
    TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE)

Constructs a TwoDGrid object.
"""
struct TwoDGrid{T<:AbstractFloat, Ta<:AbstractArray, Tfft, Trfft} <: AbstractGrid{T, Ta}
  nx::Int
  ny::Int
  nk::Int
  nl::Int
  nkr::Int

  dx::T
  dy::T
  Lx::T
  Ly::T

  x::Ta
  y::Ta
  k::Ta
  l::Ta
  kr::Ta
  Ksq::Ta
  invKsq::Ta
  Krsq::Ta
  invKrsq::Ta

  fftplan::Tfft
  rfftplan::Trfft

  # Range objects that access the aliased part of the wavenumber range
  kalias::UnitRange{Int}
  kralias::UnitRange{Int}
  lalias::UnitRange{Int}
end

function TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS,
                  effort=FFTW.MEASURE, T=Float64, dealias=1/3, ArrayType=Array)

  dx = Lx/nx
  dy = Ly/ny

  nk = nx
  nl = ny
  nkr = Int(nx/2+1)

  # Physical grid
  x = ArrayType{T}(reshape(range(x0, step=dx, length=nx), (nx, 1)))
  y = ArrayType{T}(reshape(range(y0, step=dy, length=ny), (1, ny)))

  # Wavenubmer grid
  i₁ = 0:Int(nx/2)
  i₂ = Int(-nx/2+1):-1
  j₁ = 0:Int(ny/2)
  j₂ = Int(-ny/2+1):-1

   k = ArrayType{T}(reshape(2π/Lx*cat(i₁, i₂, dims=1), (nk, 1)))
   l = ArrayType{T}(reshape(2π/Ly*cat(j₁, j₂, dims=1), (1, nl)))
  kr = ArrayType{T}(reshape(2π/Lx*cat(i₁, dims=1), (nkr, 1)))

     Ksq = @. k^2 + l^2
  invKsq = @. 1/Ksq
  invKsq[1, 1] = 0

     Krsq = @. kr^2 + l^2
  invKrsq = @. 1/Krsq
  invKrsq[1, 1] = 0

  # FFT plans
  FFTW.set_num_threads(nthreads)
  fftplan = plan_flows_fft(ArrayType{Complex{T}, 2}(undef, nx, ny), effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 2}(undef, nx, ny), effort)

  # Index endpoints for aliasfrac i, j wavenumbers
  kalias, kralias = getaliasedwavenumbers(nk, nkr, dealias)
  lalias, _ = getaliasedwavenumbers(nl, nl, dealias)

  Ta = typeof(x)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)

  TwoDGrid{T, Ta, Tfft, Trfft}(nx, ny, nk, nl, nkr, dx, dy, Lx, Ly, x, y, k, l, kr, Ksq, invKsq, Krsq, invKrsq,
           fftplan, rfftplan, kalias, kralias, lalias)
end

OneDGrid(dev::CPU, args...; kwargs...) = OneDGrid(args...; ArrayType=Array, kwargs...)
TwoDGrid(dev::CPU, args...; kwargs...) = TwoDGrid(args...; ArrayType=Array, kwargs...)

"""
    gridpoints(g)

Returns the collocation points of the grid `g` in 2D arrays `X, Y`.
"""
function gridpoints(g::AbstractGrid{T, A}) where {T, A}
  X = [ g.x[i] for i=1:g.nx, j=1:g.ny]
  Y = [ g.y[j] for i=1:g.nx, j=1:g.ny]
  A(X), A(Y)
end

"""
    dealias!(a, g, kalias)

Dealias `a` on the grid `g` with aliased x-wavenumbers `kalias`.
"""
function dealias!(a, g::OneDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)
  nothing
end

function dealias!(a, g::OneDGrid, kalias)
  @views @. a[kalias, :] = 0
  nothing
end

function dealias!(a, g::TwoDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)
  nothing
end

function dealias!(a, g::TwoDGrid, kalias)
  @views @. a[kalias, g.lalias, :] = 0
  nothing
end

"""
    makefilter(K; order=4, innerK=0.65, outerK=1)

Returns a filter acting on the non-dimensional wavenumber K that decays exponentially
for K>innerK, thus removing high-wavenumber content from a spectrum it is multiplied with.
The decay rate is determined by order and outerK determines the outer wavenumber at which
the filter is smaller than Float64 machine precision.
"""
function makefilter(K::Array; order=4, innerK=0.65, outerK=1)
  TK = typeof(K)
  K = Array(K)
  decay = 15*log(10) / (outerK-innerK)^order # decay rate for filtering function
  filt = @. exp( -decay*(K-innerK)^order )
  filt[real.(K) .< innerK] .= 1
  TK(filt)
end

function makefilter(g::TwoDGrid; realvars=true, kwargs...)
  K = realvars ?
      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2))
  makefilter(K; kwargs...)
end

function makefilter(g::OneDGrid; realvars=true, kwargs...)
  K = realvars ? g.kr*g.dx/π : @.(abs(g.k*g.dx/π))
  makefilter(K; kwargs...)
end

makefilter(g, T, sz; kwargs...) = ones(T, sz) .* makefilter(g; realvars=sz[1]==g.nkr, kwargs...)
makefilter(eq) = makefilter(eq.grid, fltype(eq.T), eq.dims)


text:
@testset "ACEtb.SlaterKoster Kwon Test" begin

using ACEtb.SlaterKoster, JuLIP, Test, LinearAlgebra
OldSK = ACEtb.SlaterKoster.OldSK

using SlaterKoster: assembleRI

@info("Old Kwon Tests... just checking that it still runs in principle")
oldkwon = OldSK.KwonHamiltonian()
at = bulk(:Si, cubic = true, pbc = false) * 2
H, M = OldSK.hamiltonian(oldkwon, at)
println(@test M == Matrix((1.0+0.0im)*I, size(M)))

@info("Test Old Kwon versus new Kwon hamiltonians")
ats = [ rattle!(bulk(:Si, cubic = true, pbc = false) * 2, 0.1),
        rattle!(bulk(:Si, cubic = true, pbc = true) * 2, 0.1),
        rattle!(bulk(:Si, cubic = false, pbc = true) * 3, 0.1) ]
kwon = KwonTB()

for at in ats
   Hold, Mold = OldSK.hamiltonian(oldkwon, at)
   # for i = 1:size(Hold,1); Hold[i,i] = 0; end
   eigold = eigvals(Hold)

   Hri = assembleRI(at, kwon)
   Hc = Matrix(Hri)
   eignew = eigvals(Symmetric(Hc))

   println(@test sort(eigold) ≈ sort(eignew))
end

end

# # -----------------------------------------------------------------------------
# # visual check that the spectrum looks ok - need to make this a proper
# # test set
# # -----------------------------------------------------------------------------
# # E = eigvals(H)
# # using Plots
# # plot(E, lw=0, ms=5, label = "")
# # eF = 0.5 * (E[length(E)÷2] + E[length(E)÷2+1])
# # hline!([eF], style=:dot, c=:black, label ="")
# # -----------------------------------------------------------------------------


text:
using JCheck

using Documenter

DocMeta.setdocmeta!(JCheck,
                    :DocTestSetup,
                    :(using JCheck),
                    recursive = true,
                    warn = false)

makedocs(sitename="JCheck.jl",
         doctest = true,
         modules = [JCheck],
         pages = ["Home" => "index.md",
                  "Reference" => "reference.md"])

deploydocs(repo = "github.com/ps-pat/JCheck.jl.git",
           devbranch = "main")

Documenter.HTML(assets = [asset("assets/logo.svg",
                                class = :ico,
                                islocal = true)])


text:
@safetestset "isempty" begin
    using Validator

    # ignore_whitespace=false
    @test Validator.isempty("") === true
    @test Validator.isempty(" ") === false
    @test Validator.isempty("foo") === false

    # ignore_whitespace=true
    @test Validator.isempty("", true) === true
    @test Validator.isempty(" ", true) === true
    @test Validator.isempty("foo", true) === false
end


text:
export RadiationModel, NoRadiation

abstract type RadiationModel end

vars_state(::RadiationModel, FT) = @vars()
vars_aux(::RadiationModel, FT) = @vars()
vars_integrals(::RadiationModel, FT) = @vars()
vars_reverse_integrals(::RadiationModel, FT) = @vars()

function atmos_nodal_update_aux!(
    ::RadiationModel,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function preodefun!(::RadiationModel, aux::Vars, state::Vars, t::Real) end
function integral_set_aux!(
    ::RadiationModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) end
function integral_load_aux!(::RadiationModel, aux::Vars, integ::Vars) end
function reverse_integral_set_aux!(
    ::RadiationModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
) end
function reverse_integral_load_aux!(::RadiationModel, aux::Vars, integ::Vars) end
function flux_radiation!(
    ::RadiationModel,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end

struct NoRadiation <: RadiationModel end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using DataFrames
using MacroTools
using Recurrences

function loop(xs::Vector{Expr}, f::Function; init::Dict, iterations::Int = 10)
    ls = [:(Symbol($(string(v))) => $v) for v in keys(init)]
    is = [:($k = $v) for (k, v) in init]
    quote
        let $(is...)
            $f(:it=>0, $(ls...))
            for it in 1:$iterations
                $(xs...)
                $f(:it=>it, $(ls...))
            end
        end
    end
end

function loop(cs::Vector{<:ClosedForm}, f::Function; init::Dict, iterations::Int = 10)
    xs = [convert(Expr, c) for c in cs]
    xs = [MacroTools.postwalk(x->(x isa Number ? Rational(x) : x), y) for y in xs]
    vs = [Symbol(string(c.func)) for c in cs]
    is = [k in vs ? :($(initvar(k, 0)) = $v) : :($k = $v) for (k, v) in init]
    ls = [v in vs ? :(Symbol($(string(v))) => $v(it)) : :(Symbol($(string(v))) => $v) for v in keys(init)]
    quote
        let $(is...)
            $(xs...)
            for it in 0:$iterations
                $f(:it=>it, $(ls...))
            end
        end
    end
end

function collecttrace(xs; init, kwargs...)
    df = DataFrame()
    df[:it] = Int[]
    for v in keys(init)
        df[v] = Rational[]
    end
    eval(loop(xs, (x...)->push!(df, Dict(x)); init=init, kwargs...))
    df
end

function cftrace(xs::Vector{Expr}; kwargs...)
    cs = lrs_sequential(xs) |> solve
    collecttrace(cs; kwargs...)
end

function trace(xs::Vector{Expr}; kwargs...)
    vs = Base.unique(Iterators.flatten(map(Recurrences.free_symbols, xs)))
    dict = Dict(v=>Rational(rand(-10:10)) for v in vs)
    ldf = collecttrace(xs; kwargs..., init=dict)
    cdf = cftrace(xs; kwargs..., init=dict)
    ldf, cdf
end


text:
using Test, Random
using StatsBase, Distributions, DependentBootstrap
using ForecastEval

#-------------------------------------------------
#DO NOT CHANGE ANY OF THE FOLLOWING INPUTS
kernelfunctypevec = [ForecastEval.KernelEpanechnikov(), ForecastEval.KernelGaussian(), ForecastEval.KernelBartlett(), ForecastEval.KernelUniform()];
kernelfuncsymvec = [:epanechnikov, :gaussian, :bartlett, :uniform];
bandwidthvec = [-1, 0, 1, 2, 3];
bootmethodsymvec = [:iid, :stationary];
bootmethodtypevec = [DependentBootstrap.BootIID(), DependentBootstrap.BootStationary()]
blocklengthvec = [0.0, 5.0];
xmat =
[-0.371817    0.672154  -1.33346    1.3462       0.321323 ;
 0.0676651  -1.03904    1.12363    0.872374    -0.508987 ;
-1.00231    -0.842998  -1.76564    0.00652721  -1.73871  ;
-0.165286   -0.503325  -1.09233   -0.22975     -0.906645 ;
 0.248638   -0.602089   1.02422    1.24328      1.54605  ;
 0.465047    1.11979    1.09211    0.693409     0.116323 ;
 1.53365     0.968045   0.135202   0.252632    -3.08523  ;
 0.600273   -2.00259   -0.488497   1.17        -1.27601  ;
 1.36719    -0.173315   1.1944    -1.57142      2.36206  ;
 1.22798     0.264002   0.427104  -0.0986267    0.508668 ;
 0.229319    1.6132    -2.14253   -1.58257     -0.530612 ;
 0.0410861  -0.258605  -0.321659  -1.20498     -0.63559  ;
 1.915       1.3493    -0.701988   0.125084     0.683679 ;
 1.11015     1.37006    1.63997   -0.0548267   -0.634976 ;
 0.106914    0.108819   0.535538  -0.073909     1.75006  ;
-0.247781   -1.42433   -0.890601  -0.694821     2.06756  ;
 0.312066    0.885002  -1.39183    0.763129     1.52937  ;
 0.488613   -0.172177  -0.224076  -0.13952     -1.47516  ;
 0.153926    0.422725   0.804834   0.141964     0.761245 ;
-0.679447   -0.502743   1.00025   -0.734587     0.119061];
xmatabs = abs.(xmat);
#-------------------------------------------------

@testset "Diebold-Mariano tests" begin
	#Generate data
	x = xmat[:,1]
	y = xmat[:,2]
	z = xmat[:,3]
	ld = (x - z).^2 - (y - z).^2
	#Test constructor
	for kk = 1:length(kernelfunctypevec)
		for bw in bandwidthvec
			dmhac1 = ForecastEval.DMHAC(alpha=0.05, kernelfunction=kernelfunctypevec[kk], bandwidth=bw)
			@test dmhac1.alpha == 0.05
			@test dmhac1.kernelfunction == kernelfunctypevec[kk]
			@test dmhac1.bandwidth == bw
			dmhac2 = ForecastEval.DMHAC(alpha=0.05, kernelfunction=kernelfuncsymvec[kk], bandwidth=bw)
			@test dmhac1.alpha == 0.05
			@test dmhac1.kernelfunction == kernelfunctypevec[kk]
			@test dmhac1.bandwidth == bw
			dmhac3 = ForecastEval.DMHAC(ld, alpha=0.05, kernelfunction=kernelfunctypevec[kk], bandwidth=bw)
			@test dmhac3.alpha == 0.05
			@test dmhac3.kernelfunction == kernelfunctypevec[kk]
			@test dmhac3.bandwidth == bw
			dmhac4 = ForecastEval.DMHAC(ld, alpha=0.05, kernelfunction=kernelfuncsymvec[kk], bandwidth=bw)
			@test dmhac4.alpha == 0.05
			@test dmhac4.kernelfunction == kernelfunctypevec[kk]
			@test dmhac4.bandwidth == bw
		end
	end
	for kk = 1:length(bootmethodtypevec)
		for bl in blocklengthvec
			dmboot1 = ForecastEval.DMBoot(ld, alpha=0.05, bootmethod=bootmethodtypevec[kk], blocklength=bl)
			@test dmboot1.alpha == 0.05
			@test dmboot1.bootinput.bootmethod == bootmethodtypevec[kk]
			if bl <= 0.0 ; @test dmboot1.bootinput.blocklength > 0.0
			else
				if bootmethodsymvec[kk] == :iid ; @test dmboot1.bootinput.blocklength == 1.0
				else @test dmboot1.bootinput.blocklength == bl
				end
			end
			dmboot2 = ForecastEval.DMBoot(ld, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=bl)
			@test dmboot2.alpha == 0.05
			@test dmboot2.bootinput.bootmethod == bootmethodtypevec[kk]
			if bl <= 0.0 ; @test dmboot2.bootinput.blocklength > 0.0
			else
				if bootmethodsymvec[kk] == :iid ; @test dmboot1.bootinput.blocklength == 1.0
				else @test dmboot1.bootinput.blocklength == bl
				end
			end
		end
	end
	#Test output (hac)
	pvalvec = [0.11440139911128933,0.05985060000192489,0.1211016948931237,0.026592349584142416,0.14442472227131237,0.14442472227131237,0.14442472227131237,0.14442472227131237,0.14442472227131237,0.12870301609991372,0.14442472227131237,0.11825398287003397,0.12491428583032126,0.07056983889420516,0.13149963805876264,0.03576447246936558,0.07469037167648004,0.06832784498031234,0.09957810645179305,0.06041771012189759]
	pvalmat = reshape(pvalvec, (4,5))
	teststatvec = [-1.5787151898146972,-1.8818925954728518,-1.5501659301885555,-2.217450096879554,-1.4595102474128858,-1.4595102474128858,-1.4595102474128858,-1.4595102474128858,-1.4595102474128858,-1.5192363442451382,-1.4595102474128858,-1.562144367590813,-1.534469114610218,-1.808235607809477,-1.5082144973278373,-2.0995951195178453,-1.782361057698249,-1.822838564430011,-1.6469024100931644,-1.8777329254314343]
	teststatmat = reshape(teststatvec, (4,5))
	for kk = 1:length(kernelfuncsymvec)
		for nn = 1:length(bandwidthvec)
			dmout = dm(ld, :hac, alpha=0.05, kernelfunction=kernelfuncsymvec[kk], bandwidth=bandwidthvec[nn])
			@test isapprox(dmout.pvalue, pvalmat[kk,nn])
			@test isapprox(dmout.teststat, teststatmat[kk,nn])
		end
	end
	#Test output (boot)
	Random.seed!(1234)
	pvalmat = [0.094  0.088 ; 0.106  0.026]
	for kk = 1:length(bootmethodsymvec)
		for nn = 1:length(blocklengthvec)
			dmout = dm(ld, :boot, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=blocklengthvec[nn])
			@test isapprox(pvalmat[kk,nn], dmout.pvalue)
		end
	end

end

@testset "Reality Check tests" begin
	#Test constructor
	for kk = 1:length(bootmethodsymvec)
		for nn = 1:length(blocklengthvec)
			rc1 = ForecastEval.RCBoot(xmat, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=blocklengthvec[nn])
			@test rc1.alpha == 0.05
			@test rc1.bootinput.bootmethod == bootmethodtypevec[kk]
			if blocklengthvec[nn] <= 0.0 ; @test rc1.bootinput.blocklength > 0.0
			else
				if bootmethodsymvec[kk] == :iid ; @test rc1.bootinput.blocklength == 1.0
				else @test rc1.bootinput.blocklength == blocklengthvec[nn]
				end
			end
		end
	end
	#Test main routine
	Random.seed!(1234)
	pvalmat = [0.868  0.872 ; 0.879  0.873];
	for kk = 1:length(bootmethodsymvec)
		for nn = 1:length(blocklengthvec)
			rcout = rc(xmat, :boot, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=blocklengthvec[nn])
			@test isapprox(pvalmat[kk,nn], rcout.pvalue)
		end
	end
end

@testset "SPA tests" begin
	#Test constructor
	for kk = 1:length(bootmethodsymvec)
		for nn = 1:length(blocklengthvec)
			for mm = 1:length(kernelfuncsymvec)
				for jj = 1:length(bandwidthvec)
					spa1 = ForecastEval.SPABoot(xmat, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=blocklengthvec[nn], kernelfunction=kernelfuncsymvec[mm], bandwidth=bandwidthvec[jj])
					@test spa1.alpha == 0.05
					@test spa1.bootinput.bootmethod == bootmethodtypevec[kk]
					if blocklengthvec[nn] <= 0.0 ; @test spa1.bootinput.blocklength > 0.0
					else
						if bootmethodsymvec[kk] == :iid ; @test spa1.bootinput.blocklength == 1.0
						else @test spa1.bootinput.blocklength == blocklengthvec[nn]
						end
					end
					@test spa1.kernelfunction == kernelfunctypevec[mm]
					@test spa1.bandwidth == bandwidthvec[jj]
				end
			end
		end
	end
    #Test routine (note, hacvariance has already been fully covered by previous tests, as has bootinput constructor)
    Random.seed!(1234)
    spaout = spa(xmat, alpha=0.05, bootmethod=bootmethodsymvec[2], blocklength=blocklengthvec[1], kernelfunction=kernelfuncsymvec[1], bandwidth=bandwidthvec[1])
    @test isapprox(spaout.pvalue_u, 0.878)
    @test isapprox(spaout.pvalue_c, 0.845)
    @test isapprox(spaout.pvalue_l, 0.783)
    @test isapprox(spaout.pvalueauto, 0.845)
    @test isapprox(spaout.teststat, 0.2827379497687274)
    @test isapprox(sum(spaout.mu_u), -0.42402380049999994)
    @test isapprox(sum(spaout.mu_c), -0.05397999049999995)
    @test isapprox(sum(spaout.mu_l), 0.06876765000000001)
end

@testset "MCS tests" begin
    #Test constructor
    for kk = 1:length(bootmethodsymvec)
        for nn = 1:length(blocklengthvec)
            mca1 = ForecastEval.MCSBoot(xmatabs, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=blocklengthvec[nn], basecaseindex=2)
            @test mca1.alpha == 0.05
            @test mca1.bootinput.bootmethod == bootmethodtypevec[kk]
            if blocklengthvec[nn] <= 0.0 ; @test mca1.bootinput.blocklength > 0.0
            else
                if bootmethodsymvec[kk] == :iid ; @test mca1.bootinput.blocklength == 1.0
                else @test mca1.bootinput.blocklength == blocklengthvec[nn]
                end
            end
            @test mca1.basecaseindex == 2
        end
    end
    for kk = 1:length(bootmethodsymvec)
        for nn = 1:length(blocklengthvec)
            mca1 = ForecastEval.MCSBootLowRAM(xmatabs, alpha=0.05, bootmethod=bootmethodsymvec[kk], blocklength=blocklengthvec[nn], basecaseindex=2)
            @test mca1.alpha == 0.05
            @test mca1.bootinput.bootmethod == bootmethodtypevec[kk]
            if blocklengthvec[nn] <= 0.0 ; @test mca1.bootinput.blocklength > 0.0
            else
                if bootmethodsymvec[kk] == :iid ; @test mca1.bootinput.blocklength == 1.0
                else @test mca1.bootinput.blocklength == blocklengthvec[nn]
                end
            end
            @test mca1.basecaseindex == 2
        end
    end
    #Test routine (note, bootinput constructor has already been fully covered by previous tests)
	Random.seed!(1234)
    mcsout1 = mcs(xmatabs, alpha=0.05, bootmethod=bootmethodsymvec[2], blocklength=blocklengthvec[1], basecaseindex=2)
    @test mcsout1.inQF == Int[3,2,4,1]
    @test mcsout1.outQF == Int[5]
    @test isapprox(sum(mcsout1.pvalueQF), 2.162)
    @test mcsout1.inMT == Int[5,3,2,4,1]
    @test mcsout1.outMT == Int[]
    @test isapprox(sum(mcsout1.pvalueMT), 2.2039999999999997)
    Random.seed!(1234)
    mcsout2 = mcs(xmatabs, ForecastEval.MCSBootLowRAM(xmatabs, alpha=0.05, bootmethod=bootmethodsymvec[2], blocklength=blocklengthvec[1], basecaseindex=2))
    @test mcsout2.inQF == Int[3,2,4,1]
    @test mcsout2.outQF == Int[5]
    @test isapprox(sum(mcsout2.pvalueQF), 2.162)
    @test mcsout2.inMT == Int[5,3,2,4,1]
    @test mcsout2.outMT == Int[]
    @test isapprox(sum(mcsout2.pvalueMT), 2.2039999999999997)
end


text:
relu(x::AbstractArray) = max.(x, 0)
relu_grad(x::AbstractArray) = float(x .> 0)  # TODO: ensure  relu_grad works for CuArray

@diffrule relu(x::AbstractArray) x relu_grad(x)


text:
using JuMP
using Clp
using DiffOpt
using Test
using ChainRulesCore

# This script creates the JuMP model for a small unit commitment instance
# represented in a solution map function taking parameters as arguments and returning
# the optimal solution as output.

# The derivatives of this solution map can then be expressed in ChainRules semantics
# and implemented using DiffOpt

ATOL=1e-4
RTOL=1e-4

"""
Solution map of the problem using parameters:
- `load1_demand, load2_demand` for the demand of two nodes
- `gen_costs` is the vector of generator costs
- `noload_costs` is the vector of fixed activation costs of the generators,
and returning the optimal output power `p`.
"""
function unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer)))
    ## Problem data
    unit_codes = [1, 2] # Generator identifiers
    load_names = ["Load1", "Load2"] # Load identifiers
    n_periods = 4 # Number of time periods
    Pmin = Dict(1 => fill(0.5, n_periods), 2 => fill(0.5, n_periods)) # Minimum power output (pu)
    Pmax = Dict(1 => fill(3.0, n_periods), 2 => fill(3.0, n_periods)) # Maximum power output (pu)
    RR = Dict(1 => 0.25, 2 => 0.25) # Ramp rates (pu/min)
    P0 = Dict(1 => 0.0, 2 => 0.0) # Initial power output (pu)
    D = Dict("Load1" => load1_demand, "Load2" => load2_demand) # Demand (pu)
    Cp = Dict(1 => gen_costs[1], 2 => gen_costs[2]) # Generation cost coefficient ($/pu)
    Cnl = Dict(1 => noload_costs[1], 2 => noload_costs[2]) # No-load cost ($)

    ## Variables
    # Note: u represents the activation of generation units.
    # Would be binary in the typical UC problem, relaxed here to u ∈ [0,1]
    # for a linear relaxation.
    @variable(model, 0 <= u[g in unit_codes, t in 1:n_periods] <= 1) # Commitment
    @variable(model, p[g in unit_codes, t in 1:n_periods] >= 0) # Power output (pu)

    ## Constraints

    # Energy balance
    @constraint(
        model,
        energy_balance_cons[t in 1:n_periods],
        sum(p[g, t] for g in unit_codes) == sum(D[l][t] for l in load_names),
    )

    # Generation limits
    @constraint(model, [g in unit_codes, t in 1:n_periods], Pmin[g][t] * u[g, t] <= p[g, t])
    @constraint(model, [g in unit_codes, t in 1:n_periods], p[g, t] <= Pmax[g][t] * u[g, t])

    # Ramp rates
    @constraint(model, [g in unit_codes, t in 2:n_periods], p[g, t] - p[g, t - 1] <= 60 * RR[g])
    @constraint(model, [g in unit_codes], p[g, 1] - P0[g] <= 60 * RR[g])
    @constraint(model, [g in unit_codes, t in 2:n_periods], p[g, t - 1] - p[g, t] <= 60 * RR[g])
    @constraint(model, [g in unit_codes], P0[g] - p[g, 1] <= 60 * RR[g])

    # Objective
    @objective(
        model,
        Min,
        sum((Cp[g] * p[g, t]) + (Cnl[g] * u[g, t]) for g in unit_codes, t in 1:n_periods),
    )

    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    # converting to dense matrix
    return JuMP.value.(p.data)
end

@show unit_commitment([1.0, 1.2, 1.4, 1.6], [1.0, 1.2, 1.4, 1.6], [1000.0, 1500.0], [500.0, 1000.0])

# Forward differentiation rule for the solution map of the unit commitment problem
# taking in input perturbations on the input parameters and returning perturbations propagated to the result
function ChainRulesCore.frule((_, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), ::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs)
    # creating the UC model with a DiffOpt optimizer wrapper around Clp
    model = Model(() -> diff_optimizer(Clp.Optimizer))
    # building and solving the main model
    pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model)
    energy_balance_cons = model[:energy_balance_cons]

    # Setting some perturbation of the right-hand side of the energy balance constraints
    # the RHS is equal to the sum of load demands at each period.
    # the corresponding perturbation are set accordingly as the set of perturbations of the two loads
    MOI.set.(
        model,
        DiffOpt.ForwardInConstraint(), energy_balance_cons,
        AffExpr[d1 + d2 for (d1, d2) in zip(Δload1_demand, Δload2_demand)],
    )

    p = model[:p]
    u = model[:u]

    # setting the perturbation of the linear objective
    Δobj = sum(Δgen_costs ⋅ p[:,t] + Δnoload_costs ⋅ u[:,t] for t in size(p, 2))
    MOI.set(model, DiffOpt.ForwardInObjective(), Δobj)
    DiffOpt.forward(JuMP.backend(model))
    # querying the corresponding perturbation of the decision
    Δp = MOI.get.(model, DiffOpt.ForwardOutVariablePrimal(), p)
    return (pv, Δp.data)
end

load1_demand = [1.0, 1.2, 1.4, 1.6]
load2_demand = [1.0, 1.2, 1.4, 1.6]
gen_costs = [1000.0, 1500.0]
noload_costs = [500.0, 1000.0]

Δload1_demand = 0 * load1_demand .+ 0.1
Δload2_demand = 0 * load2_demand .+ 0.2
Δgen_costs = 0 * gen_costs .+ 0.1
Δnoload_costs = 0 * noload_costs .+ 0.4
@show (pv, Δpv) = ChainRulesCore.frule((nothing, Δload1_demand, Δload2_demand, Δgen_costs, Δnoload_costs), unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs)

# Reverse-mode differentiation of the solution map
# The computed pullback takes a seed for the optimal solution `̄p` and returns
# derivatives wrt each input parameter.
function ChainRulesCore.rrule(::typeof(unit_commitment), load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer)))
    # solve the forward UC problem
    pv = unit_commitment(load1_demand, load2_demand, gen_costs, noload_costs, model=model)
    function pullback_unit_commitment(pb)
        p = model[:p]
        u = model[:u]
        energy_balance_cons = model[:energy_balance_cons]

        MOI.set.(model, DiffOpt.BackwardInVariablePrimal(), p, pb)
        DiffOpt.backward(JuMP.backend(model))

        obj = MOI.get(model, DiffOpt.BackwardOutObjective())

        # computing derivative wrt linear objective costs
        dgen_costs = similar(gen_costs)
        dgen_costs[1] = sum(JuMP.coefficient.(obj, p[1,:]))
        dgen_costs[2] = sum(JuMP.coefficient.(obj, p[2,:]))

        dnoload_costs = similar(noload_costs)
        dnoload_costs[1] = sum(JuMP.coefficient.(obj, u[1,:]))
        dnoload_costs[2] = sum(JuMP.coefficient.(obj, u[2,:]))

        # computing derivative wrt constraint constant
        dload1_demand = JuMP.constant.(MOI.get.(model, DiffOpt.BackwardOutConstraint(), energy_balance_cons))
        dload2_demand = copy(dload1_demand)
        return (dload1_demand, dload2_demand, dgen_costs, dnoload_costs)
    end
    return (pv, pullback_unit_commitment)
end

(pv, pullback_unit_commitment) = ChainRulesCore.rrule(unit_commitment, load1_demand, load2_demand, gen_costs, noload_costs; model = Model(() -> diff_optimizer(Clp.Optimizer)))
@show pullback_unit_commitment(ones(size(pv)))


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/evaluator/passes.jl
# Functions used to compute forward pass of nonlinear functions which include:
# set_value_post, overwrite_or_intersect, forward_pass_kernel, associated blocks
#############################################################################

xnum_yset(b, x, y) = is_num(b, x) && !is_num(b, y)
xset_ynum(b, x, y) = !is_num(b, x) && is_num(b, y)
xy_num(b, x, y) = is_num(b, x) && is_num(b, y)
xyset(b, x, y) = !(is_num(b, x) || is_num(b, y))

function varset(::Type{MC{N,T}}, i, x_cv, x_cc, l, u) where {V,N,T<:RelaxTag}
    v = seed_gradient(i, Val(N))
    return MC{N,T}(x_cv, x_cc, Interval{Float64}(l, u), v, v, false)
end

function fprop!(t::RelaxCacheAttribute, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    if l == u
        b[k] = x
    else
        z = varset(MC{N,T}, rev_sparsity(g, i, k), x, x, l, u)
        if !first_eval(t, b)
            z = z ∩ interval(b, k)
        end
        b[k] = z
    end
    nothing
end

function fprop!(t::Relax, ex::Subexpression, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    y = first_index(g, k)
    x = dependent_subexpression_index(g, y)
    if subexpression_is_num(b, x)
        b[k] = subexpression_num(b, x)
    else
        b[k] = subexpression_set(b, x)
    end
end

for (F, f) in ((PLUS, :+), (MIN, :min), (MAX, :max), (DIV, :/), (ARH, :arh))
    @eval function fprop_2!(t::Relax, v::Val{$F}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        x = child(g, 1, k)
        y = child(g, 2, k)
        if !xy_num(b, x, y)
            if xyset(b, x, y)
                z = ($f)(set(b, x), set(b, y))
            elseif xset_ynum(b, x, y)
                z = ($f)(set(b, x), num(b, y))
            else
                z = ($f)(num(b, x), set(b, y))
            end
            b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        else
            b[k] = ($f)(num(b, x), num(b, y))
        end
    end
end

function fprop!(t::Relax, v::Val{MINUS}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    if is_binary(g, k)
        y = child(g, 2, k)
        if !xy_num(b, x, y)
            if xyset(b, x, y)
                z = set(b, x) - set(b, y)
            elseif xset_ynum(b, x, y)
                z = set(b, x) - num(b, y)
            else
                z = num(b, x) - set(b, y)
            end
            b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        else
            b[k] = num(b, x) - num(b, y)
        end
    else
        if is_num(b, x)
            b[k] = -num(b, x)
        else
            z = -set(b, x)
            b[k] = cut(-set(b, x), set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        end
    end
end

function fprop_n!(t::Relax, v::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = zero(MC{N,T})
    znum = 0.0
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum += num(b, i)
        else
            numval = false
            z += set(b, i)
        end
    end
    if numval
        b[k] = znum
    else
        z += znum
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

function fprop_n!(t::Relax, v::Val{MIN}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = Inf*one(MC{N,T})
    znum = Inf
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum = min(znum, num(b, i))
        else
            numval = false
            z = min(z, set(b, i))
        end
    end
    if numval
        b[k] = znum
    else
        z = min(z, znum)
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

function fprop_n!(t::Relax, v::Val{MAX}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = -Inf*one(MC{N,T})
    znum = -Inf
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum = max(znum, num(b, i))
        else
            numval = false
            z = max(z, set(b, i))
        end
    end
    if numval
        b[k] = znum
    else
        z = max(z, znum)
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

function fprop_2!(t::Relax, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}

    x = child(g, 1, k)
    y = child(g, 2, k)

    if !xy_num(b, x, y)
        if xyset(b, x, y)
            xv = set(b, x)
            yv = set(b, y)
            if b.use_apriori_mul
                dp = b.dp
                dP = b.dP
                p_rel = b.p_rel
                p_diam = b.p_diam
                s = sparsity(g, 1)
                xr = info(b, x)
                yr = info(b, y)
                u1max, u2max, v1nmax, v2nmax = estimator_extrema(xr, yr, s, dP)
                z = xv*yv
                wIntv = z.Intv
                if (u1max < xv.Intv.hi) || (u2max < yv.Intv.hi)
                    u1cv, u2cv, u1cvg, u2cvg = estimator_under(xv, yv, xr, yr, s, dp, dP, p_rel, p_diam)
                    za_l = McCormick.mult_apriori_kernel(xv, yv, wIntv, u1cv, u2cv, u1max, u2max, u1cvg, u2cvg)
                    z = z ∩ za_l
                end
                if (v1nmax > -xv.Intv.lo) || (v2nmax > -yv.Intv.lo)
                    v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_over(xv, yv, xr, yr, s, dp, dP, p_rel, p_diam)
                    za_u = McCormick.mult_apriori_kernel(-xv, -yv, wIntv, v1ccn, v2ccn, v1nmax, v2nmax, v1ccgn, v2ccgn)
                    z = z ∩ za_u
                end
            else
                z = xv*yv
            end
        elseif xset_ynum(b, x, y)
            z = set(b, x)*num(b, y)
        else
            z = num(b, x)*set(b, y)
        end
        b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    else
        b[k] = num(b, x)*num(b, y)
    end
end

function fprop_n!(t::Relax, ::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}

    z = one(MC{N,T})
    znum = one(Float64)
    numval = true
    if b.use_apriori_mul
        zr = one(V)
        dp = b.dp
        dP = b.dP
        s = sparsity(g, 1)
        for (q,i) in enumerate(children(g, k))
            if is_num(b, i)
                znum = znum*num(b, i)
            else
                numval = false
                x = set(b, i)
                xr = info(b, i)
                u1max, u2max, v1nmax, v2nmax = estimator_extrema(zr, xr, s, dP)
                zv = z*x
                wIntv = zv.Intv
                if (u1max < z.Intv.hi) || (u2max < x.Intv.hi)
                    u1cv, u2cv, u1cvg, u2cvg = estimator_under(zr, xr, s, dp, dP)
                    za_l = McCormick.mult_apriori_kernel(z, x, wIntv, u1cv, u2cv, u1max, u2max, u1cvg, u2cvg)
                    zv = zv ∩ za_l
                end
                if (v1nmax > -z.Intv.lo) || (v2nmax > -x.Intv.lo)
                    v1ccn, v2ccn, v1ccgn, v2ccgn = estimator_under(zr, xr, s, dp, dP)
                    za_u = McCormick.mult_apriori_kernel(-z, -x, wIntv, v1ccn, v2ccn, v1nmax, v2nmax, v1ccgn, v2ccgn)
                    zv = zv ∩ za_u
                end
                zr = zr*xr
                zv = cut(zv, zv, b.ic.v, b.ϵ_sg, sparsity(g, i), b.cut, false)
                z = zv
            end
        end
    else
        for (q,i) in enumerate(children(g, k))
            if is_num(b, i)
                znum = znum*num(b, i)
            else
                numval = false
                x = set(b, i)
                z = z*x
                z = cut(z, z, b.ic.v, b.ϵ_sg, sparsity(g, i), b.cut, false)
            end
        end
    end
    if numval
        b[k] = znum
    else
        z = z*znum
        b[k] = z
    end
end

for F in (PLUS, MULT, MIN, MAX)
    @eval function fprop!(t::Relax, v::Val{$F}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        is_binary(g, k) ? fprop_2!(Relax(), Val($F), g, b, k) : fprop_n!(Relax(), Val($F), g, b, k)
    end
end
function fprop!(t::Relax, v::Val{DIV}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    fprop_2!(Relax(), v, g, b, k)
end

function fprop!(t::Relax, v::Val{POW}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    if is_num(b, y) && isone(num(b, y))
        b[k] = set(b, x)
    elseif is_num(b,y) && iszero(num(b, y))
        b[k] = one(MC{N,T})
    elseif !xy_num(b, x, y)
        if xyset(b, x, y)
            z = set(b, x)^set(b, y)
        elseif xset_ynum(b, x, y)
            z = set(b, x)^num(b, y)
        else
            z = num(b, x)^set(b, y)
        end
        b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    else
        b[k] = num(b, x)^num(b, y)
    end
end

function fprop!(t::Relax, v::Val{USER}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    f = user_univariate_operator(g, first_index(g, k))
    x = child(g, 1, k)
    if is_num(b, x)
        b[k] = f(num(b, x))
    else
        z = f(set(b, x))
        b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g, k), b.cut, b.post)
    end
end

function fprop!(t::Relax, v::Val{USERN}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    mv = user_multivariate_operator(g, first_index(g, k))
    n = arity(g, k)
    set_input = _set_input(b, n)
    num_input = _num_input(b, n)
    anysets = false
    i = 1
    for c in children(g, k)
        if is_num(b, c)
            x = num(b, c)
            if !isinf(x)
                set_input[i] = MC{N,T}(x)
                num_input[i] = x
            end
        else
            set_input[i] = set(b, c)
            anysets = true
        end
        i += 1
    end
    if anysets
        z = MOI.eval_objective(mv, set_input)::MC{N,T}
        b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
    else
        b[k] = MOI.eval_objective(mv, num_input)
    end
end

for ft in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[ft]
    (f == :user || f == :+ || f == :-) && continue
    @eval function fprop!(t::Relax, v::Val{$ft}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        x = child(g, 1, k)
        if is_num(b, x)
            return b[k] = ($f)(num(b, x))
        else
            z = ($f)(set(b, x))
            b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
        end
    end
end

for (F, f) in ((LOWER_BND, :lower_bnd), (UPPER_BND, :upper_bnd))
    @eval function fprop!(t::Relax, v::Val{$F}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        y = child(g, 2, k)
        if is_num(b, y)
            z = set(b, child(g, 1, k))
            z = ($f)(z, num(b, y))
            b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g, k), b.cut, b.post)
        end
    end
end

function fprop!(t::Relax, v::Val{BND}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    z = set(b, child(g, 1, k))
    y = child(g, 2, k)
    r = child(g, 3, k)
    if is_num(b, y) && is_num(b, r)
        z = bnd(z, num(b, y),num(b, r))
    end
    b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
end

function f_init!(t::Relax, g::DAT, b::RelaxCache)
    for k = node_count(g):-1:1
        c = node_class(g, k)
        (c == EXPRESSION)    && fprop!(t, Expression(), g, b, k)
        (c == VARIABLE)      && fprop!(t, Variable(), g, b, k)
        (c == SUBEXPRESSION) && fprop!(t, Subexpression(), g, b, k)
        b._info[k] = set(b, k)
    end
    nothing
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
import Automa
import Automa.RegExp: @re_str

using Test

@testset "SizedMemory" begin
    # Vector{UInt8}
    mem = Automa.SizedMemory(b"bar")
    @test lastindex(mem) === length(mem) === 3
    @test mem[1] === UInt8('b')
    @test mem[2] === UInt8('a')
    @test mem[3] === UInt8('r')
    @test_throws BoundsError mem[0]
    @test_throws BoundsError mem[4]

    # String
    mem = Automa.SizedMemory("bar")
    @test lastindex(mem) === length(mem) === 3
    @test mem[1] === UInt8('b')
    @test mem[2] === UInt8('a')
    @test mem[3] === UInt8('r')
    @test_throws BoundsError mem[0]
    @test_throws BoundsError mem[4]

    # SubString
    mem = Automa.SizedMemory(SubString("xbar", 2, 4))
    @test lastindex(mem) === length(mem) === 3
    @test mem[1] === UInt8('b')
    @test mem[2] === UInt8('a')
    @test mem[3] === UInt8('r')
    @test_throws BoundsError mem[0]
    @test_throws BoundsError mem[4]
end

@testset "RegExp" begin
    @test_throws ArgumentError("invalid escape sequence: \\o") Automa.RegExp.parse("\\o")
end

@testset "DOT" begin
    re = re"[A-Za-z_][A-Za-z0-9_]*"
    re.actions[:enter] = [:enter]
    re.actions[:exit]  = [:exit]
    nfa = Automa.re2nfa(re)
    @test startswith(Automa.nfa2dot(nfa), "digraph")
    dfa = Automa.nfa2dfa(nfa)
    @test startswith(Automa.dfa2dot(dfa), "digraph")
    machine = Automa.compile(re)
    @test startswith(Automa.machine2dot(machine), "digraph")
    @test occursin(r"^Automa\.NFANode\(.*\)$", repr(nfa.start))
    @test occursin(r"^Automa\.DFANode\(.*\)$", repr(dfa.start))
    @test occursin(r"^Automa\.Node\(.*\)$", repr(machine.start))
end

@testset "Determinacy" begin
    # see https://github.com/BioJulia/Automa.jl/issues/19
    notmach(re) = Automa.machine2dot(Automa.compile(re)) != Automa.machine2dot(Automa.compile(re))
    for re in (re"0?11|0?12", re"0?12|0?1*")
        @test count(_->notmach(re), 1:1000) == 0
    end
end

include("test01.jl")
include("test02.jl")
include("test03.jl")
include("test04.jl")
include("test05.jl")
include("test06.jl")
include("test07.jl")
include("test08.jl")
include("test09.jl")
include("test10.jl")
include("test11.jl")
include("test12.jl")
include("test13.jl")
include("test14.jl")
include("test15.jl")
include("test16.jl")
include("test17.jl")
include("test18.jl")

module TestFASTA
using Test
@testset "FASTA" begin
    include("../example/fasta.jl")
    @test records[1].identifier == "NP_003172.1"
    @test records[1].description == "brachyury protein isoform 1 [Homo sapiens]"
    @test records[1].sequence[1:5] == b"MSSPG"
    @test records[1].sequence[end-4:end] == b"SPPSM"
end
end

module TestNumbers
using Test
@testset "Numbers" begin
    include("../example/numbers.jl")
    @test tokens == [(:dec,"1"),(:hex,"0x0123BEEF"),(:oct,"0o754"),(:float,"3.14"),(:float,"-1e4"),(:float,"+6.022045e23")]
    @test status == :ok
    @test startswith(Automa.machine2dot(machine), "digraph")

    notmach(re) = Automa.machine2dot(Automa.compile(re)) != Automa.machine2dot(Automa.compile(re))
    @test count(_->notmach(numbers), 1:15) == 0
end
end

module TestTokenizer
using Test
@testset "MiniJulia" begin
    include("../example/tokenizer.jl")
    @test tokens[1:14] == [
        (:identifier,"quicksort"),
        (:lparen,"("),
        (:identifier,"xs"),
        (:rparen,")"),
        (:spaces," "),
        (:equal,"="),
        (:spaces," "),
        (:identifier,"quicksort!"),
        (:lparen,"("),
        (:identifier,"copy"),
        (:lparen,"("),
        (:identifier,"xs"),
        (:rparen,")"),
        (:rparen,")")]
    @test tokens[end-5:end] == [
        (:keyword,"return"),
        (:spaces," "),
        (:identifier,"j"),
        (:newline,"\n"),
        (:keyword,"end"),
        (:newline,"\n")]
end
end

module TestStream

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos
using TranscodingStreams
using Test

# Test 1
machine = let
    line = re"[^\r\n]*"
    newline = re"\r?\n"
    Automa.compile(line * newline)
end
Automa.Stream.generate_reader(:readline, machine) |> eval
@testset "Scanning a line" begin
    for (data, state) in [
            ("\n",      :ok),
            ("foo\n",   :ok),
            ("foo\r\n", :ok),
            ("",        :incomplete),
            ("\r",      :incomplete),
            ("foo",     :incomplete),
            ("\r\r",    :error),
            ("\r\nx",   :error),
            ("foo\nx",  :error),]
        s = readline(NoopStream(IOBuffer(data)))
        if state == :ok
            @test s == 0
        elseif state == :incomplete
            @test s > 0
        else
            @test s < 0
        end
    end
end

# Test 2
machine = let
    alphanum = re"[A-Za-z0-9]+"
    alphanum.actions[:enter] = [:start_alphanum]
    alphanum.actions[:exit]  = [:end_alphanum]
    whitespace = re"[ \t\r\n]*"
    Automa.compile(whitespace * alphanum * whitespace)
end
actions = Dict(
   :start_alphanum => :(@mark; start_alphanum = @relpos(p)),
   :end_alphanum   => :(end_alphanum = @relpos(p-1)),
)
initcode = :(start_alphanum = end_alphanum = 0)
returncode = :(return cs == 0 ? String(data[@abspos(start_alphanum):@abspos(end_alphanum)]) : "")
Automa.Stream.generate_reader(:stripwhitespace, machine, actions=actions, initcode=initcode, returncode=returncode) |> eval
@testset "Stripping whitespace" begin
    for (data, value) in [
            ("x", "x"),
            (" foo ", "foo"),
            ("  \r\n123\n  ", "123"),
            ("   abc123   ", "abc123"),
            ("", ""),
            ("  12+3 ", ""),
           ]
        for bufsize in [1:5; 100]
            @test stripwhitespace(NoopStream(IOBuffer(data), bufsize=bufsize)) == value
        end
    end
end

# Three-column BED file format.
cat = Automa.RegExp.cat
rep = Automa.RegExp.rep
machine = let
    chrom = re"[^\t]+"
    chrom.actions[:exit] = [:chrom]
    chromstart = re"[0-9]+"
    chromstart.actions[:exit] = [:chromstart]
    chromend = re"[0-9]+"
    chromend.actions[:exit] = [:chromend]
    record = cat(chrom, '\t', chromstart, '\t', chromend)
    record.actions[:enter] = [:mark]
    bed = rep(cat(record, re"\r?\n"))
    Automa.compile(bed)
end
#write("bed.dot", Automa.machine2dot(machine))
#run(`dot -Tsvg -o bed.svg bed.dot`)
actions = Dict(
    :mark => :(@mark; mark = @relpos(p)),
    :chrom => :(chrom = @relpos(p-1)),
    :chromstart => :(chromstart = @relpos(p-1)),
    :chromend => :(chromend = @relpos(p-1); found = true; @escape)
)
initcode = :(cs = state; mark = chrom = chromstart = chromend = 0; found = false)
loopcode = :(found && @goto __return__)
returncode = quote
    if found
        return (
            String(data[@abspos(mark):@abspos(chrom)]),
            parse(Int, String(data[@abspos(chrom)+2:@abspos(chromstart)])),
            parse(Int, String(data[@abspos(chromstart)+2:@abspos(chromend)]))),
            cs
    else
        return ("", -1, 0), cs
    end
end
Automa.Stream.generate_reader(:readrecord!, machine, arguments=(:(state::Int),), actions=actions, initcode=initcode, loopcode=loopcode, returncode=returncode) |> eval

@testset "Three-column BED (stateful)" begin
    stream = NoopStream(IOBuffer("""chr1\t10\t200\n"""))
    state = machine.start_state
    val, state = readrecord!(stream, state)
    @test val == ("chr1", 10, 200)
    val, state = readrecord!(stream, state)
    @test val == ("", -1, 0)
    @test state == 0
    stream = NoopStream(IOBuffer("""1\t10\t200000\nchr12\t0\t21000\r\nchrM\t123\t12345\n"""))
    state = machine.start_state
    val, state = readrecord!(stream, state)
    @test val == ("1", 10, 200000)
    val, state = readrecord!(stream, state)
    @test val == ("chr12", 0, 21000)
    val, state = readrecord!(stream, state)
    @test val == ("chrM", 123, 12345)
    val, state = readrecord!(stream, state)
    @test val == ("", -1, 0)
    @test state == 0
end

# FASTA
mutable struct Record
    data::Vector{UInt8}
    identifier::UnitRange{Int}
    description::UnitRange{Int}
    sequence::UnitRange{Int}
end

function Record()
    return Record(UInt8[], 1:0, 1:0, 1:0)
end

function initialize!(record::Record)
    empty!(record.data)
    record.identifier = record.description = record.sequence = 1:0
    return record
end

machine = let re = Automa.RegExp
    newline = re"\r?\n"

    identifier = re"[!-~]*"
    identifier.actions[:enter] = [:pos]
    identifier.actions[:exit] = [:identifier]

    description = re"[!-~][ -~]*"
    description.actions[:enter] = [:pos]
    description.actions[:exit] = [:description]

    header = re.cat('>', identifier, re.opt(re" " * description))
    header.actions[:exit] = [:header]

    letters = re"[A-Za-z*-]*"
    letters.actions[:enter] = [:mark, :pos]
    letters.actions[:exit] = [:letters]

    sequence = re.cat(letters, re.rep(newline * letters))

    record = re.cat(header, newline, sequence)
    record.actions[:enter] = [:mark]
    record.actions[:exit] = [:record]

    fasta = re.rep(record)

    Automa.compile(fasta)
end

actions = Dict(
    :mark => :(@mark),
    :pos => :(pos = @relpos(p)),
    :identifier => :(record.identifier = pos:@relpos(p-1)),
    :description => :(record.description = pos:@relpos(p-1)),
    :header => :(append!(record.data, data[@markpos():p-1]); push!(record.data, UInt8('\n'))),
    :letters => quote
        if isempty(record.sequence)
            record.sequence = length(record.data)+1:length(record.data)
        end
        record.sequence = first(record.sequence):last(record.sequence)+p-@abspos(pos)
        append!(record.data, data[@abspos(pos):p-1])
    end,
    :record => quote
        found = true
        @escape
    end
)
initcode = quote
    cs = state
    pos = 0
    found = false
    initialize!(record)
end
loopcode = quote
    found && @goto __return__
end
context = Automa.CodeGenContext(generator=:goto)
Automa.Stream.generate_reader(
    :readrecord!,
    machine,
    arguments=(:(state::Int), :(record::$(Record)),),
    actions=actions,
    context=context,
    initcode=initcode,
    loopcode=loopcode,
) |> eval

@testset "Streaming FASTA" begin
    stream = NoopStream(IOBuffer(""))
    state = machine.start_state
    record = Record()
    @test readrecord!(stream, state, record) == 0

    stream = NoopStream(IOBuffer("""
    >seq1 hogehoge
    ACGT
    TGCA
    """), bufsize=10)
    state = machine.start_state
    record = Record()
    @test readrecord!(stream, state, record) == 0
    @test String(record.data[record.identifier]) == "seq1"
    @test String(record.data[record.description]) == "hogehoge"
    @test String(record.data[record.sequence]) == "ACGTTGCA"

    stream = NoopStream(IOBuffer("""
    >seq1 1st sequence
    NANANANA
    >seq2 2nd sequence
    -----AAA
    GGGGG---
    """), bufsize=10)
    state = machine.start_state
    record = Record()
    state = readrecord!(stream, state, record)
    @test state > 0
    @test String(record.data[record.identifier]) == "seq1"
    @test String(record.data[record.description]) == "1st sequence"
    @test String(record.data[record.sequence]) == "NANANANA"
    state = readrecord!(stream, state, record)
    @test state == 0
    @test String(record.data[record.identifier]) == "seq2"
    @test String(record.data[record.description]) == "2nd sequence"
    @test String(record.data[record.sequence]) == "-----AAAGGGGG---"

    stream = NoopStream(IOBuffer("""
    >seq1 1st sequence

    N
    ANANANA
    >seq2 2nd sequence
    -----AAA

    GGGG

    G---

    """), bufsize=10)
    state = machine.start_state
    record = Record()
    state = readrecord!(stream, state, record)
    @test state > 0
    @test String(record.data[record.identifier]) == "seq1"
    @test String(record.data[record.description]) == "1st sequence"
    @test String(record.data[record.sequence]) == "NANANANA"
    state = readrecord!(stream, state, record)
    @test state == 0
    @test String(record.data[record.identifier]) == "seq2"
    @test String(record.data[record.description]) == "2nd sequence"
    @test String(record.data[record.sequence]) == "-----AAAGGGGG---"
end

end


text:
"""
    TrackGenerator{M<:Mesh,Q<:AzimuthalQuadrature,T<:Real}

This is the main structure of the library, which holds all the information about the ray
tracing.
"""
struct TrackGenerator{T<:Real,M<:Mesh,BC<:BoundaryConditions,Q<:AzimuthalQuadrature}
    mesh::M
    bcs::BC

    azimuthal_quadrature::Q
    n_tracks_x::Vector{Int}
    n_tracks_y::Vector{Int}
    n_tracks::Vector{Int}
    n_total_tracks::Int

    tracks::Vector{Vector{Track}}
    tracks_by_uid::Vector{Track}

    tiny_step::T
    correct_volumes::Bool
    volumes::Vector{T}
end

function show(io::IO, t::TrackGenerator)
    @unpack azimuthal_quadrature, n_total_tracks, correct_volumes = t
    @unpack δs, ϕs = azimuthal_quadrature
    n_azim_2 = nazim2(azimuthal_quadrature)

    # println(io, typeof(t))
    println(io, "  Number of azimuthal angles in (0, π): ", n_azim_2)
    println(io, "  Azimuthal angles in (0, π): ", round.(rad2deg.(ϕs), digits=2))
    println(io, "  Effective azimuthal spacings: ", round.(δs, digits=3))
    println(io, "  Total tracks: ", n_total_tracks)
    print(io,   "  Correct volumes: ", correct_volumes)
end

origins_in_x(n_tracks_x, i, j) = j ≤ n_tracks_x[i]
origins_in_y(n_tracks_x, i, j) = !origins_in_x(n_tracks_x, i, j)

"""
    TrackGenerator(
        model::UnstructuredDiscreteModel, n_azim::Int, δ::T;
        tiny_step::T=1e-8, correct_volumes=false
    ) where {T<:Real}

Initialize a [`TrackGenerator`](@ref) using an [`UnstructuredDiscreteModel`](@ref) which
holds mesh related information, the number of azimuthal angles `n_azim` and the azimuthal
spacing `δ`. The optional attributes are `tiny_step`, which is used in the track's
segmentation routine, and `correct_volumes`, which allows for cell volume correction since
the ray tracing algorithm yields to approximate volumes.
"""
function TrackGenerator(
    model::UnstructuredDiscreteModel, n_azim::Int, δ::T;
    bcs::BoundaryConditions=BoundaryConditions(), tiny_step::T=1e-8, correct_volumes=false
) where {T<:Real}

    mesh = Mesh(model)
    Δx, Δy = width(mesh), height(mesh)

    azimuthal_quadrature = AzimuthalQuadrature(Val(n_azim), δ)
    n_azim_2 = nazim2(azimuthal_quadrature)
    n_azim_4 = nazim4(azimuthal_quadrature)

    n_tracks_x = Vector{Int}(undef, n_azim_2)
    n_tracks_y = Vector{Int}(undef, n_azim_2)
    n_tracks   = Vector{Int}(undef, n_azim_2)

    for i in right_dir(azimuthal_quadrature)
        φ = π / n_azim_2 * (i - 1 / 2)

        n_tracks_x[i] = floor(Δx / δ * abs(sin(φ))) + 1
        n_tracks_y[i] = floor(Δy / δ * abs(cos(φ))) + 1
        n_tracks[i] = n_tracks_x[i] + n_tracks_y[i]

        # supplementaries:
        j = suplementary_idx(azimuthal_quadrature, i)
        n_tracks_x[j] = n_tracks_x[i]
        n_tracks_y[j] = n_tracks_y[i]
        n_tracks[j] = n_tracks[i]
    end

    n_total_tracks = sum(n_tracks)

    tracks = Vector{Vector{Track}}(undef, n_azim_2)
    for i in both_dir(azimuthal_quadrature)
        tracks[i] = Vector{Track}(undef, n_tracks[i])
    end

    tracks_by_uid = Vector{Track}(undef, n_total_tracks)

    volumes = Vector{T}(undef, num_cells(model))

    return TrackGenerator(
        mesh, bcs, azimuthal_quadrature, n_tracks_x, n_tracks_y, n_tracks, n_total_tracks,
        tracks, tracks_by_uid, tiny_step, correct_volumes, volumes
    )
end

"""
    trace!(t::TrackGenerator)

Computes and fills both the azimuthal quadrature and cyclic tracks around the rectangular
domain using the provided azimuthal angles and spacing when defining the [`TrackGenerator`](@ref)
`t`.
"""
function trace!(t::TrackGenerator{T}) where {T}
    @unpack mesh, bcs, azimuthal_quadrature = t
    @unpack n_tracks_x, n_tracks_y, n_tracks = t
    @unpack tracks, tracks_by_uid = t
    @unpack bbmin, bbmax = mesh
    @unpack δs, ϕs, ωₐ = azimuthal_quadrature

    n_azim_2 = nazim2(azimuthal_quadrature)
    n_azim_4 = nazim4(azimuthal_quadrature)

    # effective azimuthal spacings, used for some computations but not stored
    δx = Vector{T}(undef, n_azim_2)
    δy = Vector{T}(undef, n_azim_2)

    Δx, Δy = width(mesh), height(mesh)

    for i in right_dir(azimuthal_quadrature)

        # effective azimuthal angle
        ϕ = ϕs[i] = atan((Δy * n_tracks_x[i]) / (Δx * n_tracks_y[i]))

        # effective azimuthal spacings
        δx[i] = Δx / n_tracks_x[i]
        δy[i] = Δy / n_tracks_y[i]
        δs[i] = δx[i] * sin(ϕ)

        # supplementaries:
        j = suplementary_idx(azimuthal_quadrature, i)
        ϕs[j] = π - ϕ
        δx[j] = δx[i]
        δy[j] = δy[i]
        δs[j] = δs[i]
    end

    # once we have computed all the azimuthal angles, compute weights
    init_weights!(azimuthal_quadrature)

    # mesh vertices and side segments
    p1 = bbmin
    p2 = Point2D(bbmin[1], bbmax[2])
    p3 = bbmax
    p4 = Point2D(bbmax[1], bbmin[2])
    sides = (top=Segment(p2, p3), bottom=Segment(p4, p1),
             right=Segment(p3, p4), left=Segment(p1, p2))

    uid = 1
    for i in both_dir(azimuthal_quadrature)

        # get azimuthal angle
        ϕ = ϕs[i]

        # for all tracks in a given azimuthal direction
        for j in 1:n_tracks[i]

            if origins_in_x(n_tracks_x, i, j)
                if points_right(azimuthal_quadrature, i)
                    p = Point2D(δx[i] * (n_tracks_x[i] - j + 1 / 2), 0)
                else
                    p = Point2D(δx[i] * (j - 1 / 2), 0)
                end
            else
                if points_right(azimuthal_quadrature, i)
                    p = Point2D(0, δy[i] * (j - n_tracks_x[i] - 1 / 2))
                else
                    p = Point2D(Δx, δy[i] * (j - n_tracks_x[i] - 1 / 2))
                end
            end

            # it can exit at y = Δy regardless if it points to the right or left
            m = tan(ϕ)
            q = Point2D(p[1] - (p[2] - Δy) / m, Δy)

            # keep going if `q` is not an exit point
            if !(0 ≤ q[1] ≤ Δx)

                # it can exit at x = Δx if it points to the right
                if points_right(azimuthal_quadrature, i)
                    q = Point2D(Δx, p[2] + m * (Δx - p[1]))

                # or it can exit at x = 0 if it points to the left
                else
                    q = Point2D(0, p[2] - m * p[1])
                end

                if !(0 ≤ q[2] ≤ Δy)
                    throw(DomainError("could not found track exit point."))
                end
            end

            # recalibrate coordinates and compute distance
            p += bbmin
            q += bbmin
            ℓ = norm(p - q)

            ABC = general_form(p, q)
            segments = Vector{Segment{T}}(undef, 0)

            BCFwd = boundary_condition(q, sides, bcs)
            BCBwd = boundary_condition(p, sides, bcs)

            # alternative method
            if points_right(azimuthal_quadrature, i)

                if j ≤ n_tracks_y[i]
                    BCFwd1 = bcs.right
                else
                    BCFwd1 = bcs.top
                end

                if j ≤ n_tracks_x[i]
                    BCBwd1 = bcs.bottom
                else
                    BCBwd1 = bcs.left
                end

            else

                if j ≤ n_tracks_y[i]
                    BCFwd1 = bcs.left
                else
                    BCFwd1 = bcs.top
                end

                if j ≤ n_tracks_x[i]
                    BCBwd1 = bcs.bottom
                else
                    BCBwd1 = bcs.right
                end

            end

            if !isequal(BCFwd, BCFwd1) || !isequal(BCBwd, BCBwd1)
                error("Boundaries do not match!")
            end

            if j ≤ n_tracks_y[i]
                DirNextTrackFwd = Forward
            else
                if BCFwd == Periodic
                    DirNextTrackFwd = Forward
                elseif BCFwd == Vaccum || BCFwd == Reflective
                    DirNextTrackFwd = Backward
                end
            end

            if j ≤ n_tracks_x[i]
                if BCBwd == Periodic
                    DirNextTrackBwd = Backward
                elseif BCBwd == Vaccum || BCBwd == Reflective
                    DirNextTrackBwd = Forward
                end
            else
                DirNextTrackBwd = Backward
            end

            track = Track{BCFwd,BCBwd,DirNextTrackFwd,DirNextTrackBwd}(
                uid, i, j, p, q, ϕ, ℓ, ABC, segments
            )
            tracks_by_uid[uid] = track
            tracks[i][j] = track
            uid += 1
        end
    end

    # once we have all tracks and boundary conditions, set next tracks
    next_tracks(t)

    return t
end

function next_tracks(t::TrackGenerator)
    @unpack tracks_by_uid = t

    # search for the next track in forward and backward direction
    for track in tracks_by_uid
        next_track_fwd(t, track)
        next_track_bwd(t, track)
    end

    return nothing
end

function next_track_fwd(t::TrackGenerator, track::Track)
    @unpack azimuthal_quadrature, tracks = t
    @unpack n_tracks_x, n_tracks_y, n_tracks = t
    @unpack azim_idx, track_idx = track

    i, j, k = azim_idx, track_idx, suplementary_idx(azimuthal_quadrature, azim_idx)

    BCFwd = bc_fwd(track)

    # these are the tracks that arrive to the y-axis
    if j ≤ n_tracks_y[i]
        if BCFwd == Periodic
            track.next_track_fwd = tracks[i][j + n_tracks_x[i]]
        elseif BCFwd == Vaccum || BCFwd == Reflective
            track.next_track_fwd = tracks[k][j + n_tracks_x[i]]
        end
    # these are the tracks that arrive to the top (superior x-axis)
    else
        if BCFwd == Periodic
            track.next_track_fwd = tracks[i][j - n_tracks_y[i]]
        elseif BCFwd == Vaccum || BCFwd == Reflective
            track.next_track_fwd = tracks[k][n_tracks[i] + n_tracks_y[i] - j + 1]
        end
    end

    return nothing
end

function next_track_bwd(t::TrackGenerator, track::Track)
    @unpack azimuthal_quadrature, tracks = t
    @unpack n_tracks_x, n_tracks_y, n_tracks = t
    @unpack azim_idx, track_idx = track

    i, j, k = azim_idx, track_idx, suplementary_idx(azimuthal_quadrature, azim_idx)

    BCBwd = bc_bwd(track)

    # these are the tracks that arrive to the bottom (inferior x-axis)
    if j ≤ n_tracks_x[i]
        if BCBwd == Periodic
            track.next_track_bwd = tracks[i][j + n_tracks_y[i]]
        elseif BCBwd == Vaccum || BCBwd == Reflective
            track.next_track_bwd = tracks[k][n_tracks_x[i] - j + 1]
        end
    # these are the tracks that arrive to the y-axis
    else
        if BCBwd == Periodic
            track.next_track_bwd = tracks[i][j - n_tracks_x[i]]
        elseif BCBwd == Vaccum || BCBwd == Reflective
            track.next_track_bwd = tracks[k][j - n_tracks_x[i]]
        end
    end

    return nothing
end

"""
    segmentize!(t::TrackGenerator)

Segmentize tracks, i.e. divide the tracks in segments generated by the intersections between
the cells or elements of the mesh. This function call is intended to be done after calling
`[trace!]`(@ref).
"""
function segmentize!(t::TrackGenerator{T}; k::Int=5, rtol::Real=Base.rtoldefault(T)) where {T}
    @unpack tracks_by_uid = t

    !isassigned(tracks_by_uid) && error("Segmentation is intended after tracing. Please, " *
                                        "call `trace!` first!")
    for track in tracks_by_uid
        _segmentize_track!(t, track, k, rtol)
    end

    fill_volumes(t, 1)

    return t
end

function fill_volumes(t::TrackGenerator{T}, i) where {T}
    @unpack tracks_by_uid, azimuthal_quadrature, volumes = t
    @unpack δs = azimuthal_quadrature
    n_azim_2 = nazim2(azimuthal_quadrature)

    fill!(volumes, zero(T))

    for track in tracks_by_uid
        a = track.azim_idx
        for segment in track.segments
            i = segment.element
            volumes[i] += δs[a] * segment.ℓ
        end
    end

    volumes ./= n_azim_2

    #### TODO: correct volumes by changing segment lengths #####
    @unpack mesh = t
    @unpack model, cell_nodes = mesh
    ncells = num_cells(model)
    volumes2 = Vector{T}(undef, ncells)
    fill!(volumes2, zero(T))
    for i in 1:ncells
        node_ids = cell_nodes[i]
        volumes2[i] = element_volume(mesh, node_ids)
    end

    return nothing
end

function element_volume(mesh, node_ids)
    @unpack model = mesh
    node_coordinates = get_node_coordinates(get_grid(model))

    x1 = node_coordinates[node_ids[1]]
    x2 = node_coordinates[node_ids[2]]
    x3 = node_coordinates[node_ids[3]]

    return 1 / 2 * abs((x2 - x1) × (x3 - x1))
end


text:
# fixed or variable sampling time 1/0
fixTime = 0		# default: 0 (variable time steps)

#### problem parameters ####
TsPF = 0.05
# nominal sampling time
sampleN = 3 		# down-sampling from Hybrid A* to OBCA
if fixTime == 1 	#
	Ts = 0.65/3*sampleN		# 0.65/3 must be compatible with motion resolution of Hybrid A* algorithm
else
	Ts = 0.6/3*sampleN 		# 0.6/3 must be compatible with motion resolution of Hybrid A* algorithm
end

# wheelbase
#L  = 2.7            # wheel base length
motionStep = 0.1	# step length of Hybrid A*",

# "nominal" shape of ego/controlled car,  ego object is later rotated around the car center
# center of rear wheel axis is reference point
# size of car is: (x_upper + x_lower) + (y_upper + y_lower)
#	   [x_upper, y_upper, -x_lower, -y_lower ]
ego  = [ 3.7   , 1      ,  1      ,  1       ]

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc


text:
# Libcurl Julia bindings

# TODOs: tests, Method abstraction for libcurl ccall's

module Curl

warn("Curl.jl has been deprecated. Please check out Requests.jl (http://github.com/loladiro/Requests.jl) for all your HTTP client needs.")

include("constants.jl")

const DEBUG = false
const version = "0.0.1"

# global
raw_response = { :headers => "", :text => "" }

type Response
  headers::Array{Array{String,1}}
  text::String
end

function escape(curl, str::String)
  p = ccall( (:curl_easy_escape, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Ptr{Uint8}, Int), curl, str.data, 0)
  bytestring(p)
end

function dict_to_query_params(curl, d::Dict)
  join([join([string(x), escape(curl, d[x])], "=") for x=keys(d)], "&")
end

function dict_to_query_params(curl, d::Array{Any,1})
  ""
end

function write_function(buf::Ptr{Uint8}, size::Uint32, nmemb::Uint32, userp::Ptr{Void})
  global resp_text
  gen_write_cb(:text, buf, size, nmemb, userp)
  nmemb
end

function header_function(buf::Ptr{Uint8}, size::Uint32, nmemb::Uint32, userp::Ptr{Void})
  gen_write_cb(:headers, buf, size, nmemb, userp)
  nmemb
end

function gen_write_cb(resp_type::Symbol, buf::Ptr{Uint8}, size::Uint32, nmemb::Uint32, userp::Ptr{Void})

  # To test, can call directly as
  # ccall(c_write_function, (Uint32), (Ptr{Uint32}, Uint32, Uint32, Ptr{Void}), 0x00, 0x00, 0x00, 0x00)

  if DEBUG; println("write_function"); end

  # convert size provided in hex to dec
  num_bytes = parseint(base(10, nmemb))
  arr = zeros(Char, num_bytes)
  for i = 1:num_bytes
    c = char(parseint(base(10, unsafe_load(buf, i))))
    arr[i] = c
  end
  str_val = join(arr)
  global raw_response
  raw_response[resp_type] = string(raw_response[resp_type], str_val)

  if DEBUG
    println("size: $size")
    println("nmemb $nmemb")
    println("num bytes: $num_bytes")
    println("length of str val: $(length(str_val))")
    println("length of resp data: $(length(raw_response[resp_type]))")
  end

  # just report back that we saved what was sent in as input
  return nmemb

end
c_write_function = cfunction(write_function, Uint32, (Ptr{Uint8}, Uint32, Uint32, Ptr{Void}))
c_header_function = cfunction(header_function, Uint32, (Ptr{Uint8}, Uint32, Uint32, Ptr{Void}))

function curl_version()
  curl_ver = bytestring(ccall( (:curl_version, "libcurl"), Ptr{Uint8}, ()))
  curl_ver
end

function setup_curl()

  curl_ver = curl_version()
  if DEBUG; println("curl version: $curl_ver"); end
  curl = ccall( (:curl_easy_init, "libcurl"), Ptr{Uint8}, ())

end

function setup_curlopts(curl, url)
  iostr = IOString()

  ccall( (:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Int), curl, CURLOPT_FOLLOWLOCATION, 1)

  # instruct libcurl to not include the headers in the start of the body output
  ccall( (:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Int), curl, CURLOPT_HEADER, 0)
  ccall( (:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Int), curl, CURLOPT_HEADERFUNCTION, c_header_function)

  ccall( (:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_URL, url.data)
  ccall( (:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Int), curl, CURLOPT_WRITEFUNCTION, c_write_function)
  ccall( (:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_FILE, iostr.data)
end

function parse_headers(raw_headers)
  raw_header_groups = filter((x)->length(x) > 0, split(raw_headers, "\r\n\r\n"))
  [split(group, "\r\n") for group in raw_header_groups]
end

function do_curl(curl)

  global raw_response = { :headers => "", :text => "" }
  curl_resp = ccall( (:curl_easy_perform, "libcurl"), Ptr{Uint8}, (Ptr{Uint8},), curl)

  global raw_response
  # response values are now ready
  if DEBUG; println("result: $curl_resp"); end

  headers = parse_headers(raw_response[:headers])
  response = Response(headers, raw_response[:text])

  err = ccall( (:curl_easy_strerror, "libcurl"), Ptr{Uint8}, (Ptr{Uint8},), curl_resp)
  if DEBUG; println("error: $(bytestring(err))"); end

  response
end

function cleanup_curl(curl)
  ccall( (:curl_easy_cleanup, "libcurl"), Ptr, (Ptr{Uint8},), curl)
end
# finalizer(curl, cleanup_curl)

function set_ua(curl)
  ua_str = "Curl.jl Ver. $version - $(curl_version())"
  ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_USERAGENT, ua_str.data)
end

macro run_with_block(expr)
  quote
    curl = setup_curl()
    setup_curlopts(curl, url)
    set_ua(curl)

    $expr

    response = do_curl(curl)
    response

  end
end

# Http methods
function head(url)
  response = @run_with_block begin
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Int), curl, CURLOPT_NOBODY, 1)
  end
  response
end

function get(url)
  @run_with_block begin end
end

function post(url, params)
  @run_with_block begin
    escaped_params = dict_to_query_params(curl, params)
    # curl_easy_setopt(curl, CURLOPT_POSTFIELDS, "name=daniel&project=curl");
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_POSTFIELDS, escaped_params.data)
  end
end

function post(url, params)
  @run_with_block begin
    escaped_params = dict_to_query_params(curl, params)
    # curl_easy_setopt(curl, CURLOPT_POSTFIELDS, "name=daniel&project=curl");
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_POSTFIELDS, escaped_params.data)
  end
end

function patch(url, params)
  @run_with_block begin
    escaped_params = dict_to_query_params(curl, params)
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_POSTFIELDS, escaped_params.data)
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_CUSTOMREQUEST, "PATCH".data)
  end
end

function put(url, params)
  @run_with_block begin
    escaped_params = dict_to_query_params(curl, params)
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_POSTFIELDS, escaped_params.data)
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_CUSTOMREQUEST, "PUT".data)
  end
end

function delete(url)
  @run_with_block begin
    # curl_easy_setopt(handle, CURLOPT_CUSTOMREQUEST, "DELETE");
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_CUSTOMREQUEST, "DELETE".data)
  end
end

function options(url)
  @run_with_block begin
    # curl_easy_setopt(handle, CURLOPT_CUSTOMREQUEST, "OPTIONS");
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Int), curl, CURLOPT_NOBODY, 1)
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_CUSTOMREQUEST, "OPTIONS".data)
  end
end

function trace(url)
  @run_with_block begin
    # curl_easy_setopt(handle, CURLOPT_CUSTOMREQUEST, "OPTIONS");
    ccall((:curl_easy_setopt, "libcurl"), Ptr{Uint8}, (Ptr{Uint8}, Int, Ptr{Uint8}), curl, CURLOPT_CUSTOMREQUEST, "TRACE".data)
  end
end

end


text:
#!/usr/bin/env julia

##########################################################################
# run_SNR.jl
#
# Run all SNR simulations.
##########################################################################
include("run_SNR-ibc.jl")
include("run_SNR-indoors.jl")
include("run_SNR-randomlargescale.jl")
include("run_SNR-triangularmacro.jl")
include("run_SNR-triangularhetnet.jl")


text:
# ---
# title: 979. Distribute Coins in Binary Tree
# id: problem979
# author: Indigo
# date: 2021-02-21
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/distribute-coins-in-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree with `N` nodes, each `node` in the tree has
# `node.val` coins, and there are `N` coins total.
#
# In one move, we may choose two adjacent nodes and move one coin from one node
# to another.  (The move may be from parent to child, or from child to parent.)
#
# Return the number of moves required to make every node have exactly one coin.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree1.png)**
#
#
#
#     Input: [3,0,0]
#     Output: 2
#     Explanation: From the root of the tree, we move one coin to its left child, and one coin to its right child.
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree2.png)**
#
#
#
#     Input: [0,3,0]
#     Output: 3
#     Explanation: From the left child of the root, we move two coins to the root [taking two moves].  Then, we move one coin from the root of the tree to the right child.
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree3.png)**
#
#
#
#     Input: [1,0,2]
#     Output: 2
#
#
# **Example 4:**
#
# **![](https://assets.leetcode.com/uploads/2019/01/18/tree4.png)**
#
#
#
#     Input: [1,0,0,null,3]
#     Output: 4
#
#
#
#
# **Note:**
#
#   1. `1<= N <= 100`
#   2. `0 <= node.val <= N`
#
#
## @lc code=start
using LeetCode

function distribute_coins(root::TreeNode{Int})
    dfs(::Nothing) = 0
    function dfs(root::TreeNode{Int})
        l = dfs(root.left)
        r = dfs(root.right)
        res = res + abs(l) + abs(r)
        root.val + l + r - 1
    end
    res = 0
    dfs(root)
    res
end
## @lc code=end


text:
export MapFoldl, umm!
export inorm2, idot

struct MapFoldl{FM, FF} <: Function
    m::FM
    f::FF
    function MapFoldl(m::FM, f::FF) where {FM, FF}
        !isreversible(f) && throw(InvertibilityError("second parameter $f is not reversible!"))
        new{FM, FF}(m, f)
    end
end

@i function (mf::MapFoldl)(out::T, iter) where T
    @anc anc = zero(T)
    for i=1:1:length(iter)
        anc += mf.m(iter[i])
        mf.f(out, anc)
        anc -= mf.m(iter[i])
    end
end

# NOTE: also define the multiplication between two matrices
@i function umm!(x, θ, Nin::Int, Nout::Int)
    @anc k = 0
    for j=1:Nout
        for i=Nin-1:-1:j
            k += identity(1)
            ROT(x[i], x[i+1], θ[k])
        end
    end

    # uncompute k
    for j=1:Nout
        for i=Nin-1:-1:j
            k -= identity(1)
        end
    end
end

@i function idot(out, v1::AbstractVector{T}, v2) where T
    @anc anc1 = zero(T)
    for i = 1:length(v1)
        anc1 += identity(v1[i])
        CONJ(anc1)
        out += v1[i]*v2[i]
        CONJ(anc1)
        anc1 -= identity(v1[i])
    end
end

@i function inorm2(out, vec::AbstractVector{T}) where T
    @anc anc1 = zero(T)
    for i = 1:length(vec)
        anc1 += identity(vec[i])
        CONJ(anc1)
        out += anc1*vec[i]
        CONJ(anc1)
        anc1 -= identity(vec[i])
    end
end