text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DifferentialEquations
using DynamicalSystems
using DelayEmbeddings
using DelimitedFiles

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

## This script was for testing reasons only and had no impact of any results
# shown in the paper

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we look at a chaotic system: The Roessler attractor in the funnel regime

## Integrate system and determine theiler window

# # set time interval for integration
# N = 5000 # number of samples
# dt_tra = 0.02
# tspan = (0.0, (dt_tra*N))
#
# function roessler!(du,u,p,t)
#     x,y,z = u
#     a,b,c = p
#     du[1] = -(y+z) # dx
#     du[2] = x + (a*y) # dy
#     du[3] = b+z*(x-c) # dz
# end
#
# function multiplicative_noise!(du,u,p,t)
#     du[1] = .3 * u[1]
#     du[2] = .3 * u[2]
#     du[3] = du[3]
# end
#
# a = .2
# b = .2
# c = 5.7
# p = [a,b,c]
#
# function lorenz(du,u,p,t)
#   du[1] = 10.0(u[2]-u[1])
#   du[2] = u[1]*(28.0-u[3]) - u[2]
#   du[3] = u[1]*u[2] - (8/3)*u[3]
# end
#
# # random initial condition
# u0 = [rand(1); rand(1); rand(1)]
#
# prob = SDEProblem(roessler!, multiplicative_noise!, u0, tspan, p)
# sol = solve(prob)
#
# tr = zeros(N,3)
# @inbounds for i = 1:N
#     tr[i,1] = sol.u[i][1]
#     tr[i,2] = sol.u[i][2]
#     tr[i,3] = sol.u[i][3]
# end
#
# figure()
# plot3D(tr[:,1],tr[:,2],tr[:,3])

# roe = Systems.roessler(a = a, b = b, c = c)
roe = Systems.roessler()

# set time interval for integration
N = 5000 # number of samples
dt_tra = 0.05
t = 0:dt_tra:(dt_tra*N)

tr = trajectory(roe, (N*dt_tra); dt = dt_tra, Ttr = (2000*dt_tra))
tr = regularize(tr)

σ = .1
tra = zeros(size(tr))
tra[:,1] = tr[:,1] .+ σ*randn(size(tra,1))
tra[:,2] = tr[:,2] .+ σ*randn(size(tra,1))
tra[:,3] = tr[:,3] .+ σ*randn(size(tra,1))
tr = Dataset(tra)

w1 = estimate_delay(tr[:,1], "mi_min")
w2 = estimate_delay(tr[:,2], "mi_min")

w = maximum(hcat(w1,w2))
taus = 0:100

## Perform reconstructions and compute according L-value

# L-value of reference
L_ref = uzal_cost(tr, Tw = (4*w), w = w, samplesize=1)

#Standard TDE
@time Y_tde, _ = standard_embedding_cao(tr[:,2])
L_tde = uzal_cost(Y_tde, Tw = (4*w), w = w, samplesize=1)

#MDOP embedding
println("Computation time MDOP:")
@time begin
    tw1 = mdop_maximum_delay(tr[:,1])
    tw2 = mdop_maximum_delay(tr[:,2])
    lm1 = DelayEmbeddings.findlocalminima(tw1[2])
    if lm1[1] ≤ 2
        lmm1 = try lm1[2]
        catch
            lm1[1]
        end
    else
        lmm1 = lm1[1]
    end
    lm2 = DelayEmbeddings.findlocalminima(tw2[2])
    if lm2[1] ≤ 2
        lmm2 = try lm2[2]
        catch
            lm2[1]
        end
    else
        lmm2 = lm2[1]
    end
    tau_max = maximum(hcat(lmm1, lmm2))
    taus_mdop = 0:tau_max

    Y_mdop, τ_vals_mdop, ts_vals_mdop, FNNs_mdop , βs_mdop = mdop_embedding(tr[:,1:2];
                                                        τs = taus_mdop , w = w)
    L_mdop = uzal_cost(Y_mdop, Tw = (4*w), w = w, samplesize=1)
end

#Garcia & Almeida
println("Computation time Garcia & Almeida method:")
@time Y_GA, τ_vals_GA, ts_vals_GA, FNNs_GA , ns_GA = garcia_almeida_embedding(tr[:,1:2];
                                                    τs = taus , w = w, T = w)
L_GA = uzal_cost(Y_GA, Tw = (4*w), w = w, samplesize=1)

#Pecuzal
println("Computation time pecuzal method:")
@time Y_pec, τ_vals_pec, ts_vals_pec, Ls_pec , εs_pec = pecuzal_embedding_update(tr[:,1:2];
                                                            τs = taus , w = w)
L_pec = minimum(Ls_pec)

## compute other evaluation statistics

@time mfnn1, mfnn2, mfnn3, mfnn4, f1, f2, f3, f4, RQA_ref, RQA1, RQA2,
                                            RQA3, RQA4, R_ref, R1, R2, R3, R4 =
        perform_recurrence_analysis(tr, Dataset(Y_tde), Dataset(Y_mdop),
                    Dataset(Y_GA), Dataset(Y_pec); ε = 0.08, w = w, kNN = 10)

NN = 100
figure()
subplot(231)
RR_ref = grayscale(R_ref[1:NN,1:NN])
imshow(RR_ref, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("reference")
subplot(232)
RR1 = grayscale(R1[1:NN,1:NN])
imshow(RR1, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("TDE")
subplot(233)
RR2 = grayscale(R2[1:NN,1:NN])
imshow(RR2, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("MDOP")
subplot(234)
RR3 = grayscale(R3[1:NN,1:NN])
imshow(RR3, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("GA")
subplot(235)
RR4 = grayscale(R4[1:NN,1:NN])
imshow(RR4, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("PECUZAL")

# figure()
# plot3D(tr[:,1], tr[:,2], tr[:,3])
# title("reference")
# figure()
# plot3D(Y_tde[:,1], Y_tde[:,2], Y_tde[:,3])
# title("TDE")
# figure()
# plot3D(Y_GA[:,1], Y_GA[:,2], Y_GA[:,3])
# title("GA")
# figure()
# plot3D(Y_mdop[:,1], Y_mdop[:,2], Y_mdop[:,3])
# title("MDOP")
# figure()
# plot3D(Y_pec[:,1], Y_pec[:,2], Y_pec[:,3])
# title("PECUZAL")

display("mfnn_tde: $mfnn1")
display("mfnn_mdop: $mfnn2")
display("mfnn_GA: $mfnn3")
display("mfnn_pec: $mfnn4")

display("*****")

display("f_tde: $f1")
display("f_mdop: $f2")
display("f_GA: $f3")
display("f_pec: $f4")

display("*****")

display("DET_ref: $(abs(RQA_ref.DET - RQA_ref.DET)/RQA_ref.DET)")
display("DET_tde: $(abs(RQA_ref.DET - RQA1.DET)/RQA_ref.DET)")
display("DET_mdop: $(abs(RQA_ref.DET - RQA2.DET)/RQA_ref.DET)")
display("DET_GA: $(abs(RQA_ref.DET - RQA3.DET)/RQA_ref.DET)")
display("DET_pec: $(abs(RQA_ref.DET - RQA4.DET)/RQA_ref.DET)")

display("*****")

display("RTE_ref: $(abs(RQA_ref.RTE - RQA_ref.RTE)/RQA_ref.RTE)")
display("RTE_tde: $(abs(RQA_ref.RTE - RQA1.RTE)/RQA_ref.RTE)")
display("RTE_mdop: $(abs(RQA_ref.RTE - RQA2.RTE)/RQA_ref.RTE)")
display("RTE_GA: $(abs(RQA_ref.RTE - RQA3.RTE)/RQA_ref.RTE)")
display("RTE_pec: $(abs(RQA_ref.RTE - RQA4.RTE)/RQA_ref.RTE)")

display("*****")

display("Lmax_ref: $(abs(RQA_ref.Lmax - RQA_ref.Lmax)/RQA_ref.Lmax)")
display("Lmax_tde: $(abs(RQA_ref.Lmax - RQA1.Lmax)/RQA_ref.Lmax)")
display("Lmax_mdop: $(abs(RQA_ref.Lmax - RQA2.Lmax)/RQA_ref.Lmax)")
display("Lmax_GA: $(abs(RQA_ref.Lmax - RQA3.Lmax)/RQA_ref.Lmax)")
display("Lmax_pec: $(abs(RQA_ref.Lmax - RQA4.Lmax)/RQA_ref.Lmax)")

display("*****")

display("ENTR_ref: $(abs(RQA_ref.ENTR - RQA_ref.ENTR)/RQA_ref.ENTR)")
display("ENTR_tde: $(abs(RQA_ref.ENTR - RQA1.ENTR)/RQA_ref.ENTR))")
display("ENTR_mdop: $(abs(RQA_ref.ENTR - RQA2.ENTR)/RQA_ref.ENTR))")
display("ENTR_GA: $(abs(RQA_ref.ENTR - RQA3.ENTR)/RQA_ref.ENTR))")
display("ENTR_pec: $(abs(RQA_ref.ENTR - RQA4.ENTR)/RQA_ref.ENTR))")

display("*****")

display("L_ref: $(abs(RQA_ref.L - RQA_ref.L)/RQA_ref.L)")
display("L_tde: $(abs(RQA_ref.L - RQA1.L)/RQA_ref.L)")
display("L_mdop: $(abs(RQA_ref.L - RQA2.L)/RQA_ref.L)")
display("L_GA: $(abs(RQA_ref.L - RQA3.L)/RQA_ref.L)")
display("L_pec: $(abs(RQA_ref.L - RQA4.L)/RQA_ref.L)")

display("*****")

display("MRT_ref: $(abs(RQA_ref.MRT - RQA_ref.MRT)/RQA_ref.MRT)")
display("MRT_tde: $(abs(RQA_ref.MRT - RQA1.MRT)RQA_ref.MRT)")
display("MRT_mdop: $(abs(RQA_ref.MRT - RQA2.MRT)RQA_ref.MRT)")
display("MRT_GA: $(abs(RQA_ref.MRT - RQA3.MRT)RQA_ref.MRT)")
display("MRT_pec: $(abs(RQA_ref.MRT - RQA4.MRT)RQA_ref.MRT)")

display("*****")

display("LAM_ref: $(abs(RQA_ref.LAM - RQA_ref.LAM)/RQA_ref.LAM)")
display("LAM_tde: $(abs(RQA_ref.LAM - RQA1.LAM)RQA_ref.LAM)")
display("LAM_mdop: $(abs(RQA_ref.LAM - RQA2.LAM)RQA_ref.LAM)")
display("LAM_GA: $(abs(RQA_ref.LAM - RQA3.LAM)RQA_ref.LAM)")
display("LAM_pec: $(abs(RQA_ref.LAM - RQA4.LAM)RQA_ref.LAM)")

display("*****")

display("RTE-MRT_ref: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA_ref.RTE/RQA_ref.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_tde: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA1.RTE/RQA1.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_mdop: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA2.RTE/RQA2.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_GA: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA3.RTE/RQA3.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_pec: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA4.RTE/RQA4.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")

display("*****")

# # using DelimitedFiles
# #
# writedlm("Duffing_Y_ref",tr)
# writedlm("Duffing_Y_tde",Y_tde)
# writedlm("Duffing_Y_mdop",Y_mdop)
# writedlm("Duffing_Y_GA",Y_GA)
# writedlm("Duffing_Y_pec",Y_pec)


text:
function _solution_symmetric_difference(a::Vector{T}, b::Vector{T}) where T
  only_in_a = T[]
  only_in_b = T[]

  for e in a
    if e in b
      continue
    end
    push!(only_in_a, e)
  end

  for e in b
    if e in a
      continue
    end
    push!(only_in_b, e)
  end

  return only_in_a, only_in_b
end

function _solution_symmetric_difference_size(a::Vector{T}, b::Vector{T}) where T
  only_in_a, only_in_b = _solution_symmetric_difference(a, b)
  return length(only_in_a) + length(only_in_b)
end

ext: jl


text:
include("./DC3Df90.jl")

# -----------------------------
# parameters
# -----------------------------
X = -98.0:4.0:98.0
Y = X
nx = length(X)
ny = length(Y)

L1 = -50.0
L2 = 50.0
W1 = -50.0
W2 = 50.0

z = -30.0
depth = 50.0
#dip = 90.0
# -----------------------------

# -----------------------------
# calc
# -----------------------------
Uf90s = zeros(Float64,ny,nx,3)
Uf90d = zeros(Float64,ny,nx,3)
Uf90t = zeros(Float64,ny,nx,3)
for i = 1:ny
for j = 1:nx
    Uf90s[i,j,:] = DC3Df90(X[j], Y[i], z, depth, 90.0, L1, L2, W1, W2, 1e2, 0.0, 0.0)[1:3] # Δu₁ = 100.0 (cm)
    Uf90d[i,j,:] = DC3Df90(X[j], Y[i], z, depth, 45.0, L1, L2, W1, W2, 0.0, 1e2, 0.0)[1:3] # Δu₂ = 100.0 (cm)
    Uf90t[i,j,:] = DC3Df90(X[j], Y[i], z, depth, 90.0, L1, L2, W1, W2, 0.0, 0.0, 1e2)[1:3] # Δu₃ = 100.0 (cm)
end
end
# -----------------------------

# -----------------------------
# Figure
# -----------------------------
using Plots
plotlyjs()
#pyplot()

xp = permutedims(repeat(X, inner=(1,ny)),[2 1]) # km
yp = repeat(Y, inner=(1,nx))

# Note: in x and y coordinates, values of the dislocation u are thousand times magnified

# Strike-slip
u = 1e-2.*Uf90s # cm -> m
x = xp .+ u[:,:,1]
y = yp .+ u[:,:,2]
z = u[:,:,3]
plts = plot(x, y, z, linetype=:surface, c=:balance, xlabel="(km)", ylabel="(km)", zlabel="(m)",
            zlims=(-0.025,0.025), clims=(-0.02,0.02), colorbar=false, guidefont=10, titlefont=font("Times New Roman",16), title="Strike-slip, δ=90°")

# Dip-slip
u = 1e-2.*Uf90d
x = xp .+ u[:,:,1]
y = yp .+ u[:,:,2]
z = u[:,:,3]
pltd = plot(x, y, z, linetype=:surface, c=:balance, xlabel="(km)", ylabel="(km)", zlabel="(m)",
            zlims=(-0.75,0.75), clims=(-0.5,0.5), colorbar=false, guidefont=10, titlefont=font("Times New Roman",16), title="Dip-slip, δ=45°", viewangle=(130,3))

# Tensile
u = 1e-2.*Uf90t
x = xp .+ u[:,:,1]
y = yp .+ u[:,:,2]
z = u[:,:,3]
pltt = plot(x, y, z, linetype=:surface, c=:balance, xlabel="(km)", ylabel="(km)", zlabel="(m)",
            zlims=(0.0,0.20), clims=(0.0,0.15), colorbar=false, guidefont=10, titlefont=font("Times New Roman",16), title="Tensile, δ=90°")

# -----------------------------


text:
using Traitor
using Traitor: extract_arg_trait
if VERSION < v"0.5.0-dev"
    using BaseTestNext
else
    using Base.Test
end

@testset "Parsing of function argument trait constraints" begin

# Simple trait & type symbols
@test extract_arg_trait(:(x::T::A)) == (:(x::T),   :(Tuple{A}))
@test extract_arg_trait(:(x::::A))  == (:(x::Any), :(Tuple{A}))
@test extract_arg_trait(:(x::T))    == (:(x::T),     Tuple{})
@test extract_arg_trait(:(x))       == (:(x::Any),   Tuple{})
# Type exprs
@test extract_arg_trait(:(x::Union{S,T}::A)) == (:(x::Union{S,T}), :(Tuple{A}))
@test extract_arg_trait(:(x::Union{S,T}))    == (:(x::Union{S,T}),   Tuple{})
# Trait exprs
@test extract_arg_trait(:(x::T::Union{A,B})) == (:(x::T),   :(Tuple{Union{A,B}}))
@test extract_arg_trait(:(x::::Union{A,B}))  == (:(x::Any), :(Tuple{Union{A,B}}))
# Trait & type exprs
@test extract_arg_trait(:(x::Union{S,T}::Union{A,B})) == (:(x::Union{S,T}), :(Tuple{Union{A,B}}))
@test extract_arg_trait(:(x::::Union{A,B}))  == (:(x::Any),      :(Tuple{Union{A,B}}))
@test extract_arg_trait(:(x::Union{S,T}))    == (:(x::Union{S,T}), Tuple{})
# Keyword args
@test extract_arg_trait(Expr(:kw,:(x::T::A),1)) == (Expr(:kw,:(x::T),1),   :(Tuple{A}))
@test extract_arg_trait(Expr(:kw,:(x::::A),1))  == (Expr(:kw,:(x::Any),1), :(Tuple{A}))
@test extract_arg_trait(Expr(:kw,:(x::T),1))    == (Expr(:kw,:(x::T),1),     Tuple{})
@test extract_arg_trait(Expr(:kw,:x,1))         == (Expr(:kw,:(x::Any),1),   Tuple{})

end

abstract Fooness

immutable FooA <: Fooness ; end
immutable FooB <: Fooness ; end

Fooness(::Any) = FooA()
Fooness(::Int16) = FooB()

abstract Size

immutable Big <: Size ; end
immutable Small <: Size ; end
immutable Medium <: Size ; end

Size(::Int16) = Small()
Size(::Int) = Small()
Size(::BigInt) = Big()
Size(::Int128) = Medium()

@traitor function f(x::::Big)
   "Huge"
end

@traitor function f(x::::Union{Small,Medium})
   "Smallish"
end

# FIXME!  Defining this breaks the thunk
@traitor function f(x::::Tuple{Small,FooB})
   "A small FooB"
end

@testset "@traitor" begin

@test f(1) == "Smallish"
@test f(Int128(1)) == "Smallish"
@test f(BigInt(1)) == "Huge"
#@test f(Int16(1)) == "A small FooB"

end


text:
using SimpleTypePrint
using Documenter

makedocs(;
    modules=[SimpleTypePrint],
    authors="Jiayi Wei <wjydzh1@gmail.com> and contributors",
    repo="https://github.com/MrVPlusOne/SimpleTypePrint.jl/blob/{commit}{path}#L{line}",
    sitename="SimpleTypePrint.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MrVPlusOne.github.io/SimpleTypePrint.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MrVPlusOne/SimpleTypePrint.jl",
)

max_issues_repo_name: MrVPlusOne/SimpleTypePrint

max_forks_repo_name: MrVPlusOne/SimpleTypePrint


text:
struct GradientData{T} <: AlgorithmCache
    G::T
end

function algorithm_cache(algo::Gradient, prob, net, par_cache)
    T = eltype(trainable_first(net))
    g = grad_cache(T, net)
    return GradientData(g)
end

function setup_algorithm!(g::GradientData, data, par_cache)
    for field=vec_data(g.G)
        T = eltype(field)
        if T<: Real
            field .= real.(data)
        else
            field .= data
        end
    end
end

precondition!(g::GradientData, params::Gradient, iter_n) = g.G


text:
using Carretera
using FactCheck

facts("Probamos las funciones auxiliares Pos") do

  C = Carretera1D(20).carretera

  Cambiar_Vehiculo!(C[5], 5, 5, 2, 0, 1)
  Cambiar_Vehiculo!(C[11], 5, 11, 2, 0, 2)
  Cambiar_Vehiculo!(C[18], 5, 18, 2, 0, 3)
  Cambiar_Vehiculo!(C[15], -1, 15, -2, 0, 0)

  @fact Pos_izq(C) --> 18
  @fact Pos_izq(C[1:17]) --> 11
  @fact Pos_izq(C[1:10]) --> 5
  @fact Pos_izq(C[1:4]) --> 0

  @fact Pos_izq9(C) --> 18
  @fact Pos_izq9(C[1:17]) --> 11
  @fact Pos_izq9(C[1:10]) --> 8

  @fact Pos_der(C) --> 5
  @fact Pos_der(C[6:20]) --> 11
  @fact Pos_der(C[12:20]) --> 18
  @fact Pos_der(C[19:20]) --> 20

  @fact Pos_siguiente(C) --> 20
  @fact Pos_siguiente(C[1:18]) --> 17
end

facts("Probamos las funciones auxiliares Vehiculo") do

  C = Carretera1D(20).carretera

  Cambiar_Vehiculo!(C[5], 5, 5, 2, 0, 1)
  Cambiar_Vehiculo!(C[11], 5, 11, 2, 0, 2)
  Cambiar_Vehiculo!(C[18], 5, 18, 2, 0, 3)
  Cambiar_Vehiculo!(C[15], -1, 15, -2, 0, 0)

  @fact Vehiculo_Adelante(12, C, 20) --> 15
  @fact Vehiculo_Adelante(19, C, 20) --> 20

  @fact Vehiculo_Atras(11, C, 20) --> 20
  @fact Vehiculo_Atras(18, C, 20) --> 11
end

#fact("Probamos las funciones básicas") do

 #  C = Carretera1D(20).carretera


text:
export check_fileheader

function check_fileheader(s::IO, fh_new::FileHeader)

    # Check text header length
    if sizeof(fh_new.th) > 3200 @error "Text Header longer than 3200 bytes" end

    # Check dsf
    if fh_new.bfh.DataSampleFormat != 5
        @warn "DataSampleFormat not supported for writing. Attempting to convert to IEEE Float32"
        fh_new.bfh.DataSampleFormat = 5
    end

    # read fileheader of the existing file
    fh = read_fileheader(s)

    # Check first section of assigned values

    for field in fieldnames(typeof(fh.bfh))[1:end]
        if(getfield(fh.bfh, field) != getfield(fh_new.bfh, field)) @error "The existing fileheaders are different from the new SeisBlock" end
    end

end


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

# Collection of sources required to build MPCBuilder
name = "FLINT"
hash = "ea9f102a0bb8ab5fd50258aa177c345593ab11df"
version = VersionNumber("2.5.3+dev0.$(hash[1:10])")
sources = [
    "https://github.com/wbhart/flint2/archive/$hash.tar.gz" =>
    "3719310c21ebb2c6c6144f571026e4e1558847eaba7464c707484a4c33ca7356",
]

# Bash recipe for building across all platforms

script = """cd flint2-$hash
"""*raw"""
./configure --prefix=$prefix --enable-shared --disable-static --with-gmp=$prefix --with-mpfr=$prefix
make -j
make install
"""
# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:aarch64, :glibc, :blank_abi),
    BinaryProvider.Windows(:i686, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:armv7l, :glibc, :eabihf),
    BinaryProvider.Windows(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:x86_64, :glibc, :blank_abi),
    BinaryProvider.MacOS(:x86_64, :blank_libc, :blank_abi),
    BinaryProvider.Linux(:i686, :glibc, :blank_abi)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libflint", :libflint)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl"
    "https://github.com/JuliaMath/MPFRBuilder/releases/download/v4.0.1/build.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: isuruf/FLINTBuilder

max_issues_repo_name: isuruf/FLINTBuilder

max_forks_repo_name: isuruf/FLINTBuilder


text:
module ModelingToolkit

using DiffEqBase, Distributed
using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays
using Latexify, Unitful, ArrayInterface
using MacroTools
using UnPack: @unpack
using DiffEqJump
using DataStructures: OrderedDict, OrderedSet
using SpecialFunctions, NaNMath
using RuntimeGeneratedFunctions
using Base.Threads
import MacroTools: splitdef, combinedef, postwalk, striplines
import Libdl
using DocStringExtensions
using Base: RefValue
import IfElse

RuntimeGeneratedFunctions.init(@__MODULE__)

using RecursiveArrayTools

import SymbolicUtils
import SymbolicUtils: Term, Sym, to_symbolic, FnType, @rule, Rewriters, substitute, similarterm

using LinearAlgebra: LU, BlasInt

import LightGraphs: SimpleDiGraph, add_edge!

import TreeViews

using Requires

export Num, Variable
"""
$(TYPEDEF)

Wrap anything in a type that is a subtype of Real
"""
struct Num <: Real
    val
end

const show_numwrap = Ref(false)

Num(x::Num) = x # ideally this should never be called
(n::Num)(args...) = Num(value(n)(map(value,args)...))
value(x) = x
value(x::Num) = x.val

using SymbolicUtils: to_symbolic
SymbolicUtils.to_symbolic(n::Num) = value(n)
SymbolicUtils.@number_methods(Num,
                              Num(f(value(a))),
                              Num(f(value(a), value(b))))

SymbolicUtils.simplify(n::Num; kw...) = Num(SymbolicUtils.simplify(value(n); kw...))

SymbolicUtils.symtype(n::Num) = symtype(n.val)

function Base.iszero(x::Num)
    _x = SymbolicUtils.to_mpoly(value(x))[1]
    return (_x isa Number || _x isa SymbolicUtils.MPoly) && iszero(_x)
end

import SymbolicUtils: <ₑ, Symbolic, Term, operation, arguments

Base.show(io::IO, n::Num) = show_numwrap[] ? print(io, :(Num($(value(n))))) : Base.show(io, value(n))

Base.promote_rule(::Type{<:Number}, ::Type{<:Num}) = Num
Base.promote_rule(::Type{<:Symbolic{<:Number}}, ::Type{<:Num}) = Num
Base.getproperty(t::Term, f::Symbol) = f === :op ? operation(t) : f === :args ? arguments(t) : getfield(t, f)
<ₑ(s::Num, x) = value(s) <ₑ value(x)
<ₑ(s, x::Num) = value(s) <ₑ value(x)
<ₑ(s::Num, x::Num) = value(s) <ₑ value(x)

for T in (Integer, Rational)
    @eval Base.:(^)(n::Num, i::$T) = Num(Term{symtype(n)}(^, [value(n),i]))
end

macro num_method(f, expr, Ts=nothing)
    if Ts === nothing
        Ts = [Any]
    else
        @assert Ts.head == :tuple
        # e.g. a tuple or vector
        Ts = Ts.args
    end

    ms = [quote
              $f(a::$T, b::$Num) = $expr
              $f(a::$Num, b::$T) = $expr
          end for T in Ts]
    quote
        $f(a::$Num, b::$Num) = $expr
        $(ms...)
    end |> esc
end

"""
    tosymbolic(a::Union{Sym,Num}) -> Sym{Real}
    tosymbolic(a::T) -> T
"""
tosymbolic(a::Num) = tosymbolic(value(a))
tosymbolic(a::Sym) = tovar(a)
tosymbolic(a) = a
@num_method Base.isless isless(tosymbolic(a), tosymbolic(b)) (Real,)
@num_method Base.:(<) (tosymbolic(a) < tosymbolic(b)) (Real,)
@num_method Base.:(<=) (tosymbolic(a) <= tosymbolic(b)) (Real,)
@num_method Base.:(>) (tosymbolic(a) > tosymbolic(b)) (Real,)
@num_method Base.:(>=) (tosymbolic(a) >= tosymbolic(b)) (Real,)
@num_method Base.isequal isequal(tosymbolic(a), tosymbolic(b)) (Number, Symbolic)
@num_method Base.:(==) tosymbolic(a) == tosymbolic(b) (Number,)

Base.hash(x::Num, h::UInt) = hash(value(x), h)

Base.convert(::Type{Num}, x::Symbolic{<:Number}) = Num(x)
Base.convert(::Type{Num}, x::Number) = Num(x)
Base.convert(::Type{Num}, x::Num) = x

Base.convert(::Type{<:Array{Num}}, x::AbstractArray) = map(Num, x)
Base.convert(::Type{<:Array{Num}}, x::AbstractArray{Num}) = x
Base.convert(::Type{Sym}, x::Num) = value(x) isa Sym ? value(x) : error("cannot convert $x to Sym")

LinearAlgebra.lu(x::Array{Num}; kw...) = lu(x, Val{false}(); kw...)

"""
$(TYPEDEF)

TODO
"""
abstract type AbstractSystem end
abstract type AbstractODESystem <: AbstractSystem end

"""
$(TYPEDSIGNATURES)

Get the set of independent variables for the given system.
"""
function independent_variables end

"""
$(TYPEDSIGNATURES)

Get the set of states for the given system.
"""
function states end

"""
$(TYPEDSIGNATURES)

Get the set of parameters variables for the given system.
"""
function parameters end

include("variables.jl")
include("context_dsl.jl")
include("differentials.jl")

include("equations.jl")
include("utils.jl")
include("linearity.jl")
include("solve.jl")
include("direct.jl")
include("domains.jl")
include("register_function.jl")

include("systems/abstractsystem.jl")

include("systems/diffeqs/odesystem.jl")
include("systems/diffeqs/sdesystem.jl")
include("systems/diffeqs/abstractodesystem.jl")
include("systems/diffeqs/first_order_transform.jl")
include("systems/diffeqs/modelingtoolkitize.jl")
include("systems/diffeqs/validation.jl")

include("systems/jumps/jumpsystem.jl")

include("systems/nonlinear/nonlinearsystem.jl")

include("systems/optimization/optimizationsystem.jl")

include("systems/control/controlsystem.jl")

include("systems/pde/pdesystem.jl")

include("systems/reaction/reactionsystem.jl")
include("systems/dependency_graphs.jl")

include("systems/reduction.jl")

include("latexify_recipes.jl")
include("build_function.jl")
include("extra_functions.jl")

export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
export JumpSystem
export ODEProblem, SDEProblem
export NonlinearProblem, NonlinearProblemExpr
export OptimizationProblem, OptimizationProblemExpr
export SteadyStateProblem, SteadyStateProblemExpr
export JumpProblem, DiscreteProblem
export NonlinearSystem, OptimizationSystem
export ControlSystem
export ode_order_lowering
export runge_kutta_discretize
export PDESystem
export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw
export Differential, expand_derivatives, @derivatives
export IntervalDomain, ProductDomain, ⊗, CircleDomain
export Equation, ConstrainedEquation
export Term, Sym
export independent_variable, states, parameters, equations, controls, pins, observed

export calculate_jacobian, generate_jacobian, generate_function
export calculate_tgrad, generate_tgrad
export calculate_gradient, generate_gradient
export calculate_factorized_W, generate_factorized_W
export calculate_hessian, generate_hessian
export calculate_massmatrix, generate_diffusion_function
export stochastic_integral_transform

export BipartiteGraph, equation_dependencies, variable_dependencies
export eqeq_dependencies, varvar_dependencies
export asgraph, asdigraph

export toexpr, rename, get_variables
export simplify, substitute
export build_function
export @register
export modelingtoolkitize
export @variables, @parameters

const HAS_DAGGER = Ref{Bool}(false)
function __init__()
    @require Dagger="d58978e5-989f-55fb-8d15-ea34adc7bf54" include("dagger.jl")
end

end # module


text:
function is_file_path_ok(path_to_file::String)::VLResult

    if (isfile(path_to_file) == false)

        # error message -
        error_message = "Ooops! $(path_to_file) does not exist."
        error_object = ArgumentError(error_message)

        # return -
        return VLResult{ArgumentError}(error_object)
    end

    # default: return nothing -
    return VLResult(nothing)
end

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
using Cropbox
using GLMakie
using Test

#### CropRootBox Start
using Cropbox
using Distributions
import Makie
import Meshing
using GeometryBasics: GeometryBasics, Mesh, Point3f, coordinates, faces, FaceView
using CoordinateTransformations: IdentityTransformation, LinearMap, Transformation, Translation, recenter
using Rotations: RotZX, RotXY, RotX, RotY
using Colors: RGBA
using StaticArrays
import UUIDs

@system Rendering

@system Container(Rendering, Controller) begin
    dist(; p::Point3f): distance => -Inf ~ call
end

@system Pot(Container) <: Container begin
    r1: top_radius => 10 ~ preserve(u"cm", parameter)
    r2: bottom_radius => 6 ~ preserve(u"cm", parameter)
    h: height => 30 ~ preserve(u"cm", parameter)
    sq: square => false ~ preserve::Bool(parameter)

    dist(nounit(r1), nounit(r2), nounit(h), sq; p::Point3f): distance => begin
        x, y, z = p
        if z < -h # below
            -z - h
        elseif 0 < z # above
            z
        else # inside: -h <= z <= 0
            w = -z / h # [0, 1]
            r = (1-w)*r1 + w*r2
            if sq
                max(abs(x), abs(y)) - r
            else
                sqrt(x^2 + y^2) - r
            end
        end
    end ~ call
end

mesh(s::Pot) = begin
    r1 = Cropbox.deunitfy(s.r1', u"cm")
    r2 = Cropbox.deunitfy(s.r2', u"cm")
    h = Cropbox.deunitfy(s.h', u"cm")
    GeometryBasics.Mesh(x -> s.dist'(x), GeometryBasics.Rect(GeometryBasics.Vec(-2r1, -2r2, -1.5h), GeometryBasics.Vec(4r1, 4r2, 3h)), Meshing.MarchingCubes(), samples=(50, 50, 50))
end

@system Rhizobox(Container) <: Container begin
    l: length => 16u"inch" ~ preserve(u"cm", parameter)
    w: width => 10.5u"inch" ~ preserve(u"cm", parameter)
    h: height => 42u"inch" ~ preserve(u"cm", parameter)

    dist(nounit(l), nounit(w), nounit(h); p::Point3f): distance => begin
        x, y, z = p
        if z < -h # below
            -z - h
        elseif 0 < z # above
            z
        else # inside: -h <= z <= 0
            d = abs(y) - w/2
            d < 0 ? abs(x) - l/2 : d
        end
    end ~ call
end

mesh(s::Rhizobox) = begin
    l = Cropbox.deunitfy(s.l', u"cm")
    w = Cropbox.deunitfy(s.w', u"cm")
    h = Cropbox.deunitfy(s.h', u"cm")
    g = GeometryBasics.Rect3(Point3f(-l/2, -w/2, 0), Point3f(l, w, -h))
    GeometryBasics.mesh(g)
end

@system Rhizobox2(Container) <: Container begin
    l: length => 16u"inch" ~ preserve(u"cm", parameter)
    w: width => 10.5u"inch" ~ preserve(u"cm", parameter)
    h: height => 42u"inch" ~ preserve(u"cm", parameter)

    θ_w: ground_incident_w => 60 ~ preserve(u"°", parameter)
    ϕ_w(θ_w): calc_angle_w => 90u"°" - θ_w ~ preserve(u"°")

    θ_l: ground_incident_l => 90 ~ preserve(u"°", parameter)
    ϕ_l(θ_l): calc_angle_l => 90u"°" - θ_l ~ preserve(u"°")

    w_scale: width_position => 0 ~ preserve(parameter, extern, min = 0, max = 1) # Scale of 0 -> 1 for position from mid to edge
    l_scale: length_position => 0 ~ preserve(parameter, extern, min = 0, max = 1) # Scale of 0 -> 1 for position from mid to edge

    RT(nounit(w), nounit(l), w_scale, l_scale, ϕ_w, ϕ_l): transformation => begin
        R = RotXY(ϕ_l, ϕ_w) |> LinearMap
        T = Translation(l_scale * -l/2 * cos(ϕ_w), w_scale * -w/2 * cos(ϕ_l), 0)
        R ∘ T
    end ~ track::Transformation

    dist(nounit(l), nounit(w), nounit(h), RT, ϕ_w, ϕ_l; p::Point3f): distance => begin
        x, y, z = inv(RT)(p)
        if z < -h # below
            -z - h
        elseif 0 < z # above
            z
        else # inside: -h <= z <= 0
            d = abs(y) - w/2
            d < 0 ? abs(x) - l/2 : d
        end
    end ~ call
end

mesh(s::Rhizobox2) = begin
    l = Cropbox.deunitfy(s.l', u"cm")
    w = Cropbox.deunitfy(s.w', u"cm")
    h = Cropbox.deunitfy(s.h', u"cm")
    g = GeometryBasics.Rect3(Point3f(-l/2, -w/2, 0), Point3f(l, w, -h))
    m = GeometryBasics.mesh(g)

    mv = GeometryBasics.decompose(GeometryBasics.Point, m) # Decompose the mesh into coordinates
    mf = GeometryBasics.decompose(GeometryBasics.TriangleFace{Int}, m) # Decompose mesh into faces
    M = s.RT'
    GeometryBasics.normal_mesh(M.(mv), mf) # Apply rotation to coordinates and connect back together using faces
end

@system SoilCore(Container) <: Container begin
    d: diameter => 5 ~ preserve(u"cm", parameter)
    l: length => 90 ~ preserve(u"cm", parameter)
    x0: x_origin => 0 ~ preserve(u"cm", parameter)
    y0: y_origin => 0 ~ preserve(u"cm", parameter)

    dist(nounit(d), nounit(l), nounit(x0), nounit(y0); p::Point3f): distance => begin
        x, y, z = p
        if z < -l # below
            -z - l
        elseif 0 < z # above
            z
        else # inside: -l <= z <= 0
            sqrt((x - x0)^2 + (y - y0)^2) - d/2
        end
    end ~ call
end

mesh(s::SoilCore) = begin
    d = Cropbox.deunitfy(s.d', u"cm")
    l = Cropbox.deunitfy(s.l', u"cm")
    x0 = Cropbox.deunitfy(s.x0', u"cm")
    y0 = Cropbox.deunitfy(s.y0', u"cm")
    g = GeometryBasics.Cylinder(Point3f(x0, y0, 0), Point3f(x0, y0, -l), Float32(d)/2)
    GeometryBasics.mesh(g)
end

@system SoilLayer(Container) <: Container begin
    d: depth => 0 ~ preserve(u"cm", parameter)
    t: thickness => 10 ~ preserve(u"cm", parameter)

    dist(nounit(d), nounit(t); p::Point3f): distance => begin
        x, y, z = p
        a = -dc
        b = a - t
        if a <= z # above
            z - a
        elseif z < b # below
            -(z - b)
        else # inside: b <= z < a
            z
        end
    end ~ call
end

@system Tropism begin
    N: tropsim_trials => 1.0 ~ preserve(parameter)
    to(; α, β): tropism_objective => 0 ~ call
end

@system Plagiotropism(Tropism) <: Tropism begin
    RT0: parent_transformation ~ hold
    to(RT0; α, β): tropism_objective => begin
        R = RotZX(β, α) |> LinearMap
        (RT0 ∘ R).linear[9] |> abs
    end ~ call
end

@system Gravitropism(Tropism) <: Tropism begin
    RT0: parent_transformation ~ hold
    to(RT0; α, β): tropism_objective => begin
        R = RotZX(β, α) |> LinearMap
        #-(RT0 ∘ R).linear[9]
        p = (RT0 ∘ R)(Point3f(0, 0, -1))
        p[3]
    end ~ call
end

@system Exotropism(Tropism) <: Tropism begin
    to(; α, β): tropism_objective => begin
        #HACK: not exact implementation, needs to keep initial heading
        abs(Cropbox.deunitfy(α))
    end ~ call
end

@system RootSegment{ # TODO: Add diameter
    Segment => RootSegment,
    Branch => RootSegment,
}(Tropism, Rendering) begin
    box ~ <:Container(override)

    ro: root_order => 1 ~ preserve::int(extern)
    zi: zone_index => 0 ~ preserve::int(extern)

    zt(lmax, la, lb, lp, ln): zone_type => begin
        if (lmax - la) <= lp
            :apical
        elseif lp < lb
            :basal
        else
            #HACK: assume apical when no lateral branching zone exists
            iszero(ln) ? :apical : :lateral
        end
    end ~ preserve::sym

    lb: length_of_basal_zone => 0.4 ~ preserve(u"cm", extern, parameter, min=0)
    la: length_of_apical_zone => 0.5 ~ preserve(u"cm", extern, parameter, min=0)
    ln: length_between_lateral_branches => 0.3 ~ preserve(u"cm", extern, parameter, min=0)
    lmax: maximal_length => 3.9 ~ preserve(u"cm", extern, parameter, min=0.1)

    r: maximum_elongation_rate => 1.0 ~ preserve(u"cm/d", extern, parameter, min=0)
    pr(r): potential_elongation_rate ~ track(u"cm/d")

    t(context.clock.time): timestamp ~ preserve(u"minute")
    Δl(Δx) ~ preserve(u"cm", max=lr)
    lp: parent_length => 0 ~ preserve(u"cm", extern)
    ls: sibling_length => 0 ~ preserve(u"cm", extern)
    l(pr): length ~ accumulate(u"cm", max=Δl, when=!im)
    lr(lp, lmax): remaining_length => lmax - lp ~ track(u"cm")
    lt(lp, l): total_length => lp + l ~ track(u"cm")
    li(ls, l): interleaving_length => ls + l ~ track(u"cm")
    ls1(li): next_sibling_length ~ track(u"cm", when=!ib)

    Δx: axial_resolution => 1 ~ preserve(u"cm", parameter)
    σ: standard_deviation_of_angle => 30 ~ preserve(u"°", parameter)
    σ_Δx(σ, nounit(Δx)): normalized_standard_deviation_of_angle => sqrt(Δx)*σ ~ preserve(u"°")

    θ: insertion_angle => 30 ~ preserve(u"°", parameter)
    pα(zi, nounit(θ), nounit(σ_Δx);): pick_angular_angle => begin
        θ = zi == 0 ? θ : zero(θ)
        rand(Normal(θ, σ_Δx))
    end ~ call(u"°")
    pβ(;): pick_radial_angle => rand(Uniform(0, 360)) ~ call(u"°")
    αN: angular_angle_trials => 20 ~ preserve::int(parameter)
    βN: radial_angle_trials => 5 ~ preserve::int(parameter)
    A(pα, pβ, to, N, dist=box.dist, np, αN, βN): angles => begin
        n = rand() < N % 1 ? ceil(N) : floor(N)
        P = [(pα(), pβ()) for i in 0:n]
        O = [to(α, β) for (α, β) in P]
        (o, i) = findmin(O)
        (α, β) = P[i]
        d = dist(np(α, β))
        for i in 1:αN
            #HACK: look around 360° rather than original 90° in case the segment already is out of boundary
            α1 = α + 360u"°" * (i-1)/αN
            for j in 1:βN
                d < 0 && break
                β1 = pβ()
                d1 = dist(np(α1, β1))
                if d1 < d
                    d = d1
                    α, β = α1, β1
                end
            end
            d < 0 && break
        end
        (α, β)
    end ~ preserve::Tuple
    α(A): angular_angle => A[1] ~ preserve(u"°")
    β(A): radial_angle => A[2] ~ preserve(u"°")

    RT0: parent_transformation ~ track::Transformation(override)
    pp(RT0): parent_position => RT0(Point3f(0, 0, 0)) ~ preserve::Point3f
    np(RT0, nounit(Δl); α, β): new_position => begin
        R = RotZX(β, α) |> LinearMap
        (RT0 ∘ R)(Point3f(0, 0, -Δl))
    end ~ call::Point3f
    RT(nounit(l), α, β): local_transformation => begin
        # put root segment at parent's end
        T = Translation(0, 0, -l)
        # rotate root segment
        R = RotZX(β, α) |> LinearMap
        R ∘ T
    end ~ track::Transformation
    RT1(RT0, RT): global_transformation => RT0 ∘ RT ~ track::Transformation
    cp(RT1): current_position => RT1(Point3f(0, 0, 0)) ~ track::Point3f

    # Attempt at dynamic radius
    # ISSUE: Radius growth is happening backwards, thicker at the bottom
    # SOLUTION: Somehow trace back to previous parent roots and update based on time step?
    # Very intensive in terms of processing/time
    # OR: Could hard code a simulation stop time into the config that's passed here. Can then do
    # (stop time - t) as the time calculation to get diameter to work
    ri: radius_initial => 0.04 ~ preserve(u"cm", extern, parameter, min = 0.01)
    thr: radius_threshold => 0.1 ~ preserve(u"cm", extern, parameter, min = 0.01)
    gr: radius_growth_rate => 0.0005 ~ preserve(u"mm/hr", extern, parameter, min = 0.0000001)
    fl(a, thr): flag_variable => a < thr ~ flag
    a(gr): radius ~ accumulate(u"mm", init = ri, when = fl, min=0.01)

    c: color => RGBA(1, 1, 1, 1) ~ preserve::RGBA(parameter)

    n: name ~ hold
    T: transition ~ hold
    nb(T, name;): next_branch => begin
        find(r) = begin
            d = T[name]
            for (k, v) in d
                r < v ? (return k) : (r -= v)
            end
            :nothing
        end
        find(rand())
    end ~ call::sym

    ms(l, Δl, lt, lmax): may_segment => (l >= Δl && lt < lmax) ~ flag
    S(n, box, ro, zi, r, lb, la, ln, lmax, lt, ls1, wrap(RT1), a): segment => begin
        #HACK: keep lb/la/ln/lmax parameters same for consecutive segments
        produce(eval(n); box, ro, zi=zi+1, r, lb, la, ln, lmax, lp=lt, ls=ls1, RT0=RT1, a)
    end ~ produce::Segment(when=ms)

    mb(zt, li, ln): may_branch => (zt == :lateral && li >= ln) ~ flag
    B(nb, box, ro, wrap(RT1)): branch => begin
        #HACK: eval() for Symbol-based instantiation based on tabulate-d matrix
        produce(eval(nb()); box, ro=ro+1, RT0=RT1)
    end ~ produce::Branch(when=mb)

    ii(cp; c::Container): is_inside => (c.dist'(cp) <= 0) ~ call::Bool
    im(l, Δl, lt, lmax): is_mature => (l >= Δl || lt >= lmax) ~ flag

    is(S): is_segmented => !isnothing(S) ~ flag
    ib(B): is_branched => !isnothing(B) ~ flag
end

mesh(s::RootSegment) = begin
    l = Cropbox.deunitfy(s.l', u"cm")
    a = Cropbox.deunitfy(s.a', u"cm")
    (iszero(l) || iszero(a)) && return nothing
    r = a/2
    g = GeometryBasics.Rect3(Point3f(-r, -r, 0), Point3f(a, a, l))
    m = GeometryBasics.mesh(g)

    #HACK: reconstruct a mesh with transformation applied
    mv = GeometryBasics.decompose(GeometryBasics.Point, m)
    mf = GeometryBasics.decompose(GeometryBasics.TriangleFace{Int}, m)
    M = s.RT1'
    m = GeometryBasics.normal_mesh(M.(mv), mf)

    c = s.color'
    n = length(GeometryBasics.coordinates(m))
    GeometryBasics.pointmeta(m; color=fill(c, n))
end

#TODO: provide @macro / function to automatically build a series of related Systems
@system BaseRoot(RootSegment) <: RootSegment begin
    T: transition ~ tabulate(rows=(:PrimaryRoot, :FirstOrderLateralRoot, :SecondOrderLateralRoot, #= :ThirdOrderLateralRoot =#), parameter)
end
@system ThirdOrderLateralRoot{
    Segment => ThirdOrderLateralRoot,
}(BaseRoot, Gravitropism) <: BaseRoot begin
    n: name => :ThirdOrderLateralRoot ~ preserve::sym
end
@system SecondOrderLateralRoot{
    Segment => SecondOrderLateralRoot,
    #= Branch => ThirdOrderLateralRoot, =#
}(BaseRoot, Gravitropism) <: BaseRoot begin
    n: name => :SecondOrderLateralRoot ~ preserve::sym
end
@system FirstOrderLateralRoot{
    Segment => FirstOrderLateralRoot,
    Branch => SecondOrderLateralRoot,
}(BaseRoot, Gravitropism) <: BaseRoot begin
    n: name => :FirstOrderLateralRoot ~ preserve::sym
end
@system PrimaryRoot{
    Segment => PrimaryRoot,
    Branch => FirstOrderLateralRoot,
}(BaseRoot, Gravitropism) <: BaseRoot begin
    n: name => :PrimaryRoot ~ preserve::sym
end

@system RootArchitecture(Controller) begin
    box(context) ~ <:Container(override)
    minB: minimum_number_of_basal_roots => 1 ~ preserve(parameter)
    maxB: number_of_basal_roots => 1 ~ preserve(parameter, min=minB)
    RT0: initial_transformation => IdentityTransformation() ~ track::Transformation
    roots(roots, box, maxB, wrap(RT0)) => begin
        [produce(PrimaryRoot; box, RT0) for i in (length(roots)+1):maxB]
    end ~ produce::PrimaryRoot[]
end

render(s::RootArchitecture; soilcore=nothing, resolution=(500, 500)) = begin
    #HACK: comfortable default size when using WGLMakie inside Jupyter Notebook
    scene = Makie.Scene(; resolution)
    Makie.mesh!(scene, mesh(s))
    #HACK: customization for container
    Makie.mesh!(scene, mesh(s.box), color=(:black, 0.02), transparency=true, shading=false)
    !isnothing(soilcore) && Makie.mesh!(scene, mesh(soilcore), color=(:purple, 0.1), transparency=true, shading=false)
    scene
end

mesh(s::RootArchitecture; container=nothing) = begin
    meshes = GeometryBasics.Mesh[]
    gather!(s, Rendering; store=meshes, callback=render!, kwargs=(; container))
    isempty(meshes) ? nothing : merge(meshes)
end
render!(g::Gather, r::RootSegment, ::Val{:Rendering}; container=nothing) = begin
    m = isnothing(container) || r.ii'(container) ? mesh(r) : nothing
    !isnothing(m) && push!(g, m)
    visit!(g, r; container)
end
render!(g::Gather, a...; k...) = visit!(g, a...; k...)

gatherbaseroot!(g::Gather, s::BaseRoot, ::Val{:BaseRoot}) = (push!(g, s); visit!(g, s))
gatherbaseroot!(g::Gather, a...) = visit!(g, a...)

#TODO: implement a simpler generic gather interface
gather_primaryroot!(g::Gather, s::PrimaryRoot, ::Val{:BaseRoot}) = (push!(g, s); visit!(g, s))
gather_primaryroot!(g::Gather, s::BaseRoot, ::Val{:BaseRoot}) = visit!(g, s)
gather_primaryroot!(g::Gather, a...) = visit!(g, a...)

gather_firstorderlateralroot!(g::Gather, s::FirstOrderLateralRoot, ::Val{:BaseRoot}) = (push!(g, s); visit!(g, s))
gather_firstorderlateralroot!(g::Gather, s::BaseRoot, ::Val{:BaseRoot}) = visit!(g, s)
gather_firstorderlateralroot!(g::Gather, a...) = visit!(g, a...)

gather_secondorderlateralroot!(g::Gather, s::SecondOrderLateralRoot, ::Val{:BaseRoot}) = (push!(g, s); visit!(g, s))
gather_secondorderlateralroot!(g::Gather, s::BaseRoot, ::Val{:BaseRoot}) = visit!(g, s)
gather_secondorderlateralroot!(g::Gather, a...) = visit!(g, a...)

make_gather_root(::Type{PrimaryRoot}) = gather_primaryroot!
make_gather_root(::Type{FirstOrderLateralRoot}) = gather_firstorderlateralroot!
make_gather_root(::Type{SecondOrderLateralRoot}) = gather_secondorderlateralroot!
make_gather_root(::Type{BaseRoot}) = gatherbaseroot!

gathergeom!(g::Gather,  r::BaseRoot, ::Val{:BaseRoot}) = begin
    r.zi' == 0 && push!(g, (
        point=[r.pp', r.cp', r.S["**"].cp'...],
        radius=[r.a', r.a', r.S["**"].a'...] |> Cropbox.deunitfy,
        timestamp=[r.t', r.t', r.S["**"].t'...] |> Cropbox.deunitfy,
    ))
    visit!(g, r)
end
gathergeom!(g::Gather, a...) = visit!(g, a...)

using WriteVTK
gathervtk(name::AbstractString, s::System) = begin
    L = gather!(s, BaseRoot; callback=gathergeom!)
    P = Float32[]
    C = MeshCell[]
    i = 0
    for l in L
        lp = l.point
        [append!(P, p) for p in lp]
        n = length(lp)
        I = collect(1:n) .+ i
        i += n
        c = MeshCell(VTKCellTypes.VTK_POLY_LINE, I)
        push!(C, c)
    end
    P3 = reshape(P, 3, :)
    g = vtk_grid(name, P3, C)
    for k in (:radius, :timestamp)
        D = [l[k] for l in L] |> Iterators.flatten |> collect
        isempty(D) && continue
        g[string(k), VTKPointData()] = D
    end
    g
end
writevtk(name::AbstractString, s::System) = vtk_save(gathervtk(name, s))
writepvd(name::AbstractString, S::Type{<:System}; kwargs...) = begin
    pvd = paraview_collection(name)
    path = mkpath("$name-pvd")
    i = 0
    simulate(S; kwargs...) do D, s
        pvd[i] = gathervtk("$path/$name-$i", s)
        i += 1
    end
    vtk_save(pvd)
end

using MeshIO
using FileIO
writestl(name::AbstractString, s::System) = writestl(name, mesh(s))
writestl(name::AbstractString, m::Mesh) = save(File{format"STL_BINARY"}(name), m)
writestl(name::AbstractString, ::Nothing) = @warn "no mesh available for writing $name"

#### CropRootBox End

# CONFIGS

### Figure 1-4 root config
root_maize = @config(
    :RootArchitecture => :maxB => 5,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 1 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.1 ± 0.01,
        la = 18.0 ± 1.8,
        ln = 0.6 ± 0.06,
        lmax = 89.7 ± 7.4,
        r = 6.0 ± 0.6,
        Δx = 0.5,
        σ = 10,
        θ = 80 ± 8,
        N = 1.5,
        ri = 0.04 ± 0.004,
        gr = 0.005,
        thr = 0.1,
        color = RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.2 ± 0.04,
        la = 0.4 ± 0.04,
        ln = 0.4 ± 0.03,
        lmax = 0.6 ± 1.6,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 15,
        N = 1,
        ri = 0.03 ± 0.003,
        gr = 0.003,
        thr = 0.07,
        color = RGBA(0, 1, 0, 1),
    ),
    :SecondOrderLateralRoot => (;
        lb = 0,
        la = 0.4 ± 0.02,
        ln = 0,
        lmax = 0.4,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 10,
        N = 2,
        ri = 0.02 ± 0.002,
        gr = 0.001,
        thr = 0.03,
        color = RGBA(0, 0, 1, 1),
    )
)

# Figure 1
container_rhizobox = :Rhizobox => (;
    l = 16u"inch",
    w = 10.5u"inch",
    h = 42u"inch",
)

# Figure 2
container_rhizobox = :Rhizobox2 => (;
    l = 5u"inch",
    w = 10u"inch",
    h = 40u"inch",
    θ_w = 70u"°",
    θ_l = 90u"°",
    w_scale = 0,
    l_scale = 0,
)

# Figure 3
container_rhizobox = :Rhizobox2 => (;
    l = 5u"inch",
    w = 10u"inch",
    h = 40u"inch",
    θ_w = 70u"°",
    θ_l = 90u"°",
    w_scale = 0.25,
    l_scale = 1,
)

# Figure 4
container_rhizobox = :Rhizobox2 => (;
    l = 5u"inch",
    w = 10u"inch",
    h = 40u"inch",
    θ_w = 70u"°",
    θ_l = 80u"°",
    w_scale = 0,
    l_scale = 0,
)

# Figure 5
root_maize = @config(
    :RootArchitecture => :maxB => 5,
    :BaseRoot => :T => [
        # P F S
          0 1 0 ; # P
          0 0 1 ; # F
          0 0 0 ; # S
    ],
    :PrimaryRoot => (;
        lb = 0.1 ± 0.01,
        la = 18.0 ± 1.8,
        ln = 0.6 ± 0.06,
        lmax = 89.7 ± 7.4,
        r = 6.0 ± 0.6,
        Δx = 0.5,
        σ = 10,
        θ = 80 ± 8,
        N = 1.5,
        ri = 0.04 ± 0.004,
        gr = 0.005,
        thr = 1.0,
        color = RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.2 ± 0.04,
        la = 0.4 ± 0.04,
        ln = 0.4 ± 0.03,
        lmax = 0.6 ± 1.6,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 15,
        N = 1,
        ri = 0.03 ± 0.003,
        gr = 0.003,
        thr = 0.7,
        color = RGBA(0, 1, 0, 1),
    ),
    :SecondOrderLateralRoot => (;
        lb = 0,
        la = 0.4 ± 0.02,
        ln = 0,
        lmax = 0.4,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 10,
        N = 2,
        ri = 0.02 ± 0.002,
        gr = 0.001,
        thr = 0.3,
        color = RGBA(0, 0, 1, 1),
    )
)

container_rhizobox = :Rhizobox2 => (;
    l = 5u"inch",
    w = 10u"inch",
    h = 40u"inch",
    θ_w = 70u"°",
    θ_l = 90u"°",
    w_scale = 0,
    l_scale = 0,
)

# SIMULATE AND RUN CONFIGS
b = instance(Rhizobox2, config = container_rhizobox)
s = instance(RootArchitecture; config = root_maize, options = (; box = b), seed = 0)
r = simulate!(s, stop = 100u"d") #(to see diameter effect, reduce simulation length to 10d)

scn = render(s)


text:
@testset "indicator" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, x, CS.Integers([-3, 1, 2, 3]))
    @variable(m, y, CS.Integers([-3, 1, 2, 3]))
    @variable(m, b, Bin)
    @constraint(m, b => {x + y + 1 == 5})
    optimize!(m)
    com = CS.get_inner_model(m)
    constraint = get_constraints_by_type(com, CS.IndicatorConstraint)[1]

    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 2])
    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 3])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [0, 2, 3])

    constr_indices = constraint.indices
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[2],
        -3,
    )
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        1,
    )
    # not actually feasible but will not be tested fully here
    CS.fix!(com, com.search_space[constr_indices[1]], 1)
    # will be tested when setting the next
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        -3,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1
    # now setting it to 1 should be feasible
    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[3],
        -3,
    )

    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[2:3]
        @test sort(CS.values(com.search_space[ind])) == [-3, 1, 2, 3]
    end
    @test sort(CS.values(com.search_space[1])) == [0, 1]
    # feasible but remove -3
    @test CS.fix!(com, com.search_space[constr_indices[1]], 1)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[2:3]
        @test sort(CS.values(com.search_space[ind])) == [1, 2, 3]
    end
    CS.values(com.search_space[constr_indices[1]]) == [1]
end

@testset "indicator is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b, Bin)
    @variable(m, -5 <= x[1:5] <= 5, Int)
    @constraint(m, b => {x in CS.GeqSet()})
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.remove_above!(
        com,
        variables[constraint.indices[2]],
        3;
        check_feasibility = false,
    )
    @test CS.remove_below!(
        com,
        variables[constraint.indices[3]],
        4;
        check_feasibility = false,
    )
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b, Bin)
    @variable(m, -5 <= x[1:5] <= 5, Int)
    @constraint(m, b => {x in CS.GeqSet()})
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.remove_above!(
        com,
        variables[constraint.indices[2]],
        3;
        check_feasibility = false,
    )
    @test CS.remove_below!(
        com,
        variables[constraint.indices[3]],
        3;
        check_feasibility = false,
    )
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b, Bin)
    @variable(m, -5 <= x[1:5] <= 5, Int)
    @constraint(m, b => {x in CS.GeqSet()})
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 0; check_feasibility = false)
    @test CS.remove_above!(
        com,
        variables[constraint.indices[2]],
        3;
        check_feasibility = false,
    )
    @test CS.remove_below!(
        com,
        variables[constraint.indices[3]],
        4;
        check_feasibility = false,
    )
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end


text:
module TestBipartiteGraphs
using Test

using Catlab.Graphs.BasicGraphs, Catlab.Graphs.BipartiteGraphs

# Undirected bipartite graphs
#############################

g = UndirectedBipartiteGraph()
add_vertices₁!(g, 2)
add_vertices₂!(g, 3)
@test (nv₁(g), nv₂(g)) == (2,3)
@test (vertices₁(g), vertices₂(g)) == (1:2, 1:3)
@test nv(g) == (2,3)
@test vertices(g) == (1:2, 1:3)
@test g == UndirectedBipartiteGraph(2, 3)

add_edge!(g, 1, 1)
add_edges!(g, [2,2], [2,3])
@test ne(g) == 3
@test edges(g) == 1:3
@test (src(g), tgt(g)) == ([1,2,2], [1,2,3])
@test inneighbors(g, 1) == [1]
@test outneighbors(g, 2) == [2,3]

rem_edge!(g, 1, 1)
@test ne(g) == 2
rem_vertex₁!(g, 1)
rem_vertex₂!(g, 1)
@test nv(g) == (1,2)
@test ne(g) == 2

# Bipartite graphs
##################

g = BipartiteGraph(2, 3)
@test nv(g) == (2, 3)

add_edge₁₂!(g, 1, 2)
add_edge₂₁!(g, 1, 1)
add_edges₁₂!(g, [2,2], [3,3])
add_edges₂₁!(g, [2,3], [1,1])
@test (ne₁₂(g), ne₂₁(g)) == (3,3)
@test (edges₁₂(g), edges₂₁(g)) == (1:3, 1:3)
@test ne(g) == (3,3)
@test edges(g) == (1:3, 1:3)
@test (src₁(g), tgt₂(g)) == ([1,2,2], [2,3,3])
@test (src₂(g), tgt₁(g)) == ([1,2,3], [1,1,1])

rem_edge₁₂!(g, 1, 2)
@test (src₁(g), tgt₂(g)) == ([2,2], [3,3])
rem_edge₂₁!(g, 1, 1)
@test (src₂(g), tgt₁(g)) == ([3,2], [1,1])
rem_vertex₁!(g, 1)
@test nv(g) == (1,3)
@test ne(g) == (2,0)
rem_vertex₂!(g, 3)
@test nv(g) == (1,2)
@test ne(g) == (0,0)

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Small sanity tests to ensure changing the rounding of float functions work
using Base: Test, MathConstants

## Float64 checks
# a + b returns a number exactly between prevfloat(1.) and 1., so its
# final result depends strongly on the utilized rounding direction.
a = prevfloat(0.5)
b = 0.5
c = 0x1p-54
d = prevfloat(1.)

# Default rounding direction, RoundNearest
@test a + b === 1.
@test - a - b === -1.
@test a - b === -c
@test b - a === c

# RoundToZero
setrounding(Float64,RoundToZero) do
    @test a + b === d
    @test - a - b === -d
    @test a - b === -c
    @test b - a === c
end

# Sanity check to see if we have returned to RoundNearest
@test a + b === 1.
@test - a - b === -1.
@test a - b == -c
@test b - a == c

# RoundUp
setrounding(Float64,RoundUp) do
    @test a + b === 1.
    @test - a - b === -d
    @test a - b === -c
    @test b - a === c
end

# RoundDown
setrounding(Float64,RoundDown) do
    @test a + b === d
    @test - a - b === -1.
    @test a - b === -c
    @test b - a === c
end

## Float32 checks

a32 = prevfloat(0.5f0)
b32 = 0.5f0
c32 = (1.f0 - prevfloat(1.f0))/2
d32 = prevfloat(1.0f0)

# Default rounding direction, RoundNearest
@test a32 + b32 === 1.0f0
@test - a32 - b32 === -1.0f0
@test a32 - b32 === -c32
@test b32 - a32 === c32

# RoundToZero
setrounding(Float32,RoundToZero) do
    @test a32 + b32 === d32
    @test - a32 - b32 === -d32
    @test a32 - b32 === -c32
    @test b32 - a32 === c32
end

# Sanity check to see if we have returned to RoundNearest
@test a32 + b32 === 1.0f0
@test - a32 - b32 === -1.0f0
@test a32 - b32 == -c32
@test b32 - a32 == c32

# RoundUp
setrounding(Float32,RoundUp) do
    @test a32 + b32 === 1.0f0
    @test - a32 - b32 === -d32
    @test a32 - b32 === -c32
    @test b32 - a32 === c32
end

# RoundDown
setrounding(Float32,RoundDown) do
    @test a32 + b32 === d32
    @test - a32 - b32 === -1.0f0
    @test a32 - b32 === -c32
    @test b32 - a32 === c32
end

# convert with rounding
for v = [sqrt(2),-1/3,nextfloat(1.0),prevfloat(1.0),nextfloat(-1.0),
         prevfloat(-1.0),nextfloat(0.0),prevfloat(0.0)]
    pn = Float32(v,RoundNearest)
    @test pn == convert(Float32,v)
    pz = Float32(v,RoundToZero)
    @test pz == setrounding(()->convert(Float32,v), Float64, RoundToZero)
    pd = Float32(v,RoundDown)
    @test pd == setrounding(()->convert(Float32,v), Float64, RoundDown)
    pu = Float32(v,RoundUp)
    @test pu == setrounding(()->convert(Float32,v), Float64, RoundUp)

    @test pn == pd || pn == pu
    @test v > 0 ? pz == pd : pz == pu
    @test pu - pd == eps(pz)
end

for T in [Float32,Float64]
    for v in [sqrt(big(2.0)),-big(1.0)/big(3.0),nextfloat(big(1.0)),
              prevfloat(big(1.0)),nextfloat(big(0.0)),prevfloat(big(0.0)),
              pi,ℯ,eulergamma,catalan,golden,
              typemax(Int64),typemax(UInt64),typemax(Int128),typemax(UInt128),0xa2f30f6001bb2ec6]
        pn = T(v,RoundNearest)
        @test pn == convert(T,BigFloat(v))
        pz = T(v,RoundToZero)
        @test pz == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundToZero)
        pd = T(v,RoundDown)
        @test pd == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundDown)
        pu = T(v,RoundUp)
        @test pu == setrounding(()->convert(T,BigFloat(v)), BigFloat, RoundUp)

        @test pn == pd || pn == pu
        @test v > 0 ? pz == pd : pz == pu
        @test isinf(pu) || pu - pd == eps(pz)
    end
end

#fenv
@test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundNearest)) == RoundNearest
@test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundToZero)) == RoundToZero
@test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundUp)) == RoundUp
@test Base.Rounding.from_fenv(Base.Rounding.to_fenv(RoundDown)) == RoundDown
@test_throws ArgumentError Base.Rounding.from_fenv(-99)

badness = 1//0
@test_throws DivideError round(Int64,badness,RoundNearestTiesAway)
@test_throws DivideError round(Int64,badness,RoundNearestTiesUp)

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
using DrWatson
@quickactivate "ChaosThroughBilliards"
include(srcdir("style.jl"))

psize = 2.0

bd = billiard_stadium(1.0, 1.0)

frames = 1800
dt = 0.0001
speed = 200

# billiard_video(
#     videodir("1_1particle.mp4"), bd, deepcopy(ps);
#     frames, backgroundcolor = :black, colors = [PURPLE],
#     plot_particles = true, tailwidth = 4, particle_size = psize, res = MAXRES,
#     dt, speed, tail = 20000, # this makes ultra fine temporal resolution
#     framerate = 60,
# )

# input:
colors = [PURPLE, TEAL]
f(p) = p.pos[2] # the function that obtains the data from the particle
total_span = 10.0

ps = particlebeam(1.0, 0.8, 0, 2, 0.0001, nothing)
ylim = (0, 1)
ylabel = "y"

billiard_video_timeseries(
    videodir("timeseries.mp4"), bd, ps, f;
    displayfigure = true, total_span, colors,
    frames, backgroundcolor = :black,
    plot_particles = true, tailwidth = 4, particle_size = psize, res = MAXRES,
    dt, speed, tail = 20000, # this makes ultra fine temporal resolution
    framerate = 60, ylabel
)

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
module SampleChainsAbstractMCMC

using SampleChains
using AbstractMCMC

using Reexport
@reexport using SampleChains
using ConcreteStructs
using Random
using TupleVectors

using TupleVectors: chainvec

export AbstractMCMCChain

@concrete struct  AbstractMCMCChain{T} <: AbstractChain{T}
    samples     # :: AbstractVector{T}
    info        # Per-sample metadata, type depends on sampler used
    meta        # Metadata associated with the sample as a whole
    state
    transform
end

function AbstractMCMCChain(samples::AbstractVector{T}, info, meta, state, transform) where {T}
    return AbstractMCMCChain{T}(samples, info, meta, state, transform)
end

TupleVectors.summarize(ch::AbstractMCMCChain) = summarize(samples(ch))

function gettransform(chain::AbstractMCMCChain)
    return getfield(chain, :transform)
end

function SampleChains.pushsample!(chain::AbstractMCMCChain, val, state)
    tr = getfield(chain, :transform)
    samplemap = tr.samplemap
    statemap = tr.statemap

    push!(samples(chain), samplemap(val))
    push!(info(chain), statemap(state))
end

function SampleChains.step!(chain::AbstractMCMCChain)
    st = getfield(chain, :state)
    mt = getfield(chain, :meta)
    rng = st.rng
    state = st.state[]
    model = mt.model
    sampler = mt.sampler
    val, state = AbstractMCMC.step(rng, model, sampler, state)
end

function SampleChains.initialize!(rng::Random.AbstractRNG, ::Type{AbstractMCMCChain}, model::AbstractMCMC.AbstractModel, spl::AbstractMCMC.AbstractSampler, samplemap, statemap)
    (val, state) = AbstractMCMC.step(rng, model, spl)

    sample0 = samplemap(val)
    info0 = statemap(state)

    samples = chainvec(sample0)
    info = chainvec(info0)
    meta = (model = model, sampler = spl)
    state = (rng=rng, state=Ref(state))
    transform = (samplemap=samplemap, statemap=statemap)

    AbstractMCMCChain(samples, info, meta, state, transform)
end

function SampleChains.initialize!(::Type{AbstractMCMCChain}, model::AbstractMCMC.AbstractModel, spl::AbstractMCMC.AbstractSampler, samplemap, statemap)
    rng = Random.GLOBAL_RNG
    return initialize!(rng, AbstractMCMCChain, model, spl, samplemap, statemap)
end

function SampleChains.drawsamples!(chain::AbstractMCMCChain, n::Int=1000)
    @cleanbreak for j in 1:n
        val, state = step!(chain)
        pushsample!(chain, val, state)
    end
    return chain
end

end


text:
# Displaying color swatches
# -------------------------

const max_width = 180
const max_height = 150
const max_pixel_size = 25

function Base.show(io::IO, ::MIME"image/svg+xml", c::Color)
    write(io,
        """
        <?xml version="1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="25mm" height="25mm" viewBox="0 0 1 1">
             <rect width="1" height="1"
                   fill="#$(hex(c))" stroke="none"/>
        </svg>
        """)
end

function Base.show(io::IO, ::MIME"image/svg+xml",
                   cs::AbstractVector{T}) where T <: Color
    show(io, MIME"image/svg+xml"(), Base.ReshapedArray(cs, (1, length(cs)), ()))
end

function Base.show(io::IO, ::MIME"image/svg+xml",
                           cs::AbstractMatrix{T}) where T <: Color
    m, n = size(cs)
    apsect_ratio = clamp(n / m,
                         max_pixel_size / max_height,
                         max_width / max_pixel_size)
    pixel_aspect = apsect_ratio / (n / m)
    scale_factor = min(max_width / (n * max_pixel_size * pixel_aspect),
                       max_height / (m * max_pixel_size),
                       1)
    ysize = max_pixel_size * scale_factor
    xsize = ysize * pixel_aspect

    xpad = n > 50 ? 0 : 1
    ypad = m > 28 ? 0 : 1

    write(io,
        """
        <?xml version="1.0" encoding="UTF-8"?>
        <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
         "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
        <svg xmlns="http://www.w3.org/2000/svg" version="1.1"
             width="$(n*xsize)mm" height="$(m*ysize)mm"
             shape-rendering="crispEdges">
        """)

    for i in 1:m, j in 1:n
        c = ndims(cs) == 2 ? cs[i, j] : cs[j]
        write(io,
            """
            <rect x="$((j-1)*xsize)mm" y="$((i-1)*ysize)mm"
                  width="$(xsize - xpad)mm" height="$(ysize - ypad)mm"
                  fill="#$(hex(c))" stroke="none" />
            """)
    end

    write(io, "</svg>")
end


text:
module CIFAR10

import ..Datasets.unpack

function getdata(dir::String)
    mkpath(dir)
    url = "https://www.cs.toronto.edu/~kriz/cifar-10-binary.tar.gz"
    println("Downloading $url...")
    path = download(url)
    run(unpack(path,dir,".gz",".tar"))
end

function readdata(data::Vector{UInt8})
    n = Int(length(data)/3073)
    x = Matrix{Float64}(3072, n)
    y = Vector{Int}(n)
    for i = 1:n
        k = (i-1) * 3073 + 1
        y[i] = Int(data[k])
        x[:,i] = data[k+1:k+3072] / 255
    end
    x = reshape(x, 32, 32, 3, n)
    x, y
end

function traindata(dir::String)
    files = [joinpath(dir,"cifar-10-batches-bin","data_batch_$i.bin") for i=1:5]
    all(isfile, files) || getdata(dir)
    data = UInt8[]
    for file in files
        append!(data, open(read,file))
    end
    readdata(data)
end

function testdata(dir::String)
    file = joinpath(dir,"cifar-10-batches-bin","test_batch.bin")
    isfile(file) || getdata(dir)
    readdata(open(read,file))
end

end


text:
using SimpleGraphs, SimpleGraphAlgorithms, ProgressMeter

"""
    Trees(n)
This type is used for iterating over all trees with vertex set `{1,2,...,n}`, like this:

```
for T ∈ Trees(6)
    some_action(T)
end
```

**WARNING**: There are `n^(n-2)` trees on `n` vertices. With `n=10`, that's
one hundred million.
"""
struct Trees
    n::Int
end

function Base.iterate(T::Trees, state = 0)
    N = T.n
    NN = Int128(N)^(N - 2)
    if state == NN
        return nothing
    end

    code = Int.(digits(state, base = N, pad = N - 2) .+ 1)
    G = prufer_restore(code)
    name(G, "Tree")
    return (G, state + 1)
end

Base.length(T::Trees) = Int128(T.n)^(T.n - 2)


text:
using Base.Test
# bank1d.jl
import WaveletScattering: Bank1D
# spec1d.jl
import WaveletScattering: Spec1D

W = Bank1D(Spec1D(n_filters_per_octave=4, n_octaves=8))
@test ndims(W) == 1