text:
import PolyhedralRelaxations

const PR = PolyhedralRelaxations

struct InputData
    f::Function
    f_name::String
    base_partition::Vector{Float64}
end

struct ErrorResults
    num_vars_for_errs::Vector{Int64}
    errs_for_num_vars::Vector{Float64}
end

function get_error_results(input_data::InputData,err_tols::Vector{Float64},num_vars::Vector{Int64})::ErrorResults
    num_vars_for_errs = Int64[]
    err_var_dict = Dict{Float64,Int64}
    for err_tol in err_tols
        milp_data, _ = PR.construct_milp_relaxation(input_data.f,
            input_data.base_partition, error_tolerance=err_tol)
        push!(num_vars_for_errs, PR.get_num_binary_variables(milp_data))
    end

    errs_for_num_vars = Float64[]
    for num in num_vars
        milp_data, _ = PR.construct_milp_relaxation(input_data.f, input_data.base_partition,
            num_additional_binary_variables=num)
        push!(errs_for_num_vars, PR.get_error_bound(milp_data))
    end

    return ErrorResults(num_vars_for_errs, errs_for_num_vars)
end

function generate_error_table()
    inputs = [
        InputData(sin, "sin", [0,π,2*π]),
        InputData(x -> x^3, "x^3", [-1.0,0.0,1.0]),
        InputData(x -> x * abs(x), "x|x|", [-2.0,0.0,2.0]),
        InputData(x -> 1 / (1+exp(-x)), "logistic(x)", [-5.0,0.0,5.0])
    ]
    err_tols = [NaN64,0.1,0.01]
    num_vars = [0,50,100]
    results = ErrorResults[]
    for input_data in inputs
        push!(results, get_error_results(input_data,err_tols,num_vars))
    end

    headers = [
        "name",
        "base_partition",
        "num_partitions_error_0",
        "num_partitions_error_0.1",
        "num_partitions_error_0.01",
        "strength_0_vars",
        "strength_50_vars",
        "strength_100_vars",
    ]

    open("out/error_results.csv", "w") do io
        write(io, get_string_from_list(headers))
        write(io, "\n")
        for i in 1:length(results)
            input_data = inputs[i]
            error_results = results[i]

            row = String[]
            push!(row, input_data.f_name)
            push!(row, string(input_data.base_partition))
            for j in 1:length(err_tols)
                push!(row, string(error_results.num_vars_for_errs[j]))
            end

            for j in 1:length(num_vars)
                push!(row, string(error_results.errs_for_num_vars[j]))
            end

            write(io, get_string_from_list(row))
            write(io, "\n")
        end
    end
end

function get_string_from_list(l::Vector{String})::String
    s = l[1]
    for i in 2:length(l)
        s = s * ";"
        s = s * l[i]
    end
    return s
end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences1D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 1)
else
    @init_parallel_stencil(Threads, Float64, 1)
end
using Plots, Printf, LinearAlgebra

@parallel function compute_flux!(qH, H, D, dx)
    @all(qH) = -D*@d(H)/dx
    return
end

@parallel function compute_rate!(ResH, dHdt, H, Hold, qH, dt, damp, dx)
    @all(ResH) = -(@inn(H) - @inn(Hold))/dt -@d(qH)/dx
    @all(dHdt) = @all(ResH) + damp*@all(dHdt)
    return
end

@parallel function compute_update!(H, dHdt, dtau)
    @inn(H) = @inn(H) + dtau*@all(dHdt)
    return
end

@views function diffusion_1D()
    # Physics
    lx     = 10.0       # domain size
    D      = 1.0        # diffusion coefficient
    ttot   = 0.6        # total simulation time
    dt     = 0.1        # physical time step
    # Numerics
    nx     = 256        # numerical grid resolution
    tol    = 1e-6       # tolerance
    itMax  = 1e5        # max number of iterations
    damp   = 1-41/nx    # damping (this is a tuning parameter, dependent on e.g. grid resolution)
    # Derived numerics
    dx     = lx/nx      # grid size
    dtau   = (1.0/(dx^2/D/2.1) + 1.0/dt)^-1 # iterative timestep
    xc     = LinRange(dx/2, lx-dx/2, nx)
    # Array allocation
    qH     = @zeros(nx-1)
    dHdt   = @zeros(nx-2)
    ResH   = @zeros(nx-2)
    # Initial condition
    H0     = Data.Array( exp.(-(xc.-lx/2).^2) )
    Hold   = @ones(nx).*H0
    H      = @ones(nx).*H0
    t = 0.0; it = 0; ittot = 0
    # Physical time loop
    while t<ttot
        iter = 0; err = 2*tol
        # Pseudo-transient iteration
        while err>tol && iter<itMax
            @parallel compute_flux!(qH, H, D, dx)
            @parallel compute_rate!(ResH, dHdt, H, Hold, qH, dt, damp, dx)
            @parallel compute_update!(H, dHdt, dtau)
            iter += 1; err = norm(ResH)/length(ResH)
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
    end
    # Analytic solution
    Hana = 1/sqrt(4*(ttot+1/4)) * exp.(-(xc.-lx/2).^2 /(4*(ttot+1/4)))
    @printf("Total time = %1.2f, time steps = %d, iterations tot = %d, error vs analytic = %1.2e \n", round(ttot, sigdigits=2), it, ittot, norm(Array(H)-Hana))
    # Visualise
    plot(xc, Array(H0), linewidth=3); display(plot!(xc, Array(H), legend=false, framestyle=:box, linewidth=3, xlabel="lx", ylabel="H", title="damped diffusion (nt=$it, iters=$ittot)"))
    return
end

diffusion_1D()

ext: jl


text:
include("required.jl")

## Conditional arguments
# outcome = "raw" # desired outcome (required)
# outcome = "bart"
# save_figures = true # should figures be overwritten (optional)

# Make sure "outcome" is defined
if !(@isdefined outcome) || !(outcome in ["raw", "bart"])
    @error "'outcome' must be either 'raw' or 'bart'"
else
    @info "'outcome' currently set to '$(outcome)'"
end

## Load presence-absence data for all species

# Load species names
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name

# Load distributions
distributions = [
    geotiff(
        SimpleSDMPredictor, joinpath("data", "raster", "distributions_$(outcome).tif"), i
    ) for i in eachindex(spenames)
]

## Y matrix (site-by-species community data table)

# Create matrix Y
Y = Ymatrix(distributions; transform=false)

# Create matrix of observed sites only
Yobs = Ymatrix(distributions; transform=false, observed=true)

## Richness

# Get richness
richness_full = richness(Y, distributions[1])

# Plot richness
richness_plot = plot_layer(
    richness_full;
    c=:cividis,
    colorbar_title="Species richness",
    size=(650, 400),
    dpi=200,
)

## LCBD

# Get relative LCBD values
lcbd_rel = lcbd(Y, distributions[1]; transform=true, relative=true)

# Get non-relative values
lcbd_full = lcbd(Y, distributions[1]; transform=true, relative=false)

# Get total beta
beta_full = beta_total(Y)

# Plot relative values
lcbd_rel_plot = plot_layer(
    lcbd_rel;
    c=:viridis,
    colorbar_title="Relative LCBD value",
    size=(650, 400),
    dpi=200,
)

# Plot absolute values (scaled for a nicer result)
scaling_factor = lcbd_full |> maximum |> log10 |> abs |> ceil |> Int
scaling_value = 10^scaling_factor
lcbd_resc = rescale(lcbd_full, extrema(lcbd_full) .* scaling_value)
lcbd_plot = plot_layer(
    lcbd_resc;
    c=:viridis,
    colorbar_title="LCBD value (x $(format(scaling_value, commas=true)))",
    size=(650, 400),
    dpi=200,
)

## Relationship

# Functions to add total beta value
rectangle(w, h, x, y) = Shape(x .+ [0, w, w, 0], y .+ [0, 0, h, h])
function rectangle!(w, h, x, y, textstring, textsize)
    return plot!(
        rectangle(w, h, x, y);
        color=:white,
        legend=:none,
        annotations=(x + (w / 2), y + (h / 2), text(textstring, textsize, :center)),
    )
end

# Plot relationship as histogram2d
lmin, lmax = extrema(lcbd_resc)
lrange = lmax - lmin

rel2d_plot = histogram2d(
    richness_full,
    lcbd_resc;
    c=:cividis,
    bins=40,
    xlabel="Richness",
    ylabel="LCBD value (x 100,000)",
    colorbar_title="Number of sites",
    xlim=(1.0, Inf),
    size=(650, 400),
    dpi=200,
)
vline!([median(richness_full)]; label=:none, linestyle=:dash, c=:grey)
hline!([median(lcbd_resc)]; label=:none, linestyle=:dash, c=:grey)
rectangle!(
    16.0,
    0.15 * lrange,
    33.0,
    lmax - 0.2 * lrange,
    "BDtot = $(round(beta_full; digits=3))",
    10,
)

## Export figures

# save_figures = true # should figures be overwritten (optional)
if (@isdefined save_figures) && save_figures == true
    @info "Figures saved ($(outcome))"
    savefig(richness_plot, joinpath("fig", outcome, "04_$(outcome)_richness.png"))
    savefig(lcbd_plot, joinpath("fig", outcome, "04_$(outcome)_lcbd.png"))
    savefig(rel2d_plot, joinpath("fig", outcome, "04_$(outcome)_relationship.png"))
else
    @info "Figures not saved ($(outcome))"
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using JSON

struct HypersliceSegment
  fp   :: PointND
  simplex :: Simplex
  d1   :: Dim
  d2   :: Dim
  p1d1 :: Float64
  p1d2 :: Float64
  p2d1 :: Float64
  p2d2 :: Float64
end

function Base.isapprox(x::HypersliceSegment, y::HypersliceSegment)
  all([
    x.fp ≈ y.fp,
    x.simplex ≈ y.simplex,
    x.d1 == y.d1,
    x.d2 == y.d2,
    x.p1d1 ≈ y.p1d1,
    x.p1d2 ≈ y.p1d2,
    x.p2d1 ≈ y.p2d1,
    x.p2d2 ≈ y.p2d2
  ])
end

function JSON.lower(hs::HypersliceSegment)
  Dict("focusPoint" => hs.fp,
       "simplex" => hs.simplex,
       "d1" => hs.d1,
       "d2" => hs.d2,
       "p1d1" => hs.p1d1,
       "p1d2" => hs.p1d2,
       "p2d1" => hs.p2d1,
       "p2d2" => hs.p2d2
  )
end

struct SliceSet2D
  problemSpec :: ProblemSpec
  slices :: Array{HypersliceSegment}
end

function savejson(filename::String, hs::SliceSet2D)
  open(filename, "w") do io
    JSON.print(io, hs)
  end
end

function readjson(filename::String)
  j = JSON.parsefile(filename, dicttype=OrderedDict)
  ps = ProblemSpec(x[1] => Tuple(x[2]) for x in j["problemSpec"])
  slices = [HypersliceSegment(s["focusPoint"], s["simplex"], s["d1"], s["d2"], s["p1d1"], s["p1d2"], s["p2d1"], s["p2d2"]) for s in j["slices"]]
  SliceSet2D(ps, slices)
end


text:
# This is an image comparison via frequency domain
# correlation. The image correlations are differenced
# and weighted against the original signal's
# autocorrelation.
#
# author: GRAM 5/1/16

function compare(s_p, s_n, s_o, sei, i, j)

	# correlate signals
	cor_p = ifft(conj(fft(s_p)) .* fft(sei[:,j,i]))
	cor_n = ifft(conj(fft(s_n)) .* fft(sei[:,j,i]))
	cor_o = ifft(conj(fft(s_o)) .* fft(sei[:,j,i]))
	acor = ifft(conj(fft(sei[:,j,i])) .* fft(sei[:,j,i]))

	# convolution builds in neglible reals
	# instead of zeros, so elim trailing imaginaries
	cor_p = real(cor_p)
	cor_n = real(cor_n)
	cor_o = real(cor_o)
	acor = real(acor)

	# normalize amplitudes
	cor_p = cor_p ./ maximum(cor_p)
	cor_n = cor_n ./ maximum(cor_n)
	cor_o = cor_o ./ maximum(cor_o)
	acor = acor ./ maximum(acor)

	# difference and weight
	dif_p = abs(acor - cor_p)
	dif_n = abs(acor - cor_n)
	dif_o = abs(acor - cor_o)

	# sum the difference vectors
	return [sum(dif_p),sum(dif_n),sum(dif_o)]

end


text:
"""
    function most_frequent(f::Function, iterable)

Return the most frequent value of a given iterable.
"""
function most_frequent(f::Function, iterable)
    d = Dict{Int, Int}()
    for element ∈ iterable
        v = f(element)
        if haskey(d, v)
            d[v] += 1
        else
            d[v] = 1
        end
    end
    candidate = 0
    key = 0
    for (k, v) in d
        if v > candidate
            key = k
            candidate = v
        end
    end
    return key
end


text:
abstract type WeightLearner end
abstract type WeightLearnerFit end

Base.@kwdef struct IHW{M<:PValueAdjustment,WL<:WeightLearner,F,T}
    multiple_testing_method::M = BenjaminiHochberg()
    folds::F = 5
    weight_learner::WL
    α::T = 0.1
end

# interface for weight learner
function learn_weights(
    weight_learner::WeightLearner,
    Ps_train,
    Xs_train,
    Xs_test,
    α,
    multiple_testing_method,
) end

struct UnitWeightLearner <: WeightLearner end

function learn_weights(
    ::UnitWeightLearner,
    Ps_train,
    Xs_train,
    Xs_test,
    α,
    multiple_testing_method,
)
    ones(size(Xs_test, 1)), nothing
end

struct IHWResult{
    I<:IHW,
    XT,
    P,
    KF<:FoldsView, #split into folds
    WF, #modelfits
    W<:AbstractWeights,
}
    pvals::P
    Xs::XT
    weights::W
    adjusted_pvals::P
    kfolds::KF
    weighting_fits::WF #Vector{}...}
    ihw_options::I
end

function fit(ihw::IHW, Ps, Xs)
    @unpack folds, multiple_testing_method, weight_learner, α = ihw

    m = length(Ps)
    kf = kfolds(shuffleobs(1:m), folds)

    ws = PriorityWeights(ones(m))
    weighting_fits = []

    # Initialize MLJ machine
    for (train_idx, test_idx) in kf
        if isa(Xs, AbstractVector)
            Xs_train_view = view(Xs, train_idx)
            Xs_test_view = view(Xs, test_idx)
        else
            Xs_train_view = view(Xs, train_idx, :)
            Xs_test_view = view(Xs, test_idx, :)
        end
        tmp_wts, tmp_fit = learn_weights(
            weight_learner,
            view(Ps, train_idx),
            Xs_train_view,
            Xs_test_view,
            α,
            multiple_testing_method,
        )

        push!(weighting_fits, tmp_fit)
        ws[test_idx] = tmp_wts .* length(test_idx) ./ sum(tmp_wts)
    end

    adj_p = MultipleTesting.adjust(Ps, ws, multiple_testing_method)

    IHWResult(Ps, Xs, ws, adj_p, kf, weighting_fits, ihw)
end

adjust(ihwres::IHWResult) = ihwres.adjusted_pvals
weights(ihwres::IHWResult) = ihwres.weights


text:
function _get_raster(::Type{WorldClim}, ::Type{Elevation}, layer::Integer, resolution=10.0)
    res = Dict(0.5 => "30s", 2.5 => "2.5m", 5.0 => "5m", 10.0 => "10m")

    path = joinpath(SimpleSDMLayers._layers_assets_path, _rasterpath(WorldClim), _rasterpath(Elevation), res[resolution])
    isdir(path) || mkpath(path)

    output_file = joinpath(path, "wc2.1_$(res[resolution])_elev.tif")
    zip_file = joinpath(path, "wc2.1_$(res[resolution])_elev.zip")

    if !isfile(path)
        if !isfile(zip_file)
            root = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/"
            stem = "wc2.1_$(res[resolution])_elev.zip"
            Downloads.download(root * stem, zip_file)
        end
        zf = ZipFile.Reader(zip_file)
        file_to_read =
            first(filter(f -> joinpath(path, f.name) == output_file, zf.files))

        if !isfile(joinpath(path, file_to_read.name))
            write(joinpath(path, file_to_read.name), read(file_to_read))
        end
        close(zf)
    end

    return joinpath(path, file_to_read.name)
end

function _get_raster(::Type{WorldClim}, ::Type{BioClim}, layer::Integer, resolution=10.0)
    1 ≤ layer ≤ 19 || throw(ArgumentError("The layer must be between 1 and 19"))

    res = Dict(0.5 => "30s", 2.5 => "2.5m", 5.0 => "5m", 10.0 => "10m")

    path = joinpath(SimpleSDMLayers._layers_assets_path, _rasterpath(WorldClim), _rasterpath(BioClim), res[resolution])
    isdir(path) || mkpath(path)

    output_file = joinpath(path, "wc2.1_$(res[resolution])_bio_$(layer).tif")
    zip_file = joinpath(path, "bioclim_2.1_$(res[resolution]).zip")

    if !isfile(path)
        if !isfile(zip_file)
            root = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/"
            stem = "wc2.1_$(res[resolution])_bio.zip"
            Downloads.download(root * stem, zip_file)
        end
        zf = ZipFile.Reader(zip_file)
        file_to_read =
            first(filter(f -> joinpath(path, f.name) == output_file, zf.files))

        if !isfile(joinpath(path, file_to_read.name))
            write(joinpath(path, file_to_read.name), read(file_to_read))
        end
        close(zf)
    end

    return joinpath(path, file_to_read.name)
end

function _get_raster(::Type{WorldClim}, ::Type{BioClim}, mod::CMIP6, fut::SharedSocioeconomicPathway, resolution=10.0, year="2021-2040")
    res = Dict(2.5 => "2.5m", 5.0 => "5m", 10.0 => "10m")

    path = joinpath(SimpleSDMLayers._layers_assets_path, _rasterpath(WorldClim), _rasterpath(BioClim), _rasterpath(mod), _rasterpath(fut), year, res[resolution])
    isdir(path) || mkpath(path)

    zip_file = joinpath(path, "$(res[resolution])_$(_rasterpath(mod))_$(_rasterpath(fut)).zip")
    if !isfile(path)
        if !isfile(zip_file)
            root = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/fut/"
            stem = "$(res[resolution])/wc2.1_$(res[resolution])_bioc_$(_rasterpath(mod))_$(_rasterpath(fut))_$(year).zip"
            Downloads.download(root * stem, zip_file)
        end
        zf = ZipFile.Reader(zip_file)
        file_to_read = only(zf.files)
        if !isfile(joinpath(path, "stack.tif"))
            write(joinpath(path, "stack.tif"), read(file_to_read))
        end
        close(zf)
    end

    return joinpath(path, "stack.tif")
end


text:
using MeshCat, MeshCatMechanisms, RigidBodyDynamics

urdf_path = joinpath(pwd(), "models/a1_quadruped/temp/a1.urdf")
a1 = MeshCatMechanisms.parse_urdf(urdf_path, remove_fixed_tree_joints = true)
a1_visuals = MeshCatMechanisms.URDFVisuals(joinpath(pwd(), "models/a1_quadruped/temp/a1.urdf"))

state = MechanismState(a1)
state_cache = StateCache(a1)
result = DynamicsResult(a1)
result_cache = DynamicsResultCache(a1)

vis = Visualizer()
mvis = MechanismVisualizer(a1, a1_visuals, vis[:base])
render(vis)


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
struct Position
    x::Number
    y::Number
end

Base.:+(a::Position, b::Tuple{Number, Number}) = Position(a.x + b[1], a.y + b[2])

Base.:-(a::Position, b::Position) = (a.x - b.x, a.y - b.y)

mutable struct Ship1
    pos::Position
    θ::Float64 # The angle measured in normalized degrees (range 0<=θ<360), with 0º = east
end

mutable struct Ship2
    pos::Position
    waypoint::Position # position of the waypoint relative to the ship
end

function parseline(line::AbstractString)
    return (line[1], parse(Int, (line[2:end])))
end

function getinput(filename)
    parseline.(readlines(filename))
end

function executeinstruction(state::Ship1, instr)
    code, num = instr
    if code == 'N'
        state.pos = state.pos + (0, num)
    elseif code == 'E'
        state.pos = state.pos + (num, 0)
    elseif code == 'S'
        state.pos = state.pos + (0, -num)
    elseif code == 'W'
        state.pos = state.pos + (-num, 0)
    elseif code == 'L'
        state.θ += num
    elseif code == 'R'
        state.θ -= num
    elseif code == 'F'
        state.pos = state.pos + (round(Int, cosd(state.θ) * num), round(Int, sind(state.θ) * num))
    else
        throw(ErrorException("Invalid code \"" * code * "\""))
    end
end

function part1(filename="day12input.txt")
    instructions = getinput(filename)
    ship = Ship1(Position(0, 0), 0)
    for instr in instructions
        executeinstruction(ship, instr)
    end
    return abs(ship.pos.x) + abs(ship.pos.y)
end

function executeinstruction2(state::Ship2, instr)
    code, num = instr
    if code == 'N'
        state.waypoint = state.waypoint + (0, num)
    elseif code == 'E'
        state.waypoint = state.waypoint + (num, 0)
    elseif code == 'S'
        state.waypoint = state.waypoint + (0, -num)
    elseif code == 'W'
        state.waypoint = state.waypoint + (-num, 0)
    elseif code == 'L'
        x, y = state.waypoint.x, state.waypoint.y
        magnitude = sqrt(x^2 + y^2)
        θ = atand(y, x)
        θ += num
        while θ >= 360
            θ -= 360
        end
        while θ < 0
            θ += 360
        end
        x = round(Int, cosd(θ) * magnitude)
        y = round(Int, sind(θ) * magnitude)
        state.waypoint = Position(x, y)
    elseif code == 'R'
        x, y = state.waypoint.x, state.waypoint.y
        magnitude = sqrt(x^2 + y^2)
        θ = atand(y, x)
        θ -= num
        while θ >= 360
            θ -= 360
        end
        while θ < 0
            θ += 360
        end
        x = round(Int, cosd(θ) * magnitude)
        y = round(Int, sind(θ) * magnitude)
        state.waypoint = Position(x, y)
    elseif code == 'F'
        dx, dy = state.waypoint.x, state.waypoint.y
        total_move = (dx*num, dy*num)
        state.pos += total_move
    else
        throw(ErrorException("Invalid code \"" * code * "\""))
    end
end

function part2(filename="day12input.txt")
    instructions = getinput(filename)
    ship = Ship2(Position(0, 0), Position(10, 1))
    for instr in instructions
        executeinstruction2(ship, instr)
    end
    return abs(ship.pos.x) + abs(ship.pos.y)
end

println(part1())
println(part2())


text:
immutable LennardJones <: Perturbation
  ϵ::Array{Float64,1}
  σ::Array{Float64,1}
end

LennardJones(ϵ::Real, σ::Real) = LennardJones([ϵ], [σ])

ncomp(lj::LennardJones)::Int = length(lj.ϵ)

function coreradius(lj::LennardJones)::Array{Float64,2}
  0.5 * pairpotential_minimizer(lj)
end

function cutoffradius(lj::LennardJones)::Array{Float64,2}
  5 * pairpotential_minimizer(lj)
end

function pairpotential(lj::LennardJones)::Array{Function,2}
  @attach(lj, ϵ, σ)

  N::Int = ncomp(lj)
  ret = Array{Function,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    ϵᵢⱼ = √(ϵ[i]*ϵ[j])
    σᵢⱼ = (σ[i] + σ[j]) / 2

    u(r::Real)::Float64 = 4ϵᵢⱼ * ((σᵢⱼ/r)^12 - (σᵢⱼ/r)^6)

    ret[i,j] = ret[j,i] = u
  end

  return ret
end

function pairpotential_minimizer(lj::LennardJones)::Array{Float64,2}
  @attach(lj, ϵ, σ)

  N::Int = ncomp(lj)

  [ 2^(1/6) * (σ[i] + σ[j]) / 2 for i in 1:N, j in 1:N ]
end

function pairpotential_derivative(lj::LennardJones)::Array{Function,2}
  @attach(lj, ϵ, σ)

  N::Int = ncomp(lj)
  ret = Array{Function,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    ϵᵢⱼ = √(ϵ[i]*ϵ[j])
    σᵢⱼ = (σ[i] + σ[j]) / 2

    u′(r::Real)::Float64 = 4ϵᵢⱼ * (-12σᵢⱼ^12 / r^13 + 6σᵢⱼ^6 / r^7)

    ret[i,j] = ret[j,i] = u′
  end

  return ret
end


text:
module GeometryBasics

    using StaticArrays, Tables, StructArrays, IterTools, LinearAlgebra

    using Base: @propagate_inbounds

    include("fixed_arrays.jl")
    include("offsetintegers.jl")
    include("basic_types.jl")
    include("metadata.jl")
    include("viewtypes.jl")
    include("geometry_primitives.jl")
    include("rectangles.jl")
    include("triangulation.jl")
    include("meshes.jl")
    include("lines.jl")
    include("boundingboxes.jl")

    export AbstractGeometry, GeometryPrimitive
    export Mat, Point, Vec
    export LineFace, Polytope, Line, NgonFace, convert_simplex
    export LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon
    export Simplex, connect, Triangle, NSimplex, Tetrahedron
    export QuadFace, metafree, coordinates, TetrahedronFace
    export TupleView, SimplexFace, Mesh, meta
    export Triangle, TriangleP
    export AbstractFace, TriangleFace, QuadFace, GLTriangleFace
    export OffsetInteger, ZeroIndex, OneIndex, GLIndex
    export FaceView, SimpleFaceView
    export AbstractPoint, PointMeta, PointWithUV
    export decompose, coordinates, faces, normals, decompose_uv, decompose_normals
    export GLTriangleFace, GLNormalMesh3D, GLPlainTriangleMesh, GLUVMesh3D, GLUVNormalMesh3D
    export AbstractMesh, Mesh, TriangleMesh
    export GLNormalMesh2D, PlainTriangleMesh

    # all the different predefined mesh types
    # Note: meshes can contain arbitrary meta information,
    export AbstractMesh, TriangleMesh, PlainMesh, GLPlainMesh, GLPlainMesh2D, GLPlainMesh3D
    export UVMesh, GLUVMesh, GLUVMesh2D, GLUVMesh3D
    export NormalMesh, GLNormalMesh, GLNormalMesh2D, GLNormalMesh3D
    export NormalUVMesh, GLNormalUVMesh, GLNormalUVMesh2D, GLNormalUVMesh3D
    export NormalUVWMesh, GLNormalUVWMesh, GLNormalUVWMesh2D, GLNormalUVWMesh3D

    # mesh creation functions
    export triangle_mesh, triangle_mesh, uv_mesh
    export uv_mesh, normal_mesh, uv_normal_mesh

    export height, origin, radius, width, widths, xwidth, yheight
    export HyperSphere, Circle, Sphere
    export Cylinder, Cylinder2, Cylinder3, Pyramid, extremity
    export Rect, Rect2D, Rect3D, IRect, IRect2D, IRect3D, FRect, FRect2D, FRect3D
    export before, during, isinside, isoutside, meets, overlaps, intersects, finishes
    export centered, direction, area, update
    export max_dist_dim, max_euclidean, max_euclideansq, min_dist_dim, min_euclidean
    export min_euclideansq, minmax_dist_dim, minmax_euclidean, minmax_euclideansq
    export self_intersections, split_intersections

end # module


text:
using YAML
using Statistics
using Cambrian
import Cambrian.selection
import Cambrian.mutate

cfg = get_config("../cfg/ga.yaml")
selection(pop::Array{<:Individual}) = Cambrian.tournament_selection(pop, cfg.tournament_size)
mutate(i::Individual) = Cambrian.mutate(i, cfg.m_rate)

function test_ga_evo(fitness::Function, d_fitness::Int)
    cfg = get_config("../cfg/ga.yaml"; d_fitness=d_fitness, id="test")
    e = GAEvo{Cambrian.FloatIndividual}(cfg, fitness)

    step!(e)
    @test length(e.population) == cfg.population.size
    for i in e.population
        @test i.fitness == fitness(i)
    end
    best = sort(e.population)[end]
    @test e.gen == 1

    print("GA step ", string(fitness))
    @timev step!(e)

    run!(e)
    @test length(e.population) == cfg.population.size
    for i in e.population
        @test i.fitness == fitness(i)
    end
    new_best = sort(e.population)[end]
    @test !(new_best < best)
    @test e.gen == cfg.n_gen
end

@testset "Genetic Algorithm" begin
    @testset "Sum of genes" begin
        function sum_eval(i::Individual)
            [sum(i.genes)]
        end
        test_ga_evo(sum_eval, 1)
    end

    @testset "Rosenbrock" begin
        function rosenbrock(i::Individual)
            x = i.genes
            y = -(sum([(1.0 - x[i])^2 + 100.0 * (x[i+1] - x[i]^2)^2
                      for i in 1:(length(x)-1)]))
            [y]
        end
        test_ga_evo(rosenbrock, 1)
    end

    @testset "Multi-objective" begin
        function moo(i::Individual)
            [mean(i.genes), std(i.genes), minimum(i.genes)]
        end
        test_ga_evo(moo, 3)
    end
end


text:
"""
	legend(cmd0::String="", arg1=nothing; kwargs...)

Make legends that can be overlaid on maps. It reads specific legend-related information from input or file file.

Full option list at [`legend`]($(GMTdoc)legend.html)

Parameters
----------

- $(GMT.opt_B)
- $(GMT.opt_J)
- $(GMT.opt_R)
- **C** | **clearance** :: [Type => Str]

    Sets the clearance between the legend frame and the internal items [4p/4p].
    ($(GMTdoc)legend.html#c)
- **D** | **pos** | **position** :: [Type => Str]  `Arg=[g|j|J|n|x]refpoint+wwidth[/height][+jjustify][+lspacing][+odx[/dy]]`

    Defines the reference point on the map for the legend using one of four coordinate systems.
    ($(GMTdoc)legend.html#d)
- **F** | **box** :: [Type => Str | Number]   `Arg=[+cclearances][+gfill][+i[[gap/]pen]][+p[pen]][+r[radius]][+s[[dx/dy/][shade]]]`

    Without further options, draws a rectangular border around the legend using *MAP_FRAME_PEN*.
    ($(GMTdoc)legend.html#f)
- **M** :: [Type => Bool]

    Modern mode only:
    ($(GMTdoc)legend.html#m)
- **S** | **scale** :: [Type => Number]

    Scale all symbol sizes by a common scale
    ($(GMTdoc)legend.html#s)
- **T** | **leg_file** :: [Type => Str]

    Modern mode only: Write hidden legend specification file to fname.
    ($(GMTdoc)legend.html#t)
- $(GMT.opt_Jz)
- $(GMT.opt_P)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_p)
- $(GMT.opt_q)
- $(GMT.opt_t)
- $(GMT.opt_savefig)
"""
function legend(cmd0::String="", arg1=nothing; first=true, kwargs...)

    gmt_proggy = (IamModern[1]) ? "legend "  : "pslegend "
	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic(gmt_proggy, cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd, _, _, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd, = parse_common_opts(d, cmd, [:F :c :p :q :t :JZ :UVXY :params], first)

	cmd = parse_type_anchor(d, cmd, [:D :pos :position],
	                        (map=("g", arg2str, 1), outside=("J", arg2str, 1), inside=("j", arg2str, 1), norm=("n", arg2str, 1), paper=("x", arg2str, 1), width=("+w", arg2str), justify="+j", spacing="+l", offset=("+o", arg2str)), 'j')
	cmd  = parse_these_opts(cmd, d, [[:C :clearance], [:M], [:S :scale], [:T :leg_file]])

	show_kwargs[1] && legend_help()
	isa(arg1, NamedTuple) && (arg1 = text_record(mk_legend(arg1)))
	r = finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1)
	gmt("destroy")
	return r
end

#=
GMT.mk_legend(gap="-0.1i", header=(text="My Map Legend", font=(24,"Times-Roman")), hline=(pen=1, offset="0.2i"), ncolumns=2, vline=(pen=1, offset=0), symbol=(marker=:circ, size="0.15i", dx_left="0.1i", fill="p300/12", dx_right="0.3i", text="This circle is hachured"), Symbol1=(marker=:ellipse, size="0.15i", dx_left="0.1i", fill=:yellow, dx_right="0.3i", text="This ellipse is yellow"), Symbol2=(marker=:wedge, size="0.15i", dx_left="0.1i", fill=:green, dx_right="0.3i", text="This wedge is green"), hline2=(pen=1, offset="0.2i"), ncolumns2=2, map_scale=(lon=5,lat=5,length="600+u+f"), gap2="0.05i", image=(width="3i", fname="@SOEST_block4.png",justify=:CT), gap3="0.05i", label=(txt="Smith et al., @%5%J. Geophys. Res., 99@%%", justify=:R, font=(9, "Times-Roman")), gap4="0.1i", text1="Let us just try some simple text that can go on a few lines.", text2="There is no easy way to predetermine how many lines may be required")
=#

function legend_help(test::Bool=false)
	io = (!test) ? stdout : IOBuffer()
	println(io, "\n\tLegend Codes\n")
	println(io, "  header=(text=?(a text) [, font=?(text font)])  (H) --> Plots a centered text string")
	println(io, "  cmap|cpt=?(CPT file name)  (A) --> Symbol or cell color fills may be given indirectly via a z-value")
	println(io, "  colorbar=(name=?(CPT file name), offset=?, height=? [, extra=?(extra opts)])  (B) --> Plot horizontal colorbar")
	println(io, "  textcolor=?(color)  (C) --> Color with which the remaining text is to be printed via z=value")
	println(io, "  hline=(pen=?(pen)[, offset=?])  (D) --> Horizontal line with specified pen across the legend")
	println(io, "  fill=?(color))  --> (F) Fill (color of pattern) for cells")
	println(io, "  gap|vspace=?(amount))  --> (G) Specifies a vertical gap of the given length")
	println(io, "  image=(name=?(img file name), width=?, justify=?(2 char code))  --> (I) Place image justified relative to the current point")
	println(io, "  label=(text=?(a text), justify=?(L|R|C)[, font=?(text font)])  --> (L) Plots a (L)eft, (C)entered, or (R)ight-justified text")
	println(io, "  map_scale=([lon=?,] lat|y=?(lat of scale), length=? [box=?(map panel), region=?, proj=?])  --> (M) Place a map scale")
	println(io, "  ncol[xx]=?(number of columns)  --> (N) Change the number of columns in the legend")
	println(io, "  paragraph=?(Bool or a String)  --> (P) Start a new text paragraph")
	println(io, "  symbol[xx]=([dx_left=?,] marker=?(symbol name), size=?(symbol size), [fill=?(color), pen=?(outline)] [, dx_right=?, label=?(text)])  --> (S) Plots the selected symbol")
	println(io, "  text[xx]=?(the text)  (T) --> One or more of these records with paragraph-text")
	println(io, "  vline=(pen=?(pen)[, offset=?])  (V) --> Draws a vertical line between columns")
end

mk_legend(nt::NamedTuple) = mk_legend(; nt...)
function mk_legend(; kwargs...)
	leg::Vector{String} = Vector{String}(undef, length(kwargs))
	c = zeros(Bool, length(kwargs))
	ky = keys(kwargs)

	function check_unused(d::Dict, opt::String, symb=nothing)
		(symb !== nothing) && del_from_dict(d, [symb])
		(length(d) > 0) && println("\tThe following options were not consumed in the legend's '$(opt)' option => ", keys(d))
	end

	for n = 1:length(ky)
		code = isa(kwargs[n], NamedTuple) ? kwargs[n] : NamedTuple([ky[n]] .=> [kwargs[n]])
		kw_str = lowercase(string(ky[n]))
		k = keys(code)
		if (kw_str == "H" || startswith(kw_str, "header"))		# code = (header=txt, font=?)	H 24p,Times-Roman My Map Legend
			d = nt2dict(code)
			f::String = ((val = find_in_dict(d, [:font])[1]) === nothing) ? "-" : font(val)
			leg[n] = "H " * f * " " * string(d[Symbol.(keys(d))[1]])
			check_unused(d, kw_str, Symbol.(keys(d))[1])

		elseif (kw_str == "A" || startswith(kw_str, "cpt") || startswith(kw_str, "cmap"))
			leg[n] = "A " * string(code[1])

		elseif (kw_str == "B" || startswith(kw_str, "colorbar"))	# code = (name="tt.cpt", offset=0.5, height=0.5)
			d = nt2dict(code)
			((val = find_in_dict(d, [:offset])[1]) === nothing) && error("Must specify the 'offset' in 'colorbar'")
			f = string(val)
			((val = find_in_dict(d, [:height])[1]) === nothing) && error("Must specify the 'height' in 'colorbar'")
			f *= " " * string(val)
			((val = find_in_dict(d, [:extra :options])[1]) !== nothing) && (f *= " " * string(val))
			leg[n] = "B " * string(d[Symbol.(keys(d))[1]]) * " " * f
			check_unused(d, kw_str, Symbol.(keys(d))[1])

		elseif (kw_str == "C" || startswith(kw_str, "textcolor"))
			leg[n] = "C " * get_color(code[1])

		elseif (kw_str == "D" || startswith(kw_str, "hline"))	# code = (hline=pen, offset=?)	D 0.2i 1p
			d = nt2dict(code)
			f = ((val = find_in_dict(d, [:offset])[1]) === nothing) ? "" : string(val)
			leg[n] = "D " * f * " " * add_opt_pen(d, [:pen], "", true)
			check_unused(d, kw_str)

		elseif (kw_str == "F" || startswith(kw_str, "fill"))
			leg[n] = "F " *  join([@sprintf("%s ", get_color(x)) for x in code])

		elseif (kw_str == "G" || startswith(kw_str, "vspace") || startswith(kw_str, "gap"))	# code = (vspace=val)	G -0.1i
			leg[n] = "G " * string(code[1])

		elseif (kw_str == "I" || startswith(kw_str, "image"))	# code = (image=fname)	I @SOEST_block4.png 3i CT
			d = nt2dict(code)
			((val = find_in_dict(d, [:width])[1]) === nothing) && error("Must specify the 'width' in 'image'")
			f = string(val)
			((val = find_in_dict(d, [:justify :justification])[1]) === nothing) && error("Must specify the 'justify' in 'image'")
			f *= " " *string(val)
			leg[n] = "I " * string(d[Symbol.(keys(d))[1]]) * " " * f
			check_unused(d, kw_str, Symbol.(keys(d))[1])

		elseif (kw_str == "L" || startswith(kw_str, "label"))	# code = L 9p,Times-Roman R Smith et al., @%5%J. Geophys. Res., 99@%%, 2000
			d = nt2dict(code)
			((val = find_in_dict(d, [:justify :justification])[1]) === nothing) && error("Must specify the 'justify' in 'label'")
			f = string(val)
			f = (((val = find_in_dict(d, [:font])[1]) === nothing) ? "-" : font(val)) * " " * uppercase(f[1])
			leg[n] = "L " * f * " " * string(d[Symbol.(keys(d))[1]])
			check_unused(d, kw_str, Symbol.(keys(d))[1])

		elseif (kw_str == "M" || occursin("scale", kw_str))		# code = (map_scale=)	M 5 5 600+u+f
			d = nt2dict(code)
			f = ((val = find_in_dict(d, [:lon :x])[1]) === nothing) ? "-" : string(val)
			((val = find_in_dict(d, [:lat :y])[1]) === nothing) && error("Must specify the 'lat or y' in map_scale")
			leg[n] = "M " * f * " " * string(val)
			((val = find_in_dict(d, [:length])[1]) === nothing) && error("Must specify the 'length' in map_scale")
			leg[n] *= " " * string(val)
			opt_R = parse_R(d, "", false, false)[1]
			opt_J = parse_J(d, "", " ", true, false, false)[1]
			opt_F = parse_F(d, "")
			leg[n] *= opt_F * opt_R * opt_J
			check_unused(d, kw_str)

		elseif (kw_str == "N" || startswith(kw_str, "ncol"))		# code = (ncolumns=?)	N 2
			leg[n] = "N " * string(code[1])

		elseif (kw_str == "P" || startswith(kw_str, "parag"))
			# NOTE: Pargraph also accepts options like pstext -M but that's too complicated to parse
			t = (isa(code[1], Bool) && code[1]) ? "" : string(code[1])
			leg[n] = "P " * t

		elseif (kw_str == "S" || startswith(kw_str, "symb"))
			# code = (symbol=:circ, size=0.15, dx_left=0.1, fill="p300/12", dx_righ=0.3, text="This circ")
			# S [dx1 symbol size fill pen [ dx2 text ]]
			d = nt2dict(code)
			marca::String = get_marker_name(d, nothing, [:symbol, :marker], false, true)[1]
			f = ((val = find_in_dict(d, [:dx_left])[1]) === nothing) ? "- " : string(val, " ");	dx1 = f
			f = ((val = find_in_dict(d, [:fill])[1]) === nothing) ? "- " : get_color(val);			fill = f
			f = add_opt_pen(d, [:pen], "", true);	pen = (f == "") ? " -" : " " * f		# TRUE to also seek (lw,lc,ls)
			((val = find_in_dict(d, [:size])[1]) === nothing) && error("Must specify the 'size' in 'symbol'")
			_size = arg2str(val)
			f = ((val = find_in_dict(d, [:label :text])[1]) === nothing) ? "" : string(val);	label = f
			if (label != "")
				f = ((val = find_in_dict(d, [:dx_right])[1]) === nothing) ? "- " : string(val, " ");	dx2 = f
				label = " " * dx2 * label
			end
			leg[n] = "S " * dx1 * marca * " " * _size * " " * fill * pen * label
			check_unused(d, kw_str)

		elseif (kw_str == "T" || startswith(kw_str, "text"))
			leg[n] = "T " * string(code[1])

		elseif (kw_str == "V" || startswith(kw_str, "vline"))		# code = (vline=pen, offset=?)	V 0 1p
			d = nt2dict(code)
			f = ((val = find_in_dict(d, [:offset])[1]) === nothing) ? "" : string(val)
			leg[n] = "V " * f * " " * add_opt_pen(d, [:pen], "", true)
			check_unused(d, kw_str)

		else
			@warn("Unrecognizable option $(ky[n]) in 'legend'.")
			c[n] = true
		end
	end
	any(c) && (leg = leg[.!c])		# Remove entries corresponding to bad options to not let go undefineds
	leg
end

# ---------------------------------------------------------------------------------------------------
legend!(cmd0::String="", arg1=nothing; kw...) = legend(cmd0, arg1; first=false, kw...)
legend(arg1; kw...)  = legend("", arg1; first=true, kw...)
legend!(arg1; kw...) = legend("", arg1; first=false, kw...)

const pslegend  = legend			# Alias
const pslegend! = legend!			# Alias


text:
using Parameters, ComponentArrays

# include("Gillespie.jl")
# include("HybridModel.jl")

function ctDNAmodel(du,u,p,t)

	# Grab dependent variables from soln array
	N,C,B,Z,I = u

	in_stochastic_region, params = p.in_stochastic_region,p.modelparams

	@unpack rN,λC,b,a,kN,kC,kB,λB,γB,θ,δZ,δB,α,β,ϕ,τ = params

    # Tumor-killing function
	γC = γB*C/(kB + C)

	T = N + C
	rC = λC + b*(T - kN)^2/(a*T^2 + (T - kN)^2)

	# RHS
	du[1] = -rN*N*log(T/kN)
	du[2] = -rC*C*log(T/kC)
	du[3] = in_stochastic_region ? 0.0 : (λB - δB)*B - γC*B
	du[4] = θ*(α*δB + β*γC)*B - δZ*Z		# Normal T cells is in there δz
	du[5] = ϕ*γC*B - τ*I

end

function compute_rates(t,u,params)

    N,C,B,Z,I = u
    @unpack rN,λC,b,a,kN,kC,kB,λB,γB,θ,δZ,δB,α,β,ϕ,τ = params

    # Get total Tcell population
    T = N+C

    # Tumor-killing function
	γC = γB*C/(kB + C)

	T = N + C
	rC = λC + b*(T - kN)^2/(a*T^2 + (T - kN)^2)

    # du[1] = -rN*N*log(T/kN)
	# du[2] = -rC*C*log(T/kC)
	# du[3] = (λB - δB)*B - γC*B
	# du[4] = θ*(α*δB + β*γC)*B - δZ*Z		# Normal T cells is in there δz
	# du[5] = ϕ*γC*B - τ*I

    return (
        λB*B,
        (δB + γC)*B
        )

end

function main()

	params=(rN = 0.16,	# normal T cell net growth rate
	kN = 500.0,			# normal T cell carrying capacity
	λC = 0.037,			# basal CAR growth rate
	kC = 119.0, #139.0,			# CAR T cell carrying capacity
	a  = 0.423,			# signal inefficiency factor
	b  = 0.525,			# immune reconstition factor
	λB  = 0.25,			# tumor birth rate (0.03 - 0.28) Wilkins 2000
	γB = 1.15,			# tumor killing rate (0.85-1.15) Kimmel
	kB = 4.05,			# Half-maximal tumor-killing rate Kimmel
	δZ = 0.1,
	δB = 0.1,
	α = 0.1,
	β = 0.1,
	ϕ = 0.1,
	θ = 10.0,
	τ = 0.1)

	sto_mat = [[0,0,1,0,0],[0,0,-1,0,0]]

	t_current = 0.0
	t_final = 200.0
	u0 = [6.0;0.36;1e6;1e2;0]

	opts = Options()

	u_current = u0

	dmparams = (modelparams=params,in_stochastic_region=false)

	# Build the deterministic model step
	deterministic_model = DeterministicModel(ctDNAmodel,u_current,dmparams,100,1.2*u0[3],(t_current,t_final),opts)

	reaction_network = ReactionNetwork(sto_mat,compute_rates,params)

	sol = simulate(deterministic_model,reaction_network)

	return sol
end

max_stars_repo_name: gregorykimmel/HybridModel

max_issues_repo_name: gregorykimmel/HybridModel

max_forks_repo_name: gregorykimmel/HybridModel


text:
#-------------------------------------------------------------------------------------------------------#
#											Input Types
#-------------------------------------------------------------------------------------------------------#
abstract type StandardArray{T,N} <: AbstractArray{T,N} end
abstract type AbstractInput{N} <: StandardArray{Float64,N} end
struct RandInput{N} <: AbstractInput{N}
	data::Array{<:Any,N}
	stages::Vector{Int}
	input_bool::Vector{Bool}
	reward_bool::Vector{Bool}
	punishment_bool::Vector{Bool}
	cum_stages::Vector{Int}
end
struct SparseInput{N} <: AbstractInput{N}
	data::Array{<:Any,N}
	stages::Vector{Int}
	input_bool::Vector{Bool}
	reward_bool::Vector{Bool}
	punishment_bool::Vector{Bool}
	cum_stages::Vector{Int}
end
struct RestInput{N} <: AbstractInput{N}
	data::Array{<:Any,N}
	stages::Vector{Int}
	input_bool::Vector{Bool}
	reward_bool::Vector{Bool}
	punishment_bool::Vector{Bool}
	cum_stages::Vector{Int}
end
struct SparseRandInput{N} <: AbstractInput{N}
	data::Array{<:Any,N}
	stages::Vector{Int}
	input_bool::Vector{Bool}
	reward_bool::Vector{Bool}
	punishment_bool::Vector{Bool}
	cum_stages::Vector{Int}
end
struct ColorInput{N} <: AbstractInput{N}
	data::Array{<:Any,N}
	stages::Vector{Int}
	input_bool::Vector{Bool}
	reward_bool::Vector{Bool}
	punishment_bool::Vector{Bool}
	cum_stages::Vector{Int}
	image::Image
end

Inputs = Union{RandInput, SparseInput, RestInput, SparseRandInput, ColorInput}

#-------------------------------------------------------------------------------------------------------#
#											Network Types
#-------------------------------------------------------------------------------------------------------#

struct NeuronLayer{T,N} <: AbstractArray{T,N}
	data::Array{T,1}
	dims::NTuple{1,Int}
end

struct NeuronLayers{T,N}
	layers::Array{NeuronLayer{T,N},1}
end

struct SynapseLayer{T,N} <: AbstractArray{T,N}
	data::Array{T,N}
	dims::NTuple{N,Int}
end

struct MBONLayer{T,N} <: AbstractArray{T,N}
	data::Array{T,N}
	dims::NTuple{N,Int}
end

struct SynapseLayers{T,N}
	layers::Array{SynapseLayer{T,N},1}
end

abstract type AbstractSparseConnection{T,Int,N} <: AbstractSparseArray{T,Int,N}
end

abstract type AbstractDenseConnection{T,N} <: AbstractArray{T,N}
end

struct SparseConnection{T,N} <: AbstractSparseConnection{T,Int,N}
	data::AbstractSparseArray{T,Int,N}
	dims::NTuple{N,Int}
end

struct DenseConnection{T,N} <: AbstractDenseConnection{T,N}
	data::Array{T,N}
	dims::NTuple{N,Int}
end

Connection = Union{SparseConnection,DenseConnection}
NormProduct = Union{Number, AbstractSparseArray}
ConnectionLayer = Union{NeuronLayer, SynapseLayer, MBONLayer}

struct DiverseLayers{T}
	layers::Array{T,1}
end

Iterable = Union{Array, Tuple}
BrainTypes = Union{NeuronLayers, SynapseLayers, DiverseLayers}

mutable struct Dopamine <: AbstractArray{Float64, 1}
	reward::Float64
	punishment::Float64
end

size(da::Dopamine) = (2,)
function getindex(da::Dopamine, i::Int)
	if i == 1
		return da.reward
	elseif i == 2
		return da.punishment
	else
		error("Tried to access a dopamine type with index != 1 or 2")
	end
end
function setindex!(da::Dopamine, f, i::Int)
	if i == 1
		da.reward = f
	elseif i == 2
		da.punishment = f
	else
		error("Tried to access a dopamine type with index != 1 or 2")
	end
end
Dopamine(num::Number) = Dopamine(num,num)


text:
module DiscreteFactor
    importall Base

    type DF
        var::Vector{ASCIIString}
        card::Vector{Int}
        value::Vector{Float64}

        function DF(var::Vector{ASCIIString}, card::Vector{Int}, value::Vector{Float64})
            if prod(card) != length(value)
                error("Cardinalities and values do not match!")
            elseif length(var) != length(card)
                error("Variables and cardinalities do not match!")
            else
                new(var, card, value)
            end
        end
    end

    include("dfactor.jl")

    export
        locateat,
        *, # infix for DFProduct
        DF,
        DFMarginDrop,
        DFMarginKeep,
        DFNormalize,
        DFPermute,
        DFProduct,
        DFReduce
end


text:
using NeXLUncertainties
using BenchmarkTools

inputs = uvs(
    Dict(
        nl"I[low,std]" => poisson(20),
        nl"I[high,std]" => poisson(15),
        nl"I[peak,std]" => poisson(1000),
        nl"t[low,std]" => uv(5.0, 0.01),
        nl"t[high,std]" => uv(8.0, 0.01),
        nl"t[peak,std]" => uv(5.6, 0.01),
        nl"i[low,std]" => uv(10.0, 0.1),
        nl"i[high,std]" => uv(10.1, 0.1),
        nl"i[peak,std]" => uv(9.9, 0.1),
        nl"R[low,std]" => uv(100.0, 0.01),
        nl"R[high,std]" => uv(110.0, 0.01),
        nl"R[peak,std]" => uv(106.0, 0.01),
        nl"I[low,unk]" => uv(30.0, sqrt(30.0)),
        nl"I[high,unk]" => uv(25.0, sqrt(25.0)),
        nl"I[peak,unk]" => uv(5000.0, sqrt(5000.0)),
        nl"t[low,unk]" => uv(3.0, 0.01),
        nl"t[high,unk]" => uv(6.0, 0.01),
        nl"t[peak,unk]" => uv(2.6, 0.01),
        nl"i[low,unk]" => uv(10.0, 0.1),
        nl"i[high,unk]" => uv(10.1, 0.1),
        nl"i[peak,unk]" => uv(9.9, 0.1),
        nl"R[low,unk]" => uv(102.0, 0.01),
        nl"R[high,unk]" => uv(108.0, 0.01),
        nl"R[peak,unk]" => uv(106.0, 0.01),
    ),
)

struct NormI <: MeasurementModel
    position::String
    sample::String
end

function NeXLUncertainties.compute(ni::NormI, inputs::LabeledValues, withJac::Bool)
    lI, lt, li = label("I[$(ni.position),$(ni.sample)]"),
    label("t[$(ni.position),$(ni.sample)]"),
    label("i[$(ni.position),$(ni.sample)]")
    I, t, i = inputs[lI], inputs[lt], inputs[li]
    labels = [label("NI[$(ni.position),$(ni.sample)]")]
    results = [I / (t * i)]
    jac = withJac ? zeros(Float64, 1, length(inputs)) : missing
    if withJac
        jac[1, indexin(inputs, lI)] = results[1] / I
        jac[1, indexin(inputs, lt)] = -results[1] / t
        jac[1, indexin(inputs, li)] = -results[1] / i
    end
    return (LabeledValues(labels, results), jac)
end

struct IChar <: MeasurementModel
    sample::String
end

function NeXLUncertainties.compute(ic::IChar, inputs::LabeledValues, withJac::Bool)
    lNIl, lNIp, lNIh =
        label.(("NI[low,$(ic.sample)]", "NI[peak,$(ic.sample)]", "NI[high,$(ic.sample)]"))
    lRl, lRp, lRh =
        label.(("R[low,$(ic.sample)]", "R[peak,$(ic.sample)]", "R[high,$(ic.sample)]"))
    NIl, NIp, NIh = inputs[lNIl], inputs[lNIp], inputs[lNIh]
    Rl, Rp, Rh = inputs[lRl], inputs[lRp], inputs[lRh]
    labels = [label("Ichar[$(ic.sample)]")]
    results = [NIp - ((Rp - Rl) * NIh + (Rh - Rp) * NIl) / (Rh - Rl)]
    jac = withJac ? zeros(Float64, 1, length(inputs)) : missing
    if withJac
        jac[1, indexin(inputs, lNIp)] = 1.0
        jac[1, indexin(inputs, lNIl)] = (Rh - Rp) / (Rh - Rl)
        jac[1, indexin(inputs, lNIh)] = (Rl - Rp) / (Rh - Rl)
        jac[1, indexin(inputs, lRp)] = (NIl - NIh) / (Rh - Rl)
        jac[1, indexin(inputs, lRl)] = (NIh - NIl) * (Rh - Rp) / ((Rh - Rl)^2)
        jac[1, indexin(inputs, lRh)] = (NIl - NIh) * (Rl - Rp) / ((Rh - Rl)^2)
    end
    return (LabeledValues(labels, results), jac)
end

struct KRatioModel <: MeasurementModel
    id::String
end

function NeXLUncertainties.compute(kr::KRatioModel, inputs::LabeledValues, withJac::Bool)
    lIstd, lIunk = nl"Ichar[std]", nl"Ichar[unk]"
    labels = [label("k[$(kr.id)]")]
    results = [inputs[lIunk] / inputs[lIstd]]
    jac = withJac ? zeros(Float64, 1, length(inputs)) : missing
    if withJac
        jac[1, indexin(inputs, lIunk)] = results[1] / inputs[lIunk]
        jac[1, indexin(inputs, lIstd)] = -results[1] / inputs[lIstd]
    end
    return (LabeledValues(labels, results), jac)
end

retain = missing # AllInputs()
norm2char(id) = MaintainInputs([label("R[$i,$id]") for i in ("low", "peak", "high")])
ICharModel(id::String) =
    IChar(id) ∘ (NormI("low", id) | NormI("peak", id) | NormI("high", id) | norm2char(id))
TotalModel(id::String) =
    (KRatioModel(id) | retain) ∘ (ICharModel("std") | ICharModel("unk") | retain)

res = TotalModel("K-L3")(inputs)
resmc = mcpropagate(TotalModel("K-L3"), inputs, 1000)

@btime mcpropagate(TotalModel("K-L3"), inputs, 1000)
@btime TotalModel("K-L3")(inputs)

lbl = nl"k[K-L3]"
println("res   = ", value(res[lbl]), " ± ", σ(res[lbl]))
println("resmc = ", value(resmc[lbl]), " ± ", σ(resmc[lbl]))


text:
using Test
using UncertainData
using Distributions
using StaticArrays
using StatsBase
using KernelDensity

#############################################
# UncertainValues module
#############################################
@testset "Uncertain values" begin
    @testset "Assign distributions" begin
        include("uncertain_values/test_assign_distributions.jl")

        @testset "Uncertain values" begin
            include("uncertain_values/test_uncertain_values.jl")
        end

        @testset "minmax" begin
            include("uncertain_values/test_minmax.jl")
        end

        @testset "CertainValue" begin
            include("uncertain_values/test_certain_values.jl")
        end

        @testset "Merging" begin
            include("uncertain_values/test_merging.jl")
        end
    end
end

@testset "Populations" begin
    include("uncertain_values/populations/test_UncertainScalarPopulation.jl")
    include("uncertain_values/populations/test_ConstrainedUncertainScalarPopulation.jl")
end

#############################################
# UncertainDatasets module
#############################################
@testset "Uncertain datasets" begin
    include("uncertain_datasets/test_uncertain_datasets.jl")
end

#############################################
# Resampling uncertain values
#############################################
@testset "Uncertain values" begin
    include("resampling/test_resampling_uncertain_values.jl")
end
#############################################
# Resampling uncertain tuples
#############################################
@testset "Uncertain tuples" begin
    include("resampling/test_resampling_uncertain_tuples.jl")
end

#############################################
# Resampling uncertain datasets
#############################################

# Resampling schemes
@testset "Resampling" begin

    @testset "BinnedResampling" begin
        include("resampling/resampling_schemes/test_BinnedResampling.jl")
    end

    @testset "BinnedWeightedResampling" begin
        include("resampling/resampling_schemes/test_BinnedWeightedResampling.jl")
    end

    @testset "BinnedMeanResampling" begin
        include("resampling/resampling_schemes/test_BinnedMeanResampling.jl")
    end

    @testset "BinnedMeanWeightedResampling" begin
        include("resampling/resampling_schemes/test_BinnedMeanWeightedResampling.jl")
    end

    @testset "ConstrainedIndexValueResampling" begin
        include("resampling/resampling_schemes/test_ConstrainedIndexValueResampling.jl")
    end

    @testset "ContrainedValueResampling" begin
        include("resampling/resampling_schemes/test_ConstrainedValueResampling.jl")
    end

    @testset "SequentialResampling" begin
        include("resampling/resampling_schemes/test_SequentialResampling.jl")
    end

    @testset "SequentialInterpolatedResampling" begin
        include("resampling/resampling_schemes/test_SequentialInterpolatedResampling.jl")
    end

    @testset "RandomSequences" begin
        include("resampling/resampling_schemes/test_RandomSequences.jl")
    end

    # Define an array of uncertain values `uvals` that we can construct datasets from.
    include("resampling/define_uncertain_values.jl")

    @testset "Vectors of uncertain values" begin

        # Resampling vectors of uncertain values
        #---------------------------------------
        include("resampling/uncertain_vectors/test_resampling_vectors.jl")
        include("resampling/uncertain_vectors/test_resampling_vectors_apply_funcs.jl")
        include("resampling/uncertain_vectors/test_resampling_vectors_constraints.jl")
    end

    # Resampling uncertain datasets
    #-------------------------------
    @testset "UncertainDataset" begin
        include("resampling/uncertain_datasets/test_resampling_datasets.jl")
    end

    # Resampling uncertain value datasets
    #-------------------------------------
    @testset "UncertainValueDataset" begin
        include("resampling/uncertain_datasets/test_resampling_datasets_uncertainvaluedataset.jl")
        include("resampling/uncertain_datasets/test_resampling_datasets_uncertainvaluedataset_apply_funcs.jl")
        include("resampling/uncertain_datasets/test_resampling_datasets_uncertainvaluedataset_constraints.jl")
    end
    # Resampling uncertain index datasets
    #-------------------------------------
    @testset "UncertainIndexDataset" begin
        include("resampling/uncertain_datasets/test_resampling_datasets_uncertainindexdataset.jl")
        include("resampling/uncertain_datasets/test_resampling_datasets_uncertainindexdataset_apply_funcs.jl")
        include("resampling/uncertain_datasets/test_resampling_datasets_uncertainindexdataset_constraints.jl")
    end

    # Resampling uncertain index-value datasets
    #-------------------------------------
    @testset "UncertainIndexValueDataset" begin
        include("resampling/uncertain_datasets/test_resampling_uncertainindexvaluedataset.jl")
        include("resampling/uncertain_datasets/test_resampling_with_schemes.jl")
    end

    # Special resampling constraints
    #-----------------------------
    @testset "Special resampling constraints" begin
        @testset "Sequential" begin
            include("resampling/uncertain_datasets/sequential/test_resampling_sequential_increasing.jl")
            include("resampling/uncertain_datasets/sequential/test_resampling_sequential_decreasing.jl")
        end
    end

    # Resampling inplace.
    #-----------------------------
    @testset "Inplace resampling" begin
        include("resampling/test_resampling_inplace.jl")
    end

    #############################################
    # Resampling uncertain datasets element-wise
    #############################################
    @testset "Element-wise" begin
        include("resampling/uncertain_datasets/test_resampling_abstractuncertainvaluedataset_elwise.jl")
    end
end

############################################
# Interpolation and binning
#############################################
@testset "Interpolation/binning" begin
    include("generic_interpolation/test_findall_nan_chunks.jl")
    include("generic_interpolation/test_interpolate_nans.jl")
    include("resampling/uncertain_datasets/test_interpolation.jl")
    include("binning/test_binning.jl")

    # Interpolation with resampling schemes
    @testset "Intp w/ resampling scheme" begin
    include("resampling/binning/test_bin_BinnedResampling.jl")
    include("resampling/binning/test_bin_BinnedWeightedResampling.jl")
    end
end

#############################
# Mathematics
#############################
@testset "Mathematics" begin
    include("mathematics/test_mathematics.jl")
end

############################################
# Uncertain statistics
#############################################
@testset "Statistics" begin
    include("uncertain_statistics/uncertain_values/test_core_stats_values_point_estimates.jl")
    include("uncertain_statistics/uncertain_values/test_core_stats_values_pair_estimates.jl")

    include("uncertain_statistics/uncertain_datasets/test_core_stats_datasets_single_dataset_estimates.jl")
    include("uncertain_statistics/uncertain_datasets/test_core_stats_datasets_pairwise_estimates.jl")

    include("uncertain_statistics/test_hypothesistests.jl")
end
#include("uncertain_statistics/test_hypothesistests_timeseries.jl")

#########################################
# Resampling with constraints and models
#########################################

# #####################
# # Resampling schemes
# #####################
# include("resampling/resampling_schemes/test_ConstrainedResampling.jl")

# #######################
# # Sampling constraints
# #######################
@testset "Resampling with constraints" begin
    include("sampling_constraints/test_sampling_constraints.jl")
    include("sampling_constraints/test_constrain_certainvalue.jl")
    include("sampling_constraints/test_constrain_population.jl")
    include("sampling_constraints/test_constrain_uncertainvalues.jl")
    include("sampling_constraints/test_constrain_uncertainvalues_kde.jl")
    include("sampling_constraints/test_constrain_uncertaindatasets.jl")
    include("sampling_constraints/test_constrain_uncertainvaluedatasets.jl")
    include("sampling_constraints/test_constrain_uncertainindexdatasets.jl")

    include("sampling_constraints/test_constrain_with_schemes.jl")

end


text:
module Test_Collision

using Modia3D
using Modia3D.StaticArrays
using Modia3D.ModiaMath
using Modia3D.Unitful

# Materials
massProperties = "Aluminium"
vmat1   = Modia3D.Material(color="LightBlue", transparency=0.5, wireframe=false)
vblue   = Modia3D.Material(color="LightBlue" , transparency=0.5)
vred    = Modia3D.Material(color="Red", transparency=0.5)
vgreen  = Modia3D.Material(color="Green", transparency=0.5)
vblack  = Modia3D.Material(color="Black", transparency=0.5)
vyellow = Modia3D.Material(color="Yellow", transparency=0.5)
vpink   = Modia3D.Material(color="Pink", transparency=0.5)
cmat    = "Steel"
vMatAABB = Modia3D.Material(color="Grey", transparency=0.95)

# dimensions
Lx  = 0.7
Dx = 0.5
Ly    = 0.2
Lz   = 1.0

groundWidth  = 1.0
groundHeight = 0.01
cyl  = Modia3D.Cylinder(1.2*Ly, Ly/2)
Tcyl = ModiaMath.rot1(90u"°")
TGround = ModiaMath.rot1(-90u"°")

@assembly CollisionSolids(world) begin
world_f1 = Modia3D.Object3D(world , r=[0.5*Lx, 0.0, groundWidth/2])

blau1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vblue;contactMaterial = cmat),    [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
blau2     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vblue; contactMaterial = cmat),   [ [0.0, 0.0, -Lz/2], [ 0.0, 0.0, Lz/2] ] )

schwarz1  = Modia3D.Part(Modia3D.Solid(Modia3D.SolidCone(Dx,Lz; relativeTipDiameter=0.0), nothing, vblack; contactMaterial = cmat),  [ [0.0, 0.0, 0.0]] )
rosa1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vpink; contactMaterial = cmat),   [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
green1    = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vgreen; contactMaterial = cmat),     [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
green2    = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), nothing, vgreen; contactMaterial = cmat),     [ [0.0, 0.0, -Lz/2]] )
gelb1     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vyellow; contactMaterial = cmat), [ [0.0, 0.0, -Lz/2],[ 0.0, 0.0, Lz/2] ] )
gelb2     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidEllipsoid(Lx,Ly,Lz), "Aluminium", vyellow; contactMaterial = cmat), [ [0.0, 0.0, -Lz/2] ] )
gelb3     = Modia3D.Part(Modia3D.Solid(Modia3D.SolidSphere(Dx), "Aluminium", vyellow; contactMaterial = cmat),    [ [0.0, 0.0, -Lz/2] ] )
rotM1      = Modia3D.Part(Modia3D.Solid(Modia3D.SolidPipe(Dx,Lz), nothing, vred; contactMaterial = cmat),          [[0.0, 0.0, -Lz/2] , [ 0.0, 0.0,Lz/2]] )
rotM2      = Modia3D.Part(Modia3D.Solid(Modia3D.SolidCylinder(Dx, Lz), nothing, vred; contactMaterial = cmat),     [[0.0, 0.0, -Lz/2], [ 0.0, 0.0,Lz/2]])

Modia3D.connect(world_f1, blau1.frames[1]; R=ModiaMath.rot2(160u"°"))
Modia3D.connect(world_f1, blau2.frames[1]; R=ModiaMath.rot2(45u"°"))

rev2 = Modia3D.Revolute(blau2.frames[2], rosa1.frames[1]; phi_start =-pi/2)

rev5 = Modia3D.Revolute(blau1.frames[2], schwarz1.frames[1]; phi_start =-pi/2)

rev3 = Modia3D.Revolute(rosa1.frames[2], green1.frames[1])

Modia3D.connect(green1.frame0, green2.frames[1]; R=ModiaMath.rot2(45u"°"))
rev4 = Modia3D.Revolute(green1.frames[2], gelb1.frames[1])

Modia3D.connect(gelb1.frames[2],gelb2.frames[1]; R=ModiaMath.rot2(130u"°"))
Modia3D.connect(gelb1.frames[2], gelb3.frames[1]; R=ModiaMath.rot2(90u"°"))

Modia3D.connect(world, rotM1.frames[1]; fixed=false, R=ModiaMath.rot3(90u"°")) #
Modia3D.connect(rotM1.frames[2], rotM2.frames[1]; R=ModiaMath.rot3(45u"°"))

end

# Bounding Boxes AABB
@assembly BoundingBoxes(world) begin
  blau1AABB     = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  blau2AABB     = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  rotM1AABB      = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  rotM2AABB      = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  schwarz1AABB  = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  rosa1AABB     = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  green1AABB    = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  green2AABB    = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  gelb1AABB     = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  gelb2AABB     = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
  gelb3AABB     = Modia3D.Object3D(world,  Modia3D.Box(0.9,0.9,0.9      , material=vMatAABB))
end

@assembly CollisionTest() begin
  world = Modia3D.Object3D(Modia3D.CoordinateSystem(1.0))

  collisionSolids = CollisionSolids(world)
  boundingBoxes   = BoundingBoxes(world)
end

collisionTest = CollisionTest(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3,nz_max = 80))

#collisionTest = CollisionTest()
Modia3D.initAnalysis!(collisionTest)

AABB = collisionTest._internal.scene.options.contactDetection.contactPairs.AABB

tStart=0.0
tEnd  =1.0

for time = range(tStart, stop=tEnd, length=101)
  s = Modia3D.linearMovement(2*Lx, tStart, tEnd, time)
  delta_phi = Modia3D.linearMovement(pi/3, tStart, tEnd, time)

  Modia3D.setAngle!(collisionTest.collisionSolids.rev2,-pi/2 + delta_phi)
  Modia3D.setAngle!(collisionTest.collisionSolids.rev3, delta_phi)
  Modia3D.setAngle!(collisionTest.collisionSolids.rev4,-delta_phi)
  Modia3D.setAngle!(collisionTest.collisionSolids.rev5, delta_phi)

  Modia3D.set_r!(collisionTest.collisionSolids.rotM1.frames[1], [0, 0, -s])

  Modia3D.updatePosition!(collisionTest)
  Modia3D.selectContactPairsWithEvent!(collisionTest)

  collisionTest.boundingBoxes.blau1AABB.data.Lx = abs(AABB[1][1].x_max - AABB[1][1].x_min)
  collisionTest.boundingBoxes.blau1AABB.data.Ly = abs(AABB[1][1].y_max - AABB[1][1].y_min)
  collisionTest.boundingBoxes.blau1AABB.data.Lz = abs(AABB[1][1].z_max - AABB[1][1].z_min)
  collisionTest.boundingBoxes.blau1AABB.r_abs = collisionTest.collisionSolids.blau1.frame0.r_abs

  collisionTest.boundingBoxes.blau2AABB.data.Lx = abs(AABB[1][2].x_max - AABB[1][2].x_min)
  collisionTest.boundingBoxes.blau2AABB.data.Ly = abs(AABB[1][2].y_max - AABB[1][2].y_min)
  collisionTest.boundingBoxes.blau2AABB.data.Lz = abs(AABB[1][2].z_max - AABB[1][2].z_min)
  collisionTest.boundingBoxes.blau2AABB.r_abs = collisionTest.collisionSolids.blau2.frame0.r_abs

  collisionTest.boundingBoxes.rotM1AABB.data.Lx = abs(AABB[2][1].x_max - AABB[2][1].x_min)
  collisionTest.boundingBoxes.rotM1AABB.data.Ly = abs(AABB[2][1].y_max - AABB[2][1].y_min)
  collisionTest.boundingBoxes.rotM1AABB.data.Lz = abs(AABB[2][1].z_max - AABB[2][1].z_min)
  collisionTest.boundingBoxes.rotM1AABB.r_abs = collisionTest.collisionSolids.rotM1.frame0.r_abs

  collisionTest.boundingBoxes.rotM2AABB.data.Lx = abs(AABB[2][2].x_max - AABB[2][2].x_min)
  collisionTest.boundingBoxes.rotM2AABB.data.Ly = abs(AABB[2][2].y_max - AABB[2][2].y_min)
  collisionTest.boundingBoxes.rotM2AABB.data.Lz = abs(AABB[2][2].z_max - AABB[2][2].z_min)
  collisionTest.boundingBoxes.rotM2AABB.r_abs = collisionTest.collisionSolids.rotM2.frame0.r_abs

  collisionTest.boundingBoxes.schwarz1AABB.data.Lx = abs(AABB[3][1].x_max - AABB[3][1].x_min)
  collisionTest.boundingBoxes.schwarz1AABB.data.Ly = abs(AABB[3][1].y_max - AABB[3][1].y_min)
  collisionTest.boundingBoxes.schwarz1AABB.data.Lz = abs(AABB[3][1].z_max - AABB[3][1].z_min)
  collisionTest.boundingBoxes.schwarz1AABB.r_abs =  [(AABB[3][1].x_max+AABB[3][1].x_min)/2,(AABB[3][1].y_max+AABB[3][1].y_min)/2,(AABB[3][1].z_max+AABB[3][1].z_min)/2 ]

  collisionTest.boundingBoxes.rosa1AABB.data.Lx = abs(AABB[4][1].x_max - AABB[4][1].x_min)
  collisionTest.boundingBoxes.rosa1AABB.data.Ly = abs(AABB[4][1].y_max - AABB[4][1].y_min)
  collisionTest.boundingBoxes.rosa1AABB.data.Lz = abs(AABB[4][1].z_max - AABB[4][1].z_min)
  collisionTest.boundingBoxes.rosa1AABB.r_abs = collisionTest.collisionSolids.rosa1.frame0.r_abs

  collisionTest.boundingBoxes.green1AABB.data.Lx = abs(AABB[5][1].x_max - AABB[5][1].x_min)
  collisionTest.boundingBoxes.green1AABB.data.Ly = abs(AABB[5][1].y_max - AABB[5][1].y_min)
  collisionTest.boundingBoxes.green1AABB.data.Lz = abs(AABB[5][1].z_max - AABB[5][1].z_min)
  collisionTest.boundingBoxes.green1AABB.r_abs = collisionTest.collisionSolids.green1.frame0.r_abs

  collisionTest.boundingBoxes.green2AABB.data.Lx = abs(AABB[5][2].x_max - AABB[5][2].x_min)
  collisionTest.boundingBoxes.green2AABB.data.Ly = abs(AABB[5][2].y_max - AABB[5][2].y_min)
  collisionTest.boundingBoxes.green2AABB.data.Lz = abs(AABB[5][2].z_max - AABB[5][2].z_min)
  collisionTest.boundingBoxes.green2AABB.r_abs = collisionTest.collisionSolids.green2.frame0.r_abs

  collisionTest.boundingBoxes.gelb1AABB.data.Lx = abs(AABB[6][1].x_max - AABB[6][1].x_min)
  collisionTest.boundingBoxes.gelb1AABB.data.Ly = abs(AABB[6][1].y_max - AABB[6][1].y_min)
  collisionTest.boundingBoxes.gelb1AABB.data.Lz = abs(AABB[6][1].z_max - AABB[6][1].z_min)
  collisionTest.boundingBoxes.gelb1AABB.r_abs = collisionTest.collisionSolids.gelb1.frame0.r_abs

  collisionTest.boundingBoxes.gelb2AABB.data.Lx = abs(AABB[6][2].x_max - AABB[6][2].x_min)
  collisionTest.boundingBoxes.gelb2AABB.data.Ly = abs(AABB[6][2].y_max - AABB[6][2].y_min)
  collisionTest.boundingBoxes.gelb2AABB.data.Lz = abs(AABB[6][2].z_max - AABB[6][2].z_min)
  collisionTest.boundingBoxes.gelb2AABB.r_abs = collisionTest.collisionSolids.gelb3.frame0.r_abs

  collisionTest.boundingBoxes.gelb3AABB.data.Lx = abs(AABB[6][3].x_max - AABB[6][3].x_min)
  collisionTest.boundingBoxes.gelb3AABB.data.Ly = abs(AABB[6][3].y_max - AABB[6][3].y_min)
  collisionTest.boundingBoxes.gelb3AABB.data.Lz = abs(AABB[6][3].z_max - AABB[6][3].z_min)
  collisionTest.boundingBoxes.gelb3AABB.r_abs = collisionTest.collisionSolids.gelb2.frame0.r_abs

  Modia3D.setComputationFlag(collisionTest)
  Modia3D.visualize!(collisionTest,time)

#=
  time += 0.01
  s = Modia3D.linearMovement(2*Lx, tStart, tEnd, time)
  Modia3D.set_r!(rotM1, [0, 0, s])
  Modia3D.updatePosition!(world)
  Modia3D.selectContactPairsWithEvent!(world)
  Modia3D.setComputationFlag(world)
  Modia3D.visualize!(world,time)
=#
  time += 0.01
end

Modia3D.closeAnalysis!(collisionTest)
println("... success of Test_Collision.jl!")
end


text:
function fdist(x::Vector{Float64},V::Float64,Q::Float64,Vs::Float64,vs::Float64,
    ts::Float64,V0::Float64,rho::Float64,beta::Float64,C::Float64,W1end::Float64,
    c0::Float64,A0::Float64,h::Float64)

    f = zeros(2);

    f[1] = x[2] .- V0./Vs .- beta.*C./Vs.*(vs.^2.*(2.*rho./beta).*((W1end./vs .- x[1])./8 .+
        c0./vs).^2 .- sqrt.(A0));
    f[2] = (x[2] .- V./Vs)./h.*ts .- ts./Vs.*(vs.^5.*(2.*rho./beta).^2.*((W1end./vs .- x[1])./8 .+
        c0./vs).^4.*0.5.*(W1end./vs .+ x[1])) .+ ts/Vs.*Q;

    return f
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
"""
    supported_constraints(con_types::Vector{Tuple{DataType, DataType}})

Returns `true` if `Function-in-Set` is supported for Dualization and throws an error if it is not.
"""
function supported_constraints(con_types::Vector{Tuple{DataType, DataType}})
    for (F, S) in con_types
        if !supported_constraint(F, S)
            error("Constraints of funtion ", F, " in the Set ", S," are not implemented")
        end
    end
    return
end

# General case
supported_constraint(::DataType, ::DataType) = false
# List of supported constraints
# SVF - Linear
supported_constraint(::Type{SVF}, ::Type{<:MOI.GreaterThan}) = true
supported_constraint(::Type{SVF}, ::Type{<:MOI.LessThan}) = true
supported_constraint(::Type{SVF}, ::Type{<:MOI.EqualTo}) = true
# SAF - Linear
supported_constraint(::Type{SAF{T}}, ::Type{MOI.GreaterThan{T}}) where T = true
supported_constraint(::Type{SAF{T}}, ::Type{MOI.LessThan{T}}) where T = true
supported_constraint(::Type{SAF{T}}, ::Type{MOI.EqualTo{T}}) where T = true
# VVF - Linear
supported_constraint(::Type{VVF}, ::Type{MOI.Nonpositives}) = true
supported_constraint(::Type{VVF}, ::Type{MOI.Nonnegatives}) = true
supported_constraint(::Type{VVF}, ::Type{MOI.Zeros}) = true
# VAF - Linear
supported_constraint(::Type{<:VAF}, ::Type{MOI.Nonpositives}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.Nonnegatives}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.Zeros}) = true
# SOC
supported_constraint(::Type{VVF}, ::Type{MOI.SecondOrderCone}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.SecondOrderCone}) = true
# RotatedSOC
supported_constraint(::Type{VVF}, ::Type{MOI.RotatedSecondOrderCone}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.RotatedSecondOrderCone}) = true
# SDP Triangle
supported_constraint(::Type{VVF}, ::Type{MOI.PositiveSemidefiniteConeTriangle}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.PositiveSemidefiniteConeTriangle}) = true
# ExponentialCone
supported_constraint(::Type{VVF}, ::Type{MOI.ExponentialCone}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.ExponentialCone}) = true
# DualExponentialCone
supported_constraint(::Type{VVF}, ::Type{MOI.DualExponentialCone}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.DualExponentialCone}) = true
# PowerCone
supported_constraint(::Type{VVF}, ::Type{<:MOI.PowerCone}) = true
supported_constraint(::Type{VAF{T}}, ::Type{MOI.PowerCone{T}}) where T = true
# DualPowerCone
supported_constraint(::Type{VVF}, ::Type{<:MOI.DualPowerCone}) = true
supported_constraint(::Type{VAF{T}}, ::Type{MOI.DualPowerCone{T}}) where T = true
# NormOneCone
supported_constraint(::Type{VVF}, ::Type{MOI.NormOneCone}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.NormOneCone}) = true
# NormInfinityCone
supported_constraint(::Type{VVF}, ::Type{MOI.NormInfinityCone}) = true
supported_constraint(::Type{<:VAF}, ::Type{MOI.NormInfinityCone}) = true

"""
    supported_objective(primal_model::MOI.ModelLike)

Returns `true` if `MOI.ObjectiveFunctionType()` is supported for Dualization and throws an error if it is not.
"""
function supported_objective(primal_model::MOI.ModelLike)
    obj_func_type = MOI.get(primal_model, MOI.ObjectiveFunctionType())
    if !supported_obj(obj_func_type)
        error("Objective functions of type $obj_func_type are not implemented")
    end
    return
end

# General case
supported_obj(::Any) = false
# List of supported objective functions
supported_obj(::Type{SVF}) = true
supported_obj(::Type{<:SAF}) = true
supported_obj(::Type{<:SQF}) = true


text:
type LogisticGrowth <: AbstractLotkaVolterra
  r::Real
  K::Real
  LogisticGrowth(r,K) = K > 0 ? new(r,K) : error("must have positive carrying capacity")
end
LogisticGrowth() = LogisticGrowth(1.,1.)

#LogisticGrowth specific data
nspecies(Z::LogisticGrowth) = 1

#values here are scalar, as opposed to vector
specificrate(Z::LogisticGrowth,x::Real) = Z.r*(1-x/Z.k)
rate(Z::LogisticGrowth,x::Real) = x*specificrate(Z)

#existence and evaluation of fixed points
somefixedpoint(Z::LogisticGrowth) = true
isfixedpoint(Z::LogisticGrowth) = true
isfixedpoint(Z::LogisticGrowth,sector::Vector{Int64}) = issector(Z,sector)
fixedpoint(Z::LogisticGrowth) = fixedpoint(Z,[1])
function fixedpoint(Z::LogisticGrowth,sector::Vector{Int64})
  if !issector(Z,sector)
    error("Not a valid sector for Logististic Growth fixed point")
  else
    return sector == [] ? 0. : Z.K
  end
end
fixedpoints(Z::LogisticGrowth) = [0.,fixedpoint(Z::LogisticGrowth)]
function fixedpoints(Z::LogisticGrowth,kind::Symbol)
  if kind == :feasible
    return [Z.K]
  elseif kind == :stable
    return isstable(Z::LogisticGrowth) ? [Z.K] : [0.]
  elseif kind == :viable
    return isstable(Z::LogisticGrowth) ? [Z.K] : []
  end
end

#feasibility, stability
somefeasible(Z::LogisticGrowth) = true
isfeasible(Z::LogisticGrowth) = true
function isfeasible(Z::LogisticGrowth,sector::Vector{Int64})
  if !issector(Z,sector)
    error("Not a valid sector for LogististiGrowth fixed point")
  else
    return sector == [1]
  end
end

somestable(Z::LogisticGrowth) = true
isstable(Z::LogisticGrowth) = isstable(Z,[1])
function isstable(Z::LogisticGrowth,sector::Vector{Int64})
  if !issector(Z,sector)
    error("Not a valid sector for LogististiGrowth fixed point")
  else
    return (sector == []) $ (Z.r > 0)
  end
end

function odeint(Z::LogisticGrowth,x0::Real,times)
  K,r = Z.K,Z.r
  ρ = x0/K
  expterm = [exp(-r*t)for t in times]
  return K*ρ./((1-ρ)*expterm + ρ)
end


text:
# # An Intro to DifferentialEquations.jl
#
# ## Basic Introduction Via Ordinary Differential Equations
#
# This notebook will get you started with DifferentialEquations.jl by introducing you to the functionality for solving ordinary differential equations (ODEs). The corresponding documentation page is the [ODE tutorial](http://docs.juliadiffeq.org/latest/tutorials/ode_example.html). While some of the syntax may be different for other types of equations, the same general principles hold in each case. Our goal is to give a gentle and thorough introduction that highlights these principles in a way that will help you generalize what you have learned.
#
# ### Background
#
# If you are new to the study of differential equations, it can be helpful to do a quick background read on [the definition of ordinary differential equations](https://en.wikipedia.org/wiki/Ordinary_differential_equation). We define an ordinary differential equation as an equation which describes the way that a variable $u$ changes, that is
#
# $$ u' = f(u,p,t) $$
#
# where $p$ are the parameters of the model, $t$ is the time variable, and $f$ is the nonlinear model of how $u$ changes. The initial value problem also includes the information about the starting value:
#
# $$ u(t_0) = u_0 $$
#
# Together, if you know the starting value and you know how the value will change with time, then you know what the value will be at any time point in the future. This is the intuitive definition of a differential equation.

# ### First Model: Exponential Growth
#
# Our first model will be the canonical exponential growth model. This model says that the rate of change is proportional to the current value, and is this:
#
# $$ u' = au $$
#
# where we have a starting value $u(0)=u_0$. Let's say we put 1 dollar into Bitcoin which is increasing at a rate of $98\%$ per year. Then calling now $t=0$ and measuring time in years, our model is:
#
# $$ u' = 0.98u $$
#
# and $u(0) = 1.0$. We encode this into Julia by noticing that, in this setup, we match the general form when
#
# $$ f(u,p,t) = 0.98u $$
#
# with $ u_0 = 1.0 $. If we want to solve this model on a time span from `t=0.0` to `t=1.0`, then we define an `ODEProblem` by specifying this function `f`, this initial condition `u0`, and this time span as follows:

using DifferentialEquations
f(u,p,t) = 0.98u
u0 = 1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
#----------------------------------------------------------------------------

# To solve our `ODEProblem` we use the command `solve`.

sol = solve(prob)
#----------------------------------------------------------------------------

# and that's it: we have succesfully solved our first ODE!
#
# #### Analyzing the Solution
#
# Of course, the solution type is not interesting in and of itself. We want to understand the solution! The documentation page which explains in detail the functions for analyzing the solution is the [Solution Handling](http://docs.juliadiffeq.org/latest/basics/solution.html) page. Here we will describe some of the basics. You can plot the solution using the plot recipe provided by [Plots.jl](http://docs.juliaplots.org/latest/):

using Plots; gr()
plot(sol)
#----------------------------------------------------------------------------

# From the picture we see that the solution is an exponential curve, which matches our intuition. As a plot recipe, we can annotate the result using any of the [Plots.jl attributes](http://docs.juliaplots.org/latest/attributes/). For example:

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") ## legend=false
#----------------------------------------------------------------------------

# Using the mutating `plot!` command we can add other pieces to our plot. For this ODE we know that the true solution is $u(t) = u_0 exp(at)$, so let's add some of the true solution to our plot:

plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")
#----------------------------------------------------------------------------

# In the previous command I demonstrated `sol.t`, which grabs the array of time points that the solution was saved at:

sol.t
#----------------------------------------------------------------------------

# We can get the array of solution values using `sol.u`:

sol.u
#----------------------------------------------------------------------------

# `sol.u[i]` is the value of the solution at time `sol.t[i]`. We can compute arrays of functions of the solution values using standard comprehensions, like:

[t+u for (u,t) in tuples(sol)]
#----------------------------------------------------------------------------

# However, one interesting feature is that, by default, the solution is a continuous function. If we check the print out again:

sol
#----------------------------------------------------------------------------

# you see that it says that the solution has a 4th order interpolation, meaning that it is a continuous function of 4th order accuracy. We can call the solution as a function of time `sol(t)`. For example, to get the value at `t=0.45`, we can use the command:

sol(0.45)
#----------------------------------------------------------------------------

# #### Controlling the Solver
#
# DifferentialEquations.jl has a common set of solver controls among its algorithms which can be found [at the Common Solver Options](http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html) page. We will detail some of the most widely used options.
#
# The most useful options are the tolerances `abstol` and `reltol`. These tell the internal adaptive time stepping engine how precise of a solution you want. Generally, `reltol` is the relative accuracy while `abstol` is the accuracy when `u` is near zero. These tolerances are local tolerances and thus are not global guarantees. However, a good rule of thumb is that the total solution accuracy is 1-2 digits less than the relative tolerances. Thus for the defaults `abstol=1e-6` and `reltol=1e-3`, you can expect a global accuracy of about 1-2 digits. If we want to get around 6 digits of accuracy, we can use the commands:

sol = solve(prob,abstol=1e-8,reltol=1e-8)
#----------------------------------------------------------------------------

# Now we can see no visible difference against the true solution:

plot(sol)
plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")
#----------------------------------------------------------------------------

# Notice that by decreasing the tolerance, the number of steps the solver had to take was `9` instead of the previous `5`. There is a trade off between accuracy and speed, and it is up to you to determine what is the right balance for your problem.
#
# Another common option is to use `saveat` to make the solver save at specific time points. For example, if we want the solution at an even grid of `t=0.1k` for integers `k`, we would use the command:

sol = solve(prob,saveat=0.1)
#----------------------------------------------------------------------------

# Notice that when `saveat` is used the continuous output variables are no longer saved and thus `sol(t)`, the interpolation, is only first order. We can save at an uneven grid of points by passing a collection of values to `saveat`. For example:

sol = solve(prob,saveat=[0.2,0.7,0.9])
#----------------------------------------------------------------------------

# By default it always saves the first and last values, but we can turn this off as well:

sol = solve(prob,saveat=[0.2,0.7,0.9],save_start = false, save_end = false)
#----------------------------------------------------------------------------

# If we need to reduce the amount of saving, we can also turn off the continuous output directly via `dense=false`:

sol = solve(prob,dense=false)
#----------------------------------------------------------------------------

# and to turn off all intermediate saving we can use `save_everystep=false`:

sol = solve(prob,save_everystep=false)
#----------------------------------------------------------------------------

# More advanced saving behaviors, such as saving functionals of the solution, are handled via the `SavingCallback` in the [Callback Library](http://docs.juliadiffeq.org/latest/features/callback_library.html#SavingCallback-1) which will be addressed later in the tutorial.

# #### Choosing Solver Algorithms
#
# There is no best algorithm for numerically solving a differential equation. When you call `solve(prob)`, DifferentialEquations.jl makes a guess at a good algorithm for your problem, given the properties that you ask for (the tolerances, the saving information, etc.). However, in many cases you may want more direct control. A later notebook will help introduce the various *algorithms* in DifferentialEquations.jl, but for now let's introduce the *syntax*.
#
# The most crucial determining factor in choosing a numerical method is the stiffness of the model. Stiffness is roughly characterized by a Jacobian `f` with large eigenvalues. That's quite mathematical, and we can think of it more intuitively: if you have big numbers in `f` (like parameters of order `1e5`), then it's probably stiff. Or, as the creator of the MATLAB ODE Suite, Lawrence Shampine, likes to define it, if the standard algorithms are slow, then it's stiff. We will go into more depth about diagnosing stiffness in a later tutorial, but for now note that if you believe your model may be stiff, you can hint this to the algorithm chooser via `alg_hints = [:stiff]`.

sol = solve(prob,alg_hints=[:stiff])
#----------------------------------------------------------------------------

# Stiff algorithms have to solve implicit equations and linear systems at each step so they should only be used when required.
#
# If we want to choose an algorithm directly, you can pass the algorithm type after the problem as `solve(prob,alg)`. For example, let's solve this problem using the `Tsit5()` algorithm, and just for show let's change the relative tolerance to `1e-6` at the same time:

sol = solve(prob,Tsit5(),reltol=1e-6)
#----------------------------------------------------------------------------

# ### Systems of ODEs: The Lorenz Equation
#
# Now let's move to a system of ODEs. The [Lorenz equation](https://en.wikipedia.org/wiki/Lorenz_system) is the famous "butterfly attractor" that spawned chaos theory. It is defined by the system of ODEs:
#
# $$ \frac{dx}{dt} = \sigma (y - x) $$
# $$ \frac{dy}{dt} = x (\rho - z) -y $$
# $$ \frac{dz}{dt} = xy - \beta z $$
#
# To define a system of differential equations in DifferentialEquations.jl, we define our `f` as a vector function with a vector initial condition. Thus, for the vector `u = [x,y,z]'`, we have the derivative function:

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end
#----------------------------------------------------------------------------

# Notice here we used the in-place format which writes the output to the preallocated vector `du`. For systems of equations the in-place format is faster. We use the initial condition $u_0 = [1.0,0.0,0.0]$ as follows:

u0 = [1.0,0.0,0.0]
#----------------------------------------------------------------------------

# Lastly, for this model we made use of the parameters `p`. We need to set this value in the `ODEProblem` as well. For our model we want to solve using the parameters $\sigma = 10$, $\rho = 28$, and $\beta = 8/3$, and thus we build the parameter collection:

p = (10,28,8/3) ## we could also make this an array, or any other type!
#----------------------------------------------------------------------------

# Now we generate the `ODEProblem` type. In this case, since we have parameters, we add the parameter values to the end of the constructor call. Let's solve this on a time span of `t=0` to `t=100`:

tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)
#----------------------------------------------------------------------------

# Now, just as before, we solve the problem:

sol = solve(prob)
#----------------------------------------------------------------------------

# The same solution handling features apply to this case. Thus `sol.t` stores the time points and `sol.u` is an array storing the solution at the corresponding time points.
#
# However, there are a few extra features which are good to know when dealing with systems of equations. First of all, `sol` also acts like an array. `sol[i]` returns the solution at the `i`th time point.

sol.t[10],sol[10]
#----------------------------------------------------------------------------

# Additionally, the solution acts like a matrix where `sol[j,i]` is the value of the `j`th variable at time `i`:

sol[2,10]
#----------------------------------------------------------------------------

# We can get a real matrix by performing a conversion:

A = convert(Array,sol)
#----------------------------------------------------------------------------

# This is the same as sol, i.e. `sol[i,j] = A[i,j]`, but now it's a true matrix. Plotting will by default show the time series for each variable:

plot(sol)
#----------------------------------------------------------------------------

# If we instead want to plot values against each other, we can use the `vars` command. Let's plot variable `1` against variable `2` against variable `3`:

plot(sol,vars=(1,2,3))
#----------------------------------------------------------------------------

# This is the classic Lorenz attractor plot, where the `x` axis is `u[1]`, the `y` axis is `u[2]`, and the `z` axis is `u[3]`. Note that the plot recipe by default uses the interpolation, but we can turn this off:

plot(sol,vars=(1,2,3),denseplot=false)
#----------------------------------------------------------------------------

# Yikes! This shows how calculating the continuous solution has saved a lot of computational effort by computing only a sparse solution and filling in the values! Note that in vars, `0=time`, and thus we can plot the time series of a single component like:

plot(sol,vars=(0,2))
#----------------------------------------------------------------------------

# ### A DSL for Parameterized Functions
#
# In many cases you may be defining a lot of functions with parameters. There exists the domain-specific language (DSL) defined by the `@ode_def` macro for helping with this common problem. For example, we can define the Lotka-Volterra equation:
#
# $$ \frac{dx}{dt} = ax - bxy $$
# $$ \frac{dy}{dt} = -cy + dxy $$
#
# as follows:

function lotka_volterra!(du,u,p,t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = -p[3]*u[2] + p[4]*u[1]*u[2]
end
#----------------------------------------------------------------------------

# However, that can be hard to follow since there's a lot of "programming" getting in the way. Instead, you can use the `@ode_def` macro:

lv! = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
#----------------------------------------------------------------------------

# We can then use the result just like an ODE function from before:

u0 = [1.0,1.0]
p = (1.5,1.0,3.0,1.0)
tspan = (0.0,10.0)
prob = ODEProblem(lv!,u0,tspan,p)
sol = solve(prob)
plot(sol)
#----------------------------------------------------------------------------

# Not only is the DSL convenient syntax, but it does some magic behind the scenes. For example, further parts of the tutorial will describe how solvers for stiff differential equations have to make use of the Jacobian in calculations. Here, the DSL uses symbolic differentiation to automatically derive that function:

lv!.Jex
#----------------------------------------------------------------------------

# The DSL can derive many other functions; this ability is used to speed up the solvers. An extension to DifferentialEquations.jl, [Latexify.jl](https://korsbo.github.io/Latexify.jl/latest/tutorials/parameterizedfunctions.html), allows you to extract these pieces as LaTeX expressions.

# ## Internal Types
#
# The last basic user-interface feature to explore is the choice of types. DifferentialEquations.jl respects your input types to determine the internal types that are used. Thus since in the previous cases, when we used `Float64` values for the initial condition, this meant that the internal values would be solved using `Float64`. We made sure that time was specified via `Float64` values, meaning that time steps would utilize 64-bit floats as well. But, by simply changing these types we can change what is used internally.
#
# As a quick example, let's say we want to solve an ODE defined by a matrix. To do this, we can simply use a matrix as input.

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
#----------------------------------------------------------------------------

# There is no real difference from what we did before, but now in this case `u0` is a `4x2` matrix. Because of that, the solution at each time point is matrix:

sol[3]
#----------------------------------------------------------------------------

# In DifferentialEquations.jl, you can use any type that defines `+`, `-`, `*`, `/`, and has an appropriate `norm`. For example, if we want arbitrary precision floating point numbers, we can change the input to be a matrix of `BigFloat`:

big_u0 = big.(u0)
#----------------------------------------------------------------------------

# and we can solve the `ODEProblem` with arbitrary precision numbers by using that initial condition:

prob = ODEProblem(f,big_u0,tspan)
sol = solve(prob)
#----------------------------------------------------------------------------

sol[1,3]
#----------------------------------------------------------------------------

# To really make use of this, we would want to change `abstol` and `reltol` to be small! Notice that the type for "time" is different than the type for the dependent variables, and this can be used to optimize the algorithm via keeping multiple precisions. We can convert time to be arbitrary precision as well by defining our time span with `BigFloat` variables:
#

prob = ODEProblem(f,big_u0,big.(tspan))
sol = solve(prob)
#----------------------------------------------------------------------------

# Let's end by showing a more complicated use of types. For small arrays, it's usually faster to do operations on static arrays via the package [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl). The syntax is similar to that of normal arrays, but for these special arrays we utilize the `@SMatrix` macro to indicate we want to create a static array.

using StaticArrays
A  = @SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0]
u0 = @SMatrix rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)
#----------------------------------------------------------------------------

sol[3]
#----------------------------------------------------------------------------

# ## Conclusion
#
# These are the basic controls in DifferentialEquations.jl. All equations are defined via a problem type, and the `solve` command is used with an algorithm choice (or the default) to get a solution. Every solution acts the same, like an array `sol[i]` with `sol.t[i]`, and also like a continuous function `sol(t)` with a nice plot command `plot(sol)`. The Common Solver Options can be used to control the solver for any equation type. Lastly, the types used in the numerical solving are determined by the input types, and this can be used to solve with arbitrary precision and add additional optimizations (this can be used to solve via GPUs for example!). While this was shown on ODEs, these techniques generalize to other types of equations as well.


text:
"""
    solve(
        equation::AdvectionEquation{PeriodicIntervalDomain{T},F},
        u,
        tlist,
        M,
        N;
        method = "filterfirst",
        subspacedim = 10,
        solver = QNDF(),
        reltol = 1e-4,
        abstol = 1e-6,
        C̄ = nothing,
    ) where {T,F}

Solve `equation` from `tlist[1]` to `tlist[2]` with initial conditions `u` and a
discretization of `M` and `N` points. If `method` is `"filterfirst"`, the equation is filtered then
discretized. If `method` is `"discretizefirst"`, the equation is discretized then filtered.
The parameter `subspacedim` controls the accuracy of the linear exponential timestepping.
"""
function solve(
    equation::AdvectionEquation{PeriodicIntervalDomain{T},F},
    u,
    tlist,
    M,
    N;
    method = "filterfirst",
    subspacedim = 10,
    solver = QNDF(),
    reltol = 1e-4,
    abstol = 1e-6,
    degmax = 10,
    λ = 1e-6,
    λ_ridge = 1e-8,
    C̄ = nothing,
) where {T,F}
    @unpack domain, filter = equation

    x = discretize(domain, M)
    ξ = discretize(domain, N)

    ū = apply_filter(u, filter, domain)
    if filter isa TopHatFilter
        h = filter.width
        dh(x) = derivative(h, x)
        α(x) = 1 / 3 * dh(x) * h(x)
    end

    # Get matrices
    C = advection_matrix(domain, N)
    D = diffusion_matrix(domain, N)

    if method == "filterfirst"
        F === TopHatFilter ||
            error("Method \"filterfirst\" is only implemented for TopHatFilter")
        ūₕ = ū.(ξ)
        A = spdiagm(α.(ξ))
        J = DiffEqArrayOperator(-C + A * D)
        p = (; J)
        Mdu_taylor!(du, u, p, t) = mul!(du, p.J, u)
        odefunction =
            ODEFunction(Mdu_taylor!, jac = (J, u, p, t) -> (J .= p.J), jac_prototype = p.J)
        problem = ODEProblem(odefunction, ūₕ, tlist, p)
        solution = OrdinaryDiffEq.solve(problem, solver; reltol, abstol)
    elseif method == "filterfirst_exp"
        F === TopHatFilter ||
            error("Method \"filterfirst\" is only implemented for TopHatFilter")
        # M == N || @warn "Method \"filterfirst\" will only use M"
        A = spdiagm(α.(ξ))
        J = DiffEqArrayOperator(-C + A * D)
        ūₕ = ū.(ξ)
        problem = ODEProblem(J, ūₕ, tlist)
        solution = OrdinaryDiffEq.solve(
            problem,
            LinearExponential(krylov = :simple, m = subspacedim),
        )
    elseif method == "discretizefirst"
        # W = filter_matrix(filter, domain, M, N; degmax, λ)
        # R = reconstruction_matrix(filter, domain, M, N; degmax, λ)
        W, R = get_W_R(filter, domain, M, N; degmax, λ)
        ūₕ = ū.(x)
        # ūₕ = W * u.(ξ)
        p = (; J = -W * C * R)
        Mdu!(du, u, p, t) = mul!(du, p.J, u)
        odefunction = ODEFunction(
            Mdu!,
            jac = (J, u, p, t) -> (J .= p.J),
            jac_prototype = p.J,
            # mass_matrix = W * R,
        )
        problem = ODEProblem(odefunction, ūₕ, tlist, p)
        solution = OrdinaryDiffEq.solve(problem, solver; reltol, abstol)
    elseif method == "discretizefirst-without-R"
        W = filter_matrix(filter, domain, M, N; degmax, λ = λ_ridge)
        # println(W[M ÷ 2, :])
        # A = factorize(W'W + λ * I)
        A = lu(W'W + λ * I)
        J = -W * C
        ūₕ = ū.(x)
        # ūₕ = W * u.(ξ)
        p = (; u☆ = zeros(N), utmp = zeros(N), J, A)
        function Mdu!_without_R(dū, ū, p, t)
            @unpack u☆, utmp, J, A = p
            mul!(utmp, W', ū)
            ldiv!(u☆, A, utmp)
            mul!(dū, J, u☆)
        end
        odefunction = ODEFunction(
            Mdu!_without_R,
            # jac = (J, u, p, t) -> (J .= p.J),
            # jac_prototype = p.J,
            # mass_matrix = W * R,
        )
        problem = ODEProblem(odefunction, ūₕ, tlist, p)
        solution = OrdinaryDiffEq.solve(problem, solver; reltol, abstol)
    elseif method == "discretizefirst-fit-Cbar"
        ūₕ = ū.(x)
        p = (; J = -C̄)
        Mdu!_C̄(du, u, p, t) = mul!(du, p.J, u)
        odefunction = ODEFunction(
            Mdu!_C̄,
            jac = (J, u, p, t) -> (J .= p.J),
            jac_prototype = p.J,
        )
        problem = ODEProblem(odefunction, ūₕ, tlist, p)
        solution = OrdinaryDiffEq.solve(problem, solver; reltol, abstol)
    else
        error("Unknown method")
    end

    solution
end