text:
"""
```julia
lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3) -> ds
```
```math
\\begin{aligned}
\\dot{X} &= \\sigma(Y-X) \\\\
\\dot{Y} &= -XZ + \\rho X -Y \\\\
\\dot{Z} &= XY - \\beta Z
\\end{aligned}
```
The famous three dimensional system due to Lorenz [1], shown to exhibit
so-called "deterministic nonperiodic flow". It was originally invented to study a
simplified form of atmospheric convection.

Currently, it is most famous for its strange attractor (occuring at the default
parameters), which resembles a butterfly. For the same reason it is
also associated with the term "butterfly effect" (a term which Lorenz himself disliked)
even though the effect applies generally to dynamical systems.
Default values are the ones used in the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : E. N. Lorenz, J. atmos. Sci. **20**, pp 130 (1963)
"""
function lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(loop, u0, [σ, ρ, β], loop_jac)
end
function loop(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        return SVector{3}(du1, du2, du3)
    end
end
function loop_jac(u, p, t)
    @inbounds begin
        σ, ρ, β = p
        J = @SMatrix [-σ  σ  0;
        ρ - u[3]  (-1)  (-u[1]);
        u[2]   u[1]  -β]
        return J
    end
end

function lorenz_iip(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(liip, u0, [σ, ρ, β], liip_jac)
end
function liip(du, u, p, t)
    @inbounds begin
        σ = p[1]; ρ = p[2]; β = p[3]
        du[1] = σ*(u[2]-u[1])
        du[2] = u[1]*(ρ-u[3]) - u[2]
        du[3] = u[1]*u[2] - β*u[3]
        return nothing
    end
end
function liip_jac(J, u, p, t)
    @inbounds begin
    σ, ρ, β = p
    J[1,1] = -σ; J[1, 2] = σ; J[1,3] = 0
    J[2,1] = ρ - u[3]; J[2,2] = -1; J[2,3] = -u[1]
    J[3,1] = u[2]; J[3,2] = u[1]; J[3,3] = -β
    return nothing
    end
end

"""
```julia
chua(u0 = [0.7, 0.0, 0.0]; a = 15.6, b = 25.58, m0 = -8/7, m1 = -5/7)
```
```math
\\begin{aligned}
\\dot{x} &= \\alpha (y - h(x))\\\\
\\dot{y} &= x - y+z \\\\
\\dot{z} &= \\beta y
\\end{aligned}
```
where h(x) is defined by
```math
h(x) = m_1 x + \\frac 1 2 (m_0 - m_1)(|x + 1| - |x - 1|)
```
This is a 3D continuous system that exhibits chaos.

Chua designed an electronic circuit with the expressed goal of exhibiting
chaotic motion, and this system is obtained by rescaling the circuit units
to simplify the form of the equation. [1]

The parameters are a, b, m0 and m1. Setting a = 15.6, m0 = -8/7 and m1 = -5/7,
and varying the parameter b from b = 25 to b = 51, one observes a classic
period-doubling bifurcation route to chaos. [2]

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : Chua, Leon O. "The genesis of Chua's circuit", 1992.

[2] : [Leon O. Chua (2007) "Chua circuit", Scholarpedia, 2(10):1488.](http://www.scholarpedia.org/article/Chua_circuit)

"""
function chua(u0 = [0.7, 0.0, 0.0]; a = 15.6, b = 25.58, m0 = -8/7, m1 = -5/7)
    return CDS(chua_rule, u0, [a, b, m0, m1], chua_jacob)
end
function chua_rule(u, p, t)
    @inbounds begin
    a, b, m0, m1 = p
    du1 = a * (u[2] - u[1] - chua_element(u[1], m0, m1))
    du2 = u[1] - u[2] + u[3]
    du3 = -b * u[2]
    return SVector{3}(du1, du2, du3)
    end
end
function chua_jacob(u, p, t)
    a, b, m0, m1 = p
    return @SMatrix[-a*(1 + chua_element_derivative(u[1], m0, m1)) a 0;
                    1 -1 1;
                    0 -b 0]
end
# Helper functions for Chua's circuit.
function chua_element(x, m0, m1)
    return m1 * x + 0.5 * (m0 - m1) * (abs(x + 1.0) - abs(x - 1.0))
end
function chua_element_derivative(x, m0, m1)
    return m1 + 0.5 * (m0 - m1) * (-1 < x < 1 ? 2 : 0)
end

"""
```julia
roessler(u0=[1, -2, 0.1]; a = 0.2, b = 0.2, c = 5.7)
```
```math
\\begin{aligned}
\\dot{x} &= -y-z \\\\
\\dot{y} &= x+ay \\\\
\\dot{z} &= b + z(x-c)
\\end{aligned}
```
This three-dimensional continuous system is due to Rössler [1].
It is a system that by design behaves similarly
to the `lorenz` system and displays a (fractal)
strange attractor. However, it is easier to analyze qualitatively, as for example
the attractor is composed of a single manifold.
Default values are the same as the original paper.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : O. E. Rössler, Phys. Lett. **57A**, pp 397 (1976)
"""
function roessler(u0=[1, -2, 0.1]; a = 0.2, b = 0.2, c = 5.7)
    return CDS(roessler_rule, u0, [a, b, c], roessler_jacob)
end
function roessler_rule(u, p, t)
    @inbounds begin
    a, b, c = p
    du1 = -u[2]-u[3]
    du2 = u[1] + a*u[2]
    du3 = b + u[3]*(u[1] - c)
    return SVector{3}(du1, du2, du3)
    end
end
function roessler_jacob(u, p, t)
    a, b, c = p
    return @SMatrix [0.0 (-1.0) (-1.0);
                     1.0 a 0.0;
                     u[3] 0.0 (u[1]-c)]
end

"""
    double_pendulum(u0 = [π/2, 0, 0, 0.5];
                    G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
Famous chaotic double pendulum system (also used for our logo!). Keywords
are gravity (G), lengths of each rod and mass of each ball (all assumed SI units).

The variables order is [θ1, dθ1/dt, θ2, dθ2/dt].

Jacobian is created automatically (thus methods that use the Jacobian will be slower)!

(please contribute the Jacobian and the e.o.m. in LaTeX :smile:)

The parameter container has the parameters in the same order as stated in this
function's documentation string.
"""
function double_pendulum(u0=[π/2, 0, 0, 0.5]; G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
    return CDS(doublependulum_rule, u0, [G, L1, L2, M1, M2])
end
@inbounds function doublependulum_rule(state, p, t)
    G, L1, L2, M1, M2 = p

    du1 = state[2]
    del_ = state[3] - state[1]
    den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
    du2 = (M2*L1*state[2]*state[2]*sin(del_)*cos(del_) +
               M2*G*sin(state[3])*cos(del_) +
               M2*L2*state[4]*state[4]*sin(del_) -
               (M1 + M2)*G*sin(state[1]))/den1

    du3 = state[4]

    den2 = (L2/L1)*den1
    du4 = (-M2*L2*state[4]*state[4]*sin(del_)*cos(del_) +
               (M1 + M2)*G*sin(state[1])*cos(del_) -
               (M1 + M2)*L1*state[2]*state[2]*sin(del_) -
               (M1 + M2)*G*sin(state[3]))/den2
    return SVector{4}(du1, du2, du3, du4)
end

"""
    henonheiles(u0=[0, -0.25, 0.42081,0])
```math
\\begin{aligned}
\\dot{x} &= p_x \\\\
\\dot{y} &= p_y \\\\
\\dot{p}_x &= -x -2 xy \\\\
\\dot{p}_y &= -y - (x^2 - y^2)
\\end{aligned}
```

The Hénon–Heiles system [1] is a conservative dynamical system and was introduced as a simplification of the motion
of a star around a galactic center. It was originally intended to study the
existence of a "third integral of motion" (which would make this 4D system integrable).
In that search, the authors encountered chaos, as the third integral existed
for only but a few initial conditions.

The default initial condition is a typical chaotic orbit.
The function `Systems.henonheiles_ics(E, n)` generates a grid of
`n×n` initial conditions, all having the same energy `E`.

[1] : Hénon, M. & Heiles, C., The Astronomical Journal **69**, pp 73–79 (1964)
"""
function henonheiles(u0=[0, -0.25, 0.42081, 0]#=; conserveE::Bool = true=#)
    i = one(eltype(u0))
    o = zero(eltype(u0))
    J = zeros(eltype(u0), 4, 4)
    return CDS(hhrule!, u0, nothing, hhjacob!, J)
end
function hhrule!(du, u, p, t)
    @inbounds begin
        du[1] = u[3]
        du[2] = u[4]
        du[3] = -u[1] - 2u[1]*u[2]
        du[4] = -u[2] - (u[1]^2 - u[2]^2)
        return nothing
    end
end
function hhjacob!(J, u, p, t)
    @inbounds begin
    o = 0.0; i = 1.0
    J[1,:] .= (o,    o,     i,    o)
    J[2,:] .= (o,    o,     o,    i)
    J[3,:] .= (-i - 2*u[2],   -2*u[1],   o,   o)
    J[4,:] .= (-2*u[1],  -1 + 2*u[2],  o,   o)
    return nothing
    end
end
_hhenergy(x,y,px,py) = 0.5(px^2 + py^2) + _hhpotential(x,y)
_hhpotential(x, y) = 0.5(x^2 + y^2) + (x^2*y - (y^3)/3)
function henonheiles_ics(E, n=10)
    ys = range(-0.4, stop = 1.0, length = n)
    pys = range(-0.5, stop = 0.5, length = n)
    ics = Vector{Vector{Float64}}()
    for y in ys
        V = _hhpotential(0.0, y)
        V ≥ E && continue
        for py in pys
            Ky = 0.5*(py^2)
            Ky + V ≥ E && continue
            px = sqrt(2(E - V - Ky))
            ic = [0.0, y, px, py]
            push!(ics, [0.0, y, px, py])
        end
    end
    return ics
end

"""
    qbh([u0]; A=1.0, B=0.55, D=0.4)

A conservative dynamical system with rule
```math
\\begin{aligned}
\\dot{q}_0 &= A p_0 \\\\
\\dot{q}_2 &= A p_2 \\\\
\\dot{p}_0 &= -A q_0 -3 \\frac{B}{\\sqrt{2}} (q_2^2 - q_1^2) - D q_1 (q_1^2 + q_2^2) \\\\
\\dot{p}_2 &= -q_2 (A + 3\\sqrt{2} B q_1 + D (q_1^2 + q_2^2)) (x^2 - y^2)
\\end{aligned}
```

This dynamical rule corresponds to a Hamiltonian used in nuclear
physics to study the quadrupole vibrations of the nuclear surface [1,2].

```math
H(p_0, p_2, q_0, q_2) = \\frac{A}{2}\\left(p_0^2+p_2^2\\right)+\\frac{A}{2}\\left(q_0^2+q_2^2\\right)
			 +\\frac{B}{\\sqrt{2}}q_0\\left(3q_2^2-q_0^2\\right) +\\frac{D}{4}\\left(q_0^2+q_2^2\\right)^2
```

The Hamiltonian has a similar structure with the Henon-Heiles one, but it has an added fourth order term
and presents a nontrivial dependence of chaoticity with the increase of energy [3].
The default initial condition is chaotic.

[1]: Eisenberg, J.M., & Greiner, W., Nuclear theory 2 rev ed. Netherlands: North-Holland pp 80 (1975)

[2]: Baran V. and Raduta A. A., International Journal of Modern Physics E, **7**, pp 527--551 (1998)

[3]: Micluta-Campeanu S., Raportaru M.C., Nicolin A.I., Baran V., Rom. Rep. Phys. **70**, pp 105 (2018)
"""
function qbh(u0=[0., -2.5830294658973876, 1.3873470962626937, -4.743416490252585];  A=1., B=0.55, D=0.4)
    return CDS(qrule, u0, [A, B, D])
end
function qrule(z, p, t)
    @inbounds begin
        A, B, D = p
        p₀, p₂ = z[1], z[2]
        q₀, q₂ = z[3], z[4]

        return SVector{4}(
            -A * q₀ - 3 * B / √2 * (q₂^2 - q₀^2) - D * q₀ * (q₀^2 + q₂^2),
            -q₂ * (A + 3 * √2 * B * q₀ + D * (q₀^2 + q₂^2)),
            A * p₀,
            A * p₂
        )
    end
end

"""
    lorenz96(N::Int, u0 = rand(M); F=0.01)

```math
\\frac{dx_i}{dt} = (x_{i+1}-x_{i-2})x_{i-1} - x_i + F
```

`N` is the chain length, `F` the forcing. Jacobian is created automatically.
(parameter container only contains `F`)
"""
function lorenz96(N::Int, u0 = range(0; length = N, step = 0.1); F=0.01)
    @assert N ≥ 4 "`N` must be at least 4"
    lor96 = Lorenz96{N}() # create struct
    return CDS(lor96, u0, [F])
end
struct Lorenz96{N} end # Structure for size type
function (obj::Lorenz96{N})(dx, x, p, t) where {N}
    F = p[1]
    # 3 edge cases
    @inbounds dx[1] = (x[2] - x[N - 1]) * x[N] - x[1] + F
    @inbounds dx[2] = (x[3] - x[N]) * x[1] - x[2] + F
    @inbounds dx[N] = (x[1] - x[N - 2]) * x[N - 1] - x[N] + F
    # then the general case
    for n in 3:(N - 1)
      @inbounds dx[n] = (x[n + 1] - x[n - 2]) * x[n - 1] - x[n] + F
    end
    return nothing
end

"""
    duffing(u0 = [0.1, 0.25]; ω = 2.2, f = 27.0, d = 0.2, β = 1)
The (forced) duffing oscillator, that satisfies the equation
```math
\\ddot{x} + d\\cdot\\dot{x} + β*x + x^3 = f\\cos(\\omega t)
```
with `f, ω` the forcing strength and frequency and `d` the dampening.

The parameter container has the parameters in the same order as stated in this
function's documentation string.
"""
function duffing(u0 = [0.1, 0.25]; ω = 2.2, f = 27.0, d = 0.2, β = 1)
    J = zeros(eltype(u0), 2, 2)
    J[1,2] = 1
    return CDS(duffing_rule, u0, [ω, f, d, β], duffing_jacob, J)
end
@inbounds function duffing_rule(x, p, t)
    ω, f, d, β = p
    dx1 = x[2]
    dx2 = f*cos(ω*t) - β*x[1] - x[1]^3 - d * x[2]
    return SVector(dx1, dx2)
end
@inbounds function duffing_jacob(u, p, t)
    ω, f, d, β = p
    return @SMatrix [0 1 ;
    (-β - 3u[1]^2) -d]
end

"""
    shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)
Shinriki oscillator with all other parameters (besides `R1`) set to constants.
*This is a stiff problem, be careful when choosing solvers and tolerances*.
"""
function shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)
    # # Jacobian caller for Shinriki:
    # shinriki_rule(::Type{Val{:jac}}, J, u, p, t) = (shi::Shinriki)(t, u, J)
    return CDS(shinriki_rule, u0, [R1])
end
shinriki_voltage(V) = 2.295e-5*(exp(3.0038*V) - exp(-3.0038*V))
function shinriki_rule(u, p, t)
    R1 = p[1]

    du1 = (1/0.01)*(
    u[1]*(1/6.9 - 1/R1) - shinriki_voltage(u[1] - u[2]) - (u[1] - u[2])/14.5
    )

    du2 = (1/0.1)*(
    shinriki_voltage(u[1] - u[2]) + (u[1] - u[2])/14.5 - u[3]
    )

    du3 = (1/0.32)*(-u[3]*0.1 + u[2])
    return SVector{3}(du1, du2, du3)
end

"""
```julia
gissinger(u0 = [3, 0.5, 1.5]; μ = 0.119, ν = 0.1, Γ = 0.9)
```
```math
\\begin{aligned}
\\dot{Q} &= \\mu Q - VD \\\\
\\dot{D} &= -\\nu D + VQ \\\\
\\dot{V} &= \\Gamma -V + QD
\\end{aligned}
```
A continuous system that models chaotic reversals due to Gissinger [1], applied
to study the reversals of the magnetic field of the Earth.

The parameter container has the parameters in the same order as stated in this
function's documentation string.

[1] : C. Gissinger, Eur. Phys. J. B **85**, 4, pp 1-12 (2012)
"""
function gissinger(u0 = [3, 0.5, 1.5]; μ = 0.119, ν = 0.1, Γ = 0.9)
    return CDS(gissinger_rule, u0, [μ, ν, Γ], gissinger_jacob)
end
function gissinger_rule(u, p, t)
    μ, ν, Γ = p
    du1 = μ*u[1] - u[2]*u[3]
    du2 = -ν*u[2] + u[1]*u[3]
    du3 = Γ - u[3] + u[1]*u[2]
    return SVector{3}(du1, du2, du3)
end
function gissinger_jacob(u, p, t)
    μ, ν, Γ = p
    return @SMatrix [μ -u[3] -u[2];
                     u[3] -ν u[1];
                     u[2] u[1] -1]
end

"""
```julia
rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)
```
```math
\\begin{aligned}
\\dot{x} &= -\\mu x +yz \\\\
\\dot{y} &= -\\mu y +x(z-\\alpha) \\\\
\\dot{z} &= 1 - xz
\\end{aligned}
```
Rikitake's dynamo is a system that tries to model the magnetic reversal events
by means of a double-disk dynamo system.

[1] : T. Rikitake Math. Proc. Camb. Phil. Soc. **54**, pp 89–105, (1958)
"""
function rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)
    return CDS(rikitake_rule, u0, [μ, α], rikitake_jacob)
end
function rikitake_rule(u, p, t)
    μ, α = p
    x,y,z = u
    xdot = -μ*x + y*z
    ydot = -μ*y + x*(z - α)
    zdot = 1 - x*y
    return SVector{3}(xdot, ydot, zdot)
end
function rikitake_jacob(u, p, t)
    μ, α = p
    x,y,z = u
    xdot = -μ*x + y*z
    ydot = -μ*y + x*(z - α)
    zdot = 1 - x*y
    return @SMatrix [-μ  z  y;
                     z-α -μ x;
                     -y  -x 0]
end

"""
```julia
nosehoover(u0 = [0, 0.1, 0])
```
```math
\\begin{aligned}
\\dot{x} &= y \\\\
\\dot{y} &= yz - x \\\\
\\dot{z} &= 1 - y^2
\\end{aligned}
```
Three dimensional conservative continuous system, discovered in 1984 during
investigations in thermodynamical chemistry by Nosé and Hoover, then
rediscovered by Sprott during an exhaustive search as an extremely simple
chaotic system. [1]

See Chapter 4 of "Elegant Chaos" by J. C. Sprott. [2]

[1] : Hoover, W. G. (1995). Remark on ‘‘Some simple chaotic flows’’. *Physical Review E*, *51*(1), 759.

[2] : Sprott, J. C. (2010). *Elegant chaos: algebraically simple chaotic flows*. World Scientific.
"""
nosehoover(u0 = [0, 0.1, 0]) = CDS(nosehoover_rule, u0, nothing, nosehoover_jacob)
function nosehoover_rule(u, p, t)
    x,y,z = u
    xdot = y
    ydot = y*z - x
    zdot  = 1.0 - y*y
    return SVector{3}(xdot, ydot, zdot)
end
function nosehoover_jacob(u, p, t)
    x,y,z = u
    return @SMatrix [0 1 0;
                     -1 z y;
                     0 2y 0]
end

"""
    antidots([u]; B = 1.0, d0 = 0.3, c = 0.2)
An antidot "superlattice" is a Hamiltonian system that corresponds to a
smoothened periodic Sinai billiard with disk diameter `d0` and smooth
factor `c` [1].

This version is the two dimensional
classical form of the system, with quadratic dynamical rule and
a perpendicular magnetic field. Notice that the dynamical rule
is with respect to the velocity instead of momentum, i.e.:
```math
\\begin{aligned}
\\dot{x} &= v_x \\\\
\\dot{y} &= v_y \\\\
\\dot{v_x} &= B*v_y - U_x \\\\
\\dot{v_y} &= -B*v_x - U_X \\\\
\\end{aligned}
```
with ``U`` the potential energy:
```math
U = \\left(\\tfrac{1}{c^4}\\right) \\left[\\tfrac{d_0}{2} + c - r_a\\right]^4
```
if ``r_a = \\sqrt{(x \\mod 1)^2 + (y \\mod 1)^2} < \\frac{d_0}{2} + c`` and 0
otherwise. I.e. the potential is periodic with period 1 in both ``x, y`` and
normalized such that for energy value of 1 it is a circle of diameter ``d0``.
The magnetic field is also normalized such that for value `B=1` the cyclotron
diameter is 1.

[1] : G. Datseris *et al*, [New Journal of Physics 2019](https://iopscience.iop.org/article/10.1088/1367-2630/ab19cc/meta)
"""
function antidots(u0 = [0.5, 0.5, 0.25, 0.25];
    d0 = 0.5, c = 0.2, B = 1.0)
    return CDS(antidot_rule, u0, [B, d0, c], antidot_jacob)
end

function antidot_rule(u, p, t)
    B, d0, c = p
    x, y, vx, vy = u
    # Calculate quadrant of (x,y):
    U = Uy = Ux = 0.0
    xtilde = x - round(x);  ytilde = y - round(y)
    ρ = sqrt(xtilde^2 + ytilde^2)
    # Calculate derivatives and potential:
    if ρ < 0.5*d0 + c
        sharedfactor = -(4*(c + d0/2 - ρ)^(3))/(c^4*ρ)
        Ux = sharedfactor*xtilde # derivatives
        Uy = sharedfactor*ytilde
    end
    Br = 2*√(2)*B # magnetic field with prefactor
    return SVector{4}(vx, vy, Br*vy - Ux, -vx*Br - Uy)
end
function antidot_jacob(u, p, t)
    B, d0, c = p
    x, y, vx, vy = u
    xtilde = x - round(x);  ytilde = y - round(y)
    ρ = sqrt(xtilde^2 + ytilde^2)
    # Calculate derivatives and potential:
    if ρ < 0.5*d0 + c
        Uxx, Uyy, Uxy = antidot_secondderv(xtilde, ytilde, p)
    else
        Uxx, Uyy, Uxy = 0.0, 0.0, 0.0
    end
    Br = 2*√(2)*B # magnetic field with prefactor
    return SMatrix{4, 4}(0.0, 0, -Uxx, -Uxy, 0, 0, -Uxy, -Uyy,
                         1, 0, 0, -Br,        0, 1, Br, 0)
end

function antidot_potential(x::Real, y::Real, p)
    B, d0, c = p
    δ = 4
    # Calculate quadrant of (x,y):
    xtilde = x - round(x);  ytilde = y - round(y)
    ρ = sqrt(xtilde^2 + ytilde^2)
    # Check distance:
    pot = ρ > 0.5*d0 + c ? 0.0 : (1/(c^δ))*(0.5*d0 + c - ρ)^δ
end

function antidot_secondderv(x, y, p)
    B, d0, c = p
    r = sqrt(x^2 + y^2)
    Uxx = _Uxx(x, y, c, d0, r)
    Uyy = _Uxx(y, x, c, d0, r)
    Uxy =  (2c + d0 - 2r)^2 * (2c + d0 + 4r)*x*y / (2c^4*r^3)
    return Uxx, Uyy, Uxy
end

function _Uxx(x, y, c, d0, r)
    Uxx =  (2c + d0 - 2r)^2 * (-2c*y^2 - d0*y^2 + 2*r*(3x^2 + y^2)) /
    (2 * (c^4) * r^3)
end

"""
```julia
ueda(u0 = [3.0, 0]; k = 0.1, B = 12.0)
```
```math
\\ddot{x} + k \\dot{x} + x^3 = B\\cos{t}
```
Nonautonomous Duffing-like forced oscillation system, discovered by Ueda in
1961. It is one of the first chaotic systems to be discovered.

The stroboscopic plot in the (x, ̇x) plane with period 2π creates a "broken-egg
attractor" for k = 0.1 and B = 12. Figure 5 of [1] is reproduced by

```julia
using Plots
ds = Systems.ueda()
a = trajectory(ds, 2π*5e3, dt = 2π)
scatter(a[:, 1], a[:, 2], markersize = 0.5, title="Ueda attractor")
```

For more forced oscillation systems, see Chapter 2 of "Elegant Chaos" by
J. C. Sprott. [2]

[1] : [Ruelle, David, ‘Strange Attractors’, The Mathematical Intelligencer, 2.3 (1980), 126–37](https://doi.org/10/dkfd3n)

[2] : Sprott, J. C. (2010). *Elegant chaos: algebraically simple chaotic flows*. World Scientific.
"""
function ueda(u0 = [3.0, 0]; k = 0.1, B = 12.0)
    return CDS(ueda_rule, u0, [k, B], ueda_jacob)
end
function ueda_rule(u, p, t)
    x,y = u
    k, B = p
    xdot = y
    ydot = B*cos(t) - k*y - x^3
    return SVector{2}(xdot, ydot)
end
function ueda_jacob(u, p, t)
    x,y = u
    k, B = p
    return @SMatrix [0      1;
                     -3*x^2 -k]
end

struct MagneticPendulum
    magnets::Vector{SVector{2, Float64}}
end
mutable struct MagneticPendulumParams
    γs::Vector{Float64}
    d::Float64
    α::Float64
    ω::Float64
end

function (m::MagneticPendulum)(u, p, t)
    x, y, vx, vy = u
    γs::Vector{Float64}, d::Float64, α::Float64, ω::Float64 = p.γs, p.d, p.α, p.ω
    dx, dy = vx, vy
    dvx, dvy = @. -ω^2*(x, y) - α*(vx, vy)
    for (i, ma) in enumerate(m.magnets)
        δx, δy = (x - ma[1]), (y - ma[2])
        D = sqrt(δx^2 + δy^2 + d^2)
        dvx -= γs[i]*(x - ma[1])/D^3
        dvy -= γs[i]*(y - ma[2])/D^3
    end
    return SVector(dx, dy, dvx, dvy)
end

"""
    magnetic_pendulum(u=[0.7,0.7,0,0]; d=0.3, α=0.2, ω=0.5, N=3, γs=fill(1.0,N))

Create a pangetic pendulum with `N` magnetics, equally distributed along the unit circle,
with dynamical rule
```math
\\begin{aligned}
\\ddot{x} &= -\\omega ^2x - \\alpha \\dot{x} - \\sum_{i=1}^N \\frac{\\gamma_i (x - x_i)}{D_i^3} \\\\
\\ddot{y} &= -\\omega ^2y - \\alpha \\dot{y} - \\sum_{i=1}^N \\frac{\\gamma_i (y - y_i)}{D_i^3} \\\\
D_i &= \\sqrt{(x-x_i)^2  + (y-y_i)^2 + d^2}
\\end{aligned}
```
where α is friction, ω is eigenfrequency, d is distance of pendulum from the magnet's plane
and γ is the magnetic strength.
"""
function magnetic_pendulum(u = [sincos(0.12553*2π)..., 0, 0];
    γ = 1.0, d = 0.3, α = 0.2, ω = 0.5, N = 3, γs = fill(γ, N))
    m = MagneticPendulum([SVector(cos(2π*i/N), sin(2π*i/N)) for i in 1:N])
    p = MagneticPendulumParams(γs, d, α, ω)
    ds = ContinuousDynamicalSystem(m, u, p)
end

"""
    fitzhugh_nagumo(u = 0.5ones(2); a=3.0, b=0.2, ε=0.01, I=0.0)
Famous excitable system which emulates the firing of a neuron, with rule
```math
\\begin{aligned}
\\dot{v} &= av(v-b)(1-v) - w + I \\\\
\\ddot{w} &= \\varepsilon(v - w)
\\end{aligned}
```

More details in the [Scholarpedia](http://www.scholarpedia.org/article/FitzHugh-Nagumo_model) entry.
"""
function fitzhugh_nagumo(u = 0.5ones(2); a=3.0, b=0.2, ε=0.01, I=0.0)
    ds = ContinuousDynamicalSystem(fitzhugh_nagumo_rule, u, [a, b, ε, I])
end
function fitzhugh_nagumo_rule(x, p, t)
    u, w = x
    a, b, ε, I = p
    return SVector(a*u*(u-b)*(1. - u) - w + I, ε*(u - w))
end

"""
    more_chaos_example(u = rand(3))
A three dimensional chaotic system introduced in [^Sprott2020] with rule
```math
\\begin{aligned}
\\dot{x} &= y \\\\
\\dot{y} &= -x - sign(z)y \\\\
\\dot{z} &= y^2 - \\exp(-x^2)
\\end{aligned}
```
It is noteworthy because its strange attractor is multifractal with fractal dimension ≈ 3.

[^Sprott2020]: Sprott, J.C. 'Do We Need More Chaos Examples?', Chaos Theory and Applications 2(2),1-3, 2020
"""
more_chaos_example(u = [0.0246, 0.79752, 0.3535866]) =
ContinuousDynamicalSystem(more_chaos_rule, u, nothing)
function more_chaos_rule(u, p, t)
    x, y, z = u
    dx = y
    dy = -x - sign(z)*y
    dz = y^2 - exp(-x^2)
    return SVector(dx, dy, dz)
end

"""
    thomas_cyclical(u0 = [1.0, 0, 0]; b = 0.2)
```math
\\begin{aligned}
\\dot{x} &= \\sin(y) - bx\\\\
\\dot{y} &= \\sin(z) - by\\\\
\\dot{z} &= \\sin(x) - bz
\\end{aligned}
```
Thomas' cyclically symmetric attractor is a 3D strange attractor originally proposed
by René Thomas[^Thomas1999]. It has a simple form which is cyclically symmetric in the
x,y, and z variables and can be viewed as the trajectory of a frictionally dampened
particle moving in a 3D lattice of forces.
For more see the [Wikipedia page](https://en.wikipedia.org/wiki/Thomas%27_cyclically_symmetric_attractor).

Reduces to the labyrinth system for `b=0`, see
See discussion in Section 4.4.3 of "Elegant Chaos" by J. C. Sprott.

[^Thomas1999]: Thomas, R. (1999). *International Journal of Bifurcation and Chaos*, *9*(10), 1889-1905.
"""
thomas_cyclical(u0 = [1.0, 0, 0]; b = 0.2) = CDS(thomas_rule, u0, [b], thomas_jacob)
labyrinth(u0 = [1.0, 0, 0]) = CDS(thomas_rule, u0, [0.0], thomas_jacob)

function thomas_rule(u, p, t)
    x,y,z = u
    b = p[1]
    xdot = sin(y) - b*x
    ydot = sin(z) - b*y
    zdot = sin(x) - b*z
    return SVector{3}(xdot, ydot, zdot)
end
function thomas_jacob(u, p, t)
    x,y,z = u
    b = p[1]
    return @SMatrix [b cos(y) 0;
                     0 b cos(z);
                     cos(x) 0 b]
end

"""
    stommel_thermohaline(u = [0.3, 0.2]; η1 = 3.0, η2 = 1, η3 = 0.3)
Stommel's box model for Atlantic thermohaline circulation
```math
\\begin{aligned}
 \\dot{T} &= \\eta_1 - T - |T-S|\\cdot T \\\\
 \\dot{S} &= \\eta_2 - \\eta_3S - |T-S|\\cdot S
\\end{aligned}
```
Here ``T, S`` are variables standing for dimensionless temperature and salinity
differences between the boxes (polar and equitorial ocean basins) and ``\\eta_i``
are parameters.

[^Stommel1961]: Stommel, Thermohaline convection with two stable regimes offlow. Tellus, 13(2)
"""
function stommel_thermohaline(u = [0.3, 0.2]; η1 = 3.0, η2 = 1, η3 = 0.3)
    ds = ContinuousDynamicalSystem(stommel_thermohaline_rule, u, [η1, η2, η3],
    stommel_thermohaline_jacob)
end
function stommel_thermohaline_rule(x, p, t)
    T, S = x
    η1, η2, η3 = p
    q = abs(T-S)
    return SVector(η1 - T -q*T, η2 - η3*S - q*S)
end
function stommel_thermohaline_jacob(x, p, t)
    T, S = x
    η1, η2, η3 = p
    q = abs(T-S)
    if T ≥ S
        return @SMatrix [(-1 - 2T + S)  (T);
                         (-S)  (-η3 - T + 2S)]
    else
        return @SMatrix [(-1 + 2T - S)  (-T);
                         (+S)  (-η3 + T - 2S)]
    end
end


text:
using DifferentialEquations
using ProgressMeter
using PyCall
using MATLAB
using Printf

# pygui(:qt)
using PyPlot
include("PSS_Definitions_Dimer_X.jl")
# pygui(true);

H=.1
Energy=H

# H=range(.1505,stop=0.2505,length=9)

# @showprogress 1 "Computing..." for Energy in H
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/"
h=replace(@sprintf("%.15f",Energy),"."=>"_")
h_title=@sprintf("%.8f",Energy)
file_name=location*"DIMER_PSS_X"*h*".fig"
mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
mat"xlabel('Q')"
mat"ylabel('P')"
mat"title($h_title)"
mat"axis([-3, 3, -3, 3])"
# Q1=Q1_find_dimer(.2,.4,H)
# t_end=1000.0;
# q0,p0=[zeros(2) for i in 1:2]
# q0[1]=Q1;
# q0[2]=.2;
# p0[1]=0;
# p0[2]=.4;
# Hamiltonian_Dimer(q0,p0,1)
# prob= HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
# @time t,A=solve(prob, Vern9(),reltol=1e-13,abstol=1e-15,maxiters=1e15);

# q1=A[1,:]
# q2=A[2,:]
# p1=A[3,:]
# p2=A[4,:];

# Q1=(q1+q2)/sqrt(2)
# Q2=(q1-q2)/sqrt(2)
# P1=(p1+p2)/sqrt(2)
# P2=(p1-p2)/sqrt(2);

# figure()
# plot(q1,p1,"b",linewidth=.1)
# plot(q2,p2,"r",linewidth=.1);

# figure()
# plot(Q1,P1,"b",linewidth=.1)
# plot(Q2,P2,"r",linewidth=.1);

t_end=1e5;
# figure()
# @time Q_PSS,P_PSS=PSS_function( q0[2], p0[2], H, t_end,count);
# plot(Q_PSS,P_PSS,".", markersize=4, c
#     ="k");

N_iter_Q=200;#50
Q_start=.5
Q_end=3
N_iter_P=20
P_start=0.2
P_end=1.4
ArrP=range(P_start,stop=P_end,length=N_iter_P)
ArrQ=range(Q_start,stop=Q_end,length=N_iter_Q)
Brandons_Colors=["#393b79" ,"#5254a3","#6b6ecf","#9c9ede" ,"#637939","#8ca252" ,"#b5cf6b" ,"#cedb9c" ,"#8c6d31","#bd9e39" ,"#e7ba52","#e7cb94","#843c39","#ad494a" ,"#d6616b","#e7969c" ,"#7b4173" ,"#a55194","#ce6dbd" ,"#de9ed6"];

# figure()
@showprogress for k=1:N_iter_Q
    Q_PSS,P_PSS=PSS_function(ArrQ[k], 0, Energy, t_end);
    current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
    if Q_PSS!=0
        mat"plot($Q_PSS,$P_PSS,'.','MarkerSize',3,'color',$current_color); hold on;"
        mat"plot(-$Q_PSS,-$P_PSS,'.','MarkerSize',3,'color',$current_color); hold on;"

        # plot(Q_PSS,P_PSS,".", markersize=1,c=current_color);
    end
end
@showprogress for k=1:N_iter_P
    Q_PSS,P_PSS=PSS_function(0,ArrP[k], Energy, t_end);
    current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
    if Q_PSS!=0
        mat"plot($Q_PSS,$P_PSS,'.','MarkerSize',3,'color',$current_color); hold on;"
        mat"plot(-$Q_PSS,-$P_PSS,'.','MarkerSize',3,'color',$current_color); hold on;"
    end
end

mat"savefig($file_name)"

# end

ext: jl


text:
import VimBindings: yank, cut

@testset "yanking" begin
    buf = IOBuffer("0123456789")

    @test cut(buf, Motion(0, 3)) == "012"
    @test yank(buf, Motion(0, 3)) == "012"
    @test VB.vim.registers['"'] == "012"

    VB.vim.register = 'a'
    yank(buf, Motion(3, 6))
    @test VB.vim.registers['a'] == "345"

    VB.vim.register = '_'
    @test yank(buf, Motion(3, 6)) === nothing

end


text:
# Removing t = 0, such that Σ is invertible
n = 500.0;
t = Vector(0.1:1/n:1)

# Creating a test matrix Σ = tril(UV') + triu(VU',1) that is PSD
p = 2;
Ut, Vt = spline_kernel(t', p)
K  = SymEGRQSMatrix(Ut,Vt,ones(size(Ut,2)))
Σ    = Matrix(K)
chol = cholesky(Σ)
L = cholesky(K)
x = randn(size(K,1))

# Testing inverses (Using Cholesky factorizations)
B = randn(length(t),10);
@test L\B  ≈ chol.L\B
@test L'\B ≈ chol.U\B
@test L*B  ≈ chol.L*B
@test L'*B ≈ chol.U*B

# Testing logdet
@test logdet(L) ≈ logdet(chol.L)
@test det(L) ≈ det(chol.L)

# Testing traces and norm
M = SymEGRSSMatrix(Ut,Vt);
@test isapprox(tr(L,M), tr(chol\Matrix(M)), atol=1e-6)
@test trinv(L) ≈ tr(chol\Diagonal(ones(size(L,1))))
@test SymEGRSSMatrices.fro_norm_L(L) ≈ norm(chol.L[:])^2

# Testing show
@test L.L ≈ tril(Ut'*L.Wt,-1) + Diagonal(L.d)
@test L.U ≈ triu(L.Wt'*Ut,1) + Diagonal(L.d)
@test Matrix(L) ≈ tril(L.Ut'*L.Wt,-1) + Diagonal(L.d)
@test L[3,1] ≈ chol.L[3,1]
@test L[2,2] ≈ chol.L[2,2]
@test L[1,3] ≈ chol.L[1,3]

# Testing explicit-implicit-inverse
Yt, Zt = SymEGRSSMatrices.dss_create_yz(L.Ut,L.Wt,L.d)
@test tril(Yt'*Zt,-1) + Diagonal(L.d.^(-1)) ≈ inv(chol.L)
@test L*(tril(Yt'*Zt,-1) + Diagonal(L.d.^(-1))) ≈ I
@test L'*(triu(Zt'*Yt,1) + Diagonal(L.d.^(-1))) ≈ I
@test SymEGRSSMatrices.squared_norm_cols(Yt,Zt,L.d.^(-1)) ≈ sum(inv(chol.L).^2,dims=1)'

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
module StringMatch

include("sais.jl")
include("lcparray.jl")
include("suffixarray.jl")

end


text:
# ---
# title: 1344. Angle Between Hands of a Clock
# id: problem1344
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/angle-between-hands-of-a-clock/description/>
# hidden: true
# ---
#
# Given two numbers, `hour` and `minutes`. Return the smaller angle (in degrees)
# formed between the `hour` and the `minute` hand.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/26/sample_1_1673.png)
#
#
#
#     Input: hour = 12, minutes = 30
#     Output: 165
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/12/26/sample_2_1673.png)
#
#
#
#     Input: hour = 3, minutes = 30
#     Output: 75
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/12/26/sample_3_1673.png)**
#
#
#
#     Input: hour = 3, minutes = 15
#     Output: 7.5
#
#
# **Example 4:**
#
#
#
#     Input: hour = 4, minutes = 50
#     Output: 155
#
#
# **Example 5:**
#
#
#
#     Input: hour = 12, minutes = 0
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= hour <= 12`
#   * `0 <= minutes <= 59`
#   * Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Select

using GaPLAC
using GaPLAC.CSV
using GaPLAC.DataFrames
using GaPLAC.StatsBase
using GaPLAC.AbstractGPs

function run(args)
    @info "running 'select'"
    @info args

    !isempty(args["chains"]) && !isempty(args["formulae"]) && throw(ArgumentError("'select' can only take one of '--formulae' or '--chains', not both"))

    if !isempty(args["chains"])
        lp1df = CSV.read(args["chains"][1], DataFrame)
        lp1 = log2(harmmean([BigFloat(2) ^ x for x in lp1df[!, :lp]]))
        lp2df = CSV.read(args["chains"][2], DataFrame)
        lp2 = log2(harmmean([BigFloat(2) ^ x for x in lp2df[!, :lp]]))
        bayes = log2(BigFloat(2)^lp1 / BigFloat(2)^lp2)
    elseif !isempty(args["formulae"])
        spec1 = GaPLAC.gp_formula(args["formulae"][1])
        spec2 = GaPLAC.gp_formula(args["formulae"][2])

        @debug "GP formulae" GaPLAC.formula(spec1) GaPLAC.formula(spec2)

        eq1, vars1 = GaPLAC._apply_vars(GaPLAC.formula(spec1))
        kernels1 = GaPLAC._walk_kernel(eq1)
        length(vars1) == length(kernels1) || error("Something went wrong with equation parsing, number of variables should == number of kernels")

        eq2, vars2 = GaPLAC._apply_vars(GaPLAC.formula(spec2))
        kernels2 = GaPLAC._walk_kernel(eq2)
        length(vars2) == length(kernels2) || error("Something went wrong with equation parsing, number of variables should == number of kernels")

        gp1 = GP(eq1)
        gp2 = GP(eq2)

        @debug "GPs:" gp1 gp2

        df = CSV.read(args["data"], DataFrame)
        df = disallowmissing(df[completecases(df),:])

        y1 = df[!, GaPLAC.response(spec1)]
        x1 = Matrix(df[!, vars1])
        pr1 = AbstractGPs.FiniteGP(gp1, x1, 0.1, obsdim=1)

        y2 = df[!, GaPLAC.response(spec2)]
        x2 = Matrix(df[!, vars2])
        pr2 = AbstractGPs.FiniteGP(gp2, x2, 0.1, obsdim=1)

        lp1 = logpdf(pr1, y1)
        lp2 = logpdf(pr2, y2)
        pst1 = posterior(pr1, y1)
        pst2 = posterior(pr2, y2)

        bayes = log2(BigFloat(2)^lp1 / BigFloat(2)^lp2)
    else
        throw(ArgumentError("'select' command requires either '--chains' or '--formulae' arguments"))
    end

    @info """
        **Log2 Bayes**: $(round(Float64(bayes), digits=3))

        - **Log(pdf)** - model 1: $(round(lp1, digits=4))
        - **Log(pdf)** - model 2: $(round(lp2, digits=4))

        _Note_ - Positive values indicate more evidence for model 1
        """

end

end

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
module StructuredMechModels

using Flux
using LinearAlgebra
using TrajectoryOptimization
using DocStringExtensions

import TrajectoryOptimization: Model

include("utils.jl")

abstract type AbstractRigidBody end

"""
$(TYPEDEF)
"""
struct Env
    sys
    x0
    thetamask
end

Env(sys, x0) = Env(sys, x0, zeros(size(x0)))

"""
Helper function to interface parameterized rigidbody dynamics functions with TrajectoryOptimization
"""
function rigid_body_dynamics!(ẋ, x, u, p)
    res = p.f(x, u)
    ẋ[:] = res
end

function rigid_body_jacobian!(Z, x, u, p)
    error("Shouldn't be called! TrajOpt uses ForwardDiff to generate the jacobian")
end

"""
Converts parameterized rigidbody dynamics functions to `Model` useful for `TrajectoryOptimization`

$(SIGNATURES)
"""
function Model(nn::AbstractRigidBody)
    p = (f=nn,)
    return AnalyticalModel{Nominal,TrajectoryOptimization.Continuous}(
        (xdot, x, u)->rigid_body_dynamics!(xdot, x, u, p),
        (Z,x,u)->rigid_body_jacobian!(Z, x, u, p), 2*nn.qdim, nn.udim, 0, p, Dict{Symbol,Any}())
end

include("naive.jl")

"""
Models the rigidbody as a Black-Box Neural Network

$(SIGNATURES)

**Arguments**:

- `qdim:: Int`: Dimensions of position configuration
- `udim:: Int`: Dimensions of active control
- `thetamask`: Binary Vector with ones corresponding to angles in the configuration space
"""
function Naive(qdim::Int, udim::Int, thetamask; activation=Tanh(), hidden_sizes=(128, 128, 128))
    nn = f64(NaiveRigidBody(qdim, udim, hidden_sizes, thetamask))
    return nn
end

include("delan.jl")

"""
Models the rigidbody as a Structured Mechanical Model with no constraints on the external forces.

$(SIGNATURES)

**Arguments**:

- `qdim`: Dimensions of position configuration
- `udim`: Dimensions of active control
- `thetamask`: Binary Vector with ones corresponding to angles in the configuration space

Useful when no prior knowledge about the active forces on a mechanical system.
"""
function DeLan(qdim, udim, thetamask; activation=Tanh(), hidden_sizes=(32, 32, 32), bias=10.0)
    nn = f64(DelanRigidBody(qdim, udim, hidden_sizes, thetamask, activation=activation, bias=bias))
    return nn
end

"""
Models the rigidbody as a Structured Mechanical Model with the constraints that active forces on the system are control affine.

$(SIGNATURES)

**Arguments**:

- `qdim`: Dimensions of position configuration
- `udim`: Dimensions of active control
- `thetamask`: Binary Vector with ones corresponding to angles in the configuration space

Almost all mechanical systems are control affine in the input forces for control.
"""
function ControlAffine(qdim, udim, thetamask; activation=Tanh(), hidden_sizes=(32, 32, 32), bias=10.0)
    mm = DelanMassMatrix(qdim, hidden_sizes, activation)
    pot = DelanPotential(qdim, hidden_sizes, activation)
    ff = GeneralizedForces([ControlAffineForce(qdim, udim, hidden_sizes), QVForce(qdim, udim, hidden_sizes)])
    nn = f64(DelanRigidBody(mm, pot, ff, qdim, udim, bias, thetamask))
    return nn
end

include("dataset.jl")

include("datagen.jl")

#include("supervised.jl")

include("policy.jl")

include("pytorch.jl")
using .Torch

module Dynamics
    include(joinpath("dynamics", "furuta.jl"))
    include(joinpath("dynamics", "doublecartpole.jl"))
    include(joinpath("dynamics", "triplecartpole.jl"))
end
using .Dynamics

export Model

end # module


text:
module Walkers

const version = "2.1.0"

using GLMakie, GLFW, FileIO
using LinearAlgebra, Colors
using Random: seed!
import Base./

include("matutils.jl")

const NaNPoint3f0 = Point3f0(NaN32, NaN32, NaN32)

@enum Laws position=0 velocity=1 newtonlinear=2 newton=3 cyclical=4

@enum Relations onetoone=0 sparse=1 manytomany=2 electronic=3

"Display an observable each time it is changed."
monitor = obs::Observable -> on(display, obs)

"Load an image file from disk and return it to be used as icon in GLFW."
function load_icon(path)
    icon = load(path)
    reinterpret(NTuple{4, UInt8}, icon)
end

"Return the values of multiple labeled sliders."
function ls_values(sliders...)
    [s.slider.value for s in sliders]
end

"Return the layouts of multiple labeled sliders."
function ls_layouts(sliders...)
    [s.layout for s in sliders]
end

"Defines broadcast on real by point division."
function /(x::Real, pt::Point{})
    x ./ pt
end

"Maps `x` from the range `[0, 1]` to `[avg-width, avg+width]`."
function scatter(x::Any, avg=0::Real, width=1::Real)
    width * (2x .- 1) .+ avg
end

"Builds a linear random hue colormap of `steps` elements."
function huecolmap(steps=2; s=1.0, l=0.5, a=1.0)
    hues = range(rand(0:359), stop=rand(0:359), length=steps)
    cmap = Array{RGBA}(steps)
    for i = 1:steps
        color = HSLA(hues[i], s, l, a)
        cmap[i] = convert(RGBA, color)
    end
    cmap
end

"Computes walkers states."
function walk(law::Laws, n::Int, pos::Array, rels::Array)
    states = Point3f0[]
    vel = fill(Vec3f0(0, 0, 0), length(pos)) # Null initial velocity
    for i = 1:n
        append!(states, pos)
        if law == position
            # A specific part of the distance between a walker and the others
            # will be added to its position.
            pos += rels * pos - sum(rels .* pos, dims=2)
        elseif law == cyclical
            # Alg that does not uses attraction directly but variances between
            # attraction values as position modulation.
            pos += rels * pos - transpose(transpose(pos) * rels)
        else
            forces = diffs(pos) # Distance between all walkers locations
            if law == newton || law == newtonlinear
                # Newton's and Coulomb's forces norms are of form k/d².
                # newtonlinear is variation of form k/d
                # normalize(f) returns a direction vector of norm 1, while
                # norm.(f).^pow is the norm of the new force vector.
                pow = law == newton ? -2 : -1
                forces = 10^4 * nulldiag(normalize(forces) .* norm.(forces) .^ pow)
            end
            # Forces are modulated by walkers relations, then instant velocity
            # is incremented by the resulting acceleration. Position is then
            # incremented by instant velocity.
            vel += sum(.1rels .* forces, dims=2)
            pos += vel
        end
    end
    states
end

function app()
    set_window_config!(
        title = "Walkers Alpha",
        decorated = true,
        # renderloop = renderloop,
    )

    set_theme!(
        fontsize = 14,
        colgap = 0,
        rowgap = 0
    )

    app_scene = Scene(camera=campixel!)

    app_layout = GridLayout(
        app_scene,
        2, 2,
        colsizes = [Fixed(400), Auto()]
    )

    view_scene = LScene(app_scene, camera=cam3d!, raw=false)
    view_scene.scene[:show_axis] = false

    count_ls    = labelslider!(app_scene, "Walkers count",       2:40;                       sliderkw=Dict(:startvalue=>5))
    spread_ls   = labelslider!(app_scene, "Walkers spread",      LinRange(0f0, 100f0, 101);  sliderkw=Dict(:startvalue=>50f0))
    rel_avg_ls  = labelslider!(app_scene, "Average attraction",  LinRange(-.1f0, .1f0, 101); sliderkw=Dict(:startvalue=>0f0))
    rel_var_ls  = labelslider!(app_scene, "Attraction variance", LinRange(0f0, 1f0, 101);    sliderkw=Dict(:startvalue=>0f0))
    iters_ls    = labelslider!(app_scene, "Iterations",          2:1000;                     sliderkw=Dict(:startvalue=>10))
    rotspeed_ls = labelslider!(app_scene, "Rotation speed",      LinRange(0f0, 1f0, 101);    sliderkw=Dict(:startvalue=>.2f0))

    settings_layout = GridLayout(
        tellheight = false,
        alignmode = Outside(10),
    )

    settings_layout[1:6, 1] = ls_layouts(count_ls, spread_ls, rel_avg_ls, rel_var_ls, iters_ls, rotspeed_ls)

    # Don't render background for noow because it prevents interaction with sliders
    # app_layout[1, 1] = Box(app_scene, color=(:black, .04), strokevisible=false)
    app_layout[1, 1] = settings_layout
    app_layout[1, 2] = view_scene
    app_layout[2, 1:2] = Box(app_scene, color=RGBA(0, 0, 0, .08), strokevisible=false)
    app_layout[2, 1:2] = Label(app_scene, "Welcome to Walkers Alpha!", alignmode=Outside(3))

    #=
        # GUI parameters
        speed_gui,      speed_s      = textslider(0:.1:10,            editscreen)
        # cmap_gui,       cmap_s       = widget(default_cmap,               editscreen)
        # path_color_gui, path_color_s = widget(Signal(default_path_color), editscreen)
        # law_gui,        law_s        = widget(Signal(position),           editscreen)
        # relation_gui,   relation_s   = widget(Signal(onetoone),           editscreen)
        center_gui,     center_s     = button("⛶",                        editscreen)
        regen_gui,      regen_s      = button("↻",                        editscreen)
        params = Pair[
            "Rotation speed" => speed_gui,
            "Color map"      => cmap_gui,
            "Paths color"    => path_color_gui,
            "Relation rel_model" => relation_gui,
            "Dynamics law"   => law_gui,
            "Center"         => center_gui,
            "Regenerate"     => regen_gui
        ]

        # DEBUG
        # speed_s      = Signal(1)
        # cmap_s       = Signal()
        # center_s     = Signal(false)
        # regen_s      = Signal(false)

        # Tray button actions
        preserve(map(center_s, regen_s) do center, regen # preserve() avoids GC
            if center
                center!(viewscreen)
            elseif regen
                seed = abs(rand(Int32)) # Uses the global RNG
                push!(seed_s, seed)
            end
            nothing
        end)

        # Save stuff
        key_pressed = false
        preserve(map(window.inputs[:keyboard_buttons]) do ksam # preserve() avoids GC
        key, scancode, action, mods = ksam
        if key == GLFW.KEY_S && !key_pressed
            name = rand(0000:9999)
            screenshot(viewscreen, path="$name.png",       channel=:color)
            screenshot(viewscreen, path="$name-depth.png", channel=:depth)
            key_pressed = true
        else
            key_pressed = false
        end
        nothing
        end)

        # Time and rotation signal mappings
        time_s = loop(linspace(0f0, Float32(typemax(Int32)), typemax(Int32)))
        rot_s = map(time_s, speed_s) do t, s
            rotationmatrix_z(Float32(t * s/360)) # -> 4x4 Float32 rotation matrix
        end
    =#

    seed = 13121312
    law = position
    rel_model = onetoone
    rings_colors_stops::Array{RGBAf0, 1} = [
        RGBA(.27, .01, .33, .7),
        RGBA(.13, .56, .55, .6),
        RGBA(.99, .91, .14, .4)
    ]
    paths_color = RGBA(.1, .1, .1, .8)

    # RNG
    rng = MersenneTwister(0)

    # Init points
    points = lift(ls_values(count_ls, spread_ls)...) do count, spread
        seed!(rng, seed)
        scatter.(rand(rng, Point3f0, count), 0, spread)
    end

    # Relations matrix mask
    rels_mask = lift(ls_values(count_ls)...) do n
        seed!(rng, seed + 1)
        if rel_model == onetoone
            # One walker is in relation with another
            mask = offsetcols(collect(I(n)))
        elseif rel_model == sparse
            # Sparse is ManyToMany with only 25% of relations.
            # As the diagonal will be nulled, compensation is needed.
            mask = rand(rng, Float32, n, n) .< .25f0 + 1f0/n
        # elseif rel_model == electronic
        #     # Relation matrix as if each walker behaves as a +/- charged particule
        #     loads = repeat(rand(rng, n), 1, n)
        #     rel = -transpose(loads) .* loads
        else
            # Each walker is in relation with all others
            mask = trues(n, n)
        end
        nulldiag(mask)
    end

    # Relations matrix
    relations = lift(ls_values(count_ls, rel_avg_ls, rel_var_ls)..., rels_mask) do n, avg, var, mask
        seed!(rng, seed + 2)
        rels = scatter.(rand(rng, Float32, n, n), avg, var)
        rels .* mask
    end

    # System states
    last_states::Union{Array{Point3f0, 2}, Nothing} = nothing
    states = lift(ls_values(count_ls, iters_ls)..., points, relations) do n, iters, points, relations
        if length(points) == size(relations, 1)
            states = walk(law, iters, points, relations)
            states = reshape(states, n, iters)
            last_states = states
            return states
        else
            return last_states
        end
    end

    # Geometry
    rings_vertex = @lift begin
        rings = vcat($states, $states[1:1, :])
        rings = vcat(rings, fill(NaNPoint3f0, (1, size(rings, 2))))
        reshape(rings, length(rings))
    end
    paths_vertex = @lift begin
        paths = hcat($states, fill(NaNPoint3f0, size($states, 1)))
        reshape(permutedims(paths), length(paths))
    end

    # Colors
    rings_color = @lift to_colormap(rings_colors_stops, length($rings_vertex))

    # Plotting
    rings = lines!(view_scene, rings_vertex, color=rings_color, transparency=true, linewidth=2)
    paths = lines!(view_scene, paths_vertex, color=paths_color, transparency=true, linewidth=3)

    # Controls
    controls = cameracontrols(view_scene.scene)
    controls.rotationspeed[]= .02

    # Run app & set icon
    display(app_scene)
    window = app_scene.current_screens[1].glscreen
    GLFW.SetWindowIcon(window, load_icon(joinpath("rsc", "logo.png")))

    # Autorotate
    angle = lift(rotspeed_ls.slider.value) do speed
        speed / 5f1
    end
    while isopen(view_scene.scene)
        rotate_cam!(view_scene.scene, to_value(angle), 0f0, 0f0)
        sleep(.01)
    end
end

app() # For tesing purposes, should be removed when compiling

end

max_stars_repo_name: LeMinaw/walkers

max_issues_repo_name: LeMinaw/walkers

max_forks_repo_name: LeMinaw/Walkers


text:
using SparseArrays
using LinearAlgebra

const theta_double = [ 2.22045e-16, 2.58096e-8, 1.38635e-5, 0.000339717, 0.00240088,
0.00906566, 0.0238446, 0.0499123, 0.0895776, 0.144183, 0.214236,
0.299616, 0.399778, 0.513915, 0.641084, 0.780287, 0.930533, 1.09086,
1.26038, 1.43825, 1.62372, 1.81608, 2.01471, 2.21905, 2.42858, 2.64285,
2.86145, 3.084, 3.31017, 3.53967, 3.77221, 4.00756, 4.2455, 4.48582,
4.72835, 4.97292, 5.21938, 5.46759, 5.71744, 5.9688, 6.22158, 6.47568,
6.73102, 6.9875, 7.24507, 7.50365, 7.76317, 8.02359, 8.28485, 8.5469,
8.80969, 9.07319, 9.33734, 9.60212, 9.8675, 10.1334, 10.3999, 10.6669,
10.9343, 11.2022, 11.4705, 11.7392, 12.0084, 12.2778, 12.5477, 12.8178,
13.0883, 13.3591, 13.6302, 13.9016, 14.1732, 14.4451, 14.7172, 14.9896,
15.2622, 15.535, 15.808, 16.0812, 16.3546, 16.6281, 16.9019, 17.1758,
17.4498, 17.724, 17.9984, 18.2729, 18.5476, 18.8223, 19.0972, 19.3723,
19.6474, 19.9227, 20.198, 20.4735, 20.7491, 21.0248, 21.3005, 21.5764,
21.8523, 22.1284]

const theta_single = [1.19209e-7, 0.000597886, 0.0112339, 0.0511662, 0.130849, 0.249529,
0.401458, 0.580052, 0.779511, 0.995184, 1.22348, 1.46166, 1.70765,
1.95985, 2.21704, 2.47828, 2.74282, 3.01007, 3.27956, 3.55093,
3.82386, 4.09811, 4.37347, 4.64978, 4.9269, 5.20471, 5.48311, 5.76201,
6.04136, 6.32108, 6.60113, 6.88146, 7.16204, 7.44283, 7.7238, 8.00493,
8.2862, 8.56759, 8.84908, 9.13065, 9.4123, 9.69402, 9.97579, 10.2576,
10.5394, 10.8213, 11.1032, 11.3852, 11.6671, 11.949, 12.231, 12.5129,
12.7949, 13.0769, 13.3588, 13.6407, 13.9227, 14.2046, 14.4865, 14.7684]

const theta_half = [0.00195058, 0.0744366, 0.266455, 0.524205, 0.810269, 1.10823,
1.41082, 1.71493, 2.01903, 2.32247, 2.62492, 2.9263, 3.22657, 3.52578,
3.82398, 4.12123, 4.41759, 4.71314, 5.00792, 5.302, 5.59543, 5.88825,
6.18051, 6.47224, 6.76348, 7.05427, 7.34464, 7.6346, 7.92419, 8.21342,
8.50232, 8.79091, 9.07921, 9.36722, 9.65497, 9.94247, 10.2297, 10.5168,
10.8036, 11.0902, 11.3766, 11.6629, 11.9489, 12.2348, 12.5205, 12.8061,
13.0915, 13.3768, 13.6619, 13.9469, 14.2318, 14.5165, 14.8012, 15.0857,
15.3702, 15.6547, 15.9391, 16.2235, 16.5078, 16.792, 17.0761, 17.3599,
17.6437, 17.9273, 18.2108, 18.4943, 18.7776, 19.0609, 19.3442, 19.6273,
19.9104, 20.1935, 20.4764, 20.7593, 21.0422, 21.325, 21.6077, 21.8904,
22.173, 22.4556, 22.7381, 23.0206, 23.303, 23.5854, 23.8677, 24.15,
24.4322, 24.7144, 24.9966, 25.2787, 25.5608, 25.8428, 26.1248, 26.4068,
26.6887, 26.9706, 27.2524, 27.5342, 27.816, 28.0978]

function select_taylor_degree(A,
                              b_columns,
                              opnorm = LinearAlgebra.opnorm;
                              m_max = 55,
                              p_max = 8,
                              precision = :double,
                              shift = DEFAULT_SHIFT,
                              # bal,
                              force_estm = true,
                              check_positive = false)

    #SELECT_TAYLOR_DEGREE   Select degree of Taylor approximation.
    #   [M,MV,alpha,unA] = SELECT_TAYLOR_DEGREE(A,m_max,p_max) forms a matrix M
    #   for use in determining the truncated Taylor series degree in EXPMV
    #   and EXPMV_TSPAN, based on parameters m_max and p_max.
    #   MV is the number of matrix-vector products with A or A^* computed.

    #   Reference: A. H. Al-Mohy and N. J. Higham, Computing the action of
    #   the matrix exponential, with an application to exponential
    #   integrators. MIMS EPrint 2010.30, The University of Manchester, 2010.

    #   Awad H. Al-Mohy and Nicholas J. Higham, October 26, 2010.

    n = size(A, 2)

    if p_max < 2 || m_max > 60 || m_max + 1 < p_max*(p_max - 1)
        error("Invalid p_max or m_max.")
    end

    theta =
      if precision == :double
          theta_double
      elseif precision == :single
          theta_single #load theta_taylor_single
      elseif precision == :half
          theta_half
      end

    if shift && hasmethod(tr, Tuple{typeof(A)})
        mu = tr(A)/n
        A = diagshift!(A, mu) #A -= mu * I
    end

    if !force_estm
        @assert hasmethod(opnorm, Tuple{typeof(A), typeof(1)}) "opnorm not defined for $typeof(A)"
        normA = opnorm(A,1)
    end

    if !force_estm && normA <= 4*theta[m_max]*p_max*(p_max + 3)/(m_max*b_columns)
        unA = 1
        c = normA
        alpha = c*ones(p_max-1,1)
    else
        unA = 0
        eta = zeros(p_max, 1)
        alpha = zeros(p_max-1, 1)
        for p = 1:p_max
            c = normAm(A, p + 1; check_positive = check_positive)
            c = c^(1 / (p + 1))
            eta[p] = c
        end
        for p = 1:p_max-1
            alpha[p] = max(eta[p],eta[p+1])
        end
    end

    M = zeros(m_max,p_max-1)
    for p = 2:p_max
        for m = p*(p-1)-1 : m_max
            M[m,p-1] = alpha[p-1]/theta[m]
        end
    end
    return (M,alpha,unA)

end


text:
### A Pluto.jl notebook ###
# v0.16.3

using Markdown
using InteractiveUtils

# ╔═╡ 19288130-32bf-11ec-043f-295f73ca2da1
using CSV, DataFrames, NearestNeighbors

# ╔═╡ 51ae3c5b-8e06-4df3-8672-2c2b96bdc5b2
df = CSV.read("marketing_campaign.csv", DataFrame);

# ╔═╡ fcebb5ba-abcb-4f2e-8de4-a3cd4d15028d
dropmissing!(df);

# ╔═╡ 4fe253b8-e4d9-4a1d-afd2-77d11e77262d
df[1:6, :]

# ╔═╡ 2b21fe56-ac7a-48b4-8645-7175f81e07cc
names(df)

# ╔═╡ acde051d-bf58-4230-9ab1-1fd4b9c193ee
df₁ = df[:, [:Year_Birth, :Income]]; df₁[1:6, :]

# ╔═╡ 745910f8-f2c0-43da-aa5a-5b61023aaa28
md"The ```KDTree``` function will calculate the nearest neighbors for each element in this ```DataFrame```. However, it requires us to present the data as a series of 2216 column vectors, each representing a point in our data.

This can be accomplished easily by just rotating the ```DataFrame``` 90° to the left, which gives the new column indices a 1:1 correspondence with the old row indices.

We also need to convert the ```DataFrame``` to a Julia ```Array``` to work with ```KDTree```, and we might as well convert the data types to ```Float32``` while we're at it."

# ╔═╡ f7b96787-13c3-4752-9cce-f6a2918b6afe
df₂ = Array{Float32}(df₁) |> rotl90

# ╔═╡ 56d167a9-26b6-424b-af56-1b2d726bc22a
kdtree = KDTree(df₂; leafsize=10)

# ╔═╡ b049d762-ca44-4815-afbb-46efaba5ee4c
md"Let's just pick a random point in our data and find its nearest neighbors."

# ╔═╡ 4ee18911-dc24-43be-beb6-79a174ff0c1f
point = df₂[:, end]

# ╔═╡ 370d39ad-4b4e-4643-92fd-3b42fda004d9
idxs, dists = knn(kdtree, point, 7, true)

# ╔═╡ 64594f15-87ae-411a-842c-de37c53234cf
idxs

# ╔═╡ 331a1a4e-9768-4f0e-bf6b-14160586abf6
md"Let's see if these guys seem like they might be similar, I dunno.

(Wow, they actually are *super* similar...)"

# ╔═╡ 51f46b06-520e-4f81-b002-d9848b19caf9
df[idxs[1:end], :]

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"

[compat]
CSV = "~0.9.8"
DataFrames = "~1.2.2"
NearestNeighbors = "~0.4.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "29728b5bf89047611c189f412f3325fff993711b"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.9.8"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "09d9eaef9ef719d2cd5d928a191dc95be2ec8059"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.5"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "7fb0eaac190a7a68a56d2407a6beff1142daf844"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.12"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "19cb49649f8c41de7fea32d089d37de917b553da"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.0.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NearestNeighbors]]
deps = ["Distances", "StaticArrays"]
git-tree-sha1 = "16baacfdc8758bc374882566c9187e785e85c2f0"
uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
version = "0.4.9"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "f19e978f81eca5fd7620650d7dbea58f825802ee"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═19288130-32bf-11ec-043f-295f73ca2da1
# ╠═51ae3c5b-8e06-4df3-8672-2c2b96bdc5b2
# ╠═fcebb5ba-abcb-4f2e-8de4-a3cd4d15028d
# ╠═4fe253b8-e4d9-4a1d-afd2-77d11e77262d
# ╠═2b21fe56-ac7a-48b4-8645-7175f81e07cc
# ╠═acde051d-bf58-4230-9ab1-1fd4b9c193ee
# ╟─745910f8-f2c0-43da-aa5a-5b61023aaa28
# ╠═f7b96787-13c3-4752-9cce-f6a2918b6afe
# ╠═56d167a9-26b6-424b-af56-1b2d726bc22a
# ╟─b049d762-ca44-4815-afbb-46efaba5ee4c
# ╠═4ee18911-dc24-43be-beb6-79a174ff0c1f
# ╠═370d39ad-4b4e-4643-92fd-3b42fda004d9
# ╠═64594f15-87ae-411a-842c-de37c53234cf
# ╟─331a1a4e-9768-4f0e-bf6b-14160586abf6
# ╠═51f46b06-520e-4f81-b002-d9848b19caf9
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: benharmonics/ModelZoo

max_issues_repo_name: benharmonics/ModelZoo

max_forks_repo_name: benharmonics/ModelZoo


text:
using Iris
using Test

@testset "Self Energies" begin
    @testset "1D Symmetric" begin
        shape = Interval(-1,1)
        bnd = Boundary(shape,DirichletBC)
        lat = Lattices.Cartesian(.1)
        dom = LatticeDomain(bnd,lat)
        Σ = @test_nowarn Iris.Common.SelfEnergies.SelfEnergy{Iris.Symmetric}(dom,ones(ComplexF64,length(dom.x)+1))
        for i ∈ 1:20, j ∈ 1:20
            if i==1 && j==1
                @test Σ.Σ0[1][i,j] ≈ 100
                @test Σ.Σ1[1][1,1] ≈ 10
            elseif i==20 && j==20
                @test Σ.Σ0[2][20,20] ≈ 100
                @test Σ.Σ1[2][20,20] ≈ -10
            else
                @test Σ.Σ1[1][i,j] ≈ 0
                @test Σ.Σ1[2][i,j] ≈ 0
            end
        end
    end
    @testset "2D Symmetric" begin
        shape = Square(1,0,0)
        bnd = Boundary(shape,MatchedBC{1}(in=[1]),MatchedBC{2}(out=[1]),DirichletBC{3}(),DirichletBC{4}(),PML;depth=.3)
        lat = Lattices.Cartesian(.05,.05)
        dom = LatticeDomain(bnd,lat)
        s = @test_nowarn Iris.Common.SelfEnergies.SelfEnergy{Iris.Symmetric}(dom,ones(ComplexF64,20),ones(ComplexF64,20))
    end
end


text:
# Just a simple set of tests to ensure that the basis of any given quaternionic
# type actually obeys the standard quaternion rules.

@testset verbose=true "Basis" begin
    @testset "$T" for T in Types
        for Q in [Quaternion, Rotor]
            # Define basis elements
            u = Q{T}(1)
            i = Q{T}(𝐢)
            j = Q{T}(𝐣)
            k = Q{T}(𝐤)

            # Standard expressions
            @test u * u ≈ u atol=eps(T)
            @test i * i ≈ -u atol=eps(T)
            @test j * j ≈ -u atol=eps(T)
            @test k * k ≈ -u atol=eps(T)
            @test i * j * k ≈ -u atol=eps(T)

            # Full multiplication table
            @test u * u ≈ u atol=eps(T)
            @test u * i ≈ i atol=eps(T)
            @test u * j ≈ j atol=eps(T)
            @test u * k ≈ k atol=eps(T)
            @test i * u ≈ i atol=eps(T)
            @test i * i ≈ -u atol=eps(T)
            @test i * j ≈ k atol=eps(T)
            @test i * k ≈ -j atol=eps(T)
            @test j * u ≈ j atol=eps(T)
            @test j * i ≈ -k atol=eps(T)
            @test j * j ≈ -u atol=eps(T)
            @test j * k ≈ i atol=eps(T)
            @test k * u ≈ k atol=eps(T)
            @test k * i ≈ j atol=eps(T)
            @test k * j ≈ -i atol=eps(T)
            @test k * k ≈ -u atol=eps(T)
        end

        for Q in [Quaternion, QuatVec]
            # Define basis elements
            u = Q{T}(1)
            i = Q{T}(𝐢)
            j = Q{T}(𝐣)
            k = Q{T}(𝐤)
            basis = [u, i, j, k]

            # Basic self-addition/subtraction
            @test u + u == 2u
            @test u - u == 0u
            @test i + i == 2i
            @test i - i == 0i
            @test j + j == 2j
            @test j - j == 0j
            @test k + k == 2k
            @test k - k == 0k

            # Full addition/subtraction table
            for i1 in 1:4
                for i2 in 1:4
                    a = zeros(T, 4)
                    s = zeros(T, 4)
                    a[i1] += one(T)
                    a[i2] += one(T)
                    s[i1] += one(T)
                    s[i2] -= one(T)
                    @test basis[i1] + basis[i2] == Q{T}(a...)
                    @test basis[i1] - basis[i2] == Q{T}(s...)
                end
            end

            # Normalization
            for q in basis
                n = normalize(q)
                @test typeof(n) === Q{float(T)}
                if Q === QuatVec && q == u
                    continue
                end
                @test q == n
                @test q == normalize(2n)
            end
        end

        let Q = Rotor
            # Define basis elements
            u = Q{T}(1)
            i = Q{T}(𝐢)
            j = Q{T}(𝐣)
            k = Q{T}(𝐤)
            basis = [u, i, j, k]

            # Normalization
            for q in basis
                n = normalize(q)
                @test typeof(n) === Q{float(T)}
                @test q == n
                @test q == normalize(2n)
                twoq = Q{T}(2*(q.components))  # Won't normalize
                @test q == normalize(twoq)
            end
        end

        let Q = QuatVec
            i = Q{T}(𝐢)
            j = Q{T}(𝐣)
            k = Q{T}(𝐤)

            @test i × i == zero(i)
            @test j × j == zero(j)
            @test k × k == zero(k)
            @test i × j == k
            @test j × k == i
            @test k × i == j
            @test j × i == -k
            @test k × j == -i
            @test i × k == -j

            @test i ×̂ i == zero(i)
            @test j ×̂ j == zero(j)
            @test k ×̂ k == zero(k)
            @test i ×̂ j == k
            @test j ×̂ k == i
            @test k ×̂ i == j
            @test j ×̂ i == -k
            @test k ×̂ j == -i
            @test i ×̂ k == -j

            @test i ×̂ 2j == k
            @test j ×̂ 2k == i
            @test k ×̂ 2i == j
            @test j ×̂ 2i == -k
            @test k ×̂ 2j == -i
            @test i ×̂ 2k == -j

            @test 2i ×̂ j == k
            @test 2j ×̂ k == i
            @test 2k ×̂ i == j
            @test 2j ×̂ i == -k
            @test 2k ×̂ j == -i
            @test 2i ×̂ k == -j
        end

    end
end


text:
X = [-1.3046   -0.1606
   -1.4341   -0.3372
   -1.3475   -0.0421
   -1.3426    0.0746
   -1.2433   -0.0451
   -1.3477    0.0140
   -1.2695    0.0159
   -1.2560   -0.0217
   -1.4525   -0.1600
   -1.4813    0.2405
   -1.2533    0.1410
   -1.1160    0.1336
   -1.2506    0.0134
   -1.4179    0.1585
   -1.1427    0.1662
   -1.2654    0.2749
   -1.3093    0.4972
   -0.8894    0.5117
   -1.2147    0.4365
   -1.0589    0.5416
   -1.0974    0.3831
   -1.1807    0.5650
   -1.0185    0.5226
   -1.1234    0.5162
   -0.9260    0.5602
   -1.1131    0.6486
   -1.0704    0.4117
   -1.3090    0.6600
   -0.9298    0.6564
   -0.7473    0.7469
   -0.8765    0.7687
   -1.0327    0.6897
   -0.8800    0.6879
   -0.7897    0.8371
   -0.9968    0.5561
   -0.7329    0.8403
   -0.9754    0.4914
   -0.5870    0.8492
   -0.7836    0.9945
   -0.5976    0.7718
   -0.6082    1.1435
   -0.5593    0.8088
   -0.4512    0.6593
   -0.4470    0.9503
   -0.5153    0.8144
   -0.4882    0.8137
   -0.3907    0.9800
   -0.3310    0.9054
   -0.4595    0.8447
   -0.2568    0.8148
   -0.1615    0.6493
   -0.1068    0.9994
   -0.3376    0.9333
   -0.0439    0.8333
   -0.0312    0.8441
   -0.1402    0.8472
   -0.0961    1.0671
    0.0436    0.9985
   -0.2503    0.7038
   -0.1422    0.9170
   -0.1054    0.9390
   -0.0536    0.8220
    0.3137    0.8682
    0.0720    0.6508
    0.0782    0.8938
   -0.0973    0.6904
    0.1411    0.7367
    0.1296    0.7934
    0.2079    0.5609
    0.2143    0.8637
    0.0133    0.8896
    0.2598    0.6927
    0.1868    0.9497
    0.1563    0.3887
    0.2798    0.5154
    0.2081    0.6446
    0.3619    0.4604
    0.4530    0.2400
    0.3655    0.4805
    0.4878    0.3051
    0.3224    0.4864
    0.4199    0.3399
    0.4536    0.2371
    0.4703    0.3144
    0.3250    0.3584
    0.4466    0.2982
    0.6724    0.3708
    0.5432    0.2119
    0.6769    0.1165
    0.5262    0.2860
    0.8319    0.2159
    0.6931   -0.0227
    0.3347   -0.0805
    0.6928    0.0936
    0.4681    0.0717
    0.6455   -0.0896
    0.6603   -0.0498
    0.6617   -0.1023
    0.5885   -0.0414
    0.5709   -0.1069
   -0.5597    0.1367
   -0.6311    0.2434
   -0.4894    0.0448
   -0.6578    0.1441
   -0.4873    0.2604
   -0.5392   -0.0516
   -0.6547   -0.0502
   -0.6288   -0.1486
   -0.7275   -0.1015
   -0.6298    0.0501
   -0.3436   -0.2539
   -0.6768   -0.1921
   -0.6052   -0.2849
   -0.4088   -0.1798
   -0.5129   -0.2064
   -0.3433   -0.3513
   -0.5111   -0.3768
   -0.4578   -0.3235
   -0.3441   -0.5589
   -0.4199   -0.3592
   -0.3887   -0.4941
   -0.6470   -0.5287
   -0.4376   -0.4311
   -0.2540   -0.5804
   -0.5833   -0.5855
   -0.2180   -0.5837
   -0.3327   -0.4320
   -0.4364   -0.6702
   -0.4299   -0.6884
   -0.4017   -0.7903
   -0.2743   -0.5594
   -0.2300   -0.5676
   -0.0773   -0.6538
   -0.0194   -0.7473
   -0.1527   -0.7179
   -0.2595   -0.5159
   -0.0402   -0.8188
    0.1271   -0.7987
    0.0809   -0.5309
    0.1404   -0.8573
    0.1402   -0.6199
   -0.0565   -0.6858
    0.0977   -0.8960
    0.1566   -0.7853
    0.2867   -0.9838
    0.1906   -1.0191
    0.0526   -0.7277
   -0.0401   -0.8271
    0.3977   -0.7768
    0.5240   -0.9277
    0.3055   -0.9099
    0.5350   -0.8659
    0.5637   -0.8476
    0.4200   -0.9719
    0.5199   -0.7543
    0.6421   -0.7453
    0.6160   -0.7585
    0.4763   -0.8627
    0.5213   -0.8194
    0.6220   -0.7689
    0.6079   -0.7382
    0.8823   -0.8498
    0.6742   -0.6028
    1.0800   -0.5724
    0.7631   -0.7106
    0.8309   -0.7879
    0.9476   -0.7258
    0.6503   -0.6596
    0.9792   -0.6558
    0.9327   -0.5266
    0.9087   -0.4398
    1.2249   -0.5257
    1.0020   -0.6705
    0.8749   -0.6966
    0.9947   -0.4664
    1.2899   -0.7225
    1.2173   -0.3204
    1.1954   -0.5276
    1.2815   -0.4800
    1.2609   -0.4442
    1.2198   -0.5365
    1.2043   -0.3410
    1.2503   -0.3602
    1.2035   -0.3879
    1.3811   -0.0959
    1.2714   -0.1807
    1.4145   -0.0326
    1.1280   -0.3075
    1.3946   -0.0686
    1.1405   -0.2190
    1.2109   -0.1448
    1.4998    0.0787
    1.1156   -0.0481
    1.5878    0.1745
    1.2817    0.4193
    1.3417   -0.0501
    1.2290   -0.1015
    1.4039    0.2343
    1.2912    0.1612
	1.2759    0.1342];

transpose(X);

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
@inline function loopheader!(integrator::SDDEIntegrator)
    # Apply right after iterators / callbacks

    # Accept or reject the step
    if integrator.iter > 0
      if ((integrator.opts.adaptive && integrator.accept_step) || !integrator.opts.adaptive) && !integrator.force_stepfail
        integrator.success_iter += 1
        # StochasticDiffEq.apply_step!(integrator)
        StochasticDiffEq.apply_step!(integrator)
      elseif integrator.opts.adaptive && !integrator.accept_step
        if integrator.isout
          integrator.dtnew = integrator.dt*integrator.opts.qmin
        elseif !integrator.force_stepfail
          step_reject_controller!(integrator,getalg(integrator.alg))
        end
        StochasticDiffEq.choose_algorithm!(integrator,integrator.cache)
        StochasticDiffEq.fix_dtnew_at_bounds!(integrator)
        StochasticDiffEq.modify_dtnew_for_tstops!(integrator)
        reject_step!(integrator,integrator.dtnew)
        integrator.dt = integrator.dtnew
        integrator.sqdt = sqrt(abs(integrator.dt))
      end
    end

    integrator.iter += 1
    integrator.force_stepfail = false
  end

@inline function loopfooter!(integrator::SDDEIntegrator)
  ttmp = integrator.t + integrator.dt
  if integrator.force_stepfail
    if integrator.opts.adaptive
      integrator.dtnew = integrator.dt/integrator.opts.failfactor
    elseif integrator.last_stepfail
      return
    end
    integrator.last_stepfail = true
    integrator.accept_step = false
  elseif integrator.opts.adaptive
    stepsize_controller!(integrator,getalg(integrator.alg))
    integrator.isout = integrator.opts.isoutofdomain(integrator.u,integrator.p,ttmp)
    integrator.accept_step = (!integrator.isout && accept_step_controller(integrator, integrator.opts.controller)) || (integrator.opts.force_dtmin && integrator.dt <= integrator.opts.dtmin)
    if integrator.accept_step # Accepted
      step_accept_controller!(integrator,getalg(integrator.alg))
      integrator.last_stepfail = false
      integrator.tprev = integrator.t
      if typeof(integrator.t)<:AbstractFloat && !isempty(integrator.opts.tstops)
        tstop = integrator.tdir * first(integrator.opts.tstops)
        @fastmath abs(ttmp - tstop) < 10eps(integrator.t) ? (integrator.t = tstop) : (integrator.t = ttmp)
      else
        integrator.t = ttmp
      end
      StochasticDiffEq.calc_dt_propose!(integrator)
      StochasticDiffEq.handle_callbacks!(integrator)
    end
  else # Non adaptive
    integrator.tprev = integrator.t
    if typeof(integrator.t)<:AbstractFloat && !isempty(integrator.opts.tstops)
      tstop = integrator.tdir * first(integrator.opts.tstops)
      # For some reason 10eps(integrator.t) is slow here
      # TODO: Allow higher precision but profile
      @fastmath abs(ttmp - tstop) < 10eps(max(integrator.t,tstop)) ? (integrator.t = tstop) : (integrator.t = ttmp)
    else
      integrator.t = ttmp
    end
    integrator.last_stepfail = false
    integrator.accept_step = true
    integrator.dtpropose = integrator.dt
    StochasticDiffEq.handle_callbacks!(integrator)
  end
  if integrator.opts.progress && integrator.iter%integrator.opts.progress_steps==0
    @logmsg(-1,
    integrator.opts.progress_name,
    _id = :StochasticDelayDiffEq,
    message=integrator.opts.progress_message(integrator.dt,integrator.u,integrator.p,integrator.t),
    progress=integrator.t/integrator.sol.prob.tspan[2])
  end
end

@inline function postamble!(integrator::SDDEIntegrator)
  StochasticDiffEq.solution_endpoint_match_cur_integrator!(integrator)
  resize!(integrator.sol.t,integrator.saveiter)
  resize!(integrator.sol.u,integrator.saveiter)
  if integrator.opts.progress
    @logmsg(-1,
    integrator.opts.progress_name,
    _id = :StochasticDiffEq,
    message=integrator.opts.progress_message(integrator.dt,integrator.u,integrator.p,integrator.t),
    progress="done")
  end
  return nothing
end

function DiffEqBase.auto_dt_reset!(integrator::SDDEIntegrator)
  # @unpack f,g, u, t, tdir, opts, sol, destats = integrator
  @unpack f,g, u, t, tdir, opts, sol = integrator
  @unpack prob = sol
  @unpack abstol, reltol, internalnorm = opts

  # determine maximal time step
  if has_constant_lags(prob)
    dtmax = tdir * min(abs(opts.dtmax), minimum(abs, prob.constant_lags))
  else
    dtmax = opts.dtmax
  end

  # determine initial time step
  sde_prob = SDEProblem(f,g, prob.u0, prob.tspan, prob.p)
  integrator.dt = StochasticDiffEq.sde_determine_initdt(integrator.u, integrator.t,
      integrator.tdir, integrator.opts.dtmax, integrator.opts.abstol, integrator.opts.reltol,
      integrator.opts.internalnorm, sde_prob, StochasticDiffEq.get_current_alg_order(getalg(integrator.alg), integrator.cache), integrator)
  # update statistics
  # destats.nf += 2
  # destats.nf2 += 2

  nothing
end

DiffEqBase.has_reinit(integrator::SDDEIntegrator) = false # TODO true - syncronize with DDE!
function DiffEqBase.reinit!(integrator::SDDEIntegrator,u0 = integrator.sol.prob.u0;
  t0 = integrator.sol.prob.tspan[1], tf = integrator.sol.prob.tspan[2],
  erase_sol = true,
  tstops = integrator.opts.tstops_cache,
  saveat = integrator.opts.saveat_cache,
  d_discontinuities = integrator.opts.d_discontinuities_cache,
  order_discontinuity_t0 = t0 == integrator.sol.prob.tspan[1] && u0 == integrator.sol.prob.u0 ? integrator.order_discontinuity_t0 : 0,
  reinit_cache = true,reinit_callbacks = true,
  initialize_save = true,
  reset_dt = (integrator.dtcache == zero(integrator.dt)) && integrator.opts.adaptive)

  if DiffEqBase.isinplace(integrator.sol.prob)
    recursivecopy!(integrator.u,u0)
    recursivecopy!(integrator.uprev,integrator.u)
  else
    integrator.u = u0
    integrator.uprev = integrator.u
  end

  integrator.t = t0
  integrator.tprev = t0

  tType = typeof(integrator.t);
  maximum_order = StochasticDiffEq.alg_order(getalg(integrator.alg))
  tstops_internal, saveat_internal, d_discontinuities_internal =
    tstop_saveat_disc_handling(tstops, saveat, d_discontinuities, (tType(t0), tType(tf)),
                                              order_discontinuity_t0, maximum_order,
                                              integrator.sol.prob.constant_lags,
                                              integrator.sol.prob.neutral)

  integrator.opts.tstops = tstops_internal
  integrator.opts.saveat = saveat_internal
  integrator.opts.d_discontinuities = d_discontinuities_internal
  integrator.order_discontinuity_t0 = order_discontinuity_t0

  if erase_sol
    if integrator.opts.save_start
      resize_start = 1
    else
      resize_start = 0
    end
    resize!(integrator.sol.u,resize_start)
    resize!(integrator.sol.t,resize_start)
    if integrator.sol.u_analytic !== nothing
      resize!(integrator.sol.u_analytic,0)
    end
    if typeof(getalg(integrator.alg)) <: StochasticDiffEq.StochasticDiffEqCompositeAlgorithm
      resize!(integrator.sol.alg_choice,resize_start)
    end
    integrator.saveiter = resize_start

    if order_discontinuity_t0 ≤ maximum_order
      resize!(integrator.tracked_discontinuities, 1)
      integrator.tracked_discontinuities[1] = Discontinuity(integrator.tdir * integrator.t, Rational{Int}(order_discontinuity_t0))
    else
      resize!(integrator.tracked_discontinuities, 0)
    end
  end

  integrator.iter = 0
  integrator.success_iter = 0

  # full re-initialize the PI in timestepping
  integrator.qold = integrator.opts.qoldinit
  integrator.q11 = typeof(integrator.t)(1)
  integrator.u_modified = false

  if reset_dt
    DiffEqBase.auto_dt_reset!(integrator)
  end

  if reinit_callbacks
    StochasticDiffEq.initialize_callbacks!(integrator, initialize_save)
  end

  if reinit_cache
    StochasticDiffEq.initialize!(integrator.integrator, integrator.cache)
    # StochasticDiffEq.initialize!(integrator,integrator.cache)
  end

  reinit!(integrator.W,integrator.dt)
  nothing
end

@inline function DiffEqBase.get_du(integrator::SDDEIntegrator)
  (integrator.u - integrator.uprev) / integrator.dt
end

@inline function DiffEqBase.get_du!(out,integrator::SDDEIntegrator)
  DiffEqBase.@.. out = (integrator.u - integrator.uprev) / integrator.dt
end

@inline function add_tstop!(integrator::SDDEIntegrator,t)
  t < integrator.t && error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
  push!(integrator.opts.tstops, integrator.tdir * t)
end

DiffEqBase.change_t_via_interpolation!(integrator::SDDEIntegrator,t,modify_save_endpoint::Type{Val{T}}=Val{false}) where T = StochasticDiffEq.change_t_via_interpolation!(integrator, t, modify_save_endpoint)

# update integrator when u is modified by callbacks
function StochasticDiffEq.handle_callback_modifiers!(integrator::SDDEIntegrator)
    # copied from StochasticDiffEq
    if integrator.P !== nothing && integrator.opts.adaptive
        if integrator.cache isa StochasticDiffEqMutableCache
            oldrate = integrator.P.cache.currate
            integrator.P.cache.rate(oldrate, integrator.u, integrator.p, integrator.t)
        else
            integrator.P.cache.currate = integrator.P.cache.rate(
                integrator.u, integrator.p, integrator.t
            )
        end
    end

    # update heap of discontinuities
    # discontinuity is assumed to be of order 0, i.e. solution x is discontinuous
    push!(integrator.opts.d_discontinuities, Discontinuity(integrator.tdir * integrator.t, 0//1))
end

function DiffEqBase.u_modified!(integrator::SDDEIntegrator, bool::Bool)
    integrator.u_modified = bool
    nothing
end

get_proposed_dt(integrator::SDDEIntegrator) = integrator.dtpropose
set_proposed_dt!(integrator::SDDEIntegrator,dt::Number) = (integrator.dtpropose = dt)

function set_proposed_dt!(integrator::SDDEIntegrator,integrator2::SDDEIntegrator)
  integrator.dtpropose = integrator2.dtpropose
  integrator.qold = integrator2.qold
end

@inline function postamble!(integrator::HistorySDEIntegrator)
  if integrator.saveiter == 0 || integrator.sol.t[integrator.saveiter] != integrator.t
    integrator.saveiter += 1
    copyat_or_push!(integrator.sol.t, integrator.saveiter, integrator.t)
    if integrator.opts.save_idxs === nothing
      copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
    else
      copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u[integrator.opts.save_idxs],Val{false})
    end
  end

  resize!(integrator.sol.t,integrator.saveiter)
  resize!(integrator.sol.u,integrator.saveiter)
end

function DiffEqBase.postamble!(integrator::SDDEIntegrator)
  # clean up solution of the SDE integrator
  postamble!(integrator.integrator)

  # # clean solution of the SDDE integrator
  # StochasticDiffEq._postamble!(integrator)
end

@inline function DiffEqBase.savevalues!(integrator::SDDEIntegrator,force_save=false)::Tuple{Bool,Bool}
  saved, savedexactly = false, false
  !integrator.opts.save_on && return saved, savedexactly
  tdir_t = integrator.tdir * integrator.t
  while !isempty(integrator.opts.saveat) && first(integrator.opts.saveat) <= tdir_t # Perform saveat
    integrator.saveiter += 1; saved = true
    curt = integrator.tdir * pop!(integrator.opts.saveat)
    if curt!=integrator.t # If <t, interpolate
      Θ = (curt - integrator.tprev)/integrator.dt
      val = StochasticDiffEq.sde_interpolant(Θ,integrator,integrator.opts.save_idxs,Val{0}) # out of place, but force copy later
      save_val = val
      copyat_or_push!(integrator.sol.t,integrator.saveiter,curt)
      copyat_or_push!(integrator.sol.u,integrator.saveiter,save_val,Val{false})
      if typeof(integrator.alg) <: StochasticDiffEq.StochasticDiffEqCompositeAlgorithm
        copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
      end
    else # ==t, just save
      savedexactly = true
      copyat_or_push!(integrator.sol.t,integrator.saveiter,integrator.t)
      if integrator.opts.save_idxs === nothing
        copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
      else
        copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u[integrator.opts.save_idxs],Val{false})
      end
      if typeof(integrator.alg) <: Union{StochasticDiffEq.StochasticDiffEqCompositeAlgorithm,StochasticDiffEq.StochasticDiffEqRODECompositeAlgorithm}
        copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
      end
    end
  end
  if force_save || integrator.opts.save_everystep
    integrator.saveiter += 1; saved, savedexactly = true, true
    if integrator.opts.save_idxs === nothing
      copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u)
    else
      copyat_or_push!(integrator.sol.u,integrator.saveiter,integrator.u[integrator.opts.save_idxs],Val{false})
    end
    copyat_or_push!(integrator.sol.t,integrator.saveiter,integrator.t)
    if typeof(integrator.alg) <: Union{StochasticDiffEq.StochasticDiffEqCompositeAlgorithm,StochasticDiffEq.StochasticDiffEqRODECompositeAlgorithm}
      copyat_or_push!(integrator.sol.alg_choice,integrator.saveiter,integrator.cache.current)
    end
  end
  return saved, savedexactly
end

@inline function DiffEqNoiseProcess.setup_next_step!(integrator::SDDEIntegrator)
  !isnothing(integrator.W) && DiffEqNoiseProcess.setup_next_step!(integrator.W,integrator.u,integrator.p)
  !isnothing(integrator.P) && DiffEqNoiseProcess.setup_next_step!(integrator.P,integrator.u,integrator.p)
end

@inline function DiffEqNoiseProcess.reject_step!(integrator::SDDEIntegrator,dtnew = integrator.dtnew)
  !isnothing(integrator.W) && reject_step!(integrator.W,dtnew,integrator.u,integrator.p)
  !isnothing(integrator.P) && reject_step!(integrator.P,dtnew,integrator.u,integrator.p)
end

@inline function DiffEqNoiseProcess.accept_step!(integrator::SDDEIntegrator,setup)
  !isnothing(integrator.W) && accept_step!(integrator.W,integrator.dt,integrator.u,integrator.p,setup)
  !isnothing(integrator.P) && accept_step!(integrator.P,integrator.dt,integrator.u,integrator.p,setup)
end

@inline function DiffEqNoiseProcess.save_noise!(integrator::SDDEIntegrator)
  !isnothing(integrator.W) && DiffEqNoiseProcess.save_noise!(integrator.W)
  !isnothing(integrator.P) && DiffEqNoiseProcess.save_noise!(integrator.P)
end


text:
function reversestring(string)
    gnirts = ""
    for i in 1:length(string)
        gnirts = gnirts * string[length(string)-i+1]
    end
    return gnirts
end
print(reversestring("Hello, World!"))

max_stars_repo_name: ellygaytor/Projects

max_issues_repo_name: ellygaytor/Projects

max_forks_repo_name: ellygaytor/Projects


text:
"""
    struct GeneralSpace{𝕜} <: ElementarySpace{𝕜}

A finite-dimensional space over an arbitrary field `𝕜` without additional structure. It is
thus characterized by its dimension, and whether or not it is the dual and/or conjugate
space. For a real field `𝕜`, the space and its conjugate are the same.
"""
struct GeneralSpace{𝕜} <: ElementarySpace{𝕜}
    d::Int
    dual::Bool
    conj::Bool
    function GeneralSpace{𝕜}(d::Int, dual::Bool = false, conj::Bool = false) where {𝕜}
        d >= 0 ||
            throw(ArgumentError("Dimension of a vector space should be bigger than zero"))
        if 𝕜 isa Field
            new{𝕜}(d, dual, (𝕜 ⊆ ℝ) ? false : conj)
        else
            throw(ArgumentError("Unrecognised scalar field: $𝕜"))
        end
    end
end

dim(V::GeneralSpace) = V.d
Base.axes(V::GeneralSpace) = Base.OneTo(dim(V))

dual(V::GeneralSpace{𝕜}) where {𝕜} = GeneralSpace{𝕜}(V.d, !V.dual, V.conj)
isdual(V::GeneralSpace) = V.dual
Base.conj(V::GeneralSpace{𝕜}) where {𝕜} = GeneralSpace{𝕜}(V.d, V.dual, !V.conj)

function Base.show(io::IO, V::GeneralSpace{𝕜}) where {𝕜}
    if V.conj
        print(io, "conj(")
    end
    print(io, "GeneralSpace{", 𝕜, "}(", V.d, ")")
    if V.dual
        print(io, "'")
    end
    if V.conj
        print(io, ")")
    end
end
Base.show(io::IO, ::Type{GeneralSpace}) = print(io, "GeneralSpace")


text:
mutable struct ShaderModule
    mVkDevice::vk.VkDevice
    mHandle::vk.VkShaderModule
    mStage::vk.VkShaderStageFlagBits

    function ShaderModule(device::vk.VkDevice,
                            code::Ptr{UInt8}, codeSize::Int64,
                            stage::vk.VkShaderStageFlagBits = vk.VK_SHADER_STAGE_ALL)
        this = new(device, VkExt.createShaderModule(device, code, codeSize), stage)
    end
end

@class ShaderModule

function destroy(this::ShaderModule)
    vk.vkDestroyShaderModule(this.mVkDevice, this.mHandle, C_NULL)
    this.mHandle = C_NULL
end

function handle(this::ShaderModule)::vk.VkShaderModule
    return this.mHandle
end

function getStage(this::ShaderModule)::vk.VkShaderStageFlagBits
    return this.mStage
end


text:
module InplaceDHMC

export
    # kinetic energy
    GaussianKineticEnergy,
    # NUTS
    TreeOptionsNUTS,
    # reporting
    NoProgressReport, LogProgressReport,
    # mcmc
    TuningNUTS, mcmc_with_warmup, threaded_mcmc, default_warmup_stages

using DocStringExtensions: FIELDS, FUNCTIONNAME, SIGNATURES, TYPEDEF
using LinearAlgebra
using LinearAlgebra: Diagonal, Symmetric
# using LogDensityProblems: capabilities, LogDensityOrder, dimension, logdensity_and_gradient
# import NLSolversBase, Optim # optimization step in mcmc
using UnPack: @unpack
using Random
using Statistics: cov, mean, median, middle, quantile, var

using VectorizedRNG: AbstractPCG, PtrPCG
using VectorizationBase, LoopVectorization, VectorizedRNG, StackPointers, PaddedMatrices, QuasiNewtonMethods, SIMDPirates, Mmap
using QuasiNewtonMethods: dimension, logdensity, ∂logdensity!

# copy from StatsFuns.jl
function logaddexp(x, y)
    isfinite(x) && isfinite(y) || return x > y ? x : y # QuasiNewtonMethods.nanmax(x,y)
    x > y ? x + log1p(exp(y - x)) : y + log1p(exp(x - y))
end

const MMAP = Ref{Matrix{UInt8}}()
const STACK_POINTER_REF = Ref{StackPointer}()
const LOCAL_STACK_SIZE = Ref{Int}()
const GLOBAL_PCGs = Vector{PtrPCG{4}}(undef,0)
const NTHREADS = Ref{Int}()

include("hamiltonian.jl")
include("tree.jl")
include("kinetic_energy.jl")
include("stepsize.jl")
include("NUTS.jl")
include("reporting.jl")
include("warmup.jl")
include("mcmc.jl")
include("rng.jl")
include("diagnostics.jl")
include("precompile.jl")
_precompile_()

function __init__()
    NTHREADS[] = Threads.nthreads()
    # Note that 1 GiB == 2^30 == 1 << 30 bytesy
    # Allocates 0.5 GiB per thread for the stack by default.
    # Can be controlled via the environmental variable PROBABILITY_MODELS_STACK_SIZE
    LOCAL_STACK_SIZE[] = if "PROBABILITY_MODELS_STACK_SIZE" ∈ keys(ENV)
        parse(Int, ENV["PROBABILITY_MODELS_STACK_SIZE"])
    else
        1 << 30
    end + VectorizationBase.REGISTER_SIZE - 1 # so we have at least the indicated stack size after REGISTER_SIZE-alignment
    MMAP[] = Mmap.mmap(Matrix{UInt8}, LOCAL_STACK_SIZE[], NTHREADS[])
    STACK_POINTER_REF[] = StackPointer(
        VectorizationBase.align(Base.unsafe_convert(Ptr{Cvoid}, pointer(MMAP[])))
    )
    STACK_POINTER_REF[] = threadrandinit!(STACK_POINTER_REF[], GLOBAL_PCGs)
    _precompile_()
end
function realloc_stack(new_local_stack_size::Integer)
    @warn """You must redefine all probability models. The stack pointers get dereferenced at compile time, and the stack has just been reallocated.
Re-evaluating densities without first recompiling them will likely crash Julia!"""
    LOCAL_STACK_SIZE[] = new_local_stack_size
    MMAP[] = Mmap.mmap(Matrix{UInt8}, LOCAL_STACK_SIZE[], NTHREADS[])
    STACK_POINTER_REF[] = PaddedMatrices.StackPointer(
        VectorizationBase.align(Base.unsafe_convert(Ptr{Cvoid}, pointer(MMAP[])))
    )
    STACK_POINTER_REF[] = threadrandinit!(STACK_POINTER_REF[], GLOBAL_PCGs)
end
stackpointer() = STACK_POINTER_REF[] + (Threads.threadid() - 1) * LOCAL_STACK_SIZE[]

end # module


text:
# Tests for DirichletCategorical distribution

using  Distributions
using Test, Random, LinearAlgebra

Random.seed!(34567)

rng = MersenneTwister(123)

@testset "Testing DirichletCategorical with $key" for (key, func) in
    Dict("rand(...)" => [rand, rand] ) #,
         #"rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)])

# Constructors
d1 = DirichletCategorical([0.1, 0.2, 0.3], [1, 2, 3])
d2 = DirichletCategorical([4, 2, 3], [1, 1, 0])
d3 = DirichletCategorical(5)
d4 = DirichletCategorical([1.3, 4.3, 2.6, 3.9])

@testset "Constructors" for d in (d1, d2, d3, d4)
    @test typeof(d) == DirichletCategorical{Float64,Int64}
end

@testset "Constructor throws" begin
    @test_throws ArgumentError DirichletCategorical([-1, 0], [0, 0])
    @test_throws ArgumentError DirichletCategorical([0, 0], [-1, 0])
    @test_throws ArgumentError DirichletCategorical([0], [0, 0])
    @test_throws ArgumentError DirichletCategorical([0, 0], [1, 1])
    @test_throws ArgumentError DirichletCategorical([-1, 2], [0, 0])
    @test_throws ArgumentError DirichletCategorical([1.0, 1.0], [2, -1])
end

@testset "ncategories" for (d,n) in [(d1, 3), (d2, 3), (d3, 5), (d4, 4)]
    @test ncategories(d) == n
end

@testset "length" for (d,n) in [(d1, 3), (d2, 3), (d3, 5), (d4, 4)]
    @test length(d) == n
end

@testset "params" begin
    @test params(d1) == ([0.1, 0.2, 0.3], [1, 2, 3])
    @test params(d2) == ([4, 2, 3], [1, 1, 0])
    @test params(d3) == (ones(5), zeros(Int, 5))
    @test params(d4) == ([1.3, 4.3, 2.6, 3.9], zeros(Int, 4))
end

@testset "partype" for d in (d1, d2, d3, d4)
    @test partype(d) == Float64
end

@testset "rand" begin
    @testset "single" for (d,n) in [(d1, 1), (d2, 1), (d3, 1), (d4, 1)]
        @test rand(d) == n
    end

    @testset "multiple" for (d,n) in [
            (d1, [1, 2, 1, 1, 2, 2, 3, 3, 3, 1]),
            (d2, [1, 1, 1, 1, 1, 1, 3, 2, 3, 1]),
            (d3, [1, 1, 1, 1, 3, 2, 5, 3, 5, 1]),
            (d4, [1, 2, 2, 1, 2, 2, 4, 3, 4, 1])]
        @test rand(d, 10) == n
    end

    @testset "10k" for (d,n) in [
            (d1, [1613, 3334, 5053]),
            (d2, [4504, 2720, 2776]),
            (d3, [1932, 2024, 1978, 2055, 2011]),
            (d4, [1058, 3526, 2147, 3269])]
        @test begin
            r = rand(d, 10000)
            counts = map(x -> sum(r .== x), 1:maximum(d))
            counts == n
        end
    end

end

@testset "pdf" for (d,x,p) in [
        (d1, 1, 0.16666666666666669),
        (d2, 2, 0.2727272727272727),
        (d3, 3, 0.2),
        (d4, 4, 0.32231404958677684),
        (d4, -1, 0),
        (d4, 7, 0)]
    @test pdf(d, x) == p
end

@testset "logpdf" for (d,x,p) in [
        (d1, 1, -1.7917594692280547),
        (d2, 2, -1.2992829841302609),
        (d3, 3, -1.6094379124341003),
        (d4, 4, -1.1322288994670948),
        (d4, -1, -Inf),
        (d4, 9, -Inf)]
    @test logpdf(d, x) == p
end

@testset "cpdf" for (d,x,p) in [
        (d1, 1, 0.16666666666666669),
        (d2, 2, 0.7272727272727273),
        (d3, 3, 0.6000000000000001),
        (d4, 4, 1.0),
        (d4, -1, 0),
        (d3, 5, 1.0)]
    @test cdf(d, x) == p
end

@testset "quantile" for (d,x,p) in [(d1, 1, 0.1), (d2, 1, 0.4), (d3, 3, 0.6), (d4, 4, 0.8)]
     @test quantile(d, p) == x

end

@testset "quantile throws" begin
     @test_throws DomainError quantile(d1, -1.0)
     @test_throws DomainError quantile(d1, 2.0)
end

@testset "minimum" for d in (d1, d2, d3, d4)
    @test minimum(d) == 1
end

@testset "maximum" for (d,n) in [(d1, 3), (d2, 3), (d3, 5), (d4, 4)]
    @test maximum(d) == n
end

@testset "insupport" for (d,s) in [
        (d1, [0, 1, 1, 1, 0, 0, 0]),
        (d2, [0, 1, 1, 1, 0, 0, 0]),
        (d3, [0, 1, 1, 1, 1, 1, 0]),
        (d4, [0, 1, 1, 1, 1, 0, 0])]
    @test map(x -> insupport(d, x), collect(0:6)) == s
end

@testset "mode" for (d,n) in [(d1, 3), (d2, 1), (d4, 2)]
    @test mode(d) == n
end

@testset "mode throws" begin
    @test_throws ErrorException mode(d3)
end

@testset "update!" begin
    @testset "single" for (d,o,n) in [
            (d1, 1, [2, 2, 3]),
            (d2, 2, [1, 2, 0]),
            (d3, 3, [0, 0, 1, 0, 0]),
            (d4, 4, [0, 0, 0, 1])]
        @test begin
            d = deepcopy(d)
            Distributions.update!(d, o)

            d.n == n
        end
    end

    @testset "multiple" for (d,o,n) in [
            (d1, [1, 1, 1, 3, 3, 3, 2, 1, 3, 3], [5, 3, 8]),
            (d2, [3, 2, 3, 2, 3, 2, 3, 3, 3, 3], [1, 4, 7]),
            (d3, [3, 4, 5, 4, 4, 4, 5, 5, 4, 4], [0, 0, 1, 6, 3]),
            (d4, [4, 4, 1, 4, 4, 4, 2, 1, 3, 2], [2, 2, 1, 5])]
        @test begin
            d = deepcopy(d)
            Distributions.update!(d, o)

            d.n == n
        end
    end

    @testset "throws" begin
        @test_throws ArgumentError Distributions.update!(d1, -1)
        @test_throws ArgumentError Distributions.update!(d1, 5)
    end

    @testset "rand" begin
        d = DirichletCategorical([10, 0.0001, 0.0001, 0.0001])
        a = rand(d, 100)
        Distributions.update!(d, ones(Int, 10000) .* 4)
        b = rand(d, 100)

        all(a .== 1) && all(b .== 4)
    end

end

end