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
const lorenz63 = lorenz
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} &= a [y - h(x)]\\\\
\\dot{y} &= x - y+z \\\\
\\dot{z} &= b 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``, ``m_0``, and ``m_1``. Setting ``a = 15.6``, ``m_0 = -8/7``
and ``m_1 = -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 (`L1` and `L2`) and mass of each ball (`M1` and `M2`).
Everything is assumed in SI units.

The variables order is ``[θ₁, ω₁, θ₂, ω₂]`` and they satisfy:

```math
\\begin{aligned}
θ̇₁ &= ω₁ \\\\
ω̇₁ &= [M₂ L₁ ω₁² \\sin φ \\cos φ + M₂ G \\sin θ₂ \\cos φ +
       M₂ L₂ ω₂² \\sin φ - (M₁ + M₂) G \\sin θ₁] / (L₁ Δ) \\\\
θ̇₂ &= ω₂ \\\\
ω̇₂ &= [-M₂ L₂ ω₂² \\sin φ \\cos φ + (M₁ + M₂) G \\sin θ₁ \\cos φ -
         (M₁ + M₂) L₁ ω₁² \\sin φ - (M₁ + M₂) G \\sin Θ₂] / (L₂ Δ)
\\end{aligned}
```
where ``φ = θ₂-θ₁`` and ``Δ = (M₁ + M₂) - M₂ \\cos² φ``.

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

(please contribute the Jacobian 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(u, p, t)
    G, L1, L2, M1, M2 = p

    du1 = u[2]

    φ = u[3] - u[1]
    Δ = (M1 + M2) - M2*cos(φ)*cos(φ)

    du2 = (M2*L1*u[2]*u[2]*sin(φ)*cos(φ) +
               M2*G*sin(u[3])*cos(φ) +
               M2*L2*u[4]*u[4]*sin(φ) -
               (M1 + M2)*G*sin(u[1]))/(L1*Δ)

    du3 = u[4]

    du4 = (-M2*L2*u[4]*u[4]*sin(φ)*cos(φ) +
               (M1 + M2)*G*sin(u[1])*cos(φ) -
               (M1 + M2)*L1*u[2]*u[2]*sin(φ) -
               (M1 + M2)*G*sin(u[3]))/(L2*Δ)

    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 \\dot{x} + β x + x^3 = f \\cos(\\omega t)
```
with `f, ω` the forcing strength and frequency and `d` the damping.

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_y \\\\
\\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. That is, 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 ``d_0``.
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 \\\\
\\dot{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 - \\textrm{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| T \\\\
 \\dot{S} &= \\eta_2 - \\eta_3S - |T-S| S
\\end{aligned}
```
Here ``T, S`` denote the dimensionless temperature and salinity differences respectively
between the boxes (polar and equatorial ocean basins) and ``\\eta_i`` are parameters.

[^Stommel1961]: Stommel, Thermohaline convection with two stable regimes of flow. 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

"""
    lorenz84(u = [0.1, 0.1, 0.1]; F=6.846, G=1.287, a=0.25, b=4.0)
Lorenz-84's low order atmospheric general circulation model
```math
\\begin{aligned}
\\dot x = − y^2 − z^2 − ax + aF, \\\\
\\dot y = xy − y − bxz + G, \\\\
\\dot z = bxy + xz − z. \\\\
\\end{aligned}
```

This system has interesting multistability property in the phase space. For the default
parameter set we have four coexisting attractors that gives birth to interesting fractalized
phase space as shown in [^Freire2008]. One can see this by doing:

```
ds = Systems.lorenz84(rand(3))
xg = yg = range(-1.0, 2.0; length=300)
zg = range(-1.5, 1.5; length=30)
bsn, att = basins_of_attraction((xg, yg, zg), ds; mx_chk_att=4)
```

[^Freire2008]: J. G. Freire *et al*,  Multistability, phase diagrams, and intransitivity in the Lorenz-84 low-order atmospheric circulation model, Chaos 18, 033121 (2008)
"""
function lorenz84(u = [0.1, 0.1, 0.1]; F=6.846, G=1.287, a=0.25, b=4.0)
    return ContinuousDynamicalSystem(lorenz84_rule, u, [F, G, a, b],
    lorenz84_rule_jacob)
end
@inline @inbounds function lorenz84_rule(u, p, t)
    F, G, a, b = p
    x, y, z = u
    dx = -y^2 -z^2 -a*x + a*F
    dy = x*y - y - b*x*z + G
    dz = b*x*y + x*z - z
    return SVector{3}(dx, dy, dz)
end
function lorenz84_rule_jacob(u, p, t)
    F, G, a, b = p
	x, y, z = u
    return @SMatrix [(-a)  (-2y)  (-2z);
                     y-b*z  x-1  (-b*x);
                     b*y+z  b*x   x-1]
end

"""
    lorenzdl(u = [0.1, 0.1, 0.1]; R=4.7)
Diffusionless Lorenz system: it is *probably* the simplest rotationnaly invariant
chaotic flow.
```math
\\begin{aligned}
\\dot x = y − x, \\\\
\\dot y = -xz, \\\\
\\dot z = xy - R. \\\\
\\end{aligned}
```

For `R=4.7` this system has two coexisting Malasoma strange attractors that are
linked together as shown in [^Sprott2014]. The fractal boundary between the basins of attractor can be
visualized with a Poincaré section at `z=0`:
```julia
ds = Systems.lorenzdl()
xg = yg = range(-10.0, 10.0; length=300)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2)
bsn, att = basins_of_attraction((xg, yg), pmap)
```

[^Sprott2014]: J. C. Sprott,  Simplest Chaotic Flows with Involutional Symmetries, Int. Jour. Bifurcation and Chaos 24, 1450009 (2014)
"""
function lorenzdl(u = [0.1, 0.1, 0.1]; R=4.7)
    return ContinuousDynamicalSystem(lorenzdl_rule, u, R,
    lorenzdl_rule_jacob)
end
@inline @inbounds function lorenzdl_rule(u, p, t)
    R = p
    x, y, z = u
    dx = y - x
    dy = - x*z
    dz = x*y - R
    return SVector{3}(dx, dy, dz)
end
function lorenzdl_rule_jacob(u, p, t)
    x, y, z = u
    return @SMatrix [-1     1     0;
                     -z     0    -x;
                      y     x     0]
end

"""
    coupled_roessler(u0=[1, -2, 0, 0.11, 0.2, 0.1];
    ω1 = 0.18, ω2 = 0.22, a = 0.2, b = 0.2, c = 5.7, k1 = 0.115, k2 = 0.0)

Two coupled Rössler oscillators, used frequently in the study of chaotic synchronization.
The parameter container has the parameters in the same order as stated in this
function's documentation string.
The equations are:
```math
\\begin{aligned}
\\dot{x_1} &= -\\omega_1 y_1-z_1 \\\\
\\dot{y_1} &= \\omega_1 x+ay_1 + k_1(y_2 - y_1) \\\\
\\dot{z_1} &= b + z_1(x_1-c) \\\\
\\dot{x_2} &= -\\omega_2 y_2-z_2 \\\\
\\dot{y_2} &= \\omega_2 x+ay_2 + k_2(y_1 - y_2) \\\\
\\dot{z_2} &= b + z_2(x_2-c) \\\\
\\end{aligned}
```
"""
function coupled_roessler(u0=[1, -2, 0, 0.11, 0.2, 0.1];
    ω1 = 0.18, ω2 = 0.22, a = 0.2, b = 0.2, c = 5.7, k1 = 0.115, k2 = 0.0)
    p = [ω1, ω2, a, b, c, k1, k2]
    return ContinuousDynamicalSystem(coupled_roessler_f, u0, p)
end
function coupled_roessler_f(u,p,t)
    ω1, ω2, a, b, c, k1, k2 = p

    du1 = -ω1*u[2] - u[3]
    du2 = ω1*u[1] + a*u[2] + k1*(u[5]-u[2])
    du3 = b + u[3]*(u[1]-c)
    du4 = -ω2*u[5] - u[6]
    du5 = ω2*u[4] + a*u[5] + k2*(u[2]-u[5])
    du6 = b + u[6]*(u[4]-c)
    return SVector(du1,du2,du3,du4,du5,du6)
end

"""
    kuramoto(D = 20, u0 = range(0, 2π; length = D);
        K = 0.3, ω = range(-1, 1; length = D)
    )
The Kuramoto model[^Kuramoto1975] of `D` coupled oscillators with equation
```@math
\\dot{\\phi}_i = \\omega_i + \\frac{K}{D}\\sum_{j=1}^{D} \\sin(\\phi_j - \\phi_i)
```

[^Kuramoto1975]: Kuramoto, Yoshiki. International Symposium on Mathematical Problems in Theoretical Physics. 39.
"""
function kuramoto(D = 25, u0 = range(0, 2π; length = D);
    K = 0.3, ω = range(-1, 1; length = D))
    p = KuramotoParams(K, ω)
    return ContinuousDynamicalSystem(kuramoto_f, u0, p)
end
function kuramoto_f(du, u, p, t)
    ω = p.ω; K = p.K
    D = length(u)
    k = K/D
    @inbounds for i in 1:D
        du[i] = ω[i] + k*sum(sin(u[j] - u[i]) for j in 1:D)
    end
    return
end
mutable struct KuramotoParams{T<:Real, V<:AbstractVector{T}}
    K::T
    ω::V
end

"""
    sprott_dissipative_conservative(u0 = [1.0, 0, 0]; a = 2, b = 1, c = 1)
An interesting system due to Sprott[^Sprott2014b] where some initial conditios
such as `[1.0, 0, 0]` lead to quasi periodic motion on a 2-torus, while for
`[2.0, 0, 0]` motion happens on a (dissipative) chaotic attractor.

The equations are:
```math
\\begin{aligned}
\\dot{x} &= y + axy + xz \\\\
\\dot{y} &= 1 - 2x^2 + byz \\\\
\\dot{z_1} &= cx - x^2 - y^2
\\end{aligned}
```
In the original paper there were no parameters, which are added here for exploration purposes.

[^Sprott2014b]: J. C. Sprott. Physics Letters A, 378
"""
function sprott_dissipative_conservative(u0 = [1.0, 0, 0]; a = 2, b = 1, c = 1)
    return CDS(
        sprott_dissipative_conservative_f, u0, [a, b, c], sprott_dissipative_conservative_J
    )
end

function sprott_dissipative_conservative_f(u, p, t)
    a, b, c = p
    x, y, z = u
    dx = y + a*x*y + x*z
    dy = 1 - 2*x^2 + b*y*z
    dz = c*x - x^2 - y^2
    return SVector(dx, dy, dz)
end
function sprott_dissipative_conservative_J(u, p, t)
    a, b, c = p
    x, y, z = u
    return @SMatrix [a*y + z     1 + a*x     +x;
                    -4x     b*z    b*y;
                    (c - 2x)    (-2y)  0]
end

"""
```julia
hodgkinhuxley(u0=[-60.0, 0.0, 0.0, 0.0]; I = 12.0, Vna = 50.0, Vk = -77.0, Vl = -54.4, gna = 120.0,gk = 36.0, gl = 0.3) -> ds
```
```math
\\begin{aligned}
C_m \\frac{dV_m}{dt} = -\\overline{g}_\\mathrm{K} n^4 (V_m - V_\\mathrm{K}) - \\overline{g}_\\mathrm{Na} m^3 h(V_m - V_\\mathrm{Na}) - \\overline{g}_l (V_m - Vl) + I\\\\
\\dot{n} &= \\alpha_n(V_m)(1-n) - \\beta_n(V_m)n \\\\
\\dot{m} &= \\alpha_m(V_m)(1-m) - \\beta_m(V_m)m \\\\
\\dot{h} &= \\alpha_h(V_m)(1-h) - \\beta_h(V_m)h \\\\
\\alpha_n(V_m) = \\frac{0.01(V+55)}{1 - \\exp(\\frac{1V+55}{10})} \\quad
\\alpha_m(V_m) = \\frac{0.1(V+40)}{1 - \\exp(\\frac{V+40}{10})} \\quad
\\alpha_h(V_m) = 0.07 \\exp(-\\frac{(V+65)}{20}) \\\\
\\beta_n(V_m) = 0.125 \\exp(-\\frac{V+65}{80}) \\quad
\\beta_m(V_m) = 4 \\exp(-\\frac{V+65}{18}) \\quad
\\beta_h(V_m) = \\frac{1}{1 + \\exp(-\\frac{V+35}{10})}
\\end{aligned}
```
The Nobel-winning four-dimensional dynamical system due to Hodgkin and Huxley [^HodgkinHuxley1952], which describes the electrical spiking
activity (action potentials) in neurons. A complete description of all parameters and variables is given in [^HodgkinHuxley1952], [^Ermentrout2010], and [^Abbott2005]. The equations and default parameters used here are taken from [^Ermentrout2010][^Abbott2005]. They differ slightly from the original paper [^HodgkinHuxley1952], since they were changed to shift the resting potential to -65 mV, instead of the 0mV in the original paper.

Varying the injected current I from `I = -5`  to  `I = 12` takes the neuron from quiescent to a single spike, and to a tonic (repetitive) spiking. This is due to a subcritical Hopf bifurcation, which occurs close to `I = 9.5`.

[^HodgkinHuxley1952] : A. L. Hodgkin, A.F. Huxley J. Physiol., pp. 500-544 (1952).

[^Ermentrout2010] : G. Bard Ermentrout, and David H. Terman, "Mathematical Foundations of Neuroscience", Springer (2010).

[^Abbott2005] : L. F. Abbott, and P. Dayan, "Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems", MIT Press (2005).
"""
function hodgkinhuxley(u0=[-60.0, 0.0, 0.0, 0.0]; I = 12.0, Vna = 50.0, Vk = -77.0, Vl = -54.4, gna = 120.0,gk = 36.0, gl = 0.3)
#In Ermentrout's & Abbott's books
    return CDS(hodgkinhuxley_rule, u0, [I, Vna, Vk, Vl, gna, gk, gl])
end
function hodgkinhuxley_rule(u, p, t)
    @inbounds begin
        I, Vna, Vk, Vl, gna, gk, gl = p
        V, n, m, h = u

        αn = 0.01 * (V+55)/(1 -  exp(-(V+55)/10))
        αm = 0.1 * (V+40)/(1- exp(-(V+40)/10))
        αh = 0.07 * exp(-(V+65)/20.0 )
        βn = 0.125 * exp(-(V+65)/80.0)
        βm = 4.0 * exp(-(V+65)/18.0 )
        βh = 1.0/(1 + exp(-(V+35)/10))

        du1 = (I -n^4*gk*(V-Vk) - m^3*h*gna*(V-Vna) - gl*(V-Vl))/1.0
        du2 = αn*(1-n) - βn*n
        du3 = αm*(1-m) - βm*m
        du4 = αh*(1-h) - βh*h
        return SVector{4}(du1, du2, du3, du4)
    end
end

"""
```julia
vanderpol(u0=[0.5, 0.0]; μ=1.5, F=1.2, T=10) -> ds
```
```math
\\begin{aligned}
\\ddot{x} -\\mu (1-x^2) \\dot{x} + x = F \\cos(\\frac{2\\pi t}{T})
\\end{aligned}
```
The forced van der Pol oscillator is an oscillator with a nonlinear damping term driven by a sinusoidal forcing. It was proposed by Balthasar van der Pol, in his studies of nonlinear electrical circuits used in the first radios [^Kanamaru2007][^Strogatz2015].
The unforced oscillator (`F = 0`) has stable oscillations in the form of a limit cycle with a slow buildup followed by a sudden discharge, which van der Pol called relaxation oscillations [^Strogatz2015][^vanderpol1926].
The forced oscillator (`F > 0`) also has periodic behavior for some parameters, but can additionally have chaotic behavior.

The van der Pol oscillator is a specific case of both the FitzHugh-Nagumo neural model [^Kanamaru2007]. The default damping parameter is taken from [^Strogatz2015] and the forcing parameters are taken from [^Kanamaru2007], which generate periodic oscillations. Setting `\\mu=8.53` generates chaotic oscillations.

[^Kanamaru2007] : Takashi Kanamaru (2007) "Van der Pol oscillator", Scholarpedia, 2(1):2202.

[^Strogatz2015] : Steven H. Strogatz (2015) "Nonlinear dynamics and chaos : with applications to physics, biology, chemistry, and engineering", Boulder, CO :Westview Press, a member of the Perseus Books Group.

[^vanderpol1926] : B. Van der Pol (1926), "On relaxation-oscillations", The London, Edinburgh and Dublin Phil. Mag. & J. of Sci., 2(7), 978–992.
"""
function vanderpol(u0=[0.5, 0.0]; μ=1.5, F=1.2, T=10)
    return CDS(vanderpol_rule, u0, [μ, F, T], vanderpol_jac)
end
function vanderpol_rule(u, p, t)
    @inbounds begin
        μ, F, T = p
        du1 = u[2]
        du2 = μ*(1 - u[1]^2)*u[2] - u[1] + F*sin(2π*t/T)
        return SVector{2}(du1, du2)
    end
end
function vanderpol_jac(u, p, t)
    @inbounds begin
        μ = p[1]
        J = @SMatrix [0 1;
        (-μ*(2*u[1]-1)*u[2]-1) (μ*(1 - u[1]^2))]
        return J
    end
end

"""
```julia
lotkavolterra(u0=[10.0, 5.0]; α = 1.5, β = 1, δ=1, γ=3) -> ds
```
```math
\\begin{aligned}
\\dot{x} &= \\alpha x - \\beta xy, \\\\
\\dot{y} &= \\delta xy - \\gamma y
\\end{aligned}
```
The famous Lotka-Volterra model is a simple ecological model describing the interaction between a predator and a prey species (or also parasite and host species). It has been used independently in fields such as epidemics, ecology, and economics [^Hoppensteadt2006], and is not to be confused with the Competitive Lotka-Volterra model, which describes competitive interactions between species.

The `x` variable describes the number of prey, while `y` describes the number of predator.  The default parameters are taken from [^Weisstein], which lead to typical periodic oscillations.

[^Hoppensteadt2006] : Frank Hoppensteadt (2006) "Predator-prey model", Scholarpedia, 1(10):1563.

[^Weisstein] : Weisstein, Eric W., "Lotka-Volterra Equations." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/Lotka-VolterraEquations.html
"""
function lotkavolterra(u0=[10.0, 5.0]; α = 1.5, β = 1, δ=1, γ=3)
    return CDS(lotkavolterra_rule, u0, [α, β, δ, γ], lotkavolterra_jac)
end
function lotkavolterra_rule(u, p, t)
    @inbounds begin
        α, β, δ, γ = p
        du1 = α*u[1] - β*u[1]*u[2]
        du2 = δ*u[1]*u[2] - γ*u[2]
        return SVector{2}(du1, du2)
    end
end
function lotkavolterra_jac(u, p, t)
    @inbounds begin
        α, β, δ, γ = p
        J = @SMatrix [(α - β*u[2]) (-β*u[1]);
        (δ*u[2]) (δ*u[1] - γ)]
        return J
    end
end

"""
```julia
hindmarshrose(u0=[-1.0, 0.0, 0.0]; a=1, b=3, c=1, d=5, r=0.001, s=4, xr=-8/5, I=2.0) -> ds
```
```math
\\begin{aligned}
\\dot{x} &= y - ax^3 + bx^2 +I - z, \\\\
\\dot{y} &= c - dx^2 -y, \\\\
\\dot{z} &= r(s(x - x_r) - z)
\\end{aligned}
```

The Hindmarsh-Rose model reproduces the bursting behavior of a neuron's membrane potential, characterized by a fast sequence of spikes followed by a quiescent period. The `x` variable describes the membane potential, whose behavior can be controlled by the applied current `I`; the `y` variable describes the sodium and potassium ionic currents, and `z` describes an adaptation current [^HindmarshRose1984].

The default parameter values are taken from [^HindmarshRose1984], chosen to lead to periodic bursting.

[^HindmarshRose1984] : J. L. Hindmarsh and R. M. Rose (1984) "A model of neuronal bursting using three coupled first order differential equations", Proc. R. Soc. Lond. B 221, 87-102.
"""
function hindmarshrose(u0=[-1.0, 0.0, 0.0]; a=1, b=3, c=1, d=5, r=0.001, s=4, xr=-8/5, I=2.0)
    return CDS(hindmarshrose_rule, u0, [a,b,c,d,r,s,xr, I], hindmarshrose_jac)
end
function hindmarshrose_rule(u, p, t)
    @inbounds begin
        a,b,c,d,r,s, xr, I = p
        du1 = u[2] - a*u[1]^3 + b*u[1]^2 -u[3] + I
        du2 = c - d*u[1]^2 - u[2]
        du3 = r*(s*(u[1] - xr) - u[3])
        return SVector{3}(du1, du2, du3)
    end
end
function hindmarshrose_jac(u, p, t)
    @inbounds begin
        a,b,c,d,r,s, xr, I = p
        J = @SMatrix [(-3*a*u[1]^2 + 2*b*u[1]) 1 -1;
        -2*d*u[1] -1 0;
        r*s 0 -r]
        return J
    end
end

"""
```julia
stuartlandau_oscillator(u0=[1.0, 0.0]; μ=1.0, ω=1.0, b=1) -> ds
```
The Stuart-Landau model describes a nonlinear oscillation near a Hopf bifurcation, and was proposed by Landau in 1944 to explain the transition to turbulence in a fluid [^Landau1944].
It can be written in cartesian coordinates as [^Deco2017]
```math
\\begin{aligned}
\\dot{x} &= (\\mu -x^2 -y^2)x - \\omega y - b(x^2+y^2)y \\\\
\\dot{y} &= (\\mu -x^2 -y^2)y + \\omega x + b(x^2+y^2)x
\\end{aligned}
```

The dynamical analysis of the system is greatly facilitated by putting it in polar coordinates, where it becomes the normal form of the supercritical Hopf bifurcation) [^Strogatz2015].
```math
\\begin{aligned}
\\dot{r} &= \\mu r - r^3, \\\\
\\dot{\\theta} &= \\omega +br^2
\\end{aligned}
```
The parameter `\\mu` serves as the bifurcation parameter, `\\omega` is the frequency of infinitesimal oscillations, and `b` controls the dependence of the frequency on the amplitude.  Increasing `\\mu` from negative to positive generates the supercritical Hopf bifurcation, leading from a stable spiral at the origin to a stable limit cycle with radius `\\sqrt(\\mu)`.
[^Landau1944] : L. D. Landau, "On the problem of turbulence, In Dokl. Akad. Nauk SSSR (Vol. 44, No. 8, pp. 339-349) (1944).
[^Deco2017] : G. Deco et al "The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core",  Sci Rep 7, 3095 (2017).
[^Strogatz2015] : Steven H. Strogatz "Nonlinear dynamics and chaos : with applications to physics, biology, chemistry, and engineering", Boulder, CO :Westview Press, a member of the Perseus Books Group (2015).
"""
function stuartlandau_oscillator(u0=[1.0, 0.0]; μ=1.0, ω=1.0, b=1)
    return CDS(stuartlandau_rule, u0, [μ, ω, b], stuartlandau_jac)
end
function stuartlandau_rule(u, p, t)
    @inbounds begin
        μ, ω, b = p
        du1 = u[1]*(μ - u[1]^2 - u[2]^2) - ω*u[2] - b*(u[1]^2 + u[2]^2)*u[2]
        du2 = u[2]*(μ - u[1]^2 - u[2]^2) + ω*u[1] + b*(u[1]^2 + u[2]^2)*u[1]
        return SVector{2}(du1, du2)
    end
end
function stuartlandau_jac(u, p, t)
    @inbounds begin
        μ, ω, b = p
        J = @SMatrix [(μ - 3*u[1]^2 -u[2]^2 -2*b*u[1]*u[2]) (-2*u[1]*u[2] -ω -b*u[1]^2 -3*b*u[2]^2);
            (-2*u[1]*u[2] +ω +b*u[2]^2 +3*b*u[1]^2) (μ -u[1]^2 -3*u[2]^2 +2*b*u[1]*u[2])]
        return J
    end
end


text:
module FancyDiagnostics
    include("LineNumbers.jl")
    include("display.jl")
    include("hooks.jl")
end # module

ext: jl


text:
__precompile__()

module SONames

const SONAMES = Ref(Dict{String,String}())

"""
    SONames.read()

Reload a mapping of system shared libraries, with keys as names of the form `libXX` and
values as the fully qualified paths of the form `/path/to/libXX.so.ver`.
Full paths can then be looked up with `SONames.lookup("libXX")`.
"""
function read()
    global SONAMES
    empty!(SONAMES[])
    @static if Sys.islinux()
        cmd = `/sbin/ldconfig -p`
    else
        cmd = `/sbin/ldconfig -r`
    end
    for (n, line) in enumerate(eachline(cmd, keep=false))
        # The first two lines aren't relevant
        n > 2 || continue
        @static if Sys.islinux()
            # Idk, I'm writing this on FreeBSD
        else
            inds = findfirst(":-l", line)
            inds === nothing && continue
            dotind = findnext(==('.'), line, last(inds))
            dotind === nothing && continue
            spc = findlast(==(' '), line)
            spc === nothing && continue

            name = string("lib", line[nextind(line, last(inds)):prevind(line, dotind)])
            sopath = line[nextind(line, spc):end]
        end
        push!(SONAMES[], name => sopath)
    end
    return
end

"""
    SONames.lookup(lib)

Look up the fully qualified path for the given shared library name.

# Examples
```julia-repl
julia> SONames.lookup("libc")
"/lib/libc.so.7"
```
"""
lookup(lib::AbstractString) = (global SONAMES; get(SONAMES[], lib, ""))

function __init__()
    if !Sys.islinux() && Sys.KERNEL !== :FreeBSD
        error("SO name mapping is only available on Linux and FreeBSD")
    end
    read()
end

end # module


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

name = "CmdStanBinaries"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/stan-dev/cmdstan/releases/download/v2.27.0/cmdstan-2.27.0.tar.gz", "ff71f4d255cf26c2d8366a8173402656a659399468871305056aa3d56329c1d5")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd cmdstan-2.27.0/
echo 'CXX=c++' >> make/local
echo 'TBB_CXX_TYPE=c++' >> make/local
echo 'STAN_THREADS=true' >> make/local
cd stan/lib/stan_math/lib/tbb_2020.3/build/
cp linux.inc x86_64-linux-gnu.inc
cp linux.gcc.inc x86_64-linux-gnu.c++.inc
cd ../../../../../..
make build
mkdir ${bindir}
mkdir ${libdir}

cp bin/stanc ${bindir}/stanc
cp bin/diagnose ${bindir}/diagnose
cp bin/print ${bindir}/print
cp bin/stansummary ${bindir}/stansummary

exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc")
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("stanc", :stanc),
    ExecutableProduct("stansummary", :stansummary),
    ExecutableProduct("print", :print_deprecated),
    ExecutableProduct("diagnose", :diagnose)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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


text:
"""
`abstract ImageSource`

Description
============

Abstract supertype for all image sources. Every subtype of
`ImageSource` must implement the read-only functions of the array
interface. Additionally, they must fully implement the iterator
interface.

Methods
========

- **`rand`** : Loads and returns a random image, or an array of
images (depending on the parameters).

- **`length`** : Returns the total number of images that are
available in the image source.

- **`getindex`** : Loads and returns the image of the given index.

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

see also
=========

`rand`, `length`

`DirImageSource`
"""
abstract ImageSource

function Base.rand(s::ImageSource)
    index = rand(1:length(s))
    s[index]
end

function Base.rand(s::ImageSource, n::Integer)
    idx = rand(1:length(s), n)
    s[idx]
end

function Base.getindex(s::ImageSource, idx::AbstractVector)
    [s[i] for i in idx]
end

Base.eltype{T<:ImageSource}(::Type{T}) = Image
Base.endof(s::ImageSource) = length(s)

function Base.show(io::IO, s::ImageSource)
    println(io, typeof(s))
    print(io, "- $(length(s)) images found")
end


text:
project_root = "$(@__DIR__)/.."

include("$project_root/julia/code_gen.jl")

function run_example()
    generate("$project_root/example/event_registration.orm",
             "$project_root/example/output")
end

run_example()

max_stars_repo_name: RelationalAI-oss/kse-tools

max_issues_repo_name: RelationalAI-oss/kse-tools

max_forks_repo_name: RelationalAI-oss/kse-tools


text:
struct NewtonResults
    z::Vector
    λ::Vector
    ν::Vector

    W::SparseMatrixCSC
    w::Vector
    G::SparseMatrixCSC
    g::Vector
    D::SparseMatrixCSC
    d::Vector
    A::SparseMatrixCSC
    b::Vector
    active_set::Vector
    x_eval::Vector
    δ::Vector
end

function NewtonResults(Nz::Int,Np::Int,Nx::Int)
    z = zeros(Nz)
    λ = zeros(Np)
    ν = zeros(Nx)

    W = spzeros(Nz,Nz)
    w = zeros(Nz)

    G = spzeros(Np,Nz)
    g = zeros(Np)

    D = spzeros(Nx,Nz)
    d = zeros(Nx)

    A = spzeros(Nz+Np+Nx,Nz+Np+Nx)
    b = spzeros(Nz+Np+Nx)

    active_set = zeros(Bool,Np)
    x_eval = zeros(Nx)

    δ = zeros(Nz+Np+Nx)

    NewtonResults(z,λ,ν,W,w,G,g,D,d,A,b,active_set,x_eval,δ)
end

function NewtonResults(solver::Solver)
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    # batch problem dimensions
    nm = nn + mm
    Nz = nn*N + mm*(N-1)
    Np = p*(N-1) + p_N
    Nx = N*nn

    NewtonResults(Nz,Np,Nx)
end

function get_batch_sizes(solver::Solver)
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    # batch problem dimensions
    nm = nn + mm
    Nz = nn*N + mm*(N-1)
    Np = p*(N-1) + p_N
    Nx = N*n
    Nu = mm*(N-1) # number of control decision variables u

    return Nz,Np,Nx,Nu,nm
end

## Newton solve
function update_newton_results!(newton_results::NewtonResults,results::SolverIterResults,solver::Solver)
    # get problem dimensions
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)
    Nz,Np,Nx,Nu,nm = get_batch_sizes(solver)

    # pull out results for convenience
    X = results.X
    U = results.U
    Iμ = results.Iμ
    C = results.C
    Cx = results.Cx
    Cu = results.Cu
    fdx = results.fdx
    fdu = results.fdu
    x0 = solver.obj.x0

    # pull out newton results for convenience
    z = newton_results.z
    λ = newton_results.λ
    ν = newton_results.ν

    W = newton_results.W
    w = newton_results.w

    G = newton_results.G
    g = newton_results.g

    D = newton_results.D
    d = newton_results.d

    active_set = newton_results.active_set
    x_eval = newton_results.x_eval

    Inn = sparse(Matrix(I,nn,nn))

    # update batch matrices
    for k = 1:N
        x = results.X[k]

        if k != N
            u = results.U[k]
            Q,R,H,q,r = taylor_expansion(solver.obj.cost,x,u) .* solver.dt
        else
            Qf,qf = taylor_expansion(solver.obj.cost,x)
        end

        # Indices
        k != N ? idx = (((k-1)*nm + 1):k*nm) : idx = (((k-1)*nm + 1):Nz) # index over x and u
        k != N ? idx2 = (((k-1)*p + 1):k*p) : idx2 = (((k-1)*p + 1):Np) # index over p
        idx3 = ((k-1)*nm + 1):((k-1)*nm + nn) # index over x only
        idx4 = ((k-1)*nm + nn + 1):k*nm # index over u only
        idx5 = ((k-1)*mm + 1):k*mm # index through stacked u vector
        idx6 = ((k-1)*p + 1):((k-1)*p + pI) # index over p only inequality indices
        idx7 = ((k-1)*nn + 1):(k*nn)
        idx8 = ((k-2)*nm + 1):((k-1)*nm + nn)
        idx9 = ((k-2)*nm + 1):((k-2)*nm + nn)
        idx10 = ((k-2)*nm + nn + 1):((k-1)*nm)

        # Assemble W, w, G, g
        if k != N
            # W[idx,idx] = [Q H';H R] + [Cx[k]'*Iμ[k]*Cx[k] Cx[k]'*Iμ[k]*Cu[k]; Cu[k]'*Iμ[k]*Cx[k] Cu[k]'*Iμ[k]*Cu[k]]
            W[idx3,idx3] = Q + Cx[k]'*Iμ[k]*Cx[k]
            W[idx3,idx4] = H' + Cx[k]'*Iμ[k]*Cu[k]
            W[idx4,idx3] = H + Cu[k]'*Iμ[k]*Cx[k]
            W[idx4,idx4] = R + Cu[k]'*Iμ[k]*Cu[k]

            # w[idx] = [q;r] + [Cx[k]'*Iμ[k]*C[k]; Cu[k]'*Iμ[k]*C[k]]
            w[idx3] = q + Cx[k]'*Iμ[k]*C[k]
            w[idx4] = r + Cu[k]'*Iμ[k]*C[k]

            # G[idx2,idx] = [Cx[k] Cu[k]]
            G[idx2,idx3] = Cx[k]
            G[idx2,idx4] = Cu[k]

            g[idx2] = C[k]

            z[idx] = [x;u]
        else
            W[idx,idx] = Qf + Cx[N]'*Iμ[N]*Cx[N]
            w[idx] = qf + Cx[N]'*Iμ[N]*C[N]

            G[idx2,idx] = results.Cx[N]
            g[idx2] = results.C[N]

            z[idx] = x
        end

        # assemble D, d
        if k == 1
            D[1:nn,1:nn] = Inn
            d[1:nn] = X[1] - solver.obj.x0
        else
            # D[idx7,idx8] = [-fdx[k-1] -fdu[k-1] Inn]
            D[idx7,idx9] = -fdx[k-1]
            D[idx7,idx10] = -fdu[k-1]
            D[idx7,idx3] = Inn

            solver.fd(view(x_eval,idx7),X[k-1][1:n],U[k-1][1:m])
            d[idx7] = X[k] - x_eval[idx7]
        end

        # assemble λ, ν, active_set
        λ[idx2] = results.λ[k]
        ν[idx7] = results.s[k] # TODO is this still ok for penalty method?
        active_set[idx2] = results.active_set[k]

    end

    return nothing
end

function update_results_from_newton_results!(results::SolverIterResults,newton_results::NewtonResults,solver::Solver)
    n,m,N = get_sizes(solver)
    n̄,nn = get_num_states(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

    # batch problem dimensions
    Nz,Np,Nx,Nu,nm = get_batch_sizes(solver)

    z = newton_results.z
    λ = newton_results.λ
    active_set = newton_results.active_set

    # update results with stacked vector from NewtonResults
    for k = 1:N
        idx = ((k-1)*nm + 1):k*nm # index over x and u
        k != N ? idx2 = (((k-1)*p + 1):k*p) : idx2 = (((k-1)*p + 1):Np) # index over p
        k != N ? idx3 = (((k-1)*nm + 1):((k-1)*nm + nn)) : idx3 = (((k-1)*nm + 1):Nz)# index over x only
        idx4 = ((k-1)*nm + nn + 1):k*nm # index over u only
        idx5 = ((k-1)*mm + 1):k*mm # index through stacked u vector
        idx6 = ((k-1)*p + 1):((k-1)*p + pI) # index over p only inequality indices

        results.X[k] = z[idx3]
        k != N ? results.U[k] = z[idx4] : nothing
        results.λ[k] = λ[idx2]
    end

    update_constraints!(results,solver)
    update_jacobians!(results,solver)

    return nothing
end

function newton_step!(newton_results::NewtonResults,alpha::Float64=1.0)
    W = newton_results.W
    w = newton_results.w

    G = newton_results.G
    g = newton_results.g

    D = newton_results.D
    d = newton_results.d

    A = newton_results.A
    b = newton_results.b

    z = newton_results.z
    λ = newton_results.λ
    ν = newton_results.ν

    δ = newton_results.δ

    # get batch problem sizes
    Nz = size(W,1)
    Np = size(G,1)
    Nx = size(D,1)

    active_set = newton_results.active_set

    Np_as = sum(active_set)

    # assemble KKT matrix/vector
    _idx1 = Array(1:Nz)
    _idx2 = Array((1:Np) .+ Nz)[active_set]
    _idx3 = Array((1:Nx) .+ (Nz + Np))
    _idx4 = [_idx1;_idx2;_idx3]

    A[_idx1,_idx1] = W
    A[_idx1,_idx2] = G[active_set,:]'
    A[_idx1,_idx3] = D'
    A[_idx2,_idx1] = G[active_set,:]
    A[_idx3,_idx1] = D

    b[_idx1] = -w
    b[_idx2] = -g[active_set]
    b[_idx3] = -d

    # solve
    δ[_idx4] = A[_idx4,_idx4]\b[_idx4]

    # update
    z .+= alpha*δ[1:Nz]
    λ .+= alpha*δ[(1:Np) .+ Nz]
    ν .+= alpha*δ[(1:Nx) .+ (Nz+Np)]

    return nothing
end

function cost_newton(results::SolverIterResults,newton_results::NewtonResults,solver::Solver)
    # get problem dimensions
    n,m,N = get_sizes(solver)

    J = cost(solver,results)

    # add dynamics constraint costs
    X = results.X
    U = results.U
    ν = newton_results.ν
    x_eval = newton_results.x_eval

    for k = 1:N
        if k == 1
            J += ν[1:n]'*(X[1] - solver.obj.x0)
        else
            idx = ((k-1)*n+1):k*n
            solver.fd(view(x_eval,idx),X[k-1][1:n],U[k-1][1:m])
            J += ν[idx]'*(X[k] - x_eval[idx])
        end
    end
    return J
end

function newton_step!(results::SolverIterResults,newton_results::NewtonResults,solver::Solver,alpha::Float64=1.0)
    update_newton_results!(newton_results,results,solver)
    newton_step!(newton_results,alpha)
    update_results_from_newton_results!(results,newton_results,solver)
    return nothing
end

function newton_solve!(results::SolverIterResults,solver::Solver)
    # instantiate Newton Results
    newton_results = NewtonResults(solver)

    # copy current results
    results_new = copy(results)

    # get initial cost and max constraint violation
    update_newton_results!(newton_results,results_new,solver)
    J_prev = cost_newton(results_new,newton_results,solver)
    c_max = max_violation(results_new)

    println("Newton initialization")
    println("Cost: $J_prev")
    println("c_max: $c_max")

    α = 1.0
    max_iter = 10
    max_c = 1e-8
    iter = 1
    ls_param = 0.01

    while c_max > max_c && iter <= max_iter
        println("$α")
        newton_step!(results_new,newton_results,solver,α)
        J = cost_newton(results_new,newton_results,solver)
        c_max = max_violation(results_new)
        println("Newton step: $iter")
        println("Cost: $J \n    prev. Cost: $J_prev")
        println("c_max: $c_max")

        if J <= J_prev + ls_param*newton_results.b'*α*newton_results.δ
            results = copy(results_new)
            J_prev = copy(J)
            println("improved")
            α = 1.0
        else
            println("what??")
            results_new = copy(results)
            α = 0.5*α
        end
        iter += 1
    end
    println("Newton solver end")
    println("c_max: $c_max")
end


text:
let
    n = 5
    tf = 5
    t = collect(0:(tf/n):tf)
    f(t) = sin.(t)
    fdot(t) = cos.(t)

    x = f(t)
    spline = CubicSpline(t,x)
    interpolate(spline, t[1])
    deriv(spline, t[1])
    hess(spline, t[1])
    interpolate(spline, [t[1], t[end]])
    deriv(spline, [t[1], t[end]])
    hess(spline, [t[1], t[end]])

    dt = 1.0
    @test extrapolate(spline, t[end]+dt, 0)[1] == interpolate(spline, t[end])
    @test extrapolate(spline, t[end]+dt, 1)[1] == interpolate(spline, t[end]) + deriv(spline, t[end])*dt
    @test extrapolate(spline, t[1]-dt, 0)[1] == interpolate(spline, t[1])
    @test extrapolate(spline, t[1]-dt, 1)[1] == interpolate(spline, t[1]) - deriv(spline, t[1])*dt

    xdot = fdot(t)
    spline2 = CubicSpline(t,x,xdot) # prescribe knot tangents directly

    # plt = plot(t_vec, interpolate(spline, t_vec))
    # plot!(plt, t_vec, interpolate(spline2, t_vec))
    # plot!(plt, t_vec, f(t_vec))
    # plot!(plt, t, x, seriestype=:scatter)
end
let
    n = 10
    t = collect(0:(1/10):1)
    x = sin.(t)
    y = sin.(t)
    X = hcat(x,y)
    spline = ParametricSpline(t,X)
    interpolate(spline, t[1])
    deriv(spline, t[1])
    hess(spline, t[1])
    interpolate(spline, [t[1], t[end]])
    deriv(spline, [t[1], t[end]])
    hess(spline, [t[1], t[end]])
end
# test turn rate and tangent rate
let
    # circle in the plane
    tvec = collect(0:.01:2π)
    x = cos.(tvec)
    y = sin.(tvec)
    dx = -y
    dy = x
    ddx = -x
    ddy = -y
    z = zeros(length(tvec))
    dz = z
    ddz = z
    r = collect(hcat(x,y,z))
    dr = collect(hcat(dx,dy,dz))
    ddr = collect(hcat(ddx,ddy,ddz))
    spline = ParametricSpline(tvec,r,dr)

    # @show tangent_rate(spline,0.0)
    # @show tangent_rate(spline,π/4)
    #
    # tvec2 = (1/2)*(tvec[2:end].+tvec[1:end-1])
    # pts = interpolate(spline,tvec2)
    # plt = plot(map(pt->pt[1],pts),map(pt->pt[2],pts),aspect_ratio=:equal)
    # pts = deriv(spline,tvec2)
    # plot!(plt,map(pt->pt[1],pts),map(pt->pt[2],pts),aspect_ratio=:equal)
    # pts = hess(spline,tvec2)
    # plot(map(pt->pt[1],pts),map(pt->pt[2],pts),aspect_ratio=:equal)
end


text:
import MathOptInterfaceMosek
using Base.Test
using JuMP
const MOI = MathOptInterface

include("jump_sdp.jl")
include("jump_lp.jl")
include("jump_soc.jl")

@testset "JuMP tests" begin
    test_jump_lp(MathOptInterfaceMosek.MosekSolver(LOG=0))
    test_jump_soc(MathOptInterfaceMosek.MosekSolver(LOG=0))
    test_jump_sdp(MathOptInterfaceMosek.MosekSolver(LOG=0))
end


text:
module TestShap

using ShapML
using Test
using DataFrames
using Random
using Statistics
using Distributed

@testset "shap return values are correct." begin

    #--------------------------------------------------------------------------
    data = DataFrame(y = 1:10, x1 = 1:10, x2 = 1:10, x3 = 1:10)

    function model_mean(data)
        y_pred = Array{Float64}(undef, size(data, 1))
        for i in 1:size(data, 1)
            y_pred[i] = mean(data[i, 2:end])
        end
        return DataFrame(y_pred = y_pred)
    end

    function predict_function(model, data)
      data_pred = model(data)
      return data_pred
    end

    data_shap = ShapML.shap(explain = data[:, 2:end],
                            model = model_mean,
                            predict_function = predict_function,
                            sample_size = 10
                            )

    #--------------------------------------------------------------------------
    # Test that each instance has a Shapley value for each feature.
    @test size(data_shap, 1) == (size(data, 1) * (size(data, 2) - 1))
    #--------------------------------------------------------------------------
    # Test that the Shapley values are the same for a subset of target features.
    data_shap_feature = ShapML.shap(explain = data[:, 2:end],
                                    model = model_mean,
                                    predict_function = predict_function,
                                    target_features = ["x2"],
                                    sample_size = 10
                                    )

    @test all(data_shap[data_shap.feature_name .== "x2", [:shap_effect]].shap_effect .== data_shap_feature[data_shap_feature.feature_name .== "x2", [:shap_effect]].shap_effect)
    #--------------------------------------------------------------------------
end
#------------------------------------------------------------------------------
@testset "The NULL model returns 0 Shapley values." begin

    #--------------------------------------------------------------------------
    data = DataFrame(y = 1:10, x1 = 1:10, x2 = 1:10)

    function model_mean(data)
        return DataFrame(y_pred = repeat([0], size(data, 1)))
    end

    function predict_function(model, data)
      data_pred = model(data)
      return data_pred
    end

    data_shap = ShapML.shap(explain = data[:, 2:end],
                            model = model_mean,
                            predict_function = predict_function,
                            sample_size = 10
                            )

    @test all(data_shap[:, [:shap_effect]].shap_effect .== 0)
    @test all(data_shap[:, [:shap_effect_sd]].shap_effect_sd .== 0)
    #--------------------------------------------------------------------------
end
end  # End test module.

max_stars_repo_name: nredell/shapML

max_issues_repo_name: nredell/shapML

max_forks_repo_name: nredell/shapML


text:
export kl_divergence

# P is truth, Q is approximation
function kl_divergence(P, Q)

    kld = 0.0

    p_sum = sum(P)
    q_sum = sum(Q)

    for i = 1:length(P)
        if P[i] > 0.0
            pn = P[i] / p_sum
            qn = Q[i] / q_sum

            kld += pn * log(pn / qn)
        end
    end
    return kld
end


text:
module TiffImages

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
    @eval Base.Experimental.@compiler_options optimize=1
end

using ColorTypes
using DocStringExtensions
using FileIO
using FixedPointNumbers
using IndirectArrays
using OffsetArrays
using OrderedCollections
using PkgVersion
using ProgressMeter
using Base.Iterators

const PKGVERSION = @PkgVersion.Version 0

include("enum.jl")
include("utils.jl")
include("files.jl")
include("compression.jl")
include("tags.jl")
include("ifds.jl")
include("layout.jl")
include(joinpath("types", "common.jl"))
include(joinpath("types", "dense.jl"))
include(joinpath("types", "mmap.jl"))
include("load.jl")

@deprecate TiffFile(::Type{O}) where O<:Unsigned TiffFile{O}()

## Precompilation helper
mktemp() do fpath, _
    for t in [N0f8, N0f16, Float32, Float64]
        for c in [Gray, GrayA, RGB, RGBA]
            TiffImages.save(fpath, rand(c{t}, 2, 2))
            TiffImages.load(fpath)
        end
    end
end

end # module


text:
# development testing for multithreaded convolution

using Distributions
using IncrementalInference
using Base: Test

import IncrementalInference: getSample

struct Prior{T} <: IncrementalInference.FunctorSingleton where T <: Distribution
  z::T
end
getSample(s::Prior, N::Int=1) = (reshape(rand(s.z,N),1,:), )
struct LinearOffset{T} <: IncrementalInference.FunctorPairwise where T <: Distribution
  z::T
end
getSample(s::LinearOffset, N::Int=1) = (reshape(rand(s.z,N),1,:), )
function (s::LinearOffset)(res::Array{Float64},
                           userdata::FactorMetadata,
                           idx::Int,
                           meas::Tuple{Array{Float64, 2}},
                           X1::Array{Float64,2},
                           X2::Array{Float64,2}  )
  #
  res[1] = meas[1][idx] - (X2[1,idx] - X1[1,idx])
  nothing
end
struct MultiModalOffset <: IncrementalInference.FunctorPairwise
  z::Vector{Distribution}
  c::Categorical
end
getSample(s::MultiModalOffset, N::Int=1) = (reshape.(rand.(s.z, N),1,:)..., rand(s.c, N))
function (s::MultiModalOffset)(res::Array{Float64},
                               userdata::FactorMetadata,
                               idx::Int,
                               meas::Tuple,
                               X1::Array{Float64,2},
                               X2::Array{Float64,2}  )
  #
  res[1] = meas[meas[end][idx]][idx] - (X2[1,idx] - X1[1,idx])
  nothing
end

@testset "Basic ContinuousScalar example to ensure multithreaded convolutions work..." begin

@show Threads.nthreads()

N = 100

# Start with an empty factor graph
fg = emptyFactorGraph()

# add the first node
addNode!(fg, :x0, ContinuousScalar, N=N)

# this is unary (prior) factor and does not immediately trigger autoinit of :x0.
addFactor!(fg, [:x0], Prior(Normal(0,1)))

addNode!(fg, :x1, ContinuousScalar, N=N)
# P(Z | :x1 - :x0 ) where Z ~ Normal(10,1)
addFactor!(fg, [:x0, :x1], LinearOffset(Normal(10.0,1)))

pts = approxConv(fg, :x0x1f1, :x1, N=N)

@test 0.95*N <= sum(abs.(pts - 10.0) .< 5.0)

end

#


text:
"""
Author: Eric Ford and collaborators (see header of each file)
Created: August 2020
Contact: https://github.com/RvSpectML/RvSpectMLBase.jl
"""

"""
[**RvSpectMLBase.jl**](https://github.com/RvSpectML/RvSpectMLBase.jl) is the base package for a set of complementary Julia package's for measuring radial velocities (RVs) from stellar spectroscopic timeseries via machine learning (ML).

"""

module RvSpectMLBase

# Packages we'll use in many places, but try to keep these minimal.
using LinearAlgebra, Statistics
using DataFrames, Query
using NaNMath

# types.jl is responsible for exporting its own types and functions
include("types/types.jl")

# Load basic physics early, since may be used by instruments
include("util/physics.jl")
export calc_doppler_factor, apply_doppler_boost!
export absorption_line
export speed_of_light_mps

# Load basic CS algorithms early?
#include("util/cs_algs.jl")
#export searchsortednearest

# instruments.jl & the instruments it contains are responsible for exporting their own functions & modules
include("instruments/instruments.jl")

# Load the rest of the utility functions
include("util/util.jl")
#export searchsortednearest

end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
using AIBECS, Test
using SparseArrays, LinearAlgebra
using Unitful, UnitfulAstro
using WorldOceanAtlasTools
using DiffEqBase
using ForwardDiff
using DataFrames
# For CI, make sure the download does not hang
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

Circulation = Archer_etal_2000
@testset "Archer_etal_2000 3-box model" begin
    # Run tests with the 6-box toy model
    include("setup.jl")
end

Circulation = Primeau_2x2x2
@testset "Primeau's 6-box model" begin
    # Run tests with the 6-box toy model
    include("setup.jl")
    include("parameters.jl")
    include("bgc_functions.jl")
    include("gridtools.jl")
    include("cost_functions.jl")
    include("solvers.jl")
    include("time_steps.jl")
    include("derivatives.jl")
end

Circulation = OCIM1
@testset "OCIM1" begin
    include("setup.jl")
end

Circulation = OCIM0
@testset "OCIM0.1" begin
    include("setup.jl")
end


text:
n = parse(Int, readline())
a = parse.(Int, split(readline()))
maxnum = typemin(Int64)
minnum = typemax(Int64)
sumnum = 0
for num in a
    global maxnum = max(maxnum, num)
    global minnum = min(minnum, num)
    global sumnum += num
end
println(minnum, " ", maxnum, " ", sumnum)


text:
function isPrime(x)
	i = 2
	if x == 1
		return false
	end
	if x == 2
		return true
	end
	if x == 3
		return true
	end
	while i <= sqrt(x)
		if x % i == 0
		return false
		end
		i =i+1
	end
	return true
end

function nthPrime(n)
	primeCount = 0
	if n == 0
		return "there is no zeroth prime"
	end
	if typeof(n)==String
		n = parse(Int64,n)
	end
	x=2
	while primeCount <= n
		if isPrime(x)
			primeCount += 1
		end
		if primeCount == n
			return x
		end
		x+=1
	end
end

function main(input)
	nthPrime(input)
end

max_stars_repo_name: KrymoreAkbota/julia

max_issues_repo_name: KrymoreAkbota/julia

max_forks_repo_name: KrymoreAkbota/julia


text:
# broadcast.jl: Elementwise broadcasting binary functions for arrays and scalars.
# The entry format is (cudaname, julianame, kernelcode)
# With single name entries cudaname=julianame and kernelcode=name(xi,yi).
# I commented out functions if I don't know the Julia equivalent.

broadcast_ops = [
("add",".+","xi+yi"),
("sub",".-","xi-yi"),
("mul",".*","xi*yi"),
("div","./","xi/yi"),
("pow",".^","pow(xi,yi)"),
("max","max","(xi>yi?xi:yi)"),
("min","min","(xi<yi?xi:yi)"),
("eq",".==","xi==yi"),
("ne",".!=","xi!=yi"),
("gt",".>","xi>yi"),
("ge",".>=","xi>=yi"),
("lt",".<","xi<yi"),
("le",".<=","xi<=yi"),
# "hypot",
# "rhypot",
# "atan2",
# "frexp",
# "ldexp",
# "scalbn",
# "scalbln",
# "jn",
# "yn",
# "fmod",
# "remainder",
# "mod",
# "fdim",
("invxback","invxback","(-xi*yi*yi)"),
("reluback","reluback","(yi>0?xi:0)"),
("sigmback","sigmback","(xi*yi*(1-yi))"),
("tanhback","tanhback","(xi*(1-yi*yi))"),
("rpow","rpow","pow(yi,xi)"),   # need this for Array.^Scalar
"BGH",
]

# broadcast_op overloads a Julia function for KnetArrays.
# The corresponding kernel is defined in libknet8.
function broadcast_op(f, j=f, o...)
    J=Symbol(j)
    if isdefined(Base, J); eval(Expr(:import,:Base,J)); end
    for S in (32,64)
        T = Symbol("Float$S")
        F01 = "$(f)_$(S)_01"    # Scalar,Array->Array
        F11 = "$(f)_$(S)_11"    # Array,Array->Array (same size)
        F12 = "$(f)_$(S)_12"    # Array,Array->Array (different size)
        @eval begin
            function $J(x::$T,y::KnetArray{$T})
                z = similar(y)
                ccall(($F01,$libknet8),Void,(Cint,$T,Ptr{$T},Ptr{$T}),length(z),x,y,z)
                return z
            end
            function $J(x::KnetArray{$T},y::KnetArray{$T})
                if size(x)==size(y)
                    z = similar(x)
                    ccall(($F11,$libknet8),Void,(Cint,$Ptr{$T},Ptr{$T},Ptr{$T}),length(z),x,y,z)
                    return z
                else
                    (dz,sx,nx,sy,ny) = vbroadcast_shape(x,y)
                    z = similar(x,dz)
                    ccall(($F12,$libknet8),Void,(Cint,$Ptr{$T},Cint,Cint,Ptr{$T},Cint,Cint,Ptr{$T}),length(z),x,sx,nx,y,sy,ny,z)
                    return z
                end
            end
        end
    end
end

# vbroadcast_shape computes index/offset arguments for a broadcasting kernel call.
function vbroadcast_shape(x,y)
    nz = max(ndims(x),ndims(y))
    dz = ones(Int,nz)
    xdims = ydims = xsame = ysame = xlast = ylast = 0; zlen = 1
    for i=1:nz
        if size(x,i) > 1
            xdims += 1; xlast = i
            dz[i] = size(x,i)
        end
        if size(y,i) > 1
            ydims += 1; ylast = i
            if dz[i] == 1
                dz[i] = size(y,i)
            else
                dz[i] == size(y,i) || throw(DimensionMismatch("arrays could not be broadcast to a common size"))
            end
        end
        xsame += (dz[i] == size(x,i))
        ysame += (dz[i] == size(y,i))
        zlen *= dz[i]
    end
    xsame == nz || xdims <= 1 || error("Only vector broadcasting supported")
    ysame == nz || ydims <= 1 || error("Only vector broadcasting supported")
    if xdims == 0
        sx = zlen; nx = 1
    elseif xdims == 1
        sx = prod(dz[1:xlast-1]); nx=dz[xlast]
    elseif xsame == nz
        sx = 1; nx=zlen
    else
        error("Broadcasting error")
    end
    if ydims == 0
        sy = zlen; ny = 1
    elseif ydims == 1
        sy = prod(dz[1:ylast-1]); ny=dz[ylast]
    elseif ysame == nz
        sy = 1; ny=zlen
    else
        error("Broadcasting error")
    end
    return (tuple(dz...), sx, nx, sy, ny)
end

# Define all overloaded Julia functions for KnetArrays:

for f in broadcast_ops
    if !isa(f,Tuple); f=(f,); end
    broadcast_op(f...)
end

# Additional imports: fns in broadcast_ops are imported in broadcast_op()
import Base: +, -, *, /, \

# Here we'll just define some functions that specifically do not have broadcasting.
(+){T}(x::KnetArray{T},y::KnetArray{T})=(size(x)==size(y)||throw(DimensionMismatch("$(map(size,(x,y)))"));(.+)(x,y))
(-){T}(x::KnetArray{T},y::KnetArray{T})=(size(x)==size(y)||throw(DimensionMismatch("$(map(size,(x,y)))"));(.-)(x,y))
#(*){T}(x::KnetArray{T},y::KnetArray{T})=(.*)(x,y) # This is matmul
#(/){T}(x::KnetArray{T},y::KnetArray{T})=(./)(x,y) # This is another linalg op

# Broadcast max/min haven't been defined in Base:
# max(a::Array,b::Array)=broadcast(max,a,b)
# min(a::Array,b::Array)=broadcast(min,a,b)
# tkelman: These two methods aren't necessary, and overwrite Base. You can get this behavior via max.(a,b), with @compat needed on 0.4.

# Ambiguity fixes:
max{T<:Real,S<:Real}(a::KnetArray{T},s::S)=max(T(s),a)
max{T<:Real,S<:Real}(s::S,a::KnetArray{T})=max(T(s),a)
min{T<:Real,S<:Real}(a::KnetArray{T},s::S)=min(T(s),a)
min{T<:Real,S<:Real}(s::S,a::KnetArray{T})=min(T(s),a)
max{T<:Real,S<:Number}(a::KnetArray{T},s::S)=max(T(s),a)
max{T<:Real,S<:Number}(s::S,a::KnetArray{T})=max(T(s),a)
min{T<:Real,S<:Number}(a::KnetArray{T},s::S)=min(T(s),a)
min{T<:Real,S<:Number}(s::S,a::KnetArray{T})=min(T(s),a)
(+)(a::KnetArray{Bool},s::Bool)=(.+)(s,a)
(+)(s::Bool,a::KnetArray{Bool})=(.+)(s,a)
(-)(a::KnetArray{Bool},s::Bool)=(.+)(-s,a)
(-)(s::Bool,a::KnetArray{Bool})=(.-)(s,a)
(.^)(x::Base.Irrational{:e}, a::KnetArray)=(.^)(float(x),a)

# Scalar kernels are defined for scalar,array order only.
# For array,scalar we can get most for free.
(.+){T}(a::KnetArray{T},s::Number)=(.+)(T(s),a)
(.+){T}(s::Number,a::KnetArray{T})=(.+)(T(s),a)
(.-){T}(a::KnetArray{T},s::Number)=(.+)(T(-s),a)
(.-){T}(s::Number,a::KnetArray{T})=(.-)(T(s),a)
(.*){T}(a::KnetArray{T},s::Number)=(.*)(T(s),a)
(.*){T}(s::Number,a::KnetArray{T})=(.*)(T(s),a)
(./){T}(a::KnetArray{T},s::Number)=(.*)(T(1/s),a)
(./){T}(s::Number,a::KnetArray{T})=(./)(T(s),a)
max{T}(a::KnetArray{T},s::Number)=max(T(s),a)
max{T}(s::Number,a::KnetArray{T})=max(T(s),a)
min{T}(a::KnetArray{T},s::Number)=min(T(s),a)
min{T}(s::Number,a::KnetArray{T})=min(T(s),a)
(.^){T}(s::Number,a::KnetArray{T})=(.^)(T(s),a)
# Pow is the one exception, we need to define a separate kernel:
(.^){T}(a::KnetArray{T},s::Number)=rpow(T(s),a)

.=={T}(a::KnetArray{T},s::Number)=(T(s).==a)
.=={T}(s::Number,a::KnetArray{T})=(T(s).==a)
.!={T}(a::KnetArray{T},s::Number)=(T(s).!=a)
.!={T}(s::Number,a::KnetArray{T})=(T(s).!=a)
.>{T}(a::KnetArray{T},s::Number)=(T(s).<a)
.>{T}(s::Number,a::KnetArray{T})=(T(s).>a)
.>={T}(a::KnetArray{T},s::Number)=(T(s).<=a)
.>={T}(s::Number,a::KnetArray{T})=(T(s).>=a)
.<{T}(a::KnetArray{T},s::Number)=(T(s).>a)
.<{T}(s::Number,a::KnetArray{T})=(T(s).<a)
.<={T}(a::KnetArray{T},s::Number)=(T(s).>=a)
.<={T}(s::Number,a::KnetArray{T})=(T(s).<=a)

# familiar aliases for broadcasting operations of array & scalar (#7226):
(+){T}(a::KnetArray{T},s::Number)=(.+)(T(s),a)
(+){T}(s::Number,a::KnetArray{T})=(.+)(T(s),a)
(-){T}(a::KnetArray{T},s::Number)=(.+)(T(-s),a)
(-){T}(s::Number,a::KnetArray{T})=(.-)(T(s),a)
(*){T}(a::KnetArray{T},s::Number)=(.*)(T(s),a)
(*){T}(s::Number,a::KnetArray{T})=(.*)(T(s),a)
(/){T}(a::KnetArray{T},s::Number)=(.*)(T(1/s),a)
(\){T}(s::Number,a::KnetArray{T})=(.*)(T(1/s),a)

#(/){T}(s::Number,a::KnetArray{T})=(.*)(T(1/s),a) # not defined in base
#(^){T}(a::KnetArray{T},s::Number)=(.^)(a,T(s)) # linalg
#(^){T}(s::Number,a::KnetArray{T})=(.^)(T(s),a) # linalg


text:
@testset "get_edge_normal" begin
    mmesh = Mesh(2,2,1.0)
    out = get_edge_normal(mmesh)
    bdedge = bcedge(mmesh)
    n1 = mmesh.nodes[bdedge[:,1], :]
    n2 = mmesh.nodes[bdedge[:,2], :]
    @test out≈[ -0.0  -1.0
        0.0   1.0
        0.0   1.0
        1.0   0.0
    -1.0  -0.0
    -1.0  -0.0
    -0.0  -1.0
        1.0   0.0]
end


text:
@static if Sys.isapple()
    using Homebrew
    Homebrew.update()
end


text:
function get_population_stats(population)
    if length(population) > 1
        return maximum(population), minimum(population), mean(population), median(population), var(population)
    else
        val = population[1]
        return val, val, val, val, 0
    end
end
const _single_features_functions_graph = Dict(
    "nb_edge" => (ne, false),
    "nb_vertex" => (nv, false),
    "global_clustering_coefficient" => (global_clustering_coefficient, false),
    "density" => (density, false),
    "diameter" => (diameter, false),
    "radius" => (radius, false),
    "degree" => (degree, true)
)
const _single_features_functions_network = Dict(
    "PD" => (:bus, :PD),
    "QD" => (:bus, :QD),
    "GS" => (:bus, :GS),
    "BS" => (:bus, :BS),
    "VM" => (:bus, :VM),
    "VA" => (:bus, :VA),
    "VMAX" => (:bus, :VMAX),
    "VMIN" => (:bus, :VMIN),
    "PG" => (:gen, :PG),
    "QG" => (:gen, :QG),
    "QMAX" => (:gen, :QMAX),
    "QMIN" => (:gen, :QMIN),
    "PMAX" => (:gen, :PMAX),
    "PMIN" => (:gen, :PMIN),
    "BR_R" => (:branch, :BR_R),
    "BR_X" => (:branch, :BR_X),
    "BR_B" => (:branch, :BR_B)
)
const _single_features_functions_source_path = Dict(
    "nbus" => nbus,
    "nbranch" => nbranch,
    "nbranch_unique" => nbranch_unique,
    "ngen" => ngen_unique
)
const _feature_info_dict = Dict(
    "nb_edge" => (:graph, false),
    "nb_vertex" => (:graph, false),
    "global_clustering_coefficient" => (:graph, false),
    "density" => (:graph, false),
    "diameter" => (:graph, false),
    "radius" => (:graph, false),
    "degree" => (:graph, true),
    "PD" => (:network, true),
    "QD" => (:network, true),
    "GS" => (:network, true),
    "BS" => (:network, true),
    "VM" => (:network, true),
    "VA" => (:network, true),
    "VMAX" => (:network, true),
    "VMIN" => (:network, true),
    "PG" => (:network, true),
    "QG" => (:network, true),
    "QMAX" => (:network, true),
    "QMIN" => (:network, true),
    "PMAX" => (:network, true),
    "PMIN" => (:network, true),
    "BR_R" => (:network, true),
    "BR_X" => (:network, true),
    "BR_B" => (:network, true),
    "nbus" => (:source_path, false),
    "nbranch" => (:source_path, false),
    "nbranch_unique" => (:source_path, false),
    "ngen" => (:source_path, false)
)

# Single features
function get_single_feature_graph(feature_name::AbstractString, g::SimpleGraph)
    feature_tuple = _single_features_functions_graph[feature_name]
    feature = feature_tuple[1](g)
    return !feature_tuple[2] ? feature : get_population_stats(feature)
end
function get_single_feature_network(feature_name::AbstractString, network::PowerFlowNetwork)
    feature_tuple = _single_features_functions_network[feature_name]
    feature = getfield(network, feature_tuple[1])[!, feature_tuple[2]]
    return get_population_stats(feature)
end
function get_single_feature_source_path(feature_name::AbstractString, source_path::AbstractString)
    feature_functions = _single_features_functions_source_path[feature_name]
    feature = feature_functions(source_path)
    return feature
end

# All features
function get_features_graph(g::SimpleGraph)
    degree_max, degree_min, degree_mean, degree_median, degree_var = get_population_stats(degree(g))
    features = Dict(
        "nb_edge" => ne(g),
        "nb_vertex" => nv(g),
        "degree_max" => degree_max,
        "degree_min" => degree_min,
        "degree_mean" => degree_mean,
        "degree_median" => degree_median,
        "degree_var" => degree_var,
        "global_clustering_coefficient" => global_clustering_coefficient(g),
        "density" => density(g),
        "diameter" => diameter(g),
        "radius" => radius(g)
    )
    return features
end
function get_features_opf(network::PowerFlowNetwork)
    # Thanks VIM
    features = Dict()
    # bus
    features["PD_max"], features["PD_min"], features["PD_mean"], features["PD_median"], features["PD_var"] = get_population_stats(network.bus[!, :PD])
    features["QD_max"], features["QD_min"], features["QD_mean"], features["QD_median"], features["QD_var"] = get_population_stats(network.bus[!, :QD])
    features["GS_max"], features["GS_min"], features["GS_mean"], features["GS_median"], features["GS_var"] = get_population_stats(network.bus[!, :GS])
    features["BS_max"], features["BS_min"], features["BS_mean"], features["BS_median"], features["BS_var"] = get_population_stats(network.bus[!, :BS])
    features["VM_max"], features["VM_min"], features["VM_mean"], features["VM_median"], features["VM_var"] = get_population_stats(network.bus[!, :VM])
    features["VA_max"], features["VA_min"], features["VA_mean"], features["VA_median"], features["VA_var"] = get_population_stats(network.bus[!, :VA])
    features["VMAX_max"], features["VMAX_min"], features["VMAX_mean"], features["VMAX_median"], features["VMAX_var"] = get_population_stats(network.bus[!, :VMAX])
    features["VMIN_max"], features["VMIN_min"], features["VMIN_mean"], features["VMIN_median"], features["VMIN_var"] = get_population_stats(network.bus[!, :VMIN])

    # gen
    features["PG_max"], features["PG_min"], features["PG_mean"], features["PG_median"], features["PG_var"] = get_population_stats(network.gen[!, :PG])
    features["QG_max"], features["QG_min"], features["QG_mean"], features["QG_median"], features["QG_var"] = get_population_stats(network.gen[!, :QG])
    features["QMAX_max"], features["QMAX_min"], features["QMAX_mean"], features["QMAX_median"], features["QMAX_var"] = get_population_stats(network.gen[!, :QMAX])
    features["QMIN_max"], features["QMIN_min"], features["QMIN_mean"], features["QMIN_median"], features["QMIN_var"] = get_population_stats(network.gen[!, :QMIN])
    features["PMAX_max"], features["PMAX_min"], features["PMAX_mean"], features["PMAX_median"], features["PMAX_var"] = get_population_stats(network.gen[!, :PMAX])
    features["PMIN_max"], features["PMIN_min"], features["PMIN_mean"], features["PMIN_median"], features["PMIN_var"] = get_population_stats(network.gen[!, :PMIN])

    # branch
    features["BR_R_max"], features["BR_R_min"], features["BR_R_mean"], features["BR_R_median"], features["BR_R_var"] = get_population_stats(network.branch[!, :BR_R])
    features["BR_X_max"], features["BR_X_min"], features["BR_X_mean"], features["BR_X_median"], features["BR_X_var"] = get_population_stats(network.branch[!, :BR_X])
    features["BR_B_max"], features["BR_B_min"], features["BR_B_mean"], features["BR_B_median"], features["BR_B_var"] = get_population_stats(network.branch[!, :BR_B])

    return features
end
function get_features_instance(network::PowerFlowNetwork)
    g = SimpleGraph(network)
    graph_features = get_features_graph(g)
    opf_features = get_features_opf(network)
    features = merge(+, graph_features, opf_features)
    return features
end

function get_clique_features(clique::Vector{Vector{Int}})
    clique_size = [length(c) for c in clique]
    features = Dict()
    features["clique_size_max"], features["clique_size_min"], features["clique_size_mean"], features["clique_size_median"], features["clique_size_var"] = get_population_stats(clique_size)
    return features
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
"""
    simplecycles_hadwick_james(g)

Find circuits (including self-loops) in `g` using the algorithm
of Hadwick & James.

### References
- Hadwick & James, "Enumerating Circuits and Loops in Graphs with Self-Arcs and Multiple-Arcs", 2008
"""
function simplecycles_hadwick_james end
@traitfn function simplecycles_hadwick_james(g::::IsDirected)
    nvg = nv(g)
    T = eltype(g)
    B = Vector{Vector{T}}()
    for i in vertices(g)
        push!(B, Vector{T}())
    end
    blocked = falses(nvg)
    stack = Vector{T}()
    cycles = Vector{Vector{T}}()
    for v in vertices(g)
        circuit_recursive!(g, v, v, blocked, B, stack, cycles)
        resetblocked!(blocked)
        resetB!(B)
    end
    return cycles
end

"""
    resetB!(B)

Reset B work structure.
"""
resetB!(B) = map!(empty!, B, B)

"""
    resetblocked!(blocked)

Reset vector of `blocked` vertices.
"""
resetblocked!(blocked) = fill!(blocked, false)

"""
    circuit_recursive!(g, v1, v2, blocked, B, stack, cycles)

Find circuits in `g` recursively starting from v1.
"""
function circuit_recursive! end
@traitfn function circuit_recursive!{T<:Integer}(g::::IsDirected, v1::T, v2::T, blocked::BitArray, B::Vector{Vector{T}}, stack::Vector{T}, cycles::Vector{Vector{T}})
    f = false
    push!(stack, v2)
    blocked[v2] = true

    Av = out_neighbors(g, v2)
    for w in Av
        (w < v1) && continue
        if w == v1 # Found a circuit
            push!(cycles, copy(stack))
            f = true
        elseif !blocked[w]
            f = circuit_recursive!(g, v1, w, blocked, B, stack, cycles)
        end
    end
    if f
        unblock!(v2, blocked, B)
    else
        for w in Av
            (w < v1) && continue
            if !(v2 in B[w])
                push!(B[w], v2)
            end
        end
    end
    pop!(stack)
    return f
end

"""
    unblock!(v, blocked, B)

Unblock the value `v` from the `blocked` list and remove from `B`.
"""
function unblock!(v::T, blocked::BitArray, B::Vector{Vector{T}}) where T
    blocked[v] = false
    wPos = 1
    Bv = B[v]
    while wPos <= length(Bv)
        w = Bv[wPos]
        wPos += 1 - (length(Bv) - length(filter!(v -> v == w, Bv)))
        if blocked[w]
            unblock!(w, blocked, B)
        end
    end
    return nothing
end


text:
using SafeTestsets, Aqua, RasterDataSources

if VERSION >= v"1.5.0"
    # HTTP.jl `write` is full of ambiguities
    # Aqua.test_ambiguities([RasterDataSources, Base, Core])
    Aqua.test_unbound_args(RasterDataSources)
    Aqua.test_stale_deps(RasterDataSources)
    Aqua.test_undefined_exports(RasterDataSources)
    Aqua.test_project_extras(RasterDataSources)
    Aqua.test_deps_compat(RasterDataSources)
    Aqua.test_project_toml_formatting(RasterDataSources)
end

@time @safetestset "alwb" begin include("alwb.jl") end
@time @safetestset "awap" begin include("awap.jl") end
@time @safetestset "chelsa bioclim" begin include("chelsa-bioclim.jl") end
@time @safetestset "chelsa future" begin include("chelsa-future.jl") end
@time @safetestset "earthenv habitat heterogeneity" begin include("earthenv-heterogeneity.jl") end
@time @safetestset "earthenv landcover" begin include("earthenv-landcover.jl") end
@time @safetestset "worldclim bioclim" begin include("worldclim-bioclim.jl") end
@time @safetestset "worldclim climate" begin include("worldclim-climate.jl") end
@time @safetestset "worldclim weather" begin include("worldclim-weather.jl") end


text:
using GZip

"Where to download mnist from"
mnisturl = "http://yann.lecun.com/exdb/mnist"

"Where to download mnist to"
mnistdir = "../data/"

"""

This utility loads the [MNIST](http://yann.lecun.com/exdb/mnist)
hand-written digits dataset.  There are 60000 training and 10000 test
examples. Each input x consists of 784 pixels representing a 28x28
grayscale image.  The pixel values are converted to Float32 and
normalized to [0,1].  Each output y is a UInt8 indicating the correct
class.  10 is used to represent the digit 0.

```
# Usage:
include(Pkg.dir("Knet/data/mnist.jl"))
xtrn, ytrn, xtst, ytst = mnist()
# xtrn: 28×28×1×60000 Array{Float32,4}
# ytrn: 60000-element Array{UInt8,1}
# xtst: 28×28×1×10000 Array{Float32,4}
# ytst: 10000-element Array{UInt8,1}
```

"""
function mnist()
    global _mnist_xtrn,_mnist_ytrn,_mnist_xtst,_mnist_ytst
    if !@isdefined _mnist_xtrn
        @info "Loading MNIST..."
        _mnist_xtrn = _mnist_xdata("train-images-idx3-ubyte.gz")
        _mnist_xtst = _mnist_xdata("t10k-images-idx3-ubyte.gz")
        _mnist_ytrn = _mnist_ydata("train-labels-idx1-ubyte.gz")
        _mnist_ytst = _mnist_ydata("t10k-labels-idx1-ubyte.gz")
    end
    return _mnist_xtrn,_mnist_ytrn,_mnist_xtst,_mnist_ytst
end

"Utility to view a MNIST image, requires the Images package"
mnistview(x,i)=colorview(Gray,permutedims(x[:,:,1,i],(2,1)))

function _mnist_xdata(file)
    a = _mnist_gzload(file)[17:end]
    reshape(a ./ 255f0, (28,28,1,div(length(a),784)))
end

function _mnist_ydata(file)
    a = _mnist_gzload(file)[9:end]
    a[a.==0] .= 10
    # full(sparse(a,1:length(a),1f0,10,length(a)))
    return a
end

function _mnist_gzload(file)
    if !isdir(mnistdir)
        mkpath(mnistdir)
    end
    path = joinpath(mnistdir,file)
    if !isfile(path)
        url = "$mnisturl/$file"
        download(url, path)
    end
    f = gzopen(path)
    a = read(f)
    close(f)
    return(a)
end

nothing

max_stars_repo_name: davidssmith/KnetLayers

max_issues_repo_name: davidssmith/KnetLayers

max_forks_repo_name: davidssmith/KnetLayers


text:
"""
Reverse max
"""
function max_rev(a::Interval, b::Interval, c::Interval)  # a = max(b,c)

    B_lo = b.lo;  B_hi = b.hi;
    C_lo = c.lo;  C_hi = c.hi;

    (b.lo > c.lo) && (B_lo = max(b.lo,a.lo))
    (b.lo < c.lo) && (C_lo = max(c.lo,a.lo))
    (b.hi > c.hi) && (B_hi = min(b.hi,a.hi))
    (b.hi < c.hi) && (C_hi = min(c.hi,a.hi))

    if isempty(b)
        isempty(c) && (return a, ∅, ∅)
        return a, ∅, Interval(C_lo,C_hi)
    else
        isempty(c) && (return a, Interval(B_lo,B_hi), ∅)
        return a, Interval(B_lo,B_hi), Interval(C_lo,C_hi)
    end
end
max_rev(a,b,c) = max_rev(promote(a,b,c)...)

"""
Reverse min
"""
function min_rev(a::Interval, b::Interval, c::Interval)

    B_lo = b.lo;  B_hi = b.hi;
    C_lo = c.lo;  C_hi = c.hi;

    (b.lo > c.lo) && (B_lo = max(c.lo,a.lo))
    (b.lo < c.lo) && (C_lo = max(b.lo,a.lo))
    (b.hi > c.hi) && (B_hi = min(c.hi,a.hi))
    (b.hi < c.hi) && (C_hi = min(b.hi,a.hi))

    if isempty(b)
        isempty(c) && (return a, ∅, ∅)
        return a, ∅, Interval(C_lo,C_hi)
    else
        isempty(c) && (return a, Interval(B_lo,B_hi), ∅)
        return a, Interval(B_lo,B_hi), Interval(C_lo,C_hi)
    end
end

min_rev(a,b,c) = min_rev(promote(a,b,c)...)


text:
"""
Example illustrating how to update a model with new parameters.
"""

using CSV
using DataStructures
using DataFrames
using Agtor

# Start with the environment variable set for multi-threading
# $ JULIA_NUM_THREADS=4 ./julia

# Windows:
# $ set JULIA_NUM_THREADS=4 && julia dev.jl

data_dir = "test/data/"
zone_spec_dir::String = "$(data_dir)zones/"
zone_params::Dict{String, Dict} = load_spec(zone_spec_dir)

z1 = create(zone_params[:Zone_1])

# Grab some pre-generated parameter combinations
samples = DataFrame!(CSV.File("test/data/scenarios/sampled_params.csv"))

# Below we update the model with the sampled values

# The "value" field for each parameter should change.
@info "before update" z1.fields[1].irrigation.efficiency

# We update the model using values from two rows to illustrate
# that it is changing...
update_model!(z1, samples[1, :])
@info "after update" z1.fields[1].irrigation.efficiency

update_model!(z1, samples[2, :])
@info z1.fields[1].irrigation.efficiency

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
using SortingAlgorithms
import SortingAlgorithms.sortandperm_radix
include("0_setup.jl")

function sumby_sortperm1(by, val)
    (by, p) = sortandperm_radix(by)
    sumby_contiguous(by, val[p])
end

N  = 250_000_000;
K = 100;
by = nothing;
val = nothing;
gc();

# srand(1);
# by = rand(Int32(1):Int32(round(N/K)), N);
# val = rand(Int32(1):Int32(5), N)
# @time res1 = sumby_sortperm1(by, val)

srand(1);
by = rand(Int32(1):Int32(round(N/K)), N);
val = rand(Int32(1):Int32(5), N);
@time res2 = sumby(by, val);

function sumby_sortperm2(by, val)
    (by, p) = sortandperm_radix(by)
    sumby_contiguous(by, @view val[p])
end

srand(1);
by = rand(Int32(1):Int32(round(N/K)), N);
val = rand(Int32(1):Int32(5), N);
@time res3 = sumby_sortperm2(by, val)

srand(1);
by = rand(Int32(1):Int32(round(N/K)), N);
val = rand(Int32(1):Int32(5), N);
@time res4 = sumby_radixsort(by, val)

res2 == res3
res3 == res4