text:
@inline @inbounds function simplestChaoticFlow(X::AbstractArray, 𝐴::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    𝑥̇ = 𝑦
    𝑦̇ = 𝑧
    𝑧̇ = -𝐴(𝑡)*𝑧 + 𝑦^2 - 𝑥
    return SVector{3}(𝑥̇, 𝑦̇, 𝑧̇)
end

@inline @inbounds function simplestChaoticFlow_J(X::AbstractArray, 𝐴::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    J = @SMatrix [  1.0     0.0     0.0;
                    0.0     1.0     0.0;
                    2.0*𝑦   -𝐴(𝑡)   -1.0]
end

"""
    Sprott's simplest chaotic flow

Sprott, or Sprott1997
"""
function simplestChaoticFlow(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=simplestChaoticFlow_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

simplestChaoticFlowSim = Process(
    process = simplestChaoticFlow,
    X0 = [0.05, 0.05, 0.05], # As in Sprott's paper
    parameter_profile = ramp,
    parameter_profile_parameters = (2.02, 2.07, 0.0, 10000.0),
    # Parameters should only be between ~2.018 and ~2.082 (otherwise, unbounded)
    # Range can be different when parameters are time varying
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.5,
    tmax = 10000.0,
    alg = RK4(),
    solver_opts = Dict(:adaptive => false))
export simplestChaoticFlowSim

simplestChaoticFlowVis = Process(
    process = simplestChaoticFlow,
    X0 = [0.05, 0.05, 0.05],
    parameter_profile = ramp,
    parameter_profile_parameters = (2.03, 2.05, 0.0, 5000.0),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.01,
    tmax = 5000.0,
    alg =  RK4(),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-9))
export simplestChaoticFlowVis

polarReduce(x::Real) = x - 2π*(x÷π);
polarReduce(X::AbstractArray) = X .- 2π.*(X.÷π);

@inline @inbounds function doublePendulum(X::AbstractArray, p, t::Real)
    (θ₁, θ₂, 𝑝₁, 𝑝₂) = X
    (ℓ₁, ℓ₂, m₁, m₂) = p(t)
    𝑔 = 0.0

    C₁(θ₁, θ₂, 𝑝₁, 𝑝₂) = (𝑝₁*𝑝₂*sin(θ₁ - θ₂))/(ℓ₁*ℓ₂*(m₁ + m₂*sin(θ₁ - θ₂)^2))
    C₂(θ₁, θ₂, 𝑝₁, 𝑝₂) = ((ℓ₂^2*m₂*𝑝₁^2 + ℓ₁^2*(m₁ + m₂)*𝑝₂^2 - ℓ₁*ℓ₂*m₂*𝑝₁*𝑝₂*cos(θ₁ - θ₂))/(2*ℓ₁^2*ℓ₂^2*(m₁ + m₂*sin(θ₁ - θ₂)^2)^2))*(sin(2*(θ₁ - θ₂)))

    dθ₁ = (ℓ₂*𝑝₁ - ℓ₁*𝑝₂*cos(θ₁ - θ₂))/(ℓ₁^2*ℓ₂*(m₁ + m₂*sin(θ₁ - θ₂)^2))
    dθ₂ = (ℓ₁*(m₁ + m₂)*𝑝₂ - ℓ₂*m₂*𝑝₁*cos(θ₁ - θ₂))/(ℓ₁*ℓ₂^2*m₂*(m₁ + m₂*sin(θ₁ - θ₂)^2))

    d𝑝₁ = -(m₁ + m₂)*𝑔*ℓ₁*sin(θ₁) - C₁(θ₁, θ₂, 𝑝₁, 𝑝₂) + C₂(θ₁, θ₂, 𝑝₁, 𝑝₂)
    d𝑝₂ = -m₂*𝑔*ℓ₂*sin(θ₂) + C₁(θ₁, θ₂, 𝑝₁, 𝑝₂) - C₂(θ₁, θ₂, 𝑝₁, 𝑝₂)

    return SVector{4}(dθ₁, dθ₂, d𝑝₁, d𝑝₂)
end
function doublePendulum(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters))
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

"""
    Your classic double pendulum, with Hamiltonian equations of motion

    https://scienceworld.wolfram.com/physics/DoublePendulum.html
"""
function cartesianDoublePendulum(P::Process)
    # Should the initial conditions be set in cartesian coordinates?
    sol = doublePendulum(P(process=doublePendulum))
    p = parameterseries(P, transient=true)
    x₁ = p[1, :].*sin.(sol[1, :])
    y₁ = .-p[1, :].*cos.(sol[1, :])
    x₂ = x₁ .+ p[2, :].*sin.(sol[2, :])
    y₂ = y₁ .- p[2, :].*cos.(sol[2, :])
    sol = hcat(x₁, y₁, x₂, y₂)
end

doublePendulumSim = Process(
    process = doublePendulum,
    X0 = [π, 0.01, 0.0, 0.0], # Two angles and two momenta
    parameter_profile = (ramp, constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = ((0.0, 1.0, 0.0), (1.0,), (1.0,), (1.0,), (1.0,)), # (threshold, baseline, stepHeight)
    transient_t0 = -10.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.1,
    tmax = 1000.0,
    alg = AutoVern7(Rodas5()))
export doublePendulumSim

cartesianDoublePendulumSim = Process(
    process = cartesianDoublePendulum,
    X0 = [π/2, π/2, 0.0, 0.0], # Two angles and two momenta
    parameter_profile = (constantParameter, constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = ((1.0,), (1.0,), (1.0,), (1.0,)), # (ℓ₁, ℓ₂, m₁, m₂)
    transient_t0 = -10.0,
    t0 = 0.0,
    dt = 0.00001,
    savedt = 0.1,
    tmax = 100.0,
    alg = AutoVern9(Rodas5()),
    solver_opts =  Dict(:adaptive => true, :reltol => 1e-12, :abstol => 1e-12))
export cartesianDoublePendulumSim

cartesianDoublePendulumVis = Process(
    process = cartesianDoublePendulum,
    X0 = [5π/6, π, 0.0, 0.0], # Two angles (from the downward direction) and two momenta
    parameter_profile = (lorentzian, lorentzian, constantParameter, constantParameter),
    transient_t0 = -10.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.001,
    tmax = 50.0,
    parameter_profile_parameters = ((-1.0, 10.0, 25.0, 1.5), (-1.0, 5.0, 25.0, 2.0), (1.0,), (2.0,)),
    alg = AutoVern9(Rodas5()),
    solver_opts = Dict(:adaptive => false))
export cartesianDoublePendulumVis

cartesianDoublePendulumVis2 = Process(
    process = cartesianDoublePendulum,
    X0 = [5π/6, π, 0.0, 0.0], # Two angles (from the downward direction) and two momenta
    parameter_profile = (lorentzian, lorentzian, constantParameter, constantParameter),
    transient_t0 = -10.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.001,
    tmax = 100.0,
    parameter_profile_parameters = ((-1.0, 18.0, 50.0, 1.5), (-1.0, 9.0, 50.0, 2.0), (1.0,), (2.0,)),
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => false))
export cartesianDoublePendulumVis2

@inline @inbounds function waveDrivenHarmonic(X::AbstractArray, p::Function, t::Real)
    (ω, ε, k, Ω) = p(t)
    dX2 = -ω^2.0*X[1] + ε*sin(k*X[1] - Ω*t)
    dX1 = X[2]
    return SVector{2}(dX1, dX2)
end

@inline @inbounds function waveDrivenHarmonic_J(X::AbstractArray, p::Function, t::Real)
    (ω, ε, k, Ω) = p(t)
    J = @SMatrix [0.0 1.0; ε*k*cos(k*X[1] - Ω*t)-ω^2 0.0]
end

"""
    An harmonic oscillator driven by a plane wave propagating in the x-direction

Sprott, or Chernikov1988
"""
function waveDrivenHarmonic(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=waveDrivenHarmonic_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

waveDrivenHarmonicSim = Process(# parameters:  (ω, ε, k, Ω)
    process = waveDrivenHarmonic,
    X0 = [0.1, 0.1],
    parameter_profile = (constantParameter, constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = ((1π,), (5π,), (1.0π,), (3π,)),
    #parameter_profile_parameters = ((1π,), (10.0,), (1.0π,), (0.9π)),
    transient_t0 = -10.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.01,
    tmax = 100.0,
    alg = RK4())
export waveDrivenHarmonicSim

# @inline @inbounds function pulseDrivenHarmonic(X::AbstractArray, p::Function, t::Real)
#     (ω, ε, k, Ω) = p(t)
#     dX2 = -ω^2.0*X[1] + ε*sin(k*X[1] - Ω*t)^2
#     dX1 = X[2]
#     return SVector{2}(dX1, dX2)
# end

# @inline @inbounds function pulseDrivenHarmonic_J(X::AbstractArray, p::Function, t::Real)
#     (ω, ε, k, Ω) = p(t)
#     J = @SMatrix [0.0 1.0; 2*ε*k*cos(k*X[1] - Ω*t)-ω^2 0.0]
# end

# """
#     An harmonic oscillator driven by directional 'pulses'
# """
# function pulseDrivenHarmonic(P::Process)
#     seed(P.solver_rng)
#     prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters))#, jac=waveDrivenHarmonic_J)
#     sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
# end

# pulseDrivenHarmonicSim = Process(# parameters:  (ω, ε, k, Ω)
#     process = pulseDrivenHarmonic,
#     X0 = [0.0, 0.0],
#     parameter_profile = (constantParameter, constantParameter, constantParameter, constantParameter),
#     parameter_profile_parameters = ((1π,), (10.0,), (1π,), (1.5π,)),
#     transient_t0 = -10.0,
#     t0 = 0.0,
#     dt = 0.001,
#     savedt = 0.01,
#     tmax = 100.0,
#     alg = RK4())
# export pulseDrivenHarmonicSim

@inline @inbounds function thomasCyclicallySymmetric(X::AbstractArray, 𝑏::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    𝑥̇ = -𝑏(𝑡)*𝑥 + sin(𝑦)
    𝑦̇ = -𝑏(𝑡)*𝑦 + sin(𝑧)
    𝑧̇ = -𝑏(𝑡)*𝑧 + sin(𝑥)
    return SVector{3}(𝑥̇, 𝑦̇, 𝑧̇)
end

@inline @inbounds function thomasCyclicallySymmetric_J(X::AbstractArray, 𝑏::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    J = @SMatrix [  -𝑏(𝑡)   cos(𝑦)   0.0;
                    0.0     -𝑏(𝑡)    cos(𝑧);
                    cos(𝑥)  0.0      -𝑏(𝑡)]
end

"""
    Thomas' cyclically symmetric attractor

From Sprott, or Thomas1999
"""
function thomasCyclicallySymmetric(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=thomasCyclicallySymmetric_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

thomasCyclicallySymmetricSim = Process(
    process = thomasCyclicallySymmetric,
    X0 = [0.1, 0.0, 0.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (0.18,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 1000.0,
    alg = RK4(),
    solver_opts = Dict(:adaptive => true))
export thomasCyclicallySymmetricSim

thomasCyclicallySymmetricVis = Process(
    process = thomasCyclicallySymmetric,
    X0 = [0.1, 0.0, 0.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (0.18,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.05,
    tmax = 5000.0,
    alg = RadauIIA3(),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-12))
export thomasCyclicallySymmetricVis

@inline @inbounds function doubleScroll(X::AbstractArray, 𝑎::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    𝑥̇ = 𝑦
    𝑦̇ = 𝑧
    𝑧̇ = -𝑎(𝑡)*(𝑧 + 𝑦 + 𝑥 - sign(𝑥))
    return SVector{3}(𝑥̇, 𝑦̇, 𝑧̇)
end
function 𝛿(x)
    if x == 0
        return Inf
    else
        return 0
    end
end
@inline @inbounds function doubleScroll_J(X::AbstractArray, 𝑎::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    J = @SMatrix [  0.0   1.0   0.0;
                    0.0   0.0   1.0;
                    -𝑎(𝑡)*(1.0-2.0*𝛿(x))  -𝑎(𝑡)  -𝑎(𝑡)] # Need to worry about integrating delta? Hopefully we never hit 0.0...
end

"""
    Double scroll attractor

From Sprott, or Elwakil2001
"""
function doubleScroll(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=doubleScroll_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

doubleScrollSim = Process(
    process = doubleScroll,
    X0 = [0.01, 0.01, 0.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (0.8,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 1000.0,
    alg = RK4(),
    solver_opts = Dict(:adaptive => true))
export doubleScrollSim

doubleScrollVis = Process(
    process = doubleScroll,
    X0 = [0.01, 0.01, 0.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (0.8,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.01,
    tmax = 5000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-12))
export doubleScrollVis

@inline @inbounds function lorenz(dX, X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    (𝜎, 𝑟, 𝑏) = p(𝑡)
    dX[1] = 𝜎*(𝑦 - 𝑥)
    dX[2] = -𝑥*𝑧 + 𝑟*𝑥 - 𝑦
    dX[3] = 𝑥*𝑦 - 𝑏*𝑧
end

@inline @inbounds function lorenz_J(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    (𝜎, 𝑟, 𝑏) = p(𝑡)
    J = @SMatrix [  -𝜎      𝜎       0.0;
                    𝑟-𝑧     -1.0    -𝑥;
                    𝑦       𝑥       -𝑏]
end

"""
    Lorenz system

Sprott, Shin1998, Balcerzak2018
"""
function lorenz(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=lorenz_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

lorenzSim = Process(
    process = lorenz,
    X0 = [0.0, -0.01, 9.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = (10.0, 28.0, 8/3), # Sprott's recomendation
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 500.0,
    alg = AutoVern9(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-10, :abstol => 1e-10, :maxiters=>1e7))
export lorenzSim

lorenzVis= Process(
    process = lorenz,
    X0 = [0.0, 0.01, 9.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = (10.0, 28.0, 8/3),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.0005,
    savedt = 0.001,
    tmax = 1000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-15))
export lorenzVis

@inline @inbounds function diffusionlessLorenz(X::AbstractArray, 𝑎::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    𝑥̇ = 𝑦 - 𝑥
    𝑦̇ = -𝑥*𝑧
    𝑧̇ = 𝑥*𝑦 - 𝑎(𝑡)
    return SVector{3}(𝑥̇, 𝑦̇, 𝑧̇)
end
@inline @inbounds function diffusionlessLorenz_J(X::AbstractArray, 𝑎::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    J = @SMatrix [  -1.0   1.0   0.0;
                    -𝑧     0.0   -𝑥;
                    𝑦      𝑥     0.0]
end

"""
    Diffusionless Lorenz system

See Li2014, Sprott2010a (Elegant Chaos), and Schrier2000
"""
function diffusionlessLorenz(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=diffusionlessLorenz_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

diffusionlessLorenzSim = Process(
    process = diffusionlessLorenz,
    X0 = [1.0, 0.0, 1.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (1.0,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 1000.0,
    alg = RK4(),
    solver_opts = Dict(:adaptive => false))
export diffusionlessLorenzSim

diffusionlessLorenzVis = Process(
    process = diffusionlessLorenz,
    X0 = [1.0, 0.0, 1.0],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (1.0,),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.025,
    tmax = 20000.0,
    alg = AutoVern7(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-8))
export diffusionlessLorenzVis

@inline @inbounds function piecewiseLinearHyperchaos(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧, 𝑢) = X
    (𝑎, 𝑏) = p(𝑡)

    𝑥̇ = 𝑦 - 𝑥
    𝑦̇ = -𝑧*sign(𝑥) + 𝑢
    𝑧̇ = abs(𝑥) - 𝑎
    𝑢̇ = -𝑏*𝑦

    return SVector{4}(𝑥̇, 𝑦̇, 𝑧̇, 𝑢̇)
end
@inline @inbounds function piecewiseLinearHyperchaos_J(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧, 𝑢) = X
    (𝑎, 𝑏) = p(𝑡)
    J = @SMatrix [  -1.0        1.0        0.0          0.0;
                    -2*𝛿(𝑧)     0.0        -sign(𝑥)     1.0; # ? Again, hope we miss the zero
                    sign(𝑥)     0.0        0.0          0.0;
                    0.0         -𝑏         0.0          0.0]
end

"""
    A simple, piecewise-linear, hyperchaotic flow

Li2014
"""
function piecewiseLinearHyperchaos(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=piecewiseLinearHyperchaos_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

piecewiseLinearHyperchaosSim = Process(
    process = piecewiseLinearHyperchaos,
    X0 = [0.2, 0.0, 1.0, 0.0],
    parameter_profile = (constantParameter, constantParameter),
    parameter_profile_parameters = (1.0, 0.25), # * 𝑎 is an amplitude, so set to 1 without loss of generality
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 1000.0,
    alg = Vern9(),
    solver_opts = Dict(:adaptive => false, :reltol => 1e-10, :abstol => 1e-10))
export piecewiseLinearHyperchaosSim

piecewiseLinearHyperchaosVis = Process(
    process = piecewiseLinearHyperchaos,
    X0 = [0.01, 0.01, 0.01, 0.01],
    parameter_profile = (constantParameter, constantParameter),
    parameter_profile_parameters = (1.0, 0.25),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 15000.0,
    alg = AutoVern9(Rodas5()),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-14))
export piecewiseLinearHyperchaosVis

@inline @inbounds function simplifiedLorenz4D(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧, 𝑢) = X
    (𝑎, 𝑏) = p(𝑡)

    𝑥̇ = 𝑦 - 𝑥
    𝑦̇ = -𝑥*𝑧 + 𝑢
    𝑧̇ = 𝑥*𝑦 - 𝑎
    𝑢̇ = -𝑏*𝑦

    return SVector{4}(𝑥̇, 𝑦̇, 𝑧̇, 𝑢̇)
end
@inline @inbounds function simplifiedLorenz4D_J(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧, 𝑢) = X
    (𝑎, 𝑏) = p(𝑡)
    J = @SMatrix [  -1.0        1.0        0.0          0.0;
                    -𝑧          0.0        -𝑥           1.0;
                    𝑦           𝑥          0.0          0.0;
                    0.0         -𝑏         0.0          0.0]
end

"""
    The simplified 4D Lorenz system, whose attractors are exclusively hidden and can be hyperchaotic

Li2014a and Gao2006
"""
function simplifiedLorenz4D(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=simplifiedLorenz4D_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

simplifiedLorenz4DSim = Process(
    process = simplifiedLorenz4D,
    X0 = [2.0, 4.0, 0.0, 0.0],
    parameter_profile = (constantParameter, constantParameter),
    parameter_profile_parameters = (2.6, 0.44),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.05,
    tmax = 1000.0,
    alg = RK4(),
    solver_opts = Dict(:adaptive => false))
export simplifiedLorenz4DSim

simplifiedLorenz4DVis = Process(
    process = simplifiedLorenz4D,
    X0 = [2.0, 4.0, 0.0, 0.0],
    parameter_profile = (constantParameter, constantParameter),
    parameter_profile_parameters = (2.6, 0.45),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.05,
    tmax = 7500.0,
    alg = Vern9(),
    solver_opts = Dict(:adaptive => false))#, :reltol => 1e-12))
export simplifiedLorenz4DVis

@inline @inbounds function chensSystem(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    (𝑎, 𝑏, 𝑐) = p(𝑡)

    𝑥̇ = 𝑎*(𝑦 - 𝑥)
    𝑦̇ = (𝑐 - 𝑎)*𝑥 - 𝑥*𝑧 + 𝑐*𝑦
    𝑧̇ = 𝑥*𝑦 - 𝑏*𝑧

    return SVector{3}(𝑥̇, 𝑦̇, 𝑧̇)
end
@inline @inbounds function chensSystem_J(X::AbstractArray, p::Function, 𝑡::Real)
    (𝑥, 𝑦, 𝑧) = X
    (𝑎, 𝑏, 𝑐) = p(𝑡)
    J = @SMatrix [-𝑎        𝑎       0.0;
                  𝑐-𝑎-𝑧     𝑐       -𝑥;
                  𝑦         𝑥       -𝑏]
end

"""
    Yet another chaotic attractor

Chen1999, Song2004
"""
function chensSystem(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=chensSystem_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end
export chensSystem
yetAnotherChaoticAttractor = chensSystem
export yetAnotherChaoticAttractor

chensSystemSim = Process(
    process = chensSystem,
    X0 = [-3.0, 2.0, 20.0],
    parameter_profile = (constant, constant, constant),
    parameter_profile_parameters = ((44,), (5,), (30,)),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.02,
    tmax = 200.0,
    alg = AutoVern9(Rodas5()),
    solver_opts = Dict(:adaptive => true, :abs_tol=>1e-10, :rel_tol=>1e-10))
export chensSystemSim

chensSystemVis = Process(
    process = chensSystem,
    X0 = [-3.0, 2.0, 20.0],
    parameter_profile = (constantParameter, constantParameter, constantParameter),
    parameter_profile_parameters = (35.0, 8.0/3.0, 28.0),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.0001,
    savedt = 0.001,
    tmax = 1000.0,
    alg = Vern9(),
    solver_opts = Dict(:adaptive => true, :reltol => 1e-9))
export chensSystemVis


text:
# This file was generated, do not modify it. # hide
X = transform(fit!(machine(Standardizer(), X)), X);

ext: jl


text:
function merge_nodes!(ds, comp_sizes, i, j)
	iroot = find_root(ds, i)
	jroot = find_root(ds, j)
	if iroot != jroot
		root = root_union!(ds, iroot, jroot)
		comp_sizes[root] = comp_sizes[iroot] + comp_sizes[jroot]
		return root
	else
		return iroot
	end
end

function component_neighbors!(nodes_in_comp, mask, compos, label, present, g, u)
    Q = Queue{Int}()
    label[u] = u
	i = 1
	nodes_in_comp[i] = u
    enqueue!(Q, u)
	ncomp = 0
    while !isempty(Q)
        src = dequeue!(Q)
		i += 1
		nodes_in_comp[i] = src
        for vertex in all_neighbors(g, src)
			if present[vertex]
	            if label[vertex] == 0
	                enqueue!(Q, vertex)
	                label[vertex] = u
	            end
			else
				if !mask[vertex]
					mask[vertex] = true
					ncomp += 1
					compos[ncomp] = vertex
				end
			end
        end
    end
	for j in 1:i
		label[nodes_in_comp[j]] = 0
	end
	for j in 1:ncomp
		mask[compos[j]] = false
	end
    return ncomp
end

function reverse_greedy!(present, g, threshold)
	threshold = max(1, round(Int, threshold))
	n = nv(g)
	ds = IntDisjointSets(n)
	reinsert_nodes = Int[]
	comp_sizes = ones(Int, n) # size of the component where the node in
	S = zeros(Int, n)
	Svec = [Set{Int}() for _ in 1:threshold]
	label = zeros(Int, n)
	nodes_in_comp = zeros(Int, n)
	mask = fill(false, n) # auxiliary vector
	compos = zeros(Int, n) # auxiliary vector
	max_comp_size = 1 # record giant component size

	# add edges between present nodes and record component sizes
	for i in vertices(g)
		if present[i]
			for j in neighbors(g,i)
				if present[j]
					newroot = merge_nodes!(ds, comp_sizes, i, j)
					if comp_sizes[newroot] > max_comp_size
						max_comp_size = comp_sizes[newroot]
					end
				end
			end
		end
	end
	smin = threshold
	for i in vertices(g)
		if !present[i]
			s, ncomp = 1, 0
			for j in neighbors(g,i)
				if present[j]
					jroot = find_root(ds,j)
					if !mask[jroot]
						mask[jroot] = true
						ncomp += 1
						compos[ncomp] = jroot
						s += comp_sizes[jroot]
					end
				end
			end
			# reset the auxiliary vector mask
			for k in 1:ncomp
				mask[compos[k]] = false
			end
			s > threshold && continue
			S[i] = s
			push!(Svec[s], i)
			if s < smin
				smin = s
			end
		end
	end

	while smin <= threshold
		if isempty(Svec[smin])
			smin += 1
			continue
		end
		ibest = rand(Svec[smin])
		for j in neighbors(g, ibest)
			if present[j]
				newroot = merge_nodes!(ds, comp_sizes, ibest, j)
				if comp_sizes[newroot] > max_comp_size
					max_comp_size = comp_sizes[newroot]
				end
			end
		end
		max_comp_size > threshold && break
		delete!(Svec[smin], ibest)
		push!(reinsert_nodes, ibest)
		S[ibest] = 0
		present[ibest] = true

		ncomp = component_neighbors!(nodes_in_comp, mask, compos, label, present, g, ibest)
		for index in 1:ncomp
			i = compos[index]
			if 0 < S[i] <= threshold
				s, ncomp = 1, 0
				for j in neighbors(g,i)
					if present[j]
						jroot = find_root(ds,j)
						if !mask[jroot]
							mask[jroot] = true
							ncomp += 1
							compos[ncomp] = jroot
							s += comp_sizes[jroot]
						end
					end
				end
				# reset the auxiliary vector mask
				for k in 1:ncomp
					mask[compos[k]] = false
				end
				if S[i] != s
					delete!(Svec[S[i]],i)
					s <= threshold && push!(Svec[s], i)
					S[i] = s
				end
			end
		end
	end
	reinsert_nodes
end


text:
#######################################################################
# Structural equality
#######################################################################

# (T, T) → Bool
# Checks `e1` and `e2` for structural equality,
#   i.e. compares all the fields of `e1` and `e2`
# Returns `true` if values are equal
# ASSUMTION: `e1` and `e2` have the same run-time type
# NOTE: Relies on metaprogramming
@generated structEqual(e1, e2) = begin
    # if there are no fields, we can simply return true
    if fieldcount(e1) == 0
        return :(true)
    end
    mkEq    = fldName -> :(e1.$fldName == e2.$fldName)
    # generate individual equality checks
    eqExprs = map(mkEq, fieldnames(e1))
    # construct &&-expression for chaining all checks
    mkAnd   = (expr, acc) -> Expr(:&&, expr, acc)
    # no need in initial accumulator because eqExprs is not empty
    foldr(mkAnd, eqExprs)
end

# (Any, Any) → Bool
# Checks `e1` and `e2` of arbitrary run-time types
# for structural equality and returns `true` if values are equal
genericStructEqual(e1, e2) =
    # if types are different, expressions are not equal
    typeof(e1) != typeof(e2) ?
        false :
    # othewise we need to perform a structural check
        structEqual(e1, e2)

max_stars_repo_name: julbinb/julia-pkgsrc-downloader

max_issues_repo_name: julbinb/julia-pkgsrc-downloader

max_forks_repo_name: julbinb/julia-pkgsrc-downloader


text:
using Clang
using Clang.LibClang.Clang_jll

"""
    generate_ceed_wrapper(ceed_path)

Generate Julia bindings for the libCEED library, where `ceed_path` is the path
of the libCEED directory. The generated bindings are used (with some manual
modifications) to create the low-level C interface LibCEED.C for the LibCEED.jl
package.
"""
function generate_ceed_wrapper(ceed_path)
    headers = ["ceed/ceed.h", "ceed/cuda.h", "ceed/backend.h"]
    ceed_include = joinpath(ceed_path, "include")
    ceed_headers = [joinpath(ceed_include, header) for header in headers]

    wc = init(;
        headers=ceed_headers,
        output_file=joinpath(@__DIR__, "libceed_api_gen.jl"),
        common_file=joinpath(@__DIR__, "libceed_common_gen.jl"),
        clang_includes=[ceed_include, CLANG_INCLUDE],
        clang_args=map(x -> "-I"*x, find_std_headers()),
        header_wrapped=(root, current) -> root == current,
        header_library=x -> "libceed",
        clang_diagnostics=true,
    )
    run(wc, false)
end

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


text:
abstract type AbstractReporter end

struct TextReporter
end

function report(r::TextReporter, rbm::AbstractRBM,
                epoch::Int, epoch_time::Float64, score::Float64)
    println("[Epoch $epoch] Score: $score [$(epoch_time)s]")
end

function report(r::TextReporter, dbn::DBN, epoch::Int, layer::Int)
    println("[Layer $layer] Starting epoch $epoch")
end


text:
import Celeste: BoundingBox
import Celeste.DECALSIO: DECALSDataSet, _get_overlapping_ccds

@testset "decalsio" begin

    # ensure files downloaded
    DATADIR = joinpath(Pkg.dir("Celeste"), "test", "data", "decam")
    wd = pwd()
    cd(DATADIR)
    make_output = readstring(`make survey-ccds-dr5-0.0-0.5-0.0-0.5.fits`)
    if !startswith(make_output, "make: Nothing to be done for")
        print(make_output)  # only print output if something was done
    end
    cd(wd)

    @testset "overlapping ccds" begin
        dataset = DECALSDataSet(DATADIR,
                                "survey-ccds-dr5-0.0-0.5-0.0-0.5.fits")
        box = BoundingBox(0.0, 0.5, 0.0, 0.5)
        idx = _get_overlapping_ccds(dataset, box)
        @test length(idx) == 134  # all CCDs in FITS file should overlap.

        small_box = BoundingBox(0.2, 0.4, 0.2, 0.4)
        idx2 = _get_overlapping_ccds(dataset, small_box)
        @test length(idx2) == 82  # determined from running code!
                                  # (But should be less than 134)
    end
end


text:
# ****************************************************************
#                      1D heat conduction
# ****************************************************************

"""
- `A_B`   : Bond albedo
- `A_TH`  : Albedo at thermal radiation wavelength
- `k`     : Thermal conductivity
- `ρ`     : Density [kg/m³]
- `Cₚ`    : Heat capacity [J/kg/K]
- `ϵ`     : Emissivity
- `P`     : Rotation period [s]
- `l`     : Thermal skin depth
- `Γ`     : Thermal inertia
- `Δt`    : Time step
- `t_bgn` : Start time of the simulation
- `t_end` : End time of the simulation
- `Nt`    : Number of time step
- `Δz`    : Depth step
- `z_max` : Maximum depth for themal simualtion
- `Nz`    : Number of depth step
- `λ`     : Non-dimensional coefficient for heat diffusion equation
"""
struct ParamsThermo{T}
    A_B::T
    A_TH::T
    k::T
    ρ::T
    Cₚ::T
    ϵ::T
    P::T
    l::T
    Γ::T
    Δt::T
    t_bgn::T
    t_end::T
    Nt::Int64
    Δz::T
    z_max::T
    Nz::Int64
    λ::T
end

function ParamsThermo(; A_B, A_TH, k, ρ, Cₚ, ϵ, P, Δt, t_bgn, t_end, Δz, z_max)
    l = thermal_skin_depth(P, k, ρ, Cₚ)
    Γ = thermal_inertia(k, ρ, Cₚ)

    Δt /= P
    t_bgn /= P
    t_end /= P
    Nt = length(t_bgn:Δt:t_end)

    Δz /= l
    z_max /= l
    Nz = length(0:Δz:z_max)

    λ = 1/4π * (Δt/Δz^2)
    λ > 0.5 && println("λ should be smaller than 0.5 for convergence.")

    ParamsThermo(A_B, A_TH, k, ρ, Cₚ, ϵ, P, l, Γ, Δt, t_bgn, t_end, Nt, Δz, z_max, Nz, λ)
end

function Base.show(io::IO, params_thermo::ParamsThermo)
    @unpack A_B, A_TH, k, ρ, Cₚ, ϵ, P, l, Γ, Δt, t_bgn, t_end, Nt, Δz, z_max, Nz, λ = params_thermo

    println(io, "Thermophysical parameters")
    println("-------------------------")

    println("A_B   : ", A_B)
    println("A_TH  : ", A_TH)
    println("k     : ", k)
    println("ρ     : ", ρ)
    println("Cₚ    : ", Cₚ)
    println("ϵ     : ", ϵ)
    println("P     : ", P)
    println("l     : ", l)
    println("Γ     : ", Γ)
    println("Δt    : ", Δt)
    println("t_bgn : ", t_bgn)
    println("t_end : ", t_end)
    println("Nt    : ", Nt)
    println("Δz    : ", Δz)
    println("z_max : ", z_max)
    println("Nz    : ", Nz)
    println("λ     : ", λ)
end

"""
    thermal_skin_depth(P, k, ρ, Cₚ) -> l_2π

# Arguments
- `P`  :
- `k`  :
- `ρ`  :
- `Cₚ` :

# Return
`l_2π` : Thermal skin depth
"""
thermal_skin_depth(P, k, ρ, Cₚ) = √(4π * P * k / (ρ * Cₚ))
thermal_skin_depth(params) = thermal_skin_depth(params.P, params.k, params.ρ, params.Cₚ)

"""
    thermal_inertia(k, ρ, Cₚ) -> Γ

# Arguments
- `k`  :
- `ρ`  :
- `Cₚ` :

# Return
`Γ` : Thermal inertia
"""
thermal_inertia(k, ρ, Cₚ) = √(k * ρ * Cₚ)
thermal_inertia(params) = thermal_inertia(params.k, params.ρ, params.Cₚ)

"""
    update_temperature!(Tⱼ, Tⱼ₊₁, F, params_thermo)

Update temerature profie based on 1-D heat diffusion

# Arguments
- `Tⱼ`            : Temperatures
- `Tⱼ₊₁`          : Temperatures at the next time step
- `F`             : Energy flux to the surface
- `params_thermo` : Thermophysical parameters

i : index of depth
j : index of time step

for i in 2:length(Tⱼ)-1
    Tⱼ₊₁[i] = (1-2λ)*Tⱼ[i] + λ*(Tⱼ[i+1] + Tⱼ[i-1])
end
"""
function update_temperature!(Tⱼ, Tⱼ₊₁, F, params_thermo)
    @unpack λ = params_thermo

    @. Tⱼ₊₁[begin+1:end-1] = @views (1-2λ)*Tⱼ[begin+1:end-1] + λ*(Tⱼ[begin+2:end] + Tⱼ[begin:end-2])

    update_surface_temperature!(Tⱼ₊₁, F, params_thermo)  # Sourface boundary condition (Radiation)
    Tⱼ₊₁[end] = Tⱼ₊₁[end-1]                              # Internal boundary condition (Insulation)

    Tⱼ .= Tⱼ₊₁
end

function update_temperature!(shape, params_thermo)
    @unpack A_B, A_TH = params_thermo

    for facet in shape.facets
        F_sun  = facet.flux.sun
        F_scat = facet.flux.scat
        F_rad  = facet.flux.rad

        F_total = (1 - A_B)*(F_sun + F_scat) + (1 - A_TH)*F_rad

        update_temperature!(facet.Tz, shape.Tz⁺, F_total, params_thermo)
    end
end

"""
    update_surface_temperature!(T, F, params_thermo)

Solve Newton's method to get surface temperature
"""
function update_surface_temperature!(T, F, params_thermo)
    @unpack Δz, Γ, P, ϵ = params_thermo

    for _ in 1:20
        T_pri = T[begin]

        f = F + Γ / √(4π * P) * (T[begin+1] - T[begin]) / Δz - ϵ*σ_SB*T[begin]^4
        df = - Γ / √(4π * P) / Δz - 4*ϵ*σ_SB*T[begin]^3
        T[begin] -= f / df

        err = abs(1 - T_pri / T[begin])
        err < 1e-10 && return
    end
end

# ****************************************************************
#
# ****************************************************************

"""
    intensity(λ, T) -> I

Intensity of radiation at a wavelength λ and tempertature T
according to the Planck function
"""
function intensity(λ, T)
    h = 6.62607015e-34  # Planck constant [J⋅s]
    k = 1.380649e-23    # Boltzmann's constant [J/K]

    I = 2 * h * c^2 / λ^5 / (exp(h * c₀ / (λ * k * T)) - 1)
end

ν2λ(ν) = c₀ / ν
λ2ν(λ) = c₀ / λ

max_stars_repo_name: MasanoriKanamaru/Astroshaper

max_issues_repo_name: MasanoriKanamaru/Astroshaper

max_forks_repo_name: MasanoriKanamaru/Astroshaper


text:
function Green(xr,xz,yr,yz; derivative=0)
    #=
    Fundamental solution (Green function) and derivatives of the Grad-Shafranov Delta^* operator
    -Delta_x^* G(x,y)=delta(x-y)
    y=(yr,yz) is the source point (a filament with unit current)
    x=(xr,xz) is the evaluation point
    Derivatives from Itagaki, Engineering Analysis with Boundary Elements, 33 (2009) 845-857

    derivative=0 only Green function value
    derivative=1 + gradient wrt to x and y
    derivative=2 + second order cross derivatives (see return outputs below)
    =#

    mu0=4*pi*10^-7
    #modulus at points x and y
    k_xy=sqrt(4*xr*yr/((xr+yr)^2+(xz-yz)^2))

    #Complete elliptic integrals of first and second kind
    K_xy,E_xy=Elliptic.ellipke(k_xy^2)

    #Green function
    G_xy = (mu0/pi*((xr*yr)^(1/2))/k_xy)*((1-0.5*k_xy^2)*K_xy - E_xy)

    #println("derivative=",derivative)

    if(derivative==0)
        return (G_xy=G_xy,)
        #return G_xy
    end

    if(derivative>0)
        # gradient with respect to y
        dyrG_xy = (mu0/(2*pi))*((yr   )/(sqrt( (xr+yr)^2 + (xz-yz)^2 )))*(K_xy - ((yr^2-xr^2+(xz-yz)^2 )/( (xr-yr)^2 + (xz-yz)^2))*E_xy)
        dyzG_xy = (mu0/(2*pi))*((yz-xz)/(sqrt( (xr+yr)^2 + (xz-yz)^2 )))*(K_xy - ((yr^2+xr^2+(xz-yz)^2 )/( (xr-yr)^2 + (xz-yz)^2))*E_xy)

        # gradient with respect to x
        dxrG_xy = (mu0/(2*pi))*((xr   )/(sqrt( (xr+yr)^2 + (xz-yz)^2 )))*(K_xy - ((xr^2-yr^2+(xz-yz)^2 )/( (xr-yr)^2 + (xz-yz)^2))*E_xy)
        dxzG_xy = (mu0/(2*pi))*((xz-yz)/(sqrt( (xr+yr)^2 + (xz-yz)^2 )))*(K_xy - ((xr^2+yr^2+(xz-yz)^2 )/( (xr-yr)^2 + (xz-yz)^2))*E_xy)

        if(derivative==1)
            #return G_xy, dyrG_xy, dyzG_xy, dxrG_xy, dxzG_xy
            return (G_xy=G_xy, dyrG_xy=dyrG_xy, dyzG_xy=dyzG_xy, dxrG_xy=dxrG_xy, dxzG_xy=dxzG_xy)
        end

    end

    if(derivative==2)

        # second order cross derivatives
        dxrdyrG_xy = -((xr+yr)/((xr+yr)^2+(xz-yz)^2))*dxrG_xy +
            (mu0/2*pi)*xr*(1/((xr+yr)^2+(xz-yz)^2)^(1/2))*(1/((xr-yr)^2+(xz-yz)^2))*
            ( ((xr^2-yr^2+(xz-yz)^2)/((xr+yr)^2+(xz-yz)^2))*((xr-yr)*K_xy+(xr+yr)*E_xy) -2*xr*E_xy + ((4*yr*(xz-yz)^2)/((xr-yr)^2+(xz-yz)^2))*E_xy )

        dxzdyrG_xy = -((xr+yr)/((xr+yr)^2+(xz-yz)^2))*dxzG_xy +
            (mu0/2*pi)*xr*(xz-yz)*(1/((xr-yr)^2+(xz-yz)^2)^(1/2))*(1/((xr-yr)^2+(xz-yz)^2))*
            ( ((xr^2-yr^2+(xz-yz)^2)/((xr+yr)^2+(xz-yz)^2))*(K_xy+E_xy) -2*E_xy - ((4*yr*(xr-yr))/((xr-yr)^2+(xz-yz)^2))*E_xy )

        dxrdyzG_xy = ((xz-yz)/((xr+yr)^2+(xz-yz)^2))*dxrG_xy +
            (mu0/pi)*yr*xr*(xz-yz)*
            (1/((xr+yr)^2+(xz-yz)^2)^(1/2))*(1/((xr-yr)^2+(xz-yz)^2))*
            ( (1/((xr+yr)^2+(xz-yz)^2))*((xr-yr)*K_xy+(xr+yr)*E_xy) +
            ((-2*(xr-yr))/((xr-yr)^2+(xz-yz)^2))*E_xy )

        dxzdyzG_xy = (((-1)/(xz-yz))+((xz-yz)/((xr+yr)^2+(xz-yz)^2)))*dxzG_xy +
            (mu0/pi)*yr*xr*(xz-yz)^2 *
            (1/((xr+yr)^2+(xz-yz)^2)^(1/2))*(1/((xr-yr)^2+(xz-yz)^2))*
            ( (1/((xr+yr)^2+(xz-yz)^2))*(K_xy+E_xy) +
            (-2/((xr-yr)^2+(xz-yz)^2))*E_xy )

        #return G_xy, dyrG_xy, dyzG_xy, dxrG_xy, dxzG_xy, dxrdyrG_xy, dxzdyrG_xy, dxrdyzG_xy, dxzdyzG_xy
        return (G_xy=G_xy, dyrG_xy=dyrG_xy, dyzG_xy=dyzG_xy, dxrG_xy=dxrG_xy, dxzG_xy=dxzG_xy, dxrdyrG_xy=dxrdyrG_xy, dxzdyrG_xy=dxzdyrG_xy, dxrdyzG_xy=dxrdyzG_xy, dxzdyzG_xy=dxzdyzG_xy)
    end

end


text:
# Forming adjacency graphs.

# Read CSV file into associative array.
E = ReadCSV(joinpath(Base.source_dir(),"Edge.csv"))

# Get vertices and convert to numbers.
Ev = logical( E[:, StartsWith("V,")] )

# Compute vertex adjacency graph.
Av = sqIn(Ev)
printFull(Av)

# Compute edge adjacency graph.
Ae = sqOut(Ev)
printFull(Ae)


text:
function n_queens(n, ::Val{:JuMP})
    model = JuMP.Model(CBLS.Optimizer)

    @variable(model, 1 ≤ Q[1:n] ≤ n, Int)
    @constraint(model, Q in AllDifferent())

    for i in 1:n, j in i + 1:n
        @constraint(model, [Q[i],Q[j]] in Predicate(x -> x[1] != x[2]))
        @constraint(model, [Q[i],Q[j]] in Predicate(x -> x[1] != x[2] + i - j))
        @constraint(model, [Q[i],Q[j]] in Predicate(x -> x[1] != x[2] + j - i))
    end

    return model, Q
end

"""
    n_queens(n; modeler = :JuMP)

Create a model for the n-queens problem with `n` queens. The `modeler` argument accepts :JuMP (default), which refer to the JuMP model.
"""
n_queens(n; modeler=:JuMP) = n_queens(n, Val(modeler))


text:
""" geometric predicates for 2vectors, cartesian points (px,py), homogeneuos points (pw,px,py), lines (lw,lx,ly), and segments x[1:2],y[1:2].
	Julia can inline these.
"""
Vector2 = Tuple{Real,Real}
Point2C = Tuple{Real,Real}
Line2 = Tuple{Real,Real,Real}
Point2H = Tuple{Real,Real,Real}
Point2 = Union{Point2C,Point2H}

toPoint2H((x,y)::Point2C)::Point2H = (one(x),x,y)

# some special handling for floats
dot((ux,uy)::Tuple{AbstractFloat,AbstractFloat}, (vx,vy)::Tuple{AbstractFloat,AbstractFloat}) = muladd(ux, vx, 1.0uy*vy)
dot((ux,uy)::Vector2, (vx,vy)::Vector2) = muladd(ux, vx, uy*vy) # dot product for 2d vectors. (consider muladd)
perpdot((ux,uy)::Tuple{AbstractFloat,AbstractFloat}, (vx,vy)::Tuple{AbstractFloat,AbstractFloat}) = muladd(ux, vy, -1.0uy*vx)
perpdot((ux,uy)::Vector2, (vx,vy)::Vector2) = muladd(ux, vy, -uy*vx) # dot product of ccw perp(u) with v. (consider muladd)

# basic cartesian and homogeneous orientation determinants
ccw(p::Point2C,q::Point2C,r::Point2C) = perpdot(q.-p,r.-p) # orientation determinant
ccw((pw,px,py)::Point2H,(qw,qx,qy)::Point2H,(rw,rx,ry)::Point2H) = pw*perpdot((qx,qy),(rx,ry)) -qw*perpdot((px,py),(rx,ry)) + rw*perpdot((px,py),(qx,qy))# orientation determinant
ccw(p::Point2C,q::Point2C,(rw,rx,ry)::Point2H) = perpdot(q.-p,(rx,ry))+rw*perpdot(p,q.-p) # orientation determinant

# mixed cartesian and homogeneuos
ccw(p::Point2C,q::Point2H,r::Point2H) = ccw(toPoint2H(p),q,r)
ccw(p::Point2H,q::Point2C,r::Point2H) = ccw(p,toPoint2H(q),r)
ccw(p::Point2H,q::Point2H,r::Point2C) = ccw(p,q,toPoint2H(r))
ccw(p::Point2H,q::Point2C,r::Point2C) = ccw(q,r,p)
ccw(p::Point2C,q::Point2H,r::Point2C) = ccw(r,p,q)

## Geometric predicates based on ccw
leftturn(p::Point2,q::Point2,r::Point2)::Bool = ccw(p,q,r)>zero(p[1]) # p to q to r is a left turn
rightturn(p::Point2,q::Point2,r::Point2)::Bool = ccw(p,q,r)<zero(p[1]) # p to q to r is a right turn

## leftAhead halps solve on-the-llne cases consistent with infinitesimal rotation
leftAhead(p::Point2C,q::Point2C,r::Point2C)::Bool = (d = ccw(p,q,r)) == zero(d) ? dot(p.-q,r.-q)<zero(d) : d>zero(d) # should assert ≉ 0 for precision checking
rightAhead(p::Point2C,q::Point2C,r::Point2C)::Bool = (d = ccw(p,q,r)) == zero(d) ? dot(p.-q,r.-q)<zero(d) : d<zero(d) # should assert ≉ 0 for precision checking
leftOnAhead(p::Point2C,q::Point2C,r::Point2C)::Bool = (d = ccw(p,q,r)) == zero(d) ? dot(q.-p,r.-p)>zero(d) : d>zero(d) # should assert ≉ 0 for precision checking
rightOnAhead(p::Point2C,q::Point2C,r::Point2C)::Bool = (d = ccw(p,q,r)) == zero(d) ? dot(q.-p,r.-p)>zero(d) : d<zero(d) # should assert ≉ 0 for precision checking

## orientation with drected segments
#  This should probably be expanded.
leftofseg(q::Point2, xs, ys)::Bool = @inbounds leftturn((xs[1],ys[1]),(xs[2],ys[2]), q) # point q left of segment defined by pairs of xs and ys
rightofseg(q::Point2, xs, ys)::Bool = @inbounds rightturn((xs[1],ys[1]),(xs[2],ys[2]), q) # point q right of segment defined by pairs of xs and ys

leftward(s0::Point2C,s1::Point2C, t0::Point2C,t1::Point2C)::Bool = perpdot(t1.-t0, s1.-s0)>zero(s0[1])   # vector s0->s1 is in the direction from right to left of the line t0->t1.
rightward(s0::Point2C,s1::Point2C, t0::Point2C,t1::Point2C)::Bool = perpdot(t1.-t0, s1.-s0)<zero(s0[1])   # vector s0->s1 is in the direction from left to right of the line t0->t1.

# computing homogeneous line equations
line(p::Point2C,q::Point2C)::Line2 = (perpdot(p,q), last(p)-last(q), first(q)-first(p)) # line equation
orientptln((x,y)::Point2C, (lw,lx,ly)::Line2) = muladd(y,ly,muladd(x,lx,lw))
orientptln((w,x,y)::Point2H, (lw,lx,ly)::Line2) = muladd(y,ly,muladd(x,lx,w*lw))
leftof(q::Point2, l::Line2)::Bool = orientptln(q,l)>zero(l[1]) # point left of (consistent with) oriented line
rightof(q::Point2, l::Line2)::Bool = orientptln(q,l)<zero(l[1])

max_stars_repo_name: snoeyink/CGinJulia

max_issues_repo_name: snoeyink/CGinJulia

max_forks_repo_name: snoeyink/CGinJulia


text:
using LinearAlgebra

function RCCSDpT(rhf::RHF, Alg::RpTa)
    aoints = IntegralHelper()
    moints = IntegralHelper(orbitals=rhf.orbitals)
    ccsd = RCCSD(moints, aoints)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(Alg::RpTa)
    aoints = IntegralHelper{Float64}()
    rhf = RHF(aoints)
    moints = IntegralHelper(orbitals=rhf.orbitals)
    ccsd = RCCSD(moints, aoints)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(ccsd::RCCSD, moints::IntegralHelper{T,E,O}, Alg::RpTa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    output("\n   • Perturbative Triples Started\n")

    T1 = ccsd.T1.data
    T2 = ccsd.T2.data

    Vovvv = moints["OVVV"].data
    Vooov = moints["OOOV"].data
    Vovov = moints["OVOV"].data

    o,v = size(T1)

    fo = moints["Fii"].data
    fv = moints["Faa"].data

    # Pre-allocate Intermediate arrays
    Ws  = [Array{T}(undef, v,v,v) for i = 1:Threads.nthreads()]
    Vs  = [Array{T}(undef, v,v,v) for i = 1:Threads.nthreads()]
    #Evals = zeros(T, Threads.nthreads())
    Evals = zeros(T, o)

    output("Computing energy contribution from occupied orbitals:")
    LinearAlgebra.BLAS.set_num_threads(1)
    t = @elapsed begin
    Threads.@threads for i in 1:o
        id = Threads.threadid()
        W = Ws[id]
        V = Vs[id]
        @views Vovvv_1i = Vovvv[i,:,:,:]

        @views T2_1i    = T2[i,:,:,:]
        @views T1_1i    = T1[i, :]
        for j in 1:i
            Dij = fo[i] + fo[j]

            @views Vovvv_1j    = Vovvv[j,:,:,:]
            @views Vovov_1i_3j = Vovov[i,:,j,:]

            @views Vooov_2i_3j = Vooov[:,i,j,:]
            @views Vooov_2j_3i = Vooov[:,j,i,:]

            @views T2_1i_2j    = T2[i,j,:,:]
            @views T2_1j_2i    = T2[j,i,:,:]
            @views T2_1j       = T2[j,:,:,:]
            @views T1_1j       = T1[j, :]
            δij = i == j
            for k in 1:j
                Dijk = Dij + fo[k]

                @views Vovvv_1k    = Vovvv[k,:,:,:]

                @views Vooov_2j_3k = Vooov[:,j,k,:]
                @views Vooov_2k_3j = Vooov[:,k,j,:]
                @views Vooov_2k_3i = Vooov[:,k,i,:]
                @views Vooov_2i_3k = Vooov[:,i,k,:]

                @views T2_1k_2j = T2[k,j,:,:]
                @views T2_1k_2i = T2[k,i,:,:]
                @views T2_1i_2k = T2[i,k,:,:]
                @views T2_1j_2k = T2[j,k,:,:]

                @views T2_1k = T2[k,:,:,:]
                @views T1_1k = T1[k, :]

                @views Vovov_1j_3k = Vovov[j,:,k,:]
                @views Vovov_1i_3k = Vovov[i,:,k,:]

                # Build intermediates W and V for the set i j k
                @tensoropt begin
                    W[a,b,c]  =  Vovvv_1i[a,b,d]*T2_1k_2j[c,d] - Vooov_2j_3k[l,c]*T2_1i[l,a,b]  # ijk abc
                    W[a,b,c] +=  Vovvv_1i[a,c,d]*T2_1j_2k[b,d] - Vooov_2k_3j[l,b]*T2_1i[l,a,c]  # ikj acb
                    W[a,b,c] +=  Vovvv_1k[c,a,d]*T2_1j_2i[b,d] - Vooov_2i_3j[l,b]*T2_1k[l,c,a]  # kij cab
                    W[a,b,c] +=  Vovvv_1k[c,b,d]*T2_1i_2j[a,d] - Vooov_2j_3i[l,a]*T2_1k[l,c,b]  # kji cba
                    W[a,b,c] +=  Vovvv_1j[b,c,d]*T2_1i_2k[a,d] - Vooov_2k_3i[l,a]*T2_1j[l,b,c]  # jki bca
                    W[a,b,c] +=  Vovvv_1j[b,a,d]*T2_1k_2i[c,d] - Vooov_2i_3k[l,c]*T2_1j[l,b,a]  # jik bac

                    V[a,b,c]  = W[a,b,c] + Vovov_1j_3k[b,c]*T1_1i[a] + Vovov_1i_3k[a,c]*T1_1j[b] + Vovov_1i_3j[a,b]*T1_1k[c]
                end

                # Compute Energy contribution
                δjk = j == k
                for a in 1:v
                    Dijka = Dijk - fv[a]
                    for b in 1:a
                        Dijkab = Dijka - fv[b]
                        δab = a == b
                        for c in 1:b
                            Dd = Dijkab - fv[c]
                            δbc = b == c
                            @inbounds begin
                                X = (W[a,b,c]*V[a,b,c] + W[a,c,b]*V[a,c,b] + W[b,a,c]*V[b,a,c] + W[b,c,a]*V[b,c,a] + W[c,a,b]*V[c,a,b] + W[c,b,a]*V[c,b,a])
                                Y = (V[a,b,c] + V[b,c,a] + V[c,a,b])
                                Z = (V[a,c,b] + V[b,a,c] + V[c,b,a])
                                Ef = (Y - 2*Z)*(W[a,b,c] + W[b,c,a] + W[c,a,b]) + (Z - 2*Y)*(W[a,c,b]+W[b,a,c]+W[c,b,a]) + 3*X
                            end
                            Evals[i] += Ef*(2-δij-δjk)/(Dd*(1+δab+δbc))
                        end
                    end
                end
            end
        end
        output("  Orbital {} ✔", i)
    end
    end

    Et = sum(Evals)
    output("Finished in {:5.5f} s", t)
    output("Final (T) contribution: {:15.10f}", Et)
    output("CCSD(T) energy:         {:15.10f}", Et+ccsd.energy)
    return RCCSDpT{T}(ccsd, Et+ccsd.energy, Et)
end


text:
"""
    gershgorin_disc(i, row)

Return the [Gershgorin
disc](https://en.wikipedia.org/wiki/Gershgorin_circle_theorem) for the
`i`th `row` of a matrix.
"""
function gershgorin_disc(i, row::AbstractVector)
    aᵢᵢ = row[i]
    aᵢᵢ, (sum(abs, row) - abs(aᵢᵢ))
end

"""
    gershgorin_discs(A)

Compute all [Gershgorin
discs](https://en.wikipedia.org/wiki/Gershgorin_circle_theorem) of a
matrix `A`.
"""
gershgorin_discs(A::AbstractMatrix) =
    map(((i,row),) -> gershgorin_disc(i,row), enumerate(eachrow(A)))

"""
    gershgorin_bounds(A)

For a Hermitian matrix `A`, compute the bounds on the real axis of the
eigenspectrum. NB that this bound will only be useful for diagonally
dominant `A`s.
"""
function gershgorin_bounds(A::AbstractMatrix)
    ishermitian(A) ||
        throw(ArgumentError("Not implemented for non-Hermitian matrices"))
    bounds = Inf..(-Inf)
    for (i,row) in enumerate(eachrow(A))
        a,r = gershgorin_disc(i,row)
        bounds = bounds ∪ ((a-r)..(a+r))
    end
    bounds
end


text:
include("../src/DelaunayTessellationFieldEstimator.jl")
using DelaunayTessellationFieldEstimator
using Plots

Coords = rand(0.0:500.0,10,2)

tess = DelaunayTessellation()

C = Coordinates(Coords)

P = BuildDelaunay!(tess,C.GeometryCoordinates)
x,y = VoronoiDelaunay.getplotxy(delaunayedges(tess))
plot(x,y)
scatter!(C.GeometryCoordinates[:,1],C.GeometryCoordinates[:,2])

t = FindTriangles(tess,P[1])
A = ContiguousVoronoiCellArea(t)

# yet to get a triangle counter
(A*GeomToDataAreaScalingFactor(C))

u=UniquePoints(t)
DelaunayVertexDensity(tess,P[1])

length(u)/(A*GeomToDataAreaScalingFactor(C))

DelaunayVertexDensity(tess,C)


text:
using GeometryBasics
export visualize

"""
    visualize(pcloud::PointCloud; kwargs...)

Visualize PointCloud `pcloud` at `index`.

Dimension of points in PointCloud `pcloud` must be 3.

### Optional Arguments:
- color (Symbol)       - Color of the marker, default `:blue`
- markersize (Number)  - Size of the marker, default `npoints(pcloud)/5000`

"""
function visualize(p::PointCloud, index::Number = 1; kwargs...)
    points = cpu(p[index])
    size(points, 1) == 3 || error("dimension of points in PointCloud must be 3.")

    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :lightgreen)
    get!(kwargs, :markersize, npoints(p) / 5000)

    meshscatter(points[1, :], points[2, :], points[3, :]; kwargs...)
end

"""
    visualize(m::TriMesh, index::Int=1; kwargs...)

Visualize mesh at `index` in TriMesh `m`.

### Optional Arguments:
- color (Symbol)       - Color of the marker, default `:red`

"""
function visualize(m::GeometryBasics.Mesh; kwargs...) where {T,R}
    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :orange)

    mesh(GeometryBasics.normal_mesh(m); kwargs...)
end

visualize(m::TriMesh, index::Int = 1; kwargs...) = visualize(GBMesh(m, index); kwargs...)

"""
    visualize(v::VoxelGrid, index::Int=1; kwargs...)

Visualize voxel at `index` in VoxelGrid `v`.

### Optional Arguments:
- color (Symbol)       - Color of the marker, default `:red`

"""
function visualize(
    v::VoxelGrid,
    index::Int = 1,
    thresh::Number = 0.49f0;
    algo = :Exact,
    kwargs...,
)
    algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets] ||
        error("given algo: $(algo) is not supported. Accepted algo are
              {:Exact,:MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets}.")
    kwargs = convert(Dict{Symbol,Any}, kwargs)
    get!(kwargs, :color, :violet)
    method = algo == :Exact ? _voxel_exact : _voxel_algo
    v, f = method(cpu(v[index]), Float32(thresh), algo)

    m = GBMesh(v, f)
    mesh(GeometryBasics.normal_mesh(m); kwargs...)
end

visualize(v::Dataset.AbstractDataPoint; kwargs...) = visualize(v.data; kwargs...)

visualize(v::AbstractCustomObject; kwargs...) =
    error("Define visualize function for custom type: $(typeof(v)).
            Use `import Flux3D.visualize` and define function
            `visualize(v::$(typeof(v)); kwargs...)`")


text:
using Revise
using DataFrames, CSV, Statistics, Distributions, Gadfly, Cairo, Fontconfig

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end
using CTCopy

const PLOT_NAME = "plot"
const TAX_FIELD = :interesttax
const CORP_FIELD = :corplt
const MUNI_FIELD = :muni20
const TREASURY_FIELD = :tnote20
const N_SIMS = 10_000

const LHS = [:creditspread, :taxtnote20, TAX_FIELD]
const RHS = (s->Symbol(s,"L1")).(LHS)
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}}

#gets the file
readPlot(; plotName=PLOT_NAME)::DataFrame = CSV.read("$plotName.csv")

#make the intermediate variables
function getvar!(df; taxfield = TAX_FIELD, corpfield=CORP_FIELD, munifield = MUNI_FIELD,
  treasuryfield = TREASURY_FIELD)::VAR

  df[:creditspread] = df[treasuryfield] .- df[corpfield]
  df[:munispread] = df[treasuryfield] .- df[munifield]
  df[:taxtnote20] = df[treasuryfield] .* df[taxfield]
  df[:resid1t] = df[:munispread] .- (df[:creditspread] .* (1.0 .- df[taxfield]) .+ df[:taxtnote20])
  df[:resid1] = df[:munispread] .- (df[:creditspread] .+ df[:taxtnote20])

  #labels for the VAR

  (s->lagDF!(df,s)).(LHS)

  var::VAR = (VAREstimate(df, LHS, RHS)).var
  display(var.ϕ₀)
  display(var.ϕ₁)

  return var
end

#simulate the VAR once
function simulateVAR!(var::VAR, σε::Float64, n::Int; initial=zeros(length(var.ϕ₀)), k::Int = length(var.ϕ₀),
    outmatrix = Matrix{Float64}(undef, k, n),
    simdf = DataFrame(
      creditspread = Vector{MFloat64}(missing, n),
      taxtnote20 = Vector{MFloat64}(missing, n),
      taxfield = Vector{MFloat64}(missing, n),
      residuals=Vector{MFloat64}(missing, n),
      munispread=Vector{MFloat64}(missing, n)))

  k::Int = length(var.ϕ₀)

  outmatrix::Matrix{Float64} = Matrix{Float64}(undef, k, n)
  outmatrix[:,1] .= initial

  shocks::Vector{Float64} = Vector{Float64}(undef, k)

  for t ∈ 1:(n-1)
    shocks .= rand(Normal(), k)
    outmatrix[:,t+1] = var.ϕ₀ .+ var.ϕ₁ * outmatrix[:,t] + var.L*shocks
  end

  simdf[:creditspread] .= outmatrix[1,:]
  simdf[:taxtnote20] .= outmatrix[2,:]
  simdf[:taxfield] .= outmatrix[3,:]
  simdf[:residuals] .= rand(Normal(0.0, σε), n)
  simdf[:munispread1t] .= simdf[:creditspread] .* (1.0 .- simdf[:taxfield]) .+
    simdf[:taxtnote20] .+ simdf[:residuals]
  simdf[:munispread] .= simdf[:creditspread] .+ simdf[:taxtnote20] .+ simdf[:residuals]

  return simdf
end

#runs many regressions
function runRegressions(df::DataFrame, var::VAR; residualfield=:resid1, nsims::Int = N_SIMS)
  σε::Float64 = std(skipmissing(df[residualfield]))
  n::Int = sum(completecases(df[RHS]))
  k::Int = length(var.ϕ₀)

  #do some pre-allocation
  outmatrix::Matrix{Float64} = Matrix{Float64}(undef, k, n)
  simdf = DataFrame(
    creditspread = Vector{MFloat64}(missing, n),
    taxtnote20 = Vector{MFloat64}(missing, n),
    taxfield = Vector{MFloat64}(missing, n),
    residuals=Vector{MFloat64}(missing, n),
    munispread1t=Vector{MFloat64}(missing, n),
    munispread=Vector{MFloat64}(missing, n))
  resultsdf = DataFrame(b0 = Vector{MFloat64}(missing, nsims),
    bcreditspread = Vector{MFloat64}(missing, nsims),
    btaxtnote20 = Vector{MFloat64}(missing, nsims))

  #set the regression spec
  xspec::CTExpr = Meta.parse("creditspread + taxtnote20")
  xnames::Vector{Symbol} = [:intercept, :creditspread, :taxtnote20]
  yspec::Symbol = :munispread

  #results
  for i ∈ 1:nsims
     simdf = simulateVAR!(var, σε, n, outmatrix=outmatrix, simdf=simdf)
     reg::CTLM = CTLM(simdf, xspec, yspec)
     resultsdf[i,:b0] = reg.β[1]
     resultsdf[i,:bcreditspread] = reg.β[2]
     resultsdf[i,:btaxtnote20] = reg.β[3]
   end

   return resultsdf
end

function makehistograms(resultsdf::DataFrame)

  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()

  push!(plotNames, "intercept")
  push!(plots, plot(resultsdf,
    x=:b0,
    Guide.title("$(plotNames[end])"),
    Guide.xlabel("value"), Guide.ylabel("count"),
    Geom.histogram(bincount=22, density=false)))

  push!(plotNames, "betacreditspread")
  push!(plots, plot(resultsdf,
    x=:bcreditspread,
    Guide.title("$(plotNames[end])"),
    Guide.xlabel("value"), Guide.ylabel("count"),
    Geom.histogram(bincount=22, density=false)))

  push!(plotNames, "betatautnote20")
  push!(plots, plot(resultsdf,
    x=:btaxtnote20,
    Guide.title("$(plotNames[end])"),
    Guide.xlabel("value"), Guide.ylabel("count"),
    Geom.histogram(bincount=22, density=false)))

    for i ∈ 1:length(plots) #write the graphs
      draw(PDF("output\\$(plotNames[i]).pdf", 9inch, 7inch), plots[i])
      #println("Graph $(plotNames[i]) written.")
    end

  show(describe(resultsdf, stats=[:mean, :std, :min, :q25, :median, :q75, :max]))

  println("\nGraphs written.")

  return nothing
end

function runPlacebo()::Nothing
  df::DataFrame = readPlot()
  var::VAR = getvar!(df)
  resultsdf = runRegressions(df, var)
  makehistograms(resultsdf)
  return nothing
end

@time runPlacebo()

max_stars_repo_path: Fall 17 Ivo Project 2.2/histtaxes/TaxPlacebo.jl

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

function _clean(O::Operation)
    @assert isa(O.op, Variable)
    return O.op
end
_clean(x::Variable) = x
function test_diffeq_inference(name, sys, iv, dvs, ps)
    @testset "ODESystem construction: $name" begin
        @test independent_variables(sys) == Set([_clean(iv)])
        @test dependent_variables(sys)   == Set(_clean.(dvs))
        @test parameters(sys)            == Set(_clean.(ps))
    end
end
function test_nlsys_inference(name, sys, vs, ps)
    @testset "NonlinearSystem construction: $name" begin
        @test dependent_variables(sys) == Set(vs)
        @test parameters(sys)          == Set(_clean.(ps))
    end
end

# Define a differential equation
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

ModelingToolkit.simplified_expr.(eqs)[1]
:(derivative(x(t), t) = σ * (y(t) - x(t))).args
de = ODESystem(eqs)
test_diffeq_inference("standard", de, t, (x, y, z), (σ, ρ, β))
generate_function(de, [x,y,z], [σ,ρ,β])
jac_expr = generate_jacobian(de)
jac = calculate_jacobian(de)
jacfun = eval(jac_expr)
# iip
f = ODEFunction(de, [x,y,z], [σ,ρ,β])
Wfact, Wfact_t = ModelingToolkit.calculate_factorized_W(de)
fw, fwt = map(eval, ModelingToolkit.generate_factorized_W(de))
du = zeros(3)
u  = collect(1:3)
p  = collect(4:6)
f(du, u, p, 0.1)
@test du == [4, 0, -16]
J = zeros(3, 3)
jacfun(J, u, p, t)
FW = zeros(3, 3)
FWt = zeros(3, 3)
fw(FW, u, p, 0.2, 0.1)
fwt(FWt, u, p, 0.2, 0.1)
# oop
f = ODEFunction(de, [x,y,z], [σ,ρ,β])
fw, fwt = map(eval, ModelingToolkit.generate_factorized_W(de))
du = @SArray zeros(3)
u  = SVector(1:3...)
p  = SVector(4:6...)
@test f(u, p, 0.1) === @SArray [4, 0, -16]
Sfw = fw(u, p, 0.2, 0.1)
@test Sfw.L ≈ UnitLowerTriangular(FW)
@test Sfw.U ≈ UpperTriangular(FW)
sol = Sfw \ @SArray ones(3)
@test sol isa SArray
@test sol ≈ -(I - 0.2*J)\ones(3)
Sfw_t = fwt(u, p, 0.2, 0.1)
@test Sfw_t.L ≈ UnitLowerTriangular(FWt)
@test Sfw_t.U ≈ UpperTriangular(FWt)
sol = Sfw_t \ @SArray ones(3)
@test sol isa SArray
@test sol ≈ -(I/0.2 - J)\ones(3)

@testset "time-varying parameters" begin
    @parameters σ′(t-1)
    eqs = [D(x) ~ σ′*(y-x),
           D(y) ~ x*(ρ-z)-y,
           D(z) ~ x*y - β*z]
    de = ODESystem(eqs)
    test_diffeq_inference("global iv-varying", de, t, (x, y, z), (σ′, ρ, β))
    @test begin
        f = eval(generate_function(de, [x,y,z], [σ′,ρ,β]))
        du = [0.0,0.0,0.0]
        f(du, [1.0,2.0,3.0], [x->x+7,2,3], 5.0)
        du ≈ [11, -3, -7]
    end

    @parameters σ(..)
    eqs = [D(x) ~ σ(t-1)*(y-x),
           D(y) ~ x*(ρ-z)-y,
           D(z) ~ x*y - β*z]
    de = ODESystem(eqs)
    test_diffeq_inference("single internal iv-varying", de, t, (x, y, z), (σ, ρ, β))
    @test begin
        f = eval(generate_function(de, [x,y,z], [σ,ρ,β]))
        du = [0.0,0.0,0.0]
        f(du, [1.0,2.0,3.0], [x->x+7,2,3], 5.0)
        du ≈ [11, -3, -7]
    end

    eqs = [D(x) ~ x + 10σ(t-1) + 100σ(t-2) + 1000σ(t^2)]
    de = ODESystem(eqs)
    test_diffeq_inference("many internal iv-varying", de, t, (x,), (σ,))
    @test begin
        f = eval(generate_function(de, [x], [σ]))
        du = [0.0]
        f(du, [1.0], [t -> t + 2], 5.0)
        du ≈ [27561]
    end
end

# Conversion to first-order ODEs #17
@derivatives D3'''~t
@derivatives D2''~t
@variables u(t) u_tt(t) u_t(t) x_t(t)
eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
       D2(x) ~ D(x) + 2]
de = ODESystem(eqs)
de1 = ode_order_lowering(de)
lowered_eqs = [D(u_tt) ~ 2u_tt + u_t + x_t + 1
               D(x_t)  ~ x_t + 2
               D(u_t)  ~ u_tt
               D(u)    ~ u_t
               D(x)    ~ x_t]
@test de1 == ODESystem(lowered_eqs)
test_diffeq_inference("first-order transform", de1, t, [u_tt, x_t, u_t, u, x], [])
du = zeros(5)
ODEFunction(de1, [u_tt, x_t, u_t, u, x], [])(du, ones(5), nothing, 0.1)
@test du == [5.0, 3.0, 1.0, 1.0, 1.0]

# Internal calculations
a = y - x
eqs = [D(x) ~ σ*a,
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
de = ODESystem(eqs)
generate_function(de, [x,y,z], [σ,ρ,β])
jac = calculate_jacobian(de)
f = ODEFunction(de, [x,y,z], [σ,ρ,β])

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z])
test_nlsys_inference("standard", ns, (x, y, z), (σ, ρ, β))
@test begin
    f = eval(generate_function(ns, [x,y,z], [σ,ρ,β]))
    du = [0.0, 0.0, 0.0]
    f(du, [1,2,3], [1,2,3])
    du ≈ [1, -3, -7]
end

@derivatives D'~t
@parameters A B C
_x = y / C
eqs = [D(x) ~ -A*x,
       D(y) ~ A*x - B*_x]
de = ODESystem(eqs)
@test begin
    f = eval(generate_function(de, [x,y], [A,B,C]))
    du = [0.0,0.0]
    f(du, [1.0,2.0], [1,2,3], 0.0)
    du ≈ [-1, -1/3]
    du ≈ f([1.0,2.0], [1,2,3], 0.0)
end

# Now nonlinear system with only variables
@variables x y z
@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z])
jac = calculate_jacobian(ns)
@testset "nlsys jacobian" begin
    @test isequal(jac[1,1], σ * -1)
    @test isequal(jac[1,2], σ)
    @test isequal(jac[1,3], 0)
    @test isequal(jac[2,1], ρ - z)
    @test isequal(jac[2,2], -1)
    @test isequal(jac[2,3], x * -1)
    @test isequal(jac[3,1], y)
    @test isequal(jac[3,2], x)
    @test isequal(jac[3,3], -1 * β)
end
nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
jac_func = generate_jacobian(ns)
f = @eval eval(nlsys_func)

# Intermediate calculations
# Define a nonlinear system
eqs = [0 ~ σ*a,
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z])
nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
jac = calculate_jacobian(ns)
jac = generate_jacobian(ns)

function lotka(u,p,t)
  x = u[1]
  y = u[2]
  [p[1]*x - p[2]*x*y,
  -p[3]*y + p[4]*x*y]
end

prob = ODEProblem(ODEFunction{false}(lotka),[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])
de, vars, params = modelingtoolkitize(prob)
ODEFunction(de, vars, params)(similar(prob.u0), prob.u0, prob.p, 0.1)

function lotka(du,u,p,t)
  x = u[1]
  y = u[2]
  du[1] = p[1]*x - p[2]*x*y
  du[2] = -p[3]*y + p[4]*x*y
end

prob = ODEProblem(lotka,[1.0,1.0],(0.0,1.0),[1.5,1.0,3.0,1.0])

de, vars, params = modelingtoolkitize(prob)
ODEFunction(de, vars, params)(similar(prob.u0), prob.u0, prob.p, 0.1)


text:
#=
DataFrames are used to represent factors
https://en.wikipedia.org/wiki/Factor_graph

    :potential is the column containing the probabilities, ::Float64
    Each variable has its own column corresponding to its assignments and named with its name

These can be obtained using the table() function
=#

mutable struct Table
    potential::DataFrame
end

Base.convert(::Type{DataFrame}, t::Table) = t.potential

for f in [:names, :unique, :size, :eltype, :setindex!]
    @eval (Base.$f)(t::Table, x...) = $f(t.potential, x...)
end
Base.getindex(t::Table, x...) = getindex(t.potential, !, x...)

for s in [:(==), :(!=)]
    @eval (Base.$s)(t::Table, f::DataFrame) = $s(t.potential, f)
    @eval (Base.$s)(f::DataFrame, t::Table) = $s(f, t.potential)
    @eval (Base.$s)(t1::Table, t2::Table) = $s(t1.potential, t2.potential)
end

DataFrames.nrow(t::Table) = nrow(t.potential)

"""
Table multiplication
"""
function Base.:*(t1::Table, t2::Table)
    f1 =t1.potential
    f2 =t2.potential

    onnames = setdiff(intersect(propertynames(f1), propertynames(f2)), [:potential])
    finalnames = vcat(setdiff(union(propertynames(f1), propertynames(f2)), [:potential]), :potential)

    if isempty(onnames)
        j = join(f1, f2, kind=:cross, makeunique=true)
    else
        j = outerjoin(f1, f2, on=onnames, makeunique=true)
    end

    j[!,:potential] = broadcast(*, j[!,:potential], j[!,:potential_1])

    return Table(j[!,finalnames])
end

"""
    sumout(t, v)

Table marginalization
"""
function sumout(t::Table, v::NodeNameUnion)
    f = t.potential

    # vcat works for single values and vectors alike (magic?)
    remainingvars = setdiff(propertynames(f), vcat(v, :potential))

    if isempty(remainingvars)
        # they want to remove all variables except for prob column
        # uh ... 'singleton' table?
        return Table(DataFrame(potential = sum(f[!,:potential])))
    else
        # note that this will fail miserably if f is too large (~1E4 maybe?)
        #  nothing I can do; there is a github issue
        return Table(combine(df -> DataFrame(potential = sum(df[!,:potential])), DataFrames.groupby(f, remainingvars)))
    end
end

"""
Table normalization
Ensures that the `:potential` column sums to one
"""
function LinearAlgebra.normalize!(t::Table)
    t.potential[!,:potential] ./= sum(t.potential[!,:potential])

    return t
end

LinearAlgebra.normalize(t::Table) = normalize!(deepcopy(t))

"""
Given a Table, extract the rows which match the given assignment
"""
function Base.partialsort(t::Table, a::Assignment)
    f = t.potential

    commonNames = intersect(propertynames(f), keys(a))
    mask = trues(size(f, 1))
    for s in commonNames
        # mask &= (f[!,s] .== a[s])
        vals = (f[!,s] .== a[s])
        for (i,v) in enumerate(vals)
            mask[i] = mask[i] & v
        end
    end

    return Table(f[mask, :])
end

"""
takes a list of observations of assignments represented as a DataFrame
or a set of data samples (without :potential),
takes the unique assignments,
and estimates the associated probability of each assignment
based on its frequency of occurrence.
"""
function Distributions.fit(::Type{Table}, f::DataFrame)
    w = ones(size(f, 1))
    t = f
    if hasproperty(f, :potential)
        t = f[:, propertynames(t) .!= :potential]
        w = f[!,:potential]
    end
    # unique samples
    tu = unique(t)
    # add column with probabilities of unique samples
    tu[!,:potential] = Float64[sum(w[Bool[tu[j,:] == t[i,:] for i = 1:size(t,1)]]) for j = 1:size(tu,1)]
    tu[!,:potential] /= sum(tu[!,:potential])

    return Table(tu)
end

# """
# TODO: what is this for?
# """
# function estimate_convergence(t::Table, a::Assignment)
#     f = t.potential
#
#     n = size(f, 1)
#     p = zeros(n)
#     w = ones(n)
#     if hasproperty(f, :potential)
#         w = f[!,:potential]
#     end
#
#     dfindex = find([hasproperty(a, n) for n in names(f)])
#     dfvalues = [a[n] for n in names(f)[!,dfindex]]'
#
#     cumWeight = 0.0
#     cumTotalWeight = 0.0
#     for i in 1:n
#         if convert(Array, f[i, dfindex]) == dfvalues
#             cumWeight += w[i]
#         end
#         cumTotalWeight += w[i]
#         p[i] = cumWeight / cumTotalWeight
#     end
#
#     return p
# end


text:
using JuMP
# using CPLEX
using Gurobi
# using Suppressor
using XLSX
# using JLD
# using Printf
# using Clustering
# using Distributed
using DataFrames
# using Plots
# pyplot()
# Plots.PyPlotBackend()
# using StatsPlots
# using LaTeXStrings
# gurobi_env = @suppress Gurobi.Env()


text:
alphabet = "0123456789abcdefghijklmnopqrstuvwxyz" # Only supports thru base 36
intfloor(x) = convert(Int,floor(x))

"
    convertintrep(n::Int, base::Int)

Converts the input `n` ≧ 0 into its string representation with respect to
`base`.
"
function convertintrep(n::Int, base::Int)
    result = ""
    while n > 1
        result *= "$(alphabet[1 + (n % base)])"
        n = intfloor(n / base)
    end

    result[end:-1:1]
end

"
    convertfracrep(n::Real, base::Int, precision)

Converts the input 0 ≦ `n` < 1 into its string representation with respect to
`base` for `precision` number of places.
"
function convertfracrep(n::Real, base::Int, precision)
    result = "."
    for _ in 1:precision
        n *= base
        result *= "$(alphabet[1 + intfloor(n)])"
        n -= intfloor(n)
    end

    result
end

"
    convertfloat(n::Real, base::Int, precision)

Converts a real number `n` into its string representation with respect to
`base` by a fixed `precision`. This is a concatenation of `convertintrep` and
`convertfracrep`.
"
function convertfloat(n::Real, base::Int, precision)
    result = ""
    if n < 0
        result *= "-"
        n *= -1
    end
    intpart = intfloor(n)
    fracpart = n - floor(n)
    result *= convertintrep(intpart, base) * convertfracrep(fracpart, base, precision)
end

int2dec(n::Int) = convertintrep(n, 10)
int2bin(n::Int) = "0b" * convertintrep(n, 2)
int2hex(n::Int) = "0x" * convertintrep(n, 16)


text:
struct Query
    xdoc::XMLDocument
    query::XMLElement
    dataset::XMLElement
end

function Query(dataset::String)

    # create an empty XML document
    xdoc = XMLDocument()

    # create & attach a root node
    query = create_root(xdoc, "Query")
    set_attribute(query, "virtualSchemaName", "default")
    set_attribute(query, "formatter", "TSV")
    set_attribute(query, "header", "1")
    set_attribute(query, "uniqueRows", "1")
    set_attribute(query, "count", "")
    set_attribute(query, "datasetConfigVersion", "0.6")

    # create the first child
    ds = new_child(query, "Dataset")
    set_attributes(ds, name = dataset, interface = "default")

    Query(xdoc, query, ds)
end

Query(dataset::Dataset) = Query(dataset.dataset)

"""
    Query(dataset::Dataset, args...)

Build a `Query` for the given `Dataset` and combination
of `Filters` and `Attributes`. An `Interval` from `GenomicFeatures` can also
be provided to filter by genomic position.

A `Query` build the request in `XML` format but does not execute it. To do so you can
call `BioMart.execute` or directly call the `Query` object.

To build and execute a `Query` at once use `BioMart.query`.

Example:

    q = BioMart.Query(
        BioMart.Dataset("maj_gene_ensembl"),
        Interval("1", 1, 200000),
        BioMart.Attribute("external_gene_name"),
        BioMart.Attributes("ensembl_gene_id", "ensembl_transcript_id", "strand"),
        BioMart.Filter(strand = "1")
    )
    q() #BioMart.execute(q)
"""
function Query(dataset::Dataset, args...)
    q = Query(dataset)
    for el in args
        push!(q, el)
    end
    q
end

Base.show(io::IO, q::Query) = show(q.xdoc)

"""
    execute(q::Query)

Exectute a `Query` and returns a `DataFrame` with the results.
"""
function execute(q::Query)
    query = HTTP.URIs.escapeuri(string(q.xdoc))
    r = HTTP.get(string(BIOMART_URL, "?query=", query))
    CSV.read(r.body, delim = '\t')
end

(q::Query)() = execute(q)

"""
    query(args...)

Build a `Query` using `args`, exectute it immediately and returns the results.

"""
query(args...) = Query(args...)()


text:
struct EquilibResult
    sys::System
    T::Real # temperature
    X::Vector{Float64} # molar fractions of components
    Y::Vector{Float64} # molar amount of phases
    x::Array{Float64,2} # molar fractions of components in each phase
    y::Array{Float64,3} # molar fractions of constituents in each sublattice
end

function Base.display(res::EquilibResult)
    elems = res.sys.elems
    conss = res.sys.conss
    phass = res.sys.phass
    I = length(elems)
    J = length(conss)
    K = length(phass)
    @assert size(res.x) == (K,I)
    for k in 1:K
        phas = phass[k]
        res.Y[k] < 1e-4 && continue
        @printf "%s; %s: %.4f\n" phas.name constitutionstring(phas) res.Y[k] / sum(res.Y)
        !any(latt -> length(latt) > 1, phas.cons) && continue
        for i in 1:I
            res.x[k,i] < 1e-4 && continue
            @printf "\t%s: %.4f\n" elems[i].name res.x[k,i]
        end
        cons = filter(latt -> !isempty(latt), phas.cons)
        S = length(cons)
        S < 2 && length(cons[1]) == length(elems) && continue
        for s in 1:S
            length(phas.cons[s]) < 2 && continue
            if length(filter(l -> l ≠ [], phas.cons)) > 1
                println("\t" * constitutionstring(phas, s))
            end
            for j in 1:J
                consname = conss[j]
                if consname in cons[s]
                    res.y[k,s,j] < 1e-4 && continue
                    @printf "\t\t%s: %.4f\n" subscript(consname) res.y[k,s,j]
                end
            end
        end
    end
end

Base.print(res::EquilibResult) = Base.display(res)

function equilib(sys::System, X = equiatom(sys), T = 298.15)
    # JuMP model
    model = sys.model

    I = sys.nelem # number of components
    J = sys.ncons # number of constituents
    K = sys.nphas # number of phases
    S = sys.nlatt # maximum number of sites in a sublattice
    n = sys.nsite # number of sites in each sublattice

    if length(X) ≠ I
        error("Invalid length of molar amount vector of components")
    end

    # Variables
    _T = model[:T] # temperature
    _X = model[:X] # molar amount of components
    Y = model[:Y] # moalr amount of phases
    x = model[:x] # molar fraction of components in each phase
    y = model[:y] # site fraction in each sublattice

    # Fix temperature, T
    # TODO: Using fix() results in calling f(T) at T = 0 for some reason
    Tl, Tu = sys.temp
    if T < Tl || T > Tu
        error("Temperature out of range")
    end
    set_lower_bound(_T, T)
    set_upper_bound(_T, T)

    # Fix X[i]
    # TODO: Using fix() results in X[i] = 0 for some reason
    for i in 1:I
        set_lower_bound(_X[i], X[i])
        set_upper_bound(_X[i], X[i])
    end

    # Determine maximum values of Y[k]
    Y_max = [ sum( X[i] for i in 1:I ) / sum( n[k,s] for s in 1:S ) for k in 1:K ]
    for k in 1:K
        set_upper_bound(Y[k], Y_max[k])
        for s in 1:S, j in 1:J
            is_fixed(y[k,s,j]) && continue
        end
    end

    @debug model

    optimize!(model)

    return EquilibResult(sys, T, X, [ value(Y[k]) for k in 1:K ],
                                    [ value(x[k,i]) for k in 1:K, i in 1:I ],
                                    [ value(y[k,s,j]) for k in 1:K, s in 1:S, j in 1:J ])
end


text:
#### randomized SVD (from Jiahao Chen, based on http://arxiv.org/pdf/0909.4061.pdf)

import Base.LinAlg.SVD

#The simplest possible randomized svd
#Inputs
#   A: input matrix
#   n: Number of singular value/vector pairs to find
#   p: Number of extra vectors to include in computation
function rsvd(A, n, p=0)
    Q = rrange(A, n, p=p)
    rsvd_direct(A, Q)
end

#Algorithm 4.4: randomized subspace iteration
#A must support size(A), multiply and transpose multiply
#p is the oversampling parameter
#q controls the accuracy of the subspace found; it is the "number of power iterations"
#A good heuristic is that when the original scheme produces a basis whose
#approximation error is within a factor C of the optimum, the power scheme produces
#an approximation error within C^(1/(2q+1)) of the optimum.
function rrange(A, l::Integer; p::Integer=5, q::Integer=3)
    p≥0 || error()
    m, n = size(A)
    l <= m || error("Cannot find $l linearly independent vectors of $m x $n matrix")
    Ω = randn(n, l+p)
    Q = q_from_qr(A*Ω)
    for t=1:q
        Q = q_from_qr(A'*Q)
        Q = q_from_qr(A*Q)
    end
    Q = p==0 ? Q : Q[:,1:l]
end

function q_from_qr(Y, l::Integer=-1)
    Q = full(qrfact!(Y)[:Q])
    Q = l<0 ? Q : Q[:,1:l]
end

#Algorithm 5.1: direct SVD
#More accurate
function rsvd_direct(A, Q)
    B=Q'A
    S=svdfact!(B)
    SVD(Q*S[:U], S[:S], S[:Vt])
end


text:
"""
三角形有关的功能
"""

using Printf

struct AbstractTriangle{T} <: AbstractPolygon{T}
    vertex :: Array{Point2D, 1}
    edges :: Array{Segment, 1}
    center :: Point2D
    area :: Float64
end

"""
构造一个普通的三角形
"""
function Triangle(pt1, pt2, pt3)
    verts = [pt1, pt2, pt3]
    edges = Array{Segment, 1}(undef, 3)
    edges[1] = Segment(pt1, pt2)
    edges[2] = Segment(pt2, pt3)
    edges[3] = Segment(pt3, pt1)
    center = Point2D(
        sum([pt1.x, pt2.x, pt3.x]) / 3.0,
        sum([ver.y for ver in verts]) / 3.0
    )
    vec1 = pt2 - pt1
    vec2 = pt3 - pt1
    tarea =  0.5abs(vec1.x * vec2.y - vec1.y * vec2.x)
    return AbstractTriangle{:COM}(verts, edges, center, tarea)
end

"""
构造一个普通的三角形
"""
function Triangle(verts::Array{Point2D, 1})
    @assert length(verts) == 3
    return Triangle(verts[1], verts[2], verts[3])
end

"""
构造一个直角三角形，第一个角必须是直角
"""
function RtTriangle(pt1, pt2, pt3)
    verts = [pt1, pt2, pt3]
    edges = Array{Segment, 1}(undef, 3)
    edges[1] = Segment(pt1, pt2)
    edges[2] = Segment(pt2, pt3)
    edges[3] = Segment(pt3, pt1)
    center = middle_point(
        pt1,
        middle_point(pt2, pt3),
        sc1=1., sc2=1.414
    )
    vec1 = pt2 - pt1
    vec2 = pt3 - pt1
    tarea =  0.5abs(vec1.x * vec2.y - vec1.y * vec2.x)
    return AbstractTriangle{:RT}(verts, edges, center, tarea)
end

"""
构造一个直角三角形，会检查第一个角是否是直角
"""
function RtTriangle(verts::Array{Point2D, 1})
    @assert length(verts) == 3
    vec1 = verts[2] - verts[1]
    vec2 = verts[3] - verts[1]
    @assert abs(vec1.x * vec2.x + vec1.y * vec2.y) < 1e-6
    return RtTriangle(verts[1], verts[2], verts[3])
end

"""
构造一个等边三角形
"""
function EqTriangle(pt1, pt2, pt3)
    verts = [pt1, pt2, pt3]
    edges = Array{Segment, 1}(undef, 3)
    edges[1] = Segment(pt1, pt2)
    edges[2] = Segment(pt2, pt3)
    edges[3] = Segment(pt3, pt1)
    center = Point2D(
        sum([pt1.x, pt2.x, pt3.x]) / 3.0,
        sum([ver.y for ver in verts]) / 3.0
    )
    vec1 = pt2 - pt1
    vec2 = pt3 - pt1
    tarea =  0.5abs(vec1.x * vec2.y - vec1.y * vec2.x)
    return AbstractTriangle{:EQ}(verts, edges, center, tarea)
end

"""
构造一个等边三角形，会检查边长
"""
function EqTriangle(verts::Array{Point2D, 1})
    @assert length(verts) == 3
    ret = EqTriangle(verts[1], verts[2], verts[3])
    @assert abs(ret.edges[1].length - ret.edges[2].length) < 1e-6
    @assert abs(ret.edges[3].length - ret.edges[2].length) < 1e-6
    return ret
end

"""
返回一个三角形的面积
"""
function area(tri::AbstractTriangle{N}) where N
    vec1 = tri.vertex[2] - tri.vertex[1]
    vec2 = tri.vertex[3] - tri.vertex[1]
    return 0.5abs(vec1.x * vec2.y - vec1.y * vec2.x)
end

#display函数在中间加::MIME"text/plain"
function Base.show(io::IO, tri::AbstractTriangle{N}) where N
    name = Dict(
        :COM => "Triangle",
        :RT => "RTTriangle",
        :EQ => "EqTriangle"
    )
    @printf(io, "%s\n center: %s\n vertex: %s",
    name[N], tri.center, tri.vertex)
end


text:
using Test, LinearAlgebra, Distributions, AdvancedHMC, Random, ForwardDiff
using AdvancedHMC.Adaptation: WelfordVar, NaiveVar, WelfordCov, NaiveCov, get_estimation, get_estimation, reset!

# Check that the estimated variance is approximately correct.
@testset "Online v.s. naive v.s. true var/cov estimation" begin
    D = 10
    T = Float64
    sz = (D,)
    n_samples = 100_000

    var_welford = WelfordVar{T}(sz)
    var_naive = NaiveVar{T}(sz)
    var_estimators = [var_welford, var_naive]
    cov_welford = WelfordCov{T}(sz)
    cov_naive = NaiveCov{T}(sz)
    cov_estimators = [cov_welford, cov_naive]
    estimators = [var_estimators..., cov_estimators...]

    for dist in [MvNormal(zeros(D), ones(D)), Dirichlet(D, 1)]
        for _ = 1:n_samples
            s = rand(dist)
            for estimator in estimators
                push!(estimator, s)
            end
        end

        @test get_estimation(var_welford) ≈ get_estimation(var_naive) atol=0.1D
        for estimator in var_estimators
            @test get_estimation(estimator) ≈ var(dist) atol=0.1D
        end

        @test get_estimation(cov_welford) ≈ get_estimation(cov_naive) atol=0.1D^2
        for estimator in cov_estimators
            @test get_estimation(estimator) ≈ cov(dist) atol=0.1D^2
        end

        for estimator in estimators
            reset!(estimator)
        end
    end
end

@testset "MassMatrixAdaptor constructors" begin
    θ = [0.0, 0.0, 0.0, 0.0]
    pc1 = MassMatrixAdaptor(UnitEuclideanMetric) # default dim = 2
    pc2 = MassMatrixAdaptor(DiagEuclideanMetric)
    pc3 = MassMatrixAdaptor(DenseEuclideanMetric)

    # Var adaptor dimention should be increased to length(θ) from 2
    AdvancedHMC.adapt!(pc1, θ, 1.)
    AdvancedHMC.adapt!(pc2, θ, 1.)
    AdvancedHMC.adapt!(pc3, θ, 1.)
    @test AdvancedHMC.Adaptation.getM⁻¹(pc2) == ones(length(θ))
    @test AdvancedHMC.Adaptation.getM⁻¹(pc3) == LinearAlgebra.diagm(0 => ones(length(θ)))
end

@testset "Stan HMC adaptors" begin
    θ = [0.0, 0.0, 0.0, 0.0]

    adaptor1 = StanHMCAdaptor(
        MassMatrixAdaptor(UnitEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    adaptor2 = StanHMCAdaptor(
        MassMatrixAdaptor(DiagEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    adaptor3 = StanHMCAdaptor(
        MassMatrixAdaptor(DenseEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    for a in [adaptor1, adaptor2, adaptor3]
        AdvancedHMC.initialize!(a, 1_000)
        @test a.state.window_start == 76
        @test a.state.window_end == 950
        @test a.state.window_splits == [100, 150, 250, 450, 950]
        AdvancedHMC.adapt!(a, θ, 1.)
    end
    @test AdvancedHMC.Adaptation.getM⁻¹(adaptor2) == ones(length(θ))
    @test AdvancedHMC.Adaptation.getM⁻¹(adaptor3) == LinearAlgebra.diagm(0 => ones(length(θ)))

    @test_deprecated StanHMCAdaptor(
        1_000,
        MassMatrixAdaptor(DiagEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )

    @testset "buffer > `n_adapts`" begin
        AdvancedHMC.initialize!(
            StanHMCAdaptor(
                MassMatrixAdaptor(DenseEuclideanMetric),
                NesterovDualAveraging(0.8, 0.5)
            ),
            100
        )
    end
end

let D=10
    function runnuts(ℓπ, metric; n_samples=3_000)
        n_adapts = 1_500

        θ_init = rand(D)

        h = Hamiltonian(metric, ℓπ, ForwardDiff)
        prop = NUTS(Leapfrog(find_good_stepsize(h, θ_init)))
        adaptor = StanHMCAdaptor(
            MassMatrixAdaptor(metric),
            StepSizeAdaptor(0.8, prop.integrator)
        )
        samples, stats = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; verbose=false)
        return (samples=samples, stats=stats, adaptor=adaptor)
    end

    @testset "Adapted mass v.s. true variance" begin
        n_tests = 5

        @testset "DiagEuclideanMetric" begin
            for _ in 1:n_tests
                Random.seed!(1)

                # Random variance
                σ = ones(D) + abs.(randn(D))

                # Diagonal Gaussian
                target = MvNormal(zeros(D), σ)
                ℓπ = θ -> logpdf(target, θ)

                res = runnuts(ℓπ, DiagEuclideanMetric(D))
                @test res.adaptor.pc.var ≈ σ .^ 2 rtol=0.2

                res = runnuts(ℓπ, DenseEuclideanMetric(D))
                @test res.adaptor.pc.cov ≈ diagm(0 => σ .^ 2) rtol=0.25
            end
        end

        @testset "DenseEuclideanMetric" begin
            for _ in 1:n_tests
                # Random covariance
                m = randn(D, D)
                Σ = m' * m

                # Correlated Gaussian
                target = MvNormal(zeros(D), Σ)
                ℓπ = θ -> logpdf(target, θ)

                res = runnuts(ℓπ, DiagEuclideanMetric(D))
                @test res.adaptor.pc.var ≈ diag(Σ) rtol=0.2

                res = runnuts(ℓπ, DenseEuclideanMetric(D))
                @test res.adaptor.pc.cov ≈ Σ rtol=0.25
            end
        end

    end

    @testset "Initialisation adaptor by metric" begin
        target = MvNormal(zeros(D), 1)
        ℓπ = θ -> logpdf(target, θ)

        mass_init = fill(0.5, D)
        res = runnuts(ℓπ, DiagEuclideanMetric(mass_init); n_samples=1)
        @test res.adaptor.pc.var == mass_init

        mass_init = diagm(0 => fill(0.5, D))
        res = runnuts(ℓπ, DenseEuclideanMetric(mass_init); n_samples=1)
        @test res.adaptor.pc.cov == mass_init
    end
end

@testset "Deprecation" begin
    dim = 10
    @test_deprecated Preconditioner(UnitEuclideanMetric(dim))
    @test_deprecated Preconditioner(DiagEuclideanMetric(dim))
    @test_deprecated Preconditioner(DenseEuclideanMetric(dim))
    @test_deprecated Preconditioner(UnitEuclideanMetric)
    @test_deprecated Preconditioner(DiagEuclideanMetric)
    @test_deprecated Preconditioner(DenseEuclideanMetric)
    for T in [Float32, Float64]
        @test_deprecated Preconditioner(T, UnitEuclideanMetric)
        @test_deprecated Preconditioner(T, DiagEuclideanMetric)
        @test_deprecated Preconditioner(T, DenseEuclideanMetric)
    end
    @test_deprecated NesterovDualAveraging(0.8, Leapfrog(0.1))
    @test_deprecated StanHMCAdaptor(100, MassMatrixAdaptor(UnitEuclideanMetric(dim)), StepSizeAdaptor(0.8, Leapfrog(0.1)))
end


text:
println("Times for 8 argument example using all macros")
println("Load Time (using ArgMacros)")
@time using ArgMacros
using Test

# Commands Set 1
let
    empty!(ARGS)
    append!(ARGS, ["TEST STRING F", "-deeee", "30", "3.14", "-b=6.28", "--cc", "ArgMacros", "-a", "2"])

    println("Inline Arguments Time with Precompile")
    let
        @time @inlinearguments begin
            @argumentrequired Int a "-a" "--aa"
            @argtest a (<(10))
            @argumentdefault Float64 10 b "-b"
            @argumentoptional Symbol c "--cc"
            @argumentflag d "-d"
            @argumentcount e "-e"

            @positionalrequired String f
            @positionaldefault Int 20 g
            @argtest g (x -> x % 10 == 0)
            @positionaloptional Float64 h
        end
    end

    println("Inline Arguments Time")
    @time @inlinearguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    println("Struct Arguments Time")
    @time begin
        @structarguments false ArgsStruct begin
            @argumentrequired Int a "-a" "--aa"
            @argtest a (<(10))
            @argumentdefault Float64 10 b "-b"
            @argumentoptional Symbol c "--cc"
            @argumentflag d "-d"
            @argumentcount e "-e"

            @positionalrequired String f
            @positionaldefault Int 20 g
            @argtest g (x -> x % 10 == 0)
            @positionaloptional Float64 h
        end

        argsstruct = ArgsStruct()
    end

    println("Tuple Arguments Time")
    @time argstuple = @tuplearguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    println("Dict Arguments Time")
    @time argsdict = @dictarguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    @testset "Correct, All Arguments" begin
        @testset "Inline" begin
            @test a == 2
            @test b == 6.28
            @test c == :ArgMacros
            @test d
            @test e == 4
            @test f == "TEST STRING F"
            @test g == 30
            @test h  == 3.14
        end

        @testset "Struct" begin
            @test argsstruct.a == 2
            @test argsstruct.b == 6.28
            @test argsstruct.c == :ArgMacros
            @test argsstruct.d
            @test argsstruct.e == 4
            @test argsstruct.f == "TEST STRING F"
            @test argsstruct.g == 30
            @test argsstruct.h  == 3.14
        end

        @testset "Tuple" begin
            @test argstuple.a == 2
            @test argstuple.b == 6.28
            @test argstuple.c == :ArgMacros
            @test argstuple.d
            @test argstuple.e == 4
            @test argstuple.f == "TEST STRING F"
            @test argstuple.g == 30
            @test argstuple.h  == 3.14
        end

        @testset "Dict" begin
            @test argsdict[:a] == 2
            @test argsdict[:b] == 6.28
            @test argsdict[:c] == :ArgMacros
            @test argsdict[:d]
            @test argsdict[:e] == 4
            @test argsdict[:f] == "TEST STRING F"
            @test argsdict[:g] == 30
            @test argsdict[:h]  == 3.14
        end
    end
end

# Commands Set 2
let
    empty!(ARGS)
    append!(ARGS, ["OTHER TEST STRING F", "--aa=5"])

    @inlinearguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    @structarguments false ArgsStruct begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    argsstruct = ArgsStruct()

    argstuple = @tuplearguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    argsdict = @dictarguments begin
        @argumentrequired Int a "-a" "--aa"
        @argtest a (<(10))
        @argumentdefault Float64 10 b "-b"
        @argumentoptional Symbol c "--cc"
        @argumentflag d "-d"
        @argumentcount e "-e"

        @positionalrequired String f
        @positionaldefault Int 20 g
        @argtest g (x -> x % 10 == 0)
        @positionaloptional Float64 h
    end

    @testset "Correct, Minimal Arguments" begin
        @testset "Inline" begin
            @test a == 5
            @test b == 10.0
            @test isnothing(c)
            @test !d
            @test e == 0
            @test f == "OTHER TEST STRING F"
            @test g == 20
            @test isnothing(h)
        end

        @testset "Struct" begin
            @test argsstruct.a == 5
            @test argsstruct.b == 10.0
            @test isnothing(argsstruct.c)
            @test !argsstruct.d
            @test argsstruct.e == 0
            @test argsstruct.f == "OTHER TEST STRING F"
            @test argsstruct.g == 20
            @test isnothing(argsstruct.h)
        end

        @testset "Tuple" begin
            @test argstuple.a == 5
            @test argstuple.b == 10.0
            @test isnothing(argstuple.c)
            @test !argstuple.d
            @test argstuple.e == 0
            @test argstuple.f == "OTHER TEST STRING F"
            @test argstuple.g == 20
            @test isnothing(argstuple.h)
        end

        @testset "Dict" begin
            @test argsdict[:a] == 5
            @test argsdict[:b] == 10.0
            @test isnothing(argsdict[:c])
            @test !argsdict[:d]
            @test argsdict[:e] == 0
            @test argsdict[:f] == "OTHER TEST STRING F"
            @test argsdict[:g] == 20
            @test isnothing(argsdict[:h])
        end
    end
end

# Should Create non-zero exit codes
@testset "Incorrect Arguments Rejected" begin
    # Make sure the test script is set up and working with correct args
    @test success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" -deeee 30 3.14 -b=6.28 --cc ArgMacros -a 2`)
    # These should fail
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- 30 3.14 -d -eeee -b 6.28 --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- TEST STRING F 30 3.14 -d -eeee -b 6.28 --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" 23.3 3.14 -d -eeee -b 6.28 --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" -eeee 25 3.14 -d -b 6.28 --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" 30 cat -deeee -b=6.28 --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" 30 3.14 -deeee -b=dog --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" 30 3.14 -deeee -b=6.28 --cc ArgMacros -a bird`)
    @test !success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" 30 3.14 -deeee -b=6.28 --cc ArgMacros`)

end

@testset "Help Triggered Properly" begin
    # Make sure the test script is set up and working with correct args
    @test  success(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" -deeee 30 3.14 -b=6.28 --cc ArgMacros -a 2`)
    @test !success(`julia --project=.. get_args_fail_at_exit.jl    -- "TEST STRING F" -deeee 30 3.14 -b=6.28 --cc ArgMacros -a 2`)
    # These will fail on purpose if not exited early with code 0
    # We want help to trigger, so they should successfully exit early and not fail
    @test success(`julia --project=.. get_args_fail_at_exit.jl -- "TEST STRING F" 30 3.14 --help -deeee -b=6.28 --cc ArgMacros -a 2`)
    @test success(`julia --project=.. get_args_fail_at_exit.jl -- "TEST STRING F" 30 -h 3.14 -deeee -b=dog --cc ArgMacros -a 2`)
    @test success(`julia --project=.. get_args_fail_at_exit.jl -- --help`)
end

println("Time for 8 argument script cold launch and parsing including Julia start")
@time run(`julia --project=.. get_args_no_fail_at_exit.jl -- "TEST STRING F" -deeee 30 3.14 -b=6.28 --cc ArgMacros -a 2`)
println("Time for Julia start only")
@time run(`julia --project=.. blank_script.jl`)


text:
function compute_inttscale(mw, maxtime::Integer)
    #window = length(mw)
    #_autocor = StatsBase.autocor(mw, 1:window-1)
    #_autocor = StatsBase.autocor(mw)
    _autocor = StatsBase.autocor(mw, 0:maxtime)

    if isnan(sum(_autocor))
        return 0.0
    end

    # summation autocor
    x = 1:length(_autocor)
    NumericalIntegration.integrate(x, _autocor, SimpsonEven())
end

function compute_inttscale(x::AbstractVector, _autocor::AbstractVector)
    maxtime = 0

    function findmaxint(_autocor)
        for i in 1:length(_autocor) - 1
            if _autocor[i] * _autocor[i + 1] < 0
                return i
            end
        end

        return length(_autocor)
    end

    maxtime = findmaxint(_autocor)

    NumericalIntegration.integrate(x[1:maxtime], _autocor[1:maxtime], SimpsonEven())
end

function diff_mean(mw)
    #window = length(mw)
    #_autocor = StatsBase.autocor(mw, 1:window-1)
    #_autocor = StatsBase.autocor(mw)
    _autocor = StatsBase.autocor(mw, 0:50)

    if isnan(sum(_autocor))
        return 0.0
    end

    # summation autocor
    x = 1:length(_autocor)
    integrate(x, _autocor, SimpsonEven())
end

function mean_aucotor(ta::TimeArray, window::Integer; padding::Bool = false)
    A = values(ta)
    # moving window
    mw = map(1:length(ta)-window+1) do i
        A[i:i+(window-1)]
    end
    padding && (mw = [fill(zeros(window), window - 1); mw])

    _autocors = []
    for w in mw
        #_auto = StatsBase.autocor(w, 1:window-1)
        #_auto = StatsBase.autocor(w)
        _auto = StatsBase.autocor(w, 0:50)

        if !isnan(sum(_auto))
            push!(_autocors, _auto)
        end
    end

    mean(hcat(_autocors...), dims=[2])
end

"""
    pdf(data)

estimate probability density function
    by kernel density estimation
"""
function pdf(data::AbstractVector, npts::Integer; method::Symbol=:kernel)
    if method == :kernel
        _pdfKDE = KernelDensity.kde(data, npoints = npts)
        x = _pdfKDE.x
        y = _pdfKDE.density
    elseif method == :linear
        _pdfKDE = KernelDensity.kde_lscv(data, npoints = npts)
        x = _pdfKDE.x
        y = _pdfKDE.density
    elseif method == :gibbs
        #bd = kde!(float.(data))
        #bmin, bmax = getKDERange(bd)
        #x = range(bmin, stop=bmax, length=npts)
        #y = evaluateDualTree(bd, x)
    elseif method == :manual
        # Freedman-Diaconis’s Rule
        binsize = 2.0 * iqr(data) * length(data)^(-1//3) * 3
        # Scott’s Rule
        # binsize = 3.49 * std(data) * length(data)^(-1//3) * 3
        # Rice's Rule
        # binsize = length(data)^(-1//3) * 2
        @show binsize
        minx = float.(minimum(data))
        maxx = float.(maximum(data))
        nbins = Int(round((maxx - minx) / float(binsize)))

        _x = minx:binsize:maxx
        x = minx:binsize:(maxx - binsize)
        y = map(1:length(x)) do i
            count(_i -> (_x[i] <= _i < _x[i + 1]), data) / float(length(data))
        end
    end

    x, y
end


text:
function call_expr(signature)
    if signature.head == :call
        return signature
    end
    return call_expr(signature.args[1])
end

function argnames(args::Array)
    tmpcount = 0
    out = []
    for a in args
        name = argname(a)
        if name == nothing
            tmpcount += 1
            name = gensym("_$tmpcount")
        end
        push!(out, name)
    end
    out
end
argname(x::Symbol) = (x)
function argname(e::Expr)
    @assert e.head == Symbol("::")
    return length(e.args) == 2 ? (e.args[1]) : nothing
end
argtypes(args::Array) = [argtype(a) for a in args]
argtype(x::Symbol) = :(typeof($x))
function argtype(e::Expr)
    @assert e.head == Symbol("::")
    e.args[end]
end

function paramexprs(arr::Array)
    names = argnames(arr)
    return [paramexpr(names[i], arr[i]) for i in 1:length(arr)]
end
paramexpr(name, e::Symbol) = (e)
paramexpr(name, e::Expr) = :($(name)::$(argtype(e)))