text:
# Set the IPP library path
@static if Sys.isapple()
    ipplib = ENV["IPPROOT"] * "/lib"
end

@static if Sys.islinux()
    ipplib = ENV["IPPROOT"] * "/lib/intel64"
end

@static if Sys.iswindows()
    error("Unsupported OS: Windows")
end

const ippC1 = 1
const ippBorderRepl = 1

const ippNearest = 1
const ippCubic = 6
const ippSuper = 8
const ippLanczos = 16

mutable struct IppiSize
    width::Cint
    height::Cint
end

mutable struct IppiPoint
    x::Cint
    y::Cint
end

mutable struct IppiBorderSize
    borderLeft::Cint
    borderTop::Cint
    borderRight::Cint
    borderBottom::Cint
end

function resizeCubic32fC1R(src::Matrix{Float32}, width::Integer, height::Integer)
    # a destination buffer
    dst = Matrix{Float32}(undef, width, height)

    dims = size(src)
    srcSize = IppiSize(dims[1], dims[2])
    dstSize = IppiSize(width, height)

    specSize = Ref{Cint}(0)
    initSize = Ref{Cint}(0)

    # Spec and init buffer sizes
    status = ccall(
        (:ippiResizeGetSize_32f, ipplib * "/libippi.so"),
        Cint,
        (IppiSize, IppiSize, Cint, Cint, Ref{Cint}, Ref{Cint}),
        srcSize,
        dstSize,
        ippCubic,
        0,
        specSize,
        initSize,
    )

    if status != 0
        error("ippiResizeGetSize_32f::$status")
    end

    # Memory allocation
    pInitBuf =
        ccall((:ippsMalloc_8u, ipplib * "/libipps.so"), Ptr{Cvoid}, (Cint,), initSize[])

    pSpec = ccall((:ippsMalloc_8u, ipplib * "/libipps.so"), Ptr{Cvoid}, (Cint,), specSize[])

    if pInitBuf == C_NULL || pSpec == C_NULL
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pInitBuf)
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)
        error("NULL pInitBuf or pSpec")
    end

    valueB = 0.0f0
    valueC = 0.0f0

    # Filter initialization
    status = ccall(
        (:ippiResizeCubicInit_32f, ipplib * "/libippi.so"),
        Cint,
        (IppiSize, IppiSize, Cfloat, Cfloat, Ptr{Cvoid}, Ptr{Cvoid}),
        srcSize,
        dstSize,
        valueB,
        valueC,
        pSpec,
        pInitBuf,
    )
    ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pInitBuf)

    if status != 0
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)
        error("ippiResizeCubicInit_32f::$status")
    end

    borderSize = Ref{IppiBorderSize}(IppiBorderSize(0, 0, 0, 0))

    status = ccall(
        (:ippiResizeGetBorderSize_32f, ipplib * "/libippi.so"),
        Cint,
        (Ptr{Cvoid}, Ref{IppiBorderSize}),
        pSpec,
        borderSize,
    )

    if status != 0
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)
        error("ippiResizeGetBorderSize_32f::$status")
    end

    bufSize = Ref{Cint}(0)

    # General transform function
    status = ccall(
        (:ippiResizeGetBufferSize_32f, ipplib * "/libippi.so"),
        Cint,
        (Ptr{Cvoid}, IppiSize, Cint, Ref{Cint}),
        pSpec,
        dstSize,
        ippC1,
        bufSize,
    )

    if status != 0
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)
        error("ippiResizeGetBufferSize_32f::$status")
    end

    pBuffer =
        ccall((:ippsMalloc_8u, ipplib * "/libipps.so"), Ptr{Cvoid}, (Cint,), bufSize[])

    if pBuffer != C_NULL
        srcOffset = IppiPoint(0, 0)
        dstOffset = IppiPoint(0, 0)

        status = ccall(
            (:ippiResizeGetSrcRoi_32f, ipplib * "/libippi.so"),
            Cint,
            (Ptr{Cvoid}, IppiPoint, IppiSize, Ref{IppiPoint}, Ref{IppiSize}),
            pSpec,
            dstOffset,
            dstSize,
            srcOffset,
            srcSize,
        )

        if status == 0
            # finally resize the image
            srcStep = strides(src)[2] * sizeof(Float32)
            dstStep = strides(dst)[2] * sizeof(Float32)

            status = ccall(
                (:ippiResizeCubic_32f_C1R, ipplib * "/libippi.so"),
                Cint,
                (
                    Ptr{Cfloat},
                    Cint,
                    Ptr{Cfloat},
                    Cint,
                    IppiPoint,
                    IppiSize,
                    Cint,
                    Ptr{Cfloat},
                    Ptr{Cvoid},
                    Ptr{Cvoid},
                ),
                src,
                srcStep,
                dst,
                dstStep,
                dstOffset,
                dstSize,
                ippBorderRepl,
                C_NULL,
                pSpec,
                pBuffer,
            )

        end

        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pBuffer)
    end

    # Release memory
    ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)

    return dst
end

function resizeNearest8uC1R(src::Matrix{UInt8}, width::Integer, height::Integer)
    # a destination buffer
    dst = Matrix{UInt8}(undef, width, height)

    dims = size(src)
    srcSize = IppiSize(dims[1], dims[2])
    dstSize = IppiSize(width, height)

    specSize = Ref{Cint}(0)
    initSize = Ref{Cint}(0)

    # Spec and init buffer sizes
    status = ccall(
        (:ippiResizeGetSize_8u, ipplib * "/libippi.so"),
        Cint,
        (IppiSize, IppiSize, Cint, Cint, Ref{Cint}, Ref{Cint}),
        srcSize,
        dstSize,
        ippNearest,
        0,
        specSize,
        initSize,
    )

    if status != 0
        error("ippiResizeGetSize_8u::$status")
    end

    # Memory allocation
    pSpec = ccall((:ippsMalloc_8u, ipplib * "/libipps.so"), Ptr{Cvoid}, (Cint,), specSize[])

    if pSpec == C_NULL
        error("NULL pSpec")
    end

    # Filter initialization
    status = ccall(
        (:ippiResizeNearestInit_8u, ipplib * "/libippi.so"),
        Cint,
        (IppiSize, IppiSize, Ptr{Cvoid}),
        srcSize,
        dstSize,
        pSpec,
    )

    if status != 0
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)
        error("ippiResizeNearestInit_8u::$status")
    end

    bufSize = Ref{Cint}(0)

    # General transform function
    status = ccall(
        (:ippiResizeGetBufferSize_8u, ipplib * "/libippi.so"),
        Cint,
        (Ptr{Cvoid}, IppiSize, Cint, Ref{Cint}),
        pSpec,
        dstSize,
        ippC1,
        bufSize,
    )

    if status != 0
        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)
        error("ippiResizeGetBufferSize_8u::$status")
    end

    pBuffer =
        ccall((:ippsMalloc_8u, ipplib * "/libipps.so"), Ptr{Cvoid}, (Cint,), bufSize[])

    if pBuffer != C_NULL
        srcOffset = IppiPoint(0, 0)
        dstOffset = IppiPoint(0, 0)

        status = ccall(
            (:ippiResizeGetSrcRoi_8u, ipplib * "/libippi.so"),
            Cint,
            (Ptr{Cvoid}, IppiPoint, IppiSize, Ref{IppiPoint}, Ref{IppiSize}),
            pSpec,
            dstOffset,
            dstSize,
            srcOffset,
            srcSize,
        )

        if status == 0
            # finally resize the image
            srcStep = strides(src)[2] * sizeof(UInt8)
            dstStep = strides(dst)[2] * sizeof(UInt8)

            status = ccall(
                (:ippiResizeNearest_8u_C1R, ipplib * "/libippi.so"),
                Cint,
                (Ptr{UInt8}, Cint, Ptr{UInt8}, Cint, IppiPoint, IppiSize, Cint, Ptr{Cvoid}),
                src,
                srcStep,
                dst,
                dstStep,
                dstOffset,
                dstSize,
                pSpec,
                pBuffer,
            )

        end

        ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pBuffer)
    end

    # Release memory
    ccall((:ippsFree, ipplib * "/libipps.so"), Cvoid, (Ptr{Cvoid},), pSpec)

    return dst
end


text:
fname = :ARCSpectral_abs

c = ARCTR.Combi(ARCTR.hessian_dense,ARCTR.PDataSpectral,ARCTR.solve_modelARCDiagAbs,ARCTR.preprocessSpectral,ARCTR.decreaseFact,ARCTR.Tparam())

include("Template.jl")

ext: jl


text:
### Combination methods for p-values ###

"""
    combine(PValues, <:PValueCombination)

Combine p-values

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, Fisher())
0.007616871850449092

julia> combine(pvals, Stouffer())
0.00709832618126593
```

# See also

`PValueCombination`s:

[`Fisher`](@ref)
[`Logit`](@ref)
[`Stouffer`](@ref)
[`Tippett`](@ref)
[`Simes`](@ref)
[`Wilkinson`](@ref)
[`Minimum`](@ref)

"""
function combine end

function combine(pValues::AbstractVector{T}, method::M)::T where {T <: AbstractFloat,M <: PValueCombination}
    return combine(PValues(pValues), method)
end

## Fisher combination ##

"""
Fisher's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, Fisher())
0.007616871850449092
```

# References

Fisher, R.A. (1925).
Statistical methods for research workers
(Genesis Publishing Pvt Ltd).

"""
struct Fisher <: PValueCombination
end

function combine(pValues::PValues{T}, method::Fisher)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0
        return NaN
    end
    x = -2 * sum(log.(pValues))
    p = ccdf(Chisq(2n), Float64(x))
    return p
end

## Logit combination ##

"""
Logit p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, Logit())
0.006434494635148462
```

# References

Mudholkar, G.S., and George, E.O. (1977).
The Logit Statistic for Combining Probabilities - An Overview
(Rochester University NY, Dept of Statistics).

"""
struct Logit <: PValueCombination
end

function combine(pValues::PValues{T}, method::Logit)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0 || maximum(pValues) == 1
        return NaN
    end
    c = sqrt((5n + 2) * n * pi^2 / ((5n + 4) * 3))
    x = -sum(log.(pValues ./ (1 .- pValues))) / c
    p = ccdf(TDist(5n + 4), x)
    return p
end

## Stouffer combination ##

"""
Stouffer's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, Stouffer())
0.00709832618126593

julia> weights = [1.0, 2.0, 0.4, 1.5];

julia> combine(pvals, weights, Stouffer())
0.007331653763696742
```

# References

Stouffer, S.A. (1949).
The American soldier. Vol. 1: Adjustment during army life
(Princeton University Press).

Liptak, T. (1958).
On the combination of independent tests.
Magyar Tud Akad Mat Kutato Int Kozl 3, 171–197.

"""
struct Stouffer <: PValueCombination
end

function combine(pValues::PValues{T}, method::Stouffer)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0 || maximum(pValues) == 1
        return NaN
    end
    z = cquantile.(Ref(Normal()), pValues)
    z = sum(z) ./ sqrt(n)
    p = ccdf(Normal(), z)
    return p
end

function combine(pValues::PValues{T}, weights::AbstractVector{R}, method::Stouffer)::T where {T <: AbstractFloat,R <: Real}
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0 || maximum(pValues) == 1
        return NaN
    end
    z = cquantile.(Ref(Normal()), pValues) .* weights
    z = sum(z) ./ sqrt(sum(abs2, weights))
    p = ccdf(Normal(), z)
    return p
end

function combine(pValues::AbstractVector{T}, weights::AbstractVector{R}, method::Stouffer)::T where {T <: AbstractFloat,R <: Real}
    return combine(PValues(pValues), weights, method)
end

## Tippett combination ##

"""
Tippett's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, Tippett())
0.039403990000000055
```

# References

Tippett, L.H.C. (1931). The Methods of Statistics. An introduction mainly for
workers in the biological sciences.

"""
struct Tippett <: PValueCombination
end

function combine(pValues::PValues{T}, method::Tippett)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    p = 1 - (1 - minimum(pValues))^n
    return p
end

## Simes combination ##

"""
Simes's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, Simes())
0.04
```

# References

Simes, R.J. (1986). An improved Bonferroni procedure for multiple tests of
significance. Biometrika 73, 751–754.

"""
struct Simes <: PValueCombination
end

function combine(pValues::PValues{T}, method::Simes)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    pValues = sort(pValues)
    p = n * minimum(pValues ./ (1:n))
    return p
end

## Wilkinson combination ##

"""
Wilkinson's p-value combination

# Examples

```jldoctest
julia> pv = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pv, Wilkinson(1))  # combination with rank 1
0.03940399000000003

julia> combine(pv, Wilkinson(4))  # combination with rank 4
0.0625
```

# References

Wilkinson, B. (1951). A statistical consideration in psychological research.
Psychological Bulletin 48, 156.

"""
struct Wilkinson <: PValueCombination
    rank::Integer

    function Wilkinson(rank)
        if rank < 1
            throw(ArgumentError("Rank must be positive."))
        end
        return new(rank)
    end
end

function combine(pValues::PValues{T}, method::Wilkinson)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    rank = method.rank
    if rank < 1 || rank > n
        throw(ArgumentError("Rank must be in 1,..,$(n)"))
    end
    p_rank = partialsort(pValues, rank)
    p = cdf(Beta(rank, n - rank + 1), Float64(p_rank))
    return p
end

## Generalised minimum combination ##

"""
Minimum of adjusted p-value combination

# Examples

```jldoctest
julia> pv = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pv, Minimum(BenjaminiHochberg()))
0.04

julia> combine(pv, Minimum(ForwardStop()))
0.01005033585350145
```

"""
struct Minimum <: PValueCombination
    adjustment::PValueAdjustment
end

function combine(pValues::PValues{T}, method::Minimum)::T where T <: AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    padj = adjust(pValues, method.adjustment)
    p = minimum(padj)
    return p
end


text:
module DeeplyNested

  using Boot

  include_folder(DeeplyNested, @__FILE__)

end


text:
using QuantLib

function main()
  euribor3m = euribor_index(QuantLib.Time.TenorPeriod(Dates.Month(3)))

  todaysDate = Date(2006, 5, 23)
  set_eval_date!(settings, todaysDate)
  calendar = euribor3m.fixingCalendar
  fixingDays = euribor3m.fixingDays
  settlementDate = QuantLib.Time.advance(Dates.Day(fixingDays), calendar, todaysDate)

  threeMonthFraQuote = zeros(10)
  threeMonthFraQuote[2] = 0.030
  threeMonthFraQuote[3] = 0.031
  threeMonthFraQuote[4] = 0.032
  threeMonthFraQuote[7] = 0.033
  threeMonthFraQuote[10] = 0.034

  ## QUOTES ##
  h1x4 = Quote(threeMonthFraQuote[2])
  h2x5 = Quote(threeMonthFraQuote[3])
  h3x6 = Quote(threeMonthFraQuote[4])
  h6x9 = Quote(threeMonthFraQuote[7])
  h9x12 = Quote(threeMonthFraQuote[10])

  ## RATE HELPERS ##
  fraDayCount = euribor3m.dc
  convention = euribor3m.convention
  endOfMonth = euribor3m.endOfMonth

  fra1x4 = FraRateHelper(h1x4, 1, 4, fixingDays, calendar, convention, endOfMonth, fraDayCount)
  fra2x5 = FraRateHelper(h2x5, 2, 5, fixingDays, calendar, convention, endOfMonth, fraDayCount)
  fra3x6 = FraRateHelper(h3x6, 3, 6, fixingDays, calendar, convention, endOfMonth, fraDayCount)
  fra6x9 = FraRateHelper(h6x9, 6, 9, fixingDays, calendar, convention, endOfMonth, fraDayCount)
  fra9x12 = FraRateHelper(h9x12, 9, 12, fixingDays, calendar, convention, endOfMonth, fraDayCount)

  fraInstruments = FraRateHelper[fra1x4, fra2x5, fra3x6, fra6x9, fra9x12]
  termStructureDC = QuantLib.Time.ISDAActualActual()
  tol = 1.0e-15

  fraTermStructure = PiecewiseYieldCurve(settlementDate, fraInstruments, termStructureDC, QuantLib.Math.LogLinear(), Discount(), tol, IterativeBootstrap())

  euribor3m = update_termstructure(euribor3m, fraTermStructure)

  fraCalendar = euribor3m.fixingCalendar
  fraBusinessDayConvention = euribor3m.convention
  fraFwdType = LongPosition()
  fraNotional = 100.0
  monthsToStart = [1, 2, 3, 6, 9]
  fraTermMonths = 3
  fraValueDate = QuantLib.Time.advance(Dates.Month(monthsToStart[1]), fraCalendar, settlementDate, fraBusinessDayConvention)
  fraMaturityDate = QuantLib.Time.advance(Dates.Month(fraTermMonths), fraCalendar, fraValueDate, fraBusinessDayConvention)
  fraStrikeRate = threeMonthFraQuote[monthsToStart[1] + 1]
  myFRA = ForwardRateAgreement(fraValueDate, fraMaturityDate, fraFwdType, fraStrikeRate, fraNotional, euribor3m, fraTermStructure)

  println(forward_rate(myFRA).rate)
  println(spot_value(myFRA))
  println(forward_value(myFRA))
  println(implied_yield(myFRA, spot_value(myFRA), forward_value(myFRA), settlementDate, SimpleCompounding(), fraDayCount).rate)
end


text:
# functions to build MNA from parsed circuit
#

function system_matrix_size(x::MNAbuilder)
    mx(a) = length(a) == 0 ? 0 : maximum(a) # allow for size 0
    maximum((mx(x.Gi), mx(x.Gj),
            mx(x.Hi),
            mx(x.Di), mx(x.Dj),
            mx(x.Si),
            mx(keys(x.s)) )
            )
end

function mna(pc::ParsedCircuit, group2 = Set{Symbol}())
    x = MNAbuilder()
    for (key,value) in pc.nodedict
        x.nodedict[key] = value
    end
    for i in eachindex(pc.netlist)
        element = pc.netlist[i]
        forcegroup2 =  element.name ∈ group2
        process!(x, element, forcegroup2, pc)
    end
    MNA(x,pc.max_node)
end

# retain currents for group 2

process!(x::MNAbuilder{T}, e::Component, forcegroup2::Bool, pc) where T = println("unknown element $(e.name)")
function process!(x::MNAbuilder{T}, e::Resistor, forcegroup2::Bool, pc) where T
    if forcegroup2
        pc.max_node += 1
        x.currentdict[e.name] = pc.max_node
        pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node1,     pc.max_node,    one(T))
        pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node2,     pc.max_node,    -one(T))
        pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, e.node1,        one(T))
        pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, e.node1,        -one(T))
        pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, pc.max_node,    e.value)
    else
        conductance = 1/e.value
        pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node1, e.node1, conductance)
        pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node2, e.node2, conductance)
        pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node1, e.node2, -conductance)
        pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node2, e.node1, -conductance)
    end
end
function process!(x::MNAbuilder{T}, e::Inductor, forcegroup2::Bool, pc) where T
end
function process!(x::MNAbuilder{T}, e::Capacitor, forcegroup2::Bool, pc) where T
end
function process!(x::MNAbuilder{T}, e::VoltageSource, ::Bool, pc)  where T # group 2 only
    pc.max_node += 1
    x.currentdict[e.name] = pc.max_node
    pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, e.node1,    one(T))
    pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, e.node2,    -one(T))
    pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node1,     pc.max_node, one(T))
    pushvalue!(x.Gi,x.Gj,x.Gvalue, e.node2,     pc.max_node, -one(T))
    pushvalue!(x.Si,x.Sj,x.Svalue, pc.max_node, 1,           e.value)
end
function process!(x::MNAbuilder{T}, e::CurrentSource, forcegroup2::Bool, pc) where T
    if forcegroup2
        pc.max_node += 1
        x.currentdict[e.name] = pc.max_node
        pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, e.node1     , one(T))
        pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, e.node2     , -one(T))
        pushvalue!(x.Gi,x.Gj,x.Gvalue, pc.max_node, pc.max_node , one(T))
        pushvalue!(x.Si,x.Sj,x.Svalue, ps.max_node, 1           , e.value)
    else
        pushvalue!(x.Si,x.Sj,x.Svalue, e.node1, 1, -e.value)
        pushvalue!(x.Si,x.Sj,x.Svalue, e.node2, 1, e.value)
    end
end

function pushvalue!(i_array, j_array, value_array, i, j, value)
    if i>0 && j>0
        push!(i_array, i)
        push!(j_array, j)
        push!(value_array, value)
    end
end


text:
using Tabben: TabularDataset

function random_classification_dataset(n, m, k, c)
    input = rand(m, n)
    output = rand(0:c-1, k, n)
    metadata = Dict(
        "data_location" => "https://url/to/the/npz/file.npz",
        "extras_location" => "https://url/to/the/extras/file.json",
        "task" => "classification",
        "classes" => c,
    )
    extras = Dict(
        "license" => "test license",
        "bibtex" => "test bibtex",
        "column-names-attributes" => ["A$i" for i in 1:m],
        "column-names-target" => ["T$i" for i in 1:k]
    )

    return TabularDataset(input, output, metadata, extras)
end


text:
# we are making a SBML model if this gets called, so lets load the SBML strategy -
include("./strategy/SBMLStrategy.jl")

function make_sbml_model(path_to_model_file::String, path_to_output_dir::String;
    defaults_file_name::String="Defaults.toml", default_model_file_name::String="Model.sbml")

    # start -
    time_start = time();

    # Fix: Need to generate a Default.toml project first #2
    # check, do we have the Defaults file?
    path_to_defaults_file = "$(dirname(path_to_model_file))/$(defaults_file_name)"
    if (isfile(path_to_defaults_file) == false)
        throw(ArgumentError("Ooops! We are missing the Defaults.toml file. Please generate a default project using the generate_default_project function. For details see the generate_default_project help in the REPL"))
    end

    @info "Loading the $(defaults_file_name)"

    # Load the defaults file -
    tmp_dict =  TOML.parsefile(path_to_defaults_file)
    default_parameter_dictionary = tmp_dict["default_parameter_values"]
    author_metadata_dictionary = tmp_dict["author_metadata"]

    # legit path to model file?
    is_file_path_ok(path_to_model_file)

    # before we get too far along, we need to check if the user already has code in the location that we want to generate at -
    # if they do, then move it -
    if (isdir(path_to_output_dir) == true)

        @info "Ooops! We found some code where you wanted to generate your model. Moving ..."

        # ok, looks like we may have a conflict - mv the offending code
        if (move_existing_project_at_path(path_to_output_dir) == false)

            # Something happend ... the world is ending ...
            throw(ArgumentError("automatic directory conflict resolution failed. Unable to move existing directory $(path_to_output_dir)"))
        end
    end

    # let the uset know what is going on -
    @info "Starting a new code genration job: reading model file at $(path_to_model_file)"

    # initialize -
    src_component_set = Set{VLProgramComponent}()
    config_component_set = Set{VLProgramComponent}()
    root_component_set = Set{VLProgramComponent}()

    # Load the JSON configuration file -
    path_to_configuration_file = "$(path_to_package)/configuration/Configuration.json"
    config_dict = JSON.parsefile(path_to_configuration_file)

    # let the uset know what is going on -
    @info "Loaded the configuration file $(path_to_configuration_file)"

    # parse the model file -
    problem_dictionary = parse_vff_model_file(path_to_model_file, config_dict, default_parameter_dictionary)
    problem_dictionary["configuration_dictionary"] = config_dict
    problem_dictionary["default_parameter_dictionary"] = default_parameter_dictionary
    problem_dictionary["author_metadata_dictionary"] = author_metadata_dictionary

    # let the uset know what is going on -
    @info "Model file parsing complete"

    # build the SMBL program component -
    sbml_file_program_component = build_sbml_model_program_component(problem_dictionary, default_model_file_name)
    push!(root_component_set, sbml_file_program_component)

    # stop timer -
    elapsed_time = (time() - time_start)    # in ns
    problem_dictionary["elapsed_time"] = elapsed_time
    problem_dictionary["path_to_model_file"] = path_to_model_file
    problem_dictionary["path_to_output_dir"] = path_to_output_dir

    # generate model report file (useful for debugging) -
    program_component_report = build_model_generation_report(problem_dictionary)
    push!(root_component_set, program_component_report)

    # Dump the component_set to disk -
    write_program_components_to_disk("$(path_to_output_dir)", root_component_set)
    @info "Done! elapsed_time: $(time() - time_start)s"
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
using BandedMatrices, LinearAlgebra, ArrayLayouts, Test
import Base: BroadcastStyle
import Base.Broadcast: broadcasted
import BandedMatrices: BandedStyle, BandedRows

@testset "broadcasting" begin
    @testset "general" begin
        n = 1000
        A = brand(n,n,1,1)
        B = Matrix{Float64}(undef, n,n)
        B .= exp.(A)
        @test B == exp.(Matrix(A)) == exp.(A)

        @test exp.(A) isa BandedMatrix
        @test bandwidths(exp.(A)) == (n-1,n-1)

        C = similar(A)
        @test_throws BandError C .= exp.(A)

        @test identity.(A) == A
        @test identity.(A) isa BandedMatrix
        @test bandwidths(identity.(A)) == bandwidths(A)

        @test (z -> exp(z)-1).(A) == (z -> exp(z)-1).(Matrix(A))
        @test (z -> exp(z)-1).(A) isa BandedMatrix
        @test bandwidths((z -> exp(z)-1).(A)) == bandwidths(A)

        @test A .+ 1 == Matrix(A) .+ 1
        @test A .+ 1 isa BandedMatrix
        @test (A .+ 1) == Matrix(A) .+ 1

        @test A ./ 1 == Matrix(A) ./ 1
        @test A ./ 1 isa BandedMatrix
        @test bandwidths(A ./ 1) == bandwidths(A)

        @test 1 .+ A == 1 .+ Matrix(A)
        @test 1 .+ A isa BandedMatrix
        @test (1 .+ A) == 1 .+ Matrix(A)

        @test 1 .\ A == 1 .\ Matrix(A)
        @test 1 .\ A isa BandedMatrix
        @test bandwidths(1 .\ A) == bandwidths(A)

        A = brand(10,1,1,1)
        @test A[:,1] isa Vector
        @test norm(A .- A[:,1]) == 0
        @test A ≈ A[:,1]
    end

    @testset "identity" begin
        n = 100
        A = brand(n,n,2,2)
        A.data[1,1] = NaN
        B = brand(n,n,2,2)
        B.data[1,1] = NaN
        A .= B
        @test A == B
        B = brand(n,n,1,1)
        B.data[1,1] = NaN
        A .= B
        B = brand(n,n,3,3)
        B[band(3)] .= 0
        B[band(-3)] .= 0
        B.data[1,1] = NaN
        B.data[end,end] = NaN
        A .= B
        @test A == B

        B = brand(n,n,0,3)
        B[band(3)] .= 0
        B.data[1,1] = NaN
        A .= B
        @test A == B

        B = brand(n,n,3,0)
        B[band(-3)] .= 0
        B.data[end,end] = NaN
        A .= B
        @test A == B

        B = brand(n,n,-1,1)
        A .= B
        @test A == B

        B = brand(n,n,1,-1)
        A .= B
        @test A == B
    end

    @testset "lmul!/rmul!" begin
        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= (-).(A)
        @test -A isa BandedMatrix
        @test (-).(A) isa BandedMatrix
        @test bandwidths(A) == bandwidths(-A) == bandwidths((-).(A))
        @test B == -A == (-).(A)
        @test A-I isa BandedMatrix
        @test I-A isa BandedMatrix
        @test bandwidths(A) == bandwidths(A-I) == bandwidths(I-A)

        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= 2.0.*A

        @test B ==  2A == 2.0.*A
        @test 2A isa BandedMatrix
        @test 2.0.*A isa BandedMatrix
        @test bandwidths(2A) == bandwidths(2.0.*A) == bandwidths(A)

        A .= 2.0.*A
        @test A == B

        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= A.*2.0

        @test B ==  A*2 == A.*2.0
        @test A*2 isa BandedMatrix
        @test A .* 2.0 isa BandedMatrix
        @test bandwidths(A*2) == bandwidths(A.*2.0) == bandwidths(A)
        A .= A.*2.0
        @test A == B

        n = 1000
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= A ./ 2.0

        @test B == A/2 == A ./ 2.0
        @test A/2 isa BandedMatrix
        @test A ./ 2.0 isa BandedMatrix
        @test bandwidths(A/2) == bandwidths(A ./ 2.0) == bandwidths(A)

        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B .= 2.0 .\ A

        @test B == A/2 == A ./ 2.0
        @test 2\A isa BandedMatrix
        @test 2.0 .\ A isa BandedMatrix
        @test bandwidths(2\A) == bandwidths(2.0 .\ A) == bandwidths(A)

        A = brand(5,5,1,1)
        A.data .= NaN
        lmul!(0.0,A)
        @test isnan(norm(A)) == isnan(norm(lmul!(0.0,[NaN])))

        lmul!(false,A)
        @test norm(A) == 0.0

        A = brand(5,5,1,1)
        A.data .= NaN
        rmul!(A,0.0)
        @test isnan(norm(A)) == isnan(norm(rmul!([NaN],0.0)))

        rmul!(A,false)
        @test norm(A) == 0.0

        n = 100
        A = brand(n,n,2,2)
        B = brand(n,n,1,1)
        A[band(2)] .= A[band(-2)] .= 0
        B .= A ./ 2.0
        @test B == A / 2.0 == Matrix(A)/2.0

        B .= 2.0 .\ A
        @test B == 2.0 \ A == 2.0 \ Matrix(A)

        n = 100
        A = brand(n,n,2,2)
        B = brand(n,n,1,3)
        A[band(-2)] .= 0
        B .= A ./ 2.0
        @test B == A / 2.0 == Matrix(A)/2.0

        B .= 2.0 .\ A
        @test B == 2.0 \ A == 2.0 \ Matrix(A)

        A = brand(n,n,2,2)
        B = brand(n,n,3,1)
        A[band(2)] .= 0
        B .= A ./ 2.0
        @test B == A / 2.0 == Matrix(A)/2.0

        B .= 2.0 .\ A
        @test B == 2.0 \ A == 2.0 \ Matrix(A)

        @testset "trans-adj" begin
            A = brand(5,5,1,1)
            Ã = copy(A)

            lmul!(2.0, A')
            @test A == 2Ã
            lmul!(2.0, transpose(A))
            @test A == 4Ã
            lmul!(2.0, Symmetric(A))
            @test A == 8Ã
            lmul!(2.0, Hermitian(A))
            @test A == 16Ã
            rmul!(Hermitian(A), 1/2)
            @test A == 8Ã
            rmul!(Symmetric(A), 1/2)
            @test A == 4Ã
            rmul!(transpose(A), 1/2)
            @test A == 2Ã
            rmul!(A', 1/2)
            @test A == Ã
        end
    end

    @testset "axpy!" begin
        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)
        C .= A .+ B
        @test C == A + B == A .+ B  == Matrix(A) + Matrix(B)
        @test A + B isa BandedMatrix
        @test A .+ B isa BandedMatrix
        @test bandwidths(A+B) == bandwidths(A.+B) == (2,2)
        B .= A .+ B
        @test B == C

        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)
        C .= A .* B
        @test C == A .* B  == Matrix(A) .* Matrix(B)
        @test A .* B isa BandedMatrix
        @test bandwidths(A.*B) == (1,1)
        @time B .= A .* B
        @test B == C

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B[band(2)] .= B[band(-2)] .= 0
        C = brand(n,n,1,1)
        C .= A .+ B
        @test C == A + B == A .+ B == Matrix(A) + Matrix(B)

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B[band(-2)] .= 0
        C = brand(n,n,1,2)
        C .= A .+ B
        @test C == A + B == A .+ B == Matrix(A) + Matrix(B)

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        B[band(2)] .= 0
        C = brand(n,n,2,1)
        C .= A .+ B
        @test C == A + B == A .+ B == Matrix(A) + Matrix(B)

        n = 100
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)

        C .= 2.0 .* A .+ B
        @test C == 2A+B == 2.0.*A .+ B
        @test all(BLAS.axpy!(2.0, A, copy(B)) .=== C)

        @test 2A + B isa BandedMatrix
        @test 2.0.*A .+ B isa BandedMatrix
        @test bandwidths(2A+B) == bandwidths(2.0.*A .+ B) == (2,2)
        B .= 2.0 .* A .+ B
        @test B == C
    end

    @testset "gbmv!" begin
        n = 100
        x = randn(n)
        A = brand(n,n,1,1)
        y = similar(x)
        mul!(y,A,x)
        @test Matrix(A)*x ≈ y
        @test all(BLAS.gbmv!('N', n, 1, 1, 1.0, A.data, x, 0.0, copy(y)) .===
                    Broadcast.materialize!(MulAdd(1.0,A,x,0.0,similar(y))) .=== y)
        z = similar(y)
        z .= MulAdd(2.0,A,x,3.0,y)
        @test 2Matrix(A)*x + 3y ≈ z
        @test all(BLAS.gbmv!('N', n, 1, 1, 2.0, A.data, x, 3.0, copy(y)) .=== z)
        @test MemoryLayout(typeof(A')) == BandedRows{DenseColumnMajor}()
        mul!(y,A',x)
        @test Matrix(A')*x ≈ Matrix(A)'*x ≈ y
        @test all(BLAS.gbmv!('T', n, 1, 1, 1.0, A.data, x, 0.0, copy(y)) .=== y)
        z = similar(y)
        z .= MulAdd(2.0,A',x,3.0,y)
        @test 2Matrix(A')*x + 3y ≈ z
        @test all(BLAS.gbmv!('T', n, 1, 1, 2.0, A.data, x, 3.0, copy(y)) .=== z)
        @test MemoryLayout(typeof(transpose(A))) == BandedRows{DenseColumnMajor}()
        muladd!(1.0,transpose(A),x,0.0,y)
        @test Matrix(A')*x ≈ Matrix(A)'*x ≈ y
        @test all(BLAS.gbmv!('T', n, 1, 1, 1.0, A.data, x, 0.0, copy(y)) .=== y)

        n = 100
        x = randn(ComplexF64,n)
        A = brand(ComplexF64,n,n,1,1)
        y = similar(x)
        mul!(y,A,x)
        @test Matrix(A)*x ≈ y
        @test all(BLAS.gbmv!('N', n, 1, 1, 1.0+0.0im, A.data, x, 0.0+0.0im, copy(y)) .=== y)
        @test MemoryLayout(typeof(A')) == ConjLayout{BandedRows{DenseColumnMajor}}()
        z = similar(y)
        z .= MulAdd(2.0+0.0im,A,x,3.0+0.0im,y)
        @test 2Matrix(A)*x + 3y ≈ z
        @test all(BLAS.gbmv!('N', n, 1, 1, 2.0+0.0im, A.data, x, 3.0+0.0im, copy(y)) .=== z)
        mul!(y,A',x)
        @test Matrix(A')*x ≈ Matrix(A)'*x ≈ y
        @test all(BLAS.gbmv!('C', n, 1, 1, 1.0+0.0im, A.data, x, 0.0+0.0im, copy(y)) .=== y)
        @test MemoryLayout(typeof(transpose(A))) == BandedRows{DenseColumnMajor}()
        mul!(y,transpose(A),x)
        @test Matrix(transpose(A))*x ≈ transpose(Matrix(A))*x ≈ y
        @test all(BLAS.gbmv!('T', n, 1, 1, 1.0+0.0im, A.data, x, 0.0+0.0im, copy(y)) .=== y)
    end

    @testset "non-blas matvec" begin
        n = 10
        A = brand(Int,n,n,1,1)
        x = rand(Int,n)
        y = similar(x)
        mul!(y,A,x)
        @test y == Matrix(A)*x
        mul!(y,A',x)
        @test y == Matrix(A')*x
        mul!(y,transpose(A),x)
        @test y == Matrix(transpose(A))*x

        A = brand(Int,n,n,1,1) + im*brand(Int,n,n,1,1)
        x = rand(eltype(A),n)
        y = similar(x)
        mul!(y,A,x)
        @test y == Matrix(A)*x
        mul!(y,A',x)
        @test y == Matrix(A')*x
        mul!(y,transpose(A),x)
        @test y == Matrix(transpose(A))*x
    end

    @testset "gbmm!" begin
        n = 10
        A = brand(n,n,1,1)
        B = brand(n,n,2,2)
        C = brand(n,n,3,3)
        mul!(C,A,B)
        @test Matrix(C) ≈ Matrix(A)*Matrix(B)
        mul!(C,A',B)
        @test Matrix(C) ≈ Matrix(A')*Matrix(B)
        mul!(C,A,B')
        @test Matrix(C) ≈ Matrix(A)*Matrix(B')
        C = brand(n,n,4,4)
        mul!(C,A,B)
        @test Matrix(C) ≈ Matrix(A)*Matrix(B)
        B = randn(n,n)
        C = similar(B)

        mul!(C,A,B)
        @test C ≈ Matrix(A)*Matrix(B)
        mul!(C,B,A)
        @test C ≈ Matrix(B)*Matrix(A)
    end

    @testset "Subarray" begin
        n = 10
        A = brand(n,n,1,1)
        B = view(brand(2n,2n,2,0),Base.OneTo(n),1:n)
        @test BroadcastStyle(typeof(B)) == BandedStyle()

        @test A .+ B isa BandedMatrix
        @test A + B isa BandedMatrix
        @test bandwidths(A .+ B) == (2,1)
        @test A .+ B == Matrix(A) + Matrix(B)
    end

    @testset "adjtans" begin
        n = 10
        A = brand(n,n,1,1)

        @test BroadcastStyle(typeof(A')) isa BandedStyle
        @test BroadcastStyle(typeof(transpose(A))) isa BandedStyle

        @test A' .+ A isa BandedMatrix
        @test transpose(A) .+ A isa BandedMatrix
        @test A' .+ A == transpose(A) .+ A == Matrix(A)' .+ A
    end

    @testset "vector and matrix broadcastring" begin
        n = 10
        A = brand(n,n,1,2)
        b = randn(n)

        @test b .* A == b .* Matrix(A)
        @test b .* A isa BandedMatrix
        @test bandwidths(b .* A) == bandwidths(A)
        @test A .* b == Matrix(A) .* b
        @test A .* b isa BandedMatrix
        @test bandwidths(A .* b) == bandwidths(A)
        @test b .\ A == b .\ Matrix(A)
        @test b .\ A isa BandedMatrix
        @test bandwidths(b .\ A) == bandwidths(A)
        @test A ./ b == Matrix(A) ./ b
        @test A ./ b isa BandedMatrix
        @test bandwidths(A ./ b) == bandwidths(A)

        @test bandwidths(Broadcast.broadcasted(/, b, A)) == (9,9)
        @test isinf((b ./ A)[4,1])
        @test bandwidths(Broadcast.broadcasted(\, A,b)) == (9,9)
        @test isinf((A .\ b)[4,1])

        @test A .* b == Matrix(A) .* b
        @test bandwidths(A .* b) == bandwidths(A)
        @test A ./ b == Matrix(A) ./ b
        @test bandwidths(A ./ b) == bandwidths(A)
        @test isinf((A .\ b)[4,1])

        @test reshape(b,10,1) .* A isa BandedMatrix
        @test reshape(b,10,1) .* A == A .* reshape(b,10,1) == A .* b
        @test bandwidths(reshape(b,10,1) .* A) == bandwidths(A)
        @test bandwidths(A .* reshape(b,10,1)) == bandwidths(A)

        @test b .+ A == b .+ Matrix(A) == A .+ b

        @test bandwidths(broadcasted(+, A, b')) == (9,9)
        @test A .+ b' == b' .+ A == Matrix(A) .+ b'
        @test bandwidths(A .+ b') == bandwidths(b' .+ A)  == (9,9)
        @test A .* b' == b' .* A == Matrix(A) .* b'
        @test bandwidths(A .* b') == bandwidths(b' .* A) == bandwidths(A)

        @testset "nested broadcast" begin
            @test bandwidths((b ./ 2) .* A) == (1,2)
            @test (b ./ 2) .* A == (b ./ 2) .* Matrix(A)
        end
    end

    @testset "views" begin
        A = BandedMatrix(Ones(10,10),(1,1))
        B = 2A
        view(A,8:10,:) .= view(B,8:10,:)
        @test A[1,1] == 1.0
        @test A[6,7] == 1.0
        @test A[7,7] == 1.0
        @test A[8,7] == 2.0

        A = BandedMatrix{Int}(undef,(10,10),(2,2));  vec(A.data) .= (1:length(A.data));
        B = BandedMatrix{Int}(undef,(10,10),(1,1));  vec(B.data) .= (1:length(B.data));
        Ã = Matrix(A); B̃ = Matrix(B)
        view(A,5:10,:) .= view(B,5:10,:)
        view(Ã,5:10,:) .= view(B̃,5:10,:)
        @test Ã == A

        A = BandedMatrix{Int}(undef,(10,10),(2,1));  vec(A.data) .= (1:length(A.data));
        B = BandedMatrix{Int}(undef,(10,10),(1,2));  vec(B.data) .= (1:length(B.data));
        @test_throws BandError view(A,5:10,:) .= view(B,5:10,:)

        A = BandedMatrix{Int}(undef,(10,10),(1,2));  vec(A.data) .= (1:length(A.data));
        B = BandedMatrix{Int}(undef,(10,10),(2,1));  vec(B.data) .= (1:length(B.data));
        @test_throws BandError view(A,5:10,:) .= view(B,5:10,:)

        A = BandedMatrix{Int}(undef,(10,10),(1,2));  vec(A.data) .= (1:length(A.data));
        B = BandedMatrix{Int}(undef,(10,10),(2,1));  vec(B.data) .= (1:length(B.data));
        view(A,3:4,2:6) .= view(B,8:9,1:5)
        @test A[3:4,2:6] == zeros(Int,2,5)
        @test A[5,4] == 16

        A = BandedMatrix{Int}(undef,(10,10),(1,2));  vec(A.data) .= (1:length(A.data));
        B = BandedMatrix{Int}(undef,(10,10),(2,1));  vec(B.data) .= (1:length(B.data));
        view(A,3:4,2:6) .= view(B,1:2,4:8)
        @test A[3:4,:] == zeros(Int,2,10)

        A = BandedMatrix{Int}(undef,(1,1),(1,1)); vec(A.data) .= 2*(1:length(A.data));
        B = BandedMatrix{Int}(undef,(1,1),(1,1)); vec(B.data) .= (1:length(B.data));
        view(A,2:1,2:1) .= view(B,2:1,2:1)
        @test A[1,1] == 4
    end

    @testset "special broasdcast" begin
        A = brand(5,5,1,1)
        b = brand(5,1,1,1)
        B = brand(1,5,1,1)
        @test A .+ b == b .+ A == Matrix(A) .+ Matrix(b)
        @test A .+ B' == B' .+ A == Matrix(A) .+ Matrix(B)'
        @test A .+ b' == b' .+ A == Matrix(A) .+ Matrix(b)'
        @test A .+ B == B .+ A == Matrix(A) .+ Matrix(B)

        @test bandwidths(A .* b) == bandwidths(A .* b') == bandwidths(b .* A) == bandwidths(b' .* A) == (1,1)
        @test A .* b == b .* A == Matrix(A) .* Matrix(b)
        @test A .* b' == b' .* A == Matrix(A) .* Matrix(b')
        @test A .* B == B .* A == Matrix(A) .* Matrix(B)
    end

    @testset "broadcast with Int (#154)" begin
        b = BandedMatrix(8Ones(10,10), (2,1))
        @test BandedMatrices._isweakzero(round,Int,b)
        @test round.(Int, b) == round.(Int,Matrix(b))
    end

    @testset "degnerate bands" begin
		A = BandedMatrix(Zeros(1,1), (0,-3))
		C = BandedMatrix(Ones(1,1), (0,0))
		C .= 0.0 .- A
		@test C == zeros(1,1)
		C .= A .- 0.0
        @test C == zeros(1,1)

        A = BandedMatrix(Eye(3,4))
        B = randn(3,4)
        bc = broadcasted(*, A, B)

        @test copyto!(similar(bc, Float64), bc) == A .* B == B .* A ==
            BandedMatrix(B, (0,0)) == Matrix(A) .* B
    end

    @testset "adding degenerate" begin
        A = BandedMatrix(1 => 1:9)
        B = BandedMatrix(-1 => 1:9)
        C = BandedMatrix(Fill(-1,10,10),(1,1))
        C .= A .+ B
        @test diag(C) == Zeros{Int}(10)
        C = BandedMatrix(Fill(-1,10,10),(1,1))
        C .= A .+ B
        @test diag(C) == Zeros{Int}(10)
    end

    @testset "nested broadcast" begin
        A = brand(5,4,2,1)
        x = randn(5)
        B = Base.Broadcast.broadcasted(*, Base.Broadcast.broadcasted(+, 2, x), A)
        @test bandwidths(B) == bandwidths(A) == bandwidths(materialize(B)) == (2,1)
        @test materialize(B) == (2 .+ x) .* Matrix(A)

        B = Base.Broadcast.broadcasted(+, Base.Broadcast.broadcasted(+, A, A), A)
        @test bandwidths(B) == bandwidths(A)
        @test (A .+ A) .+ A == 3A
    end
end


text:
# methods for estimating dtw #

abstract type DTWDistance{D <: Union{Function, Distances.PreMetric}} <: Distances.SemiMetric end

"""
    struct DTW{D,N} <: DTWDistance{D}

# Keyword arguments:
- `radius`: The maximum allowed deviation of the matching path from the diagonal
- `dist`: Inner distance
- `transportcost` If >1, an additional penalty factor for non-diagonal moves is added.
- `normalizer`: defaults to `Nothing`

If the two time-series are of equal length, [`dtw_cost`](@ref) is called, if not, [`dtwnn`](@ref) is called.
"""
struct DTW{D,N,FK} <: DTWDistance{D}
    "The maximum allowed deviation of the matching path from the diagonal"
    radius::Int
    dist::D
    "If >1, an additional penalty factor for non-diagonal moves is added."
    transportcost::Float64
    postprocess::FK
end
DTW(r,dist::D=SqEuclidean(); transportcost=1, postprocess::FK=nothing, normalizer::Type{N}=Nothing) where {D,N,FK} = DTW{D, N, FK}(r,dist,transportcost,postprocess)
DTW(;radius,dist=SqEuclidean(), transportcost=1, postprocess::FK=nothing, normalizer::Type{N}=Nothing) where {N,FK} = DTW{typeof(dist), N, FK}(radius,dist,transportcost,postprocess)

"""
    struct SoftDTW{D, T} <: DTWDistance{D}

# Arguments:
- `γ`: smoothing parameter default to 1. Smaller value makes the distance closer to standard DTW.
- `dist`: Inner distance, defaults to `SqEuclidean()`.
- `transportcost`
"""
Base.@kwdef struct SoftDTW{D,T,R} <: DTWDistance{D}
    γ::T
    "The maximum allowed deviation of the matching path from the diagonal"
    dist::D = SqEuclidean()
    "If >1, an additional penalty factor for non-diagonal moves is added."
    transportcost::Float64 = 1.0
    radius::R = nothing
    SoftDTW(γ, dist=SqEuclidean(),transportcost=1,radius=nothing) = new{typeof(dist), typeof(γ), typeof(radius)}(γ,dist,transportcost,radius)
end

"""
    struct FastDTW{D} <: DTWDistance{D}

- `radius`
- `dist` inner distance
"""
Base.@kwdef struct FastDTW{D} <: DTWDistance{D}
    radius::Int
    dist::D = SqEuclidean()
    FastDTW(r, dist=SqEuclidean()) = new{typeof(dist)}(r, dist)
end

function Distances.evaluate(d::DTW{<:Any,N}, x, y; kwargs...) where N
    if lastlength(x) == lastlength(y)
        x, y = normalize(N, x), normalize(N, y)
        return dtw_cost(x, y, d.dist, d.radius; transportcost = d.transportcost, kwargs...)
    end
    if lastlength(x) > lastlength(y)
        x, y = y, x
    end
    dtwnn(
        x,
        y,
        d.dist,
        d.radius,
        N;
        prune_endpoints = false,
        transportcost = d.transportcost,
        kwargs...,
    ).cost
end

Distances.evaluate(d::SoftDTW, x, y) = soft_dtw_cost(x, y, d.dist, γ=d.γ, radius=d.radius)
Distances.evaluate(d::FastDTW, x, y) =
    fastdtw(x, y, d.dist, d.radius)[1]

distpath(d::DTW, x, y) = dtw(x, y, d.dist; transportcost=d.transportcost, postprocess = d.postprocess)
distpath(d::DTW, x, y, i2min::AbstractVector, i2max::AbstractVector; transportcost=d.transportcost) =
    dtw(x, y, i2min, i2max, d.dist)
distpath(d::FastDTW, x, y) = fastdtw(x, y, d.dist, d.radius)

(d::DTWDistance)(x,y;kwargs...) = Distances.evaluate(d,x,y;kwargs...)

"""
    distance_profile(d::DTWDistance, Q::AbstractArray{S}, T::AbstractArray{S}; kwargs...) where S

Optimized method for computing the distance profile using DTW distances. kwargs are sent to [`dtwnn`](@ref).
"""
function SlidingDistancesBase.distance_profile(d::DTW{<:Any,N}, Q::AbstractArray{S}, T::AbstractArray{S}; kwargs...) where {N,S}
    m = lastlength(Q)
    n = lastlength(T)
    n >= m || throw(ArgumentError("Q cannot be longer than T"))
    l = n-m+1
    res = dtwnn(Q, T, d.dist, d.radius, N; saveall=true, kwargs...)
    res.dists
end


text:
"""
Kibaek Kim - ANL MCS 2016
  Updated in 2020
Farmer example from Birge and Louveaux book.

This example writes farmer problem in a non-stochastic form
by explictly defining nonanticipativity constraints on x.
This form particularly represents a block-angular structure of the model.
In this form, each scenario index is interpreted to represent a sub-block.
All variables should be declared at the parent model.

To use MPI, add
- MPI.Init()
- DSPopt.parallelize(comm)
- MPI.Finalize()
"""

using MPI
using StructJuMP
using DSPopt

MPI.Init()

# Initialize DSPopt.jl with the communicator.
DSPopt.parallelize(MPI.COMM_WORLD)

NS = 3;                        # number of scenarios
probability = [1/3, 1/3, 1/3]; # probability

CROPS  = 1:3 # set of crops (wheat, corn and sugar beets, resp.)
PURCH  = 1:2 # set of crops to purchase (wheat and corn, resp.)
SELL   = 1:4 # set of crops to sell (wheat, corn, sugar beets under 6K and those over 6K)

Cost     = [150 230 260]   # cost of planting crops
Budget   = 500             # budget capacity
Purchase = [238 210];      # purchase price
Sell     = [170 150 36 10] # selling price
Yield    = [3.0 3.6 24.0;
            2.5 3.0 20.0;
            2.0 2.4 16.0]
Minreq   = [200 240 0]     # minimum crop requirement

m = StructuredModel(num_scenarios = NS)

@variable(m, x[i=CROPS,s=1:NS] >= 0, Int)
@variable(m, y[j=PURCH,s=1:NS] >= 0)
@variable(m, w[k=SELL,s=1:NS] >= 0)
@objective(m, Min,
	  sum(probability[s] * Cost[i] * x[i,s] for i=CROPS for s=1:NS)
    + sum(probability[s] * Purchase[j] * y[j,s] for j=PURCH for s=1:NS)
	- sum(probability[s] * Sell[k] * w[k,s] for k=SELL for s=1:NS))
@constraint(m, nonant[i=CROPS,s=2:NS], x[i,s-1] - x[i,s] == 0)

# TODO: Distributing blocks is not supported yet.
# This is because of the DSPopt implementation.
for s in 1:NS
    blk = StructuredModel(parent = m, id = s)
	@objective(blk, Min, 0.)
	@constraint(blk, const_budget, sum(x[i,s] for i=CROPS) <= Budget)
    @constraint(blk, const_minreq[j=PURCH], Yield[s,j] * x[j,s] + y[j,s] - w[j,s] >= Minreq[j])
    @constraint(blk, const_minreq_beets, Yield[s,3] * x[3,s] - w[3,s] - w[4,s] >= Minreq[3])
    @constraint(blk, const_aux, w[3,s] <= 6000)
end

status = optimize!(m,
    is_stochastic = false, # Needs to indicate that the model is NOT a stochastic program.
    solve_type = DSPopt.DW, # see instances(DSPopt.Methods) for other methods
    param = "examples/params.txt" # This path assumes running from the one-level upper directory (i.e., ../).
    )

if DSPopt.myrank() == 0 && status == MOI.OPTIMAL
    @show objective_value(m)
    @show dual_objective_value(m)
    @show value.(x)
    @show value.(y)
    @show value.(w)
    @show dual() # This is available only for solve_type = DSPopt.Dual.
end


text:
####
#### Load DaemonMode
####

# Load the Julia package manager, which is itself a Julia package
import Pkg

# Load the DaemonMode package
try
    @eval using DaemonMode
catch e
    println("* Installing DaemonMode package for DECAES server")
    Pkg.add("DaemonMode"; io = devnull)
    @eval using DaemonMode
end

# Server port number. This can be changed to any valid port
const PORT = DaemonMode.PORT

####
#### Utilities
####

# Julia command for starting server.
# Note: inherits current Julia environment - make sure number of threads
#       are set before running this script!
julia_cmd =
    ```
    $(Base.julia_cmd())
    --project=$(Base.active_project())
    --threads=$(Threads.nthreads())
    --startup-file=no
    --optimize=3
    --quiet
    ```;

# Convenience macro to write a Julia expression into a temporary script
macro mktempscript(ex = nothing)
    quote
        local fname = tempname() * ".jl"
        open(fname; write = true) do io
            println(io, "const PORT = $PORT")
            println(io, $(string(ex)))
        end
        fname
    end
end

####
#### Start DECAES server, if necessary
####

# Check for kill server command
if !isempty(ARGS) && lowercase(ARGS[1]) == "--kill"
    try
        sendExitCode(PORT)
        println("* DECAES server killed")
        exit()
    catch e
        if e isa Base.IOError
            println("* DECAES server inactive; nothing to kill")
            exit()
        else
            rethrow(e)
        end
    end
end

# Ping server by trying to delete dummy file
function ping()
    daemon_script = @mktempscript begin
        using DaemonMode
        redirect_stdout(devnull) do
            redirect_stderr(devnull) do
                runargs(PORT)
            end
        end
    end
    ping_script = @mktempscript rm(ARGS[1]; force = true)
    ping_file = touch(tempname())
    try
        run(`$(julia_cmd) $(daemon_script) $(ping_script) $(ping_file)`)
        return !isfile(ping_file)
    catch e
        if e isa ProcessFailedException
            return false
        else
            rethrow(e)
        end
    finally
        rm.([ping_file, daemon_script, ping_script]; force = true)
    end
end

if !ping()
    # Server is not started
    println("* Starting DECAES server")

    server_script = @mktempscript begin
        using DaemonMode
        serve(PORT)
    end
    server_cmd = `$(julia_cmd) $(server_script) \&`
    run(detach(server_cmd); wait = false)

    while !ping()
        sleep(1)
    end
end

####
#### Run DECAES on server
####

decaes_script = @mktempscript begin
    import Pkg
    try
        @eval using DECAES
    catch e
        println("* Installing DECAES")
        # Pkg.add("DECAES"; io = devnull)
        Pkg.develop("DECAES"; io = devnull) #TODO
        @eval using DECAES
    end

    # Call the command line interface entrypoint function
    main()
end

pushfirst!(ARGS, decaes_script)
runargs(PORT)


text:
using Test
import DynamicPolynomials: @polyvar
# using Revise
using PolynomialBasis

PB = PolynomialBasis

@polyvar x

function allequal(v1,v2,tol)
    return allequal(v1,v2;tol=tol)
end

function allequal(v1,v2;tol=1e3eps())
    np = length(v1)
    f = length(v2) == np
    return f && all([isapprox(v1[i],v2[i],atol=tol) for i = 1:np])
end

basis = PB.CubicHermiteBasis()
@test allequal(basis.points,[-1.0,1.0])
@test allequal(basis(-1.0),[1.0,0.,0.,0.])
@test allequal(basis(1.0),[0.0,0.,1.,0.])
@test allequal(PB.derivative(basis,-1.0),[0.,1.,0.,0.])
@test allequal(PB.derivative(basis,1.0),[0.,0.,0.,1.])
@test allequal(PB.gradient(basis,-1.0),[0.,1.,0.,0.])
@test allequal(PB.gradient(basis,1.0),[0.,0.,0.,1.])

hess = PB.hessian(basis,-1.0)
@test allequal(hess,[-1.5,-2.0,1.5,-1.0])
hess = PB.hessian(basis,1.0)
@test allequal(hess,[1.5,1,-1.5,2])


text:
RpcApi.SimpleOrganism(0x00000000000009f4,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100294,Dict{Int16,Int16}(),[136,85])
RpcApi.SimpleOrganism(0x00000000000009f5,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,99994,Dict{Int16,Int16}(),[128,45])
RpcApi.SimpleOrganism(0x00000000000009f6,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100594,Dict{Int16,Int16}(),[112,170])
RpcApi.SimpleOrganism(0x00000000000009f7,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,99394,Dict{Int16,Int16}(),[119,7])
RpcApi.SimpleOrganism(0x00000000000009f8,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepLeft(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,99994,Dict{Int16,Int16}(),[87,117])
RpcApi.SimpleOrganism(0x00000000000009f9,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_2(var_11::Bool=true,var_12::Int64=104785325214659113)
            Creature.eatUp(o)
            return var_11
        end
        Creature.eatUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepRight(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        var_8 = var_8 * var_8
        Creature.eatUp(o)
        Creature.eatDown(o)
        Creature.eatDown(o)
        Creature.eatDown(o)
        Creature.eatUp(o)
        Creature.eatDown(o)
        Creature.eatLeft(o)
        var_8 = Bool(abs(var_8 - var_8))
        var_8 = func_1(var_8,var_8)
        var_3 = var_3 * var_3
        o.mem[var_9] = var_9
        var_8 = Bool(abs(var_8 - var_8))
        Creature.stepDown(o)
        local var_8::Bool = false
        Creature.stepDown(o)
        local var_9::Int16 = -8672
        o.mem[var_2] = var_2
        if var_8 !== var_8
            Creature.eatUp(o)
            Creature.eatDown(o)
            var_8 = Bool(abs(var_8 - var_8))
        end
    end),15,59392,63,100794,Dict(-5900=>-5900),[300,55])
RpcApi.SimpleOrganism(0x00000000000009fa,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100594,Dict{Int16,Int16}(),[105,117])
RpcApi.SimpleOrganism(0x00000000000009fb,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100194,Dict{Int16,Int16}(),[56,176])
RpcApi.SimpleOrganism(0x00000000000009fc,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,96194,Dict{Int16,Int16}(),[90,175])
RpcApi.SimpleOrganism(0x00000000000009fd,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100494,Dict{Int16,Int16}(),[135,135])
RpcApi.SimpleOrganism(0x00000000000009fe,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepUp(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,88094,Dict{Int16,Int16}(),[133,28])
RpcApi.SimpleOrganism(0x00000000000009ff,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,98094,Dict{Int16,Int16}(),[104,9])
RpcApi.SimpleOrganism(0x0000000000000a00,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_2(var_5::Bool=false,var_6::Int8=-81,var_7::Int16=28853,var_8::Int64=-2414801820824050205)
            var_8 = Int(Creature.getEnergy(o,Int(var_7),Int(var_7)))
            return var_5
        end
        Creature.stepDown(o)
        local var_2::Int16 = -5900
        Creature.stepDown(o)
        var_2 = if haskey(o.mem,var_2)
                o.mem[var_2]
            else
                var_2
            end
        Creature.eatUp(o)
        Creature.stepUp(o)
    end),1,1000,1,28894,Dict{Int16,Int16}(),[17,1])
RpcApi.SimpleOrganism(0x0000000000000a01,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            Creature.stepRight(o)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,99894,Dict{Int16,Int16}(),[91,167])
RpcApi.SimpleOrganism(0x0000000000000a02,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_2(var_7::Int8=-65,var_8::Int8=112,var_9::Int8=40)
            Creature.stepUp(o)
            return var_7
        end
        local var_6::Int16 = -31006
        function func_1(var_5::Int16=25063)
            Creature.stepUp(o)
            o.mem[var_5] = var_5
            Creature.stepUp(o)
            return var_5
        end
        o.mem[var_6] = var_6
        Creature.eatRight(o)
        Creature.eatUp(o)
        Creature.stepRight(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
        Creature.stepLeft(o)
        local var_3::Int8 = 60
        Creature.stepUp(o)
    end),31,10616,19,83094,Dict(-31006=>-31006),[22,139])
RpcApi.SimpleOrganism(0x0000000000000a03,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        o.mem[var_2] = var_2
        local var_2::Int16 = -5900
        Creature.eatRight(o)
        Creature.stepDown(o)
        Creature.eatLeft(o)
        Creature.stepUp(o)
        Creature.stepDown(o)
        Creature.eatUp(o)
    end),1,1000,1,100093,Dict{Int16,Int16}(),[24,55])
RpcApi.SimpleOrganism(0x0000000000000a04,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100093,Dict{Int16,Int16}(),[77,143])
RpcApi.SimpleOrganism(0x0000000000000a05,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        Creature.stepDown(o)
        local var_2::Int16 = -5900
        Creature.stepDown(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepUp(o)
    end),1,1000,1,101893,Dict{Int16,Int16}(),[17,2])
RpcApi.SimpleOrganism(0x0000000000000a06,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,31693,Dict{Int16,Int16}(),[115,189])
RpcApi.SimpleOrganism(0x0000000000000a07,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100293,Dict{Int16,Int16}(),[98,57])
RpcApi.SimpleOrganism(0x0000000000000a08,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,99693,Dict{Int16,Int16}(),[123,200])
RpcApi.SimpleOrganism(0x0000000000000a09,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        o.mem[var_2] = var_2
        local var_2::Int16 = -5900
        Creature.eatRight(o)
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepDown(o)
        Creature.eatUp(o)
        Creature.eatUp(o)
    end),1,1000,1,100093,Dict{Int16,Int16}(),[24,56])
RpcApi.SimpleOrganism(0x0000000000000a0a,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        Creature.eatLeft(o)
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,97993,Dict{Int16,Int16}(),[120,3])
RpcApi.SimpleOrganism(0x0000000000000a0b,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        Creature.eatUp(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        function func_1(var_4::Bool=true,var_5::Bool=false)
            begin  # g:\my\projects\jevo\src\Code.jl, line 68:
                local i::Int8 # g:\my\projects\jevo\src\Code.jl, line 68:
                for i = 1:var_7
                    var_7 = var_7 / var_7
                    Creature.eatUp(o)
                    var_7 = var_7 - var_7
                    var_5 = var_5 | var_4
                    Creature.eatLeft(o)
                end
            end
            local var_7::Int8 = 72
            Creature.eatLeft(o)
            Creature.stepDown(o)
            if var_7 !== var_7
                Creature.stepUp(o)
            end
            var_7 = var_7 * var_7
            Creature.eatLeft(o)
            Creature.eatUp(o)
            Creature.eatLeft(o)
            var_7 = var_7 + var_7
            return var_4
        end
        Creature.eatRight(o)
        Creature.stepUp(o)
        Creature.eatDown(o)
        var_8 = var_8 | var_8
        Creature.eatDown(o)
        Creature.eatDown(o)
        Creature.stepLeft(o)
        var_3 = var_3 * var_3
        var_8 = Bool(abs(var_8 - var_8))
        Creature.stepDown(o)
        local var_8::Bool = false
        o.mem[var_2] = var_2
        if var_8 !== var_8
        end
    end),15,59392,63,99993,Dict(-5900=>-5900),[42,6])
false
RpcApi.SimpleOrganism(0x0000000000000a0d,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        Creature.eatUp(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        function func_1(var_4::Bool=true,var_5::Bool=false)
            var_7 = var_7 / var_7
            local var_7::Int8 = 72
            Creature.eatLeft(o)
            Creature.stepDown(o)
            Creature.stepDown(o)
            if var_7 !== var_7
                Creature.eatLeft(o)
                local var_9::Bool = true
            end
            var_7 = var_7 * var_7
            Creature.eatLeft(o)
            Creature.eatUp(o)
            Creature.eatLeft(o)
            var_7 = var_7 + var_7
            return var_4
        end
        Creature.eatRight(o)
        Creature.stepUp(o)
        Creature.eatDown(o)
        var_8 = var_8 | var_8
        Creature.eatDown(o)
        Creature.eatDown(o)
        Creature.stepLeft(o)
        var_3 = var_3 * var_3
        var_8 = Bool(abs(var_8 - var_8))
        Creature.stepDown(o)
        local var_8::Bool = false
        Creature.eatRight(o)
    end),15,59392,63,99993,Dict(-5900=>-5900),[42,28])
RpcApi.SimpleOrganism(0x0000000000000a0e,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        o.mem[var_2] = var_2
        local var_2::Int16 = -5900
        Creature.eatRight(o)
        Creature.stepDown(o)
        Creature.eatLeft(o)
        Creature.stepUp(o)
        Creature.stepDown(o)
        Creature.eatUp(o)
    end),1,1000,1,100093,Dict{Int16,Int16}(),[23,56])
RpcApi.SimpleOrganism(0x0000000000000a0f,[5,50,100,4,1,1,1],:(function (o::Creature.Organism,)
        function func_1(var_3::Int8=18)
            return var_3
        end
        Creature.stepDown(o)
        Creature.stepUp(o)
        Creature.stepLeft(o)
        Creature.stepUp(o)
        Creature.stepRight(o)
        Creature.eatUp(o)
        Creature.stepDown(o)
        Creature.stepDown(o)
    end),1,1000,1,100093,Dict{Int16,Int16}(),[169,8])

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# debug.jl

@time versioninfo()

println("\n")

@time whos()

max_stars_repo_name: anthony-wang/julia-scripts

max_issues_repo_name: anthony-wang/julia-scripts

max_forks_repo_name: anthony-wang/julia-scripts


text:
using LightApp: App, State, Button, h, serve

app = App()

app.state.x = 1

app.layout = h.div(
    h.h1("Hello World!"),

    h.p("My app has state: ", State("x")),

    Button(children="Increment X", state="x") do state, val
        state.x += 1
    end
)

serve(app)


text:
# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking, CmdStan
#gr(size=(600,600));

# CmdStan uses a tmp directory to store the output of cmdstan

ProjDir = rel_path("..", "scripts", "04")
cd(ProjDir)

# ### snippet 4.7

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

# ### snippet 4.8

# Use only adults

df2 = filter(row -> row[:age] >= 18, df);

# Show first 5 rows of DataFrame df

first(df2, 5)

# ### Snippet 4.14

# Generate approximate probabilities

struct Post
  mu::Float64
  sigma::Float64
  ll:: Float64
  prod::Float64
  prob::Float64
end

mu_list = repeat(range(140, 160, length=200), 200);
sigma_list = repeat(range(4, 9, length=200), inner=200);

ll = zeros(40000);
for i in 1:40000
    d1 = Normal(mu_list[i], sigma_list[i])
    ll[i] = sum(log.(pdf.(d1, df2[:height])))
end

d2 = Normal(178.0, 20.0)
d3 = Uniform(0, 50)
prod = ll + log.(pdf.(d2, mu_list)) + log.(pdf.(d3, sigma_list))
prob = exp.(prod .- maximum(prod))
post = DataFrame(mu=mu_list, sigma=sigma_list, ll=ll, prod=prod, prob=prob)
first(post, 10)

# ### Snippet 4.15

# Sample post

samples = post[sample(1:size(post, 1), Weights(post[:prob]), 10000, replace=true), :]

# ### Snippet 4.19

# Density of mu

density(samples[:mu])

# Density of sigma

density(samples[:sigma])

# ### Snippet 4.20

# Hdp mu

#hpd(samples[:mu])

# Hdp sigma

#hpd(samples[:sigma])

# End of `clip-14-20.jl`


text:
export NLNS, NSAR2, AR1, randomwalk, SNLST, random_cycles

"""
    AR1(n_steps, x₀, k)

Simple AR(1) model with no static transformation[^1].

## Equations

The system is given by the following map:
```math
x(t+1) = k x(t) + a(t),
```
where ``a(t)`` is a draw from a normal distribution with zero mean and unit
variance. `x₀` sets the initial condition and `k` is the tunable parameter in
the map.

# References

[^1]: Lucio et al., Phys. Rev. E *85*, 056202 (2012). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)
"""
function AR1(n_steps, x₀, k)
    a = rand(Normal(), n_steps)
    x = Vector{Float64}(undef, n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = k*x[i-1] + a[i]
    end
    x
end

"""
    NSAR2(n_steps, x₀, x₁)

Cyclostationary AR(2) process[^1].

## References

[^1]: Lucio et al., Phys. Rev. E *85*, 056202 (2012). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)
"""
function NSAR2(n_steps, x₀, x₁)
    T₀ = 50.0
    τ = 10.0
    M = 5.5
    Tmod = 250
    σ₀ = 1.0
    a_1_0 = 2*cos(2*π/T₀)*exp(-1/τ)
    a₂ = -exp(-2 / τ)
    a = rand(Normal(), n_steps)
    x = Vector{Float64}(undef, n_steps)
    x[1], x[2] = x₀, x₁

    for i = 3:n_steps
        T_t = T₀ + M*sin(2*π*i / Tmod)
        a₁_t = 2*cos(2*π/T_t) * exp(-1 / τ)

        σ = σ₀^2/(1 - a_1_0^2 - a₂^2 - (2*(a_1_0^2) * a₂)/(1 - a₂)) *
            (1 - a₁_t - a₂^2 - (2*(a₁_t^2) * a₂)/(1 - a₂))
        x[i] = a₁_t * x[i-1] + a₂*x[i-2] + rand(Normal(0, sqrt(σ)))
    end
    x
end

"""
    randomwalk(n_steps, x₀)

Linear random walk (AR(1) process with a unit root)[^1].
This is an example of a nonstationary linear process.

# References

[^1]: Lucio et al., Phys. Rev. E *85*, 056202 (2012). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)
"""
function randomwalk(n_steps, x₀)
    a = rand(Normal(), n_steps)

    x = Vector{Float64}(undef, n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = x[i-1] + a[i]
    end
    x
end

"""
    SNLST(n_steps, x₀, k)

Dynamically linear process transformed by a strongly nonlinear static
transformation (SNLST)[^1].

## Equations
The system is by the following map:
```math
x(t) = k x(t-1) + a(t)
```
with the transformation ``s(t) = x(t)^3``.

# References

[^1]: Lucio et al., Phys. Rev. E *85*, 056202 (2012). [https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.85.056202)
"""
function SNLST(n_steps, x₀, k)
    a = rand(Normal(), n_steps)
    x = Vector{Float64}(undef, n_steps)
    x[1] = x₀
    for i = 2:n_steps
        x[i] = k*x[i-1] + a[i]
    end
    x.^3
end

# Keyword versions of the functions
AR1(;n_steps = 500, x₀ = rand(), k = rand())    = AR1(n_steps, x₀, k)
SNLST(;n_steps = 500, x₀ = rand(), k = rand())  = SNLST(n_steps, x₀, k)
randomwalk(;n_steps = 500, x₀ = rand())         = randomwalk(n_steps, x₀)
NSAR2(;n_steps = 500, x₀ = rand(), x₁ = rand()) = NSAR2(n_steps, x₀, x₁)

"""
    random_cycles(; periods=10 dt=π/20, σ = 0.05, frange = (0.8, 2.0))
Make a timeseries that is composed of `period` full sine wave periods, each with a
random frequency in the range given by `frange`, and added noise with std `σ`.
The sampling time is `dt`.
"""
function random_cycles(; periods=10, dt=π/20, σ = 0.05, frange = (1.0, 1.6))
    dt = π/20
    x = Float64[]

    for i in 1:periods
        f = (frange[1]-frange[2])*rand() + frange[1]
        T = 2π/f
        t = 0:dt:T
        append!(x, sin.(f .* t))
    end

    N = length(x)
    x .+= randn(N)/20
    return x
end


text:
using .AxisIndices: AbstractAxis,AxisIndices

valfromaxis(ax::AbstractAxis) = keys(ax)

getdata(a::AxisIndices.AxisArray) = parent(a)

yaxcreate(::Type{<:AxisIndices.AxisArray}, data, dnames, dvals, atts) =
  AxisIndices.AxisArray(data, map(i->dvals[i], 1:ndims(data))...)


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

import LinearAlgebra: factorize

const sbas = [
    :conj,
    :abs,
    :factorial,
    :floor,
    :max,
    :min,
    :round,
    :sign,
    :acos,
    :acosh,
    :acot,
    :acoth,
    :acsc,
    :acsch,
    :asec,
    :asech,
    :asin,
    :asinh,
    :atan,
    :atanh,
    #:atan2,
    :cos,
    :cosh,
    :cot,
    :coth,
    :csc,
    :csch,
    :exp,
    :hypot,
    :log,
    :log10,
    :sec,
    :sech,
    :sin,
    :sinh,
    :sqrt,
    :tan,
    :tanh,
]

const sdep = [
    :gamma,
    :beta,
    :besseli,
    :besselj,
    :besselk,
    :bessely,
    :polygamma,
    :zeta
]

const sfun = [
    :ibeta,
    :igamma,
    :ln,
    :psi,
    :bernoulli,
    :continued_fraction,
    :ci, #
    :dilog,
    :ei,
    :si,
    :airy_ai,
    :airy_aiprime,
    :airy_bi,
    :airy_biprime,
    :hanekl1,
    :hankel2,
    :kummerm,
    :kummeru,
    :lommel1,
    :lommel2,
    :struveh,
    :struvel,
    :whittakerm,
    :whittakeru,
    :solidharmonicy,
    :sphericalharmonicy,
    :expand_cases,
    :arglength,
    :decompose,
    :num,
    :den,
    :mainvar,
    :precision,
    :setmod,
    :rootval,
    :showrules,
    :reverse,
    :lhs,
    :rhs,
    :saveas,
    :root_val,
]

const snan = [
    :rlet,
    :clearrules,
    :scientific_notation,
    :optimize,
]

const snum = [
    :ceiling,
    :fix,
]

const scom = [
    :impart,
    :repart,
]

const sint = [
    :nextprime,
    :euler,
    :fibonacci,
    :motzkin,
]

const sran = [
    :random,
    :random_new_seed
]

const sbat = [
    :det,
    :trace,
    :nullspace,
    :rank
]

const smat = [
    :tp,
]

Expr(:block,[:($i(r)=Base.$i(r)) for i ∈ [sbas[2:end];sdep;sbat;[:length]]]...) |> eval
#Expr(:toplevel,[:(import Base: $i) for i ∈ [sbas;sdep;sbat;[:length]]]...) |> eval
:(export $([sbas;sdep;sfun;snan;snum;scom;sint;sran;sbat;smat;[:length]]...)) |> eval
#:(export $(Symbol.("@",[sbas;sdep;sfun;snum;scom;sint])...)) |> eval

for fun in [sbas;sdep;sfun;snum;scom;sint;sran;sbat;smat;[:length]]
    Reduce.parsegen(fun,:unary) |> eval
end

for fun in [sbas;sdep;sfun;snum;scom;sint]
    @eval begin
        $(Reduce.unfoldgen(fun,:unary))
        function $fun(expr::String;be=0)
            convert(String, $fun(RExpr(expr);be=be))
        end
        #=macro $fun(expr)
            :($$(QuoteNode(fun))($(esc(expr))))
        end=#
    end
end

for fun in [sbat;smat]
    @eval begin
        function $fun(expr::Union{Array{T,2},T}) where T <: ExprSymbol
            $fun(RExpr(expr)) |> parse |> mat
        end
        function $fun(expr::Array{Any,2})
            $fun(RExpr(expr)) |> parse |> mat
        end
    end
end

"""
    tp(exprn)

Syntax:
```
        tp(EXPRN:matrix_expression):matrix.
```
This operator takes a single matrix argument and returns its transpose.
"""
tp(r::Array{T,1}) where T <: ExprSymbol = r |> RExpr |> tp |> parse |> mat
tp(r::Union{Vector,Adjoint}) = r |> RExpr |> tp |> parse |> mat

"""
    length(r)

`length` is a generic operator for finding the length of various objects in the system. The meaning depends on the type of the object. In particular, the length of an algebraic expression is the number of additive top-level terms its expanded representation.

*Examples:*
```Julia
julia> length(:(a+b))
2

julia> length(2)
1
```
Other objects that support a length operator include arrays, lists and matrices. The explicit meaning in these cases is included in the description of these objects.
"""
length(r::Expr) = length(r |> RExpr) |> parse |> eval

for fun in snan
    nfun = fun ≠ :rlet ? fun : :let
    @eval begin
        $fun(r::RExpr) = string($(string(nfun)),"(",string(r),")") |> rcall |> RExpr
        $fun(r) = $fun(RExpr(r)) |> parse
    end
end

"""
    rlet(::Union{Dict,Pair},expr)

The simplest use of the `let` statement is in the form
```
R"let ⟨substitution list⟩"
```
where `⟨substitution list⟩` is a list of rules separated by commas, each of the form:
```
⟨variable⟩ = ⟨expression⟩
```
or
```
⟨prefix operator⟩(⟨argument⟩,…,⟨argument⟩) = ⟨expression⟩
```
or
```
⟨argument⟩⟨infix operator⟩,…,⟨argument⟩ = ⟨expression⟩
```
For example,
```
        let {x => y^2,
             h(u,v) => u - v,
             cos(pi/3) => 1/2,
             a*b => c,
             l+m => n,
             w^3 => 2*z - 3,
             z^10 => 0}
```
The list brackets can be left out if preferred. The above rules could also have been entered as seven separate `let` statements.

After such `let` rules have been input, `x` will always be evaluated as the square of `y`, and so on. This is so even if at the time the `let` rule was input, the variable `y` had a value other than `y`. (In contrast, the assignment `R"x:=y^2"` will set `x` equal to the square of the current value of `y`, which could be quite different.)

The rule `let a*b=c` means that whenever `a` and `b` are both factors in an expression their product will be replaced by `c`. For example, `a^5*b^7*w` would be replaced by `c^5*b^2*w`.

The rule for `l+m` will not only replace all occurrences of `l+m` by `n`, but will also normally replace `l` by `n-m`, but not `m` by `n-l`. A more complete description of this case is given in Section 11.2.5.

The rule pertaining to `w^3` will apply to any power of `w` greater than or equal to the third.

Note especially the last example, `let z^10=0`. This declaration means, in effect: ignore the tenth or any higher power of `z`. Such declarations, when appropriate, often speed up a computation to a considerable degree. (See Section 11.4 for more details.)

Any new operators occurring in such `let` rules will be automatically declared `operator` by the system, if the rules are being read from a file. If they are being entered interactively, the system will ask `Declare… Operator?`. Answer `Y` or `N` and hit `<Return>`.

In each of these examples, substitutions are only made for the explicit expressions given; i.e., none of the variables may be considered arbitrary in any sense. For example, the command
```Julia
julia> Algebra.rlet( :(h(u,v)) => :(u - v) )
```
will cause `h(u,v)` to evaluate to `u - v`, but will not affect `h(u,z)` or `h` with any arguments other than precisely the symbols `u,v`.

These simple `let` rules are on the same logical level as assignments made with the `:=` operator. An assignment `R"x := p+q"` cancels a rule `rlet( :x => :(y^2) )` made earlier, and vice versa.
"""
rlet(r::Dict{String,String}) = rlet(sub_list(r))
rlet(s::Dict{<:Any,<:Any}) = rlet(Dict([=>(string.(RExpr.([b[1],b[2]]))...) for b ∈ collect(s)]...))
rlet(s::Pair{<:Any,<:Any}) = rlet(Tuple([s]))
rlet(s::T) where T <: Tuple = rlet(RExpr(s)) |> parse
rlet(s::Array{<:Pair{<:Any,<:Any},1}) = rlet(Tuple(s))

"""
    Algebra.scientific_notation(::Union{Number,Tuple,Vector})

The declaration `scientific_notation` controls the output format of floating point numbers. At the default settings, any number with five or less digits before the decimal point is printed in a fixed-point notation, e.g., `12345.6`. Numbers with more than five digits are printed in scientific notation, e.g., `1.234567E+5`. Similarly, by default, any number with eleven or more zeros after the decimal point is printed in scientific notation. To change these defaults, `scientific_notation` can be used in one of two ways.
```Julia
julia> Algebra.scientific_notation(m);
```
where `m` is a positive integer, sets the printing format so that a number with more than `m` digits before the decimal point, or `m` or more zeros after the decimal point, is printed in scientific notation.
```Julia
julia> Algebra.scientific_notation(m,n);
```
with `m` and `n` both positive integers, sets the format so that a number with more than `m` digits before the decimal point, or `n` or more zeros after the decimal point is printed in scientific notation.
"""
scientific_notation(r::T) where T <: Tuple = scientific_notation(RExpr(r)) |> parse
scientific_notation(r::Union{Vector,Adjoint}) = scientific_notation(list(r)) |> parse
scientific_notation(r::T...) where T <: Number = scientific_notation(list(r)) |> parse

for fun in [sint;sran]
    @eval begin
        function $fun(n::T) where T <: Integer
            convert(T, $fun(RExpr(n)) |> parse |> eval)
        end
    end
end

for fun in snum
    @eval begin
        function $fun(x::T) where T <: Real
            $fun(RExpr(x)) |> parse |> eval
        end
    end
end

for fun in scom
    @eval begin
        function $fun(x::T) where T <: Number
            "$x" |> Meta.parse |> RExpr |> $fun |> parse |> eval
        end
    end
end

function bernoulli(n::T) where T <: Integer
    bernoulli(RExpr(n)) |> parse |> eval
end

@doc """
    lhs(::Union{Expr,RExpr})

Returns the left-hand side of an equation.

## Examples
```Julia
julia> Algebra.lhs(R"a+b=c")

a + b

```
""" Reduce.Algebra.lhs

@doc """
    rhs(::Union{Expr,RExpr})

Returns the right-hand side of an equation.

## Examples
```Julia
julia> Algebra.rhs(R"a+b=c")

c

```
""" Reduce.Algebra.rhs

@doc """
    abs(r)

`abs` returns the absolute value of its single argument, if that argument has a numerical value. A non-numerical argument is returned as an absolute value, with an overall numerical coefficient taken outside the absolute value operator. For example:

```Julia
julia> Algebra.abs(-3/4)
0.75

julia> Algebra.abs(:(2a))
:(2 * abs(a))

julia> Algebra.abs(im)
1.0

julia> Algebra.abs(:(-x))
:(abs(x))
```
""" Reduce.Algebra.abs

@doc """
    ceiling(r)

This operator returns the ceiling (i.e., the least integer greater than the given argument) if its single argument has a numerical value. A non-numerical argument is returned as an expression in the original operator. For example:

```Julia
julia> Algebra.ceiling(-5/4)
-1

julia> Algebra.ceiling(:(-a))
:(ceiling(-a))
```
""" Reduce.Algebra.ceiling

@doc """
    conj(r)

This operator returns the ceiling (i.e., the least integer greater than the given argument) if its single argument has a numerical value. A non-numerical argument is returned as an expression in the original operator. For example:

```Julia
julia> Algebra.conj(1+im)
1 - 1im

julia> Algebra.conj(:(a+im*b))
:(repart(a) - ((impart(a) + repart(b)) * im + impart(b)))
```
""" Reduce.Algebra.conj

@doc """
    factorial(r)

If the single argument of `factorial` evaluates to a non-negative integer, its factorial is returned. Otherwise an expression involving `factorial` is returned. For example:

```Julia
julia> Algebra.factorial(5)
120

julia> Algebra.factorial(:a)
:(factorial(a))
```
""" Reduce.Algebra.factorial

@doc """
    fix(r)

This operator returns the fixed value (i.e., the integer part of the given argument) if its single argument has a numerical value. A non-numerical argument is returned as an expression in the original operator. For example:

```Julia
julia> Algebra.fix(-5/4)
-1

julia> Algebra.fix(:a)
:(fix(a))
```
""" Reduce.Algebra.fix

@doc """
    floor(r)

This operator returns the floor (i.e., the greatest integer less than the given argument) if its single argument has a numerical value. A non-numerical argument is returned as an expression in the original operator. For example:

```Julia
julia> Algebra.floor(-5/4)
-2.0

julia> Algebra.floor(:a)
:(floor(a))
```
""" Reduce.Algebra.floor

@doc """
    impart(r)

This operator returns the imaginary part of an expression, if that argument has an numerical value. A non-numerical argument is returned as an expression in the operators `repart` and `impart`. For example:

```Julia
julia> Algebra.impart(1+im)
1

julia> Algebra.impart(:(a+im*b))
:(impart(a) + repart(b))
```
""" Reduce.Algebra.impart

@doc """
    nextprime(r)

`nextprime` returns the next prime greater than its integer argument, using a probabilistic algorithm. A type error occurs if the value of the argument is not an integer. For example:

```Julia
julia> Algebra.nextprime(5)
7

julia> Algebra.nextprime(-2)
2

julia> Algebra.nextprime(-7)
-5

julia> Algebra.nextprime(1000000)
1000003
```
whereas `Algebra.nextprime(:a)` gives a type error.
""" Reduce.Algebra.nextprime

@doc """
    random(r)

`random(n)` returns a random number ``r`` in the range ``0 ≤ r < n``. A type error occurs if the value of the argument is not a positive integer in algebraic mode, or positive number in symbolic mode. For example:

```Julia
julia> Algebra.random(5)
3

julia> Algebra.random(1000)
191
```
whereas `Algebra.random(:a)` gives a type error.
""" Reduce.Algebra.random

@doc """
    random_new_seed(r)

`random_new_seed(n)` reseeds the random number generator to a sequence determined by the integer argument `n`. It can be used to ensure that a repeatable pseudo-random sequence will be delivered regardless of any previous use of `random`, or can be called early in a run with an argument derived from something variable (such as the time of day) to arrange that different runs of a REDUCE program will use different random sequences. When a fresh copy of REDUCE is first created it is as if `random_new_seed(1)` has been obeyed.

A type error occurs if the value of the argument is not a positive integer.
""" Reduce.Algebra.random_new_seed

@doc """
    repart(r)

This returns the real part of an expression, if that argument has an numerical value. A non-numerical argument is returned as an expression in the operators `repart` and `impart`. For example:

```Julia
julia> Algebra.repart(1+im)
1

julia> Algebra.repart(:(a+im*b))
:(repart(a) - impart(b))
```
""" Reduce.Algebra.repart

@doc """
    round(r)

This operator returns the rounded value (i.e, the nearest integer) of its single argument if that argument has a numerical value. A non-numeric argument is returned as an expression in the original operator. For example:

```Julia
julia> Algebra.round(-5/4)
-1.0

julia> Algebra.round(:a)
:(round(a))
```
""" Reduce.Algebra.round

@doc """
    sign(r)

`sign` tries to evaluate the sign of its argument. If this is possible `sign` returns one of `1`, `0` or `-1`. Otherwise, the result is the original form or a simplified variant. For example:

```Julia
julia> Algebra.sign(-5)
-1

julia> Algebra.sign(:(-a^2*b))
:(-(sign(b)))
```
Note that even powers of formal expressions are assumed to be positive only as long as the switch `complex` is off.
""" Reduce.Algebra.sign

@doc """
    bernoulli(n)

The unary operator `bernoulli` provides notation and computation for Bernoulli numbers. `bernoulli(n)` evaluates to the `n`th Bernoulli number; all of the odd Bernoulli numbers, except `bernoulli(1)`, are zero.

The algorithms are based upon those by Herbert Wilf, presented by Sandra Fillebrown [?]. If the `rounded` switch is off, the algorithms are exactly those; if it is on, some further rounding may be done to prevent computation of redundant digits. Hence, these functions are particularly fast when used to approximate the Bernoulli numbers in rounded mode.
""" Reduce.Algebra.bernoulli

@doc """
    fibonacci(n)

The unary operator `fibonacci` provides notation and computation for Fibonacci numbers. `fibonacci(n)` evaluates to the `n`th Fibonacci number. If `n` is a positive or negative integer, it will be evaluated following the definition:

\$F_0 = 0; F_1 = 1; F_n = F_{n-1} + F_{n-2}\$
""" Reduce.Algebra.fibonacci

@doc """
    motzkin(n)

A Motzkin number \$M_n\$ (named after Theodore Motzkin) is the number of different ways of drawing non-intersecting chords on a circle between n points. For a non-negative integer `n`, the operator `motzkin(n)` returns the `n`th Motzkin number, according to the recursion formula

\$M_0 =  1;  M_1 = 1;  M_{n+1}  =  \\frac{2n+3}{n+3}M_n + \\frac{3n}{n+3}M_{n-1}.\$
""" Reduce.Algebra.motzkin

@doc """
    saveas(expr)

If the user wishes to assign the workspace to a variable or expression for later use, the `saveas` statement can be used. It has the syntax
```Julia
R"saveas ⟨expression⟩"
```
For example, after the differentiation in the last example, the workspace holds the expression `2*x+2*y`. If we wish to assign this to the variable `z` we can now say
```Julia
julia> Algebra.saveas(:z)
```
If the user wishes to save the expression in a form that allows him to use some of its variables as arbitrary parameters, the `for all` command can be used.

*Example:*
```Julia
R"for all x saveas h(x)"
```
with the above expression would mean that `h(z)` evaluates to `2*y+2*z`.
""" Reduce.Algebra.saveas

@doc """
    decompose(p)

The `decompose` operator takes a multivariate polynomial as argument, and returns an expression and a list of equations from which the original polynomial can be found by composition. Its syntax is:
```Julia
R"decompose(EXPRN:polynomial)"
```
For example:
```Julia
julia> Algebra.decompose(:(x^8-88*x^7+2924*x^6-43912*x^5+263431*x^4-218900*x^3+65690*x^2-7700*x+234))
(:(u ^ 2 + 35u + 234), :(u = v ^ 2 + 10v), :(v = x ^ 2 - 22x))

julia> Algebra.decompose(:(u^2+v^2+2u*v+1))
(:(w ^ 2 + 1), :(w = u + v))
```
Users should note however that, unlike factorization, this decomposition is not unique.
""" Reduce.Algebra.decompose

den(r::T) where T <: Number = den(RExpr(r)) |> parse

@doc """
    den(r)

This is used with the syntax:
```Julia
R"den(EXPRN:rational)"
```
It returns the denominator of the rational expression `EXPRN`. If `EXPRN` is a polynomial, 1 is returned.

*Examples:*
```Julia
julia> Algebra.den(:(x/y^2))
:(y ^ 2)

julia> Algebra.den(100//6)
3               [since 100/6 is first simplified to 50/3]

julia> Algebra.den(:(a/4+b/6))
12

julia> Algebra.den(:(a+b))
1
```
""" Reduce.Algebra.den

mainvar(r::T) where T <: Number = mainvar(RExpr(r)) |> parse

@doc """
    mainvar(exprn)

Syntax:
```Julia
R"mainvar(EXPRN:polynomial)"
```
Returns the main variable (based on the internal polynomial representation) of `EXPRN`. If `EXPRN` is a domain element, 0 is returned.

*Examples:* Assuming `a` has higher kernel order than `b`, `c`, or `d`:
```Julia
julia> Algebra.on(:exp)

julia> Algebra.mainvar(:((a+b)*(c+2*d)^2))
:a

julia> Algebra.mainvar(2)
0
```
""" Reduce.Algebra.mainvar

num(r::T) where T <: Number = num(RExpr(r)) |> parse

@doc """
    num(exprn)

Syntax:
```Julia
R"num(EXPRN:rational)"
```
Returns the numerator of the rational expression `EXPRN`. If `EXPRN` is a polynomial, that polynomial is returned.

*Examples:*
```Julia
julia> Algebra.num(:(x/y^2))
:x

julia> Algebra.num(100//6)
50

julia> Algebra.num(:(a/4+b/6))
:(3a + 2b)

julia> Algebra.num(:(a+b))
:(a + b)
```
""" Reduce.Algebra.num

setmod(r::Integer) = setmod(RExpr(r)) |> parse

@doc """
    setmod(::Integer)

REDUCE includes facilities for manipulating polynomials whose coefficients are computed modulo a given base. To use this option, two commands must be used; `R"setmod ⟨integer⟩"`, to set the prime modulus, and `on(:modular)` to cause the actual modular calculations to occur. For example, with `R"setmod 3"` and `R"on modular"`, the polynomial `(a+2*b)^3` would become `a^3+2*n^3`.

The argument of `setmod` is evaluated algebraically, except that non-modular (integer) arithmetic is used. Thus the sequence
```Julia
R"setmod 3; on modular; setmod 7"
```
will correctly set the modulus to 7.
""" Reduce.Algebra.setmod

@doc """
    root_val(exprn)

The `root_val` operator takes a single univariate polynomial as argument, and returns a list of root values at system precision (or greater if required to separate roots). It is used with the syntax
```Julia
R"root_val(EXPRN:univariate polynomial)"
```
For example, the sequence
```Julia
reduce> on rounded; root_val(x^3-x-1);
```
gives the result
```
{0.562279512062*I - 0.662358978622, - 0.562279512062*I

  - 0.662358978622,1.32471795724}
```
""" Reduce.Algebra.root_val

@doc """
    clearrules(r)

`clearrules` has the syntax
```Julia
R"clearrules <rule list>|<name of rule list>(,...)"
```
""" Reduce.Algebra.clearrules

@doc """
    showrules(r)

The operator `showrules` takes a single identifier as argument, and returns in rule-list form the operator rules associated with that argument. For example:
```Julia
reduce> showrules log;

{log(e) => 1,

 log(1) => 0,

      ~x
 log(e  ) => ~x,

                    1
 df(log(~x),~x) => ----}
                    ~x
```
Such rules can then be manipulated further as with any list. For example `R"rhs first ws"` has the value `1`. Note that an operator may have other properties that cannot be displayed in such a form, such as the fact it is an odd function, or has a definition defined as a procedure.
""" Reduce.Algebra.showrules

@doc """
    det(exprn)

Syntax:
```
        det(EXPRN:matrix_expression):algebraic.
```
The operator `det` is used to represent the determinant of a square matrix expression. E.g.,
```
Algebra.det(:(y^2))
```
is a scalar expression whose value is the determinant of the square of the matrix `y`, and
```
Algebra.det([:a :b :c; :d :e :f; :g :h :j])
```
is a scalar expression whose value is the determinant of the matrix
```Julia
3×3 Array{Symbol,2}:
 :a  :b  :c
 :d  :e  :f
 :g  :h  :j
```
Determinant expressions have the instant evaluation property. In other words, the statement
```
        let det mat((a,b),(c,d)) = 2;
```
sets the value of the determinant to `2`, and does not set up a rule for the determinant itself.
""" Reduce.Algebra.det

@doc """
    trace(exprn)

Syntax:
```
        TRACE(EXPRN:matrix_expression):algebraic.
```
The operator TRACE is used to represent the trace of a square matrix.
""" Reduce.Algebra.trace

@doc """
    nullspace(exprn)

Syntax:
```
        NULLSPACE(EXPRN:matrix_expression):list
```
`nullspace` calculates for a matrix `a` a list of linear independent vectors (a basis) whose linear combinations satisfy the equation \$Ax = 0\$. The basis is provided in a form such that as many upper components as possible are isolated.

Note that with `b := nullspace a` the expression `length b` is the *nullity* of `a`, and that `second length a - length b` calculates the *rank* of `a`. The rank of a matrix expression can also be found more directly by the `rank` operator described below.

*Example:* The command
```
        nullspace mat((1,2,3,4),(5,6,7,8));
```
gives the output
```
        {
         [ 1  ]
         [ 0  ]
         [ - 3]
         [ 2  ]
         ,
         [ 0  ]
         [ 1  ]
         [ - 2]
         [ 1  ]
         }
```
In addition to the REDUCE matrix form, `nullspace` accepts as input a matrix given as a list of lists, that is interpreted as a row matrix. If that form of input is chosen, the vectors in the result will be represented by lists as well. This additional input syntax facilitates the use of `nullspace` in applications different from classical linear algebra.
""" Reduce.Algebra.nullspace

@doc """
    rank(exprn)

Syntax:
```
        RANK(EXPRN:matrix_expression):integer
```
`rank` calculates the rank of its argument, that, like `nullspace` can either be a standard matrix expression, or a list of lists, that can be interpreted either as a row matrix or a set of equations.

*Example:*
```Julia
Algebra.rank([:a :b :c; :d :e :f])
```
returns the value `2`.
""" Reduce.Algebra.rank