text:
# # [Linear Time Continuous System](@id linear_continuous)
#
# Similar to the [`linear time discrete example`](@ref linear_discrete), we will now estimate a linear time continuous system ``\partial_t u = A u``.
# We simulate the correspoding system using `OrdinaryDiffEq.jl` and generate a [`ContinuousDataDrivenProblem`](@ref DataDrivenProblem) from the simulated data.

using DataDrivenDiffEq
using ModelingToolkit
using LinearAlgebra
using OrdinaryDiffEq
#md using Plots

A = [-0.9 0.2; 0.0 -0.2]
u0 = [10.0; -10.0]
tspan = (0.0, 10.0)

f(u,p,t) = A*u

sys = ODEProblem(f, u0, tspan)
sol = solve(sys, Tsit5(), saveat = 0.05);

# We could use the `DESolution` to define our problem, but here we want to use the data for didactic purposes.
# For a [`ContinuousDataDrivenProblem`](@ref DataDrivenProblem), we need either the state trajectory and the timepoints or the state trajectory and its derivate.

X = Array(sol)
t = sol.t
prob = ContinuousDataDrivenProblem(X, t)

# And plot the problems data.

#md plot(prob)

# We can see that the derivative has been automatically added via a [`collocation`](@ref collocation) method, which defaults to a `LinearInterpolation`.
# We can do a visual check and compare our derivatives with the interpolation of the `ODESolution`.

#md DX = Array(sol(t, Val{1}))
#md scatter(t, DX', label = ["Solution" nothing], color = :red, legend = :bottomright)
#md plot!(t, prob.DX', label = ["Linear Interpolation" nothing], color = :black)

# Since we have a linear system, we can use `gDMD`, which approximates the generator of the dynamics

#md res = solve(prob, DMDSVD())
#md println(res)

# We see that the system has been recovered correctly, indicated by the small error and high AIC score of the result. We can confirm this by looking at the resulting [`Basis`](@ref)

#md system = result(res)
#md println(system)

# And also plot the prediction of the recovered dynamics

#md plot(res)

# Or a have a look at the metrics of the result

#md metrics(res)

# And check the parameters of the result

#md parameters(res)

# or the generator of the system

#md Matrix(generator(system))

# to see that the operator is slightly off, but within expectations.
# In a real example, this could have many reasons, e.g. noisy data, insufficient time samples or missing states.

# Sticking to the same procedure as earlier, we now use a linear sparse regression to solve the problem

using ModelingToolkit

@parameters t
@variables x[1:2](t)

basis = Basis(x, x, independent_variable = t, name = :LinearBasis)
#md print(basis) #hide

# Afterwards, we simply `solve` the already defined problem with our `Basis` and a `SparseOptimizer`

sparse_res = solve(prob, basis, STLSQ(1e-1))
#md println(sparse_res)

# Which holds the same equations
sparse_system = result(sparse_res)
#md println(sparse_system)

# Again, we can have a look at the result

#md plot(
#md     plot(prob), plot(sparse_res), layout = (1,2)
#md )

# Both results can be converted into an `ODESystem`

@named sys = ODESystem(
    equations(sparse_system),
    get_iv(sparse_system),
    states(sparse_system),
    parameters(sparse_system)
    );

# And simulated using `OrdinaryDiffEq.jl` using the (known) initial conditions and the parameter mapping of the estimation.

x0 = [x[1] => u0[1], x[2] => u0[2]]
ps = parameter_map(sparse_res)

ode_prob = ODEProblem(sys, x0, tspan, ps)
estimate = solve(ode_prob, Tsit5(), saveat = prob.t);

# And look at the result
#md plot(sol, color = :black)
#md plot!(estimate, color = :red, linestyle = :dash)

#md # ## [Copy-Pasteable Code](@id linear_continuous_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

@test all(aic(sparse_res) .>= 1e3) #src
@test all(l2error(sparse_res) .<= 5e-1) #src
@test all(determination(sparse_res) .>= 0.97) #src
@test Array(sol) ≈ Array(estimate) rtol = 5e-2 #src


text:
using Test
using MPI: mpiexec

# run tests on Travis CI in parallel
const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all")
const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3)
const TRIXI_NTHREADS   = clamp(Sys.CPU_THREADS, 2, 3)

@time @testset "Trixi.jl tests" begin
  # This is placed first since tests error out otherwise if `TRIXI_TEST == "all"`,
  # at least on some systems.
  @time if TRIXI_TEST == "all" || TRIXI_TEST == "mpi"
    # Do a dummy `@test true`:
    # If the process errors out the testset would error out as well,
    # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
    @test true

    # There are spurious test failures of Trixi.jl with MPI on Windows, see
    # https://github.com/trixi-framework/Trixi.jl/issues/901
    # To reduce their impact, we do not test MPI with coverage on Windows.
    # This reduces the chance to hit a spurious test failure by one half.
    cmd = string(Base.julia_cmd())
    coverage = occursin("--code-coverage", cmd) && !occursin("--code-coverage=none", cmd)
    if !(coverage && Sys.iswindows())
      mpiexec() do cmd
        run(`$cmd -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes $(abspath("test_mpi.jl"))`)
      end
    end
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "threaded"
    # Do a dummy `@test true`:
    # If the process errors out the testset would error out as well,
    # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
    @test true

    run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none $(abspath("test_threaded.jl"))`)
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part1"
    include("test_tree_1d.jl")
    include("test_tree_2d_part1.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part2"
    include("test_tree_2d_part2.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part3"
    include("test_tree_2d_part3.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part4"
    include("test_tree_3d_part1.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part5"
    include("test_tree_3d_part2.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "tree_part6"
    include("test_tree_3d_part3.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "structured"
    include("test_structured_1d.jl")
    include("test_structured_2d.jl")
    include("test_structured_3d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "p4est_part1"
    include("test_p4est_2d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "p4est_part2"
    include("test_p4est_3d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "unstructured_dgmulti"
    include("test_unstructured_2d.jl")
    include("test_dgmulti_1d.jl")
    include("test_dgmulti_2d.jl")
    include("test_dgmulti_3d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part1"
    include("test_unit.jl")
    include("test_visualization.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc_part2"
    include("test_special_elixirs.jl")
    include("test_performance_specializations.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "paper_self_gravitating_gas_dynamics"
    include("test_paper_self_gravitating_gas_dynamics.jl")
  end
end

ext: jl


text:
abstract type AbstractHyperViscosity end

#statsheader(a::AbstractForcing) = ""

struct NoHyperViscosity <: AbstractHyperViscosity end

    statsheader(a::AbstractHyperViscosity) = ""

    stats(a::AbstractHyperViscosity,s::AbstractSimulation) = (),()

    msg(a::NoHyperViscosity) = "\nHyper viscosity: no hyperviscosity\n\n"

    @inline nuh(::Type{NoHyperViscosity}) = nothing
    @inline nuh(a::AbstractHyperViscosity) = nuh(typeof(a))

    @inline get_hyperviscosity_exponent(::Type{NoHyperViscosity}) = nothing
    @inline get_hyperviscosity_exponent(a::AbstractHyperViscosity) = get_hyperviscosity_exponent(typeof(a))

struct HyperViscosity{νh,M} <: AbstractHyperViscosity
end

    @inline nuh(::Type{<:HyperViscosity{n,M}}) where {n,M} = n
    @inline get_hyperviscosity_exponent(::Type{<:HyperViscosity{n,M}}) where {n,M} = M

    statsheader(a::HyperViscosity) = "hdissh,hdissv,hdiss"

    stats(a::HyperViscosity,s::AbstractSimulation) = !STAT_FIL ? (hyperviscosity_stats(s.reductionh,s.reductionv,s.u,s),(0.0,0.0)) : hyperviscosity_stats(s.reductionh,s.reductionv,s.reductionh_fil,s.reductionv_fil,s.u,s)

    msg(a::HyperViscosity{nh,M}) where {nh,M} = "\nHyper viscosity: νh = $(nh), m = $(M)\n\n"

struct SpectralBarrier{initp,endp,F} <: AbstractHyperViscosity
    func::F
    function SpectralBarrier(func,initk,endpk)
        @assert endpk >= initk
        return new{initk,endpk,typeof(func)}(func)
    end
end

SpectralBarrier(initk,endpk) = SpectralBarrier(x->log(esp(x,initk,endpk,Val{1/4}())),initk,endpk)

    statsheader(a::SpectralBarrier) = "hdissh,hdissv,hdiss"

    stats(a::SpectralBarrier,s::AbstractSimulation) = !STAT_FIL ? (spectral_barrier_stats(s.reductionh,s.reductionv,s.u,a), (0.0,0.0)) : spectral_barrier_stats(s.reductionh,s.reductionv,s.reductionh_fil,s.reductionv_fil,s.u,a)

    msg(a::SpectralBarrier{initp,endp}) where {initp,endp} = "\nSpectral barrier: init = $(initp), cutoff = $(endp)\n\n"

@inline esp(k,i,p,::Val{n}) where n = abs(k) <= i ? 1.0 : abs(k) < p ? exp(-((abs(k) - i)^2)/((p)^n - abs(k)^n)) : 0.0
#@inline esp(k,i,p) = abs(k) <= i ? 1.0 : abs(k) < p ? exp(-((abs(k) - i)^2)/((p)^2 - k^2)) : 0.0

function spectral_barrier_stats(reductionh,reductionv,u::VectorField{T},hv::SpectralBarrier{ini,cut,F}) where {T,F,ini,cut}
    isrealspace(u) && fourier!(u)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            f = hv.func
            eeh = zero(T)
            eev = zero(T)
            ii = Threads.threadid()
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    k = sqrt(k2)
                    uh = u[i,j,l]
                    fk = f(k)
                    fk = ifelse(fk == -Inf, 0.0,fk)
                    eeh += (1 + (i>1)) * fk*(proj(uh.x,uh.x) + proj(uh.y,uh.y))
                    eev += (1 + (i>1)) * fk*proj(uh.z,uh.z)
                end
            resulth[ii] += eeh
            resultv[ii] += eev
        end
    end
    dissh = -sum(resulth)
    dissv = -sum(resultv)
    return dissh,dissv,dissh+dissv
end

function spectral_barrier_stats(reductionh,reductionv,reductionh_fil,reductionv_fil,u::VectorField{T},hv::SpectralBarrier{ini,cut,F}) where {T,F,ini,cut}
    isrealspace(u) && fourier!(u)
    resulth = fill!(reductionh,zero(T))
    resultv = fill!(reductionv,zero(T))
    resulth_fil = fill!(reductionh_fil,zero(T))
    resultv_fil = fill!(reductionv_fil,zero(T))
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        l = ZRANGE[kk]
        j = YRANGE[jj]
        @inbounds begin
            f = hv.func
            eeh = zero(T)
            eev = zero(T)
            eeh_fil = zero(T)
            eev_fil = zero(T)
            ii = Threads.threadid()
            kz2 = KZ[l]*KZ[l]
                kyz2 = muladd(KY[j], KY[j], kz2)
                for i in XRANGE
                    k2 = muladd(KX[i], KX[i], kyz2)
                    k = sqrt(k2)
                    uh = u[i,j,l]
                    fk = f(k)
                    fk = ifelse(fk == -Inf, 0.0,fk)
                    leeh = (1 + (i>1)) * fk*(proj(uh.x,uh.x) + proj(uh.y,uh.y))
                    leev = (1 + (i>1)) * fk*proj(uh.z,uh.z)
                    eeh += leeh
                    eev += leev

                    G2 = Gaussfilter(SPEC_FIL_D2,k2)^2

                    eeh_fil += leeh*G2
                    eev_fil += leev*G2
                end
            resulth[ii] += eeh
            resultv[ii] += eev
            resulth_fil[ii] += eeh_fil
            resultv_fil[ii] += eev_fil
        end
    end
    dissh = -sum(resulth)
    dissv = -sum(resultv)
    dissh_fil = -sum(resulth_fil)
    dissv_fil = -sum(resultv_fil)
    return (dissh,dissv,dissh+dissv),(dissh_fill,dissv_fill,dissh_fill+dissv_fill)
end


text:
# After manual inspection of assembly code...
using Primes
x = parse(Int, last(split(readline("input.txt"))))
println((x - 2)^2)
println(sum(map(i -> !isprime(x * 100 + 100_000 + 17 * i), 0:1000)))


text:
module IntrinsicPersistence

using PersistenceBarcodes

using Distances
using LightGraphs, SimpleWeightedGraphs
using LightGraphs: DijkstraState
using NearestNeighbors
using StaticArrays
using RecipesBase

include("TriangleIterator.jl")
include("result_types.jl")
include("persistence.jl")

export persistence

end


text:
function updatejs(version = "0.8.0")
    file = "$version.zip"
    if (version != "master")
        file = "v$file"
    end

    download("https://github.com/rohitvarkey/three-js/archive/$file", file)
    run(`unzip -qq $file`)
    mv("three-js-$version", joinpath(dirname(@__FILE__), "..", "assets", "bower_components", "three-js"), remove_destination=true)
end

# Install the three-js webcomponents
mktempdir(dir -> cd(updatejs, dir))


text:
# Nonlinear curve fitting.

"""
# Nonlinear multi-variable least squares

Uses a Newton like algorithm to compute the least squares fit
of a model

## Parameters:

 * `x` an array where each colum is a variable
 * `fun` The function that should be fitted (more later on)
 * `a0` An initial guess for each fitting parameter
 * `eps` Convergence criterion
 * `maxiter` Maximum number of iterations

## Specifying the model

The model is specified in argument `fun` that should be callable according
to the following signature:

`r = fun(x, a)`

where both `x` and `a` are vectors

The way the algorithm is implemented allows for the function to be implicit. This means
that we are trying to fit the relationship `fun(x, a) = 0`.

## Initial guess

The initial guess is important and should be as close as possible to the expeted values. It should, at least, give an order of magnitude of each parameter. Zero values are not recommended since in this
case no order of magnitude is available and it is assumed. The initial step is assumed to be `a0/10`.

## Convergence criterion

The initial guess `a0` will probably not fit the data very well. The largest value
of `maxerr = maxabs(fun(x, a0))` (maximum error) is used as a reference and
every time a new approximation is computed, if the change in paramaters is compared to this
reference error. If it is small, `maxabs(da) < eps*maxerr`, the algorithm converged.

## Return values

A tuple containing:

 * The estimated coefficients
 * A `Bool` stating whether the algorithm converged
 * Number of iterations it took to converge.

"""
function nonlinear_fit(x, fun, a0, eps=1e-8, maxiter=200)

    na = length(a0)
    np = size(x, 1)
    nv = size(x, 2)
    for i in 1:na
        if a0[i] == 0
            a0[i] = 0.01
        end
    end

    da = a0 / 10
    a1 = zeros(na)
    a = zeros(na)
    for k = 1:na
        a1[k] = a0[k] - da[k]
    end
    xp = zeros(nv)
    A = zeros(np, na)
    r0 = zeros(np)
    r1 = zeros(np)
    for p = 1:np
        for k = 1:nv
            xp[k] = x[p,k]
        end
        r0[p] = fun(xp, a0)
        r1[p] = fun(xp, a1)
    end
    maxerr = maximum(abs, [maximum(abs, r0), maximum(abs, r1)])
    iter = 1
    convergence = false
    for i = 1:maxiter
        iter = i
        for p = 1:np
            for k = 1:nv
                xp[k] = x[p,k]
            end
            for k = 1:na
                aa = a1[k]
                a1[k] = a1[k] + da[k]
                r = fun(xp, a1)
                a1[k] = aa

                A[p,k] = -(r1[p] - r) / da[k]
            end
        end
        da = A \ r1
        for k = 1:na
            a1[k] = a1[k] - da[k]
        end
        for p = 1:np
            r0[p] = r1[p]
            for k = 1:nv
                xp[k] = x[p,k]
            end
            r1[p] = fun(xp, a1)
        end
        if maximum(abs, r1) < eps * maxerr
            convergence = true
            break
        end
    end

    return a1, convergence, iter
end

function nonlinear_fit0(x, fun, dflst, a0, eps=1e-8, maxiter=200)

    na = length(a0)
    np = size(x, 1)
    for i in 1:na
        if a0[i] == 0
            a0[i] = 0.01
        end
    end
    aref = eps .* abs(a0)

    A = zeros(Float64, np, na)
    r = zeros(Float64, np)
    iter = 1
    convergence = false
    for i = 1:maxiter
        iter = i
        for p = 1:np
            xp = x[p,:]
            r[p] = -fun(xp, a0)
            for k = 1:na
                A[p,k] = dflst(k, xp, a0)
            end
        end
        da = least_squares(A, r)
        a0 = a0 + da
        if maxabs(da/aref) < eps
            convergence = true
            break
        end
    end

    return a0, convergence, iter
end

function makeDerivFun(fun, na, da=1e-7)

    if length(da)==1
        da = [da for i in 1:na]
    end

    dfun(k, x, a) = begin
        a[k] = a[k] + da[k]
        x1 = fun(x, a)
        a[k] = a[k] - 2*da[k]
        x2 = fun(x, a)
        a[k] = a[k] + da[k]
        (x1 - x2) / (2*da[k])
    end

    return dfun
end


text:
struct CollaborativePredatorPreyHexWorldDecPOMDP
    hexes::Vector{Tuple{Int, Int}}
    hexWorldDiscreteMDP::DiscreteMDP
end

n_agents(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP) = 2

ordered_states(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, i::Int) = vec(collect(1:length(decpomdp.hexes)))
ordered_states(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP) = vec(collect(Iterators.product([ordered_states(decpomdp, i) for i in 1:(n_agents(decpomdp) + 1)]...)))

ordered_actions(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, i::Int) = vec(collect(1:n_actions(decpomdp.hexWorldDiscreteMDP)))
ordered_joint_actions(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP) = vec(collect(Iterators.product([ordered_actions(decpomdp, i) for i in 1:n_agents(decpomdp)]...)))

n_actions(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, i::Int) = length(ordered_actions(decpomdp, i))
n_joint_actions(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP) = length(ordered_joint_actions(decpomdp))

ordered_observations(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, i::Int) = [true, false]
ordered_joint_observations(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP) = vec(collect(Iterators.product([ordered_observations(decpomdp, i) for i in 1:n_agents(decpomdp)]...)))

n_observations(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, i::Int) = length(ordered_observations(decpomdp, i))
n_joint_observations(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP) = length(ordered_joint_observations(decpomdp))

function transition(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, s, a, s′)
    prob = 0.0

    predator1Hex = decpomdp.hexes[s[1]]
    predator2Hex = decpomdp.hexes[s[2]]
    preyHex = decpomdp.hexes[s[3]]
    preyHex′ = decpomdp.hexes[s′[3]]

    preyValidHexes = vcat([preyHex],
                          [hex for hex in hex_neighbors(preyHex)
                           if hex in decpomdp.hexes && hex != predator1Hex && hex != predator2Hex])

    # When a prey is caught (new prey born), it teleports to a random location and the predator remains (eating).
    # Otherwise, both transition following HexWorldMDP.
    if s[1] == s[3] && s[2] == s[3]
        prob = Float64(s′[1] == s[1] && s′[2] == s[2]) / length(decpomdp.hexes)
    elseif s[1] == s[3] && s[2] != s[3]
        prob = (Float64(s′[1] == s[1])
                * decpomdp.hexWorldDiscreteMDP.T[s[2], a[2], s′[2]]
                / length(decpomdp.hexes))
    elseif s[1] != s[3] && s[2] == s[3]
        prob = (decpomdp.hexWorldDiscreteMDP.T[s[1], a[1], s′[1]]
                * Float64(s′[2] == s[2])
                / length(decpomdp.hexes))
    else
        prob = (decpomdp.hexWorldDiscreteMDP.T[s[1], a[1], s′[1]]
                * decpomdp.hexWorldDiscreteMDP.T[s[2], a[2], s′[2]]
                * Float64(preyHex′ in preyValidHexes) / length(preyValidHexes))
    end
    return prob
end

function joint_observation(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, a, s′, o)
    prob = 1.0

    predator1Hex = decpomdp.hexes[s′[1]]
    predator2Hex = decpomdp.hexes[s′[2]]
    preyHex = decpomdp.hexes[s′[3]]

    if predator1Hex == preyHex || preyHex in hex_neighbors(predator1Hex)
        if o[1] == true
            prob *= 1.0
        else
            prob *= 0.0
        end
    else
        if o[1] == false
            prob *= 1.0
        else
            prob *= 0.0
        end
    end

    if predator2Hex == preyHex || preyHex in hex_neighbors(predator2Hex)
        if o[2] == true
            prob *= 1.0
        else
            prob *= 0.0
        end
    else
        if o[2] == false
            prob *= 1.0
        else
            prob *= 0.0
        end
    end

    return prob
end

function reward(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP, s, a)
    r = 0.0

    # Predators get -1 for moving and 10 for catching the prey.
    if s[1] == s[3] || s[2] == s[3]
        return 10.0
    else
        return -1.0
    end

    return r
end

function DecPOMDP(decpomdp::CollaborativePredatorPreyHexWorldDecPOMDP)
    return DecPOMDP(
        decpomdp.hexWorldDiscreteMDP.γ,
        vec(collect(1:n_agents(decpomdp))),
        ordered_states(decpomdp),
        [ordered_actions(decpomdp, i) for i in 1:n_agents(decpomdp)],
        [ordered_observations(decpomdp, i) for i in 1:n_agents(decpomdp)],
        (s, a, s′) -> transition(decpomdp, s, a, s′),
        (a, s′, o) -> joint_observation(decpomdp, a, s′, o),
        (s, a) -> reward(decpomdp, s, a)
    )
end

function CollaborativePredatorPreyHexWorldDecPOMDP(hexes::Vector{Tuple{Int,Int}},
                                r_bump_border::Float64,
                                p_intended::Float64,
                                γ::Float64)
    hexWorld = HexWorldMDP(hexes,
                           r_bump_border,
                           p_intended,
                           Dict{Tuple{Int64,Int64},Float64}(),
                           γ)
    mdp = hexWorld.mdp
    return CollaborativePredatorPreyHexWorldDecPOMDP(hexes, mdp)
end

# const SimpleCollaborativePredatorPreyHexWorld = CollaborativePredatorPreyHexWorldDecPOMDP(
#     [
#          (0, 1),
#      (0, 0), (1, 0),
#      ],
#     HexWorldRBumpBorder,
#     HexWorldPIntended,
#     0.95
# )

function SimpleCollaborativePredatorPreyHexWorld()
    SimpleCollaborativePredatorPreyHexWorld = CollaborativePredatorPreyHexWorldDecPOMDP(
        [
             (0, 1),
         (0, 0), (1, 0),
         ],
        HexWorldRBumpBorder,
        HexWorldPIntended,
        0.95
    )
    return SimpleCollaborativePredatorPreyHexWorld
end

# const CircleCollaborativePredatorPreyHexWorld = CollaborativePredatorPreyHexWorldDecPOMDP(
#     [
#          (-1, 2), (0, 2),
#      (-1, 1),         (1, 1),
#          (0, 0), (1, 0),
#      ],
#     HexWorldRBumpBorder,
#     HexWorldPIntended,
#     0.95
# )

function CircleCollaborativePredatorPreyHexWorld()
    CircleCollaborativePredatorPreyHexWorld = CollaborativePredatorPreyHexWorldDecPOMDP(
        [
             (-1, 2), (0, 2),
         (-1, 1),         (1, 1),
             (0, 0), (1, 0),
         ],
        HexWorldRBumpBorder,
        HexWorldPIntended,
        0.95
    )
    return CircleCollaborativePredatorPreyHexWorld
end

# const CollaborativePredatorPreyHexWorld = CollaborativePredatorPreyHexWorldDecPOMDP(
#     [
#      (-1, 2), (0, 2), (1, 2),
#      (-1, 1), (1, 1), (3, 1), (4, 1),
#      (0, 0), (1, 0), (2, 0), (3, 0), (4, 0),
#      ],
#     HexWorldRBumpBorder,
#     HexWorldPIntended,
#     0.95
# )

function CollaborativePredatorPreyHexWorld()
    CollaborativePredatorPreyHexWorld = CollaborativePredatorPreyHexWorldDecPOMDP(
        [
         (-1, 2), (0, 2), (1, 2),
         (-1, 1), (1, 1), (3, 1), (4, 1),
         (0, 0), (1, 0), (2, 0), (3, 0), (4, 0),
         ],
        HexWorldRBumpBorder,
        HexWorldPIntended,
        0.95
    )
    return CollaborativePredatorPreyHexWorld
end


text:
include("types/location.jl")
include("types/lintmessage.jl")

# in fact our parent and Base.parent have the same meaning: as if climbing a
# tree toward the root
import Base: parent

function Typeof(x::ANY)
    if x === Vararg
        typeof(x)
    elseif isa(x, Type)
        Type{x}
    else
        typeof(x)
    end
end

@compat abstract type AdditionalVarInfo end
extractobject(_::AdditionalVarInfo) = Nullable()

"""
A struct with information about a variable.
"""
type VarInfo
    location::Location
    typeactual::Type

    "The number of times the variable has been used."
    usages::Int

    """
    The source of the variable. Currently possible values are `:defined`,
    `:used`, and `:imported`.
    """
    source::Symbol

    "Additional known information about the particular object."
    extra::Nullable{AdditionalVarInfo}

    VarInfo(loc::Location = UNKNOWN_LOCATION, t::Type = Any;
            source::Symbol = :defined) =
        new(loc, t, 0, source, Nullable())
end

VarInfo(vi::VarInfo; source::Symbol = :defined) =
    VarInfo(location(vi), vi.typeactual)

location(vi::VarInfo) = vi.location
registeruse!(vi::VarInfo) = (vi.usages += 1; vi)
usages(vi::VarInfo) = vi.usages
source(vi::VarInfo) = vi.source
function info!(vi::VarInfo, info::AdditionalVarInfo)
    vi.extra = Nullable(info)
end

extractobject(vi::VarInfo) =
    flatten(BROADCAST(extractobject, vi.extra))

immutable ModuleInfo <: AdditionalVarInfo
    name          :: Symbol
    globals       :: Dict{Symbol, VarInfo}
    exports       :: Set{Symbol}

    ModuleInfo(name) = new(name, Dict(), Set())
end

name(data::ModuleInfo) = data.name
export!(data::ModuleInfo, sym::Symbol) = push!(data.exports, sym)
exports(data::ModuleInfo) = data.exports
set!(data::ModuleInfo, sym::Symbol, info::VarInfo) = data.globals[sym] = info
function lookup(data::ModuleInfo, sym::Symbol)::Nullable{VarInfo}
    if sym in keys(data.globals)
        return data.globals[sym]
    end

    # check standard library
    val = stdlibobject(sym)
    if !isnull(val)
        vi = VarInfo(UNKNOWN_LOCATION, Typeof(get(val)))
        info!(vi, StandardLibraryObject(get(val)))
        return vi
    end

    return Nullable{VarInfo}()
end

immutable MethodInfo <: AdditionalVarInfo
    # signature :: ...
    location :: Location
    body     :: Any
    isstaged :: Bool
end
location(mi::MethodInfo) = mi.location

"""
The binding is known to reference a standard library object. The "standard
library" consists of `Core`, `Base`, `Compat`, and their submodules.
"""
immutable StandardLibraryObject <: AdditionalVarInfo
    object :: Any
end
# TODO: remove {typeof(x.object)} part when #21397 fixed
extractobject(x::StandardLibraryObject) =
    Nullable{typeof(x.object)}(x.object)

# TODO: currently, this is not actually used
immutable FunctionInfo <: AdditionalVarInfo
    name    :: Symbol
    methods :: Vector{MethodInfo}
end
name(data::FunctionInfo) = data.name
method!(data::FunctionInfo, mi::MethodInfo) = push!(data.methods, mi)

type PragmaInfo
    location :: Location
    used     :: Bool
end

@compat abstract type _LintContext end
istoplevel(ctx::_LintContext) = false
toplevel(ctx::_LintContext) = istoplevel(ctx) ? ctx : toplevel(parent(ctx))
pragmas(ctx::_LintContext) = Dict{String, PragmaInfo}()
function pragma!(ctx::_LintContext, pragma, location::Location)
    pragmas(ctx)[pragma] = PragmaInfo(location, false)
end
finish(ctx::_LintContext, _) = nothing
globalset!(ctx::_LintContext, sym::Symbol, info::VarInfo) =
    set!(ctx, sym, info)
localset!(ctx::_LintContext, sym::Symbol, info::VarInfo) =
    set!(ctx, sym, info)

# A special context for linting a `module` keyword
immutable ModuleContext <: _LintContext
    parent        :: Nullable{_LintContext}
    data          :: ModuleInfo
    pragmas       :: Dict{String, PragmaInfo}

    """
    Methods whose linting has been deferred until the completion of this
    context.
    """
    deferred      :: Vector{MethodInfo}

    ModuleContext(parent, data) = new(parent, data, Dict(), [])
end

isroot(mctx::ModuleContext) = isnull(mctx.parent)
pragmas(mctx::ModuleContext) = mctx.pragmas
parent(mctx::ModuleContext) = get(mctx.parent)
data(mctx::ModuleContext) = mctx.data
lookup(mctx::ModuleContext, args...; kwargs...) =
    lookup(mctx.data, args...; kwargs...)
locallookup(mctx::ModuleContext, name::Symbol) = Nullable()
set!(mctx::ModuleContext, sym::Symbol, info::VarInfo) =
    set!(mctx.data, sym, info)
function defer!(mctx::ModuleContext, mi::MethodInfo)
    push!(mctx.deferred, mi)
end
export!(mctx::ModuleContext, sym::Symbol) = export!(mctx.data, sym)
exports(mctx::ModuleContext) = exports(mctx.data)
istoplevel(mctx::ModuleContext) = true
function finish(ctx::ModuleContext, cursor)
    for x in keys(ctx.data.globals)
        vi = ctx.data.globals[x]
        if source(vi) ∉ [:imported, :used]  # allow imported/used bindings
            loc = location(vi)
            if !isnull(stdlibobject(x))
                msg(cursor, :I343, x, "global variable defined at $loc with same name as export from Base")
            end
        end
    end
    for method in ctx.deferred
        lintfunctionbody(cursor, method)
    end
end

type LocalContext <: _LintContext
    parent        :: _LintContext
    declglobs     :: Set{Symbol}
    localvars     :: Dict{Symbol, VarInfo}
    oosvars       :: Set{Symbol}
    pragmas       :: Dict{String, PragmaInfo}

    """
    Methods whose linting has been deferred until the completion of this
    context.
    """
    deferred      :: Vector{MethodInfo}
    LocalContext(parent) = new(parent, Set(), Dict(), Set(), Dict(), [])
end
parent(ctx::LocalContext) = ctx.parent
pragmas(ctx::LocalContext) = ctx.pragmas
function defer!(ctx::LocalContext, mi::MethodInfo)
    push!(ctx.deferred, mi)
end
function finish(ctx::LocalContext, cursor)
    for method in ctx.deferred
        lintfunctionbody(cursor, method)
    end
    tl = toplevel(ctx)
    nl = parent(ctx)
    for x in keys(ctx.localvars)
        loc = location(ctx.localvars[x])
        if usages(ctx.localvars[x]) == 0 && !startswith(string(x), "_")
            # TODO: a better line number
            msg(cursor, :I340, x, "unused local variable, defined at $loc")
        elseif !isnull(stdlibobject(x))
            msg(cursor, :I342, x, "local variable defined at $loc shadows export from Base")
        elseif !isnull(lookup(tl, x))
            msg(cursor, :I341, x, "local variable defined at $loc shadows global variable defined at $(location(get(lookup(tl, x))))")
        elseif !isnull(lookup(nl, x))
            msg(cursor, :I344, x, "local variable defined at $loc shadows local variable defined at $(location(get(lookup(nl, x))))")
        end
    end
end

function set!(s::LocalContext, name::Symbol, vi::VarInfo)
    # TODO: check if it's soft or hard local scope
    var = locallookup(s, name)
    if !isnull(var)
        # TODO: warn about type instability?
        get(var).typeactual = Union{get(var).typeactual, vi.typeactual}
    else
        localset!(s, name, vi)
    end
end

function localset!(s::LocalContext, name::Symbol, vi::VarInfo)
    # TODO: check if already set?
    s.localvars[name] = vi
end

function globalset!(s::LocalContext, name::Symbol, vi::VarInfo)
    # TODO: check if this global declaration is legal
    # TODO: remove from localvars if needed
    push!(s.declglobs, name)
    globalset!(parent(s), name, vi)
end

function locallookup(ctx::LocalContext, name::Symbol)::Nullable{VarInfo}
    if name in keys(ctx.localvars)
        return ctx.localvars[name]
    elseif name in ctx.declglobs
        return lookup(toplevel(ctx), name)
    else
        return locallookup(parent(ctx), name)
    end
end

function lookup(ctx::LocalContext, name::Symbol)::Nullable{VarInfo}
    var = locallookup(ctx, name)
    if isnull(var)
        lookup(toplevel(ctx), name)
    else
        var
    end
end

@auto_hash_equals immutable LintIgnore
    errorcode :: Symbol
    variable  :: String
end

const LINT_IGNORE_DEFAULT = [
    LintIgnore(:I341, ""),
    LintIgnore(:I342, ""),
    LintIgnore(:W651, "")
]

type LintContext
    file         :: String
    "Current line number."
    line         :: Int
    "Line number at which the current toplevel expression begins at."
    lineabs      :: Int
    scope        :: String # usually the function name
    path         :: String
    included     :: Array{AbstractString,1} # list of files included
    macrocallLvl :: Int
    quoteLvl     :: Int
    messages     :: Array{LintMessage, 1}
    versionreachable:: Function # VERSION -> true means this code is reachable by VERSION
    ignore       :: Array{LintIgnore, 1}
    ifdepth      :: Int
    current      :: _LintContext
    function LintContext()
        mdata = ModuleInfo(:Main)
        mctx = ModuleContext(Nullable(), mdata)
        new("none", 0, 1, "", ".", AbstractString[],
            0, 0, LintMessage[], _ -> true,
            copy(LINT_IGNORE_DEFAULT), 0, mctx)
    end
end
location(ctx::LintContext) = Location(ctx.file, ctx.line)
function location!(ctx::LintContext, loc::Location)
    ctx.file = file(loc)
    ctx.path = dirname(ctx.file)
    ctx.line = line(loc)
end

finish(cur::LintContext) = finish(cur.current, cur)

function LintContext(file::AbstractString; ignore::Array{LintIgnore, 1} = LintIgnore[])
    ctx = LintContext()
    append!(ctx.ignore, ignore)
    ctx.file = file
    if ispath(file)
        ctx.path = dirname(abspath(file))
    end
    return ctx
end

function withcontext(f, ctx::LintContext, temp::_LintContext)
    old = ctx.current
    ctx.current = temp
    f()
    finish(ctx.current, ctx)
    ctx.current = old
end

function lookup(ctx::LintContext, sym::Symbol)::Nullable{VarInfo}
    lookup(ctx.current, sym)
end

function msg(ctx::LintContext, code::Symbol, variable, str::AbstractString)
    variable = string(variable)
    m = LintMessage(location(ctx), code, ctx.scope, variable, str)
    # filter out messages to ignore
    i = max(findfirst(ctx.ignore, LintIgnore(code, variable)),
            findfirst(ctx.ignore, LintIgnore(code, "")))
    if i == 0
        push!(ctx.messages, m)
    end
end

function msg(ctx::LintContext, code::Symbol, str::AbstractString)
    msg(ctx, code, "", str)
end


text:
"""
    module TwoCollidingBalls

Simulates a ball that is first sliding, then rolling and then colliding
with another ball on a table.
"""
module TwoCollidingBalls

using ModiaLang
import Modia3D
using  Modia3D.ModiaInterface

vmatSolids = Modia3D.VisualMaterial(color="Red"  , transparency=0.0)
vmatTable  = Modia3D.VisualMaterial(color="Green", transparency=0.1)

LxTable  = 3.0
LyTable  = 0.5
LzTable  = 0.02
diameter = 0.06

collidingBalls = Model(
    world = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, 0, -1]),
                                          enableContactDetection=true)),

    table = Object3D(parent=:world,
                     translation=[0.0, 0.0, -LzTable/2],
                     feature=Solid(shape=Box(lengthX=LxTable, lengthY=LyTable, lengthZ=LzTable),
                                   solidMaterial="DryWood",
                                   visualMaterial=vmatTable,
                                   contactMaterial="BilliardTable",
                                   collision=true)),

    ball1 = Object3D(feature=Solid(shape=Sphere(diameter=diameter),
                                   solidMaterial="BilliardBall",
                                   visualMaterial=vmatSolids,
                                   collision=true)),
    joint1 = FreeMotion(obj1=:world, obj2=:ball1,
                        r=Var(init=[-1.3, 0.0, diameter/2]),
                        v=Var(init=[3.0, 0.0, 0.0]),
                        rot=Var(init=[pi/2, 0.0, 0.0])),

    ball2 = Object3D(feature=Solid(shape=Sphere(diameter=diameter),
                                   solidMaterial="BilliardBall",
                                   visualMaterial=vmatSolids,
                                   collision=true)),
    joint2 = FreeMotion(obj1=:world, obj2=:ball2,
                        r=Var(init=[0.0, 0.0, diameter/2]),
                        rot=Var(init=[pi/2, 0.0, 0.0]))
)

twoCollidingBalls = @instantiateModel(buildModia3D(collidingBalls), unitless=true, log=false, logStateSelection=false, logCode=false)

stopTime = 1.5
tolerance = 1e-8
requiredFinalStates = [0.4489974852019629, 4.971482124288075e-6, 0.029996861274764766, 0.5336096206429409, 4.775435383850775e-6, -8.566849073992856e-7, 1.5708013764495665, 1.2298784775382833e-5, -55.4441088791202, -6.408457808447407e-5, -9.894396552324208e-5, -17.78647176388158, 1.4032405177409832, -4.464837913922324e-6, 0.029997744288682248, 1.4383823207541608, -4.574293232521244e-6, -2.3103945635420407e-6, 1.570796533109225, 1.2115305686199998e-8, -44.07171252097457, 0.00015344214722251976, 2.5704782290946986e-6, -47.9463685248874]
simulate!(twoCollidingBalls, stopTime=stopTime, tolerance=tolerance, log=true, logStates=true, logEvents=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(twoCollidingBalls, ["joint1.r" "joint1.rot"; "joint1.v" "joint1.w"], figure=1)

end


text:
using GeometricIntegrators.Equations

ϑ₁(t, q) = g₁₁(t,q) * q[4] + A₁(t,q)
ϑ₂(t, q) = g₂₂(t,q) * q[5] + A₂(t,q)
ϑ₃(t, q) = g₃₃(t,q) * q[6] + A₃(t,q)

function ϑ(t, q, θ)
    θ[1] = ϑ₁(t,q)
    θ[2] = ϑ₂(t,q)
    θ[3] = ϑ₃(t,q)
    θ[4] = zero(eltype(q))
    θ[5] = zero(eltype(q))
    θ[6] = zero(eltype(q))
    nothing
end

function ω(t, q, Β)
    Β[1,1] = 0
    Β[1,2] = dϑ₁dx₂(t,q) - dϑ₂dx₁(t,q)
    Β[1,3] = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
    Β[1,4] = +1
    Β[1,5] = 0
    Β[1,6] = 0

    Β[2,1] = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)
    Β[2,2] = 0
    Β[2,3] = dϑ₂dx₃(t,q) - dϑ₃dx₂(t,q)
    Β[2,4] = 0
    Β[2,5] = +1
    Β[2,6] = 0

    Β[3,1] = dϑ₃dx₁(t,q) - dϑ₁dx₃(t,q)
    Β[3,2] = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
    Β[3,3] = 0
    Β[3,4] = 0
    Β[3,5] = 0
    Β[3,6] = +1

    Β[4,1] = -1
    Β[4,2] = 0
    Β[4,3] = 0
    Β[4,4] = 0
    Β[4,5] = 0
    Β[4,6] = 0

    Β[5,1] = 0
    Β[5,2] = -1
    Β[5,3] = 0
    Β[5,4] = 0
    Β[5,5] = 0
    Β[5,6] = 0

    Β[6,1] = 0
    Β[6,2] = 0
    Β[6,3] = -1
    Β[6,4] = 0
    Β[6,5] = 0
    Β[6,6] = 0

    nothing
end

function dϑ(t, q, dϑ)
    dϑ[1,1] = dϑ₁dx₁(t,q)
    dϑ[1,2] = dϑ₁dx₂(t,q)
    dϑ[1,3] = dϑ₁dx₃(t,q)
    dϑ[1,4] = 0
    dϑ[1,5] = 0
    dϑ[1,6] = 0

    dϑ[2,1] = dϑ₂dx₁(t,q)
    dϑ[2,2] = dϑ₂dx₂(t,q)
    dϑ[2,3] = dϑ₂dx₃(t,q)
    dϑ[2,4] = 0
    dϑ[2,5] = 0
    dϑ[2,6] = 0

    dϑ[3,1] = dϑ₃dx₁(t,q)
    dϑ[3,2] = dϑ₃dx₂(t,q)
    dϑ[3,3] = dϑ₃dx₃(t,q)
    dϑ[3,4] = 0
    dϑ[3,5] = 0
    dϑ[3,6] = 0

    dϑ[4,1] = 0
    dϑ[4,2] = 0
    dϑ[4,3] = 0
    dϑ[4,4] = 0
    dϑ[4,5] = 0
    dϑ[4,6] = 0

    dϑ[5,1] = 0
    dϑ[5,2] = 0
    dϑ[5,3] = 0
    dϑ[5,4] = 0
    dϑ[5,5] = 0
    dϑ[5,6] = 0

    dϑ[6,1] = 0
    dϑ[6,2] = 0
    dϑ[6,3] = 0
    dϑ[6,4] = 0
    dϑ[6,5] = 0
    dϑ[6,6] = 0

    nothing
end

β₁(t,q) = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
β₂(t,q) = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
β₃(t,q) = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)

β(t,q) = sqrt(β₁(t,q)^2 + β₂(t,q)^2 + β₃(t,q)^2)

hamiltonian(t,q) = 0.5 * (g₁₁(t,q) * q[4]^2 + g₂₂(t,q) * q[5]^2 + g₃₃(t,q) * q[6]^2)
hamiltonian(t,q,p) = hamiltonian(t,q)

toroidal_momentum(t,q) = ϑ₃(t,q)

dHdx₁(t, q) = zero(eltype(q))
dHdx₂(t, q) = zero(eltype(q))
dHdx₃(t, q) = zero(eltype(q))
dHdx₄(t, q) = q[4]
dHdx₅(t, q) = q[5]
dHdx₆(t, q) = q[6]

function dH(t, q, dH)
    dH[1] = dHdx₁(t, q)
    dH[2] = dHdx₂(t, q)
    dH[3] = dHdx₃(t, q)
    dH[4] = dHdx₄(t, q)
    dH[5] = dHdx₅(t, q)
    dH[6] = dHdx₆(t, q)
    nothing
end

f₁(t, q, v) = q[4]
f₂(t, q, v) = q[5]
f₃(t, q, v) = q[6]
f₄(t, q, v) = q[5] * B₃(t,q) - q[6] * B₂(t,q)
f₅(t, q, v) = q[6] * B₁(t,q) - q[4] * B₃(t,q)
f₆(t, q, v) = q[4] * B₂(t,q) - q[5] * B₁(t,q)

function charged_particle_3d_periodicity(qᵢ)
    periodicity = zeros(eltype(qᵢ), size(qᵢ,1))
    periodicity[3] = 2π
    return periodicity
end

function charged_particle_3d_pᵢ(qᵢ, tᵢ=0)
    pᵢ = zero(qᵢ)

    if ndims(qᵢ) == 1
        ϑ(tᵢ, qᵢ, pᵢ)
    else
        for i in 1:size(qᵢ,2)
            @views ϑ(tᵢ, qᵢ[:,i], pᵢ[:,i])
        end
    end
    pᵢ
end

function charged_particle_3d_v(t, q, v)
    v[1] = f₁(t, q, v)
    v[2] = f₂(t, q, v)
    v[3] = f₃(t, q, v)
    v[4] = f₄(t, q, v)
    v[5] = f₅(t, q, v)
    v[6] = f₅(t, q, v)

    nothing
end

charged_particle_3d_iode_ϑ(t, q, v, p) = ϑ(t, q, p)

function charged_particle_3d_iode_f(t, q, v, f)
    f[1] = dA₁dx₁(t,q) * v[1] + dA₂dx₁(t,q) * v[2] + dA₃dx₁(t,q) * v[3] + E₁(t,q)
    f[2] = dA₁dx₂(t,q) * v[1] + dA₂dx₂(t,q) * v[2] + dA₃dx₂(t,q) * v[3] + E₂(t,q)
    f[3] = dA₁dx₃(t,q) * v[1] + dA₂dx₃(t,q) * v[2] + dA₃dx₃(t,q) * v[3] + E₃(t,q)
    f[4] = v[1] - q[4]
    f[5] = v[2] - q[5]
    f[6] = v[3] - q[6]
    nothing
end

function charged_particle_3d_iode_g(t, q, λ, g)
    g[1] = dA₁dx₁(t,q) * λ[1] + dA₂dx₁(t,q) * λ[2] + dA₃dx₁(t,q) * λ[3]
    g[2] = dA₁dx₂(t,q) * λ[1] + dA₂dx₂(t,q) * λ[2] + dA₃dx₂(t,q) * λ[3]
    g[3] = dA₁dx₃(t,q) * λ[1] + dA₂dx₃(t,q) * λ[2] + dA₃dx₃(t,q) * λ[3]
    g[4] = λ[1]
    g[5] = λ[2]
    g[6] = λ[3]
    nothing
end

charged_particle_3d_v(t, q, p, v) = charged_particle_3d_v(t, q, v)

function charged_particle_3d_sode_fx(t, q::Vector{DT}, f::Vector{DT}, Δt) where {DT}
    x = q[1:3]
    v = q[4:6]

    f .= vcat(x .+ Δt .* v, v)

    nothing
end

function charged_particle_3d_sode_fv(t, q::Vector{DT}, f::Vector{DT}, Δt) where {DT}
    x = q[1:3]
    v = q[4:6]

    local lB₁ = B₁(t-Δt, x) * x[1]
    local lB₂ = B₂(t-Δt, x) * x[1]
    local lB₃ = B₃(t-Δt, x) / x[1]

    local Ω = zeros(DT,3,3)

    Ω[1,2] = - lB₃
    Ω[1,3] = + lB₂
    Ω[2,1] = + lB₃
    Ω[2,3] = + lB₁
    Ω[3,1] = - lB₂
    Ω[3,2] = - lB₁

    C1 = eye(DT,3) .+ (0.5*Δt) .* Ω
    C2 = eye(DT,3) .- (0.5*Δt) .* Ω
    R  = inv(C1) * C2

    f .= vcat(x, R * v)

    nothing
end

function charged_particle_3d_ode(qᵢ=qᵢ; periodic=true)
    if periodic
        ODE(charged_particle_3d_v, qᵢ; periodicity=charged_particle_3d_periodicity(qᵢ))
    else
        ODE(charged_particle_3d_v, qᵢ)
    end
end

function charged_particle_3d_sode(qᵢ=qᵢ; periodic=true)
    if periodic
        SODE((charged_particle_3d_sode_fv, charged_particle_3d_sode_fx), qᵢ; periodicity=charged_particle_3d_periodicity(qᵢ))
    else
        SODE((charged_particle_3d_sode_fv, charged_particle_3d_sode_fx), qᵢ)
    end
end

function charged_particle_3d_iode(qᵢ=qᵢ)
    IODE(charged_particle_3d_iode_ϑ, charged_particle_3d_iode_f,
            charged_particle_3d_iode_g, qᵢ, charged_particle_3d_pᵢ(qᵢ);
            v=charged_particle_3d_v, h=hamiltonian)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related with True Equator Mean Equinox (TEME) reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Remarks
# ==============================================================================
#
# As mentioned in [1](p. 233), there is not an official definition of the TEME
# frame. Hence, in this package, it is considered the definition presented in
# [1](p. 233) in which the complete form of the Equation of Equinoxes is used.
# This seems to be the case when comparing the values shown in Table 3-6 [1](p.
# 232).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export r_teme_to_tod,  r_tod_to_teme
export r_teme_to_mod,  r_mod_to_teme
export r_teme_to_gcrf, r_gcrf_to_teme
export r_teme_to_pef,  r_pef_to_teme

#                                 TEME <=> TOD
# ==============================================================================

"""
    r_teme_to_tod([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the True of Date (TOD) frame at the Julian Day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-76/FK5 theory and TEME definition in **[1]**(p. 233).
Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) [rad] and in longitude (`δΔψ_1980`) [rad] that are usually obtained
from IERS EOP Data (see [`get_iers_eop`](@ref)).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the TOD frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function r_teme_to_tod(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0)
    return r_teme_to_tod(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_teme_to_tod(
    T::Type,
    JD_TT::Number,
    δΔϵ_1980::Number = 0,
    δΔψ_1980::Number = 0
)
    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    mϵ_1980, Δϵ_1980, Δψ_1980 = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δϵ_1980 += δΔϵ_1980
    Δψ_1980 += δΔψ_1980

    # Compute the obliquity.
    ϵ_1980 = mϵ_1980 + Δϵ_1980

    # Evaluate the Delaunay parameters associated with the Moon in the interval
    # [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    T_TT = (JD_TT - JD_J2000) / 36525
    r    = 360
    Ω_m  = @evalpoly(
        T_TT,
        125.04452222,
        -5r - 134.1362608,
        0.0020708,
        2.2e-6
    )
    Ω_m = mod(Ω_m, 360) * π / 180

    # Compute the equation of Equinoxes.
    #
    # According to [2], the constant unit before `sin(2Ω_m)` is also in [rad].
    Eq_equinox1982 = Δψ_1980*cos(mϵ_1980) +
        (0.002640sin(1Ω_m) + 0.000063sin(2Ω_m)) * π / 648000

    # Compute the rotation.
    return angle_to_rot(T, -Eq_equinox1982, 0, 0, :ZYX)
end

"""
    r_tod_to_teme([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True of Date (TOD) frame with the True
Equator Mean Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-76/FK5 theory and TEME definition in **[1]**(p.
233). Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) [rad] and in longitude (`δΔψ_1980`) [rad] that are usually obtained
from IERS EOP Data (see [`get_iers_eop`](@ref)).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TOD frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function r_tod_to_teme(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0)
    return r_tod_to_teme(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_tod_to_teme(
    T::Type,
    JD_TT::Number,
    δΔϵ_1980::Number = 0,
    δΔψ_1980::Number = 0
)
    return inv_rotation(r_teme_to_tod(T, JD_TT, δΔϵ_1980, δΔψ_1980))
end

#                                 TEME <=> MOD
# ==============================================================================

"""
    r_teme_to_mod([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the Mean of Date (MOD) frame at the Julian Day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-76/FK5 theory and TEME definition in **[1]**(p. 233).
Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) [rad] and in longitude (`δΔψ_1980`) [rad] that are usually obtained
from IERS EOP Data (see [`get_iers_eop`](@ref)).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the MOD frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function r_teme_to_mod(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0)
    return r_teme_to_mod(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_teme_to_mod(
    T::Type,
    JD_TT::Number,
    δΔϵ_1980::Number = 0,
    δΔψ_1980::Number = 0
)
    # Notice that, in this case, we will not use `r_teme_to_tod`, and `rTODtoMOD`
    # because this would call the function `nutation` twice, leading to a huge
    # performance drop. Hence, the code of those two functions is almost
    # entirely rewritten here.

    # Compute the nutation in the Julian Day (Terrestrial Time) `JD_TT`.
    mϵ_1980, Δϵ_1980, Δψ_1980 = nutation_fk5(JD_TT)

    # Add the corrections to the nutation in obliquity and longitude.
    Δϵ_1980 += δΔϵ_1980
    Δψ_1980 += δΔψ_1980

    # Compute the obliquity.
    ϵ_1980 = mϵ_1980 + Δϵ_1980

    # Evaluate the Delaunay parameters associated with the Moon in the interval
    # [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    T_TT = (JD_TT - JD_J2000) / 36525
    r    = 360
    Ω_m  = @evalpoly(
        T_TT,
        125.04452222,
        -5r - 134.1362608,
        0.0020708,
        2.2e-6
    )
    Ω_m = mod(Ω_m, 360) * π / 180

    # Compute the equation of Equinoxes.
    #
    # According to [2], the constant unit before `sin(2Ω_m)` is also in [rad].
    Eq_equinox1982 = Δψ_1980*cos(mϵ_1980) +
        (0.002640sin(1Ω_m) + 0.000063sin(2Ω_m)) * π / 648000

    # Compute the rotation TEME => TOD.
    r_TOD_TEME = angle_to_rot(T, -Eq_equinox1982, 0, 0, :ZYX)

    # Compute the rotation TOD => MOD.
    r_MOD_TOD = angle_to_rot(T, ϵ_1980, Δψ_1980, -mϵ_1980, :XZX)

    # Return the full rotation.
    return compose_rotation(r_TOD_TEME, r_MOD_TOD)
end

"""
    r_mod_to_teme([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the Mean of Date (MOD) frame with the True
Equator Mean Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-76/FK5 theory and TEME definition in **[1]**(p.
233). Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) [rad] and in longitude (`δΔψ_1980`) [rad] that are usually
obtained from IERS EOP Data (see [`get_iers_eop`](@ref)).  .

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the MOD frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function r_mod_to_teme(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0)
    return r_mod_to_teme(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_mod_to_teme(
    T::Type,
    JD_TT::Number,
    δΔϵ_1980::Number = 0,
    δΔψ_1980::Number = 0
)
    return inv_rotation(r_teme_to_mod(T, JD_TT, δΔϵ_1980, δΔψ_1980))
end

#                                TEME <=> GCRF
# ==============================================================================

"""
    r_teme_to_gcrf([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the Geocentric Celestial Reference Frame (GCRF) at the Julian Day `JD_TT`
[Terrestrial Time]. This algorithm uses the IAU-76/FK5 theory and TEME
definition in **[1]**(p. 233). Notice that one can provide corrections for the
nutation in obliquity (`δΔϵ_1980`) [rad] and in longitude (`δΔψ_1980`) [rad]
that are usually obtained from IERS EOP Data (see [`get_iers_eop`](@ref)).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    The EOP data related to the nutation of the obliquity (`δΔϵ_1980`) and the
    nutation of the longitude (`δΔψ_1980`) can be omitted. In this case, the
    GCRF frame is what is usually called J2000 reference frame.

# Returns

The rotation that aligns the TEME frame with the GCRF frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function r_teme_to_gcrf(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0)
    return r_teme_to_gcrf(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_teme_to_gcrf(
    T::Type,
    JD_TT::Number,
    δΔϵ_1980::Number = 0,
    δΔψ_1980::Number = 0
)
    # Compute the rotation TEME => MOD.
    r_MOD_TEME  = r_teme_to_mod(T, JD_TT, δΔϵ_1980, δΔψ_1980)

    # Compute the rotation MOD => GCRF.
    r_GCRF_MOD = r_mod_to_gcrf_fk5(T, JD_TT)

    # Return the full rotation.
    return compose_rotation(r_MOD_TEME, r_GCRF_MOD)
end

"""
    r_gcrf_to_teme([T,] JD_TT::Number [, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0])

Compute the rotation that aligns the GCRF frame with the True Equator Mean
Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial Time]. This
algorithm uses the IAU-76/FK5 theory and TEME definition in **[1]**(p. 233).
Notice that one can provide corrections for the nutation in obliquity
(`δΔϵ_1980`) [rad] and in longitude (`δΔψ_1980`) [rad] that are usually obtained
from IERS EOP Data (see [`get_iers_eop`](@ref)).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

!!! info
    The EOP data related to the nutation of the obliquity (`δΔϵ_1980`) and the
    nutation of the longitude (`δΔψ_1980`) can be omitted. In this case, the
    GCRF frame is what is usually called J2000 reference frame.

# Returns

The rotation that aligns the GCRF frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function r_gcrf_to_teme(JD_TT::Number, δΔϵ_1980::Number = 0, δΔψ_1980::Number = 0)
    return r_gcrf_to_teme(DCM, JD_TT, δΔϵ_1980, δΔψ_1980)
end

function r_gcrf_to_teme(
    T::Type,
    JD_TT::Number,
    δΔϵ_1980::Number = 0,
    δΔψ_1980::Number = 0
)
    return inv_rotation(r_teme_to_gcrf(T, JD_TT, δΔϵ_1980, δΔψ_1980))
end

#                                 TEME <=> PEF
# ==============================================================================

"""
    r_teme_to_pef([T,] JD_TT::Number)

Compute the rotation that aligns the True Equator Mean Equinox (TEME) frame with
the Pseudo-Earth Fixed (PEF) frame at the Julian Day `JD_TT` [Terrestrial Time].
This algorithm uses the IAU-76/FK5 theory and TEME definition in **[1]**(p.
233).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the TEME frame with the PEF frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
r_teme_to_pef(JD_UT1::Number) = r_teme_to_pef(DCM, JD_UT1)

function r_teme_to_pef(T::Type, JD_UT1::Number)
    # Compute the Greenwich Mean Sidereal Time.
    θ_gmst = jd_to_gmst(JD_UT1)

    # Compute the rotation.
    return angle_to_rot(T, θ_gmst, 0, 0, :ZYX)
end

"""
    r_pef_to_teme([T,] JD_TT::Number)

Compute the rotation that aligns the Pseudo-Earth Fixed (PEF) frame with the
True Equator Mean Equinox (TEME) frame at the Julian Day `JD_TT` [Terrestrial
Time]. This algorithm uses the IAU-76/FK5 theory and TEME definition in
**[1]**(p. 233).

The rotation type is described by the optional variable `T`. If it is `DCM`,
then a DCM will be returned. Otherwise, if it is `Quaternion`, then a Quaternion
will be returned. In case this parameter is omitted, then it falls back to
`DCM`.

# Returns

The rotation that aligns the PEF frame with the TEME frame. The rotation
representation is selected by the optional parameter `T`.

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
r_pef_to_teme(JD_UT1::Number) = r_pef_to_teme(DCM, JD_UT1)
r_pef_to_teme(T::Type, JD_UT1::Number) = inv_rotation(r_teme_to_pef(T, JD_UT1))


text:
@testset "310.minimum-height-trees.jl" begin
    @test find_min_height_trees(1, Vector{Int}[]) == [1]
    @test find_min_height_trees(2, [[1, 2]]) == [1, 2]
    @test find_min_height_trees(6, [[4, 1], [4, 2], [4, 3], [4, 5], [6, 5]]) == [4, 5]
    @test find_min_height_trees(4, [[1, 2], [3, 2], [2, 4]]) == [2]
end


text:
# include all unit testing files

include("test_texpr.jl")
include("test_tmode.jl")

include("test_scalar_backend.jl")
#include("test_extensions.jl")


text:
function enableDac(rhd::FPGA,dacChannel::Int,enabled::Bool)

    #error checking goes here

    if rhd.usb3
        if (OPEN_EPHYS)
            dacEnMask = 0x0400
        else
            dacEnMask = 0x0800
        end
    else
        dacEnMask = 0x0200
    end

    if dacChannel == 0
        SetWireInValue(rhd,WireInDacSource1,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 1
        SetWireInValue(rhd,WireInDacSource2,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 2
        SetWireInValue(rhd,WireInDacSource3,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 3
        SetWireInValue(rhd,WireInDacSource4,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 4
        SetWireInValue(rhd,WireInDacSource5,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 5
        SetWireInValue(rhd,WireInDacSource6,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 6
        SetWireInValue(rhd,WireInDacSource7,(enabled ? dacEnMask : 0x0000), dacEnMask)
    elseif dacChannel == 7
        SetWireInValue(rhd,WireInDacSource8,(enabled ? dacEnMask : 0x0000), dacEnMask)
    end

    UpdateWireIns(rhd)

    nothing
end

function selectDacDataStream(rhd::FPGA,dacChannel, stream)
    #error checking goes here

    if rhd.usb3
        if (OPEN_EPHYS)
            dacStreamMask = 0x03e0
        else
            dacStreamMask = 0x07e0
        end
    else
        dacStreamMask = 0x01e0
    end

    if dacChannel == 0
         SetWireInValue(rhd,WireInDacSource1, stream << 5, dacStreamMask)
    elseif dacChannel == 1
         SetWireInValue(rhd,WireInDacSource2, stream << 5, dacStreamMask)
    elseif dacChannel == 2
         SetWireInValue(rhd,WireInDacSource3, stream << 5, dacStreamMask)
    elseif dacChannel == 3
         SetWireInValue(rhd,WireInDacSource4, stream << 5, dacStreamMask)
    elseif dacChannel == 4
         SetWireInValue(rhd,WireInDacSource5, stream << 5, dacStreamMask)
    elseif dacChannel == 5
         SetWireInValue(rhd,WireInDacSource6, stream << 5, dacStreamMask)
    elseif dacChannel == 6
         SetWireInValue(rhd,WireInDacSource7, stream << 5, dacStreamMask)
    elseif dacChannel == 7
         SetWireInValue(rhd,WireInDacSource8, stream << 5, dacStreamMask)
    end

    UpdateWireIns(rhd)

    nothing
end

function selectDacDataChannel(rhd::FPGA,dacChannel::Int, dataChannel)
    #error checking goes here

    if dacChannel == 0
        SetWireInValue(rhd,WireInDacSource1,dataChannel << 0, 0x001f)
    elseif dacChannel == 1
        SetWireInValue(rhd,WireInDacSource2,dataChannel << 0, 0x001f)
    elseif dacChannel == 2
        SetWireInValue(rhd,WireInDacSource3,dataChannel << 0, 0x001f)
    elseif dacChannel == 3
        SetWireInValue(rhd,WireInDacSource4,dataChannel << 0, 0x001f)
    elseif dacChannel == 4
        SetWireInValue(rhd,WireInDacSource5,dataChannel << 0, 0x001f)
    elseif dacChannel == 5
        SetWireInValue(rhd,WireInDacSource6,dataChannel << 0, 0x001f)
    elseif dacChannel == 6
        SetWireInValue(rhd,WireInDacSource7,dataChannel << 0, 0x001f)
    elseif dacChannel == 7
        SetWireInValue(rhd,WireInDacSource8,dataChannel << 0, 0x001f)
    end

    UpdateWireIns(rhd)

    nothing
end

#Sets DAC to 16 bit value.
setDacManual(rhd::FPGA,value)=(SetWireInValue(rhd,WireInManualTriggers,value<<15,0xffff0000);UpdateWireIns(rhd))

setDacGain(rhd::FPGA,gain)=(SetWireInValue(rhd,WireInResetRun,gain<<13,0xe000);UpdateWireIns(rhd))

setAudioNoiseSuppress(rhd::FPGA,noiseSuppress)=(SetWireInValue(rhd,WireInResetRun,noiseSuppress<<6,0x1fc0);UpdateWireIns(rhd))

function setDacThreshold(rhd::FPGA,dacChannel, threshold, trigPolarity)

    #error checking goes here

    #Set threshold level
    SetWireInValue(rhd,WireInMultiUse,threshold)
    UpdateWireIns(rhd)
    if (OPEN_EPHYS)
        ActivateTriggerIn(rhd,TrigInDacThresh, dacChannel)
    else
        ActivateTriggerIn(rhd,TrigInDacConfig, dacChannel)
    end

    #Set threshold polarity
    SetWireInValue(rhd,WireInMultiUse, (trigPolarity ? 1 : 0))
    UpdateWireIns(rhd)
    if (OPEN_EPHYS)
        ActivateTriggerIn(rhd,TrigInDacThresh, dacChannel+8)
    else
        ActivateTriggerIn(rhd,TrigInDacConfig, dacChannel+8)
    end

    nothing
end


text:
# originates from sisue #458

# using Revise

using Test
using IncrementalInference

@testset "Test clique factors, #458 Example 1" begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)
addVariable!(fg, :x4, ContinuousScalar)

addVariable!(fg, :l0, ContinuousScalar)
addVariable!(fg, :l1, ContinuousScalar)

lc = LinearRelative(Normal())
lp = Prior(Normal())
addFactor!(fg, [:x0;:x1], lc, graphinit=false)
addFactor!(fg, [:x1;:x2], lc, graphinit=false)
addFactor!(fg, [:x2;:x3], lc, graphinit=false)
addFactor!(fg, [:x3;:x4], lc, graphinit=false)

addFactor!(fg, [:x0;:l0], lc, graphinit=false)
addFactor!(fg, [:x2;:l0], lc, graphinit=false)

addFactor!(fg, [:x0;:l1], lc, graphinit=false)
addFactor!(fg, [:x2;:l1], lc, graphinit=false)

addFactor!(fg, [:x0;], lp, graphinit=false)
addFactor!(fg, [:l0;], lp, graphinit=false)

# build tree
vo = Symbol[:x2, :x0, :l0, :x3, :x1, :l1, :x4]
tree = buildTreeReset!(fg, vo)
# drawTree(tree, show=true)

# check that frontal variables only show up once
frontals = getCliqFrontalVarIds.((x->getClique(tree, x)).([:x0; :l0; :x4]))

@test intersect(frontals[1], frontals[2]) |> length == 0
@test intersect(frontals[2], frontals[3]) |> length == 0
@test intersect(frontals[1], frontals[3]) |> length == 0

# check that all variables exist as frontals
lsvars = ls(fg)
@test intersect(union(frontals...), lsvars) |> length == lsvars |> length

## Now check if factors in cliques are okay

C3 = getClique(tree, :x0)
C3_fg = buildCliqSubgraph(fg, C3)
# drawGraph(C3_fg, show=true)

C3_fcts = [:x0l0f1;:x0l1f1;:x0x1f1;:x0f1]

@test intersect(ls(C3_fg), [:x0; :x1; :l0; :l1]) |> length == 4
@test intersect(lsf(C3_fg), C3_fcts) |> length == length(C3_fcts)
@test intersect(getCliqFactorIdsAll(C3), C3_fcts) |> length == length(C3_fcts)

C2 = getClique(tree, :l0)
C2_fg = buildCliqSubgraph(fg, C2)
# drawGraph(C2_fg, show=true)

C2_fcts = [:x1x2f1; :x2x3f1; :x2l0f1; :x2l1f1; :l0f1]

@test intersect(ls(C2_fg), [:x3; :x2; :x1; :l0; :l1]) |> length == 5
@test intersect(lsf(C2_fg), C2_fcts) |> length == length(C2_fcts)
@test intersect(getCliqFactorIdsAll(C2), C2_fcts) |> length == length(C2_fcts)

C1 = getClique(tree, :x4)
C1_fg = buildCliqSubgraph(fg, C1)
# drawGraph(C1_fg, show=true)

C1_fcts = [:x3x4f1;]

@test intersect(ls(C1_fg), [:x3; :x4; :x1; :l1]) |> length == 4
@test intersect(lsf(C1_fg), C1_fcts) |> length == length(C1_fcts)
@test intersect(getCliqFactorIdsAll(C1), C1_fcts) |> length == length(C1_fcts)

# check that all factors are counted
allCliqFcts = union(C1_fcts, C2_fcts, C3_fcts)

@test length(intersect(lsf(fg), allCliqFcts)) == length(allCliqFcts)

end

@testset "Test clique factors, #458 Example 2" begin

@warn "Test for Example 2 from 458 must still be coded."

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)

addVariable!(fg, :lm0, ContinuousScalar)
addVariable!(fg, :lm3, ContinuousScalar)

lc = LinearRelative(Normal())
lp = Prior(Normal())
addFactor!(fg, [:x0;:x1], lc, graphinit=false)
addFactor!(fg, [:x1;:x2], lc, graphinit=false)
addFactor!(fg, [:x2;:x3], lc, graphinit=false)

addFactor!(fg, [:x0;:lm0], lc, graphinit=false)
addFactor!(fg, [:x1;:lm0], lc, graphinit=false)
addFactor!(fg, [:x2;:lm3], lc, graphinit=false)
addFactor!(fg, [:x3;:lm3], lc, graphinit=false)

# particular order from 458
vo = Symbol[:x0, :x2, :x1, :lm3, :lm0, :x3]
tree = buildTreeReset!(fg, vo)
# drawTree(tree, show=true)

sfg_1 = buildCliqSubgraph(fg,tree,:x3)
sfg_2 = buildCliqSubgraph(fg,tree,:x2)
sfg_3 = buildCliqSubgraph(fg,tree,:x0)
# drawGraph(sfg_1, show=true)
# drawGraph(sfg_2, show=true)
# drawGraph(sfg_3, show=true)

## WIP
C1_fcts = [:x1lm0f1; :x3lm3f1]
C2_fcts = [:x1x2f1; :x2x3f1; :x2lm3f1]
C3_fcts = [:x0x1f1; :x0lm0f1]

# check all factors are accounted for
@test union(C1_fcts, C2_fcts, C3_fcts) |> length == lsf(fg) |> length
@test intersect(C1_fcts, C2_fcts) |> length == 0
@test intersect(C1_fcts, C3_fcts) |> length == 0
@test intersect(C2_fcts, C3_fcts) |> length == 0

# clique 1
@test intersect(getCliqFactorIdsAll(tree, :x3), C1_fcts) |> length == 2
@test intersect(lsf(sfg_1), C1_fcts) |> length == 2

# clique 2
@test intersect( getCliqFactorIdsAll(tree, :x2), C2_fcts) |> length == 3
@test intersect( lsf(sfg_2), C2_fcts) |> length == 3

# clique 3
@test intersect( getCliqFactorIdsAll(tree, :x0), C3_fcts) |> length == 2
@test intersect( lsf(sfg_3), C3_fcts) |> length == 2

end

@testset "Test cliqueSubgraph frontals, seperators, potentials" begin
# clique subfg test
fg = generateCanonicalFG_lineStep(4, landmarkPriorsAt=[0,4])
# dfgplot(fg)
tree = buildTreeReset!(fg)

#check clique 1
cliq = getClique(tree, 1)

cliqfron = [:x0, :lm0, :x2]
cliqsep = Symbol[]
cliqfacs = [:lm0f1, :x0x2f1, :x0lm0f1, :x0f1, :x2lm0f1]

@test issetequal(getCliqFrontalVarIds(cliq), cliqfron)
@test issetequal(getCliqSeparatorVarIds(cliq), cliqsep)
@test issetequal(getCliqFactorIdsAll(cliq), cliqfacs)

sfg = buildCliqSubgraph(fg, cliq)
@test issetequal(ls(sfg), union(cliqsep,  cliqfron))
@test issetequal(lsf(sfg), cliqfacs)

# check clique 2
cliq = getClique(tree, 2)

cliqfron = [:x4, :lm4]
cliqsep = [:x2]
cliqfacs = [:x2lm4f1, :x2x4f1, :x4lm4f1, :lm4f1]

@test issetequal(getCliqFrontalVarIds(cliq), cliqfron)
@test issetequal(getCliqSeparatorVarIds(cliq), cliqsep)
@test issetequal(getCliqFactorIdsAll(cliq), cliqfacs)

sfg = buildCliqSubgraph(fg, cliq)
@test issetequal(ls(sfg), union(cliqsep,  cliqfron))
@test issetequal(lsf(sfg), cliqfacs)

sfg = buildCliqSubgraph(fg, tree, :x2)
# dfgplot(sfg)
end
#


text:
export
   AbstractCTQW,
   CTQW,
   matrix

"""
    AbstractCTQW

Abstract CTQW model. By default evolve according to Schrödinger equation and
performs measurmenet by taking square of absolute values of its elements. Default
representation of `AbstractCTQW` is `CTQW`.
"""
abstract type AbstractCTQW <: QWModelCont end

"""
    CTQW(graph[, matrix])

Default representation of `AbstractCTQW`. `matrix` defaults
to `:adjacency`. The Hamiltonian is a sparse matrix.
"""
struct CTQW <: AbstractCTQW
   graph::Graph
   matrix::Symbol
   CTQW(graph::Graph, matrix::Symbol=:adjacency) = new(graph, matrix)
end

"""
    CTQWDense(graph[, matrix])

Alternative representation of `AbstractCTQW`. `matrix` defaults
to `:adjacency`. The Hamiltonian is a dense matrix.
"""
struct CTQWDense <: AbstractCTQW
   graph::Graph
   matrix::Symbol
   CTQWDense(graph::Graph, matrix::Symbol=:adjacency) = new(graph, matrix)
end

"""
    matrix(ctqw)

Returns the matrix symbol defining matrix graph used.
"""
matrix(ctqw::AbstractCTQW) = ctqw.matrix

"""
    proj(type, i, n)

Return a canonic projection onto `i`-th subspace of `n`-dimensional space.
 """
function proj(::Type{T}, i::Int, n::Int) where T<:Number
   result = spzeros(T, n, n)
   result[i,i] = 1
   result
end

"""
    check_qwmodel(ctqw, parameters)

Private functions which checks the existance of `:hamiltonian`, its type and
dimensionality.
"""
function check_qwmodel(ctqw::AbstractCTQW,
                    parameters::Dict{Symbol})
   @assert :hamiltonian ∈ keys(parameters) "parameters needs to have key hamiltonian"
   @assert isa(parameters[:hamiltonian], AbstractMatrix{<:Number}) "value for :hamiltonian needs to be AbstractMatrix with numbers"
   @assert size(parameters[:hamiltonian], 1) == size(parameters[:hamiltonian], 2) == nv(ctqw.graph) "Hamiltonian needs to be square matrix of order equal to graph order"
   nothing
end

function check_qwmodel(ctqw::CTQWDense,
                    parameters::Dict{Symbol})
   @assert :hamiltonian ∈ keys(parameters) "parameters needs to have key hamiltonian"
   @assert isa(parameters[:hamiltonian], DenseMatrix{<:Number}) "value for :hamiltonian needs to be AbstractMatrix with numbers"
   @assert size(parameters[:hamiltonian], 1) == size(parameters[:hamiltonian], 2) == nv(ctqw.graph) "Hamiltonian needs to be square matrix of order equal to graph order"
   nothing
end

"""
    QWSearch([type, ]ctqw, marked[, penalty, jumpingrate])

Creates `QWSearch` according to `AbstractCTQW` model. The evolution is defined by
Hamiltonian, is the sum of the quantum walk part and the oracle. By default `type` equals
`ComplexF64`, `jumpingrate` equals largest eigenvalue of adjacency matrix of graph if
`matrix(CTQW)` outputs `:adjacency` and error otherwise, and `penalty` equals 0.
Details can be found in the original [paper](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.70.022314).

    QWSearch(qws_ctqw; marked, penalty)

Updates quantum walk search to new subset of marked elements and new penalty. By
default marked and penalty are the same as in qws.

"""
function QWSearch(::Type{T},
                  ctqw::AbstractCTQW,
                  marked::Vector{Int},
                  penalty::Real = 0.,
                  jumpingrate::T = jumping_rate(T, ctqw)) where T<:Number

   hamiltonian = jumpingrate*graph_hamiltonian(T, ctqw)
   hamiltonian -= sum(proj(T, v, nv(ctqw.graph)) for v=marked)

   parameters = Dict{Symbol,Any}(:hamiltonian => hamiltonian)

   QWSearch(ctqw, parameters, marked, penalty)
end

function QWSearch(::Type{T},
                  ctqw::CTQWDense,
                  marked::Vector{Int},
                  penalty::Real = 0.,
                  jumpingrate::T = jumping_rate(T, ctqw)) where T<:Number

   hamiltonian = jumpingrate*Matrix(graph_hamiltonian(T, ctqw))
   hamiltonian -= sum(proj(T, v, nv(ctqw.graph)) for v=marked)

   parameters = Dict{Symbol,Any}(:hamiltonian => hamiltonian)

   QWSearch(ctqw, parameters, marked, penalty)
end

function QWSearch(ctqw::AbstractCTQW,
                  marked::Vector{Int},
                  penalty::Real = 0.,
                  jumpingrate::Real = jumping_rate(Float64, ctqw))
   QWSearch(ComplexF64, ctqw, marked, penalty, ComplexF64(jumpingrate))
end

function QWSearch(qws::QWSearch{<:AbstractCTQW};
                  marked::Vector{Int}=qws.marked,
                  penalty::Real=qws.penalty)
   oldmarked = qws.marked

   hamiltonian = copy(parameters(qws)[:hamiltonian])
   hamiltonian -= sum(proj(eltype(hamiltonian), v, nv(graph(qws))) for v=marked)
   hamiltonian += sum(proj(eltype(hamiltonian), v, nv(graph(qws))) for v=oldmarked)

   QWSearch(model(qws), Dict(:hamiltonian => hamiltonian), marked, penalty)
end

"""
    check_qwdynamics(QWSearch, ctqw, parameters, marked)

Checks whetver combination of `ctqw`, `marked` and `parameters` produces valid
`QWSearch` object. It checks if `parameters` consists of key `:hamiltonian` with
real-valued matrix. Furthermore the hamiltonian
needs to be square of size equals to `graph(ctqw)` order. the hermiticity is
not checked for efficiency issue.
"""
function check_qwdynamics(::Type{QWSearch},
                          ctqw::AbstractCTQW,
                          parameters::Dict{Symbol},
                          marked::Vector{Int})
   check_qwmodel(ctqw, parameters)
end

"""
    QWEvolution([type_ctqw, ]ctqw)

Creates `QWEvolution` according to `AbstractCTQW` model. By default `type` equals
`ComplexF64`. The Hamiltonian of the system is either (normalized) Laplacian, or
minus adjacency matrix.
"""
function QWEvolution(::Type{U},
                     ctqw::AbstractCTQW) where U<:Number
   parameters = Dict{Symbol,Any}(:hamiltonian => graph_hamiltonian(U, ctqw))
   QWEvolution(ctqw, parameters)
end

function QWEvolution(ctqw::AbstractCTQW)
   QWEvolution(ComplexF64, ctqw)
end

"""
    check_qwdynamics(QWEvolution, ctqw, parameters)

Checks if combination of `ctqw` and `parameters` produces valid
`QWSearch` object. It checks if `parameters` consists of key `:hamiltonian` with
real-valued matrix. Furthermore the hamiltonian
needs to be square of size equals to `graph(ctqw)` order. The hermiticity is
not checked for efficiency issues.
"""
function check_qwdynamics(::Type{QWEvolution},
                          ctqw::AbstractCTQW,
                          parameters::Dict{Symbol})
   check_qwmodel(ctqw, parameters)
end

include("ctqw_utils.jl")
include("ctqw_evolution.jl")


text:
using ModelingToolkit, Test
@parameters t a
@variables x y
D = Differential(t)
eqs = [D(x) ~ a * x - x * y,
    D(y) ~ -3y + x * y]
f = build_function([x.rhs for x in eqs], [x, y], [a], t, expression = Val{false},
                   target = ModelingToolkit.CTarget())
f2 = eval(build_function([x.rhs for x in eqs], [x, y], [a], t)[2])
du = rand(2);
du2 = rand(2);
u = rand(2)
p = rand(1)
_t = rand()
f(du, u, p, _t)
f2(du2, u, p, _t)
@test du == du2


text:
module solver

function euler(y0::AbstractArray{Float64}, f, a::Float64, b::Float64, h::Float64)::AbstractArray{Float64}
    k = ceil(Integer, (b - a) / h)
    t_range = range(a, stop=b, length=k)
    y = y0
    for t_n = t_range
        if t_n == a
            continue
        end
        y_n = y[:,end]
        y_next = y_n + h * f(t_n, y_n)
        y = [y y_next]
    end
    return y
end

function midpoint(y0::AbstractArray{Float64}, f, a::Float64, b::Float64, h::Float64)::AbstractArray{Float64}
    k = ceil(Integer, (b - a) / h)
    t_range = range(a, stop=b, length=k)
    y = y0
    for t_n = t_range
        if t_n == a
            continue
        end
        y_n = y[:,end]
        y_next = y_n + h * f(t_n + h / 2.0, y_n + h / 2.0 * f(t_n, y_n))
        y = [y y_next]
    end
    return y
end

end # module

max_stars_repo_name: applied-math-coding/article-snippets

max_issues_repo_name: applied-math-coding/article-snippets

max_forks_repo_name: applied-math-coding/article-snippets


text:
function data2SOM!(
    data   :: AbstractArray,
    Nz_bot :: AbstractArray{Int64, 3},
)

    _, Nx, Ny = size(Nz_bot)

    for i=1:Nx, j=1:Ny

        _Nz_bot = Nz_bot[1, i, j]

        if _Nz_bot != 0
            data[1:_Nz_bot, i, j] .= data[1, i, j]
        end

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
#---------------------------------------------------------------------#
#This function computes the global Mass and Differentiation Matrices.
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function global_matrices(ψ,dψdξ,ωq,intma,coord,Npoin,Ne,Np,Nq,DFloat)

    #Initialize
    M=zeros(DFloat,Npoin,Npoin)
    L=zeros(DFloat,Npoin,Npoin)
    inode=zeros(Int64,Np)
    x=zeros(DFloat,Np)

    #Loop through elements
    for e=1:Ne
        #Store Coordinates
        for i=1:Np
            I=intma[i,e]
            inode[i]=I
            x[i]=coord[I]
        end

        #Element Length & Jacobian
        dx=x[Np]-x[1]
        dxdξ=dx/2
        dξdx=2.0/dx

        #Build Matrices
        @inbounds for k=1:Nq
            for j=1:Np
                J=inode[j]
                for i=1:Np
                    I=inode[i]
                    M[I,J]+=ωq[k]*ψ[i,k]*ψ[j,k]*dxdξ
                    L[I,J]+=ωq[k]*(dψdξ[i,k]*dξdx)*(dψdξ[j,k]*dξdx)*dxdξ
                 end #i
            end #j
        end #k
    end #e

    #Impose Dirichlet BC on L
    for J=1:Npoin
        L[1,J]=0
        L[Npoin,J]=0
    end
    L[1,1]=DFloat(1.0)
    L[Npoin,Npoin]=DFloat(1.0)

    #Return arrays
    return(M,L)
end


text:
module InterpolatingFunctions

abstract AbstractInterpolation

include("grid.jl")
include("boundary.jl")
include("linear.jl")
include("bilinear.jl")
include("cubic.jl")
include("bicubic.jl")
include("polyharmonic.jl")
include("polynomial.jl")

typealias AI AbstractInterpolation
typealias AG AbstractGrid
typealias ABC AbstractBoundaryCondition

function interpolate_gen(I,G,BC)
    bc_expr = apply_bc(I, BC, :interp, :x)
    ifunc = quote
        interp = preprocess($I, axes, f)
        function InterpolatingFunction(args...)
            if length(args) == 1
                if typeof(args[1]) <: AbstractArray
                    x = args[1]
                else
                    x = [args[1]]
                end
            else
                x = collect(args)
            end
            $(bc_expr)

            inds = get_indices($I, $G, interp, x)
            val = get_value($I, interp, x, inds)
            return val
        end
    end
    return ifunc
end

function interpolate{N,T,I<:AI}(f::AbstractArray{T}, axes::NTuple{N,AbstractArray}, ::Type{I})
    G,BC = default_settings(I,axes)
    interpolate(f,axes,I,G,BC)
end

@generated function interpolate{N,T,I<:AI,G<:AG,BC<:ABC}(f::AbstractArray{T}, axes::NTuple{N,AbstractVector},::Type{I},::Type{G},::Type{BC})
    interpolate_gen(I,G,BC)
end

export interpolate, InterpolatingFunction
export Uniform, Irregular
export Nearest, Value, Reflect, Error, None
export CubicSpline, Monotonic, NonMonotonic
export PolyharmonicSpline, polyharmonic
export Polynomial, LinearSpline
export BicubicSpline, BilinearSpline

end # module


text:
function nn_lstm(::Type{T}, vocab::Int) where T
    wordembeds = Normal(0,0.01)(T, 100, vocab)
    x = lookup(parameter(wordembeds), Node(name=:x))
    hsize = size(wordembeds, 1)
    h = LSTM(T,hsize,hsize,1,0.0,false)(x, Node(name=:dims_x))
    h = Linear(T,hsize,size(wordembeds,2))(h)
    Graph(h)
end

mutable struct LM_LSTM <: Parametric
    params
end

Merlin.parameters(f::LM_LSTM) = parameters(f.params...)
Merlin.todevice(f::LM_LSTM, device) = LM_LSTM(todevice(f.params,device))

function LM_LSTM(::Type{T}, vocab::Int) where T
    wordembeds = Normal(0,0.01)(T, 100, vocab)
    hsize = size(wordembeds, 1)
    lstm = LSTM(T, hsize, hsize, 1, 0. 0, false)
    linear = Linear(T, hsize, size(wordembeds,2))
    params = (wordembeds=wordembeds, lstm=lstm, linear=linear)
    LM_LSTM(params)
end

function (nn::LM_LSTM)(nt::NamedTuple)
    p = nn.params
    h = lookup(p.wordembeds, nt.x)
    h = p.lstm(h, nt.dims_x)
    z = p.linear(h)
    if t.train
        l = softmax_crossentropy(nt.y, z)
        gradient!(l)
        sum(Array(l.data))
    else
        throw("Not implemented yet.")
    end
end


text:
# simulation controls
alphas = PROWIMData.props_on["lift_distribution_alpha"]
ploti = 1:length(alphas)
# wing definition
wing_b = PROWIMData.wing["span"]
wing_TR = PROWIMData.wing["TR"] # 1.0 # hershey bar
wing_c = PROWIMData.wing["chord"] # 240e-3
wing_AR = wing_b / wing_c
wing_θroot = PROWIMData.wing["theta_root"]
wing_θtip = PROWIMData.wing["theta_tip"]
wing_le_sweep = PROWIMData.wing["le_sweep"]
wing_ϕ = PROWIMData.wing["dihedral"]

wing_lexs = [[0.0, 0.0]]
wing_leys = [[0.0, wing_b/2]]
wing_lezs = [[0.0, 0.0]]
wing_chords = [[wing_c, wing_c]]
wing_thetas = [[wing_θroot, wing_θtip]]
wing_phis = [[wing_ϕ, wing_ϕ]]
wing_nspanwisepanels = 50
wing_nchordwisepanels = 1

# environment
environment = AS.Environment()
ν = environment.ν
# freestream definition
Re_c = PROWIMData.props_on["Re_c"]
vinf = Re_c * ν / wing_c
vinfs = fill(vinf, length(alphas))
# rotor definition

# Epema propeller w/ 4 blades to match
nblades = [PROWIMData.rotor["blades"]]
rtip = [PROWIMData.rotor["diameter"]/2]
radii = fill(PROWIMData.rotor["radii"], length(nblades))
rhub = fill(radii[1][1], length(nblades))
J = PROWIMData.rotor["J"]
D = PROWIMData.rotor["diameter"]
n = vinf / J / D # rev/s
omega = n .* 2 * pi # rad/s
omegas = ones(length(nblades),length(alphas)) .* omega # rad/s
chords_raw = PROWIMData.rotor["chords_raw"]
chords = fill(AS.FM.linear(chords_raw[:,1], chords_raw[:,2], radii[1]), length(nblades))

twists_raw = PROWIMData.rotor["twists_raw"]
twist = AS.FM.linear(twists_raw[:,1] .* rtip[1], twists_raw[:,2], radii[1]) .* pi/180 # radians
# set collective to 30 deg. at 70\% radius
twists_07_desired = [PROWIMData.rotor["twist_07"]]
twists_07 = AS.FM.linear(radii[1], twist, 0.7 * rtip[1])
twists = Vector{Vector{typeof(twist[1])}}(undef, length(nblades))
twist_extra = -0.3 * pi/180
for (i, twistval) in enumerate(twists_07_desired)
    twists[i] = twist .+ twistval .- twists_07 .+ twist_extra
end

locations = string.(Int.(round.(rs_desired .* 1000, digits=0)))
locations[1:end-1] = "0" .* locations[1:end-1]
# uncomment below to use distinct airfoils at each section
# contourfilenames = "epema_interpolated_bspline_n30_" .* locations
# uncomment below to repeat the r/R = 0.7 airfoil across the blade
contour_filenames = fill("epema_interpolated_bspline_n30_" * locations[8],length(radii[1]))
contour_paths = fill(joinpath.(contour_directory, contour_filenames .* ".dat"),length(nblades))

polar_directory = polar_directory
plotstepi = 1:length(alphas)

Res = [5e4, 1e5, 5e5, 1e6, 1e7]
Machs = [0.0, 0.1, 0.2, 0.3]

Res_list = fill(fill(Res, length(radii[1])), length(nblades))
Ms_list = fill(fill(Machs, length(radii[1])), length(nblades))

index = [1]
rotor_positions = [[-PROWIMData.rotor["x"], PROWIMData.rotor["y"], 0.0]]
rotor_orientations = [[-1.0, 0.0, 0.0]]
spin_directions = [true]

# call template
plot_base_name = "PROWIM_props_on"
wake_developement_factors = [1.0] #[0.0, 0.5, 1.0]
swirl_recovery_factors = [0.5] #[0.0, 0.5, 1.0]
wake_developement_factor = 1.0 #[0.0, 0.5, 1.0]
swirl_recovery_factor = 0.5 #[0.0, 0.5, 1.0]

data_PROWIM_vlm_bem = AS.vlm_bem_template(vinfs, plotstepi, alphas,
    wing_lexs, wing_leys, wing_lezs, wing_chords, wing_thetas, wing_phis, wing_nspanwisepanels, wing_nchordwisepanels,
    omegas, nblades, rhub, rtip, radii, chords, twists,
    contour_paths, index,
    rotor_positions, rotor_orientations, spin_directions, Res_list, Ms_list;
        Vref = 1.0, symmetric = true, iref = 1, static_margin = 0.10, liftingline_x_over_c = 0.25,
        cambers = [fill((xc) -> 0, length(lex)) for lex in wing_lexs],
        spanwise_spacings = fill(VL.Uniform(), length(wing_lexs)),
        chordwise_spacings = fill(VL.Uniform(), length(wing_lexs)),
        polar_directory,
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames = ["wing"],
        rotor_names = ["rotor"],
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(AS.topdirectory, "data","plots",TODAY),
        plot_base_name,
        plot_extension = ".png",
        step_symbol = LS.L"\alpha "
)
objective = AS.runsimulation!(data_PROWIM_vlm_bem...)

# plot validation
waketag = "_wdf_$(round(wake_developement_factor; digits=1))_srf_$(round(swirl_recovery_factor; digits=1))"
fig_rotor = plt.figure(plot_base_name * "_rotor_sweep")
ax_rotor = fig_rotor.get_axes()[1]
ax_rotor.scatter(0.85, 0.168, c="r", label="Velduis")
ax_rotor.legend(loc="upper left", bbox_to_anchor=(1.01,1))
fig_rotor.tight_layout()
fig_rotor.savefig(joinpath(plot_directory,"PROWIM_props_on_rotor_sweep" * waketag * ".pdf"), bbox_inches="tight")
# fig_rotor.savefig(joinpath(notebookdirectory,"PROWIM_propson_rotorsweep.pdf"), bbox_inches="tight")

fig_cf = plt.figure(plot_base_name * "_cf_distribution")
# local fig_lift = plt.figure(plot_base_name * "_lift_distribution")
fig_lift = plt.figure(plot_base_name * "_lift_distribution")
axs_cf = fig_cf.get_axes()
# local ax_lift = fig_lift.get_axes()[1]
ax_lift = fig_lift.get_axes()[1]
cl_data_props_on = PROWIMData.props_on["lift_distribution"]
plot_labels = "Velduis, " .* [LS.L"\alpha = 0^\circ", LS.L"\alpha = 4^\circ", LS.L"\alpha = 10^\circ"]
# local aircraft = data_PROWIM_vlm_bem[1]
aircraft = data_PROWIM_vlm_bem[1]
b = aircraft.wing_system.lifting_line_rs[1][2,end]
for (i, data) in enumerate(cl_data_props_on)
    cratio = i / length(cl_data_props_on)
    axs_cf[3].scatter(data[:,1] .* b, data[:,2], marker="+", color=(0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label=plot_labels[i])
    ax_lift.scatter(data[:,1], data[:,2], marker="+", color=(0.05, 0.85-cratio*0.7, 0.15 + 0.75 * cratio), label=plot_labels[i])
end
axs_cf[3].legend(loc="upper left", bbox_to_anchor=(1.01,1))
fig_cf.tight_layout()
fig_cf.savefig(joinpath(plot_directory,"PROWIM_props_on_cf_distribution" * waketag * ".pdf"), bbox_inches="tight")
ax_lift.legend(loc="upper left", bbox_to_anchor=(1.01,1))
fig_lift.tight_layout()
fig_lift.savefig(joinpath(plot_directory,"PROWIM_props_on_lift_distribution" * waketag * ".pdf"), bbox_inches="tight")
# fig_cf.savefig(joinpath(notebookdirectory,"PROWIM_propson_lift_distribution" * waketag * ".pdf"), bbox_inches="tight")

fig_CF = plt.figure(plot_base_name * "_cl_alpha_sweep")
axs_CF = fig_CF.get_axes()

# Fetch data.
# local CL_data = PROWIMData.props_on["CL_alpha"]
CL_data = PROWIMData.props_on["CL_alpha"]
# local CD_data_balance = PROWIMData.props_on["CD_balance"]
CD_data_balance = PROWIMData.props_on["CD_balance"]
# local CD_data_wakesurvey = PROWIMData.props_on["CD_wake_survey"]
CD_data_wakesurvey = PROWIMData.props_on["CD_wake_survey"]

axs_CF[3].scatter(CL_data[:,1], CL_data[:,2], marker = "+", label="Velduis")
axs_CF[3].legend(loc="upper left", bbox_to_anchor=(1.01,1))
axs_CF[1].scatter(CD_data_balance[:,1], CD_data_balance[:,2], marker = "+", label="Velduis, balance")
axs_CF[1].scatter(CD_data_wakesurvey[:,1], CD_data_wakesurvey[:,2], marker = "+", label="Velduis, wake survey")
axs_CF[1].legend(loc="upper left", bbox_to_anchor=(1.01,1))
fig_CF.tight_layout()
fig_CF.savefig(joinpath(plot_directory,"PROWIM_props_on_cl_alpha_sweep" * waketag * ".pdf"), bbox_inches="tight")
# fig_CF.savefig(joinpath(notebookdirectory,"PROWIM_propson_clalphasweep" * waketag * ".pdf"), bbox_inches="tight")

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Test

@time @testset "Potentials.jl" begin
    include("lj_test.jl")
end


text:
using PkgBenchmark
using Dates

benchmarkpkg(
    dirname(@__DIR__),
    BenchmarkConfig(
        env = Dict(
            "JULIA_NUM_THREADS" => "1",
            "OMP_NUM_THREADS" => "1",
        ),
    ),
    resultfile = joinpath(@__DIR__, "results", "$(Dates.format(now(), dateformat"yyyymmddTHHMMSS")).json"),
)


text:
# This file was generated, do not modify it. # hide
df.waterfront = (df.waterfront .!= "FALSE")
coerce!(df, :waterfront => OrderedFactor);


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# test/MPBWrapper.jl
# Test the MathProgBase.jl interface for the ECOS.jl solver wrapper
#############################################################################

@testset "Test the MPB wrapper with linprog" begin
    include("MPB_linear.jl")
end

import ECOS

import MathProgBase
@static if VERSION >= v"0.7-"
    const MPB_test_path = joinpath(dirname(pathof(MathProgBase)), "..", "test")
else
    const MPB_test_path = joinpath(Pkg.dir("MathProgBase"), "test")
end

@testset "Run the conic interface test from MathProgBase.jl" begin
    include(joinpath(MPB_test_path, "conicinterface.jl"))
    coniclineartest(ECOS.ECOSSolver(), duals=true)
    conicSOCtest(ECOS.ECOSSolver(), duals=true)
    conicEXPtest(ECOS.ECOSSolver(), duals=true)
    include(joinpath(MPB_test_path, "quadprog.jl"))
    socptest(ECOS.ECOSSolver())
end


text:
using TuringModels

Turing.setadbackend(:reversediff);

d = CSV.read(joinpath(@__DIR__, "..", "..", "data", "reedfrogs.csv"), delim=';');

size(d) # Should be 48x5

# Set number of tanks

d[!, :tank] = 1:size(d,1);

# Define the Turing model

@model m12_1(density, tank, surv) = begin

    # Number of unique tanks in the data set
    N_tank = length(tank)
    a_tank ~ filldist(Normal(0, 5), N_tank)

    logitp = a_tank[tank]
    surv .~ BinomialLogit.(density, logitp)
end

# Sample

chns = sample(m12_1(Vector{Int64}(d[:, :density]), Vector{Int64}(d[:, :tank]),
    Vector{Int64}(d[:, :surv])), Turing.NUTS(0.65), 1000);

# CmdStan results

m2_1_rethinking = "
             mean   sd  5.5% 94.5% n_eff Rhat
a_tank[1]   2.49 1.16  0.85  4.53  1079    1
a_tank[2]   5.69 2.75  2.22 10.89  1055    1
a_tank[3]   0.89 0.75 -0.23  2.16  1891    1
a_tank[4]   5.71 2.70  2.21 10.85   684    1
a_tank[5]   2.52 1.14  0.92  4.42  1640    1
a_tank[6]   2.49 1.13  0.94  4.52  1164    1
a_tank[7]   5.74 2.71  2.25 10.86   777    1
a_tank[8]   2.52 1.19  0.95  4.42  1000    1
a_tank[9]  -0.46 0.69 -1.62  0.55  2673    1
a_tank[10]  2.53 1.19  0.93  4.59  1430    1
a_tank[11]  0.93 0.72 -0.17  2.11  1387    1
a_tank[12]  0.47 0.74 -0.63  1.70  1346    1
a_tank[13]  0.91 0.76 -0.25  2.30  1559    1
a_tank[14]  0.00 0.66 -1.04  1.06  2085    1
a_tank[15]  2.50 1.19  0.95  4.40  1317    1
a_tank[16]  2.50 1.14  0.98  4.31  1412    1
a_tank[17]  3.49 1.12  1.94  5.49   945    1
a_tank[18]  2.59 0.75  1.50  3.81  1561    1
a_tank[19]  2.11 0.64  1.15  3.15  1712    1
a_tank[20]  6.40 2.57  3.11 11.04   996    1
a_tank[21]  2.59 0.74  1.54  3.93  1233    1
a_tank[22]  2.63 0.79  1.49  4.01  1184    1
a_tank[23]  2.64 0.83  1.45  4.13  1379    1
a_tank[24]  1.74 0.59  0.85  2.72  1736    1
a_tank[25] -1.19 0.45 -1.90 -0.50  2145    1
a_tank[26]  0.09 0.41 -0.53  0.78  2167    1
a_tank[27] -1.75 0.56 -2.65 -0.88  1666    1
a_tank[28] -0.58 0.43 -1.25  0.08  1567    1
a_tank[29]  0.08 0.39 -0.54  0.71  3053    1
a_tank[30]  1.43 0.49  0.66  2.24  2754    1
a_tank[31] -0.79 0.44 -1.50 -0.12  1299    1
a_tank[32] -0.42 0.41 -1.12  0.23  1661    1
a_tank[33]  3.84 1.08  2.31  5.70   808    1
a_tank[34]  3.00 0.85  1.83  4.36  1038    1
a_tank[35]  2.96 0.82  1.82  4.25  1578    1
a_tank[36]  2.14 0.55  1.31  3.08  1734    1
a_tank[37]  2.12 0.56  1.31  3.04  1131    1
a_tank[38]  6.72 2.62  3.45 11.44   706    1
a_tank[39]  2.95 0.73  1.85  4.08  1509    1
a_tank[40]  2.48 0.65  1.53  3.61  1731    1
a_tank[41] -2.15 0.57 -3.11 -1.29  1231    1
a_tank[42] -0.67 0.35 -1.22 -0.14  1444    1
a_tank[43] -0.54 0.35 -1.12  0.03  1776    1
a_tank[44] -0.43 0.34 -1.00  0.10  1735    1
a_tank[45]  0.54 0.36 -0.04  1.14  1376    1
a_tank[46] -0.67 0.34 -1.25 -0.15  1619    1
a_tank[47]  2.14 0.55  1.31  3.04  1916    1
a_tank[48] -0.06 0.35 -0.61  0.50  1932    1
";

# Describe chainsd

describe(chns)

# End of `12/m12.1t.jl`


text:
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

@test sprint(MOIU.print_with_acronym, "MathOptInterface") == "MOI"
@test sprint(MOIU.print_with_acronym, "MathOptInterface.MathOptInterface") ==
      "MOI.MOI"
@test sprint(
    MOIU.print_with_acronym,
    "MathOptInterface.Utilities.MathOptInterface",
) == "MOIU.MOI"
@test sprint(MOIU.print_with_acronym, "MathOptInterfaceXXBridges") ==
      "MOIXXBridges"
@test sprint(MOIU.print_with_acronym, "MathOptInterface.BridgesXX") == "MOIBXX"
@test sprint(MOIU.print_with_acronym, "MathOptInterface.Test.x") == "MOIT.x"
@test sprint(MOIU.print_with_acronym, "MathOptInterface.x.Test") == "MOI.x.Test"
@test sprint(MOIU.print_with_acronym, "MathOptInterface.Utilities.Test") ==
      "MOIU.Test"
@test sprint(MOIU.print_with_acronym, "MathOptInterface.Utilities.Test") ==
      "MOIU.Test"


text:
@testset "OneHot primitive" begin
  using PrimitiveOneHot: OneHot
  idx = rand(5:10)
  ot = OneHot(10)
  o = OneHot{10}(idx)

  @test PrimitiveOneHot.onehotsize(o) == 10
  @test o === OneHot(10, idx)
  @test o === ot(idx)
  @test o == OneHot(11, idx)
  @test o !== OneHot(11, idx)

  @test size(o) == (10,)
  for i = 1:10
    @test o[i] == (i == idx)
  end
  @test o[:] === o
  @test_throws BoundsError o[11]

  @test UInt32(o) == UInt32(idx)
  @test UInt64(o) == UInt64(idx)
  @test Int32(o) == Int32(idx)
  @test Int64(o) == Int64(idx)

  @test one(o) < o

  @test convert(UInt32, o) == UInt32(idx)
  @test convert(UInt64, o) == UInt64(idx)
  @test convert(Int32, o) == Int32(idx)
  @test convert(Int64, o) == Int64(idx)
  @test convert(Int8, o) == Int8(idx)
  @test_throws InexactError convert(Bool, o)
  @test convert(ot, idx) == o
  @test convert(ot, UInt8(idx)) == o
  @test convert(ot, UInt16(idx)) == o
  @test convert(ot, UInt32(idx)) == o
  @test convert(ot, UInt64(idx)) == o
  @test convert(ot, UInt128(idx)) == o
  @test convert(ot, Int8(idx)) == o
  @test convert(ot, Int16(idx)) == o
  @test convert(ot, Int32(idx)) == o
  @test convert(ot, Int64(idx)) == o
  @test convert(ot, Int128(idx)) == o

  @test convert(ot, true) == OneHot(10, 1)
  @test convert(ot, false) == OneHot(10, 0)

  @test convert(OneHot(15), o) == OneHot(15, idx)
  @test convert(OneHot(idx+1), o) == OneHot(idx+1, idx)
  @test_throws PrimitiveOneHot.OneHotEncodeError convert(OneHot(idx-1), o)

  @test convert(Array, o) == begin
    z = zeros(Bool, 10)
    z[idx] = true
    z
  end
  @test reinterpret(Int32, o) == Int32(idx)

  @test one(o) === OneHot(10, 1)
  @test one(o) === one(typeof(o))
  @test isone(one(o))

  @test zero(o) === OneHot(10, 0)
  @test zero(o) === zero(typeof(o))
  @test !iszero(o)
  @test iszero(zero(o))

  @test_throws PrimitiveOneHot.OneHotEncodeError convert(ot, -1)
  @test_throws PrimitiveOneHot.OneHotEncodeError convert(ot, 12)
  @test_throws PrimitiveOneHot.OneHotEncodeError convert(ot, typemax(Int64))
  @test_throws PrimitiveOneHot.OneHotEncodeError convert(ot, typemax(UInt64))

  @test iszero(typemin(typeof(o)))
  @test Int(typemax(typeof(o))) == 10

  @test argmax(o) == argmax(collect(o))
  @test argmin(o) == argmin(collect(o))
  @test argmin(one(o)) == 2

end


text:
using COSMO, Test, Random

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "Printing" begin
for T in UnitTestFloats

   ws = COSMO.Workspace{T}()
   iter = 10
   cost = 20.
   r_prim = 1.5e-3
   r_dual = 1.2e-2
   status = :Solved
   rt = 0.7

   @testset "Printing (T = $(T))" begin
      @test COSMO.print_iteration(ws, iter, cost, r_prim, r_dual) == nothing
      @test COSMO.print_iteration(ws, ws.settings.check_termination, cost, r_prim, r_dual) == nothing
      @test COSMO.print_result(status, iter, cost, rt) == nothing
   end
end
end
nothing