text:
"""
    Difference array

# Brief:
    For any array A with N numbers, we can create difference array D with N numbers also.
    Every k-th element of D is equal to difference between A[k] and A[k-1]
        D[k] = A[k-1] - A[k]

# Complexity of some operations
    - add some value x to the m consecutive elements of array - O(1)
    - print array after any numbers of changes                - O(N)

# Functions
    - create_diff_arr(original::Array{T}) - Create difference array for array 'original'
    - calculate_arr(diff_arr::Array{T}) - Create a original array from the given difference array
    - add_to_arr(diff_arr::Array{T}, l::Int, r::Int, x::Number) - Add x to all elements with index from [l, r]

# Contributed by: [Nikola Mircic](https://github.com/Nikola-Mircic)
"""
module DifferenceArray
# Create difference array for array 'original'
# Parameters:
#   - original - an array which is used for calculating the difference array
function create_diff_arr(original::Array{T}) where {T<:Number}
    n = length(original)

    diff_arr = Array(original)

    for i in 2:n
        diff_arr[i] = original[i] - original[i-1]
    end

    return diff_arr
end

# Create a original array from the given difference array
# Parameters:
#   - diff_arr - an difference array which is used for calculating the original array
function calculate_arr(diff_arr::Array{T}) where {T<:Number}
    n = length(diff_arr)

    arr = Array(diff_arr)

    for i in 2:n
        arr[i] = diff_arr[i] + arr[i-1]
    end

    return arr
end

# Add x to all elements with index from [l, r]
# Parameters:
#   - dif_arr - a difference array of the array you want to change
#   - l - leftmost index of the affected range
#   - r - rightmost index of the affected range
#   - x - a value to be added to all elements from a given range
function add_to_arr(
    diff_arr::Array{T},
    l::Int,
    r::Int,
    x::Number,
) where {T<:Number}
    diff_arr[l] += x
    if r < length(diff_arr)
        diff_arr[r+1] -= x
    end
end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using CosmoParticles
using Documenter

DocMeta.setdocmeta!(CosmoParticles, :DocTestSetup, :(using CosmoParticles); recursive=true)

makedocs(;
    modules=[CosmoParticles],
    authors="Lucas Valenzuela",
    repo="https://github.com/lucasvalenzuela/CosmoParticles.jl/blob/{commit}{path}#{line}",
    sitename="CosmoParticles.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lucasvalenzuela.github.io/CosmoParticles.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Particles" => "particles.md",
        "Particle Collections" => "particlecollection.md",
        "Operations" => "operations.md",
        "Properties" => "properties.md",
        "Geometry" => "geometry.md",
        "Utils" => "utils.md",
    ],
)

deploydocs(;
    repo="github.com/lucasvalenzuela/CosmoParticles.jl",
    devbranch="main",
)

ext: jl


text:
import .CUDA
CUDA.libcufft()
CUDA.libcudnn()
CUDA.libcurand()
path_cublas = CUDA.libcublas()
using Libdl
if isdefined(CUDA, :libcudart)
    CUDA.libcudart()
else
    @warn """
    HACK
    CUDA.jl version not expose `libcudart`. See
    https://discourse.julialang.org/t/shared-c-dependencies-and-artifacts/68525/4?u=jw3126
    But as an implementation detail libcudart does ship with the cublas artifact
    """
    dir = splitdir(path_cublas)[1]
    path_cudart = joinpath(dir, "libcudart.so.11.0")
    dlopen(path_cudart)
end


text:
# Thanks @tlienart

using Franklin, PkgPage

using Markdown

function hfun_doc(params)
    fname = join(params[1:max(1, length(params)-2)], " ")
    head = params[end-1]
    type = params[end]
    doc = eval(Meta.parse("using FileTrees; @doc FileTrees.$fname"))
    txt = Markdown.plain(doc)
    # possibly further processing here
    body = Franklin.fd2html(txt, internal=true)
    return """

            $head

$type

$body

    """
end


text:
using CloudGraphs
using Caesar, IncrementalInference, RoME
using KernelDensityEstimate

include("../entities/roverPose.jl")

configuration = CloudGraphs.CloudGraphConfiguration("localhost", 7474, "neo4j", "neo5j", "localhost", 27017, false, "", "");
cloudGraph = connect(configuration);
conn = cloudGraph.neo4j.connection;
session = "RoverLock_" * string(Base.Random.uuid1())[1:8] #Name+SHA
# register types of interest in CloudGraphs
registerGeneralVariableTypes!(cloudGraph)
Caesar.usecloudgraphsdatalayer!()
println("Current session: $session")

fg = Caesar.initfg(sessionname=session, cloudgraph=cloudGraph)

# REF: https://github.com/dehann/Caesar.jl/blob/master/examples/wheeled/victoriapark_onserver.jl
# Tuning params - Move these out.
Podo=diagm([0.1;0.1;0.005]) # TODO ASK: Noise?
N=100
lcmode=:unimodal # TODO ASK: Solver?
lsrNoise=diagm([0.1;1.0]) # TODO ASK: ?

lastPoseVertex = initFactorGraph!(fg)
curPose = RoverPose();
curPose.timestamp = Base.Dates.datetime2unix(now())
# Bump to x2 so we don't have two x1's (x1 = base node)
curPose.poseIndex = 2
for i = 1:2
    curPose.x += 1
    curPose.y += 1
    @time lastPoseVertex, factorPose = addOdoFG!(fg, poseIndex(curPose), odoDiff(curPose), Podo, N=N, labels=["POSE"])
    curPose = RoverPose(curPose)
end


text:
function test_inner_product_layer(backend::Backend, n_input, T, eps)
  println("-- Testing InnerProductLayer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################
  # inner-product layer should be able to take blobs with
  # different num, but the input feature dimensions should
  # be the same
  batch_size   = [abs(rand(Int)) % 50 + 1 for i = 1:n_input]
  orig_dim_all = (10, 20, 3)
  orig_dim     = prod(orig_dim_all)
  target_dim   = 30

  X = [rand(T, orig_dim_all..., n) for n in batch_size]
  W = rand(T, orig_dim, target_dim)
  b = rand(T, target_dim)

  ############################################################
  # Setup
  ############################################################
  layer = InnerProductLayer(name="ip", output_dim=target_dim,
      tops=Array{Symbol}(n_input), bottoms=Array{Symbol}(n_input))
  inputs = Blob[make_blob(backend, x) for x in X]
  diffs = Blob[make_blob(backend, x) for x in X]

  println("    > Setup")
  state  = setup(backend, layer, inputs, diffs)

  @test length(state.W) == length(W)
  @test length(state.b) == length(b)
  copy!(state.W, W)
  copy!(state.b, b)

  println("    > Forward")
  forward(backend, state, inputs)

  X2 = [reshape(X[i], orig_dim, batch_size[i]) for i = 1:n_input]
  res = [W' * x2 .+ b for x2 in X2]

  res_layer = [similar(x) for x in res]
  for i = 1:n_input
    copy!(res_layer[i], state.blobs[i])
    @test all(-eps .< res_layer[i] - res[i] .< eps)
  end

  println("    > Backward")
  top_diff = [rand(T, size(state.blobs_diff[i])) for i = 1:n_input]
  for i = 1:n_input
    copy!(state.blobs_diff[i], top_diff[i])
  end
  backward(backend, state, inputs, diffs)

  top_diff = [reshape(top_diff[i], target_dim, batch_size[i]) for i = 1:n_input]
  bias_grad = similar(b)
  copy!(bias_grad, state.∇b)
  bias_grad_expected = sum([sum(top_diff[i],2) for i = 1:n_input])
  @test all(-eps .< vec(bias_grad) - vec(bias_grad_expected) .< eps)

  X_mat = [reshape(X[i], orig_dim, batch_size[i]) for i = 1:n_input]
  weight_grad = similar(W)
  copy!(weight_grad, state.∇W)
  weight_grad_expected = sum([X_mat[i]*top_diff[i]' for i = 1:n_input])
  @test all(-eps .< vec(weight_grad) - vec(weight_grad_expected) .< eps)

  back_grad = [similar(X_mat[i]) for i = 1:n_input]
  for i = 1:n_input
    copy!(back_grad[i], diffs[i])
    @test all(-eps .< vec(back_grad[i]) - vec(W * top_diff[i]) .< eps)
  end

  shutdown(backend, state)
end
function test_inner_product_layer(backend::Backend)
  test_inner_product_layer(backend, 3, Float32, 1e-3)
  test_inner_product_layer(backend, 3, Float64, 1e-10)
end

if test_cpu
  test_inner_product_layer(backend_cpu)
end
if test_gpu
  test_inner_product_layer(backend_gpu)
end


text:
"""
    QuadraticOTNewton

Semi-smooth Newton method (Algorithm 2 of Lorenz et al. 2019 [^LMM19]) for solving quadratically regularised optimal transport

[^LMM19]: Lorenz, Dirk A., Paul Manns, and Christian Meyer. [Quadratically regularized optimal transport.](https://doi.org/10.1007/s00245-019-09614-w) Applied Mathematics & Optimization 83.3 (2021): 1919-1949.

See also: [`QuadraticOTNewton`](@ref), [`quadreg`](@ref)
"""
struct QuadraticOTNewton{T<:Real,K<:Real,D<:Real} <: QuadraticOT
    θ::T
    κ::K
    δ::D
    armijo_max::Int
end

"""
    QuadraticOTNewton(;θ = 0.1, κ = 0.5, δ = 1e-5, armijo_max = 50)

Semi-smooth Newton method (Algorithm 2 of Lorenz et al. 2019) with Armijo parameters `θ` and `κ`, and conjugate gradient regularisation `δ`.
`armijo_max` sets the maximum number of Armijo step size trials.

See also: [`QuadraticOTNewton`](@ref), [`quadreg`](@ref)
"""
function QuadraticOTNewton(; θ=0.1, κ=0.5, δ=1e-5, armijo_max=50)
    return QuadraticOTNewton(θ, κ, δ, armijo_max)
end

Base.show(io::IO, ::QuadraticOTNewton) = print(io, "Semi-smooth Newton algorithm")

struct QuadraticOTNewtonCache{U,V,C,P,GT,X}
    u::U
    v::V
    δu::U
    δv::V
    σ::C
    γ::P
    G::GT
    x::X
    M::Int
    N::Int
end

function build_cache(
    ::Type{T},
    ::QuadraticOTNewton,
    μ::AbstractVector,
    ν::AbstractVector,
    C::AbstractMatrix,
    ε::Real,
) where {T}
    # create and initialize dual potentials
    u = similar(μ, T, size(μ, 1))
    v = similar(ν, T, size(ν, 1))
    fill!(u, zero(T))
    fill!(v, zero(T))
    δu = similar(u, T)
    δv = similar(v, T)
    # intermediate variables (don't need to be initialised)
    σ = similar(C, T)
    γ = similar(C, T)
    M = size(μ, 1)
    N = size(ν, 1)
    G = similar(u, T, M + N, M + N)
    fill!(G, zero(T))
    # initial guess for conjugate gradient
    x = similar(u, T, M + N)
    fill!(x, zero(T))
    return QuadraticOTNewtonCache(u, v, δu, δv, σ, γ, G, x, M, N)
end

function check_convergence(
    μ::AbstractVector,
    ν::AbstractVector,
    cache::QuadraticOTNewtonCache,
    convergence_cache::QuadraticOTConvergenceCache,
    atol::Real,
    rtol::Real,
)
    γ = cache.γ
    norm_diff = max(
        norm(vec(sum(γ; dims=2)) .- μ, Inf), norm(vec(sum(γ; dims=1) .- ν), Inf)
    )
    isconverged =
        norm_diff <
        max(atol, rtol * max(convergence_cache.norm_source, convergence_cache.norm_target))
    return isconverged, norm_diff
end

function descent_dir!(solver::QuadraticOTSolver{<:QuadraticOTNewton})
    # unpack solver
    eps = solver.eps
    C = solver.C
    μ = solver.source
    ν = solver.target
    cache = solver.cache
    # unpack cache
    u = cache.u
    v = cache.v
    δu = cache.δu
    δv = cache.δv
    σ = cache.σ
    γ = cache.γ
    G = cache.G
    x = cache.x
    M = cache.M
    N = cache.N
    # Armijo parameters
    δ = solver.alg.δ

    # setup intermediate variables
    @. γ = u + v' - C
    @. σ = γ ≥ 0
    @. γ = NNlib.relu(γ) / eps

    # setup kernel matrix G
    # G[diagind(G)[1:M]] .= vec(sum(σ; dims=2))
    G[1:M, 1:M] .= Diagonal(vec(sum(σ; dims=2)))
    # G[diagind(G)[(M + 1):end]] .= vec(sum(σ; dims=1))
    G[(M + 1):end, (M + 1):end] .= Diagonal(vec(sum(σ; dims=1)))
    G[1:M, (M + 1):end] .= σ
    G[(M + 1):end, 1:M] .= σ'
    # G[diagind(G)] .+= δ # regularise cg
    G += δ * I

    # cg step
    b = -eps * vcat(vec(sum(γ; dims=2)) .- μ, vec(sum(γ; dims=1)) .- ν)
    cg!(x, G, b)
    δu .= x[1:M]
    return δv .= x[(M + 1):end]
end

function descent_step!(solver::QuadraticOTSolver{<:QuadraticOTNewton})
    # unpack solver
    eps = solver.eps
    C = solver.C
    μ = solver.source
    ν = solver.target
    cache = solver.cache
    # unpack cache
    u = cache.u
    v = cache.v
    δu = cache.δu
    δv = cache.δv
    γ = cache.γ

    # Armijo parameters
    θ = solver.alg.θ
    κ = solver.alg.κ
    armijo_max = solver.alg.armijo_max
    armijo_counter = 0

    # dual objective
    function Φ(u, v, μ, ν, C, ε)
        return norm(NNlib.relu.(u .+ v' .- C))^2 / 2 - ε * dot(μ, u) - ε * dot(ν, v)
    end

    # compute directional derivative
    d = -eps * (dot(δu, μ) + dot(δv, ν)) + eps * dot(γ, δu .+ δv')
    t = 1
    Φ0 = Φ(u, v, μ, ν, C, eps)
    while (armijo_counter < armijo_max) &&
        (Φ(u + t * δu, v + t * δv, μ, ν, C, eps) ≥ Φ0 + t * θ * d)
        t = κ * t
        armijo_counter += 1
    end
    u .= u + t * δu
    return v .= v + t * δv
end

function solve!(solver::QuadraticOTSolver{<:QuadraticOTNewton})
    # unpack solver
    μ = solver.source
    ν = solver.target
    atol = solver.atol
    rtol = solver.rtol
    maxiter = solver.maxiter
    check_convergence = solver.check_convergence
    cache = solver.cache
    convergence_cache = solver.convergence_cache

    isconverged = false
    to_check_step = check_convergence
    for iter in 1:maxiter
        # compute descent direction
        descent_dir!(solver)
        # Newton step
        descent_step!(solver)
        # check source marginal
        # always check convergence after the final iteration
        to_check_step -= 1
        if to_check_step == 0 || iter == maxiter
            # reset counter
            to_check_step = check_convergence

            isconverged, abserror = OptimalTransport.check_convergence(
                μ, ν, cache, convergence_cache, atol, rtol
            )
            @debug string(solver.alg) *
                   " (" *
                   string(iter) *
                   "/" *
                   string(maxiter) *
                   ": absolute error of source marginal = " *
                   string(maximum(abserror))

            if isconverged
                @debug "$(solver.alg) ($iter/$maxiter): converged"
                break
            end
        end
    end

    if !isconverged
        @warn "$(solver.alg) ($maxiter/$maxiter): not converged"
    end

    return nothing
end

# interface

function plan(solver::QuadraticOTSolver{<:QuadraticOTNewton})
    cache = solver.cache
    γ = NNlib.relu.(cache.u .+ cache.v' .- solver.C) / solver.eps
    return γ
end


text:
"""
Some utilities functions routinely used in MR simulations.
"""
module utils

using ..MRphy
using ..Unitful, ..UnitfulMR
#= General =#
export CartesianLocations
"""
    CartesianLocations(dim::Dims, doShift::Bool=true)

Retuns a `(prod(dim),length(dim))` sized array of grid locations. `doShift`
shifts the locations to be consistent with `fftshift`.

# Examples
```julia-repl
julia> loc = [CartesianLocations((2,2)), CartesianLocations((2,2),false)]
2-element Array{Array{Int64,2},1}:
 [-1 -1; 0 -1; -1 0; 0 0]
 [1 1; 2 1; 1 2; 2 2]
```
"""
CartesianLocations(dim::Dims, doShift::Bool=true) = doShift ?
  Array(hcat(collect.(Tuple.(CartesianIndices(dim).-ctrSub(dim)))...)') :
  Array(hcat(collect.(Tuple.(CartesianIndices(dim)))...)')

export ctrSub
"""
    ctrSub(dim::Dims) = CartesianIndex(dim .÷ 2 .+ 1)

As a separate function, ensure consistent behaviour of getting ::CartesianIndex
to the center of a Nd-Array of size `dim`.
This `center` should match `fftshift`'s `center`.

See also: [`ctrInd`](@ref).

# Notes:
The function may be removed once julia FFT packages provides this functionality.
"""
ctrSub(dim::Dims) = CartesianIndex(dim .÷ 2 .+ 1)

export ctrInd
"""
    ctrInd(dim::Dims) = sum((dim.÷2) .* [1; cumprod([dim[1:end-1]...])])+1

As a separate fn, ensure consistent behariour of getting the linear index to the
center of a Nd-array of size `dim`.
This `center` should match `fftshift`'s `center`.

See also: [`ctrSub`](@ref).
"""
ctrInd(dim::Dims) = sum((dim.÷2) .* [1; cumprod([dim[1:end-1]...])])+1

#= MR =#
export k2g
"""
    k2g(k::TypeND(K0D,:), isTx::Bool=false; dt::T0D=4e-6u"s", γ::Γ0D=γ¹H)
Gradient, `g`, of the `TxRx` k-space, (trasmit/receive, excitation/imaging).

# Usage
*INPUTS*:
- `k::TypeND(K0D, :)` (nSteps, Nd...), Tx or Rx k-space, w/ unit u"cm^-1".
- `isTx::Bool`, if `true`, compute transmit k-space, `k`, ends at the origin.
*KEYWORDS*:
- `dt::T0D` (1,), gradient temporal step size, i.e., dwell time.
- `γ::Γ0D` (1,), gyro-ratio.
*OUTPUTS*:
- `g::TypeND(GR0D, :)` (nSteps, Nd...), gradient

# Note
The function asserts if `k` ends at the origin for `isTx==true`.

See also: [`g2k`](@ref), [`g2s`](@ref).
"""
k2g(k::TypeND(K0D,:), isTx::Bool=false; dt::T0D=4e-6u"s", γ::Γ0D=γ¹H) =
  (isTx&&any(ustrip.(selectdim(k,1,size(k,1))) .!=0)) ?
    error("Tx `k` must end at 0") : [selectdim(k,1,1:1); diff(k,dims=1)]/(γ*dt)

export g2k
"""
    g2k(g::TypeND(GR0D,:); isTx::Bool=false, dt::T0D=4e-6u"s", γ::Γ0D=γ¹H)
Compute k-space from gradient.

# Usage
*INPUTS*:
- `g::TypeND(GR0D, :)` (nSteps, Nd...), gradient
- `isTx::Bool`, if `true`, compute transmit k-space, `k`, ends at the origin.
*KEYWORDS*:
- `dt::T0D` (1,), gradient temporal step size, i.e., dwell time.
- `γ::Γ0D` (1,), gyro-ratio.
*OUTPUTS*:
- `k::TypeND(K0D, :)` (nSteps, Nd...), k-space, w/ unit u"cm^-1".

See also: [`k2g`](@ref), [`g2s`](@ref).
"""
g2k(g::TypeND(GR0D,:), isTx::Bool=false; dt::T0D=4e-6u"s", γ::Γ0D=γ¹H) =
  γ*dt*cumsum(g,dims=1) |> k->isTx ? k.-selectdim(k,1,size(k,1):size(k,1)) : k

export g2s
"""
    g2s(g::TypeND(GR0D,:); dt::T0D=4e-6u"s")
Slew rate `sl`, of the gradient, `g`.

# Usage
*INPUTS*:
- `g::TypeND(GR0D, :)` (nSteps, Nd...)
*KEYWORDS*:
- `dt::T0D` (1,), gradient temporal step size, i.e., dwell time.
*OUTPUTS*:
- `sl::TypeND(Quantity{<:Real, 𝐁/𝐋/𝐓, :)` (nSteps, Nd...), slew rate

See also: [`g2k`](@ref), [`k2g`](@ref).

# Note
No `s2g` is provided for the moment.
"""
g2s(g::TypeND(GR0D,:);dt::T0D=4e-6u"s") = [selectdim(g,1,1:1);diff(g,dims=1)]/dt

end # module utils


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/LGforcing.jl")

    add_comp!(m, LGforcing)

    set_param!(m, :LGforcing, :c_LGconcentration, readpagedata(m, "test/validationdata/$valdir/c_LGconcentration.csv"))

    # run Model
    run(m)

    forcing = m[:LGforcing,:f_LGforcing]
    forcing_compare = readpagedata(m, "test/validationdata/$valdir/f_LGforcing.csv")

    @test forcing ≈ forcing_compare rtol = 1e-3
end


text:
using PyCall

torch = pyimport("torch")
np = pyimport("numpy")

"""
Note: the below code is a minimal test of using `cvxpylayers` and `torch` to
incorporate "differentiable convex programming".
Unfortunately, differentiable convex programming seems not realised in pure Julia packages.
"""

"""
Test for cvxpylayers
"""
function main()
    py"""
    import cvxpy as cp
    import torch
    from cvxpylayers.torch import CvxpyLayer
    import torch.optim as optim
    import copy

    n, m = 2, 3
    x = cp.Variable(n)
    A = cp.Parameter((m, n))
    b = cp.Parameter(m)
    constraints = [x >= 0]
    objective = cp.Minimize(0.5 * cp.pnorm(A @ x - b, p=1))
    problem = cp.Problem(objective, constraints)
    assert problem.is_dpp()

    cvxpylayer = CvxpyLayer(problem, parameters=[A, b], variables=[x])
    A_tch = torch.randn(m, n, requires_grad=True)
    b_tch = torch.randn(m, requires_grad=True)

    # solve the problem
    solution, = cvxpylayer(A_tch, b_tch)
    print(A_tch)
    A_tch_prev = copy.deepcopy(A_tch)

    # compute the gradient of the sum of the solution with respect to A, b
    optimiser = optim.SGD((A_tch, b_tch), lr=1e3)
    optimiser.zero_grad()
    solution.sum().backward()
    print(A_tch.grad)
    optimiser.step()
    print(A_tch)
    print(A_tch == A_tch_prev)
    """
end

"""
Test for network construction and inference
"""
function main2(n)
    py"""
    import torch
    import torch.nn as nn
    import torch.nn.functional as F

    class Net(nn.Module):
        def __init__(self):
            super(Net, self).__init__()
            self.fc = nn.Linear($$n, 10)

        def forward(self, x):
            x = self.fc(x)
            return x

    """
    # a = torch.from_numpy(np.array([1, 2, 3, 4, 5.], dtype=np.float32))
    a = torch.from_numpy(PyObject(np.random.randn(n)).astype(np.float32))
    net = py"Net"()
    net(a)
end


text:
import ThreeJS
import Colors: colormap

# Run in Escher to get updatable surf and mesh plots

function surf(f::Function)
 ThreeJS.outerdiv() <<
            (ThreeJS.initscene() <<
                [
                    ThreeJS.mesh(0.0, 0.0, 0.0) <<
                    [
                        ThreeJS.parametric(100,100,-10:10, -10:10,f),
                        ThreeJS.material(Dict(:kind=>"lambert",:color=>"white",
                        :colorkind=>"vertex"));
                    ],
                    ThreeJS.ambientlight(),
                    ThreeJS.camera(30.0, 30.0, 30.0)
                ]
            )
end

function mesh(f::Function)
 ThreeJS.outerdiv() <<
            (ThreeJS.initscene() <<
                [
                    ThreeJS.meshlines(50, 50, -10:10, -10:10, f),
                    ThreeJS.ambientlight(),
                    ThreeJS.camera(30.0, 30.0, 30.0);
                ]
            )
end
main(window) =  begin
    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
    push!(window.assets, "codemirror")
    push!(window.assets, "layout2")
    inp = Signal(Dict())
    s = sampler() # A thing that lets you watch widgets/behaviors upon updates to other behaviors
    default = "(x,y) -> sin(x) * cos(y)"
    editor = watch!(s, :code, codemirror(default))
    code_cell = trigger!(s, :submit, keypress("ctrl+enter shift+enter", editor))
    t, plots = wire(
                    tabs(["Surf";"Mesh";]),
                    pages(
                    [
                        map(inp) do f
                            fn = get(f,:code,default)
                            try
                                surf(eval(parse(fn)))
                            catch
                                Elem(:div, "Something went wrong. Please check your syntax and try again. Contact rohitvarkey@gmail.com for more assistance.")
                            end
                        end;
                        map(inp) do f
                            fn = get(f,:code,default)
                            try
                                mesh(eval(parse(fn)))
                            catch
                                Elem(:div, "Something went wrong. Please check your syntax and try again. Contact rohitvarkey@gmail.com for more assistance.")
                            end
                        end;
                    ]
                    ),
                    :tabschannel,
                    :selected
                )
    vbox(
            md"""Enter an anonymous function with 2 variables.
                `ctrl+enter` or `shift+enter` to redraw the plot.
                Use the mouse the drag, zoom and pan.
                The function is plotted with x and y between -10 and 10 and
                with 100 steps and 50 steps in both axes for the surf and the
                mesh respectively. Try resizing the browser if you cant see a codebox""",
            intent(s, code_cell) >>> inp,
            t, plots
        ) |> pad(2em)
end


text:
module CatElements
export ThElements, AbstractElements, Elements, elements, inverse_elements

using DataStructures: OrderedDict

using ..CSets, ..FinSets
using ...Present, ...Theories
using ...Theories: Category, ob, hom, dom_nums, codom_nums

@present ThElements(FreeSchema) begin
  (El, Arr, Ob, Hom)::Ob
  Name::AttrType
  src::Hom(Arr, El)
  tgt::Hom(Arr, El)
  dom::Hom(Hom, Ob)
  cod::Hom(Hom, Ob)
  πₑ::Hom(El, Ob)
  πₐ::Hom(Arr, Hom)
  nameo::Attr(Ob, Name)
  nameh::Attr(Hom, Name)

  src⋅πₑ == πₐ⋅dom
  tgt⋅πₑ == πₐ⋅cod
end

@abstract_acset_type AbstractElements
@acset_type Elements(ThElements, index=[:src, :tgt, :πₑ, :πₐ]) <: AbstractElements

"""    elements(X::AbstractACSet)

construct the category of elements from an ACSet. This only correctly handles the CSet part.
This transformation converts an instance of C into a Graph homomorphism. The codomain of the
homomorphism is a graph shaped like the schema. This is one half of the isomorphism between
databases and knowledge graphs.
"""
function elements(X::StructACSet{S}) where S
  Y = Elements{Symbol}()

  obs = ob(S)
  add_parts!(Y, :Ob, length(obs), nameo=obs)
  els = map(enumerate(obs)) do (i,c)
    add_parts!(Y, :El, nparts(X, c), πₑ = i)
  end

  add_parts!(Y, :Hom, length(hom(S)), dom=dom_nums(S), cod=codom_nums(S), nameh=hom(S))
  map(enumerate(zip(hom(S), dom_nums(S), codom_nums(S)))) do (i, (f, ci, di))
    c = obs[ci]
    add_parts!(Y, :Arr, nparts(X, c), src=els[ci], tgt=view(els[di], X[f]), πₐ=i)
  end
  return Y
end

"""    elements(f::ACSetTransformation)

Apply category of elements functor to a morphism f: X->Y.
This relies on the fact `elements` of an object puts El components from the same
Ob in a contiguous index range.
"""
function elements(f::ACSetTransformation{S}) where S
  X, Y = elements.([dom(f), codom(f)])

  # Apply offset to homomorphism data
  offs = map((parts(X, :Ob))) do i
    off=findfirst(==(i), Y[:πₑ])
    isnothing(off) ? 0 : off-1
  end
  pts = vcat([collect(f[o]).+off for (o, off) in zip(ob(S), offs)]...)
  # *strict* ACSet transformation uniquely determined by its action on vertices
  return only(homomorphisms(X, Y; initial=Dict([:El=>pts])))
end

"""    inverse_elements(X::AbstractElements, typ::StructACSet)
Compute inverse grothendieck transformation on the result of `elements`.
Does not assume that the elements are ordered.
Rather than dynamically create a new ACSet type, it requires any instance of
the ACSet type that it's going to try to create

If the typed graph tries to assert conflicting values for a foreign key, fail.
If no value is specified for a foreign key, the result will have 0's.
"""
function inverse_elements(X::AbstractElements, typ::StructACSet)
  res = typeof(typ)()
  o_ids = ob_ids(X)
  for (o, is) in pairs(o_ids)
    add_parts!(res, o, length(is))
  end
  cols = [[:πₐ,:nameh], :src, :tgt, [:src,:πₑ, :nameo], [:tgt, :πₑ, :nameo]]
  for (h, s, t, so, to) in zip([X[col] for col in cols]...)
    src_ind, tgt_ind = o_ids[so][s], o_ids[to][t]
    err = "Not a discrete opfibration: $h: $so#$s sent to multiple values"
    res[src_ind, h] == 0 || error(err)
    set_subpart!(res, src_ind, h, tgt_ind)
  end
  return res
end

"""
Determine the (partial) mapping from Elements to indices for each component
E.g. [V,E,V,V,E] ⟶ [V=>{1↦1, 3↦2, 4↦3}, E=>{2↦1, 5↦2}]
"""
function ob_ids(X::AbstractElements)
  obs = Dict([o => findall(i -> X[i,[:πₑ,:nameo]] == o, parts(X,:El))
              for o in X[:nameo]])
  return Dict([o => Dict(v=>i for (i,v) in enumerate(is)) for (o, is) in obs])
end

"""    inverse_elements(X::AbstractElements, typ::StructACSet)

Compute inverse grothendieck transformation on a morphism of Elements
"""
function inverse_elements(f::ACSetTransformation, typ::StructACSet)
  iX, iY = [inverse_elements(x, typ) for x in [dom(f), codom(f)]]
  oX, oY = ob_ids.([dom(f), codom(f)])
  comps = map(dom(f)[:nameo]) do ob
    offY = oY[ob]
    x_ids = sort(collect(keys(oX[ob])))
    ob => [offY[f[:El](xi)] for xi in x_ids]
  end
  return ACSetTransformation(iX, iY; Dict(comps)...)
end

"""    presentation(X::AbstractElements)

convert a category of elements into a new schema. This is useful for generating large schemas
that are defined as the category of elements of a specified C-Set. For example, the schema for
Bipartite Graphs is the category of elements of the graph with 2 vertices and 1 edge. The 2 vertices
will get mapped to elements `V_1, V_2` and the one edge will be `E_1` and the source and target maps
will be mapped to `src_E_1, tgt_E_1`.
"""
function presentation(X::AbstractElements)
  P = Presentation(FreeSchema)
  obs = OrderedDict{Tuple{Symbol, Int}, Any}()
  for a in 1:nparts(X, :Ob)
    class = X[a, :nameo]
    for (i,j) in enumerate(incident(X, a, :πₑ))
      ob = Ob(FreeSchema, Symbol("$(class)_$i"))
      obs[(class, j)] = ob
    end
  end
  add_generators!(P, values(obs))
  homs = OrderedDict{Tuple{Symbol, Int}, Any}()
  for f in 1:nparts(X, :Hom)
    class = X[f, :nameh]
    for (i,j) in enumerate(incident(X, f, :πₐ))
      domkey = (X[j, [:src, :πₑ, :nameo]], X[j, :src])
      dom = obs[domkey]
      codkey = (X[j, [:tgt, :πₑ, :nameo]], X[j, :tgt])
      codom = obs[codkey]
      homs[(class, j)] = Hom(Symbol("$(class)_$dom"), dom, codom)
    end
  end
  add_generators!(P, values(homs))
  return P, obs, homs
end
end


text:
module ModelingToolkitStandardLibrary

include("Blocks/Blocks.jl")
include("Electrical/Electrical.jl")
#include("Magnetic/Magnetic.jl")
include("Thermal/Thermal.jl")

end


text:
# parameter settings: \nu, \lambda, \sigma, \sigma_e)=(0.001, 0.02, 0.1, 0.001,
# parameter settings: sigma_epsilon = 0.001
# start value: x_0 = 5

using Pkg
using DataFrames
import Statistics.mean
import Statistics.var
using LinearAlgebra
using CSV
using PyPlot

# load functions
include(pwd()*"/src/SDEMEM OU neuron data/ou_sdemem.jl")
include(pwd()*"/src/SDEMEM OU neuron data/mcmc.jl")

data = CSV.read("data/SDEMEM OU neuron data/OUsim.dat"; header=0)

y = data[2:end,1]

dt = 1

time_grid = 1:dt:length(y)

PyPlot.figure()
PyPlot.plot(time_grid, y)

M_subjects = 1

ϕ = zeros(M_subjects,3) # lambda nu sigma

ϕ[1,:] = log.([0.001, 0.02, 0.1])

log_σ_ϵ = log(0.001)

nbr_particles = 1*ones(Int64, M_subjects)
ρ = 0.0

# set number of obs for each data set
T_vec = zeros(Int64, M_subjects)
for i = 1:M_subjects; T_vec[i] = length(y); end

u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
u_resample_old = [randn(T_vec[1], 2)]

for i in 2:M_subjects;
    append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
    append!(u_resample_old, [randn(T_vec[i],2)])
end

m = 1

dt
# Plot paths

loglik_est, x_hat_kalman = @time kalman_filter_test(y, exp(log_σ_ϵ),  ϕ[m,:], dt, true)

PyPlot.figure()
PyPlot.plot(time_grid, y, "g")
PyPlot.plot(time_grid, x_hat_kalman, "r--")

nbr_particles_bootstrap = 10000*ones(Int64, M_subjects)

# set number of obs for each data set
T_vec = zeros(Int64, M_subjects)
for i = 1:M_subjects; T_vec[i] = length(y); end

u_prop_old_bootstrap = [randn(nbr_particles_bootstrap[1], T_vec[1]+1)]
u_resample_old_bootstrap = [randn(T_vec[1], 2)]

for i in 2:M_subjects;
    append!(u_prop_old_bootstrap, [randn(nbr_particles_bootstrap[i],T_vec[i]+1)])
    append!(u_resample_old_bootstrap, [randn(T_vec[i],2)])
end

loglik, x,nbr_resamlping = @time particlefilter_test(y, exp(log_σ_ϵ), ϕ[m,:], dt, u_prop_old_bootstrap[m], u_resample_old_bootstrap[m], nbr_particles_bootstrap[m], true, true)

PyPlot.figure()
PyPlot.plot(time_grid, x[1:100,:]', "r--")
PyPlot.plot(time_grid, y, "g")

loglik, x,nbr_resamlping = diffusionbridgeparticlefilter_test(y, exp(log_σ_ϵ), ϕ[m,:], dt, u_prop_old[m], u_resample_old[m], nbr_particles[m], true, true)

PyPlot.figure()
PyPlot.plot(time_grid, y, "g")
PyPlot.plot(time_grid, x[:], "r--")

# run multiple loglik calcs

runs = 100

loglik_matrix = zeros(3, runs)

for i in 1:runs

    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y); end

    u_prop_old = [randn(nbr_particles[1], T_vec[1]+1)]
    u_resample_old = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old, [randn(nbr_particles[i],T_vec[i]+1)])
        append!(u_resample_old, [randn(T_vec[i],2)])
    end

    nbr_particles_bootstrap = 10000*ones(Int64, M_subjects)

    # set number of obs for each data set
    T_vec = zeros(Int64, M_subjects)
    for i = 1:M_subjects; T_vec[i] = length(y); end

    u_prop_old_bootstrap = [randn(nbr_particles_bootstrap[1], T_vec[1]+1)]
    u_resample_old_bootstrap = [randn(T_vec[1], 2)]

    for i in 2:M_subjects;
        append!(u_prop_old_bootstrap, [randn(nbr_particles_bootstrap[i],T_vec[i]+1)])
        append!(u_resample_old_bootstrap, [randn(T_vec[i],2)])
    end

    m = 1

    loglik_matrix[1,i], ~ = kalman_filter_test(y, exp(log_σ_ϵ),  ϕ[m,:], dt, true)
    loglik_matrix[2,i], ~ = particlefilter_test(y, exp(log_σ_ϵ), ϕ[m,:], dt, u_prop_old_bootstrap[m], u_resample_old_bootstrap[m], nbr_particles_bootstrap[m], true, true)
    loglik_matrix[3,i], ~ = diffusionbridgeparticlefilter_test(y, exp(log_σ_ϵ), ϕ[m,:], dt, u_prop_old[m], u_resample_old[m], nbr_particles[m], true, true)

end

println("Mean:")
println(round.(mean(loglik_matrix, dims = 2),digits=2))

println("Std:")
println(round.(std(loglik_matrix, dims = 2),digits=2))

max_stars_repo_path: src/SDEMEM OU neuron data/test_particle_filters_with_simulated_data.jl


text:
#=
Created on Thursday 16 Jan 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Solution for the color transfer.
=#

using STMO
using STMO.OptimalTransport
using Plots, LaTeXStrings

using Images, Colors

# function to subsample image
subsample(image, every=8) = image[1:every:size(image,1), 1:every:size(image,2)]

# load images
image1 = load("chapters/06.OptimalTransport/Figures/butterfly3.jpg")
image2 = load("chapters/06.OptimalTransport/Figures/butterfly2.jpg")

# subsample
image1 = subsample(image1, 8)
image2 = subsample(image2, 5)

# get colors, can also using sampling
colors1 = vec(image1)
colors2 = vec(image2)

# distance matrix
C = [colordiff(c1, c2) for c1 in colors1, c2 in colors2]

n, m = size(C)

# optimal transport
P = sinkhorn(C, ones(n)/n, ones(m)/m, λ=10)

mapdistr(X, P) = Diagonal(sum(P, dims=2)[:].^-1) * P * X

# map colors
image1transf = reshape(mapdistr(colors2, P), size(image1)...)
image2transf = reshape(mapdistr(colors1, P'), size(image2)...)

plot(
plot(image1, title="image 1", aspect_ratio=:equal),
plot(image2, title="image 2", aspect_ratio=:equal),
plot(image1transf, title="image 1 transfered", aspect_ratio=:equal),
plot(image2transf, title="image 2 transfered", aspect_ratio=:equal)
)

savefig("figures/colortransfer")

max_stars_repo_name: mferrodo/STMO

max_issues_repo_name: mferrodo/STMO

max_forks_repo_name: mferrodo/STMO


text:
@testset "build equivalence list" begin
    pred1 = Predicate((x,y) -> x == y)
    indl1 = -10:10
    eqC_1 = EquivalenceClasses(pred1, indl1)
    @test typeof(eqC_1) === EquivalenceClasses
    @test eqC_1.pred === pred1
    @test all(eqC_1.indl .== indl1)
    @test all(eqC_1.classes .== collect(1:length(indl1)))
    @test typeof(eqC_1.classes) <: AbstractArray
end

@testset "complete tests" begin
    pred1 = EquivalenceClassesConstructor.Predicate((x,y) -> (x .== -y || x .== (y .+ 10)))
    test_med_ind = (i for i in -300:300);
    #test_large_ind = ([i,j,k] for i in -30:20 for j in -10:30 for k in -20:30);
    test_large_ind = ([el...] for el in Iterators.product(-30:30, -10:30, -20:30))
end


text:
# # Wiring diagrams and syntactic expressions
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/wiring_diagrams/diagrams_and_expressions.ipynb)
#
# Morphisms in a monoidal category can be represented as syntactic expressions,
# such as $f \cdot g$ and $f \otimes g$, and also as [wiring
# diagrams](wiring_basics.ipynb), aka *string diagrams*. Catlab provides
# facilities for transforming between these two representations.

using Catlab.Doctrines, Catlab.WiringDiagrams
using Catlab.Graphics

function show_diagram(d::WiringDiagram)
  to_graphviz(d, direction=:horizontal, labels=false)
end

# ## Expressions to diagrams
#
# Converting a morphism expression to a wiring diagram is conceptually and
# algorithmically simple, because every expression determines a unique diagram.
#
# As a simple example, here is the monoidal product of two generators, $f$ and
# $g$, first as an expression (displayed using LaTeX) and then as a wiring
# diagram (displayed using Graphviz).

A, B, C, D, E = Ob(FreeCartesianCategory, :A, :B, :C, :D, :E)
f = Hom(:f, A, B)
g = Hom(:g, B, C)

expr = f ⊗ g
#-
show_diagram(to_wiring_diagram(expr))

# Here is a monoidal product of compositions:

h = Hom(:h, C, D)
k = Hom(:k, D, E)

expr = (f ⋅ g) ⊗ (h ⋅ k)
#-
show_diagram(to_wiring_diagram(expr))

# ## Diagrams to expressions
#
# Converting a wiring diagram to a syntactic expression is algorithmically more
# challenging, due to the fact that a single wiring diagram generally admits
# many different representations as an expression. Thus, a particular expression
# must be singled out.
#
# To bring this out, we define a function that round-trips a morphism expression
# to a wiring diagram and then back to an expression.

function roundtrip_expr(expr::FreeCartesianCategory.Hom)
  d = to_wiring_diagram(expr)
  to_hom_expr(FreeCartesianCategory, d)
end

# We can recover the expression just considered above:

roundtrip_expr((f ⋅ g) ⊗ (h ⋅ k))

# But here is a different expression that round-trips to the same thing:

roundtrip_expr((f ⊗ h) ⋅ (g ⊗ k))

# The equality of these two expressions,
#
# ```math
# (f \cdot g) \otimes (h \cdot k) = (f \otimes h) \cdot (g \otimes k),
# ```
#
# is the *interchange law* in a monoidal category. It says that composition and
# monoidal products can be interchanged. As this example shows, the conversion
# algorithm in Catlab favors products over composition, placing products towards
# the root of the expression tree wherever possible.
# Other laws can be discovered by this procedure. Since we are working in a
# [cartesian monoidal
# category](https://ncatlab.org/nlab/show/cartesian+monoidal+category),
# operations of copying, $\Delta_A: A \to A \otimes A$, and deleting,
# $\lozenge_A: A \to I$, are available.
#
# Consider the operation of copying the product $A \otimes B$.

expr = mcopy(A ⊗ B)
#-
show_diagram(add_junctions!(to_wiring_diagram(expr)))
#-
roundtrip_expr(expr)

# The equation just witnessed,
#
# ```math
# \Delta_{A \otimes B} = (\Delta_A \otimes \Delta_B) \cdot (1_A \otimes \sigma_{A,B} \otimes 1_B),
# ```
#
# is one of the *coherence laws* for cartesian products
# ([arXiv:0908.3347](https://arxiv.org/abs/0908.3347), Table 7). Another
# coherence law for products is
#
# ```math
# \lozenge_{A \otimes B} = \lozenge_A \otimes \lozenge_B.
# ```

expr = delete(A ⊗ B)
#-
roundtrip_expr(expr)


text:
__precompile__()

module VariableHeightMomentumControl

using RigidBodyDynamics
using RigidBodyDynamics.OdeIntegrators
using JuMP
using Gurobi
using StaticArrays

import RigidBodyDynamics: GenericJoint

export
    VariableHeightMomentumController,
    control,
    initialize_state!,
    set_desired_configuration!

include("util.jl")
include("initialization.jl")
include("control.jl")

end # module

max_stars_repo_name: tkoolen/VariableHeightMomentumControl

max_issues_repo_name: tkoolen/VariableHeightMomentumControl

max_forks_repo_name: tkoolen/VariableHeightMomentumControl


text:
# Extracting an subdomain from VLSV outputs and Plotting the Poynting vectors.
#
# Usage:
#   julia -t 4 demo_poynting_mt_pyplot.jl
# or
#   JULIA_NUM_THREADS=4 julia demo_poynting_mt_pyplot.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Statistics: mean, normalize
using LinearAlgebra: ×, ⋅
using Vlasiator
using Vlasiator: μ₀, prep2d, prep2dslice
using Glob, DSP, Polyester, Printf, PyPlot

"""
    extract_EM(files, range1, range2, pArgs, ndim=2; normal=:y, verbose=true)

Extract time-series of EM fields from `files` within `range1` in the 1st dim and `range1` in
the 2nd dim.
"""
function extract_EM(files, range1, range2, pArgs, ndim=2; normal=:y, verbose=true)
   @assert ndim ∈ (2,3) "Only support extracting EM fields from 2D/3D runs."
   nfile = length(files)
   e = zeros(Float32, 3, length(range1), length(range2), nfile)
   b = zeros(Float32, 3, length(range1), length(range2), nfile)

   if ndim == 2
      @batch for i in eachindex(files)
         verbose && println("i = $i/$nfile, file = $(files[i])")
         meta = load(files[i])

         e[1,:,:,i] = prep2d(meta, "vg_e_vol", 1)[range1, range2]
         e[2,:,:,i] = prep2d(meta, "vg_e_vol", 2)[range1, range2]
         e[3,:,:,i] = prep2d(meta, "vg_e_vol", 3)[range1, range2]
         b[1,:,:,i] = prep2d(meta, "vg_b_vol", 1)[range1, range2]
         b[2,:,:,i] = prep2d(meta, "vg_b_vol", 2)[range1, range2]
         b[3,:,:,i] = prep2d(meta, "vg_b_vol", 3)[range1, range2]
      end
   else
      @batch for i in eachindex(files)
         verbose && println("i = $i/$nfile, file = $(files[i])")
         meta = load(files[i])

         e[1,:,:,i] = prep2dslice(meta, "vg_e_vol", normal, 1, pArgs)[range1, range2]
         e[2,:,:,i] = prep2dslice(meta, "vg_e_vol", normal, 2, pArgs)[range1, range2]
         e[3,:,:,i] = prep2dslice(meta, "vg_e_vol", normal, 3, pArgs)[range1, range2]
         b[1,:,:,i] = prep2dslice(meta, "vg_b_vol", normal, 1, pArgs)[range1, range2]
         b[2,:,:,i] = prep2dslice(meta, "vg_b_vol", normal, 2, pArgs)[range1, range2]
         b[3,:,:,i] = prep2dslice(meta, "vg_b_vol", normal, 3, pArgs)[range1, range2]
      end
   end
   e, b
end

"Moving box average of vector `g` with box size `n`."
function moving_average(g::AbstractVector{<:AbstractFloat}, n::Int)
   nbackward = div(n,2)
   nforward = isodd(n) ? div(n,2) : div(n,2) - 1
   len = length(g)
   g_avg = similar(g)
   @inbounds @batch for i = 1:len
       lo = max(1, i - nbackward)
       hi = min(len, i + nforward)
       g_avg[i] = mean(@view g[lo:hi])
   end
   g_avg
end

"Extract the perturbation from vector `v` with moving box size `nbox`."
function detrend(v::AbstractVector{<:AbstractFloat}; nbox=length(v)÷12)
   v̄ = moving_average(v, nbox)
   dv = v .- v̄
end

"Extract the perturbation from 4D array `v` along the last dim with moving box size `nbox`."
function detrend(v::AbstractArray{<:AbstractFloat}; nbox=size(v,4)÷12)
   v̄ = similar(v)
   for idim = 1:3
      for j in axes(v, 3), i in axes(v, 2)
         v̄[idim,i,j,:] = @views moving_average(v[idim,i,j,:], nbox)
      end
   end
   dv = v .- v̄
   dv, v̄
end

"Band pass filter for vector field `var`."
function band_pass_filter(var, responsetype, designmethod)

   filter = digitalfilter(responsetype, designmethod)

   var_filtered = zeros(eltype(var), size(var,4), 3, size(var,2), size(var,3))

   for idim = 1:3
      for j in axes(var, 3), i in axes(var, 2)
         var_filtered[:,idim,i,j] = filtfilt(filter, var[idim,i,j,:])
      end
   end
   var_filtered
end

"""
    calc_poynting(de, db, b̄, it)

Return full Poynting vector, also parallel and perpendicular Poynting vector components to
the magnetic field. The perpendicular components are in-plane.
"""
function calc_poynting(de, db, b̄, it)

   s = zeros(3, size(de,3), size(de,4))
   # parallel and perpendicular Poynting vector magnitudes
   s_par  = zeros(size(de,3), size(de,4))
   s_perp = zeros(size(de,3), size(de,4))

   @views for j in axes(de,4), i in axes(de,3)
      # Poynting vector = dE × dB
      s[:,i,j] = de[it,:,i,j] × db[it,:,i,j] / μ₀

      # Transform into parallel and perpendicular direction w.r.t. B
      b̂ = normalize(b̄[:,i,j,it])
      s_par[i,j] = s[:,i,j] ⋅ b̂
      #s_perp[i,j] = norm(s[:,i,j] .- s_par[i,j] .* b̂)
      # in-plane perpendicular component only
      s_perp[i,j] = sqrt((s[1,i,j] - s_par[i,j]*b̂[1])^2 + (s[3,i,j] - s_par[i,j]*b̂[3])^2)
   end

   s, s_par, s_perp
end

"Output figures of Poynting vectors for each snapshot."
function plot_poynting(de_filtered, db_filtered, b̄, x1, x2, frequency_range)

   nt = size(de_filtered, 1)

   norm1 = let
      sparmax = frequency_range == "high" ? 5e-7 : 1e-5
      sparmin = -sparmax
      levels = matplotlib.ticker.MaxNLocator(nbins=255).tick_values(sparmin, sparmax)
      matplotlib.colors.BoundaryNorm(levels, ncolors=256, clip=true)
   end

   norm2 = let
      sperpmax = frequency_range == "high" ? 5e-7 : 1e-5
      sperpmin = 0.0
      levels = matplotlib.ticker.MaxNLocator(nbins=255).tick_values(sperpmin, sperpmax)
      matplotlib.colors.BoundaryNorm(levels, ncolors=256, clip=true)
   end

   stride = 10 # number of strides for quivers

   s, s_par, s_perp = calc_poynting(de_filtered, db_filtered, b̄, 1)

   fig, axs = plt.subplots(1, 2; figsize=(11,6), sharex=true, sharey=true,
      constrained_layout=true)

   axs[1].set_title(L"S_\parallel"; fontsize="large")
   axs[2].set_title(L"S_\perp"; fontsize="large")

   for ax in axs
      ax.set_aspect("equal")
      ax.set_xlabel(L"x1 [R_E]"; fontsize="large")
      ax.set_ylabel(L"x2 [R_E]"; fontsize="large")
   end

   c1 = axs[1].pcolormesh(x1, x2, s_par';
      cmap=matplotlib.cm.RdBu_r,
      shading="nearest",
      norm=norm1,
      )

   cb1 = colorbar(c1; ax=axs[1], extend="both")
   cb1.ax.set_ylabel(L"[W/m^2]"; fontsize="large")

   c2 = axs[2].pcolormesh(x1, x2, s_perp';
      cmap=matplotlib.cm.turbo,
      shading="nearest",
      norm=norm2,
      )

   cb2 = colorbar(c2; ax=axs[2], extend="max")
   cb2.ax.set_ylabel(L"[W/m^2]"; fontsize="large")

   q = axs[1].quiver(x1[1:stride:end], x2[1:stride:end],
      (s[1,:,:] ./ sqrt.(s[1,:,:].^2 .+ s[3,:,:].^2))[1:stride:end, 1:stride:end]',
      (s[3,:,:] ./ sqrt.(s[1,:,:].^2 .+ s[3,:,:].^2))[1:stride:end, 1:stride:end]';
      color="k")

   qb = axs[1].quiver(x1[1:stride:end], x2[1:stride:end],
      (b̄[1,:,:,1] ./ sqrt.(b̄[1,:,:,1].^2 .+ b̄[3,:,:,1].^2))[1:stride:end, 1:stride:end]',
      (b̄[3,:,:,1] ./ sqrt.(b̄[1,:,:,1].^2 .+ b̄[3,:,:,1].^2))[1:stride:end, 1:stride:end]';
      color="tab:purple")

   if frequency_range == "high"
      fig.suptitle("t = $(t[1]) s, [0.1, 1.0] Hz";
         fontsize="xx-large")
   else
      fig.suptitle("t = $(t[1]) s, [$(responsetype.w1), $(responsetype.w2)] Hz";
         fontsize="xx-large")
   end

   for it = 1:nt
      @info "it = $it"
      local s, s_par, s_perp
      outname = "poynting/poynting_$(lpad(it, 4, '0'))_band$frequency_range.png"
      isfile(outname) && continue
      s, s_par, s_perp = calc_poynting(de_filtered, db_filtered, b̄, it)

      c1.set_array(s_par')
      c2.set_array(s_perp')

      ŝx = @. s[1,:,:] / sqrt(s[1,:,:]^2 + s[3,:,:]^2)
      ŝz = @. s[3,:,:] / sqrt(s[1,:,:]^2 + s[3,:,:]^2)

      q.set_UVC(ŝx[1:stride:end, 1:stride:end]', ŝz[1:stride:end, 1:stride:end]')

      b̂x = @. b̄[1,:,:,it] / sqrt(b̄[1,:,:,it]^2 + b̄[3,:,:,it]^2)
      b̂z = @. b̄[3,:,:,it] / sqrt(b̄[1,:,:,it]^2 + b̄[3,:,:,it]^2)

      qb.set_UVC(b̂x[1:stride:end, 1:stride:end]', b̂z[1:stride:end, 1:stride:end]')

      if frequency_range == "high"
         fig.suptitle("t = $(t[it]) s, [0.1, 1.0] Hz";
            fontsize="xx-large")
      else
         fig.suptitle("t = $(t[it]) s, [$(responsetype.w1), $(responsetype.w2)] Hz";
            fontsize="xx-large")
      end
      savefig(outname; dpi=300, bbox_inches="tight")
   end
end

########## Main

files = glob("bulk*.vlsv", ".")
nfile = length(files)

frequency_range = "low" # filtered frequency range ∈ ("low", "high")

extent = [6., 16., -7., 7.] # [RE], default full domain: [-Inf32, Inf32, -Inf32, Inf32]

normal = :none # plane normal direction for 3D data, (:none, :x, :y, :z)

fs = 2.0 # sampling frequency, [Hz]

# WARNING: t may not be exact due to round-off errors in output time stamps!
ndim, pArgs, t, range1, range2 = let meta = load(files[1])
   pArgs = Vlasiator.set_args(meta, "vg_e_vol", EARTH; normal=:none)
   x1, x2 = Vlasiator.get_axis(pArgs)

   tstart = meta.time
   tend = load(files[end]).time

   ndims(meta), pArgs,
   range(round(tstart,digits=1), round(tend, digits=1), length=nfile),
   searchsortedfirst(x1, extent[1]):searchsortedlast(x1, extent[2]),
   searchsortedfirst(x2, extent[3]):searchsortedlast(x2, extent[4])
end

@assert ndim == 3 "3D not working yet!"

e, b = extract_EM(files, range1, range2, pArgs, ndim; normal, verbose=true)

designmethod = Butterworth(5)
if frequency_range == "high"
   responsetype = Highpass(0.1; fs)
   nbox = 40
elseif frequency_range == "low"
   responsetype = Bandpass(0.02, 0.067; fs)
   nbox = 200
end

de, _ = detrend(e; nbox)
db, b̄ = detrend(b; nbox)

de_filtered = band_pass_filter(de, responsetype, designmethod)
db_filtered = band_pass_filter(db, responsetype, designmethod)

x1 = range(extent[1], extent[2], length=length(range1))
x2 = range(extent[3], extent[4], length=length(range2))

plot_poynting(de_filtered, db_filtered, b̄, x1, x2, frequency_range)


text:
# The various exports to be included in the
export Chromosome,          # The abstract Chromosome type
    process,                # Processes a Chromosome
    find_active,            # Gets the list of active Nodes
    get_genes,              # Gets the genes from a Chromosome
    clone,                  # Clones a Chromrosome
    forward_connections,    # Gets a list of the forward connections
    get_output_trace        # Gets an output trace of the Chromrosome

# An abstract chromosome, to be defined
abstract type Chromosome end

#  An abstract node to be defined
abstract type Node end

# What is a CGP node look like?
mutable struct CGPNode <: Node
    connections::Array{Int64} # An array of connections coming in
    f::Function # The explicit function for the Node
    active::Bool # Whether or not it is currently active
    p::Float64 # A value to parameterize the node and connectivity
    output::Any # What output value it gives, if any
end

# Series of constuctors for CGP Nodes

# Basic constrcutor, giving inputs and a function
function CGPNode(ins::Array{Int64}, f::Function)
    CGPNode(ins, f, false, 1.0, 0.0)
end

# Extension, allowing it to be defined as active or not
function CGPNode(ins::Array{Int64}, f::Function, active::Bool)
    CGPNode(ins, f, active, 1.0, 0.0)
end

# Extension, adding a parameter value to be added
function CGPNode(ins::Array{Int64}, f::Function, active::Bool, p::Float64)
    CGPNode(ins, f, active, p, 0.0)
end

# I think this function snaps all values into the proper input range
# and then processes them
function snap(fc::Array{Float64}, p::Array{Float64})::Array{Int64}
    map(x->indmin(abs.(p-x)), fc)
end

# Processing a chromosome (individual ?) and inputs
function process(c::Chromosome, inps::Array)::Array{Float64}
    # For all input Nodes
    for i in 1:c.nin
        # Because we know these are all input values, we set their output to
        # the input value in the full CGP graph
        c.nodes[i].output = inps[i]
    end

    # Set the max length to the maximum length any given input value
    maxlength = maximum(map(length, inps))

    # Now we need to evaluate all the Nodes
    for n in c.nodes

        # Check if the Node is an active one (else we can ignore it)
        if n.active
            # Grab the output from that active Node
            output = n.output

            # If a Node is an input Node (?)
            if n.f == Config.f_input
                # The output is some sort of indexing into it (?)
                output = Config.index_in(inps, (n.p+1.0)/2.0)
            else
                # Grab the scaling factor, set it if it doesn't exist
                p = n.p
                if ~CGP.Config.weights
                    p = 1.0
                end

                # Compute the output based on function and connections and parameter
                output = CGP.Config.scaled(p * n.f(c.nodes[n.connections[1]].output,
                                                     c.nodes[n.connections[2]].output,
                                                     n.p))
            end

            # If output doesn't exist, set it to 0
            if length(output) == 0
                output = 0.0
            # Else set it to the max length of what we can handle
            elseif length(output) > maxlength
                output = output[1:maxlength]
            end

            # Update the Node's output
            n.output = output
        end
    end

    # Convert the outputs to an array of outputs
    outs = Array{Float64}(c.nout)

    # Because outputs have to be singular values, reduce to the mean
    for i in eachindex(outs)
        outs[i] = mean(c.nodes[c.outputs[i]].output)
    end

    # Return them out
    outs
end

# Function for activating Nodes with recurrences
function recur_active!(active::BitArray, connections::Array{Int64}, ind::Int64)::Void
    # Only if that Node is not active
    if ~active[ind]
        # Activate it
        active[ind] = true

        # Look at its connections
        for i in 1:2
            recur_active!(active, connections, connections[i, ind])
        end
    end
end

# Find the active sub-graph
function find_active(nin::Int64, outputs::Array{Int64}, connections::Array{Int64})::BitArray
    # Set everything to false
    active = falses(size(connections, 2)+nin)

    # Activate all the input Nodes
    active[1:nin] = true

    # For each of the final outputs, activate its connections
    for i in eachindex(outputs)
        recur_active!(active, connections, outputs[i])
    end

    # Set the inputs to false because we don't need that anymore
    active[1:nin] = false

    # Return them out
    active
end

# Another function for getting active subgraph
# Same as above, just doesn't deal with inputs
function find_active(outputs::Array{Int64}, connections::Array{Int64})::BitArray
    active = falses(size(connections, 2))
    for i in eachindex(outputs)
        recur_active!(active, connections, outputs[i])
    end
    active
end

# Resets a chromrosome
function reset!(c::Chromosome)
    for n in c.nodes
        n.output = 0.0
    end
end

# Has to implementeed down the line
function get_positions(c::Chromosome)
    error("Must be implemented in subclass")
end

# A recursive function for tracing output
function recur_output_trace(c::Chromosome, ind::Int64, visited::Array{Int64})
    # Check if a Node has been visited
    if ~contains(==, visited, ind)

        # add it to the list of visited Nodes
        append!(visited, [ind])

        # If the Node is not an input Node and is not an Input function
        if (ind > c.nin) && (c.nodes[ind].f != Config.f_input)

            # Visit all its connections
            for i in c.nodes[ind].connections
                recur_output_trace(c, i, visited)
            end
        end
    end

    # Return out the visited Nodes
    visited
end

# Calls the output trace for a given chromosome and a specific output index
function get_output_trace(c::Chromosome, output_ind::Int64)
    # similar to decode, just return a list of node indices that determine the output
    recur_output_trace(c, c.outputs[output_ind], Array{Int64}(0))
end

# Gets a full output trace
function get_output_trace(c::Chromosome, outputs::Array{Int64})
    if length(outputs) > 0
        return unique(reduce(vcat, map(x->get_output_trace(c, x), outputs)))
    else
        return Array{Int64}(0)
    end
end

# Gets a full output trace by implicitly collecting its outputs
function get_output_trace(c::Chromosome)
    # same as indexing over active nodes
    get_output_trace(c, collect(1:c.nout))
end

# Is a subclass implementation
function node_genes(c::Chromosome)
    # number of genes per node (position, x, y, f, param)
    error("Must be implemented in subclass")
end

# Gets all the genes for a chromosome and Node
function get_genes(c::Chromosome, node_id::Int64)
    c.genes[(c.nin+c.nout)+((node_id-1-c.nin)*node_genes(c))+(1:node_genes(c))]
end

# Gets all the genes for a chromosome and an array of Nodes
function get_genes(c::Chromosome, nodes::Array{Int64})
    if length(nodes) > 0
        return reduce(vcat, map(x->get_genes(c, x), nodes))
    else
        return Array{Int64}(0)
    end
end

# Updates Genes for a chromosome at a given Node
function set_genes!(c::Chromosome, node_id::Int64, genes::Array{Float64})
    c.genes[(c.nin+c.nout)+((node_id-1-c.nin)*node_genes(c))+(1:node_genes(c))] = genes
end

# Gets all the forward connections
function forward_connections(c::Chromosome)
    # Get the list of connections for all Nodes
    connections = [[i] for i in 1:length(c.nodes)]

    #  Now iterate over each Node
    for ci in eachindex(c.nodes)
        # Iterate over its connections
        for i in 1:2
            # Get the connection
            conn = c.nodes[ci].connections[i]

            # IF it's greater than 0, we want it (it exists)
            if conn > 0
                # Check if it's already in the connection list
                # If no, include it
                if ~(contains(==, connections[ci], conn))
                    append!(connections[ci], [conn])
                end

                # For each additional connection, check if its in the list
                # and include it
                for j in eachindex(connections[conn])
                    if ~(contains(==, connections[ci], j))
                        append!(connections[ci], [j])
                    end
                end
            end
        end
    end

    # Return out all the connections
    connections
end

# Clones a chromosome
function clone(c::Chromosome)
    # call gene constructor method
    typeof(c)(deepcopy(c.genes), c.nin, c.nout)
end


text:
module SfThermo

using Unitful, Documenter

DocMeta.setdocmeta!(SfThermo, :DocTestSetup, :(using Unitful, ..SfThermo); recursive=true)

import PhysicalConstants.CODATA2018: k_B, N_A, σ
import ..SfUnits: Energy, SpecificEnergy, MolarMass

export maxwell_boltzmann_peak_speed, rms_thermal_velocity, black_body_radiant_flux, black_body_radiated_power,
	ideal_gas_kinetic_energy, ideal_gas_cooling_time, convert_energy, convert_temperature

"""
    maxwell_boltzmann_peak_speed(m::Unitful.Mass, T::Unitful.Temperature)

Compute the velocity at the peak of the Maxwell-Boltzmann speed distribution for a gas with molecular mass `m` and temperature `T`.

``\\sqrt{\\frac{2k_BT}{m}}``

where ``k_B`` is the Boltzmann constant.
"""
maxwell_boltzmann_peak_speed(m::Unitful.Mass, T::Unitful.Temperature) = sqrt((2k_B*T)/m) |> u"m/s"

"""
    rms_thermal_velocoty(m::Unitful.Mass, T::Unitful.Temperature)

Compute RMS of the Maxwell-Boltzmann speed distribution for a gas with molecular mass `m` and temperature `T`.

``\\sqrt{\\frac{3k_BT}{m}}``

where ``k_B`` is the Boltzmann constant.
"""
rms_thermal_velocity(m::Unitful.Mass, T::Unitful.Temperature) = sqrt((3k_B*T)/m) |> u"m/s"

"""
    black_body_radiant_flux(t::Unitful.Temperature, ϵ::Real = 1, amb::Unitful.Temperature = 0u"K")

Radiated power per unit area across all frequencies from a black body with emissivity `ϵ` and temperature `t`
into an environment with ambient temperature `amb`.

``j^* = \\epsilon \\sigma(T_{hot}^4 - T_{ambient}^4)``

where ``k_B`` is Boltzmann's constant, ``\\sigma`` is the Stefan-Boltzmann constant and ``\\epsilon`` is the radiator emissivity.
"""
black_body_radiant_flux(t::Unitful.Temperature, ϵ::Real = 1, amb::Unitful.Temperature = 0u"K") = upreferred(ϵ * σ * (t^4 - amb^4))

"""
    black_body_radiated_power(t::Unitful.Temperature, a::Unitful.Area, ϵ::Real = 1, amb::Unitful.Temperature = 0u"K")

Radiated power across all frequencies from a black body with radiating area `a`, emissivity `ϵ` and temperature `t`
into an environment with ambient temperature `amb`.

``P = A \\epsilon \\sigma(T_{hot}^4 - T_{ambient}^4)``

where ``k_B`` is Boltzmann's constant, ``\\sigma`` is the Stefan-Boltzmann constant, ``A`` is the radiating area and
``\\epsilon`` is the radiator emissivity.
"""
black_body_radiated_power(t::Unitful.Temperature, a::Unitful.Area, ϵ::Real = 1, amb::Unitful.Temperature = 0u"K") = upreferred(ϵ * a * σ * (t^4 - amb^4))

"""
    ideal_gas_kinetic_energy(n, t::Unitful.Temperature)

Kinetic energy of `n` particles of an ideal gas at temperature `t`.

``E = N \\frac{3}{2}k_BT``

where ``k_B`` is Boltzmann's constant.
"""
ideal_gas_kinetic_energy(n, t::Unitful.Temperature) = 3n * k_B * t / 2 |> u"J"

"""
    ideal_gas_cooling_time(n, t_initial::Unitful.Temperature, t_final::Unitful.Temperature, a::Unitful.Area, ϵ::Real = 1)

Time taken for `n` particles of an ideal gas with radiating area `a` and emissivity `ϵ` to cool from a temperature
of `t_initial` to a temperature of `t_final`.

``t_{cooling} = \\frac{Nk_B}{2\\sigma \\epsilon A}\\left[ \\frac{1}{T_{final}^3} - \\frac{1}{T_{initial}^3} \\right]``

where ``k_B`` is Boltzmann's constant, ``\\sigma`` is the Stefan-Boltzmann constant, ``A`` is the radiating area, ``\\epsilon`` is
the radiator emissivity, and ``T_initial`` and ``T_final`` and the start and end temperatures.
"""
ideal_gas_cooling_time(n, t_initial::Unitful.Temperature, t_final::Unitful.Temperature, a::Unitful.Area, ϵ::Real = 1) =
	upreferred(((n * k_B) / (2ϵ * σ * a)) * ((1/t_final^3) - (1 / t_initial^3)))

"""
    convert_temperature(e::Energy)

Convert a temperature defined as eg. energy per atom `e` to kelvin.

# Example
```jldoctest
julia> convert_temperature(1u"eV")
11604.518121550082 K
```"""
convert_temperature(e::Energy) = e / k_B |> u"K"

"""
    convert_temperature(t::Unitful.Temperature)

Convert a temperature `t` in eg. kelvin to eV per atom.

# Example
```jldoctest
julia> convert_temperature(12000u"K")
1.0340799914574215 eV
```"""
convert_temperature(t::Unitful.Temperature) = t * k_B |> u"eV"

"""
    convert_energy(e::SpecificEnergy, mn::MolarMass)

Convert specific energy `e` for given material `mn` to eV per atom.

# Example
```jldoctest
julia> convert_energy(1e10u"erg/g", 60u"g/mol")
0.6218561793757305 eV
```"""
convert_energy(e::SpecificEnergy, mn::MolarMass) = e * mn / N_A |> u"eV"

"""
    convert_energy(e::Energy, mn::MolarMass)

Convert energy per atom `e` to specific energy for a material with given molar mass `mn`.

# Example
```jldoctest
julia> convert_energy(0.6218561793757305u"eV", 60u"g/mol")
999999.9999999999 J kg^-1
```"""
convert_energy(e::Energy, mn::MolarMass) = (e * N_A) / mn |> u"J/kg"

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
@testset "O(logn)" begin
    function latex_ground(args...)
        translate(-200, -100)
        background("white")
        sethue("black")
        fontsize(50)
    end

    function latex_blend(latex_string, fsize)
        orangeblue = blend(Point(50, 70), Point(50, 120), "orange", "blue")
        setblend(orangeblue)
        fontsize(fsize)
        latex(latex_string, Point(50, 70))
    end

    video = Video(400, 200)
    Background(1:35, latex_ground)
    runtime = Object((args...) -> latex(L"\mathcal{O}(\log{n})"))
    act!(runtime, Action(1:15, appear(:draw_text)))
    act!(runtime, Action(21:35, disappear(:draw_text)))
    Object((args...) -> latex_blend(L"\mathcal{O}\left(\frac{\log{x}}{2}\right)", 20))
    render(video; tempdirectory = "images", pathname = "")
    @test_reference "refs/ologn.png" load("images/0000000016.png")
    @test_reference "refs/ologn_mid.png" load("images/0000000007.png")
    @test_reference "refs/ologn_dis_mid.png" load("images/0000000027.png")
    for i in 1:35
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "latex pos in function" begin
    function latex_ground(args...)
        background("white")
        sethue("black")
        fontsize(30)
    end

    video = Video(400, 200)
    Background(1:1, latex_ground)
    Object((args...) -> latex(L"8")) # default fontsize 50
    render(video; tempdirectory = "images", pathname = "")
    @test_reference "refs/latex_8.png" load("images/0000000001.png")
    rm("images/0000000001.png")
end

@testset "LaTeX 3x3 matrix" begin
    function latex_ground(args...)
        translate(-200, -100)
        background("white")
        sethue("black")
        fontsize(20)
    end

    video = Video(400, 200)
    Background(1:1, latex_ground)
    Object(
        (args...) ->
            latex(L"$\begin{equation}\left[\begin{array}{ccc}1 & 2 & 3 \\4 & 5 & 6 \\7 & 8 & 9 \\\end{array}\right]\end{equation}$"),
    )
    render(video; tempdirectory = "images", pathname = "")
    @test_reference "refs/latex_3x3_matrix.png" load("images/0000000001.png")
    rm("images/0000000001.png")
end

@testset "LaTeX 8" begin
    function latex_ground(args...)
        translate(-200, -100)
        background("white")
        sethue("black")
        fontsize(50)
    end

    video = Video(400, 200)
    Background(1:1, latex_ground)
    Object((args...) -> latex(L"\mathcal{O}(\log{n})"))
    render(video; tempdirectory = "images", pathname = "")
    @test_reference "refs/ologn_circ.png" load("images/0000000001.png")
    rm("images/0000000001.png")
end


text:
"DSL for Financial Contracts."
module FinancialDSL

import BusinessDays
import InterestRates
import ForwardDiff
import OptimizingIR
import Distributions
import Roots

const CAL_BRL = BusinessDays.BRSettlement()

function __init__()
    BusinessDays.initcache(CAL_BRL)
end

include("BS.jl")
include("currencies/currencies.jl")
include("marketdata.jl")
include("core/core.jl")

end # module


text:
using CSV
using DataFrames

"""
    AgentsIO.populate_from_csv!(model, filename [, agent_type, col_map]; row_number_is_id, kwargs...)

Populate the given `model` using CSV data contained in `filename`. Use `agent_type` to
specify the type of agent to create (In the case of multi-agent models) or a function
that returns an agent to add to the model. The CSV row is splatted into the `agent_type`
constructor/function.

`col_map` is a `Dict{Symbol,Int}` specifying a mapping of keyword-arguments to row number.
If `col_map` is specified, the specified data is splatted as keyword arguments.

The keyword `row_number_is_id = false` specifies whether the row number will be passed as
the first argument (or as `id` keyword) to `agent_type`.

Any other keyword arguments are forwarded to `CSV.Rows`. If the `types` keyword is not
specified and `agent_type` is a struct, then the mapping from struct field to type will be used.
`Tuple{...}` fields will be suffixed with `_1`, `_2`, ... similarly to [`AgentsIO.dump_to_csv`](@ref)

For example,
```
struct Foo <: AbstractAgent
    id::Int
    pos::NTuple{2,Int}
    foo::Tuple{Int,String}
end

model = ABM(Foo, ...)
AgentsIO.populate_from_csv!(model, "test.csv")
```
Here, `types` will be inferred to be
```
Dict(
    :id => Int,
    :pos_1 => Int,
    :pos_2 => Int,
    :foo_1 => Int,
    :foo_2 => String,
)
```
It is not necessary for all these fields to be present as columns in the CSV. Any column
names that match will be converted to the appropriate type. There should exist a constructor
for `Foo` taking the appropriate combination of fields as parameters.

If `"test.csv"` contains the following columns: `pos_1, pos_2, foo_1, foo_2`, then `model`
can be populated as `AgentsIO.populate_from_csv!(model, "test.csv"; row_number_is_id = true)`.
"""
function populate_from_csv!(
    model::ABM{S,A},
    filename,
    agent_type::B = A,
    col_map::Dict{Symbol,Int} = Dict{Symbol,Int}();
    row_number_is_id = false,
    kwargs...,
) where {A,B<:Union{Type{<:A},Function},S}
    @assert(
        agent_type isa Function || !(agent_type isa Union),
        "agent_type cannot be a Union. It must be a Function or concrete subtype of AbstractAgent"
    )
    if !haskey(kwargs, :types) && isstructtype(agent_type)
        kwargs = (
            kwargs...,
            types = Dict(
                fieldname(agent_type, i) => fieldtype(agent_type, i) for
                i in 1:fieldcount(agent_type)
            ),
        )
        for (k, v) in kwargs.types
            if v <: Tuple && isconcretetype(v)
                len = length(fieldtypes(v))
                for i in 1:len
                    kwargs.types[Symbol(k, "_$i")] = fieldtypes(v)[i]
                end
            end
        end
    end

    if isempty(col_map)
        if row_number_is_id
            for (id, row) in enumerate(CSV.Rows(read(filename); kwargs...))
                add_agent_pos!(agent_type(id, row...), model)
            end
        else
            for row in CSV.Rows(read(filename); kwargs...)
                add_agent_pos!(agent_type(row...), model)
            end
        end
    else
        if row_number_is_id
            for (id, row) in enumerate(CSV.Rows(read(filename); kwargs...))
                add_agent_pos!(agent_type(; id, (k => row[v] for (k, v) in col_map)...), model)
            end
        else
            for row in CSV.Rows(read(filename); kwargs...)
                add_agent_pos!(agent_type(; (k => row[v] for (k, v) in col_map)...), model)
            end
        end
    end
end

"""
    AgentsIO.dump_to_csv(filename, agents [, fields]; kwargs...)

Dump `agents` to the CSV file specified by `filename`. `agents` is any iterable
sequence of types, such as from [`allagents`](@ref). `fields` is an iterable sequence of
`Symbol`s specifying which fields of each agent are dumped. If not explicitly specified,
it is automatically inferred using `eltype(agents)`. All `kwargs...` are forwarded
to `CSV.write`.

All `Tuple{...}` fields are flattened to multiple columns suffixed by `_1`, `_2`...
similarly to [`AgentsIO.populate_from_csv!`](@ref)

For example,
```
struct Foo <: AbstractAgent
    id::Int
    pos::NTuple{2,Int}
    foo::Tuple{Int,String}
end

model = ABM(Foo, ...)
...
AgentsIO.dump_to_csv("test.csv", allagents(model))
```
The resultant `"test.csv"` file will contain the following columns: `id`, `pos_1`, `pos_2`,
`foo_1`, `foo_2`.
"""
function dump_to_csv(filename, agents, fields = collect(fieldnames(eltype(agents))); kwargs...)
    atype = eltype(agents)
    data = DataFrame()
    for f in fields
        ftype = fieldtype(atype, f)
        if ftype <: Tuple && isconcretetype(ftype)
            flen = length(fieldtypes(ftype))
            for i in 1:flen
                data[!, Symbol(f, "_$i")] = [getproperty(a, f)[i] for a in agents]
            end
        else
            data[!, f] = [getproperty(a, f) for a in agents]
        end
    end

    CSV.write(filename, data; kwargs...)
end


text:
function  Bistochastic_Normalization!(X::Array{Float64,2},iterations::Int)

    n,m = size(X)

    for _ in 1:iterations
        for i = 1:n
            X[i,:] ./= norm(X[i,:],1)
        end

        for j = 1:m
            X[:,j] ./= norm(X[:,j],1)
        end
    end

end

function Bistochastic_Normalization!(x::Array{Float64,1},iterations::Int,n::Int,m::Int)

    X = reshape(x,n,m)

    for _ in 1:iterations
        for i = 1:n
            X[i,:] ./= norm(X[i,:],1)
        end

        for j = 1:m
            X[:,j] ./= norm(X[:,j],1)
        end
    end

end

function build_assignment(X::Array{Float64,2};use_colsort::Bool=false,
                          return_matrix=false)
    """-------------------------------------------------------------------------
        For a HOFASM_iteration matrix result X, computing a greedy matching on
      the maximum elements of the columns sequentially. if using the option
      use_colsort, will find the matching in the columns ordered by the largest
      element of each column. If return_matrix is true, a matrix is returned of
      the same dimensions as X, but X(i,j) =  1 iff i is matched to j.

      Notes:
      ------
      * Currently implemented for square matrices
    -------------------------------------------------------------------------"""
    n, m = size(X)
    selected_vertices = Set{Int}()
    match_size = minimum((n,m))
    matching = Array{Tuple{Int,Int},1}(undef,match_size)
    index = 1

    if use_colsort
        col_order = sortperm([maximum(X[:,i]) for i in 1:m],rev=true)
    else
        col_order = 1:match_size
    end

    for j in col_order #search over columns
        max_val = -Inf
        arg_max = -1
        for i = 1:n
            if (X[i,j] > max_val) && !(i in selected_vertices)
                max_val = X[i,j]
                arg_max = i
            end
        end
        matching[index] = (j,arg_max)
        index += 1
        push!(selected_vertices,arg_max)
        if index > match_size
            break
        end
    end
    if return_matrix
        Z = zeros(n,m)
        for (i,j) in matching
            Z[i,j] = 1
        end
        return Z
    else
        return matching
    end
end

function HOFASM_iterations(Bn_indices_list::Array{Array{Int,2},1},Bn_vals_list::Array{Array{Float64,1},1},
                           Hn_indices_list::Array{Array{Int,2},1},n::Int,m::Int,max_iterations::Int=12)

    x_k::Array{Float64,1} = rand(Float64,n*m)
    x_k ./= norm(x_k)
    beta = 30 # constant from the paper

    i = 1
    x_k_1::Array{Float64,1} = zeros(n*m)

    while true

        implicit_kronecker_model_marginalization_contraction!(Bn_indices_list,Bn_vals_list,
                                                              Hn_indices_list,x_k,x_k_1,n,m)

        max_elem = maximum(x_k_1)
        inflation_normalization = x -> exp(x*beta/max_elem)
        broadcast!(inflation_normalization,x_k_1,x_k_1)

        while true

            x_normalized = copy(x_k_1)
            Bistochastic_Normalization!(x_normalized,15,n,m)
  #          print(f"{norm(X_normalized - X_k_1,'fro')}")
            if norm(x_normalized - x_k_1) < 1e-4
                x_k_1 = x_normalized
                break
            else
                x_k_1 = copy(x_normalized)
            end
        end
#        return X_k_1

        x_k_1 ./= norm(x_k_1)
        res = norm(x_k_1 - x_k) / norm(x_k)

        #print("Finished iteration $i with difference $res")
        i += 1
        if i > max_iterations
            return x_k_1
        end

        if res < 1e-4
            return x_k_1
        else
            x_k = copy(x_k_1)
        end
    end
end

function HOFASM_iterations(matrices::Array{SparseMatrixCSC{Float64,Int64},1},n::Int,m::Int,max_iterations::Int=12)

    #    X_k = np.ones((n,m))
    x_k::Array{Float64,1} = rand(Float64,n*m)
    x_k ./= norm(x_k)
    beta = 30 # constant from the paper

    i = 1

    while true

        x_k_1::Array{Float64,1} = zeros(n*m)
        for Tn in matrices
            x_k_1 .+= Tn*x_k
        end

        max_elem = maximum(x_k_1)
        inflation_normalization = x -> exp(x*beta/max_elem)
        broadcast!(inflation_normalization,x_k_1,x_k_1)

        while true

            x_normalized = copy(x_k_1)
            Bistochastic_Normalization!(x_normalized,15,n,m)
  #          print(f"{norm(X_normalized - X_k_1,'fro')}")
            if norm(x_normalized - x_k_1) < 1e-4
                x_k_1 = x_normalized
                break
            else
                x_k_1 = copy(x_normalized)
            end
        end
#        return X_k_1

        x_k_1 ./= norm(x_k_1)
        res = norm(x_k_1 - x_k) / norm(x_k)

        #print("Finished iteration $i with difference $res")
        i += 1
        if i > max_iterations
            return x_k_1
        end

        if res < 1e-4
            println("tolerance met")
            return x_k_1
        else
            x_k = copy(x_k_1)
        end
    end
end

function HOFASM_iterations(tensor_pairs::Array{NTuple{2,SparseMatrixCSC{Float64,Int64}},1},n::Int,
                           m::Int,max_iterations::Int=12,tol=1e-4)

    #    X_k = np.ones((n,m))
    X_k::Array{Float64,2} = rand(Float64,n,m)
    X_k ./= norm(X_k)
    beta = 30 # constant from the paper

    X_k_1 = Array{Float64,2}(undef,n,m)
    i = 1

    preprocessed_pairs = [(findnz(H)...,B) for (H,B) in tensor_pairs]

    while true

#        HOFASM_contraction!_test(tensor_pairs,X_k,X_k_1)
        HOFASM_contraction!(preprocessed_pairs,X_k,X_k_1)

 #       SIMD_sparse_matrix_multiplication(Hn_matrices,Bn_matrices,X_k,X_k_1)
#        packed_X_k = pack(X_k')
 #       for (Hn,Bn) in tensor_pairs

#            X_k_1 .+= Hn*(Bn*X_k')'
#            X_k_1 .+= (Bn*(Hn*X_k)' )'

#            temp = unpack(Bn*packed_X_k)'
#            X_k_1 .+= unpack(Hn*pack(temp))
           #println(typeof(Hn*(Bn*X_k')'))
           # X_k_1 .+= Hn*(X_k*Bn')
  #      end

 #       println("done")

        max_elem = maximum(X_k_1)
        inflation_normalization = x -> exp(x*beta/max_elem)
        broadcast!(inflation_normalization,X_k_1,X_k_1)

        while true

            X_normalized = copy(X_k_1)
            Bistochastic_Normalization!(X_normalized,15)
  #          println("BN")
   #         print("{norm(X_normalized - X_k_1,'fro')}")
            if norm(X_normalized - X_k_1) < 1e-4
                X_k_1 = X_normalized
                break
            else
                X_k_1 = copy(X_normalized)
            end
        end
#        return X_k_1

        X_k_1 ./= norm(X_k_1)
        res = norm(X_k_1 - X_k) / norm(X_k)

  #      println("Finished iteration $i with difference $res")
        i += 1
        if i > max_iterations
            return X_k_1
        end

        if res < tol
            return X_k_1
        else
            X_k = copy(X_k_1)
        end
    end
end

#TODO: could generalize to take a contraction function to unify all graduated assignment routines
function HOM_graduated_assignment(marg_tensor::SparseMatrixCSC{Float64,Int64},m::Int=-1,max_iterations::Int=12,tol::Float64=1e-4)
    println(m)
    println(max_iterations)
    #    X_k = np.ones((n,m))
    n = size(marg_tensor,1)
    if m == -1
        m = Int(sqrt(n))
    end
    #m = Int(sqrt(n))
    x_k::Array{Float64,1} = rand(Float64,n)
    x_k ./= norm(x_k)
    beta = 30 # constant from the paper
    i = 1

    while true

        x_k_1 = marg_tensor*x_k

        max_elem = maximum(x_k_1)
        inflation_normalization = x -> exp(x*beta/max_elem)
        broadcast!(inflation_normalization,x_k_1,x_k_1)

        while true

            x_normalized = copy(x_k_1)
            Bistochastic_Normalization!(x_normalized,10,Int(n/m),m)
  #          print(f"{norm(X_normalized - X_k_1,'fro')}")
            if norm(x_normalized - x_k_1) < 1e-4
                x_k_1 = x_normalized
                break
            else
                x_k_1 = copy(x_normalized)
            end
        end
#        return X_k_1

        x_k_1 ./= norm(x_k_1)
        res = norm(x_k_1 - x_k) / norm(x_k)

        #print("Finished iteration $i with difference $res")
        i += 1
        if i > max_iterations
            return x_k_1
        end

        if res < 1e-4
            println("tolerance met")
            return x_k_1
        else
            x_k = copy(x_k_1)
        end
    end

end

#Sequential Second-order expansion for higher order matching
function SSoEfHOM(tensor_pairs::Array{NTuple{2,SparseMatrixCSC{Float64,Int64}},1},n::Int,
                  m::Int, max_iterations::Int=12,tol::Float64=1e-4)

    X_k::Array{Float64,2} = rand(Float64,n,m)
    beta = 30 # constant from the paper

    X_k_1 = Array{Float64,2}(undef,n,m)
    B_k_1 = Array{Float64,2}(undef,n,m)
    temp2 = Array{Float64,2}(undef,n,m)
    i = 1

    preprocessed_pairs = [(findnz(H)...,B) for (H,B) in tensor_pairs]

    HOFASM_contraction!(preprocessed_pairs,X_k,X_k_1)
    S_opt = dot(X_k,X_k_1)
    x_opt = build_assignment(rand(n,m),return_matrix=true)

    while true

        #Graduated assignment routine
        max_elem = maximum(X_k_1)
        inflation_normalization = x -> exp(x*beta/max_elem)
        broadcast!(inflation_normalization,B_k_1,X_k_1)

        while true

            X_normalized = copy(B_k_1)
            Bistochastic_Normalization!(X_normalized,15)

            if norm(X_normalized - B_k_1) < 1e-6
                B_k_1 = X_normalized
                break
            else
                B_k_1 = copy(X_normalized)
            end
        end

        B_k_1 = build_assignment(B_k_1,return_matrix=true)

        diff = B_k_1-X_k_1
        println(norm(diff))
        C = dot(X_k_1,diff)
    #    println(C)

        HOFASM_contraction!(preprocessed_pairs,diff,temp2)
        D = dot(diff,temp2)
#        println(D)

        if D >= 0
            X_k = copy(B_k_1)
        else
            r = minimum((-C/D,1))
            X_k .+= r * diff
        end

        temp3 = copy(B_k_1)
        HOFASM_contraction!(preprocessed_pairs,B_k_1,temp3)
        S = dot(temp3,B_k_1)
        println(S)

        if S >= S_opt
            x_opt = copy(B_k_1)
            S_opt = S
        end

  #      println(X_k)

        i += 1
        if i >= max_iterations
            break
        else
            #tensor vector contraction
            HOFASM_contraction!(preprocessed_pairs,X_k,X_k_1)
        end
    end

    return x_opt
end

#=
#TODO: could generalize to take a contraction function to unify all graduated assignment routines
function HOM_graduated_assignment(A::ssten.COOTen,max_iterations::Int=12,tol::Float64=1e-4)

    #    X_k = np.ones((n,m))
    n = A.cubical_dimension
    m = Int(sqrt(n))
    x_k::Array{Float64,1} = rand(Float64,n)
    x_k ./= norm(x_k)
    beta = 30 # constant from the paper
    i = 1

    while true

        x_k_1 = ssten.contract_k_1(A,x_k)

        max_elem = maximum(x_k_1)
        inflation_normalization = x -> exp(x*beta/max_elem)
        broadcast!(inflation_normalization,x_k_1,x_k_1)

        while true

            x_normalized = copy(x_k_1)
            Bistochastic_Normalization!(x_normalized,10,m,m)
  #          print(f"{norm(X_normalized - X_k_1,'fro')}")
            println(norm(x_normalized - x_k_1))
            if norm(x_normalized - x_k_1) < 1e-4
                x_k_1 = x_normalized
                break
            else
                x_k_1 = copy(x_normalized)
            end
        end
#        return X_k_1

        x_k_1 ./= norm(x_k_1)
        res = norm(x_k_1 - x_k) / norm(x_k)

        #print("Finished iteration $i with difference $res")
        i += 1
        if i > max_iterations
            return x_k_1
        end

        if res < 1e-4
            println("tolerance met")
            return x_k_1
        else
            x_k = copy(x_k_1)
        end
    end

end
=#


text:
using Test
using CUDA: CUDA, devices, device, functional
using Knet.KnetArrays: KnetMems, KnetPtr, blocksize, initKnetMems, cuallocator
_cuallocator = cuallocator[]

# gc does not work reliably inside function, module, let, if, @testset etc. so some of this code is outside.

function _testkptr(kptrs, navail, nfree)
    # sump = sum(p->p.nptr, values(mem.pools))
    # @show (mem.limit, mem.bytes, mem.avail, length(mem.pools), sump)
    # klens = kptrs == nothing ? [] : (p->p.len).(kptrs)
    # @show klens
    mem = KnetMems[device().handle+1]
    blocksizes = blocksize.(2 .^ (1:10))
    @test length(KnetMems) == length(devices())
    @test length(mem.pools) == 10
    @test sort(collect(keys(mem.pools))) == blocksizes
    @test mem.limit >= mem.bytes
    @test mem.bytes == sum(blocksizes)
    @test mem.bfree == navail
    @test all(Bool[v.nptr==1 && length(v.free)==nfree for (k,v) in mem.pools])
    if nfree == 0
        @test (p->p.len).(kptrs) == blocksizes
    end
end

# Test the knet allocator, this can only be done before first alloc when KnetMems is nothing
_testingkptr = CUDA.functional() && (KnetMems === nothing)
if _testingkptr
    cuallocator[]=false
    initKnetMems()
    @testset "kptr:alloc"   begin; _testkptr(KnetPtr.(2 .^ (1:10)), 0, 0); end
end

GC.gc(true)

if _testingkptr
    @testset "kptr:gc"      begin; _testkptr(nothing, sum(blocksize.(2 .^ (1:10))), 1); end
    @testset "kptr:realloc" begin; _testkptr(KnetPtr.(2 .^ (1:10)), 0, 0); end
end

# Test the cuda allocator.
if CUDA.functional()
    cuallocator[]=true
    usedmem = (isdefined(CUDA,:pool) ? CUDA.pool[].used_memory : CUDA.used_memory)
    used = usedmem()
    @testset "kptr:cuda" begin
        @test (p = KnetPtr(128); usedmem() == used + 128)
        p = nothing
    end
end

GC.gc(true)

if CUDA.functional()
    @testset "kptr:cudagc" begin
        @test usedmem() == used
    end
end

cuallocator[]=_cuallocator

nothing


text:
# This implementation of DQN on cartpole is to verify the cartpole.jl env
using Flux, Gym, Printf
using Flux.Tracker: track, @grad, data#, gradient
using Flux.Optimise: Optimiser, _update_params!#, update!
using Statistics: mean
using DataStructures: CircularBuffer
using CuArrays

#Load game environment
env = CartPoleEnv()
reset!(env)

#ctx = Ctx(env)

#display(ctx.s)
#using Blink# when not on Juno
#body!(Blink.Window(), ctx.s)

# ----------------------------- Parameters -------------------------------------

STATE_SIZE = length(env.state)
ACTION_SIZE = length(env.action_space)
MAX_TRAIN_REWARD = env.x_threshold*env.θ_threshold_radians
SEQ_LEN = 8

# Optimiser params
η = 3f-2

# ------------------------------ Model Architecture ----------------------------

sign(x::TrackedArray) = track(sign, x)
@grad sign(x) = Base.sign.(data(x)), x̄ -> (x̄,)

model = Chain(Dense(STATE_SIZE, 24, relu),
              Dense(24, 48, relu),
              Dense(48, 1, tanh), x->sign(x)) |> gpu

opt = ADAM(η)

action(state) = model(state)

function loss(rewards)
  ep_len = size(rewards, 1)
  max_rewards = ones(Float32, ep_len) * MAX_TRAIN_REWARD |> gpu
  Flux.mse(rewards, max_rewards)
end

# ----------------------------- Helper Functions -------------------------------

function train_reward()
  state = env.state
  x, ẋ, θ, θ̇  = state[1:1], state[2:2], state[3:3], state[4:4]
  # Custom reward for training
  # Product of Triangular function over x-axis and θ-axis
  # Min reward = 0, Max reward = env.x_threshold * env.θ_threshold_radians
  x_upper = env.x_threshold .- x
  x_lower = env.x_threshold .+ x

  r_x     = max.(0f0, min.(x_upper, x_lower))

  θ_upper = env.θ_threshold_radians .- θ
  θ_lower = env.θ_threshold_radians .+ θ

  r_θ     = max.(0f0, min.(θ_upper, θ_lower))

  return r_x .* r_θ
end

function replay(rewards)
  #grads = gradient(() -> loss(rewards), params(model))
  #for p in params(model)
  #  update!(opt, p, grads[p])
  #end

  Flux.back!(loss(rewards))
  _update_params!(opt, params(model))
end

#@which _update_params!(opt, params(model))

function episode!(env, train=true)
  done = false
  total_reward = 0f0
  rewards = []
  frames = 1
  while !done && frames <= 200
    #render(env, ctx)
    #sleep(0.01)

    a = action(env.state)
    s′, r, done, _ = step!(env, a)
    total_reward += r

    train && push!(rewards, train_reward())

    if train && (frames % SEQ_LEN == 0 || done)
      rewards = vcat(rewards...)
      replay(rewards)
      rewards = []
      env.state = param(env.state.data)
    end

    frames += 1
  end
  total_reward
end

# -------------------------------- Testing -------------------------------------

function test()
  score_mean = 0f0
  for _=1:100
    reset!(env)
    total_reward = episode!(env, false)
    score_mean += total_reward / 100
  end
  return score_mean
end

# ------------------------------ Training --------------------------------------

e = 1

while true
  global e
  reset!(env)
  total_reward = episode!(env)
  print("Episode: $e | Score: $total_reward | ")

  score_mean = test()
  score_mean_str = @sprintf "%6.2f" score_mean
  print("Mean score over 100 test episodes: " * score_mean_str)

  println()

  if score_mean > 195
    println("CartPole-v0 solved!")
    break
  end
  e += 1
end


text:
function equivalence_singlevariable_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        CP.Equivalence{MOI.EqualTo{Int}, MOI.EqualTo{Int}},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x2]),
        CP.Equivalence(MOI.EqualTo(1), MOI.EqualTo(1)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
    end
end

function equivalence_scalaraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{Int},
        MOI.EqualTo{Int},
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.Equivalence{MOI.EqualTo{Int}, MOI.EqualTo{Int}},
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, 1 * x1, MOI.EqualTo(1))
    c2 = MOI.add_constraint(
        model,
        MOIU.vectorize([x1, x2]),
        CP.Equivalence(MOI.EqualTo(1), MOI.EqualTo(1)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
    end
end

const equivalencetests = Dict(
    "equivalence_singlevariable" => equivalence_singlevariable_test,
    "equivalence_scalaraffinefunction" =>
        equivalence_scalaraffinefunction_test,
)

MOIT.@moitestset equivalence


text:
@timedtestset "gaussian_process_probabilistic_programme" begin

    @timedtestset "split" begin
        x = BlockData(randn(5), randn(4))
        @testset "Vector" begin
            x1, x2 = split(x, randn(9))
            @test length(x1) == 5
            @test length(x2) == 4
        end
        @testset "Matrix" begin
            x1, x2 = split(x, randn(9, 3))
            @test size(x1) == (5, 3)
            @test size(x2) == (4, 3)
        end
    end

    # Build a toy collection of processes.
    gpc = GPC()
    f1 = wrap(GP(sin, SEKernel()), gpc)
    f2 = wrap(GP(cos, Matern52Kernel()), gpc)
    f3 = f1 + 3 * f2

    # Use them to build a programme.
    f = Stheno.GPPP((f1 = f1, f2 = f2, f3 = f3), gpc)

    # The same answers should be obtained manually or via the GPPP.
    @testset "External Consistency" begin

        x0 = GPPPInput(:f1, randn(4))
        x1 = GPPPInput(:f3, randn(3))

        @test mean(f1, x0.x) == mean(f, x0)
        @test mean(f3, x1.x) == mean(f, x1)

        @test cov(f1, x0.x) == cov(f, x0)
        @test cov(f3, x1.x) == cov(f, x1)

        @test cov(f1, f3, x0.x, x1.x) == cov(f, x0, x1)
        @test var(f3, f1, x1.x, x0.x) == var(f, x1, x0)

        y = rand(f(x1))
        @test cov(posterior(f3(x1.x), y)(x1.x)) == cov(posterior(f(x1), y)(x1))
    end

    # The GPPP must be self-consistent like any other AbstractGP.
    # This should hold for all of the various permutations of applicable input types.
    @testset "Internal Conistency ($(typeof(x0)), $(typeof(x1))" for (x0, x1) in [
        (
            GPPPInput(:f1, randn(4)),
            GPPPInput(:f3, randn(3)),
        ),
        (
            GPPPInput(:f1, randn(4)),
            BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))]),
        ),
        (
            BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))]),
            GPPPInput(:f1, randn(4)),
        ),
        (
            BlockData([GPPPInput(:f2, randn(3)), GPPPInput(:f3, randn(2))]),
            BlockData([GPPPInput(:f1, randn(6))]),
        ),
    ]

        atol=1e-9
        rtol=1e-9

        m = mean(f, x0)
        @test m isa AbstractVector{<:Real}
        @test length(m) == length(x0)

        @assert length(x0) ≠ length(x1)

        # Check that unary cov is consistent with binary cov and conforms to the API
        K_x0 = cov(f, x0)
        @test K_x0 isa AbstractMatrix{<:Real}
        @test size(K_x0) == (length(x0), length(x0))
        @test K_x0 ≈ cov(f, x0, x0) atol=atol rtol=rtol
        @test minimum(eigvals(K_x0)) > -atol
        @test K_x0 ≈ K_x0' atol=atol rtol=rtol

        # Check that single-process binary cov is consistent with single-process binary-cov
        K_x0_x1 = cov(f, x0, x1)
        @test K_x0_x1 isa AbstractMatrix{<:Real}
        @test size(K_x0_x1) == (length(x0), length(x1))

        # Check that single-process binary var is consistent.
        K_x0_x0_diag = var(f, x0, x0)
        @test K_x0_x0_diag isa AbstractVector{<:Real}
        @test length(K_x0_x0_diag) == length(x0)
        @test K_x0_x0_diag ≈ diag(cov(f, x0, x0)) atol=atol rtol=rtol

        # Check that unary var conforms to the API and is consistent with unary cov
        K_x0_diag = var(f, x0)
        @test K_x0_diag isa AbstractVector{<:Real}
        @test length(K_x0_diag) == length(x0)
        @test K_x0_diag ≈ diag(cov(f, x0)) atol=atol rtol=rtol
    end

    @testset "gppp macro" begin

        # Declare a GPPP using the helper functionality.
        f = @gppp let
            f1 = GP(SEKernel())
            f2 = GP(Matern52Kernel())
            f3 = f1 + f2
        end
    end

    # No custom rules to worry about, just need to make sure that nothing errors.
    @timedtestset "Zygote" begin
        x = GPPPInput(:f3, randn(5))
        s = 0.1
        y = rand(f(x, s))
        Zygote.gradient((x, y, f, s) -> logpdf(f(x, s), y), x, y, f, s)
    end
end


text:
# Mu' = f
"""
$(TYPEDEF)

Defines a steady state problem.

# Fields

$(FIELDS)
"""
struct SteadyStateProblem{uType,isinplace,P,F,K} <: AbstractSteadyStateProblem{uType,isinplace}
  """f: The function in the ODE."""
  f::F
  """The initial guess for the steady state."""
  u0::uType
  """Parameter values for the ODE function."""
  p::P
  kwargs::K
  @add_kwonly function SteadyStateProblem{iip}(f::AbstractODEFunction{iip},
                                               u0,p=NullParameters();
                                               kwargs...) where {iip}
    new{typeof(u0),isinplace(f),typeof(p),typeof(f),typeof(kwargs)}(f,u0,p,kwargs)
  end

  """
  $(SIGNATURES)

  Define a steady state problem using the given function.
  `isinplace` optionally sets whether the function is inplace or not.
  This is determined automatically, but not inferred.
  """
  function SteadyStateProblem{iip}(f,u0,p=NullParameters()) where iip
    SteadyStateProblem(ODEFunction{iip}(f),u0,p)
  end
end

"""
$(SIGNATURES)

Define a steady state problem using an instance of
[`AbstractODEFunction`](@ref AbstractODEFunction).
"""
function SteadyStateProblem(f::AbstractODEFunction,u0,p=NullParameters();kwargs...)
  SteadyStateProblem{isinplace(f)}(f,u0,p;kwargs...)
end

function SteadyStateProblem(f,u0,p=NullParameters();kwargs...)
  SteadyStateProblem(ODEFunction(f),u0,p;kwargs...)
end

"""
$(SIGNATURES)

Define a steady state problem from a standard ODE problem.
"""
SteadyStateProblem(prob::AbstractODEProblem) =
      SteadyStateProblem{isinplace(prob)}(prob.f,prob.u0,prob.p)


text:
module YaoQDNN

using Yao

include("layer.jl")

include("AD.jl")

end # module