text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module BellNumbers
using Nemo, Triangles, Products, NumberTheory, SeqUtils

export ModuleBellNumbers
export BellTrans, BellTriangle, BellNumberList, BellNumber
export V000110, L000110
export T137452, T264428, T137513, T104556, T001497, T132062, T039683, T203412
export T004747, T051141, T265606, T119274, T000369, T051142

"""

The Bell transform transforms an integer sequence into an integer triangle; also known as incomplete Bell polynomials. Let ``X`` be an integer sequence, then

 ``B_{n, k}(X) = \\sum_{m=1}^{n-k+1} \\binomial{n-1}{m-1} X[m] B_{n-m,k-1}(X)``

 where ``B_{0,0} = 1, B_{n,0} = 0`` for ``n≥1, B_{0,k} = 0`` for ``k≥1``.

The Bell transform is (0,0)-based and the associated triangle always has as first column 1,0,0,0,... This column is often missing in the OEIS. Other Stirling number related sequences are implemented in the module StirlingLahNumbers.

* BellTrans, BellTriangle, BellNumberList, BellNumber
* V000110, L000110, T137452, T264428, T137513, T104556, T001497, T132062, T039683, T203412, T004747, T051141, T265606, T119274, T000369, T051142
"""
const ModuleBellNumbers = ""

"""

Return a list of the first m Bell numbers (a.k.a. exponential numbers).
"""
function BellNumberList(m::Int)
    m == 0 && return fmpz[]
    R = ZArray(m)
    R[1] = 1
    m == 1 && return R
    R[2] = 1
    m == 2 && return R

    A = ZArray(m)
    A[1] = fmpz(1)
    for n ∈ 2:m-1
        A[n] = A[1]
        for k ∈ n:-1:2
            A[k-1] += A[k]
        end
        R[n+1] = A[1]
    end
    R
end

# --- much slower:
#function BellNumberList(m::Int)
#    m == 0 && return fmpz[]
#    SeqArray(m, n -> Nemo.bell(n))
#end

"""

Return the n-th Bell number. Bell numbers count the ways to partition a set of ``n`` labeled elements.

```
julia> BellNumber(10)
115975
```
"""
BellNumber(n::Int) = Nemo.bell(n)

"""

Return the n-th Bell number ``B_n``.

```
julia> V000110(11)
678570
```
"""
V000110(n::Int) = Nemo.bell(n)

"""

Return a list of Bell numbers of length len.

```
julia> L000110(10)
[1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147]
```
"""
L000110(len::Int) = BellNumberList(len)

"""

The Bell transform transforms an integer sequence into an integer triangle; also known as incomplete Bell polynomials. Let ``X`` be an integer sequence, then
``B_{n,k}(X) = \\sum_{m=1}^{n-k+1} \\binomial{n-1}{m-1} X[m] B_{n-m,k-1}(X)``
where ``B_{0,0} = 1, B_{n,0} = 0`` for ``n≥1, B_{0,k} = 0`` for ``k≥1``.
"""
function BellTrans(n::Int, k::Int, X::Array)
    global CacheBellA
    haskey(CacheBellA, (n, k, X)) && return CacheBellA[(n, k, X)]

    a, s = fmpz(1), fmpz(0)
    ((n == 0) && (k == 0)) && return a
    ((n == 0) || (k == 0)) && return s

    for m ∈ 1:n-k+1
        s += a * BellTrans(n - m, k - 1, X) * X[m]
        a = div(a * (n - m), m)
    end

    CacheBellA[(n, k, X)] = s
end

const CacheBellA = Dict{Tuple{Int,Int,Array},fmpz}()

"""

The Bell transform transforms an integer sequence into an integer triangle; also known as incomplete Bell polynomials. Let ``F`` be an integer sequence generating function, then

``B_{n,k}(F) = \\sum_{m=1}^{n-k+1} \\binomial{n-1}{m-1} F(m) B_{n-m,k-1}(F)``

where ``B_{0,0} = 1, B_{n,0} = 0`` for ``n≥1, B_{0,k} = 0`` for ``k≥1``.
"""
function BellTrans(n::Int, k::Int, F::Function)
    global CacheBellF
    haskey(CacheBellF, (n, k, F)) && return CacheBellF[(n, k, F)]

    a, s = fmpz(1), fmpz(0)
    ((n == 0) && (k == 0)) && return a
    ((n == 0) || (k == 0)) && return s

    for m ∈ 1:n-k+1
        s += a * BellTrans(n - m, k - 1, F) * F(m - 1)
        a = div(a * (n - m), m)
    end

    CacheBellF[(n, k, F)] = s
end

const CacheBellF = Dict{Tuple{Int,Int,Function},fmpz}()

"""

The Bell triangle gathers the results of the Bell transform applied to the initial segments of the input sequence.

Famously the sequence (1,1,1,...) is mapped to the triangle of the Stirling set numbers.

```
julia> ShowAsΔ(BellTriangle(5, k -> 1))
1
0 1
0 1 1
0 1 3 1
0 1 7 6 1
```
"""
function BellTriangle(n::Int, seq)
    M = ZTriangle(n)
    i = 1

    for j ∈ 0:n-1, k ∈ 0:j
        M[i] = BellTrans(j, k, seq)
        i += 1
    end
    return M
end

"""

Return the coefficients of the first ``n`` Abel polynomials.

```
julia> ShowAsΔ(T137452(5))
1
0 1
0 -2 1
0 9 -6 1
0 -64 48 -12 1
```
"""
T137452(n::Int) = BellTriangle(n, k -> (-k - 1)^k)

"""

Return the Bell transform of the Bell numbers.

```
julia> ShowAsΔ(T264428(5))
1
0 1
0 1 1
0 2 3 1
0 5 11 6 1
```
"""
T264428(n::Int) = BellTriangle(n, BellNumber)

# --- Keep the comments below.
# Stirling number related sequences are now implemented in the module StirlingLah.
# doc"""
# Return the triangle of Stirling set numbers (a.k.a. Stirling numbers of 2nd kind).
# """
# StirlingSetTriangle(n::Int) = BellTriangle(n, k -> 1)
# T048993(n::Int) = StirlingSetTriangle(n)

"""

Return the triangle of the coefficients of the Mittag-Leffler polynomials.

```
julia> ShowAsΔ(T137513(5))
1
0 2
0 0 4
0 4 0 8
0 0 32 0 16
```
"""
T137513(n::Int) = BellTriangle(n, k -> isOdd(k) ? 0 : 2*factorial(k))

"""

Return the matrix inverse of coefficients of Bessel polynomials; essentially the same as coefficients of modified Hermite polynomials T096713.

```
julia> ShowAsΔ(T104556(5))
1
0 1
0 -1 1
0 0 -3 1
0 0 3 -6 1
```
"""
T104556(n::Int) = BellTriangle(n, k -> k < 2 ? (-1)^k : 0)

"""

Return a triangle of coefficients of Bessel polynomials (better use A132062).
"""
T001497(n::Int) = BellTriangle(n, MultiFactorial(2, 1))

# Stirling number related sequences are now implemented in the module StirlingLah.
# Keep the comments below.
# doc"""
# Return the triangle of unsigned Stirling numbers of the first kind (a.k.a.
# Stirling cycle numbers).
# """
# StirlingCycleTriangle(n::Int) = BellTriangle(n, MultiFactorial(1,1))
# T132393(n::Int) = StirlingCycleTriangle(n)

"""

Return the triangle of coefficients of Bessel polynomials, also the Sheffer triangle ``(1, 1 - √(1 - 2x))`` (Cf. A001497).

```
julia> ShowAsΔ(T132062(5))
1
0 1
0 1 1
0 3 3 1
0 15 15 6 1
```
"""
T132062(n::Int) = BellTriangle(n, MultiFactorial(2, 1))

"""

Return the signed double Pochhammer triangle: expansion of ``x(x-2)(x-4)..(x-2n+2)``.

```
julia> ShowAsΔ(T039683(5))
1
0 1
0 2 1
0 8 6 1
0 48 44 12 1
```
"""
T039683(n::Int) = BellTriangle(n, MultiFactorial(2, 2))

"""

Return the Bell transform of the MultiFactorial numbers of type (3,1).

```
julia> ShowAsΔ(T203412(5))
1
0 1
0 1 1
0 4 3 1
0 28 19 6 1
```
"""
T203412(n::Int) = BellTriangle(n, MultiFactorial(3, 1))

"""

Return the Bell transform of the MultiFactorial numbers of type (3,2).

```
julia> ShowAsΔ(T004747(5))
1
0 1
0 2 1
0 10 6 1
0 80 52 12 1
```
"""
T004747(n::Int) = BellTriangle(n, MultiFactorial(3, 2))

"""

Return the triangle ``3^{n-m}S1(n, m)`` where S1 are the signed Stirling numbers of first kind.

```
julia> ShowAsΔ(T051141(5))
1
0 1
0 3 1
0 18 9 1
0 162 99 18 1
```
"""
T051141(n::Int) = BellTriangle(n, MultiFactorial(3, 3))

"""

Return the Bell transform of the quartic factorial numbers.

```
julia> ShowAsΔ(T265606(5))
1
0 1
0 1 1
0 5 3 1
0 45 23 6 1
```
"""
T265606(n::Int) = BellTriangle(n, MultiFactorial(4, 1))

"""

Return the triangle of coefficients of numerators in Pade approximation to ``e^x``.

```
julia> ShowAsΔ(T119274(5))
1
0 1
0 2 1
0 12 6 1
0 120 60 12 1
```
"""
T119274(n::Int) = BellTriangle(n, MultiFactorial(4, 2))

"""

Return the Bell transform of the MultiFactorial numbers of type (4,3).

```
julia> ShowAsΔ(T000369(5))
1
0 1
0 3 1
0 21 9 1
0 231 111 18 1
```
"""
T000369(n::Int) = BellTriangle(n, MultiFactorial(4, 3))

"""

Return the Bell transform of the MultiFactorial numbers of type (4,4).

```
julia> ShowAsΔ(T051142(5))
1
0 1
0 4 1
0 32 12 1
0 384 176 24 1
```
"""
T051142(n::Int) = BellTriangle(n, MultiFactorial(4, 4))

#START-TEST-########################################################

using Test, SeqTests

function test()

    @testset "Bell" begin
        seq = fmpz[1, 1, 2, 5, 15, 52]
        a = [BellTrans(6, k, seq) for k ∈ 0:5]
        b = [0, 52, 205, 210, 85, 15]
        @test all(a .== b)

        seq = [1, 1, 1, 1, 1, 1]
        a = BellTriangle(6, seq)
        b = fmpz[1, 0, 1, 0, 1, 1, 0, 1, 3, 1, 0, 1, 7, 6, 1, 0, 1, 15, 25]
        @test all(a[1:5] .== b[1:5])

        a = fmpz[1, 1, 2, 5, 15, 52, 203]
        b = [BellNumber(n) for n ∈ 0:6]
        @test all(a .== b)

        if data_installed()
            SeqTest(V000110, 'V')
            SeqTest([L000110], 'L')

            T = [T264428, T137452, T132062, T265606]
            SeqTest(T, 'T')

            # These triangles have an additional column 1,0,0,... on the left
            # or have different signs. But they are essentially the same.
            P = [T104556, T001497, T039683, T203412, T004747,
                 T051141, T119274, T000369, T051142, T137513]
            SeqTest(P, 'P')
        end
    end
end

function demo()
    seq = [1, 1, 2, 5, 15, 52]
    len = size(seq)[1]

    for n ∈ 0:len
        println([BellTrans(n, k, seq) for k ∈ 0:n])
    end

    M = BellTriangle(len, seq)
    ShowAsΔ(M)

    for n ∈ 0:6
        println(n, " -> ", BellNumberList(n))
    end

    T = T137452(8)
    ShowAsΔ(T)
end

"""

BellNumberList(1000) :: 0.539088 seconds (855.00 k allocations: 13.062 MiB)
BellTriangle(100, [1 for _ ∈ 1:100]) :: 0.576952 seconds (1.22 M allocations: 24.571 MB, 14.23% gc time)
"""
function perf()
    GC.gc()
    BellNumberList(5)
    @time BellNumberList(1000)
    @time BellTriangle(100, [1 for _ ∈ 1:100])
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
_strip_module(s) = split(s, '.', limit=2)[end]
_strip_type_param(s) = replace(s, r"{.+?}" => "")
_clean_name(d::Distributions.Distribution) = _strip_module(_strip_type_param(repr(d)))

function scatter(d::Distributions.ContinuousUnivariateDistribution)
    ls(a, b, c) = range(a, stop=b, length=c)
    x = ls(quantile.([d], [0.01, 0.99])..., 100)
    trace = scatter(x=x, y=pdf.([d], x), name=_clean_name(d))
end

Plot(d::Distributions.UnivariateDistribution...) = Plot(collect(map(scatter, d)))

ext: jl


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8]);F = [261*x1^4+324*x1^3*x6+744*x1^3*x7+9604*x1^2*x2^2+14896*x1^2*x2*x3-1092*x1^2*x2*x6+5776*x1^2*x3^2-1680*x1^2*x3*x8-474*x1^2*x4*x6+1296*x1^2*x5^2+729*x1^2*x6^2+3348*x1^2*x6*x7-612*x1^2*x6*x8+5693*x1^2*x7^2-10780*x1*x2*x3^2-20604*x1*x2*x4*x7-2992*x1*x2*x5*x7-4914*x1*x2*x6^2-11284*x1*x2*x6*x7-8360*x1*x3^3-9576*x1*x3*x4*x7+3040*x1*x3*x5*x7-6120*x1*x4^2*x5+1512*x1*x4*x6^2+3472*x1*x4*x6*x7-86*x1*x4*x7*x8+6552*x1*x5*x7^2-2754*x1*x6^2*x8-6324*x1*x6*x7*x8+11520*x2^2*x4^2+8281*x2^2*x6^2+9745*x2^2*x7^2+138*x2*x3*x4*x7-3936*x2*x3*x4*x8+23250*x2*x4^2*x7+192*x2*x4^2*x8-5096*x2*x4*x6^2-4416*x2*x4*x6*x7-3128*x2*x4*x7*x8+9282*x2*x6^2*x8-2484*x2*x6*x7^2-17472*x2*x7^3+3025*x3^4+9*x3^2*x4^2+6930*x3^2*x4*x7-2200*x3^2*x5*x7+4817*x3^2*x8^2+378*x3*x4^3-408*x3*x4^2*x8-324*x3*x4*x6*x7+3024*x3*x4*x6*x8-3444*x3*x4*x7*x8+3772*x3*x6*x7*x8+11194*x4^4-8568*x4^3*x8+1513*x4^2*x6^2-6804*x4^2*x6*x7-9737*x4^2*x7^2+4625*x4^2*x8^2-2520*x4*x5*x7^2-2856*x4*x6^2*x8-3864*x4*x6*x7^2+7344*x4*x6*x7*x8+400*x5^2*x7^2+5032*x6^2*x7^2+2601*x6^2*x8^2+8281*x7^4-120*x1^2*x6-18032*x1*x2*x3-13984*x1*x3^2-5074*x1*x5*x7-7912*x1*x7*x8+11328*x2*x4*x5-3648*x2*x4*x7+17664*x2*x4*x8-4508*x2*x5*x7+10120*x3^3-588*x3*x4*x5+11592*x3*x4*x7-3680*x3*x5*x7+448*x3*x6*x8+3116*x3*x7*x8-12348*x4^2*x5+13446*x4*x5*x8+216*x4*x6^2-3192*x4*x7^2+184*x4*x8^2+10584*x5*x6*x7+3496*x6*x7^2-1260*x1^2-3168*x1*x5-6880*x1*x7+10368*x2*x4+8448*x2*x7+8464*x3^2+8968*x3*x8+7480*x4^2+2268*x4*x6-4368*x4*x7+160*x4*x8+13085*x5^2+10856*x5*x8+16*x6^2+4784*x6*x7-6564*x7^2+8464*x8^2+9440*x5+336*x6+3952*x7+14720*x8+12804, 19208*x1^2*x2+14896*x1^2*x3-1092*x1^2*x6-10780*x1*x3^2-20604*x1*x4*x7-2992*x1*x5*x7-4914*x1*x6^2-11284*x1*x6*x7+23040*x2*x4^2+16562*x2*x6^2+19490*x2*x7^2+138*x3*x4*x7-3936*x3*x4*x8+23250*x4^2*x7+192*x4^2*x8-5096*x4*x6^2-4416*x4*x6*x7-3128*x4*x7*x8+9282*x6^2*x8-2484*x6*x7^2-17472*x7^3-18032*x1*x3+11328*x4*x5-3648*x4*x7+17664*x4*x8-4508*x5*x7+10368*x4+8448*x7, 14896*x1^2*x2+11552*x1^2*x3-1680*x1^2*x8-21560*x1*x2*x3-25080*x1*x3^2-9576*x1*x4*x7+3040*x1*x5*x7+138*x2*x4*x7-3936*x2*x4*x8+12100*x3^3+18*x3*x4^2+13860*x3*x4*x7-4400*x3*x5*x7+9634*x3*x8^2+378*x4^3-408*x4^2*x8-324*x4*x6*x7+3024*x4*x6*x8-3444*x4*x7*x8+3772*x6*x7*x8-18032*x1*x2-27968*x1*x3+30360*x3^2-588*x4*x5+11592*x4*x7-3680*x5*x7+448*x6*x8+3116*x7*x8+16928*x3+8968*x8, -474*x1^2*x6-20604*x1*x2*x7-9576*x1*x3*x7-12240*x1*x4*x5+1512*x1*x6^2+3472*x1*x6*x7-86*x1*x7*x8+23040*x2^2*x4+138*x2*x3*x7-3936*x2*x3*x8+46500*x2*x4*x7+384*x2*x4*x8-5096*x2*x6^2-4416*x2*x6*x7-3128*x2*x7*x8+18*x3^2*x4+6930*x3^2*x7+1134*x3*x4^2-816*x3*x4*x8-324*x3*x6*x7+3024*x3*x6*x8-3444*x3*x7*x8+44776*x4^3-25704*x4^2*x8+3026*x4*x6^2-13608*x4*x6*x7-19474*x4*x7^2+9250*x4*x8^2-2520*x5*x7^2-2856*x6^2*x8-3864*x6*x7^2+7344*x6*x7*x8+11328*x2*x5-3648*x2*x7+17664*x2*x8-588*x3*x5+11592*x3*x7-24696*x4*x5+13446*x5*x8+216*x6^2-3192*x7^2+184*x8^2+10368*x2+14960*x4+2268*x6-4368*x7+160*x8, 2592*x1^2*x5-2992*x1*x2*x7+3040*x1*x3*x7-6120*x1*x4^2+6552*x1*x7^2-2200*x3^2*x7-2520*x4*x7^2+800*x5*x7^2-5074*x1*x7+11328*x2*x4-4508*x2*x7-588*x3*x4-3680*x3*x7-12348*x4^2+13446*x4*x8+10584*x6*x7-3168*x1+26170*x5+10856*x8+9440, 324*x1^3-1092*x1^2*x2-474*x1^2*x4+1458*x1^2*x6+3348*x1^2*x7-612*x1^2*x8-9828*x1*x2*x6-11284*x1*x2*x7+3024*x1*x4*x6+3472*x1*x4*x7-5508*x1*x6*x8-6324*x1*x7*x8+16562*x2^2*x6-10192*x2*x4*x6-4416*x2*x4*x7+18564*x2*x6*x8-2484*x2*x7^2-324*x3*x4*x7+3024*x3*x4*x8+3772*x3*x7*x8+3026*x4^2*x6-6804*x4^2*x7-5712*x4*x6*x8-3864*x4*x7^2+7344*x4*x7*x8+10064*x6*x7^2+5202*x6*x8^2-120*x1^2+448*x3*x8+432*x4*x6+10584*x5*x7+3496*x7^2+2268*x4+32*x6+4784*x7+336, 744*x1^3+3348*x1^2*x6+11386*x1^2*x7-20604*x1*x2*x4-2992*x1*x2*x5-11284*x1*x2*x6-9576*x1*x3*x4+3040*x1*x3*x5+3472*x1*x4*x6-86*x1*x4*x8+13104*x1*x5*x7-6324*x1*x6*x8+19490*x2^2*x7+138*x2*x3*x4+23250*x2*x4^2-4416*x2*x4*x6-3128*x2*x4*x8-4968*x2*x6*x7-52416*x2*x7^2+6930*x3^2*x4-2200*x3^2*x5-324*x3*x4*x6-3444*x3*x4*x8+3772*x3*x6*x8-6804*x4^2*x6-19474*x4^2*x7-5040*x4*x5*x7-7728*x4*x6*x7+7344*x4*x6*x8+800*x5^2*x7+10064*x6^2*x7+33124*x7^3-5074*x1*x5-7912*x1*x8-3648*x2*x4-4508*x2*x5+11592*x3*x4-3680*x3*x5+3116*x3*x8-6384*x4*x7+10584*x5*x6+6992*x6*x7-6880*x1+8448*x2-4368*x4+4784*x6-13128*x7+3952, -1680*x1^2*x3-612*x1^2*x6-86*x1*x4*x7-2754*x1*x6^2-6324*x1*x6*x7-3936*x2*x3*x4+192*x2*x4^2-3128*x2*x4*x7+9282*x2*x6^2+9634*x3^2*x8-408*x3*x4^2+3024*x3*x4*x6-3444*x3*x4*x7+3772*x3*x6*x7-8568*x4^3+9250*x4^2*x8-2856*x4*x6^2+7344*x4*x6*x7+5202*x6^2*x8-7912*x1*x7+17664*x2*x4+448*x3*x6+3116*x3*x7+13446*x4*x5+368*x4*x8+8968*x3+160*x4+10856*x5+16928*x8+14720];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
module juliaPackage

"Function returns minumum of a function with Powell's interpolation method"

function powell(f, x0; step=1.0, ϵ=1e-5)

    c = x0 + step

    if f(x0)<f(c)
        a = x0 - step
    else
        a = x0 + 2*step
        a, c = c, a
    end

    b = x0

    i = 1
    maxIter = 10^6

    while i < maxIter

        first = (1/2)
        second = (b^2-c^2) * f(a) + (c^2-a^2) * f(b) + (a^2-b^2)f(c)
        third = (b-c) * f(a) + (c-a) * f(b) + (a-b) * f(c)

        d = first*second/third

        first = (b-c) * f(a) + (c-a) * f(b) + (a-b) * f(c)
        second = (a-b) * (b-c) * (c-a)

        if first/second < 0

            if abs(a-d)<ϵ
                return (a, f(a))
            end

            if abs(b-d)<ϵ
                return (b, f(b))
            end

            if abs(c-d)<ϵ
                return (c, f(c))
            end

            list = sort([(a,f(a)), (b,f(b)), (c,f(c))], by = x -> x[2])

            if f(d) > f(a)
                far = sort([(a,abs(a-d)), (b,abs(b-d)), (c,abs(c-d))], by = x -> x[2])
                final = sort([far[1][1], far[2][1], d])
                a = final[1]
                b = final[2]
                c = final[3]
            else
                a = d
                b = list[1][1]
                c = list[2][1]
            end
        end
    end

    throw("Couldn't find minimum")

end

export powell

end # module

max_stars_repo_name: Haloxxx/juliaPackage

max_issues_repo_name: Haloxxx/juliaPackage

max_forks_repo_name: Haloxxx/juliaPackage


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

name = "boost"
version = v"1.76.0"

# Collection of sources required to build boost
sources = [
    ArchiveSource("https://boostorg.jfrog.io/artifactory/main/release/1.76.0/source/boost_$(version.major)_$(version.minor)_$(version.patch).tar.bz2",
                  "f0397ba6e982c4450f27bf32a2a83292aba035b827a5623a14636ea583318c41"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/boost*/

./bootstrap.sh --prefix=$prefix --without-libraries=python --with-toolset="--cxx=${CXX_FOR_BUILD}"

# Patch adapted from
# https://svnweb.freebsd.org/ports/head/devel/boost-libs/files/patch-boost_math_tools_config.hpp?revision=439932&view=markup
# to be able to build long double math libraries
atomic_patch -p1 ../patches/boost_math_tools_config_hpp.patch

rm project-config.jam
toolset=gcc
targetos=linux
extraargs=

# BinaryBuilderBase compiler wrappers don't like it when we use -march=ANYTHING.
# So we patch that out. However, we cannot just insert an empty string; so we instead
# add another "harmless" option; we choose `-Wall` which is already passed anyway
sed -i "s/-march=i686/-Wall/g" tools/build/src/tools/gcc.*

if [[ $target == *apple* ]]; then
    targetos=darwin
    toolset=darwin-6.0
    extraargs="binary-format=mach-o link=shared"
    echo "using darwin : 6.0 : $CXX : <cxxflags>-stdlib=libc++ <linkflags>-stdlib=libc++ ;" > project-config.jam
    if [[ "${target}" == aarch64-* ]]; then
        # Fix error
        #     Undefined symbols for architecture arm64:
        #       "_jump_fcontext", referenced from:
        # See https://github.com/boostorg/context/issues/170#issuecomment-863669877
        extraargs="abi=aapcs ${extraargs}"
    fi
elif [[ $target == x86_64*mingw* ]]; then
    targetos=windows
    extraargs="address-model=64 binary-format=pe abi=ms link=shared"
elif [[ $target == i686*mingw* ]]; then
    targetos=windows
    extraargs="address-model=32 binary-format=pe abi=ms link=shared"
elif [[ $target == *freebsd* ]]; then
    targetos=freebsd
    toolset=clang-6.0
    extraargs="address-model=64 link=shared"
    echo "using clang : 6.0 : $CXX : <linkflags>\\"$LDFLAGS\\" ;" > project-config.jam
fi
./b2 -j${nproc} toolset=$toolset target-os=$targetos $extraargs variant=release --prefix=$prefix --without-python --layout=system --debug-configuration install

install_license LICENSE_1_0.txt
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libboost_atomic", :libboost_atomic),
    LibraryProduct("libboost_chrono", :libboost_chrono),
    LibraryProduct("libboost_container", :libboost_container),
    LibraryProduct("libboost_context", :libboost_context),
    LibraryProduct("libboost_contract", :libboost_contract),
    LibraryProduct("libboost_coroutine", :libboost_coroutine),
    LibraryProduct("libboost_date_time", :libboost_date_time),
    LibraryProduct("libboost_filesystem", :libboost_filesystem),
    LibraryProduct("libboost_graph", :libboost_graph),
    LibraryProduct("libboost_iostreams", :libboost_iostreams),
    # boost_locale segfaults on windows, see https://github.com/benlorenz/boostBuilder/issues/2
    #LibraryProduct("libboost_locale", :libboost_locale),
    LibraryProduct("libboost_log", :libboost_log),
    LibraryProduct("libboost_log_setup", :libboost_log_setup),
    LibraryProduct("libboost_prg_exec_monitor", :libboost_prg_exec_monitor),
    LibraryProduct("libboost_program_options", :libboost_program_options),
    LibraryProduct("libboost_random", :libboost_random),
    LibraryProduct("libboost_regex", :libboost_regex),
    LibraryProduct("libboost_serialization", :libboost_serialization),
    LibraryProduct("libboost_stacktrace_basic", :libboost_stacktrace_basic),
    LibraryProduct("libboost_stacktrace_noop", :libboost_stacktrace_noop),
    LibraryProduct("libboost_system", :libboost_system),
    LibraryProduct("libboost_thread", :libboost_thread),
    LibraryProduct("libboost_timer", :libboost_timer),
    LibraryProduct("libboost_type_erasure", :libboost_type_erasure),
    LibraryProduct("libboost_unit_test_framework", :libboost_unit_test_framework),
    LibraryProduct("libboost_wave", :libboost_wave),
    LibraryProduct("libboost_wserialization", :libboost_wserialization),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7", julia_compat="1.6")

max_stars_repo_name: staticfloat/BinaryBuilderRepos

max_issues_repo_name: staticfloat/BinaryBuilderRepos

max_forks_repo_name: staticfloat/BinaryBuilderRepos


text:
struct MultilevelScheme{T<:AbstractScheme}
  scheme::T
  substeps::Int64
  levels::UnitRange{Int64}
end

function schemes(s::MultilevelScheme)
  # subdivide.(s.scheme, s.substeps.^s.levels)
  [subdivide(s.scheme, s.substeps^l) for l in s.levels]
end

function npaths(s::MultilevelScheme, max_budget)
  budget = max_budget
  nlevels = length(s.levels)
  steps = s.substeps.^s.levels
  n_total = zeros(Int64, nlevels)
  for i in nlevels:-1:1
    level_budget = minimum(floor.(Int64, budget / length(s.levels[1:i]) ./ steps[1:i]) .* steps[1:i])
    n_partial = div.(level_budget, steps[1:i])
    n_total[1:i] += n_partial
    budget -= dot(n_partial, steps[1:i])
  end
  n_total
end

function cost(s::MultilevelScheme, npaths)
  steps = s.substeps.^s.levels
  s.substeps.^s.levels .* npaths
end

simulate(model, scheme::MultilevelScheme, t0, x0, nsteps, npaths::Int64) =
  simulate(model, scheme, t0, x0, nsteps, fill(npaths, length(scheme.levels)))

sample(model, scheme::MultilevelScheme, t0, x0, nsteps, npaths::Int64) =
  sample(model, scheme, t0, x0, nsteps, fill(npaths, length(scheme.levels)))

function simulate(model, scheme::MultilevelScheme, t0, x0::T, nsteps, npaths::Array{Int64,1}) where {T<:Union{Float64, SVector{D,Float64} where D}}
  nlevels = length(scheme.levels)
  ml_schemes = schemes(scheme)
  ml_steps = scheme.substeps.^scheme.levels
  x1, t1 = simulate(model, ml_schemes[1], t0, x0, nsteps*ml_steps[1], npaths[1])
  x = Array{typeof(x1)}(undef, 2*nlevels-1)
  t = Array{typeof(t1)}(undef, 2*nlevels-1)
  x[1] = x1
  t[1] = t1
  for i in 2:nlevels
    x[2i-2], x[2i-1], t[2i-2], t[2i-1] = _multilevel_simulate(model, ml_schemes[i-1], ml_schemes[i], scheme.substeps, t0, x0, nsteps*ml_steps[i-1], npaths[i])
  end
  x, t
end

function sample(model, scheme::MultilevelScheme, t0, x0::T, nsteps, npaths::Array{Int64,1}) where {T<:Union{Float64, SVector{D,Float64} where D}}
  nlevels = length(scheme.levels)
  ml_schemes = schemes(scheme)
  ml_steps = scheme.substeps.^scheme.levels
  x1 = sample(model, ml_schemes[1], t0, x0, nsteps*ml_steps[1], npaths[1])
  x = Array{typeof(x1)}(undef, 2*nlevels-1)
  x[1] = x1
  for i in 2:nlevels
    x[2i-2], x[2i-1] = _multilevel_sample(model, ml_schemes[i-1], ml_schemes[i], scheme.substeps, t0, x0, nsteps*ml_steps[i-1], npaths[i])
  end
  x
end

function _multilevel_sample(model, coarse_scheme, fine_scheme, substeps, t0, x0::T, nsteps, npaths) where T
  xc = Array{T}(undef, npaths)
  xf = Array{T}(undef, npaths)
  for i in 1:npaths
    xc[i], xf[i] = _multilevel_sample(model, coarse_scheme, fine_scheme, substeps, t0, x0, nsteps)
  end
  xc, xf
end

function _multilevel_sample(model::T1, coarse_scheme, fine_scheme, substeps, t0, x0::T2, nsteps) where {T1,T2}
  xc = x0
  xf = x0
  σ = sqrt(fine_scheme.Δt)
  for n in 1:nsteps
    Δw = zero(T2)
    tc = t0 + (n - 1) * coarse_scheme.Δt
    for i in 1:substeps
      tf = tc + (i - 1) * fine_scheme.Δt
      Δw += δw = σ * _randn(T1)
      xf = step(model, fine_scheme, tf, xf, δw)
    end
    xc = step(model, coarse_scheme, tc, xc, Δw)
  end
  xc, xf
end

function _multilevel_simulate(model, coarse_scheme, fine_scheme, substeps, t0, x0::T, nsteps, npaths::Vararg{Integer,N}) where {T,N}
  tf = t0 .+ fine_scheme.Δt * (0:nsteps*substeps)
  tc = tf[1:substeps:end]
  xf = Array{T}(undef, nsteps*substeps+1, npaths...)
  wiener!(xf, fine_scheme.Δt)
  xc = xf[1:substeps:end,:]
  simulate!(xf, model, fine_scheme, t0, x0, xf)
  simulate!(xc, model, coarse_scheme, t0, x0, xc)
  xc, xf, tc, tf
end

## wrapper methods to make sure StaticArrays are used internaly ##

function sample(model, scheme::MultilevelScheme, t0, x0::Array{Float64,1}, nsteps, npaths::Array{Int64,1})
  D = length(x0)
  x = sample(model, scheme, t0, SVector{D}(x0), nsteps, npaths)
  data = [reshape(reinterpret(Float64, vec(xx)), (D, size(xx)...)) for xx in x]
end

function simulate(model, scheme::MultilevelScheme, t0, x0::Array{Float64,1}, nsteps, npaths::Array{Int64,1})
  D = length(x0)
  x, t = simulate(model, scheme, t0, SVector{D}(x0), nsteps, npaths)
  data = [reshape(reinterpret(Float64, vec(xx)), (D, size(xx)...)) for xx in x]
  data, t
end


text:
# From node.

mutable struct FromNode <: TabularNode
    source::Union{SQLTable, Symbol, Nothing}

    FromNode(; source) =
        new(source isa AbstractString ? Symbol(source) : source)
end

FromNode(source) =
    FromNode(source = source)

"""
    From(; source)
    From(source)

`From` outputs the content of a database table.  The parameter `source`
could be a [`SQLTable`](@ref) object, a symbol (for use with [`With`](@ref)),
or `nothing`.

```sql
SELECT ...
FROM \$source
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person);

julia> print(render(q))
SELECT
  "person_1"."person_id",
  "person_1"."year_of_birth"
FROM "person" AS "person_1"
```

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> condition_occurrence =
           SQLTable(:condition_occurrence, columns = [:condition_occurrence_id,
                                                      :person_id,
                                                      :condition_concept_id]);

julia> q = From(person) |>
           Where(Fun.in(Get.person_id, From(:essential_hypertension) |>
                                       Select(Get.person_id))) |>
           With(:essential_hypertension =>
                    From(condition_occurrence) |>
                    Where(Get.condition_concept_id .== 320128));

julia> print(render(q))
WITH "essential_hypertension_1" ("person_id") AS (
  SELECT "condition_occurrence_1"."person_id"
  FROM "condition_occurrence" AS "condition_occurrence_1"
  WHERE ("condition_occurrence_1"."condition_concept_id" = 320128)
)
SELECT
  "person_1"."person_id",
  "person_1"."year_of_birth"
FROM "person" AS "person_1"
WHERE ("person_1"."person_id" IN (
  SELECT "essential_hypertension_1"."person_id"
  FROM "essential_hypertension_1"
))
```

```jldoctest
julia> q = From(nothing);

julia> print(render(q))
SELECT NULL
```
"""
From(args...; kws...) =
    FromNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(From), pats::Vector{Any}) =
    dissect(scr, FromNode, pats)

Base.convert(::Type{AbstractSQLNode}, source::SQLTable) =
    FromNode(source)

function PrettyPrinting.quoteof(n::FromNode, ctx::QuoteContext)
    source = n.source
    if source isa SQLTable
        tex = get(ctx.vars, source, nothing)
        if tex === nothing
            tex = quoteof(source, limit = true)
        end
        Expr(:call, nameof(From), tex)
    elseif source isa Symbol
        Expr(:call, nameof(From), QuoteNode(source))
    else
        Expr(:call, nameof(From), source)
    end
end

function label(n::FromNode)
    source = n.source
    if source isa SQLTable
        source.name
    elseif source isa Symbol
        source
    else
        :_
    end
end


text:
using SailRoute, Test, Dates

function circular_performance()
    twa = [a for a in LinRange(40.0, 180.0, 10)]
    tws = [10.0, 15.0]
    perf = fill(5.0, 10, 2)
    perf[:, 2] .= 7.5
    polar = SailRoute.setup_perf_interpolation(tws, twa, perf)
    performance = SailRoute.Performance(polar, 1.0, 1.0, nothing)
    return performance
end

@testset "Test straight line no current" begin
   lon1 = 0.0
   lon2 = 10.0
   lat = 0.0
   dist, theta = SailRoute.haversine(lon1, lat, lon2, lat)
   wind_spd = 10.0
   wind_dir = 0.0
   curr_spd = 0.0
   curr_dir = 0.0
   bearing = SailRoute.min_angle(theta, wind_dir)
   vt = dist/SailRoute.perf_interp(circular_performance(), bearing, wind_spd)

   n_nodes = 7
   test_route = SailRoute.Route(lon1, lon2, lat, lat, n_nodes, n_nodes)

   #create times
   start_time = Dates.DateTime(2019,7,9,0,0,0)
   times = start_time:Dates.Hour(1):Dates.DateTime(2019,7,10,0,0,0)
   t_len = length(times)
   # create weather data
   wisp = fill(wind_spd, t_len, n_nodes, n_nodes)
   widi = fill(wind_dir, t_len, n_nodes, n_nodes)
   wadi = fill(0.0, t_len, n_nodes, n_nodes)
   wahi = fill(0.0, t_len, n_nodes, n_nodes)
   cusp = fill(curr_spd, t_len, n_nodes, n_nodes)
   cudi = fill(curr_dir, t_len, n_nodes, n_nodes)

   # provide environment
   x = LinRange(0.5, 9.5, n_nodes)
   y = LinRange(-5.0, 5.0, n_nodes)
   grid_x = [i for i in x, j in y]
   grid_y = [j for i in x, j in y]

   # solve the route
   rs = SailRoute.route_solve(test_route, circular_performance(), start_time, times, grid_x, grid_y, wisp, widi, wadi, wahi, cusp, cudi)
   @test rs[1] ≈ vt
end

@testset "Straight line with current" begin
   lon1 = 0.0
   lon2 = 10.0
   lat = 0.0
   dist, theta = SailRoute.haversine(lon1, lat, lon2, lat)
   wind_spd = 10.0
   wind_dir = 0.0
   curr_spd = 1.0
   curr_dir = 270.0
   v_t = dist/SailRoute.solve_speed_given_current(wind_spd, wind_dir, curr_spd, curr_dir, theta, circular_performance())
   n_nodes = 11
   test_route = SailRoute.Route(lon1, lon2, lat, lat, n_nodes, n_nodes)

   #create times
   start_time = Dates.DateTime(2019,7,9,0,0,0)
   times = start_time:Dates.Hour(1):Dates.DateTime(2019,7,10,0,0,0)
   t_len = length(times)

   # create weather data
   wisp = fill(wind_spd, t_len, n_nodes, n_nodes)
   widi = fill(wind_dir, t_len, n_nodes, n_nodes)
   wadi = fill(0.0, t_len, n_nodes, n_nodes)
   wahi = fill(0.0, t_len, n_nodes, n_nodes)
   cusp = fill(curr_spd, t_len, n_nodes, n_nodes)
   cudi = fill(curr_dir, t_len, n_nodes, n_nodes)

   # provide environment

   x = LinRange(0.5, 9.5, n_nodes)
   y = LinRange(-5.0, 5.0, n_nodes)
   grid_x = [i for i in x, j in y]
   grid_y = [j for i in x, j in y]
   # solve the route
   rs = SailRoute.route_solve(test_route, circular_performance(), start_time, times, grid_x, grid_y, wisp, widi, wadi, wahi, cusp, cudi)
   @test v_t ≈ rs[1]
end


text:
function Neuroimaging.plot(
    p::Plots.Plot,
    dip::Union{Dipole,Coordinate};
    c = :green,
    m = (8, :rect),
    l = "",
    kwargs...,
)

    Plots.scatter!(
        p.subplots[1],
        [ustrip(dip.x) * 1000],
        [ustrip(dip.y) * 1000],
        m = m,
        c = c,
        lab = l,
        legend = false;
        kwargs...,
    )
    Plots.scatter!(
        p.subplots[2],
        [ustrip(dip.y) * 1000],
        [ustrip(dip.z) * 1000],
        m = m,
        c = c,
        lab = l,
        legend = false;
        kwargs...,
    )
    Plots.scatter!(
        p.subplots[3],
        [ustrip(dip.x) * 1000],
        [ustrip(dip.z) * 1000],
        m = m,
        c = c,
        lab = l,
        legend = true;
        kwargs...,
    )

    return p
end

function Neuroimaging.plot(p::Plots.Plot, dips::Vector{Dipole}; l = "", kwargs...)

    for dip = 1:length(dips)
        if dip == length(dips)
            p = Neuroimaging.plot(p, dips[dip], l = l; kwargs...)
        else
            p = Neuroimaging.plot(p, dips[dip], l = ""; kwargs...)
        end
    end

    return p
end


text:
# :src isn't tested; changes depending on file/web origin
function compare_SeisHdr(H1::SeisHdr, H2::SeisHdr)
  for f in fieldnames(EQLoc)
    if typeof(getfield(H1.loc, f)) <: Union{AbstractFloat, Int64}
      (f == :rms) && continue
      @test isapprox(getfield(H1.loc, f), getfield(H2.loc,f))
    end
  end
  @test getfield(H1.loc, :typ) == getfield(H2.loc, :typ)
  @test getfield(H1.loc, :src) == getfield(H2.loc, :src)
  @test H1.mag.val ≈ H2.mag.val
  @test H1.mag.gap ≈ H2.mag.gap
  @test H1.mag.src == H2.mag.src
  @test H1.mag.scale == H2.mag.scale
  @test H1.mag.nst == H2.mag.nst
  @test H1.id == H2.id
  @test H1.ot == H2.ot
  @test H1.typ == H2.typ
  return nothing
end

function compare_SeisSrc(R1::SeisSrc, R2::SeisSrc)
  @test R1.id == R2.id
  @test R1.eid == R2.eid
  @test R1.m0 ≈ R2.m0
  @test R1.mt ≈ R2.mt
  @test R1.dm ≈ R2.dm
  @test R1.gap ≈ R2.gap
  @test R1.pax ≈ R2.pax
  @test R1.planes ≈ R2.planes
  @test R1.st.desc == R2.st.desc
  @test R1.st.dur ≈ R2.st.dur
  @test R1.st.rise ≈ R2.st.rise
  @test R1.st.decay ≈ R2.st.decay
  return nothing
end

function compare_SeisData(S1::SeisData, S2::SeisData)
  sort!(S1)
  sort!(S2)
  @test S1.id == S2.id
  @test S1.name == S2.name
  @test S1.units == S2.units
  @test isapprox(S1.fs, S2.fs)
  @test isapprox(S1.gain, S2.gain)
  for i in 1:S1.n
    L1 = S1.loc[i]
    L2 = S2.loc[i]
    @test isapprox(L1.lat, L2.lat)
    @test isapprox(L1.lon, L2.lon)
    @test isapprox(L1.el, L2.el)
    @test isapprox(L1.dep, L2.dep)
    @test isapprox(L1.az, L2.az)
    @test isapprox(L1.inc, L2.inc)

    R1 = S1.resp[i]
    R2 = S2.resp[i]
    for f in fieldnames(PZResp)
      @test isapprox(getfield(R1, f), getfield(R2,f))
    end

    # Changed 2020-03-05
    t1 = t_expand(S1.t[i], S1.fs[i])
    t2 = t_expand(S2.t[i], S2.fs[i])
    ii = sortperm(t1)
    jj = sortperm(t2)
    @test isapprox(t1[ii], t2[jj])
    @test isapprox(S1.x[i][ii], S2.x[i][jj])
  end
  return nothing
end
#=
  old :t, :x tests:
  @test S1.t[i] == S2.t[i]
  @test isapprox(S1.x[i],S2.x[i])

  reason for change:
  with gaps, representations of sample times aren't unique.

  in rare cases, writing and rereading :t, :x to/from ASDF yields a different
    time matrix :t but the data and sample times haven't changed.
=#

function compare_events(Ev1::SeisEvent, Ev2::SeisEvent)
  compare_SeisHdr(Ev1.hdr, Ev2.hdr)
  compare_SeisSrc(Ev1.source, Ev2.source)
  S1 = convert(SeisData, Ev1.data)
  S2 = convert(SeisData, Ev2.data)
  compare_SeisData(S1, S2)
  return nothing
end


text:
"""
    static_analysis(assembly; kwargs...)

Perform a static analysis of the system of nonlinear beams contained in
`assembly`. Return the resulting system and a flag indicating whether the
iteration procedure converged.

# Keyword Arguments
- `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`:
       A dictionary with keys corresponding to the points at
       which prescribed conditions are applied and elements of type
       [`PrescribedConditions`](@ref) which describe the prescribed conditions
       at those points.  If time varying, this input may be provided as a
       function of time.
- `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: A dictionary
       with keys corresponding to the elements to which distributed loads are
       applied and elements of type [`DistributedLoads`](@ref) which describe
       the distributed loads at those points.  If time varying, this input may
       be provided as a function of time.
 - `linear = false`: Set to `true` for a linear analysis
 - `linearization_state`: Linearization state variables.  Defaults to zeros.
 - `update_linearization_state`: Flag indicating whether to update the linearization state
    variables for a linear analysis with the instantaneous state variables.
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.BackTracking(maxstep=1e6)`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `tvec = 0`: Time vector/value. May be used in conjunction with time varying
     prescribed conditions and distributed loads to gradually increase
    displacements/loads.
 - `reset_state = true`: Flag indicating whether the state variables should be
    reset prior to performing the analysis.  This keyword argument is only valid
    for the pre-allocated version of this function.
"""
function static_analysis(assembly;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    tvec=0.0,
    )

    static = true

    pc = typeof(prescribed_conditions) <: AbstractDict ? prescribed_conditions : prescribed_conditions(tvec[1])

    system = System(assembly, static; prescribed_points=keys(pc))

    return static_analysis!(system, assembly;
        prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads,
        linear=linear,
        linearization_state=linearization_state,
        update_linearization_state=update_linearization_state,
        method=method,
        linesearch=linesearch,
        ftol=ftol,
        iterations=iterations,
        tvec=tvec,
        reset_state=false)
end

"""
    static_analysis!(system, assembly; kwargs...)

Pre-allocated version of `static_analysis`.
"""
function static_analysis!(system, assembly;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    tvec=0.0,
    reset_state=true)

    # check to make sure system is static
    @assert system.static == true

    # reset state, if specified
    if reset_state
        reset_state!(system)
    end

    # unpack pre-allocated storage and pointers
    x = system.x
    F = system.r
    J = system.K
    force_scaling = system.force_scaling
    mass_scaling = system.mass_scaling
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem

    converged = true
    for t in tvec

        # update stored time
        system.t = t

        # current parameters
        pcond = typeof(prescribed_conditions) <: AbstractDict ?
            prescribed_conditions : prescribed_conditions(t)
        dload = typeof(distributed_loads) <: AbstractDict ?
            distributed_loads : distributed_loads(t)

        # solve the system of equations
        f! = (F, x) -> system_residual!(F, x, assembly, pcond, dload, force_scaling,
            mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

        j! = (J, x) -> system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
            mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

        if linear
            # linear analysis
            if !update_linearization_state
                if isnothing(linearization_state)
                    x .= 0
                else
                    x .= linearization_state
                end
            end
            f!(F, x)
            j!(J, x)
            x .-= safe_lu(J) \ F
        else
            # nonlinear analysis
            df = NLsolve.OnceDifferentiable(f!, j!, x, F, J)

            result = NLsolve.nlsolve(df, x,
                linsolve=(x, A, b) -> ldiv!(x, safe_lu(A), b),
                method=method,
                linesearch=linesearch,
                ftol=ftol,
                iterations=iterations)

            # update the solution
            x .= result.zero
            J .= df.DF

            # update convergence flag
            converged = result.f_converged
        end
    end

    return system, converged
end

"""
    steady_state_analysis(assembly; kwargs...)

Perform a steady-state analysis for the system of nonlinear beams contained in
`assembly`.  Return the resulting system and a flag indicating whether the
iteration procedure converged.

# Keyword Arguments
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`:
        A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: A dictionary
        with keys corresponding to the elements to which distributed loads are
        applied and elements of type [`DistributedLoads`](@ref) which describe
        the distributed loads at those points.  If time varying, this input may
        be provided as a function of time.
 - `linear = false`: Set to `true` for a linear analysis
 - `linearization_state`: Linearization state variables.  Defaults to zeros.
 - `update_linearization_state`: Flag indicating whether to update the linearization state
    variables for a linear analysis with the current state variables.
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.LineSearches.BackTracking(maxstep=1e6)`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `origin = zeros(3)`: Global frame origin vector. If time varying, this input
        may be provided as a function of time.
 - `linear_velocity = zeros(3)`: Global frame linear velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `angular_velocity = zeros(3)`: Global frame angular velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `tvec = 0.0`: Time vector/value. May be used in conjunction with time varying
    prescribed conditions, distributed loads, and global motion to gradually
    increase displacements/loads.
 - `reset_state = true`: Flag indicating whether the state variables should be
    reset prior to performing the analysis.  This keyword argument is only valid
    for the pre-allocated version of this function.
"""
function steady_state_analysis(assembly;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    tvec=0.0,
    )

    static = false

    pc = typeof(prescribed_conditions) <: AbstractDict ? prescribed_conditions : prescribed_conditions(tvec[1])

    system = System(assembly, static; prescribed_points=keys(pc))

    return steady_state_analysis!(system, assembly;
        prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads,
        linear=linear,
        linearization_state=linearization_state,
        update_linearization_state=update_linearization_state,
        method=method,
        linesearch=linesearch,
        ftol=ftol,
        iterations=iterations,
        origin=origin,
        linear_velocity=linear_velocity,
        angular_velocity=angular_velocity,
        tvec=tvec,
        reset_state=false,
        )
end

"""
    steady_state_analysis!(system, assembly; kwargs...)

Pre-allocated version of `steady_state_analysis`.
"""
function steady_state_analysis!(system, assembly;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    tvec=0.0,
    reset_state=true,
    )

    # check to make sure the simulation is dynamic
    @assert system.static == false

    # reset state, if specified
    if reset_state
        reset_state!(system)
    end

    # unpack pointers to pre-allocated storage
    x = system.x
    F = system.r
    J = system.K
    force_scaling = system.force_scaling
    mass_scaling = system.mass_scaling
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem

    converged = true
    for t in tvec

        # update stored time
        system.t = t

        # current parameters
        pcond = typeof(prescribed_conditions) <: AbstractDict ?
            prescribed_conditions : prescribed_conditions(t)
        dload = typeof(distributed_loads) <: AbstractDict ?
            distributed_loads : distributed_loads(t)
        x0 = typeof(origin) <: AbstractVector ? SVector{3}(origin) : SVector{3}(origin(t))
        v0 = typeof(linear_velocity) <: AbstractVector ? SVector{3}(linear_velocity) : SVector{3}(linear_velocity(t))
        ω0 = typeof(angular_velocity) <: AbstractVector ? SVector{3}(angular_velocity) : SVector{3}(angular_velocity(t))

        f! = (F, x) -> system_residual!(F, x, assembly, pcond,
            dload, force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1,
            irow_elem2, icol_point, icol_elem, x0, v0, ω0)

        j! = (J, x) -> system_jacobian!(J, x, assembly, pcond,
            dload, force_scaling, mass_scaling, irow_point, irow_elem, irow_elem1,
            irow_elem2, icol_point, icol_elem, x0, v0, ω0)

        # solve the system of equations
        if linear
            # linear analysis
            if !update_linearization_state
                if isnothing(linearization_state)
                    x .= 0
                else
                    x .= linearization_state
                end
            end
            f!(F, x)
            j!(J, x)
            x .-= safe_lu(J) \ F
        else
            # nonlinear analysis
            df = NLsolve.OnceDifferentiable(f!, j!, x, F, J)

            result = NLsolve.nlsolve(df, x,
                linsolve=(x, A, b) -> ldiv!(x, safe_lu(A), b),
                method=method,
                linesearch=linesearch,
                ftol=ftol,
                iterations=iterations)

            # update the solution
            x .= result.zero
            J .= df.DF

            # update the convergence flag
            convergence = result.f_converged
        end
    end

    return system, converged
end

"""
    eigenvalue_analysis(assembly; kwargs...)

Compute the eigenvalues and eigenvectors of the system of nonlinear beams
contained in `assembly`.  Return the modified system, eigenvalues, eigenvectors,
and a convergence flag indicating whether the corresponding steady-state analysis
converged.

# Keyword Arguments
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`:
        A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: A dictionary
        with keys corresponding to the elements to which distributed loads are
        applied and elements of type [`DistributedLoads`](@ref) which describe
        the distributed loads at those points.  If time varying, this input may
        be provided as a function of time.
 - `linear = false`: Set to `true` for a linear analysis
 - `linearization_state`: Linearization state variables.  Defaults to zeros.
 - `update_linearization_state`: Flag indicating whether to update the linearization state
    variables for a linear analysis with the current state variables.
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.LineSearches.BackTracking(maxstep=1e6)`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `reset_state = true`: Flag indicating whether the state variables should be
    reset prior to performing the steady-state analysis.  This keyword argument
    is only valid for the pre-allocated version of this function.
 - `find_steady_state = reset_state && !linear`: Flag indicating whether the
    steady state solution should be found prior to performing the eigenvalue analysis.
 - `origin = zeros(3)`: Global frame origin.
    If time varying, this vector may be provided as a function of time.
 - `linear_velocity = zeros(3)`: Global frame linear velocity vector.
    If time varying, this vector may be provided as a function of time.
    May be provided either as a constant or as a function of time.
 - `angular_velocity = zeros(3)`: Global frame angular velocity vector.
    If time varying, this vector may be provided as a function of time.
 - `tvec`: Time vector. May be used in conjunction with time varying
    prescribed conditions, distributed loads, and global motion to gradually
    increase displacements/loads during the steady-state analysis.
 - `nev = 6`: Number of eigenvalues to compute
"""
function eigenvalue_analysis(assembly;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    method=:newton,
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    find_steady_state=!linear,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    tvec=0.0,
    nev=6
    )

    static = false

    pc = typeof(prescribed_conditions) <: AbstractDict ? prescribed_conditions : prescribed_conditions(tvec[1])

    system = System(assembly, static; prescribed_points=keys(pc))

    return eigenvalue_analysis!(system, assembly;
        prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads,
        linear=linear,
        linearization_state=linearization_state,
        update_linearization_state=update_linearization_state,
        method=method,
        linesearch=linesearch,
        ftol=ftol,
        iterations=iterations,
        reset_state=false,
        find_steady_state=find_steady_state,
        origin=origin,
        linear_velocity=linear_velocity,
        angular_velocity=angular_velocity,
        tvec=tvec,
        nev=nev,
        )
end

"""
    eigenvalue_analysis!(system, assembly; kwargs...)

Pre-allocated version of `eigenvalue_analysis`.  Uses the state variables stored in
`system` as an initial guess for iterating to find the steady state solution.
"""
function eigenvalue_analysis!(system, assembly;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    reset_state=true,
    find_steady_state=!linear && reset_state,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    tvec=0.0,
    nev=6,
    )

    if reset_state
        reset_state!(system)
    end

    # perform steady state analysis (if nonlinear)
    if find_steady_state
        system, converged = steady_state_analysis!(system, assembly;
            prescribed_conditions=prescribed_conditions,
            distributed_loads=distributed_loads,
            linear=linear,
            linearization_state=linearization_state,
            update_linearization_state=update_linearization_state,
            method=method,
            linesearch=linesearch,
            ftol=ftol,
            iterations=iterations,
            origin=origin,
            linear_velocity=linear_velocity,
            angular_velocity=angular_velocity,
            tvec=tvec,
            )
    else
        # set linearization state variables
        if linear && !update_linearization_state
            if isnothing(linearization_state)
                system.x .= 0
            else
                system.x .= linearization_state
            end
        end
        # converged by default
        converged = true
    end

    # unpack state vector, stiffness, and mass matrices
    x = system.x # populated during steady state solution
    K = system.K # needs to be updated
    M = system.M # still needs to be populated

    # unpack scaling parameter
    force_scaling = system.force_scaling
    mass_scaling = system.mass_scaling

    # also unpack system indices
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem

    # current time
    t = system.t

    # current parameters
    pcond = typeof(prescribed_conditions) <: AbstractDict ?
        prescribed_conditions : prescribed_conditions(t)
    dload = typeof(distributed_loads) <: AbstractDict ?
        distributed_loads : distributed_loads(t)
    x0 = typeof(origin) <: AbstractVector ? SVector{3}(origin) : SVector{3}(origin(t))
    v0 = typeof(linear_velocity) <: AbstractVector ? SVector{3}(linear_velocity) : SVector{3}(linear_velocity(t))
    ω0 = typeof(angular_velocity) <: AbstractVector ? SVector{3}(angular_velocity) : SVector{3}(angular_velocity(t))

    # solve for the system stiffness matrix
    K = system_jacobian!(K, x, assembly, pcond, dload, force_scaling, mass_scaling,
        irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem, x0, v0, ω0)

    # solve for the system mass matrix
    M = system_mass_matrix!(M, x, assembly, force_scaling, mass_scaling, irow_point, irow_elem,
        irow_elem1, irow_elem2, icol_point, icol_elem)

    # construct linear map
    T = eltype(system)
    nx = length(x)
    Kfact = safe_lu(K)
    f! = (b, x) -> ldiv!(b, Kfact, M * x)
    fc! = (b, x) -> mul!(b, M', Kfact' \ x)
    A = LinearMap{T}(f!, fc!, nx, nx; ismutating=true)

    # compute eigenvalues and eigenvectors
    λ, V = partialeigen(partialschur(A; nev=min(nx, nev), which=LM())[1])

    # sort eigenvalues by magnitude
    perm = sortperm(λ, by=(λ) -> (abs(λ), imag(λ)), rev=true)
    λ .= λ[perm]
    V .= V[:,perm]

    # eigenvalues are actually -1/λ, no modification necessary for eigenvectors
    λ .= -1 ./ λ

    return system, λ, V, converged
end

"""
    initial_condition_analysis(assembly, t0; kwargs...)

Perform an analysis to obtain a consistent set of initial conditions.  Return the
final system with the new initial conditions.

# Keyword Arguments
 - `prescribed_conditions: A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads: A dictionary with keys corresponding to the elements to
        which distributed loads are applied and elements of type
        [`DistributedLoads`](@ref) which describe the distributed loads at those
        points.  If time varying, this input may be provided as a function of
        time.
 - `linear = false`: Set to `true` for a linear analysis
 - `linearization_state`: Linearization state variables.  Defaults to zeros.
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.LineSearches.BackTracking(maxstep=1e6)`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `reset_state = true`: Flag indicating whether the state variables should be
    reset prior to performing the analysis.  This keyword argument is only valid
    for the pre-allocated version of this function.
 - `origin = zeros(3)`: Global frame origin.
    If time varying, this vector may be provided as a function of time.
 - `linear_velocity = zeros(3)`: Global frame linear velocity vector.
    If time varying, this vector may be provided as a function of time.
    May be provided either as a constant or as a function of time.
 - `angular_velocity = zeros(3)`: Global frame angular velocity vector.
    If time varying, this vector may be provided as a function of time.
 - `u0=fill(zeros(3), length(assembly.elements))`: Initial displacment of each beam element,
 - `theta0=fill(zeros(3), length(assembly.elements))`: Initial angular displacement of each beam element,
 - `udot0=fill(zeros(3), length(assembly.elements))`: Initial time derivative with respect to `u`
 - `thetadot0=fill(zeros(3), length(assembly.elements))`: Initial time derivative with respect to `theta`
 - `save=1:length(tvec)`: Steps at which to save the time history
"""
function initial_condition_analysis(assembly, t0;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    u0=fill((@SVector zeros(3)), length(assembly.elements)),
    theta0=fill((@SVector zeros(3)), length(assembly.elements)),
    udot0=fill((@SVector zeros(3)), length(assembly.elements)),
    thetadot0=fill((@SVector zeros(3)), length(assembly.elements)),
    )

    static = false

    pc = typeof(prescribed_conditions) <: AbstractDict ? prescribed_conditions : prescribed_conditions(t0)

    system = System(assembly, static; prescribed_points=keys(pc))

    return initial_condition_analysis!(system, assembly, t0;
        prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads,
        linear=linear,
        linearization_state=linearization_state,
        method=method,
        linesearch=linesearch,
        ftol=ftol,
        iterations=iterations,
        reset_state=false,
        origin=origin,
        linear_velocity=linear_velocity,
        angular_velocity=angular_velocity,
        u0=u0,
        theta0=theta0,
        udot0=udot0,
        thetadot0=thetadot0,
        )
end

"""
    initial_condition_analysis!(system, assembly, t0; kwargs...)

Pre-allocated version of `initial_condition_analysis`.
"""
function initial_condition_analysis!(system, assembly, t0;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    reset_state=true,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    u0=fill((@SVector zeros(3)), length(assembly.elements)),
    theta0=fill((@SVector zeros(3)), length(assembly.elements)),
    udot0=fill((@SVector zeros(3)), length(assembly.elements)),
    thetadot0=fill((@SVector zeros(3)), length(assembly.elements)),
    )

    # check to make sure the simulation is dynamic
    @assert system.static == false

    if reset_state
        reset_state!(system)
    end

    # unpack pre-allocated storage and pointers for system
    x = system.x
    F = system.r
    J = system.K
    force_scaling = system.force_scaling
    mass_scaling = system.mass_scaling
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem
    udot = system.udot
    θdot = system.θdot
    Pdot = system.Pdot
    Hdot = system.Hdot

    nelem = length(assembly.elements)

    # set current time step
    system.t = t0

    # set current parameters
    pcond = typeof(prescribed_conditions) <: AbstractDict ?
        prescribed_conditions : prescribed_conditions(t0)
    dload = typeof(distributed_loads) <: AbstractDict ?
        distributed_loads : distributed_loads(t0)
    x0 = typeof(origin) <: AbstractVector ? SVector{3}(origin) : SVector{3}(origin(t0))
    v0 = typeof(linear_velocity) <: AbstractVector ? SVector{3}(linear_velocity) : SVector{3}(linear_velocity(t0))
    ω0 = typeof(angular_velocity) <: AbstractVector ? SVector{3}(angular_velocity) : SVector{3}(angular_velocity(t0))

    # construct residual and jacobian functions
    f! = (F, x) -> system_residual!(F, x, assembly, pcond, dload, force_scaling,
        mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, u0, theta0, udot0, thetadot0)

    j! = (J, x) -> system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
        mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
        x0, v0, ω0, u0, theta0, udot0, thetadot0)

    # solve system of equations
    if linear
        # linear analysis
        if isnothing(linearization_state)
            x .= 0
        else
            x .= linearization_state
        end
        f!(F, x)
        j!(J, x)
        x .-= safe_lu(J) \ F
    else
        # nonlinear analysis
        df = OnceDifferentiable(f!, j!, x, F, J)

        result = NLsolve.nlsolve(df, x,
            linsolve=(x, A, b) -> ldiv!(x, safe_lu(A), b),
            method=method,
            linesearch=linesearch,
            ftol=ftol,
            iterations=iterations)

        x .= result.zero
        J .= df.DF
    end

    # get convergence flag
    converged = result.f_converged

    # save states and state rates
    for ielem = 1:nelem
        icol = icol_elem[ielem]
        # extract rotation parameters for this beam element
        C = get_C(theta0[ielem])
        Cab = assembly.elements[ielem].Cab
        CtCab = C' * Cab
        # save states and state rates
        udot[ielem] = udot0[ielem]
        θdot[ielem] = thetadot0[ielem]
        Pdot[ielem] = CtCab' * SVector(x[icol], x[icol + 1], x[icol + 2]) .* mass_scaling
        Hdot[ielem] = CtCab' * SVector(x[icol + 3], x[icol + 4], x[icol + 5]) .* mass_scaling
# restore original state vector
        x[icol:icol + 2] .= u0[ielem]
        x[icol + 3:icol + 5] .= theta0[ielem]
    end

    return system, converged
end

"""
    time_domain_analysis(assembly, tvec; kwargs...)

Perform a time-domain analysis for the system of nonlinear beams contained in
`assembly` using the time vector `tvec`.  Return the final system, a post-processed
solution history, and a convergence flag indicating whether the iterations
converged for each time step.

# Keyword Arguments
 - `prescribed_conditions: A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads: A dictionary with keys corresponding to the elements to
        which distributed loads are applied and elements of type
        [`DistributedLoads`](@ref) which describe the distributed loads at those
        points.  If time varying, this input may be provided as a function of
        time.
 - `linear = false`: Set to `true` for a linear analysis
 - `linearization_state`: Linearization state variables.  Defaults to zeros.
 - `update_linearization_state`: Flag indicating whether to update the linearization state
    variables for a linear analysis with the current state variables.
 - `method = :newton`: Method (as defined in NLsolve) to solve nonlinear system of equations
 - `linesearch = LineSearches.LineSearches.BackTracking(maxstep=1e6)`: Line search used to solve nonlinear system of equations
 - `ftol = 1e-9`: tolerance for solving nonlinear system of equations
 - `iterations = 1000`: maximum iterations for solving the nonlinear system of equations
 - `reset_state = true`: Flag indicating whether the state variables should be
    reset prior to performing the analysis.  This keyword argument is only valid
    for the pre-allocated version of this function.
 - `initialize = true`: Flag indicating whether a consistent set of initial
    conditions should be found using [`initial_condition_analysis`](@ref). If
    `false`, the keyword arguments `u0`, `theta0`, `udot0` and `thetadot0` will
    be ignored and the system state vector will be used as the initial state
    variables.
 - `origin`: Global frame origin vector. If time varying, this input
        may be provided as a function of time.
 - `linear_velocity`: Global frame linear velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `angular_velocity`: Global frame angular velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `u0=fill(zeros(3), length(assembly.elements))`: Initial displacment of each beam element,
 - `theta0=fill(zeros(3), length(assembly.elements))`: Initial angular displacement of each beam element,
 - `udot0=fill(zeros(3), length(assembly.elements))`: Initial time derivative with respect to `u`
 - `thetadot0=fill(zeros(3), length(assembly.elements))`: Initial time derivative with respect to `theta`
 - `save=1:length(tvec)`: Steps at which to save the time history
"""
function time_domain_analysis(assembly, tvec;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    initialize=true,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    u0=fill((@SVector zeros(3)), length(assembly.elements)),
    theta0=fill((@SVector zeros(3)), length(assembly.elements)),
    udot0=fill((@SVector zeros(3)), length(assembly.elements)),
    thetadot0=fill((@SVector zeros(3)), length(assembly.elements)),
    save=1:length(tvec)
    )

    static = false

    pc = typeof(prescribed_conditions) <: AbstractDict ? prescribed_conditions : prescribed_conditions(tvec[1])

    system = System(assembly, static; prescribed_points=keys(pc))

    return time_domain_analysis!(system, assembly, tvec;
        prescribed_conditions=prescribed_conditions,
        distributed_loads=distributed_loads,
        linear=linear,
        linearization_state=linearization_state,
        update_linearization_state=update_linearization_state,
        method=method,
        linesearch=linesearch,
        ftol=ftol,
        iterations=iterations,
        reset_state=false,
        initialize=initialize,
        origin=origin,
        linear_velocity=linear_velocity,
        angular_velocity=angular_velocity,
        u0=u0,
        theta0=theta0,
        udot0=udot0,
        thetadot0=thetadot0,
        save=save,
        )
end

"""
    time_domain_analysis!(system, assembly, tvec; kwargs...)

Pre-allocated version of `time_domain_analysis`.
"""
function time_domain_analysis!(system, assembly, tvec;
    prescribed_conditions=Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads=Dict{Int,DistributedLoads{Float64}}(),
    linear=false,
    linearization_state=nothing,
    update_linearization_state=false,
    method=:newton,
    linesearch=LineSearches.BackTracking(maxstep=1e6),
    ftol=1e-9,
    iterations=1000,
    reset_state=true,
    initialize=true,
    origin=(@SVector zeros(3)),
    linear_velocity=(@SVector zeros(3)),
    angular_velocity=(@SVector zeros(3)),
    u0=fill((@SVector zeros(3)), length(assembly.elements)),
    theta0=fill((@SVector zeros(3)), length(assembly.elements)),
    udot0=fill((@SVector zeros(3)), length(assembly.elements)),
    thetadot0=fill((@SVector zeros(3)), length(assembly.elements)),
    save=1:length(tvec)
    )

    # check to make sure the simulation is dynamic
    @assert system.static == false

    if reset_state
        reset_state!(system)
    end

    # perform initial condition analysis
    if initialize
        system, converged = initial_condition_analysis!(system, assembly, tvec[1];
            prescribed_conditions=prescribed_conditions,
            distributed_loads=distributed_loads,
            linear=linear,
            linearization_state=linearization_state,
            method=method,
            linesearch=linesearch,
            ftol=ftol,
            iterations=iterations,
            reset_state=false,
            origin=origin,
            linear_velocity=linear_velocity,
            angular_velocity=angular_velocity,
            u0=u0,
            theta0=theta0,
            udot0=udot0,
            thetadot0=thetadot0,
            )
    else
        # converged by default
        converged = true
    end

    # unpack pre-allocated storage and pointers for system
    x = system.x
    F = system.r
    J = system.K
    force_scaling = system.force_scaling
    mass_scaling = system.mass_scaling
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem
    udot = system.udot
    θdot = system.θdot
    Pdot = system.Pdot
    Hdot = system.Hdot

    # number of beam elements
    nelem = length(assembly.elements)

    # initialize storage for each time step
    isave = 1
    history = Vector{AssemblyState{eltype(system)}}(undef, length(save))

    # add initial state to the solution history
    if isave in save
        pcond = typeof(prescribed_conditions) <: AbstractDict ?
            prescribed_conditions : prescribed_conditions(tvec[1])
        history[isave] = AssemblyState(system, assembly, prescribed_conditions=pcond)
        isave += 1
    end

    # --- Begin Time Domain Simulation --- #

    for it = 2:length(tvec)

        # update current time
        system.t = tvec[it]

        # current time step size
        dt = tvec[it] - tvec[it - 1]

        # current parameters
        pcond = typeof(prescribed_conditions) <: AbstractDict ?
            prescribed_conditions : prescribed_conditions(tvec[it])
        dload = typeof(distributed_loads) <: AbstractDict ?
            distributed_loads : distributed_loads(tvec[it])
        x0 = typeof(origin) <: AbstractVector ? SVector{3}(origin) : SVector{3}(origin(tvec[it]))
        v0 = typeof(linear_velocity) <: AbstractVector ? SVector{3}(linear_velocity) : SVector{3}(linear_velocity(tvec[it]))
        ω0 = typeof(angular_velocity) <: AbstractVector ? SVector{3}(angular_velocity) : SVector{3}(angular_velocity(tvec[it]))

        # current initialization parameters
        for ielem = 1:nelem
            icol = icol_elem[ielem]
            # get beam element states
            u = SVector(x[icol], x[icol + 1], x[icol + 2])
            θ = SVector(x[icol + 3], x[icol + 4], x[icol + 5])
            P = SVector(x[icol + 12], x[icol + 13], x[icol + 14]) .* mass_scaling
            H = SVector(x[icol + 15], x[icol + 16], x[icol + 17]) .* mass_scaling
            # extract rotation parameters
            C = get_C(θ)
            Cab = assembly.elements[ielem].Cab
            CtCab = C' * Cab
            # store `udot_init` in `udot`
            udot[ielem] = 2 / dt * u + udot[ielem]
            # store `θdot_init` in `θdot`
            θdot[ielem] = 2 / dt * θ + θdot[ielem]
            # store `CtCabPdot_init` in `Pdot`
            Pdot[ielem] = 2 / dt * CtCab * P + CtCab * Pdot[ielem]
            # store `CtCabHdot_init` in `Hdot`
            Hdot[ielem] = 2 / dt * CtCab * H + CtCab * Hdot[ielem]
        end

        # solve for the state variables at the next time step
        f! = (F, x) -> system_residual!(F, x, assembly, pcond, dload, force_scaling,
            mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
            x0, v0, ω0, udot, θdot, Pdot, Hdot, dt)

        j! = (J, x) -> system_jacobian!(J, x, assembly, pcond, dload, force_scaling,
            mass_scaling, irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem,
            x0, v0, ω0, udot, θdot, Pdot, Hdot, dt)

        # solve system of equations
        if linear
            # linear analysis
            if !update_linearization_state
                if isnothing(linearization_state)
                    x .= 0
                else
                    x .= linearization_state
                end
            end
            f!(F, x)
            j!(J, x)
            x .-= safe_lu(J) \ F
        else
            df = OnceDifferentiable(f!, j!, x, F, J)

            result = NLsolve.nlsolve(df, x,
                linsolve=(x, A, b) -> ldiv!(x, safe_lu(A), b),
                method=method,
                linesearch=linesearch,
                ftol=ftol,
                iterations=iterations)

            x .= result.zero
            J .= df.DF
        end

        # set new state rates
        for ielem = 1:nelem
            icol = icol_elem[ielem]
            # get beam element states
            u = SVector(x[icol], x[icol + 1], x[icol + 2])
            θ = SVector(x[icol + 3], x[icol + 4], x[icol + 5])
            P = SVector(x[icol + 12], x[icol + 13], x[icol + 14]) .* mass_scaling
            H = SVector(x[icol + 15], x[icol + 16], x[icol + 17]) .* mass_scaling
            # extract rotation parameters
            C = get_C(θ)
            Cab = assembly.elements[ielem].Cab
            CtCab = C' * Cab
            # save state rates
            udot[ielem] = 2 / dt * u - udot[ielem]
            θdot[ielem] = 2 / dt * θ - θdot[ielem]
            Pdot[ielem] = 2 / dt * P - CtCab' * Pdot[ielem]
            Hdot[ielem] = 2 / dt * H - CtCab' * Hdot[ielem]
        end

        # add state to history
        if it in save
            history[isave] = AssemblyState(system, assembly, prescribed_conditions=prescribed_conditions)
            isave += 1
        end

        # stop early if unconverged
        if !linear && !result.f_converged
            converged = false
            break
        end

    end

    return system, history, converged
end


text:
# Key types

struct AES256Key <: AbstractAESKey
	key::Array{UInt8, 1}
end

struct AES192Key <: AbstractAESKey
	key::Array{UInt8, 1}
end

struct AES128Key <: AbstractAESKey
	key::Array{UInt8, 1}
end

@inline convert_key(key::AbstractAESKey, key_length) = key
@inline convert_key(key::Array{UInt8,1}, ::Val{128}) = AES128Key(key)
@inline convert_key(key::Array{UInt8,1}, ::Val{192}) = AES192Key(key)
@inline convert_key(key::Array{UInt8,1}, ::Val{256}) = AES256Key(key)
@inline convert_key(key::String, key_length) = convert_key(Array{UInt8}(key), key_length)

@inline Base.getindex(k::AbstractAESKey, i) = k.key[i]

@inline get_key_length(key::Array{UInt8}) = 8 * length(key)
@inline get_key_length(key::AES128Key) = 128
@inline get_key_length(key::AES192Key) = 192
@inline get_key_length(key::AES256Key) = 256

# Cache Type
mutable struct CipherCache{kType,tmpType,modecacheType} <: AbstractAESCache
	""" registers to contain intermediate keys """
	K::kType
	""" temporary word """
	tmp::tmpType
	""" cache specific to mode """
	modecache::modecacheType
end

@inline get_intermediate_words(key::AES128Key) = @SArray[ zeros(UInt8, 4) for i in 1:4 ]
@inline get_intermediate_words(key::AES192Key) = @SArray[ zeros(UInt8, 4) for i in 1:6 ]
@inline get_intermediate_words(key::AES256Key) = @SArray[ zeros(UInt8, 4) for i in 1:4 ]

@inline get_modecache(m::Val{CTR}) = similar(Array{UInt8}, 16)
@inline get_modecache(m) = nothing

function gen_cache(key::AbstractAESKey, mode=CBC)
	K = get_intermediate_words(key)
	tmp = zeros(UInt8, 4)
	modecache = get_modecache(Val(mode))
	return CipherCache(K, tmp, modecache)
end

# Ciphertext type
mutable struct CipherText{ivType,mode}
	data::Array{UInt8}
	iv::ivType
	keylength::Int
end
CipherText(data, iv::IV, keylength, mode::MODE) where {IV} = CipherText{IV,mode}(data, iv, keylength)
const CT = CipherText

# Cipher Type
mutable struct Cipher{mode,cacheType,keyType} <: AbstractCipher
	cache::cacheType
	key::keyType
end
const AESCipher = Cipher

function Cipher(;key_length=128,mode=CBC,key=keygen(key_length))
	if !is_valid_key_length(key_length)
		error("$key_length is an invalid key length. Key length can be 128, 196 or 256.")
	end

	if get_key_length(key) !==  key_length
		error("Provided Key and Key Length do not agree.")
	end

	aes_key = convert_key(key, Val(key_length))

	cache = gen_cache(aes_key, mode)

	Cipher{mode,typeof(cache),typeof(aes_key)}(cache,aes_key)
end

@inline isecb(cipher::Cipher{m,c,k}) where {m,c,k} = m == ECB
@inline iscbc(cipher::Cipher{m,c,k}) where {m,c,k} = m == CBC
@inline iscfb(cipher::Cipher{m,c,k}) where {m,c,k} = m == CFB
@inline isofb(cipher::Cipher{m,c,k}) where {m,c,k} = m == OFB
@inline isctr(cipher::Cipher{m,c,k}) where {m,c,k} = m == CTR

@inline get_mode(cipher::Cipher{m,c,k}) where {m,c,k} = m
@inline get_key_length(cipher::Cipher) = get_key_length(cipher.key)

@inline needs_iv(cipher) = iscbc(cipher) || isctr(cipher)

# Block type

mutable struct AESBlock{dataType,offsetType} <: AbstractArray{UInt8, 1}
	data::dataType
	offset::offsetType
	function AESBlock(data, offset=1, len=16)
		new{typeof(data), typeof(offset)}(data, offset)
	end
end

function Base.getindex(block::B, i) where {B<:AESBlock}
	return UInt8(block.data[16*(block.offset - 1) + i])
end

function Base.setindex!(block::B, val, i) where {B<:AESBlock}
	block.data[16*(block.offset - 1) + i] = eltype(block.data)(val)
end

function Base.copyto!(block::B, arr::AbstractArray{T2,1}) where {B<:AESBlock, T2}
	for i in 1:16
		block.data[16*(block.offset - 1) + i] = eltype(block.data)(arr[i])
	end
end

@inline Base.size(block::AESBlock) = (16,)
@inline Base.length(block::AESBlock) = 16

increment!(block::B) where {B<:AESBlock} = block.offset += 1


text:
include("particle.jl")

function generateSwarm(numberOfParticles,numberOfOperands,minRang,maxRang)
    swarm = [];
    for i in 1:numberOfParticles
        particle = generateParticle(numberOfOperands,minRang,maxRang);
        push!(swarm,particle);
    end
    return swarm;
end

function  printSwarm(swarm)
    for i in 1:length(swarm)
        printParticle(swarm[i]);
    end
end

function updateSwarmVelocity!(swarm,globalBestParticle,targetValue)
    for i in 1:length(swarm)
        updateVelocity!(swarm[i],globalBestParticle,targetValue);

    end
end

function updateSwarmPosition!(swarm,targetValue)
    for i in 1:length(swarm)
        updatePosition!(swarm[i],targetValue);
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.7172866311917714 + 0.8530094465195428*1/(1 + exp(-(0.49975691401319544 + 0.817850285997455*1/(1 + exp(-(-0.9044530867886094 + -0.3504940249833268*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + -0.925387475816613*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))) + -0.15011340339301915*1/(1 + exp(-(-0.49866226232151334 + 0.31623027699891626*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + 0.1823838704363938*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))))))))) + -0.2302663292873568*1/(1 + exp(-(0.6056562038214484 + 0.017530899400252764*1/(1 + exp(-(-0.9044530867886094 + -0.3504940249833268*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + -0.925387475816613*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))) + 0.9141090581590978*1/(1 + exp(-(-0.49866226232151334 + 0.31623027699891626*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + 0.1823838704363938*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))))))))) + 1/(1 + exp(-(0.11802092343800519 + -0.23198020992542823*1/(1 + exp(-(0.49975691401319544 + 0.817850285997455*1/(1 + exp(-(-0.9044530867886094 + -0.3504940249833268*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + -0.925387475816613*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))) + -0.15011340339301915*1/(1 + exp(-(-0.49866226232151334 + 0.31623027699891626*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + 0.1823838704363938*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))))))))) + 0.3694431087223533*1/(1 + exp(-(0.6056562038214484 + 0.017530899400252764*1/(1 + exp(-(-0.9044530867886094 + -0.3504940249833268*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + -0.925387475816613*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))) + 0.9141090581590978*1/(1 + exp(-(-0.49866226232151334 + 0.31623027699891626*1/(1 + exp(-(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])))) + 0.1823838704363938*1/(1 + exp(-(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("vcm.jl")

# sample new sigma_g from InverseGamma with IG(a, b)
# see: supplemental of http://papers.nips.cc/paper/4579-a-nonparametric-variable-clustering-model
function random_sigma_g(B::VCMBuffer, γ_a::Float64, γ_b::Float64)

    Gm = (B.G.*B.C)

    a = γ_a + B.D  / 2.0
    b = γ_b + .5 * sum(Gm.*Gm)

    dist = InverseGamma(a, b)

    return sqrt( rand(dist) )

end

# sample new sigma_noise from InverseGamma with IG(a, b)
# see: supplemental of http://papers.nips.cc/paper/4579-a-nonparametric-variable-clustering-model
function random_sigma_noise(B::VCMBuffer, γ_a::Float64, γ_b::Float64)

    E = B.Y - (B.G.*B.C) * B.X

    a = ( γ_a + (B.D*B.N) ) / 2.0
    b = γ_b + .5 * trace(E'*E)

    dist = InverseGamma(a, b)

    return sqrt( rand(dist) )

end

# sample latent factor matrix X
# see: supplemental of http://papers.nips.cc/paper/4579-a-nonparametric-variable-clustering-model
function sample_latent_factors!(Y::Array, B::VCMBuffer)

    Gm = B.G .* B.C
    tem = Gm' / B.σ_noise^2

    sigma_x_p = try
            B.σ_x^(-2)
        catch
            1
        end

    λ = (B.σ_noise^(-2)) * sum( Gm .* Gm, 1 ) + sigma_x_p

    for n in 1:B.N

        μ = (tem * Y[:,n]) ./ λ'
        B.X[:,n] = λ'.^(-0.5) .* randn(B.K, 1) + μ

    end

    B
end

function sample_vcm_data(N::Int, D::Int, G::Array{Float64}, X::Array{Float64}, C::Array{Float64},
    mu_noise::Float64, sigma_noise::Float64)

    ϵ = mu_noise + sigma_noise * randn(D, N)

    return (G .* C) * X + ϵ

end

# sample cluster assignments c and lattent factor loadings G by integrating out g and instantiating X.
# See: Neal (2000) Algorithm 8 for details and
# supplemental of http://papers.nips.cc/paper/4579-a-nonparametric-variable-clustering-model
function gibbs_aux_assignment!(B::VCMBuffer; out_dim = -1)

    (D, N) = size(B.Y)

    if out_dim != -1
        perm = collect(out_dim)
    else
        perm = randperm(D)
    end

    auxDist = Normal(B.μ_x, B.σ_x)

    X2Cache = sum(B.X.^2, 2)

    for d in perm

        #println("out dim", out_dim)

        if out_dim == -1
            # get cluster of current item
            cluster = B.cc[d]

            #println("got cluster")

            # check if cluster is empty
            if sum(B.cc .== cluster) == 1

                #println("found empty cluster")

                B.K -= 1
                idx = find(B.cc .> cluster)
                B.cc[idx] -= 1

                #C = C[:,[1:cluster-1, cluster+1:end]]
                B.G = B.G[:,[1:cluster-1; cluster+1:end]]

                r = B.X[cluster,:]
                B.X = B.X[[1:cluster-1; cluster+1:end],:]

                #println("removed empty cluster")

            end
        else
            B.G = cat(1, B.G, zeros(1, B.K))
            B.cc = cat(1, B.cc, 0)
            #C = cat(1, C, zeros(Int, 1, K))
        end

        #println("generate cache...")

        Xm = reduce(hcat, map(m -> rand(auxDist, B.N), 1:B.m_aux))'
        X2m = sum(Xm.^2, 2)

        XX = cat(1, B.X, Xm)
        XX2Cache = cat(1, X2Cache, X2m)

        #println("prepare p..")

        # compute P(c_i = c | c_-i, y_i, ϕ_1, ..., ϕ_h ) according to Neal (2000)
        p = ones(B.K+B.m_aux) .* -Inf
        p[1:B.K] = log( map( c -> sum(B.cc .== c) / (B.α + D - 1), 1:B.K ) )
        p[B.K + 1:B.K + B.m_aux] = log( B.α / B.m_aux ) - log( B.α + D - 1 )

        Gsamps = zeros(1, B.K + B.m_aux)

        xy = XX * B.Y[d,:]'

        for k in 1:B.K + B.m_aux

            λ = XX2Cache[k] / B.σ_noise^2 + 1.0 / B.σ_g^2
            μ = ( xy[k] / B.σ_noise^2 + B.μ_g/B.σ_g^2 ) / λ

            p[k] = p[k] - .5 * ( log(λ) - λ * μ^2 )

            Gsamps[k] = μ + λ^(-.5) * randn()
        end

        #println("finished comp. of p...")

        # common normalization
        p = exp( p - maximum(p) )

        kk = rand_indices(p)

        #println("got new k")

        g = Gsamps[kk]

        #println("found g")

        if kk > B.K
            # open new cluster
            B.X = cat(1, B.X, XX[kk,:])
            #C = cat(2, C, zeros(Int, D, 1))
            B.G = cat(2, B.G, zeros(D, 1))

            B.K += 1
            kk = B.K

            X2Cache = sum(B.X.^2, 2)

        end

        if out_dim == -1
            B.G[d, kk] = g
            B.cc[d] = kk
            #C[d,:] = zeros(1, K)
            #C[d,kk] = 1
        else
            B.G[out_dim, kk] = g
            B.cc[out_dim] = kk
            #C[out_dim,:] = zeros(1, K)
            #C[out_dim,kk] = 1
        end

        #println("loop end")

    end

    B.C = spzeros(D, B.K)
    B.C[B.cc] = 1

    B
end


text:
struct Kernel
    st::Symbol
    locs::Symbol
    init::Vector{Any}
    brt
end

function Kernel(brt::BitRoutine)
    @gensym st locs
    return Kernel(st, locs, [], brt)
end


text:
#######################################################################################################################################################################################################
# Changes to this structure
# General
#     2021-Aug-04: refactor the structure with constants, variables, and temporary cache
#     2021-Aug-04: add concentrations and characteristic curves altogether
#     2021-Aug-10: add CBC and PRO supoort
#     2021-Agu-10: add constructors within the structure rather than initialize it externally
#     2021-Sep-30: rename LeafBio to LeafBiophysics to be more specific
#     2021-Oct-19: sort variable to prognostic and dignostic catergories
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

Struct that contains leaf biophysical traits used to run leaf reflection and transmittance.

# Fields
$(TYPEDFIELDS)

---
# Examples
```julia
lbio = Emerald.LeafBiophysics{FT}();
lbio = Emerald.LeafBiophysics{FT}(collect(FT,400:5:2500.1));
lbio = Emerald.LeafBiophysics{FT}(collect(FT,400:5:2500.1); opti=Emerald.OPTI_2021);
lbio = Emerald.LeafBiophysics{FT}(WaveLengthSet{FT}());
lbio = Emerald.LeafBiophysics{FT}(WaveLengthSet{FT}(); opti=Emerald.OPTI_2021);
```
"""
mutable struct LeafBiophysics{FT<:AbstractFloat}
    # parameters that do not change with time
    "Leaf fluorescence quantum efficiency (Fo standard)"
    FQE::FT
    "Dry matter content (dry leaf mass per unit area) `[g cm⁻²]`"
    LMA::FT
    "Carbon-based constituents in LMA `[g cm⁻²]`"
    CBC::FT
    "Protein content in LMA (PRO = LMA - CBC) `[g cm⁻²]`"
    PRO::FT
    "Specific absorption coefficients of anthocynanin `[-]`"
    K_ANT::Vector{FT}
    "Specific absorption coefficients of chlorophyll a and b `[-]`"
    K_CAB::Vector{FT}
    "Specific absorption coefficients of violaxanthin carotenoid `[-]`"
    K_CAR_V::Vector{FT}
    "Specific absorption coefficients of zeaxanthin carotenoid `[-]`"
    K_CAR_Z::Vector{FT}
    "Specific absorption coefficients of carbon-based constituents `[-]`"
    K_CBC::Vector{FT}
    "Specific absorption coefficients of water `[-]`"
    K_H₂O::Vector{FT}
    "Specific absorption coefficients of dry matter `[-]`"
    K_LMA::Vector{FT}
    "Specific absorption coefficients of protein `[-]`"
    K_PRO::Vector{FT}
    "Specific absorption coefficients of PS I and II `[-]`"
    K_PS::Vector{FT}
    "Specific absorption coefficients of senescent material `[-]`"
    K_SENES::Vector{FT}
    "Leaf mesophyll structural parameter that describes mesophyll reflection and transmittance"
    MESOPHYLL_N::FT
    "Doubling adding layers"
    NDUB::Int
    "Refractive index `[-]`"
    NR::Vector{FT}

    # prognostic variables that change with time
    "Anthocynanin content `[ug cm⁻²]`"
    ant::FT
    "Chlorophyll a and b content `[ug cm⁻²]`"
    cab::FT
    "Carotenoid content `[ug cm⁻²]`"
    car::FT
    "Senescent material (brown pigments) fraction `[-]`"
    f_senes::FT
    "Zeaxanthin fraction in Carotenoid (1=all Zeaxanthin, 0=all Violaxanthin) `[-]`"
    f_zeax
    "Equivalent water thickness `[cm]`"
    l_H₂O::FT

    # dignostic variables that change with time
    "Fluorescence excitation matrix backwards `[-]`"
    mat_b::Matrix{FT}
    "Fluorescence excitation matrix forwards `[-]`"
    mat_f::Matrix{FT}
    "Relative absorption by Chlorophyll `[-]`"
    α_cab::Vector{FT}
    "Relative absorption by Chlorophyll+Carotenoid `[-]`"
    α_cabcar::Vector{FT}
    "Broadband thermal reflectance, related to blackbody emittance `[-]`"
    ρ_LW::FT
    "Shortwave leaf reflectance `[-]`"
    ρ_SW::Vector{FT}
    "Broadband thermal transmission, related to blackbody emittance `[-]`"
    τ_LW::FT
    "Shortwave leaf transmission `[-]`"
    τ_SW::Vector{FT}

    # caches to speed up calculations
    "Specific absorption coefficients of all materials"
    _k_all::Vector{FT}
    "Shortwave absorption, 1 .- ρ_SW .- τ_SW  `[-]`"
    _α_SW::Vector{FT}

    # constructors
    LeafBiophysics{FT}(swl::Vector{FT}=FT.(WAVELENGTHS); opti::String=OPTI_2021) where {FT<:AbstractFloat} = LeafBiophysics{FT}(WaveLengthSet{FT}(swl; opti=opti); opti=opti)
    LeafBiophysics{FT}(wls::WaveLengthSet{FT}; opti::String=OPTI_2021) where {FT<:AbstractFloat} = (
        @unpack NΛ, NΛ_SIF, NΛ_SIFE, SΛ = wls;

        # read data from the MAT file
        _opti    = matread(opti)["optipar"];
        __nr     = _opti["nr"  ];
        __Klma   = _opti["Kdm" ];
        __Kcab   = _opti["Kab" ];
        __Kant   = _opti["Kant"];
        __Kh2o   = _opti["Kw"  ];
        __Ksenes = _opti["Ks"  ];
        __Kps    = _opti["phi" ];
        __KcaV   = _opti["KcaV"];
        __KcaZ   = _opti["KcaZ"];
        __lambda = _opti["wl"  ];

        # read protein and carbon-based constituents if exist (not in file OPTI_2017)
        if opti == OPTI_2017
            __Kpro = __Klma;
            __Kcbc = __Klma;
        else
            __Kpro = _opti["Kp"  ];
            __Kcbc = _opti["Kcbc"];
        end;

        # create data to parse
        _nr     = zeros(FT, NΛ);
        _Klma   = zeros(FT, NΛ);
        _Kcab   = zeros(FT, NΛ);
        _Kant   = zeros(FT, NΛ);
        _Kh2o   = zeros(FT, NΛ);
        _Ksenes = zeros(FT, NΛ);
        _Kps    = zeros(FT, NΛ);
        _KcaV   = zeros(FT, NΛ);
        _KcaZ   = zeros(FT, NΛ);
        _Kpro   = zeros(FT, NΛ);
        _Kcbc   = zeros(FT, NΛ);

        # fill in the data arrays
        @inbounds for _i in 1:NΛ
            _wo = findall( (__lambda.>=SΛ[_i]) .& (__lambda.<SΛ[_i+1]) );
            if length(_wo) < 1
                @warn "Warning, some wavelengths out of bounds $(string(SΛ[_i]))";
            end

            _nr[_i]     = nanmean(__nr[_wo]    );
            _Klma[_i]   = nanmean(__Klma[_wo]  );
            _Kcab[_i]   = nanmean(__Kcab[_wo]  );
            _Kant[_i]   = nanmean(__Kant[_wo]  );
            _Kh2o[_i]   = nanmean(__Kh2o[_wo]  );
            _Ksenes[_i] = nanmean(__Ksenes[_wo]);
            _Kps[_i]    = nanmean(__Kps[_wo]   );
            _KcaV[_i]   = nanmean(__KcaV[_wo]  );
            _KcaZ[_i]   = nanmean(__KcaZ[_wo]  );
            _Kpro[_i]   = nanmean(__Kpro[_wo]  );
            _Kcbc[_i]   = nanmean(__Kcbc[_wo]  );
        end;

        return new{FT}(1, 0.012, 0, 0, _Kant, _Kcab, _KcaV, _KcaZ, _Kcbc, _Kh2o, _Klma, _Kpro, _Kps, _Ksenes, 1.4, 10, _nr, 0, 40, 10, 0, 0, 0.01, zeros(FT,NΛ_SIF,NΛ_SIFE), zeros(FT,NΛ_SIF,NΛ_SIFE),
                       zeros(FT,NΛ), zeros(FT,NΛ), 0.01, zeros(FT,NΛ), 0.01, zeros(FT,NΛ), zeros(FT,NΛ), zeros(FT,NΛ))
    )
end

#######################################################################################################################################################################################################
# Changes to this structure
# General
#     2021-Aug-04: refactor the structure with BIO_PHYSICS as a field
#     2021-Oct-19: sort variable to prognostic and dignostic catergories
# To do
#     TODO: add leaf physiological parameters as a field well
#     TODO: add leaf hydraulics as a field as well
#     TODO: link leaf water content to BIO_PHYSICS.l_H₂O
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

Struct that contains leaf traits such as plant photosynthesis, hydraulics, and biophysics.

# Fields
$(TYPEDFIELDS)

---
# Examples
```julia
leaf = Leaf{FT}();
leaf = Leaf{FT}(collect(FT,400:5:2500.1));
leaf = Leaf{FT}(collect(FT,400:5:2500.1); opti=Emerald.OPTI_2021);
leaf = Leaf{FT}(WaveLengthSet{FT}());
leaf = Leaf{FT}(WaveLengthSet{FT}(); opti=Emerald.OPTI_2021);
```
"""
mutable struct Leaf{FT<:AbstractFloat}
    # parameters that do not change with time
    "Biophysical parameter structure"
    BIO_PHYSICS::LeafBiophysics{FT}

    # prognostic variables that change with time
    "Current leaf temperature"
    t::FT

    # dignostic variables that change with time

    # caches to speed up calculations
    "Last leaf temperature. If different from t, then make temperature correction"
    _t::FT

    # constructors
    Leaf{FT}(wls::WaveLengthSet{FT}; opti::String=OPTI_2021) where {FT<:AbstractFloat} = new{FT}(LeafBiophysics{FT}(wls; opti=opti), 298.15, 0.0)
    Leaf{FT}(swl::Vector{FT}=FT.(WAVELENGTHS); opti::String=OPTI_2021) where {FT<:AbstractFloat} = new{FT}(LeafBiophysics{FT}(swl; opti=opti), 298.15, 0.0)
end