text:
struct NLoptAlg <: AbstractOptimizer
    alg::Symbol
end

@params struct NLoptOptions
    nt::NamedTuple
end
function NLoptOptions(;
    ftol_rel = 1e-6,
    ftol_abs = 1e-6,
    xtol_rel = 1e-6,
    xtol_abs = 1e-6,
    kwargs...,
)
    NLoptOptions(
        merge(
            NamedTuple(kwargs),
            (
                ftol_rel = ftol_rel,
                ftol_abs = ftol_abs,
                xtol_rel = xtol_rel,
                xtol_abs = xtol_abs,
            )
        )
    )
end

@params mutable struct NLoptWorkspace <: Workspace
    model::VecModel
    problem::NLopt.Opt
    x0::AbstractVector
    options::NLoptOptions
    alg::NLoptAlg
    counter::Base.RefValue{Int}
end
function NLoptWorkspace(
    model::VecModel, optimizer::NLoptAlg,
    x0::AbstractVector = getinit(model);
    options = NLoptOptions(), kwargs...,
)
    problem, counter = get_nlopt_problem(optimizer.alg, model, copy(x0))
    return NLoptWorkspace(model, problem, copy(x0), options, optimizer, counter)
end
@params struct NLoptResult <: AbstractResult
    minimizer
    minimum
    problem
    status
    alg
    fcalls::Int
end

function optimize!(workspace::NLoptWorkspace)
    @unpack problem, options, x0, counter = workspace
    counter[] = 0
    foreach(keys(options.nt)) do k
        v = options.nt[k]
        setproperty!(problem, k, v)
    end
    minf, minx, ret = NLopt.optimize(problem, x0)
    return NLoptResult(
        copy(minx), minf, problem, ret, workspace.alg, counter[]
    )
end

function Workspace(model::VecModel, optimizer::NLoptAlg, args...; kwargs...,)
    return NLoptWorkspace(model, optimizer, args...; kwargs...)
end

function get_nlopt_problem(alg, model::VecModel, x0::AbstractVector)
    eq = if length(model.eq_constraints.fs) == 0
        nothing
    else
        model.eq_constraints
    end
    ineq = if length(model.ineq_constraints.fs) == 0
        nothing
    else
        model.ineq_constraints
    end
    obj = CountingFunction(getobjective(model))
    return get_nlopt_problem(
        alg,
        obj,
        ineq,
        eq,
        x0,
        getmin(model),
        getmax(model),
    ), obj.counter
end
function get_nlopt_problem(alg, obj, ineq_constr, eq_constr, x0, xlb, xub)
    onehot(n, i) = [zeros(i-1); 1.0; zeros(n-i)]
    local lastx, objval, objgrad
    local ineqconstrval, ineqconstrjac
    if ineq_constr !== nothing
        ineqconstrval, ineqconstrjac = nothing, nothing
    end
    local eqconstrval, eqconstrjac
    if eq_constr !== nothing
        eqconstrval, eqconstrjac = nothing, nothing
    end

    function update_cached_values(x, updategrad = true)
        lastx = copy(x)
        if updategrad
            objval, objpb = Zygote.pullback(obj, lastx)
            objgrad = objpb(1.0)[1]
        else
            objval = obj(lastx)
        end
        if ineq_constr !== nothing
            if updategrad
                ineqconstrval, ineqconstrpb = Zygote.pullback(
                    ineq_constr, lastx,
                )
                ineqconstrjac = mapreduce(vcat, 1:length(ineqconstrval)) do i
                    ineqconstrpb(onehot(length(ineqconstrval), i))[1]'
                end
            else
                ineqconstrval = ineq_constr(lastx)
            end
        end
        if eq_constr !== nothing
            if updategrad
                eqconstrval, eqconstrpb = Zygote.pullback(
                    eq_constr, lastx,
                )
                eqconstrjac = mapreduce(vcat, 1:length(eqconstrval)) do i
                    eqconstrpb(onehot(length(eqconstrval), i))[1]'
                end
            else
                eqconstrval = eq_constr(lastx)
            end
        end
    end
    update_cached_values(x0)

    function nlopt_obj(x, grad)
        if x == lastx
            if length(grad) > 0
                grad .= objgrad
            end
            return objval
        else
            update_cached_values(x, length(grad) > 0)
            if length(grad) > 0
                grad .= objgrad
            end
            return objval
        end
        return val
    end
    function nlopt_ineqconstr_gen(i)
        function (x, grad)
            if x == lastx
                if length(grad) > 0
                    grad .= ineqconstrjac[i, :]
                end
                return ineqconstrval[i]
            else
                update_cached_values(x, length(grad) > 0)
                if length(grad) > 0
                    grad .= ineqconstrjac[i, :]
                end
                return ineqconstrval[i]
            end
        end
    end
    function nlopt_eqconstr_gen(i)
        function (x, grad)
            if x == lastx
                if length(grad) > 0
                    grad .= eqconstrjac[i, :]
                end
                return eqconstrval[i]
            else
                update_cached_values(x, length(grad) > 0)
                if length(grad) > 0
                    grad .= eqconstrjac[i, :]
                end
                return eqconstrval[i]
            end
        end
    end

    problem = NLopt.Opt(alg, length(x0))
    problem.lower_bounds = xlb
    problem.upper_bounds = xub

    problem.min_objective = nlopt_obj
    if ineq_constr !== nothing
        for i in 1:length(ineqconstrval)
            NLopt.inequality_constraint!(
                problem, nlopt_ineqconstr_gen(i), 0.0,
            )
        end
    end
    if eq_constr !== nothing
        for i in 1:length(eqconstrval)
            NLopt.equality_constraint!(
                problem, nlopt_eqconstr_gen(i), 0.0,
            )
        end
    end
    return problem
end


text:
# # Getting Started
#
# This guide introduces the basic functionality of this package in a step by step manner.  This is a good starting point for learning about how to use this package.  For more details about how to use a particular model or coupling, refer to the documentation for that particular model or coupling.
#
#-
#md # !!! tip
#md #     This guide is also available as a Jupyter notebook:
#md #     [`guide.ipynb`](@__NBVIEWER_ROOT_URL__/examples/guide.ipynb).
#-
#
# ## Installing Aeroelasticity.jl
#
# If you haven't yet, now would be a good time to install this package.  It can be
# installed from the Julia REPL by typing `]` (to enter the package manager) and then
# running `add https://github.com/byuflowlab/Aeroelasticity.jl`.
#
# Now, that the package is installed we need to load it so that we can use it.

using Aeroelasticity
#!jl nothing #hide

# ## Assembling a Coupled Model
#
# For the purposes of this guide, we will be working with a two-degree-of-freedom typical
# section model, as shown in the following figure.
#
# ![](../assets/section-drawing.svg)
#
# Our goal is to create an 2D aeroelastic model which we can use to simulate the behavior
# of this system.  For the aerodynamic model, we will be using Peters' finite state model
# (see [`Peters`](@ref)).  For the structural model, we will be using a two degree of
# freedom typical section model (see [`Section`](@ref)).  To assemble an aeroelastic model
# using these two models we use the [`assemble_model`](@ref) function.

model = assemble_model(
    aerodynamic_model = Peters{4}(),
    structural_model = Section())

#!jl nothing #hide

# ## Defining Model Parameters
#
# Model parameters are defined using the [`assemble_parameters`](@ref) function.  In
# addition to the aerodynamic and structural parameters, we need to define a set of
# additional parameters for the model coupling.  A description of these parameters may
# be found in the documentation for each submodel and for the coupling.

## non-dimensional parameters
V = 1.0 # = U/(b*ωθ) (reduced velocity)
a = -1/5 # reference point normalized location
e = -1/10 # center of mass normalized location
μ = 20 # = m/(ρ*pi*b^2) (mass ratio)
r2 = 6/25 # = Iθ/(m*b^2) (radius of gyration about P)
σ = 2/5 # = ωh/ωθ (natural frequency ratio)
xθ = e - a # normalized distance from the reference point to the center of mass
a0 = 2*pi # lift curve slope
α0 = 0 # zero lift angle
cd0 = 0 # zero lift drag coefficient
cm0 = 0 # zero lift moment coefficient

## chosen dimensional parameters
b = 0.5 # semi-chord
ρ = 1 # air density
ωθ = 1 # pitch natural frequency
c = 343 # air speed of sound

## dimensionalized parameters
U = V*b*ωθ # velocity
m = μ*ρ*pi*b^2 # mass
Sθ = m*xθ*b # structural imbalance
Iθ = r2*m*b^2 # moment of inertia
ωh = σ*ωθ # plunge natural frequency
kh = m*ωh^2 # linear spring constant
kθ = Iθ*ωθ^2 # torsional spring constant

p = assemble_parameters(model;
    aerodynamic_parameters = (a=a, b=b, a0=a0, alpha0=α0, cd0=cd0, cm0=cm0),
    structural_parameters = (kh=kh, ktheta=kθ, m=m, Stheta=Sθ, Itheta=Iθ),
    additional_parameters = (U=U, rho=ρ, c=c)
)

#!jl nothing #hide

# ## Performing a Steady State Analysis
#
# We can find steady state solution for our model using
# [`DifferentialEquations`](https://github.com/SciML/DifferentialEquations.jl).

using DifferentialEquations

f = ODEFunction(model)

## initial guess for state variables
x0 = zeros(number_of_states(model))

## steady state problem
prob = SteadyStateProblem(f, x0, p)

## steady state solution
x_ss = solve(prob, SSRootfind())

#!jl nothing #hide

# Note that since our system is linear with respect to the state variables, the
# equilibrium point of our system will be at the origin.

# ## Performing a Stability Analysis
#
# The stability of a model for a given set of state variables, inputs, and parameters
# may be determined by calling the [`get_eigen`](@ref) function, which returns
# eigenvalues, left eigenvectors, and right eigenvectors.  For nonlinear systems, the
# provided state variables must correspond to an equilibrium point for the stability
# analysis to be theoretically valid.  Since our aeroelastic system is linear with
# respect to the state variables, any set of state variables will yield the same result.

K, M = linearize(model, x_ss, p)

λ, U, V = get_eigen(model, K, M)

#!jl nothing #hide

# A positive real part corresponding to any eigenvalue returned from the
# [`get_eigen`](@ref) function indicates that the system is unstable for the
# provided set of state variables, inputs, and parameters.

# ## Performing a Time Domain Simulation
#
# We can simulate the behavior of our system using
# [`DifferentialEquations.jl`](https://github.com/SciML/DifferentialEquations.jl).  Initial
# state variable values and rates may be specified using the [`assemble_states`](@ref)
# function.

using DifferentialEquations

## construct ODE function
f = ODEFunction(model)

## non-zero plunge degree of freedom
x0 = assemble_states(model;
    aerodynamic_states = (;lambda=zeros(4)),
    structural_states = (;h=0.5, theta=0, hdot=0, thetadot=0))

## simulate for 100 seconds
tspan = (0.0, 100.0)

## assemble problem
prob = DifferentialEquations.ODEProblem(f, x0, tspan, p)

## solve ODE problem
sol = DifferentialEquations.solve(prob)

#!jl nothing #hide

# We can use DifferentialEquations' built-in interface with the
# [Plots](https://github.com/JuliaPlots/Plots.jl) package to plot the simulation results.

using Plots
pyplot()

plot(sol,
    vars = [5,6,7,8],
    xlabel = "t",
    ylabel = permutedims([
        "\$h\$",
        "\$\\theta\$",
        "\$\\dot{h}\$",
        "\$\\dot{\\theta}\$",
        ]),
    label = "",
    layout = (4, 1),
    size = (600,1200)
    )

#md savefig("../assets/guide-solution.svg") #hide
#md nothing #hide

#md # ![](guide-solution.svg)

# We can also visualize the results of the typical section model using the
# [`section_coordinates`](@ref) function.

## animation parameters
a = -1/5
b = 0.5

## create animation
#md anim = @animate for t in range(tspan[1], tspan[2], length=200)
#!md @gif for t in range(tspan[1], tspan[2], length=200)
    h, θ = sol(t, idxs=7:8)
    xplot, yplot = section_coordinates(h, θ; a, b)
    plot(xplot, yplot;
        framestyle = :origin,
        grid = :false,
        xlims = (-1.0, 1.0),
        ylims = (-0.75, 0.75),
        aspect_ratio = 1.0,
        label = "t = $(round(t, digits=1))")
end

#md gif(anim, "../assets/guide-simulation.gif")
#md nothing #hide

#md # ![](guide-simulation.gif)

ext: jl


text:
using DSUtils
using Test

@testset "Data Science Utilities" begin
    @testset "ranks binning" begin
        x = [1, 1, 1, 2, 2, 2, 2, 2, 3, 3]
        xm = [1, 1, missing, 2, 2, 2, 3, 3, 3, 3]
        @test ranks(x) == [1, 1, 1, 5, 5, 5, 5, 5, 8, 8]
        @test ranks(x, rev=true) == [8, 8, 8, 4, 4, 4, 4, 4, 1, 1]
        @test ranks(x, rank=StatsBase.ordinalrank) == [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
        @test ranks(x, rank=StatsBase.competerank) == [0, 0, 0, 3, 3, 3, 3, 3, 8, 8]
        @test ranks(x, rank=StatsBase.denserank) == [0, 0, 0, 1, 1, 1, 1, 1, 2, 2]
        @test isequal(ranks(xm), [1, 1, missing, 3, 3, 3, 6, 6, 6, 6])
    end

    @testset "strstd string standardization" begin
        @test strstd("abc") == "abc"
        @test strstd("  abc def  ") == "abc_def"
        @test strstd("") == ""
        @test strstd("   ") == ""
        @test strstd(missing) == ""

        @test strstd("abc", empty="0") == "abc"
        @test strstd("  abc def  ", empty="0") == "abc_def"
        @test strstd("", empty="0") == "0"
        @test strstd("   ", empty="0") == "0"
        @test strstd(missing, empty="0") == "0"

        @test strstd("abc", empty=missing) == "abc"
        @test strstd("  abc def  ", empty=missing) == "abc_def"
        @test ismissing(strstd("", empty=missing))
        @test ismissing(strstd("   ", empty=missing))
        @test ismissing(strstd(missing, empty=missing))
    end

    @testset "onehot encoder" begin
        x = [1, 2, 3, 1, 2, 3, 4, 4, 4]
        xm = [1, 2, 3, 1, 2, 3, 4, 4, 4, missing, missing]
        s = ["a", "b", "c"]
        sm = ["a", "b", "c", missing]

        @test onehot(1, x) == [1, 0, 0, 1, 0, 0, 0, 0, 0]
        @test onehot(2, x) == [0, 1, 0, 0, 1, 0, 0, 0, 0]
        @test onehot(5, x) == [0, 0, 0, 0, 0, 0, 0, 0, 0]
        @test onehot(missing, x) == [0, 0, 0, 0, 0, 0, 0, 0, 0]

        @test onehot(1, xm) == [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]
        @test onehot(2, xm) == [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0]
        @test onehot(5, xm) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
        @test onehot(missing, xm) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]

        @test onehot("a", s) == [1, 0, 0]
        @test onehot("b", s) == [0, 1, 0]
        @test onehot(" ", s) == [0, 0, 0]
        @test onehot(missing, s) == [0, 0, 0]

        @test onehot("a", sm) == [1, 0, 0, 0]
        @test onehot("b", sm) == [0, 1, 0, 0]
        @test onehot(" ", sm) == [0, 0, 0, 0]
        @test onehot(missing, sm) == [0, 0, 0, 1]
    end

    @testset "sumxm sum with missings as 0" begin
        v1 = [1, 1, 1, 1, 1, 1]
        v2 = [2, 2, 2, 2, 2, 2]
        vm = [3, 3, 3, missing, missing, missing]

        @test sumxm(1, 2, 3) == 6
        @test sumxm(1, 2, 3, missing) == 6

        @test sumxm(v1, v2) == [3, 3, 3, 3, 3, 3]
        @test sumxm(v1, v2, vm) == [6, 6, 6, 3, 3, 3]
        @test sumxm(vm, vm) == [6, 6, 6, 0, 0, 0]
    end

    @testset "ks roc                      " begin
        t = [0, 0, 0, 1, 1, 1]
        x = [1, 1, 1, 2, 2, 2]

        @test kstest(t, x) == (n = 6, n1 = 3, n0 = 3, baserate = 0.5, ks = 1.0, ksarg = 2, ksdep = 0.5)
        @test auroc(t, x) == (conc = 9, tied = 0, disc = 0, auc = 1.0, gini = 1.0)
    end

end


text:
@testset "base" begin
    @testset "find_all_max" begin
        @test find_all_max([-Inf, -Inf, -Inf]) == (-Inf, [1, 2, 3])
        @test find_all_max([-Inf, -Inf, -Inf], [true, false, true]) == (-Inf, [1, 3])

        @test find_all_max([Inf, Inf, Inf]) == (Inf, [1, 2, 3])
        @test find_all_max([Inf, Inf, Inf], [true, true, false]) == (Inf, [1, 2])

        @test find_all_max([Inf, 0, Inf]) == (Inf, [1, 3])
        @test find_all_max([Inf, 0, Inf], [false, true, false]) == (0, [2])

        @test find_all_max([0, 1, 2, 1, 2, 1, 0]) == (2, [3, 5])
        @test find_all_max([0, 1, 2, 1, 2, 1, 0], Bool[1, 1, 0, 0, 0, 1, 1]) == (1, [2, 6])
    end

    @testset "flatten_batch" begin
        x = rand(2, 3, 4)
        y = flatten_batch(x)
        @test size(y) == (2, 12)
    end

    @testset "discount_rewards" begin
        reward, γ = [1.0], 0.5
        @test discount_rewards(reward, γ) ≈ reward
        @test discount_rewards(reward, γ; init = 2.0) ≈ [2.0]

        reward, γ = Float64[1, 2, 3], 0.5

        @test discount_rewards(reward, γ) ≈ [2.75, 3.5, 3.0]
        @test discount_rewards(reward, γ; init = 4.0) ≈ [3.25, 4.5, 5.0]

        @test discount_rewards(reward, γ; terminal = [false true false], init = 2.0) ≈
              [2.0, 2.0, 4.0]
        @test discount_rewards(reward, γ; terminal = [true false true], init = 2.0) ≈
              [1.0, 3.5, 3.0]

        # type stable
        reward, γ = [1, 2, 3], 0.5f0
        @test eltype(discount_rewards(reward, γ)) == Float32

        # 2D
        reward, γ = reshape(1:9, 3, 3), 0.5
        init = [-2.0, 0.0, 2.0]
        # for 2d rewards, the keyword argument of `dim` must be either `1` or `2`
        @test_throws MethodError discount_rewards(reward, γ)
        @test discount_rewards(reward, γ; dims = 1) ≈
              [2.75 8.0 13.25; 3.5 8.0 12.5; 3.0 6.0 9.0]
        @test discount_rewards(reward, γ; dims = 2) ≈
              [4.75 7.5 7.0; 6.5 9.0 8.0; 8.25 10.5 9.0]
        @test discount_rewards(reward, γ; init = init, dims = 1) ≈
              [2.5 8.0 13.5; 3.0 8.0 13.0; 2.0 6.0 10.0]
        @test discount_rewards(reward, γ; init = init, dims = 2) ≈
              [4.5 7.0 6.0; 6.5 9.0 8.0; 8.5 11.0 10]

        terminal = [false true false; true false true; false true false]
        @test discount_rewards(reward, γ; dims = 1, terminal = terminal) ≈
              [2.0 4.0 11.0; 2.0 8.0 8.0; 3.0 6.0 9.0]
        @test discount_rewards(reward, γ; dims = 1, terminal = terminal, init = init) ≈
              [2.0 4.0 11.0; 2.0 8.0 8.0; 2.0 6.0 10.0]
        @test discount_rewards(reward, γ; dims = 2, terminal = terminal, init = init) ≈
              [3.0 4.0 6.0; 2.0 9.0 8.0; 6.0 6.0 10.0]
    end

    @testset "discount_rewards_reduced" begin
        reward, γ = [1.0], 0.5
        @test discount_rewards_reduced(reward, γ) ≈ 1.0

        reward, γ = Float64[1, 2, 3], 0.5
        @test discount_rewards_reduced(reward, γ) ≈ 2.75
        @test discount_rewards_reduced(reward, γ; init = 4.0) ≈ 3.25
        @test discount_rewards_reduced(reward, γ; terminal = [false, true, false]) ≈ 2.0
        @test discount_rewards_reduced(
            reward,
            γ;
            terminal = [false, true, false],
            init = 4.0,
        ) ≈ 2.0

        reward, γ = reshape(1:9, 3, 3), 0.5
        init = [-2.0, 0.0, 2.0]
        terminal = [false true false; true false true; false true false]

        # for reward of 2D, `dims` must be provided
        @test_throws Exception discount_rewards_reduced(reward, γ)

        @test discount_rewards_reduced(reward, γ; dims = 1) ≈ [2.75, 8.0, 13.25]
        @test discount_rewards_reduced(reward, γ; dims = 2) ≈ [4.75, 6.5, 8.25]
        @test discount_rewards_reduced(
            reward,
            γ;
            dims = 1,
            terminal = terminal,
            init = init,
        ) ≈ [2.0, 4.0, 11.0]
        @test discount_rewards_reduced(
            reward,
            γ;
            dims = 2,
            terminal = terminal,
            init = init,
        ) ≈ [3.0, 2.0, 6.0]
    end

    @testset "generalized_advantage_estimation" begin
        reward, values, γ, λ = [1.0], [2.0, 3.0], 0.5, 0.3
        @test generalized_advantage_estimation(reward, values, γ, λ) ≈ [0.5]

        reward, values, γ, λ = [1.0, 1.0], [1, 2, 3], 0.5, 0.3
        @test generalized_advantage_estimation(reward, values, γ, λ) ≈ [1.075, 0.5]

        reward, values, γ, λ = Float64[1, 2, 3], [1, 2, 3, 4], 0.5, 0.3
        @test generalized_advantage_estimation(reward, values, γ, λ) ≈ [1.27, 1.8, 2]

        @test generalized_advantage_estimation(
            reward,
            values,
            γ,
            λ;
            terminal = [true false true],
        ) ≈ [0.0, 1.5, 0.0]

        # type stable
        reward, values, γ, λ = [1, 2, 3], [1, 2, 3, 4], 0.5f0, 0.5f0
        @test eltype(generalized_advantage_estimation(reward, values, γ, λ)) == Float32

        # 2D
        reward, values, γ, λ = reshape(1:9, 3, 3), reshape(1:12, 4, 3), 0.5, 0.3
        # for 2d rewards, the keyword argument of `dim` must be either `1` or `2`
        @test_throws MethodError generalized_advantage_estimation(reward, values, γ, λ)
        @test generalized_advantage_estimation(reward, values, γ, λ; dims = 1) ≈
              [1.27 2.4425 3.615; 1.8 2.95 4.1; 2.0 3.0 4.0]

        values = reshape(1:12, 3, 4)
        @test generalized_advantage_estimation(reward, values, γ, λ; dims = 2) ≈
              [2.6375 4.25 5.0; 3.22375 4.825 5.5; 3.81 5.4 6.0]

        reward, values, γ, λ = reshape(1:9, 3, 3), reshape(1:12, 4, 3), 0.5, 0.3
        terminal = [false true false; true false true; false true false]

        @test generalized_advantage_estimation(
            reward,
            values,
            γ,
            λ;
            dims = 1,
            terminal = terminal,
        ) ≈ [1.0 -1.0 2.7; 0.0 2.35 -2.0; 2.0 -1.0 4.0]

        values = reshape(1:12, 3, 4)
        @test generalized_advantage_estimation(reward, values, γ, λ; dims = 2) ≈
              [2.6375 4.25 5.0; 3.22375 4.825 5.5; 3.81 5.4 6.0]
    end

end


text:
using LinearAlgebra

"""
	sim_activity(J::Matrix{T}, g::T, h0::Vector{T}; dt=0.1, t_bound=1000.0,
		s=0.0, b=nothing, savestep=1, minimum_norm=0) where T <: AbstractFloat

Integrates a tanh rate network with coupling matrix `J` and strength scaling `g`
from initial condition `h0`, returning a `N x T` matrix of network unit activity.

When `minimum_norm > 0` the integration will be cut off when the network state
falls too close to the origin to prevent floating point underflow. The keyword
arguments `s` and `b` add self-coupling to `J` and add a bias term, respectively.
"""
function sim_activity(J::Matrix{T}, g::T, h0::Vector{T}; dt=0.1, t_bound=1000.0,
		s=0.0, b=nothing, savestep=1, minimum_norm=0) where T <: AbstractFloat

	N = size(J)[1]
	# Check network size and type consistency
	@assert N == size(J)[2] == length(h0)
	@assert T == typeof(dt) == typeof(t_bound) == typeof(s)

	if b isa Array{T} && (size(b) == (N,) || size(b) == (N, 1))
		use_bias = true
	else
		use_bias = false
	end

	maxstep = round(Int, t_bound / dt) # Discretized simulation length

	# Allocate memory for output
	out = Matrix{T}(undef, N, maxstep ÷ savestep)

	h = copy(h0) # Set initial condition
	r = Vector{T}(undef, N)

	J = g .* J + s * I

	step = 0
	while step < maxstep
		step += 1
		map!(tanh, r, h)
		BLAS.gemv!('N', dt, J, r, 1 - dt, h)
		if use_bias
			h .+= dt .* b
		end
		if step % savestep == 0 # Save every specified time interval
			savept = step ÷ savestep
			out[:,savept] .= h
		end
		if norm(h) < minimum_norm
			D = -I + J .* ((sech.(h)).^2)'
			if all(real(eigvals(D)) .< 0)
				break
			end
		end
	end

	lastindex = step ÷ savestep
	return out[:,1:lastindex]
end

"""
	sim_lyapunov(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T},
		Q0::Matrix{T}; s=0.0, b=nothing, dt=0.1, t_bound=1000.0, t_ONS=1.0,
		wait=0.0, return_max_LE=false, printstuff=false, scale_limit=100,
		p_factor=1e-1) where T <: AbstractFloat

Integrates a tanh network with a given weight matrix and gain parameter from a
specified initial condition, and returns a Lyapunov spectrum calculated along
the trajectory.

Optionally returns the maximum Lyapunov exponent calculated using a direct
numerical simulation of a perturbed trajectory.
"""
function sim_lyapunov(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T},
		Q0::Matrix{T}; s=0.0, b=nothing, dt=0.1, t_bound=1000.0, t_ONS=1.0,
		wait=0.0, return_max_LE=false, printstuff=false, scale_limit=100,
		p_factor=1e-1) where T <: AbstractFloat

	N = size(J)[1]

	@assert(size(J) == (length(h0), length(dh)) == size(Q0) == (N,N),
			"Dimensions of the coupling matrix, IC, initial ONS, and " *
			"perturbation must be consistent"
	)
	@assert(T == typeof(s) == typeof(dt) == typeof(t_bound) == typeof(t_ONS)
			== typeof(wait), "Make sure all scalar parameters are the same type"
	)

	maxstep = round(Int64, t_bound / dt)
	startstep = round(Int64, wait / dt)
	tONSstep = round(t_ONS / dt)
	@assert(maxstep % tONSstep == startstep % tONSstep == 0,
		"t_ONS must exactly divide t_bound and wait"
	)

	# Whether or not to include a bias term in the network dynamics
	if b isa Array{T} && size(b) == (N,)
		use_bias = true
	elseif b isa Array{T} && size(b) == (N, 1)
		use_bias = true
		b = reshape(b, (N,))
	else
		use_bias = false
	end

	# Allocate memory for output
	spectrum = zeros(T, N)

	# Allocate vector for activity and firing rate of current state
	h = copy(h0)
	r = Vector{T}(undef, N)

	# Allocate memory for Jacobian calculation and orthonormal basis
	sech2(x::T)::T = (sech(x))^2
	sech2h = Matrix{T}(undef, 1, N)
	D = Matrix{T}(undef, N, N)
	Q = copy(Q0)
	dQ = Matrix{T}(undef, N, N)

	if return_max_LE
		max_LE = 0
		p_factor = convert(T, p_factor)
		dh = copy(dh) ./ norm(dh)
		h2 = h .+ p_factor .* dh
		r2 = Vector{T}(undef, N)
	end

	# The effective coupling matrix J
	J = g .* J + s * I

	step = 0
	while step < maxstep
		step += 1

		#Evolve the trajectory and the orthonormal basis in the tangent space
		map!(tanh, r, h)
		BLAS.gemv!('N', dt, J, r, 1 - dt, h)
		if use_bias
			h .+= dt .* b
		end
		map!(sech2, sech2h, h)
		copyto!(D, J)
		D .*= sech2h
		D[diagind(D)] .-= 1
		BLAS.gemm!('N', 'N', dt, D, Q, T(0), dQ)
		Q .+= dQ

		#Evolve the perturbed trajectory and calculate how much it has grown
		if return_max_LE
			map!(tanh, r2, h2)
			BLAS.gemv!('N', dt, J, r2, 1 - dt, h2)
			if use_bias
				h2 .+= dt .* b
			end
			dh .= (h2 .- h) ./ p_factor
			scale = norm(dh)

			# If the perturbation scale has grown or shrunk too much, then
			# reset it to 1 and track the cumulative sum of log scales
			if scale > scale_limit || scale < 1 / scale_limit || step == startstep
				h2 .= h .+ dh ./ scale .* p_factor
				if step > startstep
					max_LE += log(scale)
				end
			end
		end

		# At each re-orthonormalization step, the log diagonal entries of R are
		# saved in a cumulative sum
		if step % tONSstep == 0
			printstuff && println("Time step: ", step)
			copyto!(dQ, Q)
			F = qr!(dQ)
			diags = view(F.factors, diagind(F.factors))
			if step > startstep  # forcing positive R-diagonal elements
				spectrum .+= log.(abs.(diags))
			end
			fill!(Q, 0)
			map!(sign, view(Q, diagind(Q)), diags)
			lmul!(F.Q, Q)
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	true_time = dt * (maxstep - startstep)
	if return_max_LE
		max_LE += log(norm((h2 .- h) ./ p_factor))
		return spectrum ./ true_time, max_LE / true_time
	else
		return spectrum ./ true_time
	end
end

"""
	sim_direct_max_LE(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T};
		s=0.0, b=nothing, dt=0.1, t_bound=1000.0, wait=0.0, printstuff=false,
		scale_limit=100, p_factor=1e-1) where T <: AbstractFloat

Calculates the maximum Lyapunov exponent of a tanh rate network given by
coupling matrix J, dependent on a given initial condition and initial
perturbation, using a direct method that tracks the growth of the perturbed
trajectory relative to the unperturbed one.
"""
function sim_direct_max_LE(J::Matrix{T}, g::T, h0::Vector{T}, dh::Vector{T};
		s=0.0, b=nothing, dt=0.1, t_bound=1000.0, wait=0.0, printstuff=false,
		scale_limit=100, p_factor=1e-1) where T <: AbstractFloat

	N = size(J)[1]
	@assert(size(J) == (length(h0), length(dh)) == (N,N), "Dimensions of the " *
			"coupling matrix, IC, and perturbation must be consistent"
	)

	# Whether or not to include a bias term in the network dynamics
	if b isa Array{T} && size(b) == (N,)
		use_bias = true
	elseif b isa Array{T} && size(b) == (N, 1)
		use_bias = true
		b = reshape(b, (N,))
	else
		use_bias = false
	end

	maxstep = round(Int64, t_bound / dt)
	startstep = round(Int64, wait / dt)

	# Effective coupling matrix
	J = g .* J + s * I

	# Order of the perturbation handled by the sim
	p_factor = convert(T, p_factor)
	# Allocate memory for output
	max_LE = 0

	# Set the initial condition and perturbed IC
	h = copy(h0)
	dh = copy(dh) ./ norm(dh)
	h2 = h .+ p_factor .* dh

	r = Vector{T}(undef, N)
	r2 = Vector{T}(undef, N)

	step = 0
	while step < maxstep
		step += 1

		# Evolve the trajectory and its perturbation
		map!(tanh, r, h)
		map!(tanh, r2, h2)
		BLAS.gemv!('N', dt, J, r, 1 - dt, h)
		BLAS.gemv!('N', dt, J, r2, 1 - dt, h2)
		if use_bias
			h .+= dt .* b
			h2 .+= dt .* b
		end
		dh .= (h2 .- h) ./ p_factor
		scale = norm(dh)

		# If the perturbation scale has grown or shrunk too much, then
		# reset it to 1 and track the cumulative sum of log scales
		if scale > scale_limit || scale < 1 / scale_limit || step == startstep
			printstuff && println("Time step: ", step, ", Scale: ", scale)
			h2 .= h .+ dh ./ scale .* p_factor
			if step > startstep
				max_LE += log(scale)
			end
		end
	end

	# Update the maximum Lyapunov exponent after the simulation has ended so
	# that the full simulation time is used for the calculation
	max_LE += log(norm((h2 .- h) ./ p_factor))
	true_time = dt * (maxstep - startstep)
	return max_LE / true_time
end


text:
using CausticsEngineering
using Test

"""
$(SIGNATURES)

"""
function testSquareMesh!()
    mesh = create_mesh(100, 50)

    println(mesh.nodeArray[1, 1])
    println(mesh.nodes[1])

    mesh.nodeArray[1, 1].x = 8
    println(mesh.nodeArray[1, 1])
    println(mesh.nodes[1])

    mesh.nodes[1].y += 12
    println(mesh.nodeArray[1, 1])
    println(mesh.nodes[1])
end

"""
$(SIGNATURES)

TO REFACTOR.
"""
function testSolidify!()
    println("Testing solidification")
    width = 100
    height = 100
    origMesh = create_mesh(width, height)

    for y = 1:height, x = 1:width
        x2 = (x - width / 2) / width
        y2 = (y - height / 2) / height
        value = x2 * x2 + y2 * y2
        origMesh.nodeArray[x, y].z = 15 - value * 25
    end

    save_stl!(origMesh, "./examples/testSolidify.obj")
    solidMesh = solidify(origMesh, 0)
    save_stl!(solidMesh, "./examples/testSolidify2.obj")
end

function testAreaTriangle(x1, y1, x2, y2, x3, y3)
    a = sqrt((x2 - x1)^2 + (y2 - y1)^2)
    b = sqrt((x3 - x2)^2 + (y3 - y2)^2)
    c = sqrt((x1 - x3)^2 + (y1 - y3)^2)
    s = (a + b + c) / 2.0

    surface_sq = (a + b + c) * (-a + b + c) * (a - b + c) * (a + b - c) / 16.0
    return surface_sq <= 1e-100 ? 0.0 : surface_sq
end

function doTestAreaTriangle(N = 10_000)
    X1 = 10_000 .* rand(N)
    Y1 = 10_000 .* rand(N)
    X2 = 10_000 .* rand(N)
    Y2 = 10_000 .* rand(N)
    X3 = 10_000 .* rand(N)
    Y3 = 10_000 .* rand(N)

    for x1 ∈ X1, y1 ∈ Y1, x2 ∈ X2, y2 ∈ Y2, x3 ∈ X3, y3 ∈ Y3

        _ = testAreaTriangle(x1, y1, x2, y2, x3, y3) > 0.0
    end
end

@testset "triangle area" begin
    for x1 ∈ 10_000 .* rand(10_000),
        y1 ∈ 10_000 .* rand(10_000),
        x2 ∈ 10_000 .* rand(10_000),
        y2 ∈ 10_000 .* rand(10_000),
        x3 ∈ 10_000 .* rand(10_000),
        y3 ∈ 10_000 .* rand(10_000)

        @test testAreaTriangle(x1, y1, x2, y2, x3, y3) >= 0.0
    end
end

@testset "More triangle area" begin
    @test area(0, 1, 1, 0, 0, 0)
end

v1 = CausticsEngineering.Vertex3D(0.0, 1.0, 0.0, 0.0, 0.0)
v2 = CausticsEngineering.Vertex3D(1.0, 0.0, 0.0, 0.0, 0.0)
v3 = CausticsEngineering.Vertex3D(0.0, 0.0, 0.0, 0.0, 0.0)
CausticsEngineering.area(v1, v2, v3)

v1 = CausticsEngineering.Vertex3D(0.0, 2.0, 0.0, 0.0, 0.0)
v2 = CausticsEngineering.Vertex3D(2.0, 0.0, 0.0, 0.0, 0.0)
v3 = CausticsEngineering.Vertex3D(0.0, 0.0, 0.0, 0.0, 0.0)
CausticsEngineering.area(v1, v2, v3)

fm = CausticsEngineering.FaceMesh(512, 512);
fm.toptriangles[1, 1]
CausticsEngineering.triangle3D(fm, 1, 1, :top)
CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 1, 1, :top)...)
CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 1, 1, :bottom)...)

CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 126, 13, :top)...)
CausticsEngineering.area(CausticsEngineering.triangle3D(fm, 126, 13, :bottom)...)

v1 = [7250.768146661309, 4849.994680605491]
v2 = [1725.2665442814496, 8169.783790398901]
v3 = [6560.1725812433, 5264.912964907599]

a = sqrt(sum((v2 - v1) .^ 2))
b = sqrt(sum((v3 - v2) .^ 2))
c = sqrt(sum((v1 - v3) .^ 2))
s = (a + b + c) / 2.0

s * (s - a) * (s - b) * (s - c)
(a + b + c) * (-a + b + c) * (a - b + c) * (a + b - c) / 16.0

b + c
a

10^8 * (b + c - a)


text:
module Splines2

export ns, ns_, bs, bs_, is, is_, ms, ms_, basis

using OffsetArrays
using LinearAlgebra
using Statistics

# utility functions
zeroArray(a :: Array{T}) where T<:Real = OffsetArray(a,(size(a).*0).-1)

# type definitions
abstract type AbstractSplineBasis{T<:Real} end
# Note: we have used OffsetArray for converting from C code
mutable struct SplineBasis{T<:Real} <: AbstractSplineBasis{T}
    order::Int # order of the spline
    nknots::Int # number of knots
    ncoef::Int # number of coefficients
    knots::OffsetArray{T,1} # knot vector
end
mutable struct BSplineBasis{T<:Real} <: AbstractSplineBasis{T}
    spline_basis::SplineBasis{T}
    boundary_knots::Tuple{T,T}
    interior_knots::Union{Array{T,1}, Nothing}
    intercept::Bool
    df::Int
end
mutable struct NSplineBasis{T<:Real} <: AbstractSplineBasis{T}
    b_spline_basis::BSplineBasis{T}
    qmat::Array{T,2}
    tl0::Array{T,1}
    tl1::Array{T,1}
    tr0::Array{T,1}
    tr1::Array{T,1}
end

# constructors
function SplineBasis(knots :: Array{T,1}, order :: Int = 4) where T<:Real
    SplineBasis(order,
                length(knots),
                length(knots)-order,
                zeroArray(knots))
end
function BSplineBasis(boundary_knots :: Tuple{T,T},
                      interior_knots :: Union{Array{T,1}, Nothing} = nothing,
                      order :: Int = 4,
                      intercept :: Bool = false) where T<:Real
    l_interior_knots = interior_knots == nothing ? 0 : length(interior_knots)
    df = Int(intercept) + order - 1 + l_interior_knots
    nknots = l_interior_knots + 2*order
    ncoef = nknots - order
    knots = zeros(T,nknots)
    for i=1:order
        knots[i] = boundary_knots[1]
        knots[nknots-i+1] = boundary_knots[2]
    end
    if (l_interior_knots>0)
        for i=1:l_interior_knots
            knots[i+order]=interior_knots[i]
        end
    end
    BSplineBasis(SplineBasis(knots,order), boundary_knots, interior_knots, intercept, df)
end
function NSplineBasis(boundary_knots :: Tuple{T,T},
                      interior_knots :: Union{Array{T,1}, Nothing} = nothing,
                      order :: Int = 4,
                      intercept :: Bool = false) :: NSplineBasis{T} where T<:Real
    bs = BSplineBasis(boundary_knots,interior_knots,order,intercept)
    co = basis(bs,[bs.boundary_knots[1],bs.boundary_knots[2]],2)
    qmat_ = LinearAlgebra.qr(transpose(co)).Q
    qmat = collect(transpose(qmat_[:, 3:size(qmat_,2)]))
    tl0 = qmat * basis(bs,boundary_knots[1])
    tl1 = qmat * basis(bs,boundary_knots[1],1)
    tr0 = qmat * basis(bs,boundary_knots[2])
    tr1 = qmat * basis(bs,boundary_knots[2],1)
    NSplineBasis(bs,
                 qmat,
                 tl0,tl1,tr0,tr1)
end

function basis(bs :: SplineBasis{T}, x :: T, ders :: Int = 0) where T<:Real
    function find_interval(x :: T)
        k = bs.order - 1
        n = bs.nknots - k - 1
        l = k
        while (x < bs.knots[l] && l != k)
            l = l - 1
        end
        l = l + 1
        while (x >= bs.knots[l] && l != n)
            l = l + 1
        end
        l - 1
    end
    t = bs.knots
    k = bs.order - 1
    m = ders
    hh = bs.order
    ell = find_interval(x)
    result = zeroArray(zeros(T,2*k+2))
    one = T(1)
    zero = T(0)
    result[0] = one
    for j = 1:(k-m)
        for n = 0:(j-1)
            result[hh+n] = result[n]
        end
        result[0] = zero
        for n = 1:j
            ind = ell + n
            xb = t[ind]
            xa = t[ind - j]
            if (xb == xa)
                result[n] = zero
                continue
            end
            w = result[hh + n - 1]/(xb - xa)
            result[n - 1] = result[n-1] + w*(xb - x)
            result[n] = w*(x - xa)
        end
    end
    for j = (k - m + 1):k
        for n = 0:(j-1)
            result[hh+n] = result[n]
        end
        result[0] = zero
        for n = 1:j
            ind = ell + n
            xb = t[ind]
            xa = t[ind - j]
            if (xb == xa)
                result[m] = zero
                continue
            end
            w = j*result[hh + n - 1]/(xb - xa)
            result[n - 1] = result[n - 1] - w
            result[n] = w
        end
    end
    offset = ell - k
    v = zeros(T,bs.ncoef)
    v[(1+offset):(k+1+offset)] = result[0:k]
    v
end

function basis(bs :: BSplineBasis{T}, x :: T, ders :: Int = 0) where T<:Real
    if (x < bs.boundary_knots[1] || x > bs.boundary_knots[2])
        if (x < bs.boundary_knots[1])
            k_pivot = T(0.75)*bs.boundary_knots[1]+
            T(0.25)*bs.spline_basis.knots[5-1] # 0-based
        else
            k_pivot = T(0.75)*bs.boundary_knots[2]+
            T(0.25)*bs.spline_basis.knots[length(bs.spline_basis.knots)-4-1] # 0-based
        end
        delta = x - k_pivot
        if (ders==0)
            vec = basis(bs.spline_basis,k_pivot,0) +
                basis(bs.spline_basis,k_pivot,1)*delta +
                basis(bs.spline_basis,k_pivot,2)*delta*delta/T(2.0) +
                basis(bs.spline_basis,k_pivot,3)*delta*delta*delta/T(6.0)
        elseif (ders==1)
            vec = splines.basis(bs.spline_basis,k_pivot,1) +
                splines.basis(bs.spline_basis,k_pivot,2)*delta +
                splines.basis(bs.spline_basis,k_pivot,3)*delta*delta/T(2.0)
        elseif (ders==2)
            vec = splines.basis(bs.spline_basis,k_pivot,2) +
                splines.basis(bs.spline_basis,k_pivot,3)*delta
        elseif (ders==3)
            vec = splines.basis(bs.spline_basis,k_pivot,3)
        else
            vec = k_pivot .* T(0)
        end
    else
        vec = basis(bs.spline_basis, x, ders)
    end
    if (!bs.intercept)
        vec = vec[2:length(vec)]
    end
    vec
end

function basis(ns :: NSplineBasis{T}, x::T, ders :: Int = 0) where T<:Real
    if (x < ns.b_spline_basis.boundary_knots[1])
        if (ders==0)
            vec = ns.tl0 + (x - ns.b_spline_basis.boundary_knots[1])*ns.tl1
        elseif (ders==1)
            vec = ns.tl1
        else
            vec = ns.tl1.* T(0)
        end
    elseif (x > ns.b_spline_basis.boundary_knots[2])
        if (ders==0)
            vec = ns.tr0 + (x - ns.b_spline_basis.boundary_knots[2])*ns.tr1
        elseif (ders==1)
            vsc = ns.tr1
        else
            vec = ns.tr1.* T(0)
        end
    else
        vec = ns.qmat*basis(ns.b_spline_basis, x, ders)
    end
    vec
end

function basis(bs :: AbstractSplineBasis{T}, x :: Array{T,1}, ders :: Int = 0) where T<:Real
    f(xi) = basis(bs, xi, ders)
    copy(transpose(reduce(hcat, f.(x))))
end

# utility function for processing the spline arguments
function spline_args(x :: Array{T,1};
                     boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
                     interior_knots :: Union{Array{T,1},Nothing} = nothing,
                     order :: Int = 4,
                     intercept :: Bool = false,
                     df :: Int = 3 + Int(intercept),
                     knots :: Union{Array{T,1}, Nothing} = nothing,
                     knots_offset :: Int = 0) where T<:Real
    if (interior_knots != nothing && boundary_knots != nothing)
        # pass
    elseif (knots != nothing)
        boundary_knots = extrema(knots)
        interior_knots = length(knots)==2 ? nothing : knots[2:(length(knots)-1)]
    else
        if (boundary_knots == nothing)
            boundary_knots = extrema(x)
        end
        iKnots = df - order + knots_offset + 1 - Int(intercept)
        if (iKnots>0)
            p = range(T(0), length=iKnots+2, stop=T(1))[2:(iKnots+1)]
            index = (x .>= boundary_knots[1]) .* (x .<= boundary_knots[2])
            interior_knots = Statistics.quantile(x[index], p)
        end
    end
    (boundary_knots, interior_knots)
end

"""
    bs_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for B-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

# Examples
```jldoctest
julia> Splines2.bs_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 0.384  0.096  0.008
 0.432  0.288  0.064
 0.288  0.432  0.216
 0.096  0.384  0.512
 0.0    0.0    1.0
```
"""
function bs_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int = 4,
            intercept :: Bool = false,
            df :: Int = order - 1 + Int(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing,
            centre :: Union{T,Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, interior_knots=interior_knots,
                    order=order, intercept=intercept,
                    df=df, knots=knots)
    spline = BSplineBasis(boundary_knots, interior_knots, order, intercept)
    spline
end

"""
    bs(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for B-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.bs(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 0.384  0.096  0.008
 0.432  0.288  0.064
 0.288  0.432  0.216
 0.096  0.384  0.512
 0.0    0.0    1.0
```
"""
function bs(x :: Array{T,1}; ders :: Int = 0, kwargs...) where T<:Real
    bs_(x; kwargs...)(x, ders=ders)
end

"""
    ns_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for natural B-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

# Examples
```jldoctest
julia> Splines2.ns_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
  0.0       0.0        0.0
 -0.100444  0.409332  -0.272888
  0.102383  0.540852  -0.359235
  0.501759  0.386722  -0.172481
  0.418872  0.327383   0.217745
 -0.142857  0.428571   0.714286
```
"""
function ns_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int = 4,
            intercept :: Bool = false,
            df :: Int = order - 3 + Int(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing,
            centre :: Union{T,Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, interior_knots=interior_knots,
                    order=order, intercept=intercept, df=df, knots=knots, knots_offset=2)
    spline = NSplineBasis(boundary_knots, interior_knots, order, intercept)
    spline
end

function (spline::Union{BSplineBasis, NSplineBasis})(x; ders :: Int = 0)
    b = basis(spline, x, ders)
    #if (centre != nothing && ders==0)
    #    bc = basis(spline, centre, ders)
    #    for i=1:size(b,1)
    #        b[i,:] -= bc
    #    end
    #end
    b
end

"""
    ns(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for natural B-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.ns(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
  0.0       0.0        0.0
 -0.100444  0.409332  -0.272888
  0.102383  0.540852  -0.359235
  0.501759  0.386722  -0.172481
  0.418872  0.327383   0.217745
 -0.142857  0.428571   0.714286
```
"""
function ns(x :: Array{T,1}; ders :: Int = 0, kwargs...) where T<:Real
    ns_(x; kwargs...)(x, ders=ders)
end

"""
    is_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for I-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)

# Examples
```jldoctest
julia> Splines2.is_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
 0.0     0.0     0.0
 0.1808  0.0272  0.0016
 0.5248  0.1792  0.0256
 0.8208  0.4752  0.1296
 0.9728  0.8192  0.4096
 1.0     1.0     1.0
```
"""
function is_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int = 4,
            intercept :: Bool = false,
            df :: Int = order - 1 + Int(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, order=order+1,
                    intercept=intercept, df=df+1, knots=knots)
    spline = BSplineBasis(boundary_knots, interior_knots, order+1, false)
    knots = parent(spline.spline_basis.knots)
    findj(x) = interior_knots==nothing ? (order+1) : searchsortedlast(knots,x)
    function eval(x :: Array{T,1}; ders :: Int = 0)
        b = basis(spline,x,ders)
        (nrow,ncol) = size(b)
        for i=1:nrow
            js = findj(x[i])
            for j in ncol:-1:1
                if j>js
                    b[i,j]=T(0)
                elseif j==ncol
                    b[i,j]=b[i,j]
                else
                    b[i,j]=b[i,j]+b[i,j+1]
                end
            end
            for j in ncol:-1:1
                if j < js-(order+1)
                    b[i,j]=T(1)
                end
            end
        end
        intercept ? b : b[:, 2:size(b,2)]
    end
    eval
end

"""
    is(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for I-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `ders :: Int = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.is(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
 0.0     0.0     0.0
 0.1808  0.0272  0.0016
 0.5248  0.1792  0.0256
 0.8208  0.4752  0.1296
 0.9728  0.8192  0.4096
 1.0     1.0     1.0
```
"""
function is(x :: Array{T,1}; ders :: Int = 0, kwargs...) where T<:Real
    is_(x; kwargs...)(x, ders=ders)
end

"""
    ms_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for M-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

# Examples
```jldoctest
julia> Splines2.ms_(collect(0.0:0.2:1.0), df=3)(collect(0.0:0.2:1.0))
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 1.536  0.384  0.032
 1.728  1.152  0.256
 1.152  1.728  0.864
 0.384  1.536  2.048
 0.0    0.0    4.0
```
"""
function ms_(x :: Array{T,1};
            boundary_knots :: Union{Tuple{T,T},Nothing} = nothing,
            interior_knots :: Union{Array{T,1},Nothing} = nothing,
            order :: Int = 4,
            intercept :: Bool = false,
            df :: Int = order - 1 + Int(intercept),
            knots :: Union{Array{T,1}, Nothing} = nothing,
            centre :: Union{T,Nothing} = nothing) where T<:Real
    (boundary_knots, interior_knots) =
        spline_args(x, boundary_knots=boundary_knots, interior_knots=interior_knots,
                    order=order, intercept=intercept, df=df, knots=knots)
    spline = BSplineBasis(boundary_knots, interior_knots, order, intercept)
    function eval(x :: Array{T,1}; ders :: Int = 0)
        b = basis(spline, x, ders)
        knots = parent(spline.spline_basis.knots)
        function loop(j)
            denom = knots[j+order]-knots[j]
            denom > T(0) ? order/denom : T(0)
        end
        transCoef = map(loop,1:(length(knots)-order))
        if (!intercept)
            transCoef = transCoef[2:length(transCoef)]
        end
        if (centre != nothing && ders==0)
            bc = basis(spline, centre, ders)
            for i=1:size(b,1)
                b[i,:] -= bc
            end
        end
        for i=1:size(b,1)
            b[i,:] .*= transCoef
        end
        b
    end
    eval
end

"""
    ms(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for M-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int = 0`: derivatives of the splines

# Examples
```jldoctest
julia> Splines2.ms(collect(0.0:0.2:1.0), df=3)
6×3 Array{Float64,2}:
 0.0    0.0    0.0
 1.536  0.384  0.032
 1.728  1.152  0.256
 1.152  1.728  0.864
 0.384  1.536  2.048
 0.0    0.0    4.0
```
"""
function ms(x :: Array{T,1}; ders :: Int = 0, kwargs...) where T<:Real
    ms_(x; kwargs...)(x, ders=ders)
end

end # module


text:
# It's not a model --- its a simulation!

mutable struct Simulation{M, Δ, C, I, T, W, R, D, O, P, F, Π}
                 model :: M
                    Δt :: Δ
         stop_criteria :: C
        stop_iteration :: I
             stop_time :: T
       wall_time_limit :: W
              run_time :: R
           diagnostics :: D
        output_writers :: O
              progress :: P
    iteration_interval :: F
            parameters :: Π
end

"""
    Simulation(model; Δt,
         stop_criteria = Function[iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded],
        stop_iteration = Inf,
             stop_time = Inf,
       wall_time_limit = Inf,
           diagnostics = OrderedDict{Symbol, AbstractDiagnostic}(),
        output_writers = OrderedDict{Symbol, AbstractOutputWriter}(),
              progress = nothing,
    iteration_interval = 1,
            parameters = nothing)

Construct an Oceananigans.jl `Simulation` for a `model` with time step `Δt`.

Keyword arguments
=================
- `Δt`: Required keyword argument specifying the simulation time step. Can be a `Number`
  for constant time steps or a `TimeStepWizard` for adaptive time-stepping.
- `stop_criteria`: A list of functions or callable objects (each taking a single argument,
  the `simulation`). If any of the functions return `true` when the stop criteria is
  evaluated the simulation will stop.
- `stop_iteration`: Stop the simulation after this many iterations.
- `stop_time`: Stop the simulation once this much model clock time has passed.
- `wall_time_limit`: Stop the simulation if it's been running for longer than this many
   seconds of wall clock time.
- `progress`: A function with a single argument, the `simulation`. Will be called every
  `iteration_interval` iterations. Useful for logging simulation health.
- `iteration_interval`: How often to update the time step, check stop criteria, and call
  `progress` function (in number of iterations).
- `parameters`: Parameters that can be accessed in the `progress` function.
"""
function Simulation(model; Δt,
        stop_criteria = Any[iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded],
       stop_iteration = Inf,
            stop_time = Inf,
      wall_time_limit = Inf,
          diagnostics = OrderedDict{Symbol, AbstractDiagnostic}(),
       output_writers = OrderedDict{Symbol, AbstractOutputWriter}(),
             progress = nothing,
   iteration_interval = 1,
           parameters = nothing)

   if stop_iteration == Inf && stop_time == Inf && wall_time_limit == Inf
       @warn "This simulation will run forever as stop iteration = stop time " *
             "= wall time limit = Inf."
   end

   if Δt isa TimeStepWizard && iteration_interval == 1
       @warn "You have used the default iteration_interval=1. This simulation will " *
             "recalculate the time step every iteration which can be slow."
   end

   run_time = 0.0

   return Simulation(model, Δt, stop_criteria, stop_iteration, stop_time, wall_time_limit,
                     run_time, diagnostics, output_writers, progress, iteration_interval,
                     parameters)
end

Base.show(io::IO, s::Simulation) =
    print(io, "Simulation{$(typeof(s.model).name){$(typeof(s.model.architecture)), $(eltype(s.model.grid))}}\n",
            "├── Model clock: time = $(prettytime(s.model.clock.time)), iteration = $(s.model.clock.iteration) \n",
            "├── Next time step ($(typeof(s.Δt))): $(prettytime(get_Δt(s.Δt))) \n",
            "├── Iteration interval: $(s.iteration_interval)\n",
            "├── Stop criteria: $(s.stop_criteria)\n",
            "├── Run time: $(prettytime(s.run_time)), wall time limit: $(s.wall_time_limit)\n",
            "├── Stop time: $(prettytime(s.stop_time)), stop iteration: $(s.stop_iteration)\n",
            "├── Diagnostics: $(ordered_dict_show(s.diagnostics, "│"))\n",
            "└── Output writers: $(ordered_dict_show(s.output_writers, "│"))")


text:
export AAFT
"""
    AAFT()

An amplitude-adjusted-fourier-transform surrogate[^Theiler1991].

AAFT have the same linear correlation, or periodogram, and also
preserves the amplitude distribution of the original data.

AAFT can be used to test the null hypothesis that the data come from a
monotonic nonlinear transformation of a linear Gaussian process
(also called integrated white noise)[^Theiler1991].

[^Theiler1991]: J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
"""
struct AAFT <: Surrogate end

function surrogenerator(x, method::AAFT)
    init = surrogenerator(x, RandomFourier(true))
    return SurrogateGenerator(method, x, init)
end

function (rf::SurrogateGenerator{<:AAFT})()
    x = rf.x
    xs = sort(x)
    s = rf.init()
    s[sortperm(s)] .= xs
    return s
end


text:
module MPN

import ..XInts: lenvec, vec
using ..XInts: XInt, SVec

using Base.GC: @preserve
using Base.GMP: Limb # mp_limb_t
const PLimb = Ptr{Limb}
const Size = Clong # mp_size_t

len(p::Pair) = p[2]
ptr(p::Pair) = ptr(p[1])
lenvec(p::Pair) = p[2], vec(p[1])

vec(x::Vector) = x
vec(x::SVec) = x
ptr(p::Ptr) = p
len(p::AbstractArray) = length(p)
ptr(p::AbstractArray) = pointer(p) # Vector and views thereof

len(b::BigInt) = abs(b.size)
ptr(b::BigInt) = b.d
lenvec(b::BigInt) = len(b), b

isheap(x) = false
isheap(x::XInt) = true
isheap(x::BigInt) = true
isheap(x::Vector) = true
@inline isheap(x, y) = isheap(x) & isheap(y)

add_1(rp::PLimb, s1::PLimb, n, s2::Limb) = ccall((:__gmpn_add_1, :libgmp), Limb,
                                                 (PLimb, PLimb, Size, Limb),
                                                 rp, s1, n, s2)

add(rp::PLimb, s1::PLimb, l1, s2::PLimb, l2) = ccall((:__gmpn_add, :libgmp), Limb,
                                                     (PLimb, PLimb, Size, PLimb, Size),
                                                     rp, s1, l1, s2, l2)

sub_1(rp::PLimb, s1::PLimb, n, s2::Limb) = ccall((:__gmpn_sub_1, :libgmp), Limb,
                                                 (PLimb, PLimb, Size, Limb),
                                                 rp, s1, n, s2)

sub(rp::PLimb, s1::PLimb, l1, s2::PLimb, l2) = ccall((:__gmpn_sub, :libgmp), Limb,
                                                     (PLimb, PLimb, Size, PLimb, Size),
                                                     rp, s1, l1, s2, l2)

sub_n(rp::PLimb, s1::PLimb, s2::PLimb, n) = ccall((:__gmpn_sub_n, :libgmp), Limb,
                                                  (PLimb, PLimb, PLimb, Size),
                                                  rp, s1, s2, n)

cmp(s1, s2, n=len(s1)) = @preserve s1 s2 ccall((:__gmpn_cmp, :libgmp), Cint,
                                               (Ptr{Limb}, Ptr{Limb}, Clong),
                                               ptr(s1), ptr(s2), n)

lshift(rp, s1, c) = @preserve rp s1 ccall((:__gmpn_lshift, :libgmp), Limb,
                                          (PLimb, PLimb, Size, Cuint),
                                          ptr(rp), ptr(s1), len(s1), c)

rshift(rp, s1, c) = @preserve rp s1 ccall((:__gmpn_rshift, :libgmp), Limb,
                                          (PLimb, PLimb, Size, Cuint),
                                          ptr(rp), ptr(s1), len(s1), c)

and_n(rp, s1, s2, n=len(s1)) = @preserve rp s1 s2 ccall((:__gmpn_and_n, :libgmp), Cvoid,
                                                        (PLimb, PLimb, PLimb, Size),
                                                        ptr(rp), ptr(s1), ptr(s2), n)

mul_1(rp, s1, s2limb) = @preserve rp s1 ccall((:__gmpn_mul_1, :libgmp), Limb,
                                              (PLimb, PLimb, Size, Limb),
                                              ptr(rp), ptr(s1), len(s1), s2limb)

mul(rp, s1, s2) = @preserve rp s1 s2 ccall((:__gmpn_mul, :libgmp), Limb,
                                           (PLimb, PLimb, Size, PLimb, Size),
                                           ptr(rp), ptr(s1), len(s1), ptr(s2), len(s2))

get_str(str, base, s1, n=len(s1)) = @preserve str s1 ccall((:__gmpn_get_str, :libgmp), Size,
                                                           (Ptr{Cuchar}, Cint, PLimb, Size),
                                                           ptr(str), base, ptr(s1), n)

get_d(s1, n, sign, exp=0) = @preserve s1 ccall((:__gmpn_get_d, :libgmp), Cdouble,
                                               (PLimb, Size, Size, Clong),
                                               ptr(s1), n, sign, exp)

@inline function add_1(rv::Vector, x, y::Limb)
    xl, xv = lenvec(x)
    if isheap(xv)
        @preserve rv x add_1(ptr(rv), ptr(xv), xl, y) % Bool
    else
        add_1(rv, xv, xl, y)
    end
end

@inline function add_1(rv::Vector, xv, xl, y)
    @inbounds rv[1], c =  Base.add_with_overflow(xv[1], y)
    for i=2:xl
        ri = @inbounds xv[i]
        @inbounds rv[i] = ri + c
        c &= ri == typemax(Limb)
    end
    c
end

@inline function add(rp::Vector, x, y)
    xl, xv = lenvec(x)
    yl, yv = lenvec(y)
    # native MPN is not really faster than non-inlined Julia version here,
    # but if we inline Julia version (like we do now), there is no reason
    # to use MPN version
    use_MPN = false
    if use_MPN && isheap(xv, yv)
        @preserve rp xv yv add(ptr(rp), ptr(xv), xl, ptr(yv), yl) % Bool
    else
        add(rp, xv, xl, yv, yl)
    end
end

@inline function add(rv, xv, xl, yv, yl)
    # @assert xl >= yl
    c = false
    for i=1:yl
        xi = @inbounds xv[i]
        yi = @inbounds yv[i]
        ri, c2 = Base.add_with_overflow(xi, yi)
        @inbounds rv[i] = ri+c
        c = c2 | (c & (ri == typemax(Limb)))
    end
    for i=yl+1:xl
        ri = @inbounds xv[i]
        @inbounds rv[i] = ri + c
        c &= ri == typemax(Limb)
    end
    c
end

@inline function sub_1(rv::Vector, x, y::Limb)
    xl, xv = lenvec(x)
    if isheap(xv)
        @preserve rv x sub_1(ptr(rv), ptr(xv), xl, y) % Bool
    else
        sub_1(rv, xv, xl, y)
    end
end

@inline function sub_1(rv::Vector, xv, xl, y)
    @inbounds rv[1], c =  Base.sub_with_overflow(xv[1], y)
    for i=2:xl
        ri = @inbounds xv[i]
        @inbounds rv[i] = ri - c
        c &= iszero(ri)
    end
    c
end

@inline function sub(rp::Vector, x, y)
    xl, xv = lenvec(x)
    yl, yv = lenvec(y)
    if isheap(xv, yv)
        @preserve rp xv yv sub(ptr(rp), ptr(xv), xl, ptr(yv), yl) % Bool
    else
        sub(rp, xv, xl, yv, yl)
    end
end

@inline function sub(rv, xv, xl, yv, yl)
    # @assert xl >= yl
    c = sub_n(rv, xv, yv, 1:yl, false)
    for i=yl+1:xl
        ri = @inbounds xv[i]
        @inbounds rv[i] = ri-c
        c &= iszero(ri)
    end
    c
end

@inline function sub_n(rp::Vector, x, y, n)
    xv = vec(x)
    yv = vec(y)
    if isheap(xv, yv)
        @preserve rp xv yv sub_n(ptr(rp), ptr(xv), ptr(yv), n) % Bool
    else
        sub_n(rp, xv, yv, 1:n, false)
    end
end

@inline function sub_n(rv, xv, yv, r::AbstractUnitRange, c::Bool)
    for i=r
        xi = @inbounds xv[i]
        yi = @inbounds yv[i]
        ri, c2 = Base.sub_with_overflow(xi, yi)
        @inbounds rv[i] = ri-c
        c = c2 | (c & iszero(ri))
    end
    c
end

end # MPN


text:
module FrequencyDirectionSpectra

    using LinearAlgebra, RecipesBase, Distributions

    include("typestructure.jl")
    include("common.jl")
    include("swell.jl")
    include("JONSWAP.jl")
    include("MLM.jl")
    include("MEM.jl")
    include("plotting.jl")

end


text:
using EnergyExpressions
import EnergyExpressions: OrbitalMatrixElement

include("multipole_refs.jl")

@testset "Multipole expansions" begin
    @testset "Pure state" begin
        for (a,b,c,d) in [(SpinOrbital(o"1s", 0, half(1)),
                           SpinOrbital(o"1s", 0, -half(1)),
                           SpinOrbital(o"1s", 0, half(1)),
                           SpinOrbital(o"1s", 0, -half(1))),
                          (SpinOrbital(ro"1s", half(1)),
                           SpinOrbital(ro"1s", -half(1)),
                           SpinOrbital(ro"1s", half(1)),
                           SpinOrbital(ro"1s", -half(1))),
                          (SpinOrbital(o"1s", 0, half(1)),
                           SpinOrbital(o"2p", 1, half(1)),
                           SpinOrbital(o"1s", 0, half(1)),
                           SpinOrbital(o"2p", 1, half(1))),
                          (SpinOrbital(ro"1s", half(1)),
                           SpinOrbital(ro"2p", half(3)),
                           SpinOrbital(ro"1s", half(1)),
                           SpinOrbital(ro"2p", half(3)))]
            ome = OrbitalMatrixElement([a,b], CoulombInteraction(), [c,d])
            mome = multipole_expand(ome)
            @test length(mome.terms) == 1
            t = first(mome.terms)
            @test length(t.factors) == 1
            @test first(t.factors).o.k == 0
            @test t.coeff ≈ 1
        end
    end

    @testset "Direct interaction" begin
        orbitals = sos"k[s-d]"
        g = CoulombInteraction()
        for a in orbitals
            for b in orbitals
                ome = OrbitalMatrixElement([a,b], g, [a,b])
                mome = multipole_expand(ome)
                mref = ref_direct_multipoles(a, b)
                @test mome ≈ mref atol=1e-14
            end
        end
    end

    @testset "Exchange interaction" begin
        orbitals = sos"k[s-d]"
        g = CoulombInteraction()
        for a in orbitals
            for b in orbitals
                ome = OrbitalMatrixElement([a,b], g, [b,a])
                mome = multipole_expand(ome)
                if b.m[2] == a.m[2]
                    mref = ref_exchange_multipoles(a, b)
                    @test mome ≈ mref atol=1e-14
                else
                    @test iszero(mome)
                end
            end
        end
    end
end


text:
# script for Fig.5

using MultiscaleGraphSignalTransforms, LightGraphs, Plots; gr(dpi = 200)
import WaveletsExt: wiggle

## Build Graph
N = 512; G = path_graph(N)
X = zeros(N,2); X[:, 1] = 1:N
L = Matrix(laplacian_matrix(G))
𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1,:])'
W = 1.0 * adjacency_matrix(G)

## Build NGWPs
Gstar_Sig = GraphSig(W)
G_Sig = GraphSig(W, xy = X)
GP_dual = partition_tree_fiedler(Gstar_Sig; swapRegion = false)
GP_primal = pairclustering(𝚽, GP_dual)

@time VM_NGWP = vm_ngwp(𝚽, GP_dual)

#################### Fig.5
j = 5
for k in [1, 2, 5]
   WW = sort_wavelets(VM_NGWP[GP_dual.rs[k, j]:(GP_dual.rs[k + 1, j] - 1), j, :]')
   if k == 2
      WW[:, end] *= -1
   end
   plt = wiggle(WW; sc = 0.75)
   # savefig(plt, joinpath(@__DIR__, "../paperfigs/Path512_VM_NGWP_j$(j-1)k$(k-1).png"))
end
current()


text:
# default settings
const default_chunk_size = 0
const default_verbose_level = 1

# interesting... using semicolons would make it into a 1-dimensional array
const magic = [
         b"\x00\x00\x00\x00\x00\x00\x00\x00"  ;
         b"\x00\x00\x00\x00\xc2\xea\x81\x60"  ;
         b"\xb3\x14\x11\xcf\xbd\x92\x08\x00"  ;
         b"\x09\xc7\x31\x8c\x18\x1f\x10\x11" ]

const align_1_checker_value = b"3"
const align_1_offset = 32
const align_1_length = 1
const align_1_value = 4
const u64_byte_checker_value = b"3"
const align_2_offset = 35
const align_2_length = 1
const align_2_value = 4
const endianness_offset = 37
const endianness_length = 1
const platform_offset = 39
const platform_length = 1
const encoding_offset = 70
const encoding_length = 1
const dataset_offset = 92
const dataset_length = 64
const file_type_offset = 156
const file_type_length = 8
const date_created_offset = 164
const date_created_length = 8
const date_modified_offset = 172
const date_modified_length = 8
const header_size_offset = 196
const header_size_length = 4
const page_size_offset = 200
const page_size_length = 4
const page_count_offset = 204
const page_count_length = 4
const sas_release_offset = 216
const sas_release_length = 8
const sas_server_type_offset = 224
const sas_server_type_length = 16
const os_version_number_offset = 240
const os_version_number_length = 16
const os_maker_offset = 256
const os_maker_length = 16
const os_name_offset = 272
const os_name_length = 16
const page_bit_offset_x86 = 16
const page_bit_offset_x64 = 32
const subheader_pointer_length_x86 = 12
const subheader_pointer_length_x64 = 24
const page_type_offset = 0
const page_type_length = 2
const block_count_offset = 2
const block_count_length = 2
const subheader_count_offset = 4
const subheader_count_length = 2
const page_meta_type = 0
const page_data_type = 256
const page_amd_type = 1024
const page_metc_type = 16384
const page_comp_type = -28672
const page_mix_types = [512, 640]
const page_meta_data_mix_types = vcat(page_meta_type, page_data_type, page_mix_types)
const subheader_pointers_offset = 8
const subheader_comp_uncompressed = 0
const subheader_comp_truncated = 1
const subheader_comp_compressed = 4
const text_block_size_length = 2
const row_length_offset_multiplier = 5
const row_count_offset_multiplier = 6
const col_count_p1_multiplier = 9
const col_count_p2_multiplier = 10
const row_count_on_mix_page_offset_multiplier = 15
const column_name_pointer_length = 8
const column_name_text_subheader_offset = 0
const column_name_text_subheader_length = 2
const column_name_offset_offset = 2
const column_name_offset_length = 2
const column_name_length_offset = 4
const column_name_length_length = 2
const column_data_offset_offset = 8
const column_data_length_offset = 8
const column_data_length_length = 4
const column_type_offset = 14
const column_type_length = 1
const sas_date_origin = Date(1960, 1, 1)
const sas_datetime_origin = DateTime(sas_date_origin)

const rle_compression = b"SASYZCRL"
const rdc_compression = b"SASYZCR2"

const compression_method_none = 0x00
const compression_method_rle  = 0x01
const compression_method_rdc  = 0x02

# Courtesy of Evan Miller's ReadStat project
# Source: https://github.com/WizardMac/ReadStat/blob/master/src/sas/readstat_sas.c
const encoding_names = Dict(
    20 => "UTF-8",
    28 => "US-ASCII",
    29 => "ISO-8859-1",
    30 => "ISO-8859-2",
    31 => "ISO-8859-3",
    34 => "ISO-8859-6",
    35 => "ISO-8859-7",
    36 => "ISO-8859-8",
    39 => "ISO-8859-11",
    40 => "ISO-8859-9",
    60 => "WINDOWS-1250",
    61 => "WINDOWS-1251",
    62 => "WINDOWS-1252",
    63 => "WINDOWS-1253",
    64 => "WINDOWS-1254",
    65 => "WINDOWS-1255",
    66 => "WINDOWS-1256",
    67 => "WINDOWS-1257",
    118 => "CP950",
    119 => "EUC-TW",
    123 => "BIG5-HKSCS",
    125 => "GB18030",
    126 => "CP936",
    134 => "EUC-JP",
    138 => "CP932",
    140 => "EUC-KR",
    141 => "CP949"
)

const index_rowSizeIndex = 0
const index_rowSizeIndex = 0
const index_columnSizeIndex = 1
const index_subheaderCountsIndex = 2
const index_columnTextIndex = 3
const index_columnNameIndex = 4
const index_columnAttributesIndex = 5
const index_formatAndLabelIndex = 6
const index_columnListIndex = 7
const index_dataSubheaderIndex = 8
const index_end_of_header = 100

const subheader_signature_to_index = Dict(

    b"\xF7\xF7\xF7\xF7"                 => index_rowSizeIndex,
    b"\xF7\xF7\xF7\xF7\x00\x00\x00\x00" => index_rowSizeIndex,
    b"\x00\x00\x00\x00\xF7\xF7\xF7\xF7" => index_rowSizeIndex,
    b"\xF7\xF7\xF7\xF7\xFF\xFF\xFB\xFE" => index_rowSizeIndex,

    b"\xF6\xF6\xF6\xF6"                 => index_columnSizeIndex,
    b"\x00\x00\x00\x00\xF6\xF6\xF6\xF6" => index_columnSizeIndex,
    b"\xF6\xF6\xF6\xF6\x00\x00\x00\x00" => index_columnSizeIndex,
    b"\xF6\xF6\xF6\xF6\xFF\xFF\xFB\xFE" => index_columnSizeIndex,

    b"\x00\xFC\xFF\xFF"                 => index_subheaderCountsIndex,
    b"\xFF\xFF\xFC\x00"                 => index_subheaderCountsIndex,
    b"\x00\xFC\xFF\xFF\xFF\xFF\xFF\xFF" => index_subheaderCountsIndex,
    b"\xFF\xFF\xFF\xFF\xFF\xFF\xFC\x00" => index_subheaderCountsIndex,

    b"\xFD\xFF\xFF\xFF"                 => index_columnTextIndex,
    b"\xFF\xFF\xFF\xFD"                 => index_columnTextIndex,
    b"\xFD\xFF\xFF\xFF\xFF\xFF\xFF\xFF" => index_columnTextIndex,
    b"\xFF\xFF\xFF\xFF\xFF\xFF\xFF\xFD" => index_columnTextIndex,

    b"\xFF\xFF\xFF\xFF"                 => index_columnNameIndex,
    b"\xFF\xFF\xFF\xFF\xFF\xFF\xFF\xFF" => index_columnNameIndex,

    b"\xFC\xFF\xFF\xFF"                 => index_columnAttributesIndex,
    b"\xFF\xFF\xFF\xFC"                 => index_columnAttributesIndex,
    b"\xFC\xFF\xFF\xFF\xFF\xFF\xFF\xFF" => index_columnAttributesIndex,
    b"\xFF\xFF\xFF\xFF\xFF\xFF\xFF\xFC" => index_columnAttributesIndex,

    b"\xFE\xFB\xFF\xFF"                 => index_formatAndLabelIndex,
    b"\xFF\xFF\xFB\xFE"                 => index_formatAndLabelIndex,
    b"\xFE\xFB\xFF\xFF\xFF\xFF\xFF\xFF" => index_formatAndLabelIndex,
    b"\xFF\xFF\xFF\xFF\xFF\xFF\xFB\xFE" => index_formatAndLabelIndex,

    b"\xFE\xFF\xFF\xFF"                 => index_columnListIndex,
    b"\xFF\xFF\xFF\xFE"                 => index_columnListIndex,
    b"\xFE\xFF\xFF\xFF\xFF\xFF\xFF\xFF" => index_columnListIndex,
    b"\xFF\xFF\xFF\xFF\xFF\xFF\xFF\xFE" => index_columnListIndex
)

# List of frequently used SAS date and datetime formats
# http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_intervals_sect009.htm
# https://github.com/epam/parso/blob/master/src/main/java/com/epam/parso/impl/SasFileConstants.java
const sas_date_formats = ["DATE", "DAY", "DDMMYY", "DOWNAME", "JULDAY", "JULIAN",
                    "MMDDYY", "MMYY", "MMYYC", "MMYYD", "MMYYP", "MMYYS",
                    "MMYYN", "MONNAME", "MONTH", "MONYY", "QTR", "QTRR",
                    "NENGO", "WEEKDATE", "WEEKDATX", "WEEKDAY", "WEEKV",
                    "WORDDATE", "WORDDATX", "YEAR", "YYMM", "YYMMC", "YYMMD",
                    "YYMMP", "YYMMS", "YYMMN", "YYMON", "YYMMDD", "YYQ",
                    "YYQC", "YYQD", "YYQP", "YYQS", "YYQN", "YYQR", "YYQRC",
                    "YYQRD", "YYQRP", "YYQRS", "YYQRN",
                    "YYMMDDP", "YYMMDDC", "E8601DA", "YYMMDDN", "MMDDYYC",
                    "MMDDYYS", "MMDDYYD", "YYMMDDS", "B8601DA", "DDMMYYN",
                    "YYMMDDD", "DDMMYYB", "DDMMYYP", "MMDDYYP", "YYMMDDB",
                    "MMDDYYN", "DDMMYYC", "DDMMYYD", "DDMMYYS",
                    "MINGUO"]

const sas_datetime_formats = ["DATETIME", "DTWKDATX",
                        "B8601DN", "B8601DT", "B8601DX", "B8601DZ", "B8601LX",
                        "E8601DN", "E8601DT", "E8601DX", "E8601DZ", "E8601LX",
                        "DATEAMPM", "DTDATE", "DTMONYY", "DTMONYY", "DTWKDATX",
                        "DTYEAR", "TOD", "MDYAMPM"]

const zero_space = b"\x00 "

const column_type_none    = 0x00
const column_type_decimal = 0x01
const column_type_string  = 0x02

const VENDOR_STAT_TRANSFER = 0x00
const VENDOR_SAS           = 0x01

const FALLBACK_ENCODING = "UTF-8"
const ENCODINGS_OK_WITH_BASE_TRANSCODER = [ "UTF-8" , "US-ASCII" ]

const REGULAR_STR_ARRAY(n) = Array{String}(undef, n)
const EMPTY_STRING = ""


text:
# This file implements the DotProductNode and corresponding sum-product rules

# using ForneyLab

# import ForneyLab: DeltaFactor, Univariate, Multivariate, generateId, associate!, addNode!
export DotProduct

"""
Description:

    out = β⋅in

    in1: d-dimensional vector
    in2: d-dimensional vector
    out: scalar

           in2
           |
      in1  V   out
    ----->[⋅]----->

    f(in1, in2, out) =  δ(out - in1'*in2)

Interfaces:

    1 i[:in1], 2 i[:in2], 3 i[:out]

Construction:

    DotProduct(out, in1, in2, id=:my_node)
"""
mutable struct DotProduct <: DeltaFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function DotProduct(out::Variable, in1::Variable, in2::Variable; id=generateId(DotProduct))
        self = new(id, Array{Interface}(3), Dict{Int,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:in1] = self.interfaces[2] = associate!(Interface(self), in1)
        self.i[:in2] = self.interfaces[3] = associate!(Interface(self), in2)

        return self
    end
end

slug(::Type{DotProduct}) = "dot"

# Update rule definitions

#             β
#             | ↓ N{d}
#             |
#             V
#    x ----->[⋅]-----> y
#        →        →
#       δ{d}      N

@sumProductRule(:node_type     => DotProduct,
                :outbound_type => Message{Gaussian},
                :inbound_types => (Nothing, Message{PointMass}, Message{Gaussian}),
                :name          => SPDotProductOutVPG)

function ruleSPDotProductOutVPG(msg_out::Nothing,
                            msg_in1::Message{PointMass, Multivariate},
                            msg_in2::Message{Gaussian, Multivariate})
    x = msg_in1.dist.params[:m]
    β = convert(ProbabilityDistribution{Multivariate,GaussianMeanVariance}, msg_in2.dist)
#     β = ensureParameters!(msg_in2.dist, (:m, :V))
    out_m = x' * β.m
    out_v = x' * β.V * x

    return Message(Multivariate, Gaussian, m=out_m, v=out_v)

end

#             β
#             | ↓ δ{d}
#             |
#             V
#    x ----->[⋅]-----> y
#        →        →
#       N{d}      N

@sumProductRule(:node_type     => DotProduct,
                :outbound_type => Message{Gaussian},
                :inbound_types => (Nothing, Message{Gaussian}, Message{PointMass}),
                :name          => SPDotProductOutVGP)

function ruleSPDotProductOutVGP(msg_out::Nothing,
                            msg_in1::Message{Gaussian, Multivariate},
                            msg_in2::Message{PointMass, Multivariate})
    SPDotProductOutVGP(msg_out, msg_in2, msg_in1)
end

#             β
#             | ↑ N{d}
#             |
#             V
#    x ----->[⋅]-----> y
#        →        ←
#       δ{d}      N

@sumProductRule(:node_type     => DotProduct,
                :outbound_type => Message{Gaussian},
                :inbound_types => (Message{Gaussian}, Message{PointMass}, Nothing),
                :name          => SPDotProductIn2GPV)

function ruleSPDotProductIn2GPV(msg_out::Message{Gaussian, Univariate},
                            msg_in1::Message{PointMass, Multivariate},
                            msg_in2::Nothing)
    x = msg_in1.dist.params[:m] # We'll call in1 x
    d = length(x)
    y = convert(ProbabilityDistribution{Univariate,GaussianWeightedMeanPrecision}, msg_out.dist)
#     y = ensureParameters!(msg_out.dist, (:xi, :w))
    xi = x * y.params[:xi]
    w = x * y.params[:w] * x'

    return Message(Multivariate, Gaussian, xi=xi, w=w)
end

#             β
#             | ↓ δ{d}
#             |
#             V
#    x ----->[⋅]-----> y
#        ←        ←
#       N{d}      N

@sumProductRule(:node_type     => DotProduct,
                :outbound_type => Message{Gaussian},
                :inbound_types => (Message{Gaussian}, Nothing, Message{PointMass}),
                :name          => SPDotProductIn1GVP)

function ruleSPDotProductIn1GVP(msg_out::Message{Gaussian, Univariate},
                            msg_in1::Nothing,
                            msg_in2::Message{PointMass, Multivariate})
    SPDotProductIn2GPV(msg_out, msg_in2, msg_in1)
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
# Autogenerated wrapper script for Ipopt_jll for armv7l-linux-musleabihf-libgfortran4-cxx03
export amplexe, libipopt, libipoptamplinterface

using ASL_jll
using MUMPS_seq_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Ipopt")
JLLWrappers.@declare_executable_product(amplexe)
JLLWrappers.@declare_library_product(libipopt, "libipopt.so.3")
JLLWrappers.@declare_library_product(libipoptamplinterface, "libipoptamplinterface.so.3")
function __init__()
    JLLWrappers.@generate_init_header(ASL_jll, MUMPS_seq_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        amplexe,
        "bin/ipopt",
    )

    JLLWrappers.@init_library_product(
        libipopt,
        "lib/libipopt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libipoptamplinterface,
        "lib/libipoptamplinterface.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ThermoProbe
include("IdealGas.jl")
include("Utils.jl")
using LightXML
using Printf

"""
This function read the input file specidfied
    and calculates the pure species and mixture properties at
    the specified temperature and pressure
#   Usage:
-   run(input_file, lib_dir)
-   input_file: name of the input file (String)
-   lib_dir: path to the folder where the data files are present
"""
function run(input_file::AbstractString, lib_dir::AbstractString)

    #load the xml file
    xmldoc = parse_file(input_file)
    xmlroot = root(xmldoc)

    local molwt = Array{Float64,1}()
    local mole_fracs = Array{Float64,1}()
    local gas_species = Array{String,1}()
    local thermo_obj::IdealGas.SpeciesThermoObj

    gas_species = get_collection_from_xml(xmlroot,"gasphase")
    thermo_file = lib_dir*"therm.dat"
    thermo_obj = IdealGas.create_thermo(gas_species, thermo_file)

    #Get the species composition
    mole_fracs = get_molefraction_from_xml(xmlroot,thermo_obj.molwt,gas_species)
    #Get the temperature
    T = get_value_from_xml(xmlroot,"T")
    #Get the pressire
    p = get_value_from_xml(xmlroot,"p")

    #calculation of pure species thermodynamic properties
    enthalpy = Array{Float64,1}()
    entropy = Array{Float64,1}()
    cp = Array{Float64,1}()

    enthalpy = IdealGas.H_all(thermo_obj,T)
    entropy = IdealGas.S_all(thermo_obj,T)
    cp = IdealGas.cp_all(thermo_obj,T)

    println("\nThermodynamic properties of pure species at $T(K) and $p(Pa)\n")
    println("Species \t Enthalpy(J/mol) \t Entropy(J/mol-K) \t cp(J/mol-K)\n")
    for i in eachindex(gas_species)
        @printf("%s \t\t %+.4e \t\t %+.4e \t\t %.4e \n", gas_species[i], enthalpy[i],  entropy[i],  cp[i])
    end

    #Calculation of mixture properties
    hmix = IdealGas.Hmix(thermo_obj,T,mole_fracs)
    cpmix = IdealGas.cpmix(thermo_obj,T,mole_fracs)
    gmix = IdealGas.Gmix(thermo_obj,T,p,mole_fracs)
    println("\nMixture properties\n")
    @printf("Enthalpy (J/mol): %e \n", hmix)
    @printf("Specific Heat (J/mol-K): %e \n", cpmix)
    @printf("Gibbs free energy (J/mol): %e \n ", gmix)

    return Symbol("Success")
end

end


text:
module DiffEqNoiseProcess

using ResettableStacks, DiffEqBase, RecipesBase
using RecursiveArrayTools, StaticArrays, Random, Statistics
using LinearAlgebra

import RandomNumbers: Xorshifts

import RandomNumbers, Random123

import DiffEqBase: isinplace, solve, AbstractNoiseProcess,
       DEIntegrator, AbstractNoiseProblem

import PoissonRandom, Distributions

import QuadGK, Optim

import GPUArrays

using DiffEqBase: @..

include("types.jl")
include("wiener.jl")
include("solve.jl")
include("geometric_bm.jl")
include("compoundpoisson.jl")
include("ornstein_uhlenbeck.jl")
include("rswm.jl")
include("bridges.jl")
include("noise_interfaces/simple_noise_process_interface.jl")
include("noise_interfaces/noise_process_interface.jl")
include("noise_interfaces/noise_function_interface.jl")
include("noise_interfaces/virtual_brownian_tree_interface.jl")
include("noise_interfaces/noise_grid_interface.jl")
include("noise_interfaces/noise_approximation_interface.jl")
include("noise_interfaces/noise_wrapper_interface.jl")
include("noise_interfaces/box_wedge_tail_interface.jl")
include("noise_interfaces/common.jl")
include("correlated_noisefunc.jl")
include("pCN.jl")

export RSWM

export NoiseProcess, SimpleNoiseProcess, adaptive_alg

export WienerProcess, WienerProcess!, SimpleWienerProcess, SimpleWienerProcess!

export RealWienerProcess, RealWienerProcess!

export BrownianBridge, BrownianBridge!

export GeometricBrownianMotionProcess, GeometricBrownianMotionProcess!

export CompoundPoissonProcess, CompoundPoissonProcess!

export GeometricBrownianBridge, GeometricBrownianBridge!

export CompoundPoissonBridge, CompoundPoissonBridge!

export OrnsteinUhlenbeckProcess, OrnsteinUhlenbeckProcess!

export NoiseWrapper, NoiseFunction, NoiseGrid, NoiseApproximation

export VirtualBrownianTree, VirtualBrownianTree!

export BoxWedgeTail, BoxWedgeTail!

export accept_step!, reject_step!, calculate_step!, setup_next_step!, save_noise!

export CorrelatedWienerProcess, CorrelatedWienerProcess!

export pCN, pCN!

end # module


text:
using TensorTrains

w=2
L=5
h,V = hubbard_2D(w,L;U=4)
H_tto = hV_to_mpo(h,V,tol=1e-5)
ψ_0 = half_filling(w*L)
E_als,ψ_tt = als_eigsolv(H_tto,ψ_0,sweep_schedule=[2,5],rmax_schedule=[16,32])
E_mals,ψ_tt,r = mals_eigsolv(H_tto,ψ_tt,tol=1e-6,sweep_schedule=[2,5],rmax_schedule=[32,48])

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
"""
GS step for CSC
"""

"""
gsstep!(x::Vector{Tv}, Q::SparseMatrix.SparseCSC{Tv,Ti}, b::Vector{Tv}; alpha::Tv=Tv(1), sigma::Tv=Tv(0), omega::Tv=Tv(1))::Nothing
gsstep!(x::Vector{Tv}, Q::SparseMatrix.SparseCSR{Tv,Ti}, b::Vector{Tv}; alpha::Tv=Tv(1), sigma::Tv=Tv(0), omega::Tv=Tv(1))::Nothing

GS (Gauss-Seidal) or SOR (Successive Over Relaxation) step for the following linear equation

    alpha * trans(A - sigma I) * x = b

    notrans:
        x := (D/omega + L)^(-1) (b/alpha - (U - D (1-omega)/omega - sigma I) * x)
    trans:
        x := (D/omega + tr(U))^(-1) (b/alpha - (tr(L) - D (1-omega)/omega - sigma I) * x)

    where
        A: square matrix
        x: vector (in; initial vector for the step, out; updated vector)
        b: constant vector

Note that notrans and trans are determined by the type of matrix Q.
If Q is SparseCSC, gstep! provides the step for trans.
If Q is SparseCSR, gsstep! provides the step for notrans.
"""

function gsstep!(x::Vector{Tv}, Q::SparseCSC{Tv,Ti}, b::Vector{Tv};
        alpha::Tv=Tv(1), sigma::Tv=Tv(0), omega::Tv=Tv(1))::Nothing where {Tv, Ti}
    m, n = size(Q)
    @assert m == n
    @inbounds for j = 1:n
        tmpd::Tv = 0
        tmpx::Tv = b[j] / alpha
        for z = Q.colptr[j]:Q.colptr[j+1]-1
            i = Q.rowind[z]
            if i == j
                tmpd = Q.val[z]
                tmpx += sigma * x[i]
            else
                tmpx -= Q.val[z] * x[i]
            end
        end
        x[j] = omega / tmpd * tmpx + (1 - omega) * x[j]
    end
    nothing
end

function gsstep!(x::Vector{Tv}, Q::SparseCSR{Tv,Ti}, b::Vector{Tv};
        alpha::Tv=Tv(1), sigma::Tv=Tv(0), omega::Tv=Tv(1))::Nothing where {Tv, Ti}
    m, n = size(Q)
    @assert m == n
    @inbounds for i = 1:m
        tmpd::Tv = 0
        tmpx::Tv = b[j] / alpha
        for z = Q.rowptr[i]:Q.rowptr[i+1]-1
            j = Q.colind[z]
            if i == j
                tmpd = Q.val[z]
                tmpx += sigma * x[j]
            else
                tmpx -= Q.val[z] * x[j]
            end
        end
        x[i] = omega / tmpd * tmpx + (1 - omega) * x[i]
    end
    nothing
end


text:
# an (stack-allocatable) array that errors if you try to index it
# This is useful for the idiom in C of passing a null pointer to a function
# to signal that an argument is was not supplied

import Base: size, getindex, setindex!

"""
  Special array type that errors if indexed.
  This is useful for the idiom in C of passing a null pointer to a function
  to signal that an argument is was not supplied

  This array can be stack-allocated.  It supportes the standard Array
  constructors.
"""
struct NullArray{T, N} <: DenseArray{T, N}
  dims::NTuple{N, Int}
  offset::Int  # compatability with ArrayViews

  function NullArray{T, N}(dims::NTuple{N, Int}, offset::Integer) where {T, N}
    return new(dims, offset)
  end

  function NullArray{T, N}(dims::NTuple{N, Int}) where {T, N}
    return new(dims, 0)
  end
end

size(A::NullArray) = A.dims

function getindex(A::NullArray{T, 1}, i::Int) where {T}
  throw(ErrorException("getindex! not defind for NullArray"))
end

function setindex!(A::NullArray, v, i::Int)
  throw(ErrorException("setindex! not defind for NullArray"))
end

# construct from individual dimensions
# These constructors can be called as NullArray{T}(i, j) (which is not
# documented in Julia).
function NullArray{T}(i::Integer) where {T}
  return NullArray{T, 1}( (Int(i), ))
end

function NullArray{T}(i::Integer, j::Integer) where {T}
  return NullArray{T, 2}( (Int(i), Int(j)))
end

function NullArray{T}(i::Integer, j::Integer, k::Integer) where {T}
  return NullArray{T, 3}( (Int(i), Int(j), Int(k)) )
end

function NullArray{T}(i::Integer, j::Integer, k::Integer, l::Integer
                         ) where {T}
  return NullArray{T, 4}( (Int(i), Int(j), Int(k), Int(l)) )
end

function NullArray{T}(i::Integer, j::Integer, k::Integer, l::Integer,
                          m::Integer) where {T}
  return NullArray{T, 5}( (Int(i), Int(j), Int(k), Int(l), Int(m)) )
end

function NullArray{T}(i::Integer, j::Integer, k::Integer, l::Integer,
                          m::Integer, n::Integer) where {T}
  return NullArray{T, 6}( (Int(i), Int(j), Int(k), Int(l), Int(m), Int(n)) )
end

function NullArray{T}(i::Integer, j::Integer, k::Integer, l::Integer,
                          m::Integer, n::Integer, o::Integer) where {T}
  return NullArray{T, 7}( (Int(i), Int(j), Int(k), Int(l), Int(m), Int(n), Int(o)) )
end

# NTuple constructor
function NullArray{T}(d::NTuple{T, Int})
  return NullArray{T, N}(d)
end

#------------------------------------------------------------------------------
# compatability with sview (limited to contiguous views, 2D parent array)

import ODLCommonTools.sview

function sview(A::NullArray{T, 2}, rng::UnitRange, i::Integer) where {T}
  dims = (Int(length(rng)), )
  offset = rng[1] - 1
  return NullArray{T, 1}(dims, offset)
end

function sview(A::NullArray{T, 2}, idx::Colon, i::Integer) where {T}

  dims = (Int(A.dims[1]),)
  return NullArray{T, 1}(dims)
end


text:
using Random: Sampler
using Distributions: Geometric

struct Typing_sampler{DL, DT, T} <: Sampler{Vector{Vector{T}}}
    hyperdegree::Int
    length::DL
    letter::DT
end
Typing_sampler(hyperdegree::Int, n::Int) =
    Typing_sampler{Geometric{Float64}, UnitRange{Int}, Int}(
        hyperdegree, Geometric(1/(n+1)), 1:n)
Typing_sampler(hyperdegree, space_weight, weights, letters) =
    Typing_sampler{Geometric{Float64}, AliasTable{Char, Vector{Char}}, Char}(
        hyperdegree, Geometric(space_weight), AliasTable(weights, letters))

Base.rand(rng::AbstractRNG, s::Typing_sampler) =
    [rand(rng, s.letter, rand(rng, s.length)) for _ in 1:s.hyperdegree]

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
module ModSimOpt
        using Distributions
        using Plots

        export
        ConditionalDistribution,
        GaussianConditional,
        DiscreteConditional,

        ### methods
        inconsistent_gaussian_gibbs

        ### source files
        include("incompatible_gibbs.jl")
end # module


text:
using ApproxFun, SpectralMeasures, Plots
gr(linewidth=3,legend=true,xlims=(0,100))

D = freejacobioperator() - I

## Discrete Diffusion equations
u0 = [zeros(20);ones(5);zeros(20);1.2*ones(3)]

# J1 = D
J1 = D
x1,U1 = eigen(J1)
plot(pad(u0,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Pure diffusion on the half-line")
  for t in [1,10,100,1000]
    ut1 = U1\(exp(x1*t)*(U1*u0))
    plot!(pad(ut1.coefficients,100),ylims=(0,1.5),marker=2,label="t=$(t)")
  end
  plot!()
Plots.savefig("purediffusion1.pdf")

J1n = full(J1[1:70,1:70])
plot(pad(u0,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Pure diffusion on the half-line (n=70 finite section method)")
  for t in [1,10,100,1000]
    ut1 = expm(J1n*t)*pad(u0,70)
    push!(areasfs,sum(ut1))
    push!(energiesfs,sum(ut1.^2))
    plot!(pad(ut1,100),ylims=(0,1.5),marker=2,label="t=$(t)")
  end
  plot!()
Plots.savefig("purediffusion1finitesection.pdf")

# J2 = D + diag(v) where v is a reaction term
react = [zeros(35);-1.;-1.;-1]
J2 = D + SymTriOperator(react,[0.])
x2,U2 = eigen(J2)
plot(pad(u0,100),marker=2,label="t=0",ylims=(0,1.5),xlabel="k",ylabel="u_k",title="Diffusion plus reaction on the half-line")
  for t in [1,10,100,1000]
    ut2 = U2\(exp(x2*t)*(U2*u0))
    plot!(pad(ut2.coefficients,100),ylims=(0,1.5),marker=2,label="t=$(t)")
  end
  plot!()
plot!(pad(-react,100),label="Reaction potential",color=:black,ylims=(0,1.5),linewidth=2)
Plots.savefig("reactiondiffusion.pdf")
eigs = discreteeigs(J2)
plot(U2[1,1:100],marker=2,title="Stable eigenmodes of reaction-diffusion Jacobi operator",xlabel="k",ylabel="u_k",label="v_1")
  plot!(U2[2,1:100],marker=2,label="v_2")
  plot!(pad(react,100),color=:black,label="Reaction potential",yticks=-1.5:0.5:1.5)
Plots.savefig("stableeigenmodes.pdf")

# Pure fractional diffusion α = 0.65
J4 = D
x4,U4 = eigen(J4)
α = 0.65
plot(pad(u0,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Pure fractional diffusion on the half-line, \\alpha = 0.65")
  for t in [1,10,100,1000]
    h4 = Fun(x->exp(-(abs(x)^α)*t),space(x4))
    ut4 = U4\(h4*(U4*u0))
    plot!(pad(ut4.coefficients,100),ylims=(0,1.5),marker=2,label="t=$(t)")
  end
  plot!()
Plots.savefig("purefracdiffusion1.pdf")

# Pure fractional diffusion α = 0.85
J5 = D
x5,U5 = eigen(J5)
α = 0.85
plot(pad(u0,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Pure fractional diffusion on the half-line, \\alpha = 0.85")
  for t in [1,10,100,1000]
    h5 = Fun(x->exp(-(abs(x)^α)*t),space(x5))
    ut5 = U5\(h5*(U5*u0))
    plot!(pad(ut5.coefficients,100),ylims=(0,1.5),marker=2,label="t=$(t)")
  end
  plot!()
Plots.savefig("purefracdiffusion085.pdf")

# Pure fractional diffusion α = .25
J6 = D
x6,U6 = eigen(J6)
α = .25
plot(pad(u0,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Pure fractional diffusion on the half-line, \\alpha = 0.25")
  for t in [1,10,100,1000]
    h6 = Fun(x->exp(-(abs(x)^α)*t),space(x6))
    ut6 = U6\(h6*(U6*u0))
    plot!(pad(ut6.coefficients,100),ylims=(0,1.5),marker=2,label="t=$(t)")
  end
  plot!()
Plots.savefig("purefracdiffusion025.pdf")

# Fractional diffusion with reaction α = 0.85
react = [zeros(20);1.;1.;1.]
J7 = D + SymTriOperator(react,[0.])
discreteeigs(J7)
x7,U7 = eigen(J7)
α = 0.85
plot(pad(u0,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Fractional diffusion with reaction on the half-line, \\alpha = 0.85")
#for t in [1;10;100;1000]
t = 1
h7 = Fun(x->exp(-(abs(x)^α)*t),space(x7)[3])
pt1 = space(x7)[1].points[1]
pt2 = space(x7)[2].points[1]
h7 = Fun(space(x7),[exp((pt1^α)*t);exp((pt2^α)*t);coefficients(h7)])

h7 = Fun(x->exp(-(abs(x)^α)*t),space(x7))

ut7 = U7\(h7*(U7*u0))
plot!(pad(ut7.coefficients,100),ylims=(0,1.5),marker=2,label="t=$(t)")
#end
plot!()
plot!(pad(react,100),color=:black,label="Reaction potential")
Plots.savefig("fracreactdiffusion085withpot.pdf")

SymTriPertToeplitz(V8+1.0,[.5],1.0,-.5)

# Discrete Schrodinger equations
V8 = [zeros(40);2*ones(5)]
  J8 = -D + SymTriOperator(V8,[0.])
  x8,U8 = eigen(J8)
  λ8 = discreteeigs(J8)
  n8 = length(λ8)
  Eigv8 = full(U8[1:n8,1:100])

Eigv8[1,40:100]

n8

V9 = [zeros(40);1*ones(3)]
  J9 = D + SymTriOperator(V9,[0.])
  x9,U9 = eigen(J9)
  λ9 = discreteeigs(J9)
  n9 = length(λ9)
  Eigv9 = full(U9[1:n9,1:100])

if n8 > 0
  plot(Eigv8[1,:]/4+λ8[1]-2,marker=2,title="Tunneling barrier potential and eigenstates",xlabel="k",ylabel="u_k",label="v_1")
  if n8 > 1
    for k in 2:n8
      plot!(Eigv8[k,:]/4+λ8[k]-2,marker=2,label="v_$(k)")
    end
  end
  plot!(pad(V8,100),color=:black,label="Potential")
end
Plots.savefig("tunnelbarriereig.pdf")

if n9 > 0
  plot(Eigv9[1,:]/4+λ9[1],marker=2,title="Double well potential and eigenstates (far)",xlabel="k",ylabel="u_k",label="v_1",ylims=(-.075,.4))
  if n9 > 1
    for k in 2:n9
      plot!((-1)^(k+1)*Eigv9[k,:]/4+λ9[k],marker=2,label="v_$(k)")
    end
  end
  plot!(.2-pad(V9,100),color=:black,label="Potential")
end
Plots.savefig("doublewelleigfar.pdf")

V8 = [zeros(40);2*ones(5)]
  J8 = -D + SymTriOperator(V8,[0.])
  n = 400
  uinit = pad!(exp(-(-28:.75:28).^2),n)
  Jn = full(J8[1:n,1:n])
t = 20
  ut8 = expm(im*Jn*t)*uinit
  plot(abs(ut8[1:100]),ylims=(-.5,1.5),marker=2,label="t=$(t) (abs)")
  plot!(pad(V8,100),ylims=(-.5,1.5),color=:black,label="Potential",xlabel="k",ylabel="u_k",title="No quantum tunneling on the half-line")

d,Q = eigen(Jn)

for t in [0;5;10;20;30;40]
  ut8 = expm(im*Jn*t)*uinit
  plot(abs(ut8[1:100]),ylims=(-.5,1.5),marker=2,label="t=$(t) (abs)")
  plot!(pad(V8,100),ylims=(-.5,1.5),color=:black,label="Potential",xlabel="k",ylabel="u_k",title="No quantum tunneling on the half-line")
  Plots.savefig("noquantumtunnelingt=$(t).pdf")
end
plot!()

V8 = [zeros(10);-.5*ones(6);zeros(10);-.5*ones(6)]
  J8 = D + SymTriOperator(V8,[0.])
  n = 400
  uinit = pad!([zeros(15);.45],n)
  Jn = full(J8[1:n,1:n])
  t = 50
  ut8 = expm(im*Jn*t)*uinit
  plot(abs(ut8[1:100]),ylims=(-.5,1),marker=2,label="t=$(t)")
  plot!(pad(V8,100),ylims=(-.5,1),color=:black,label="Potential",xlabel="k",ylabel="u_k",title="Quantum tunneling on the half-line")
Plots.savefig("doublewellsimclose.pdf")

n = 1000
uinit = pad!([zeros(11);.05;.09;.09;.05],n)
  Jn = full(J9[1:n,1:n])
  plot(pad(V9,100),ylims=(0,1.5),marker=2,label="t=0",xlabel="k",ylabel="u_k",title="Quantum tunneling on the half-line")
  t = 1000000
  ut9 = expm(im*Jn*t)*uinit
  plot!(abs(ut9[1:100]),ylims=(0,1.5),marker=2,label="t=$(t)")
Plots.savefig("doublewellsimfar.pdf")

u0 = [zeros(10);1.;1.;1.;1.]
plot(pad(u0,100),marker=2)
for t = [5;10;15;20;25]
  ut = SchroEvo(u0U\(exp(im*x*t)*(U*u0))
  plot!(pad(abs(ut.coefficients),100),marker=2)
end
plot!()