text:
@info "Testing whiten.jl"

@testset "whiten_data_chol!(X; dims, df)" begin
    n = 10
    p = 3
    for T in (Float32, Float64)
        # test matrix with too few rows
        @test_throws ErrorException DA.whiten_data_chol!(zeros(T,p,p), dims=1)
        @test_throws ErrorException DA.whiten_data_chol!(zeros(T,p,p), dims=2)

        # test singular matrix
        @test_throws ErrorException DA.whiten_data_chol!(zeros(T,n,p), dims=1)
        @test_throws ErrorException DA.whiten_data_chol!(zeros(T,p,n), dims=2)

        # test whitening
        X = T[diagm(0 => ones(T, p));
              rand(n-p, p) .- 0.5]
        X .= X .- mean(X, dims=1)  # Data must be centered
        Xt = copy(transpose(X))

        ### rows
        X_test = copy(X)
        W, detΣ = DA.whiten_data_chol!(X_test, dims=1)
        @test isapprox(cov(X*W, dims=1), diagm(0 => ones(T, p)))
        @test isapprox(det(cov(X, dims=1)), detΣ)

        ### cols
        Xt_test = copy(Xt)
        W, detΣ = DA.whiten_data_chol!(Xt_test, dims=2)
        @test isapprox(cov(W*Xt, dims=2), diagm(0 => ones(T, p)))
        @test isapprox(det(cov(Xt, dims=2)), detΣ)
    end
end

@testset "whiten_data_svd!(X, γ; dims, df)" begin
    n = 10
    p = 3
    for T in (Float32, Float64)
        # test limits for γ
        @test_throws ErrorException DA.whiten_data_svd!(zeros(T,p,p), T(0) - eps(T(0)), dims=1)
        @test_throws ErrorException DA.whiten_data_svd!(zeros(T,p,p), T(1) - eps(T(1)), dims=1)

        # test matrix with too few rows
        @test_throws ErrorException DA.whiten_data_svd!(zeros(T,p,p), T(0), dims=1)
        @test_throws ErrorException DA.whiten_data_svd!(zeros(T,p,p), T(0), dims=2)

        # test singular matrix
        @test_throws ErrorException DA.whiten_data_svd!(zeros(T,n,p), T(0), dims=1)
        @test_throws ErrorException DA.whiten_data_svd!(zeros(T,p,n), T(0), dims=2)

        # test whitening
        Iₚ = diagm(0 => ones(T, p))

        X = T[Iₚ; rand(n-p, p) .- 0.5]
        X .= X .- mean(X, dims=1)  # Data must be centered
        Xt = copy(transpose(X))

        S = X'X
        Σ = S/(n-1)

        for γ in range(zero(T), step=convert(T, 0.25), stop=one(T))
            S_γ = (one(T) - γ)*S + γ*(tr(S)/p)*I
            Σ_γ = (one(T) - γ)*Σ + γ*(tr(Σ)/p)*I

            ### rows
            W_test, detΣ_test = DA.whiten_data_svd!(copy(X), γ, dims=1)

            @test isapprox(det(Σ_γ), detΣ_test)
            @test isapprox(transpose(W_test)*Σ_γ*W_test, Iₚ)
            if γ == zero(T)
                @test isapprox(cov(X*W_test, dims=1), Iₚ)
            end

            W_test, detS_test = DA.whiten_data_svd!(copy(X), γ, dims=1, df=1)
            @test isapprox(det(S_γ), detS_test)
            @test isapprox(transpose(W_test)*S_γ*W_test, Iₚ)

            ### cols
            W_test, detΣ_test = DA.whiten_data_svd!(copy(Xt), γ, dims=2)

            @test isapprox(det(Σ_γ), detΣ_test)
            @test isapprox(W_test*Σ_γ*transpose(W_test), Iₚ)
            if γ == zero(T)
                @test isapprox(cov(W_test*Xt, dims=2), Iₚ)
            end

            W_test, detS_test = DA.whiten_data_svd!(copy(Xt), γ, dims=2, df=1)
            @test isapprox(det(S_γ), detS_test)
            @test isapprox(W_test*S_γ*transpose(W_test), Iₚ)
        end
    end
end


text:
# ============================================
# Hybrid flowpipe of possibly different types
# ============================================

"""
    PartitionedFlowpipe{N, D, FT<:AbstractFlowpipe, VOA<:VectorOfArray{N, D, Vector{FT}}} <: AbstractFlowpipe

Type that wraps a vector of flowpipes of possibly different types.

### Fields

- `Fk`  -- vector of flowpipes
- `ext` -- (optional, default: empty) dictionary for extensions

### Notes
"""
struct PartitionedFlowpipe{T, S<:Tuple} <: AbstractFlowpipe # TODO: ask <:AbstractFlowpipe for each element in the tuple..?
    Fk::ArrayPartition{T, S}
    ext::Dict{Symbol, Any}
end

ext: jl


text:
using Images

urlbase = "http://www.imagemagick.org/Usage/images/"

const savedir = joinpath(tempdir(), "Images")
const writedir = joinpath(savedir, "write")

function getfile(name)
    file = joinpath(savedir, name)
    if !isfile(file)
        file = download(urlbase*name, file)
    end
    file
end

if !isdir(savedir)
    mkdir(savedir)
end
if !isdir(writedir)
    mkdir(writedir)
end

# Gray
file = getfile("jigsaw_tmpl.png")
img = imread(file)
@assert colorspace(img) == "Gray"
@assert ndims(img) == 2
@assert colordim(img) == 0
@assert eltype(img) == Uint8
outname = joinpath(writedir, "jigsaw_tmpl.png")
imwrite(img, outname)
imgc = imread(outname)
@assert img.data == imgc.data

# Gray with alpha channel
file = getfile("wmark_image.png")
img = imread(file)
@assert colorspace(img) == "GrayAlpha"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert eltype(img) == Uint8
outname = joinpath(writedir, "wmark_image.png")
imwrite(img, outname)
sleep(0.2)
imgc = imread(outname)
# @assert img.data == imgc.data   # libmagick bug: doesn't write GrayAlpha properly?

# RGB
file = getfile("rose.png")
img = imread(file)
@assert colorspace(img) == "RGB"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert size(img, 1) == 3
@assert eltype(img) == Uint8
outname = joinpath(writedir, "rose.ppm")
imwrite(img, outname)
imgc = imread(outname)
T = eltype(imgc)
lim = limits(imgc)
@assert (typeof(lim[1]) == typeof(lim[2]) == T)  # issue #62
@assert img.data == imgc.data

# RGBA with 16 bit depth
file = getfile("autumn_leaves.png")
img = imread(file)
@assert colorspace(img) == "ARGB"
@assert ndims(img) == 3
@assert colordim(img) == 1
@assert size(img, 1) == 4
@assert eltype(img) == Uint16
outname = joinpath(writedir, "autumn_leaves.png")
imwrite(img, outname)
sleep(0.2)
imgc = imread(outname)
@assert img.data == imgc.data

# Indexed
file = getfile("present.gif")
img = imread(file)


text:
function restart_25(backcolor)
    empty_figure(filename = joinpath(@__DIR__, "test_25.png"))
    background(backcolor)
    sethue(PALETTE[8])
end
"Prepare Gaussian noise mask as a first iteration. Normalize it to 0..1"
function noisepic_1(A)
    no = randn(size(A));
    # That is too random (?). Let's try a more gently varying nose on the left 2/3
    rate = 0.1
    nrows = size(no)[1]
    ncols = div(size(no)[2], 3)
    for row in 1:nrows
        for col in 1:ncols
            no[row,col] = noise(row * rate * 2, col * rate * 2)
        end
        for col in (ncols+1):(2 * ncols)
            no[row,col] = noise(row * rate * 4, col * rate * 4)
        end
    end
    normalize_datarange(no)
end
"""
This is not for use in the final algorithm, just for debugging.
Plots a polyline until NaN or end of vectors.
"""
function prline_2(origo, vx, vy)
    move(origo + (vx[1], vy[1]))
    for (x, y) in zip(vx, vy)
        isnan(x) && break
        isnan(y) && break
        line(origo + (x, y))
    end
    do_action(:stroke)
end

"""
    rk4_step_1!(f, vx, vy, h, n)

Coordinates (vx[n + 1], vy[n + 1]) are updated with the estimated position,
along function f.
        f is a tuple-valued gradient in two dimensions
        vx and vy are vectors representing coordinates like (vx[n], vy[n])
        h is a step-size quantity. If the integration variable is time,
        it defines the length of time betweeen points. It can be negative.
"""
function rk4_step_1!(f, vx, vy, h, n)
    fx0, fy0 = f(vx[n], vy[n])
    x1 = vx[n] + fx0 * h * 0.5
    y1 = vy[n] + fy0 * h * 0.5
    fx1, fy1 = f(x1, y1)
    x2 = vx[n] + fx1 * h * 0.5
    y2 = vy[n] + fy1 * h * 0.5
    fx2, fy2  = f(x2, y2)
    x3 = vx[n] + fx2 * h
    y3 = vy[n] + fy2 * h
    fx3, fy3 = f(x3, y3)
    vx[n + 1] = vx[n] + 1/6 * h * ( fx0  + 2∙fx1 + 2∙fx2 + fx3 )
    vy[n + 1] = vy[n] + 1/6 * h * ( fy0  + 2∙fy1 + 2∙fy2 + fy3 )
end

"""
    rk4_steps_1!(f, vx, vy, h)

Coordinates (vx[n > 1], vy[n > 1]) are updated with the estimated positions,
along function f.
     f(x,y) returns a tuple-valued gradient in two dimensions
     vx and vy are vectors representing coordinates like (vx[n], vy[n])
     h is a step-size quantity. If the integration variable is time,
     it defines the length of time betweeen points. It can be negative
"""
function rk4_steps_1!(f, vx, vy, h)
    @assert length(vx) == length(vy)
    if !isnan(f(vx[1], vy[1])[1])
        for n in 1:(length(vx) - 1)
            rk4_step_1!(f, vx, vy, h, n)
        end
    else
        fill!(vx, NaN * vx[1])
        fill!(vy, NaN * vy[1])
    end
    nothing
end

function _convolve_contrib_1(prev, w, nxy, vx, vy, i)
    @assert i > 1
    x = vx[i]
    y = vy[i]
    x0 = vx[i - 1]
    y0 = vy[i - 1]
    if isnan(x0) || isnan(y0) || isnan(x) || isnan(y)
        return prev
    end
    k = w[i]
    ds =  hypot((x - x0, y - y0))
    n = nxy(x, y)
    n * k * ds
end

"""
    convolute_pixel_1(wf, wb, vxf, vyf, vxb, vyb, nxy, h, cutoff)
where
    wf    Window forward: A vector defining a convolution (filter) window, or kernel, with the same length as buffers.
          The second value is used for weighting the second coordinate.
    wb    Window backward: Similar, but for the backwards direction. The second value is used for weighting the second coordinate.
          Hence, it will be identical to wf in many cases.
    vxf   Vector or buffer of x coordinates, forward direction
    vyf   Vector y coordinates, forward.
    vxb   Vector x, backward
    vyb   Vector y, backwards direction
    fxy   Function, vector field (output is a tuple)
    nxy   Function, noise field
    h     Step length, duration
    cutoff Velocity, substituted to calculate out-of bounds streamlines
"""
function convolute_pixel_1(wf, wb, vxf, vyf, vxb, vyb, nxy, h, cutoff)
    @assert length(vxf) == length(vyf)
    @assert length(vxb) == length(vyb)
    x0 = vxf[1]
    y0 = vyf[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    pv = 0.0 * prevcontrib
    for i = 2:length(vxf)
        prevcontrib = _convolve_contrib_1(prevcontrib, wf, nxy, vxf, vyf, i)
        pv += prevcontrib
    end
    x0 = vxb[1]
    y0 = vyb[1]
    prevcontrib = nxy(x0, y0) * cutoff * h
    if isnan(x0) || isnan(y0)
        return 0.0 * prevcontrib
    end
    for i = 2:length(vxb)
        prevcontrib = _convolve_contrib_1(prevcontrib, wb, nxy, vxb, vyb, i)
        pv += prevcontrib
    end
    pv
end

function convolute_image_1(xs, ys, f_xy, n_xy, h, cutoff)
    # Output buffer
    M = Array{Float64}(undef, length(ys), length(xs)) # Image processing convention: a column has a horizontal line of pixels
    # Window, triangle. Experiment with FIR, Hanning
    wf = SVector{NS_1}(range(1.0, 0.1, length = NS_1));
    wb = SVector{NS_1}(range(1.0, 0.1, length = NS_1));
    cx = xs[1]::Quantity{Float64}
    cy = ys[1]
    @assert cx isa Quantity{Float64}
    @assert cy isa Quantity{Float64}
    @assert n_xy(cx, cy) isa Float64
    @assert first(ys) < last(ys)
    # prepare a mutable but fast working buffer forward, and another backward.
    vxf = similar(SVector{NS_1}(fill(cx, NS_1)))
    vyf = similar(vxf)
    vxb = similar(vxf)
    vyb = similar(vxf)
    @assert wf[1] == 1.0    # Since hardcoded below (??)
    @assert wb[1] == 1.0
    i = 1
    j = 1
    rowsy = length(ys)
    for (i::Int64, cx::Quantity{Float64}) in enumerate(xs), (j::Int64, cy::Quantity{Float64}) in enumerate(ys)
        vxf[1], vyf[1], vxb[1], vyb[1] = cx, cy, cx, cy
        # Find the streamline path
        # Put coordinates to sample in forward and aft buffers (no need to store them in memory really?
        rk4_steps_1!(f_xy, vxf, vyf, h)
        rk4_steps_1!(f_xy, vxb, vyb, -h)
        # Find the intensity for our one pixel.
        pv  = convolute_pixel_1(wf, wb, vxf, vyf, vxb, vyb, n_xy, h, cutoff)
        # Find the original indexes and update the image matrix
        ii = rowsy + 1 - j
        jj = i
        M[ii, jj] = pv / oneunit(pv)
    end
    normalize_datarange(M)
end

"Draws a few random streamlines. xs and ys are iterators covering 'all' coordinates for function f_xy"
function draw_streamlines_1(origo, xs, ys, f_xy, h)
    cx = xs[1]::Quantity{Float64}
    cy = ys[1]
    @assert cx isa Quantity{Float64}
    @assert cy isa Quantity{Float64}
    # prepare a mutable but fast working buffer forward, and another backward.
    vxf = similar(SVector{NS_1}(fill(cx, NS_1)))
    vyf = similar(vxf)
    vxb = similar(vxf)
    vyb = similar(vxf)
    for (i::Int64, cx::Quantity{Float64}) in enumerate(xs), (j::Int64, cy::Quantity{Float64}) in enumerate(ys)
        vxf[1], vyf[1], vxb[1], vyb[1] = cx, cy, cx, cy
        if rand() < 0.0002
            # Find the streamline path
            rk4_steps_1!(f_xy, vxf, vyf, h)
            rk4_steps_1!(f_xy, vxb, vyb, -h)
            # Blueish is forward
            sethue(PALETTE[1])
            prline_2(origo, vxf, vyf)
            # Reddish and darker is backward
            sethue(color_with_lumin(PALETTE[4], 30))
            prline_2(origo, vxb, vyb)
            sethue(PALETTE[8])
        end
    end
    nothing
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

using Test

@testset "FEMQuad.jl" begin
    @testset "test_glquad" begin include("test_glquad.jl") end
    @testset "test_gltri" begin include("test_gltri.jl") end
    @testset "test_gltet" begin include("test_gltet.jl") end
    @testset "test_glwed" begin include("test_glwed.jl") end
    @testset "test_glpyr" begin include("test_glpyr.jl") end
end


text:
@reexport module Systems

using Distributions
using DynamicalSystems
import DynamicalSystemsBase:
    DynamicalSystem,
    ContinuousDynamicalSystem,
    DiscreteDynamicalSystem,
    Dataset

import DifferentialEquations: @ode_def

include("discretemaps/ar1.jl")
include("discretemaps/anishchenko1.jl")
include("discretemaps/henon2.jl")
include("discretemaps/henon4.jl")
include("discretemaps/henontriple.jl")
include("discretemaps/linearmap.jl")
include("discretemaps/linearmap3d_nonlinearcoupling.jl")
include("discretemaps/nontrivial_pegiun.jl")
include("discretemaps/logistic2.jl")
include("discretemaps/logistic3.jl")
include("discretemaps/logistic4.jl")
include("discretemaps/var1.jl")
include("discretemaps/var1coupled.jl")
include("discretemaps/verdes.jl")

include("continuous_systems/chuacircuits_driven.jl")
include("continuous_systems/chuacircuit_nscroll_sine.jl")
include("continuous_systems/hindmarsh_rose.jl")
include("continuous_systems/lorenztriple.jl")
include("continuous_systems/mediated_link.jl")
include("continuous_systems/rosslerrossler.jl")
include("continuous_systems/rosslerlorenz.jl")

# We export both the equations of motion and function to generate
# the systems for every example. The equations of motion functions
# are always prepended with `eom_`.
# Some systems (involving multiple lags) do not return instances
# of DiscreteDynamicalSystem, but a Dataset with the finished
# iterated map.
export
###############
# Discrete maps
###############
eom_ar1, ar1,
eom_anishchenko1, anishchenko1,
eom_henon2, henon2,
eom_henon_triple, henon_triple,
eom_henon4, henon4,
eom_linearmap1, linearmap1,
eom_linear3d_nonlinearcoupling, linear3d_nonlinearcoupling,
eom_logistic2, logistic2,
eom_logistic3, logistic3,
eom_logistic4, logistic4,
eom_nontrivial_pegiun, nontrivial_pegiun,
eom_var1, var1,
eom_var1coupled, var1coupled,
eom_verdes, verdes

export
####################
# Continuous systems
###################
eom_chuacircuit_nscroll_sine, chuacircuit_nscroll_sine,
eom_chuacircuits_driven, chuacircuits_driven,
eom_hindmarsh_rose, hindmarsh_rose,
eom_lorenz_triple, lorenz_triple,
eom_mediated_link, mediated_link,
eom_rossler_rossler, rossler_rossler,
eom_rossler_lorenz, rossler_lorenz

# Initialise all the systems once, generating a trajectory.
#trajectory(anishchenko1(), 10)
#trajectory(henon2(), 10)
#trajectory(henon4(), 10)
#trajectory(logistic2(), 10)
#trajectory(logistic3(), 10)
#trajectory(logistic4(), 10)
#trajectory(var1(), 10)
#trajectory(var1coupled(), 10)
#trajectory(verdes(), 10)

#trajectory(rossler_rossler(), 10)
#trajectory(rossler_lorenz(), 10)
#trajectory(lorenz_triple(), 10)
#trajectory(mediated_link(), 10)
#trajectory(chuacircuits_driven(), 10)
#trajectory(chuacircuit_nscroll_sine(), 10)

end


text:
export EpisodeTurnBuffer, episode_RTSA_buffer

"""
    EpisodeTurnBuffer{names, types, Tbs} <: AbstractTurnBuffer{names, types}
    EpisodeTurnBuffer{names, types}() where {names, types}

Similar to [`CircularTurnBuffer`](@ref), but instead of using [`CircularArrayBuffer`](@ref), it uses a vector to store each element specified by `names` and `types`. And when it reaches the end of an episode, the buffer is emptied first when a new observation is pushed.

!!! note
    Notice that, before emptifying the `EpisodeTurnBuffer`, the last element of each field is exracted and then pushed at the head of the buffer. Without this step, the first transition of the new episode will be lost!

See also: [`episode_RTSA_buffer`](@ref)
"""
struct EpisodeTurnBuffer{names,types,Tbs} <: AbstractTurnBuffer{names,types}
    buffers::Tbs
    function EpisodeTurnBuffer(configs::Pair{Symbol,DataType}...)
        buffers = merge(NamedTuple(), (name, Vector{type}()) for (name, type) in configs)
        names = Tuple(name for (name, _) in configs)
        types = Tuple{(type for (_, type) in configs)...}
        new{names,types,typeof(buffers)}(buffers)
    end
end

function Base.similar(b::EpisodeTurnBuffer{names,types}) where {names,types}
    EpisodeTurnBuffer(Pair.(names, types)...)
end

"""
    episode_RTSA_buffer(;kwargs...) -> EpisodeTurnBuffer

Initialize an `EpisodeTurnBuffer` with fields of **R**eward, **T**erminal, **S**tate, **A**ction.

# Keywords

- `state_eltype::Type=Int`: the type of state.
- `action_eltype::Type=Int`: the type of action.
- `reward_eltype::Type=Float32`: the type of reward.
- `terminal_eltype::Type=Bool`: the type of terminal.
"""
function episode_RTSA_buffer(
    ;
    state_eltype = Int,
    action_eltype = Int,
    reward_eltype = Float32,
    terminal_eltype = Bool,
)
    EpisodeTurnBuffer(
        :reward => reward_eltype,
        :terminal => terminal_eltype,
        :state => state_eltype,
        :action => action_eltype,
    )
end

"if the last turn is the terminal, then the EpisodeTurnBuffer is full"
isfull(b::EpisodeTurnBuffer{RTSA}) = length(b) > 0 && convert(Bool, b.buffers.terminal[end])

function Base.push!(b::EpisodeTurnBuffer{F}; kw...) where {F}
    if isfull(b)
        for f in F
            bf = getfield(buffers(b), f)
            x = bf[end]
            empty!(bf)
            push!(bf, x)
        end
    end
    for (k, v) in kw
        hasproperty(buffers(b), k) && push!(getproperty(buffers(b), k), v)
    end
end


text:
# evaluating EntropyPenaltyFunctionals

#------------------------------------------------------------------------------
# API functions

"""
  Evaluates [`EntropyPenaltyFunctional`](@ref)s

  Note that this function may overwrite eqn.res.
"""
function _evalFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::EntropyPenaltyFunctional) where {Tmsh, Tsol, Tres}

  val = calcFunctional(mesh, sbp, eqn, opts, func)
  val = MPI.Allreduce(val, MPI.SUM, eqn.comm)

  return val
end

#------------------------------------------------------------------------------
# derivative wrt q

#TODO: EntropyPenaltyFunctional should be EntropyDissipationData?
# derivative of functional wrt q
function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol, Tres}, opts,
                           func::EntropyPenaltyFunctional,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol, Tres}

  # this is a trick to populate eqn.res (the forward sweep of reverse mode)
  calcFunctional(mesh, sbp, eqn, opts, func)

  # compute reverse mode of the contraction, take val_bar = 1
  w_j = zeros(Tsol, mesh.numDofPerNode)
  w_bar_j = zeros(Tres, mesh.numDofPerNode)
  fill!(eqn.res_bar, 0); fill!(eqn.q_bar, 0)
  A0inv = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = sview(eqn.q, :, j, i)
      q_bar_j = sview(eqn.q_bar, :, j, i)
      fill!(w_bar_j, 0)
      convertToIR(eqn.params, q_j, w_j)
      for k=1:mesh.numDofPerNode
        #val += w_j[k]*eqn.res[k, j, i]
        #-----------------------------
        # reverse sweep
        w_bar_j[k] += eqn.res[k, j, i]
        eqn.res_bar[k, j, i] += w_j[k]
      end

      getIRA0inv(eqn.params, q_j, A0inv)

      smallmatTvec_kernel!(A0inv, w_bar_j, q_bar_j, 1, 1)
    end
  end

  if func.do_face_term
    setParallelData(eqn.shared_data_bar, PARALLEL_DATA_ELEMENT)
    startSolutionExchange_rev2(mesh, sbp, eqn, opts, send_q=false)

    # do reverse mode of the face integrals
    if typeof(mesh.sbpface) <: DenseFace
      @assert opts["parallel_data"] == PARALLEL_DATA_ELEMENT

      # local part
      face_integral_functor = func.func
      flux_functor = ErrorFlux_revq()  # not used, but required by the interface
      getFaceElementIntegral_revq(mesh, sbp, eqn, face_integral_functor, flux_functor, mesh.sbpface, mesh.interfaces)

      # define anonymous function for parallel faces
      calc_func = (mesh, sbp, eqn, opts, data, data_bar) ->
        calcSharedFaceElementIntegrals_element_inner_revq(mesh, sbp, eqn, opts,
                              data, data_bar, face_integral_functor, flux_functor)
    else # SparseFace
      flux_functor_revq = func.func_sparseface_revq
      calcFaceIntegral_nopre_revq(mesh, sbp, eqn, opts, flux_functor_revq,
                                  mesh.interfaces)
      # define anonymous functions for parallel faces
      calc_func = (mesh, sbp, eqn, opts, data, data_bar) ->
        calcSharedFaceIntegrals_nopre_element_inner_revq(mesh, sbp, eqn, opts,
          data, data_bar, flux_functor_revq)

    end
  end

  if func.do_lps && opts["addStabilization"]
    addStabilization_revq(mesh, sbp, eqn, opts)
  end

  if func.do_sc && opts["addShockCapturing"]
    evalShockCapturing_revq(mesh, sbp, eqn, opts)
  end

  # do the parallel face calculations
  if func.do_face_term
    finishExchangeData_rev2(mesh, sbp, eqn, opts, eqn.shared_data, eqn.shared_data_bar, calc_func)
  end

  copy!(func_deriv_arr, eqn.q_bar)

  return nothing
end

# wrapper method for NegEntropyDissipationData
function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol, Tres}, opts,
                           func::NegEntropyDissipations,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol, Tres}

  _evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func.func, func_deriv_arr)
  scale!(func_deriv_arr, -1)

  return nothing
end

# method for EntropyDissipaiton2Data
function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol, Tres}, opts,
                           func::EntropyDissipation2Data,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol, Tres}

  _evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func.flux_functional,
                        func_deriv_arr)

  fill!(eqn.q_bar, 0)
  flux_func = func.potential_flux_revq
  val_bar = zeros(Tres, mesh.numDofPerNode); val_bar[1] = 1
  integrateFaceQuantity_revq(mesh, sbp, eqn, opts, flux_func, val_bar)

  for i=1:length(func_deriv_arr)
    func_deriv_arr[i] += eqn.q_bar[i]
  end

  return nothing
end

function _evalFunctionalDeriv_q(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol, Tres}, opts,
                           func::TotalEntropyDissipationData,
                           func_deriv_arr::Abstract3DArray) where {Tmsh, Tsol, Tres}

  fill!(eqn.res, 0)
  bcs_orig = setZeroBCs(mesh, sbp, eqn, opts, func)
  evalResidual(mesh, sbp, eqn, opts)

  # compute reverse mode of the contraction, take val_bar = 1
  w_j = zeros(Tsol, mesh.numDofPerNode)
  w_bar_j = zeros(Tres, mesh.numDofPerNode)
  fill!(eqn.res_bar, 0); fill!(eqn.q_bar, 0)
  A0inv = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = sview(eqn.q, :, j, i)
      q_bar_j = sview(eqn.q_bar, :, j, i)
      fill!(w_bar_j, 0)
      convertToIR(eqn.params, q_j, w_j)
      for k=1:mesh.numDofPerNode
        #val += w_j[k]*eqn.res[k, j, i]
        #-----------------------------
        # reverse sweep
        w_bar_j[k] += eqn.res[k, j, i]
        eqn.res_bar[k, j, i] += w_j[k]
      end

      getIRA0inv(eqn.params, q_j, A0inv)

      smallmatTvec_kernel!(A0inv, w_bar_j, q_bar_j, 1, 1)
    end
  end

  # back propigate res_bar -> q_bar
  # This partially duplicates evalResidual_revq, but that function
  # produces a vector, whereas this function requires a 3D array as output.
  # That was probably a bad design decision, but we are stuck with it now
  setParallelData(eqn.shared_data_bar, PARALLEL_DATA_ELEMENT)
  startSolutionExchange_rev2(mesh, sbp, eqn, opts, send_q=false)

  evalResidual_revq(mesh, sbp, eqn, opts, 0.0)
  assertReceivesWaited(eqn.shared_data_bar)

  copy!(func_deriv_arr, eqn.q_bar)

  resetBCs(mesh, sbp, eqn, opts, bcs_orig)

  return nothing
end

#------------------------------------------------------------------------------
# derivative wrt metrics

function _evalFunctionalDeriv_m(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::AbstractSolutionData{Tsol}, opts,
                           func::EntropyPenaltyFunctional,
                           val_bar::Number=1,
                           ) where {Tmsh, Tsol}

  # compute reverse mode of the contraction, take val_bar = 1
  w_j = zeros(Tsol, mesh.numDofPerNode)
  fill!(eqn.res_bar, 0);

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = sview(eqn.q, :, j, i)
      convertToIR(eqn.params, q_j, w_j)
      for k=1:mesh.numDofPerNode
        #val += w_j[k]*eqn.res[k, j, i]
        #-----------------------------
        # reverse sweep
        eqn.res_bar[k, j, i] += w_j[k]*val_bar
      end
    end
  end

  # do reverse mode of the face integrals
  if func.do_face_term
    if typeof(mesh.sbpface) <: DenseFace
      @assert opts["parallel_data"] == PARALLEL_DATA_ELEMENT

      face_integral_functor = func.func
      flux_functor = ErrorFlux()  # not used, but required by the interface

      # local part
      getFaceElementIntegral_revm(mesh, sbp, eqn, face_integral_functor, flux_functor, mesh.sbpface, mesh.interfaces)

      # define anonymous function for shared faces
      calc_func = (mesh, sbp, eqn, opts, data) ->
        calcSharedFaceElementIntegrals_element_inner_revm(mesh, sbp, eqn, opts,
                              data, face_integral_functor, flux_functor)
    else # SparseFace
      flux_functor_revm = func.func_sparseface_revm
      calcFaceIntegral_nopre_revm(mesh, sbp, eqn, opts, flux_functor_revm,
                                  mesh.interfaces)

      # anonymous function for shared faces
      calc_func = (mesh, sbp, eqn, opts, data) ->
        calcSharedFaceIntegrals_nopre_element_inner_revm(mesh, sbp, eqn, opts,
          data, flux_functor_revm)
    end
  end

  if func.do_lps && opts["addStabilization"]
    addStabilization_revm(mesh, sbp, eqn, opts)
  end

  if func.do_sc && opts["addShockCapturing"]
    evalShockCapturing_revm(mesh, sbp, eqn, opts)
  end

  # do the parallel part of the computation
  if func.do_face_term
    finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calc_func)
  end

  return nothing
end

# NegEntropyDissipationData method
function _evalFunctionalDeriv_m(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::AbstractSolutionData{Tsol}, opts,
                           func::NegEntropyDissipations,
                           val_bar::Number=1,
                           ) where {Tmsh, Tsol}

  _evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func.func, -val_bar)

  return nothing
end

# method for EntropyDissipaiton2Data
function _evalFunctionalDeriv_m(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::EulerData{Tsol, Tres}, opts,
                           func::EntropyDissipation2Data,
                           val_bar::Number=1) where {Tmsh, Tsol, Tres}

  # back propigate val1 contribution
  _evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func.flux_functional,
                        val_bar)

  # back propigate val2 contribution
  flux_func = func.potential_flux_revm
  vals_bar = zeros(Tres, mesh.numDofPerNode); vals_bar[1] = val_bar
  integrateFaceQuantity_revm(mesh, sbp, eqn, opts, flux_func, vals_bar)

  return nothing
end

function _evalFunctionalDeriv_m(mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator,
                           eqn::AbstractSolutionData{Tsol}, opts,
                           func::TotalEntropyDissipationData,
                           val_bar::Number=1,
                           ) where {Tmsh, Tsol}

  # compute reverse mode of the contraction, take val_bar = 1
  w_j = zeros(Tsol, mesh.numDofPerNode)
  fill!(eqn.res_bar, 0);

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = sview(eqn.q, :, j, i)
      convertToIR(eqn.params, q_j, w_j)
      for k=1:mesh.numDofPerNode
        #val += w_j[k]*eqn.res[k, j, i]
        #-----------------------------
        # reverse sweep
        eqn.res_bar[k, j, i] += w_j[k]*val_bar
      end
    end
  end

  bcs_orig = setZeroBCs(mesh, sbp, eqn, opts, func)
  evalResidual_revm(mesh, sbp, eqn, opts, 0.0)
  resetBCs(mesh, sbp, eqn, opts, bcs_orig)

  return nothing
end

#------------------------------------------------------------------------------
# Implementation for each functional

function calcFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::EntropyPenaltyFunctional) where {Tmsh, Tsol, Tres}

  fill!(eqn.res, 0.0)
  # use functions called from evalResidual to compute the dissipation, then
  # contract it with the entropy variables afterwards

  if func.do_face_term
    if typeof(mesh.sbpface) <: DenseFace
      @assert opts["parallel_data"] == PARALLEL_DATA_ELEMENT

      # local part
      face_integral_functor = func.func
      flux_functor = ErrorFlux()  # not used, but required by the interface
      getFaceElementIntegral(mesh, sbp, eqn, face_integral_functor, flux_functor, mesh.sbpface, mesh.interfaces)

      # parallel part

      # temporarily change the face integral functor
      face_integral_functor_orig = eqn.face_element_integral_func
      eqn.face_element_integral_func = face_integral_functor

      # finish data exchange/do the shared face integrals
      # this works even if the receives have already been waited on
      finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceElementIntegrals_element)

      # restore eqn object to original state
      eqn.face_element_integral_func = face_integral_functor_orig
    else  # SparseFace
      flux_functor = func.func_sparseface
      calcFaceIntegral_nopre(mesh, sbp, eqn, opts, flux_functor, mesh.interfaces)

      # parallel part
      doSparseFaceSharedFaceIntegrals(mesh, sbp, eqn, opts, flux_functor)
    end
  end

  if func.do_lps && opts["addStabilization"]
    addStabilization(mesh, sbp, eqn, opts)
  end

  if func.do_sc && opts["addShockCapturing"]
    evalShockCapturing(mesh, sbp, eqn, opts)
  end

  # compute the contraction
  val = zero(Tres)
  w_j = zeros(Tsol, mesh.numDofPerNode)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = sview(eqn.q, :, j, i)
      convertToIR(eqn.params, q_j, w_j)
      for k=1:mesh.numDofPerNode
        val += w_j[k]*eqn.res[k, j, i]

      end
    end
  end

  return val
end

function calcFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::NegEntropyDissipations) where {Tmsh, Tsol, Tres}

  return -calcFunctional(mesh, sbp, eqn, opts, func.func)
end

function calcFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::EntropyDissipation2Data) where {Tmsh, Tsol, Tres}

  # get -\int (wL - wR)^T * F(uL, uR)  (negative sign because the face term
  # has a negative sign in the weak form.
  val1 = calcFunctional(mesh, sbp, eqn, opts, func.flux_functional)

  # now compute \int (psi_L - psi_R)
  # We can't use the regular face integral machinery for this, because it
  # is conservative: it computes resL += flux, resR -= flux, so the sum is
  # always zero.  Instead we need to integral and sum a quantity at the
  # face nodes
  # this only need to work for diagonal E

  flux_func = func.potential_flux
  val = zeros(Tres, mesh.numDofPerNode)
  integrateFaceQuantity(mesh, sbp, eqn, opts, flux_func, val)

  val2 = val[1]
  return val2 + val1
end

function calcFunctional(mesh::AbstractMesh{Tmsh},
            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres}, opts,
            func::TotalEntropyDissipationData) where {Tmsh, Tsol, Tres}

  fill!(eqn.res, 0)
  # zero out all BCs except those specified
  bcs_orig = setZeroBCs(mesh, sbp, eqn, opts, func)

  evalResidual(mesh, sbp, eqn, opts)

  # reset original BCs
  resetBCs(mesh, sbp, eqn, opts, bcs_orig)

  # compute the contraction
  val = zero(Tres)
  w_j = zeros(Tsol, mesh.numDofPerNode)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = sview(eqn.q, :, j, i)
      convertToIR(eqn.params, q_j, w_j)
      for k=1:mesh.numDofPerNode
        val += w_j[k]*eqn.res[k, j, i]

      end
    end
  end

  return val
end

struct BCSet
  bndry_funcs::Vector{BCType}
  bndry_funcs_revm::Vector{BCType_revm}
  bndry_funcs_revq::Vector{BCType_revq}
end

function BCSet(n::Integer)

  bndry_funcs = Array{BCType}(n)
  bndry_funcs_revm = Array{BCType_revm}(n)
  bndry_funcs_revq = Array{BCType_revq}(n)

  return BCSet(bndry_funcs, bndry_funcs_revm, bndry_funcs_revq)
end

"""
  Sets the boundary flux functors to be ZeroFlux for all boundaries except
  those in func.bcnums.  Also does the revq and revm ones.
"""
function setZeroBCs(mesh, sbp, eqn, opts, func::TotalEntropyDissipationData)

  zfunc = ZeroFluxBC(mesh, eqn)
  zfunc_revq = ZeroFluxBC_revq(mesh, eqn)
  zfunc_revm = ZeroFluxBC_revm(mesh, eqn)

  nbcs = length(mesh.bndry_funcs)
  bcs_orig = BCSet(nbcs)

  for i=1:nbcs
    bcs_orig.bndry_funcs[i] = mesh.bndry_funcs[i]
    bcs_orig.bndry_funcs_revm[i] = mesh.bndry_funcs_revm[i]
    bcs_orig.bndry_funcs_revq[i] = mesh.bndry_funcs_revq[i]
    if !(i in func.bcnums)
      mesh.bndry_funcs[i] = zfunc
      mesh.bndry_funcs_revm[i] = zfunc_revm
      mesh.bndry_funcs_revq[i] = zfunc_revq
    end
  end

  return bcs_orig
end

function resetBCs(mesh, sbp, eqn, opts, bcs_orig::BCSet)

  nbcs = length(mesh.bndry_funcs)
  for i=1:nbcs
    mesh.bndry_funcs[i] = bcs_orig.bndry_funcs[i]
    mesh.bndry_funcs_revm[i] = bcs_orig.bndry_funcs_revm[i]
    mesh.bndry_funcs_revq[i] = bcs_orig.bndry_funcs_revq[i]
  end

  return nothing
end

"""
  Do the shared face integrals with the given flux function
"""
function doSparseFaceSharedFaceIntegrals(mesh, sbp, eqn, opts, flux_functor::FluxType)

  # figure out which parallel function to call
  if opts["parallel_data"] == PARALLEL_DATA_FACE
    pfunc = (mesh, sbp, eqn, opts, data) -> calcSharedFaceIntegrals_nopre_inner(mesh, sbp, eqn, opts, data, flux_functor)
  elseif opts["parallel_data"] == PARALLEL_DATA_ELEMENT
    pfunc = (mesh, sbp, eqn, opts, data) -> calcSharedFaceIntegrals_nopre_element_inner(mesh, sbp, eqn, opts, data, flux_functor)
  else
    error("unregonized parallel data: $(opts["parallel_data"])")
  end

  # do shared face integrals
  finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, pfunc)

  return nothing
end


text:
abstract IT2Operation <: IntervalType2FuzzySet

#############################
# Interection of 2 Membership Functions
immutable IT2IntersectOP <: IT2Operation
    t1::IntervalType2FuzzySet
    t2::IntervalType2FuzzySet
end
# export IT2IntersectOP

import Base: getindex
function getindex(op::IT2IntersectOP, value::Real)
	o1 = op.t1[value]
	o2 = op.t2[value]
	(tnorm(o1[1], o2[1]),tnorm(o1[2], o2[2]))
end
export getindex

# Union of 2 Membership Functions
immutable IT2UnionOP <: IT2Operation
    t1::IntervalType2FuzzySet
    t2::IntervalType2FuzzySet
end
# export IT2UnionOP

import Base: getindex
function getindex(op::IT2UnionOP, value::Real)
	o1 = op.t1[value]
	o2 = op.t2[value]
	(tconorm(o1[1], o2[1]),tconorm(o1[2], o2[2]))
end
export getindex

# Complement of
immutable IT2ComplementOP <: IT2Operation
    t1::IntervalType2FuzzySet
end
import Base: getindex
function getindex(op::IT2ComplementOP, value::Real)
	o = op.t1[value]
	(1.0 - o[2],1.0 - o[1])
end
export getindex
# export IT2ComplementOP

Base.complement(mf::IntervalType2FuzzySet) = IT2ComplementOP(mf)
export complement

Base.intersect(mf::IntervalType2FuzzySet, mfs::IntervalType2FuzzySet...) = intersect(IT2IntersectOP(mf,mfs[1]),mfs[2:end]...)
Base.intersect(mf::IntervalType2FuzzySet, mf1::IntervalType2FuzzySet) = IT2IntersectOP(mf,mf1)
export intersect

Base.union(mf::IntervalType2FuzzySet, mfs::IntervalType2FuzzySet...) = union(IT2UnionOP(mf,mfs[1]),mfs[2:end]...)
Base.union(mf::IntervalType2FuzzySet, mf1::IntervalType2FuzzySet) = IT2UnionOP(mf,mf1)
export union


text:
abstract type AbstractObjective end

Base.length(obj::AbstractObjective) = length(obj.cost)

"""$(TYPEDEF)
Objective: stores stage cost(s) and terminal cost functions
Constructors:
```julia
Objective(cost, N)
Objective(cost, cost_term, N)
Objective(costs::Vector{<:CostFunction}, cost_term)
```
"""
struct Objective <: AbstractObjective
    cost::CostTrajectory
end

function Objective(cost::CostFunction,N::Int)
    Objective([cost for k = 1:N])
end

function Objective(cost::CostFunction,cost_terminal::CostFunction,N::Int)
    Objective([k < N ? cost : cost_terminal for k = 1:N])
end

function Objective(cost::CostTrajectory,cost_terminal::CostFunction)
    Objective([cost...,cost_terminal])
end

import Base.getindex

getindex(obj::Objective,i::Int) = obj.cost[i]

"""```julia
cost(obj::Objective, X::Vector, U::Vector, dt::Vector)
```
Calculate cost over entire state and control trajectories
"""
function cost(obj::Objective, X::AbstractVectorTrajectory, U::AbstractVectorTrajectory, dt::Vector)
    N = length(X)
    J = 0.0
    for k = 1:N-1
        J += stage_cost(obj[k],X[k],U[k],dt[k])
    end
    J += stage_cost(obj[N],X[N])
    return J
end

"$(SIGNATURES) Compute the second order Taylor expansion of the cost for the entire trajectory"
function cost_expansion!(Q::ExpansionTrajectory{T}, obj::Objective,
        X::AbstractVectorTrajectory{T}, U::AbstractVectorTrajectory{T}, dt::Vector{T}) where T
    cost_expansion!(Q,obj.cost,X,U,dt)
end

function cost_expansion!(Q::ExpansionTrajectory{T}, c::CostTrajectory,
        X::AbstractVectorTrajectory{T}, U::AbstractVectorTrajectory{T}, dt::Vector{T}) where T
    N = length(X)
    for k = 1:N-1
        cost_expansion!(Q[k],c[k],X[k],U[k],dt[k])
    end
    cost_expansion!(Q[N],c[N],X[N])
end

function cost_expansion_sqrt!(Q::ExpansionTrajectory{T}, obj::Objective,
        X::AbstractVectorTrajectory{T}, U::AbstractVectorTrajectory{T}, dt::Vector{T}) where T
    cost_expansion_sqrt!(Q,obj.cost,X,U,dt)
end

function cost_expansion_sqrt!(Q::ExpansionTrajectory{T}, c::CostTrajectory,
        X::AbstractVectorTrajectory{T}, U::AbstractVectorTrajectory{T}, dt::Vector{T}) where T
    N = length(X)
    for k = 1:N-1
        cost_expansion!(Q[k],c[k],X[k],U[k],dt[k])

        try
            Q[k].xx .= cholesky(Q[k].xx).U
        catch PosDefException
            error("State stage cost Hessian must be PD for sqrt backward pass")
        end
        try
            Q[k].uu .= cholesky(Q[k].uu).U
        catch PosDefException
            error("Control stage cost Hessian must be PD for sqrt backward pass")
        end
    end
    cost_expansion!(Q[N],c[N],X[N])

    try
        Q[N].xx .= cholesky(Q[N].xx).U
    catch PosDefException
        error("Terminal cost Hessian must be PD for sqrt backward pass")
    end
end

@doc raw"""```julia
LQRObjective(Q, R, Qf, xf, N)
```
Create an objective of the form
`` (x_N - x_f)^T Q_f (x_N - x_f) + \sum_{k=0}^{N-1} (x_k-x_f)^T Q (x_k-x_f) + u_k^T R u_k``
"""
function LQRObjective(Q::AbstractArray, R::AbstractArray, Qf::AbstractArray, xf::AbstractVector,N::Int)
    H = zeros(size(R,1),size(Q,1))
    q = -Q*xf
    r = zeros(size(R,1))
    c = 0.5*xf'*Q*xf
    qf = -Qf*xf
    cf = 0.5*xf'*Qf*xf

    ℓ = QuadraticCost(Q, R, H, q, r, c)
    ℓN = QuadraticCost(Qf, qf, cf)

    Objective([k < N ? ℓ : ℓN for k = 1:N])
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
function assess(gen::AtDynamic{T,U,K}, args, constraints) where {T,U,K}
    (key::K, kernel_args) = args
    if length(get_internal_nodes(constraints)) + length(get_leaf_nodes(constraints)) > 1
        error("Not all constraints were consumed")
    end
    kernel_constraints = get_internal_node(constraints, key)
    subtrace::U = assess(gen.kernel, kernel_args, kernel_constraints)
    sub_call = get_call_record(subtrace)
    call = CallRecord{T}(sub_call.score, sub_call.retval, args)
    trace = AtDynamicTrace{T,U,K}(call, subtrace, key)
    trace
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Revise

using Soss
# varying intercept, varying slope
vivs = @model subj,item  begin
    N = length(subj)
    J = maximum(subj)
    K = maximum(item)

    σu ~ HalfCauchy() |> iid(2)
    σw ~ HalfCauchy() |> iid(2)
    σe ~ HalfCauchy()
    u ~ For(1:2, 1:J) do i,j
            Normal(0.0, σu[i])
        end
    w ~ For(1:2, 1:J) do i,j
            Normal(0.0, σw[i])
        end
    β ~ Cauchy() |> iid(2)
    so ~ Bernoulli() |> iid(N)
    rt ~ For(1:N) do i
        μ = (β[1] + u[1,subj[i]] + w[1,item[i]]
          + (β[2] + u[2,subj[i]] + w[2,item[i]]) * (so[i] ? 1 : -1)
            )
        LogNormal(μ,σe)
    end
end;

N = 10;
so = rand(16);
subj = repeat(1:4, inner=4);
item = repeat(1:4, outer=4);

truth = rand(vivs(so=so,subj=subj, item=item))

dynamicHMC(vivs(so=so,subj=subj, item=item), (rt = truth.rt,so=truth.so))


text:
using LinearAlgebraicRepresentation
using Plasm, SparseArrays, LinearAlgebra
Lar = LinearAlgebraicRepresentation
using Rebugger, Triangle
using JuliaInterpreter

# source data from minimal test example
v = [0.0 1.0 0.0 0.0 0.5 0.25 0.25 0.25 1.25 0.25 0.25; 0.0 0.0 1.0 0.0 0.25 0.5 0.25 0.25 0.25 1.25 0.25; 0.0 0.0 0.0 1.0 0.25 0.25 0.5 0.25 0.25 0.25 1.25]
ev = Array{Int64,1}[[1, 2], [1, 3], [2, 3], [1, 4], [2, 4], [3, 4], [5, 6], [5, 7], [6, 7], [5, 8], [5, 9], [6, 8], [6, 10], [9, 10], [7, 8], [7, 11], [9, 11], [10, 11]]
cscEV = ([1, 2, 4, 1, 3, 5, 2, 3, 6, 4, 5, 6, 7, 8, 10, 11, 7, 9, 12, 13, 8, 9, 15, 16, 10, 12, 15, 11, 14, 17, 13, 14, 18, 16, 17, 18], [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11], Int8[-1, -1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1, 1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1])
fe = Array{Int64,1}[[1, 2, 3], [1, 4, 5], [2, 4, 6], [3, 5, 6, 7, 8, 9], [7, 8, 9], [7, 10, 12], [7, 11, 13, 14], [8, 10, 15], [8, 11, 16, 17], [9, 12, 15], [9, 13, 16, 18], [14, 17, 18]]
cscFE = ([1, 2, 1, 3, 1, 4, 2, 3, 2, 4, 3, 4, 4, 5, 6, 7, 4, 5, 8, 9, 4, 5, 10, 11, 6, 8, 7, 9, 6, 10, 7, 11, 7, 12, 8, 10, 9, 11, 9, 12, 11, 12], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18], Int8[1, 1, -1, 1, 1, 1, -1, -1, 1, -1, 1, 1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, -1, 1, -1, 1, 1])
fv = Array{Int64,1}[[2, 3, 1], [4, 2, 1], [4, 3, 1], [7, 4, 2, 3, 5, 6], [7, 5, 6], [5, 8, 6], [9, 10, 5, 6], [7, 5, 8], [7, 9, 11, 5], [7, 8, 6], [7, 10, 11, 6], [9, 10, 11]]

Plasm.view(Plasm.numbering(0.25)((v,[[[k] for k=1:size(v,2)],ev,fv])))

#V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W::Lar.Points, cop_EW::Lar.ChainOp, cop_FE::Lar.ChainOp)

copEV = SparseArrays.sparse(cscEV...);
copFE = SparseArrays.sparse(cscFE...);
V = convert(Lar.Points, v');

#breakpoint(Triangle.constrained_triangulation)

function minimal_3cycles(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp)

    triangulated_faces = Array{Any, 1}(undef, FE.m)

    function face_angle(e::Int, f::Int)
        if !isassigned(triangulated_faces, f)
            vs_idxs = Array{Int64, 1}()
            edges_idxs = FE[f, :].nzind
            edge_num = length(edges_idxs)
            edges = zeros(Int64, edge_num, 2)

            for (i, ee) in enumerate(edges_idxs)
                edge = EV[ee, :].nzind
                edges[i, :] = edge
                vs_idxs = union(vs_idxs, edge)
            end

            vs = V[vs_idxs, :]

            v1 = normalize(vs[2, :] - vs[1, :])
            v2 = [0 0 0]		# added for debug
            v3 = [0 0 0]
            err = 1e-8
            i = 3
            while -err < norm(v3) < err
                v2 = normalize(vs[i, :] - vs[1, :])
                v3 = cross(v1, v2)
                i = i + 1
            end

            M = reshape([v1; v2; v3], 3, 3)

            vs = vs*M

            triangulated_faces[f] = Triangle.constrained_triangulation(
                Array{Float64,2}(vs), vs_idxs, edges, fill(true, edge_num))

        end
        edge_vs = EV[e, :].nzind

        t = findfirst(x->edge_vs[1] in x && edge_vs[2] in x, triangulated_faces[f])

        v1 = normalize(V[edge_vs[2], :] - V[edge_vs[1], :])

        if abs(v1[1]) > abs(v1[2])
            invlen = 1. / sqrt(v1[1]*v1[1] + v1[3]*v1[3])
            v2 = [-v1[3]*invlen, 0, v1[1]*invlen]
        else
            invlen = 1. / sqrt(v1[2]*v1[2] + v1[3]*v1[3])
            v2 = [0, -v1[3]*invlen, v1[2]*invlen]
        end

        v3 = cross(v1, v2)

        M = reshape([v1; v2; v3], 3, 3)

        triangle = triangulated_faces[f][t]
        third_v = setdiff(triangle, edge_vs)[1]
        vs = V[[edge_vs..., third_v], :]*M

        v = vs[3, :] - vs[1, :]
        angle = atan(v[2], v[3])
        return angle
    end

    #EF = FE'
    EF = convert(Lar.ChainOp, LinearAlgebra.transpose(FE))
x	FC = Lar.Arrangement.minimal_cycles(face_angle, true)(V, EF)

	#FC'
    return -convert(Lar.ChainOp, LinearAlgebra.transpose(FC))
end

# Juno.@enter debug_function()
function debug_function()
	minimal_3cycles(V::Lar.Points, copEV::Lar.ChainOp, copFE::Lar.ChainOp)

end
# minimal_3cycles(V::Lar.Points, copEV::Lar.ChainOp, copFE::Lar.ChainOp)


text:
"""
    get_latents(examinees::Vector{<:AbstractExaminee})

#Description

Returns the latent objects of a vector of examinees as a vector of subtypes of `AbstractLatent`.
"""
function get_latents(examinees::Vector{<:AbstractExaminee})
    map( e -> e.latent, examinees)::Vector{<:AbstractLatent}
end

"""
    get_latents_vals(examinees::Vector{<:AbstractExaminee})

#Description

Returns a matrix with latent values displayed by row.
"""
function get_latents_vals(examinees::Vector{<:AbstractExaminee})
    ret = Vector{Vector{Float64}}(undef, size(examinees, 1))
    max_length = size(get_latents_vals(examinees[1]), 1)
    e_2 = 0
    #compute maximum latent dimension
    for e in examinees
        e_2 += 1
        local latents = get_latents_vals(e)
        if size(latents, 1) == 1
            ret[e_2] = [latents]
        else
            ret[e_2] = latents
        end
        #create vectors
        if size(latents, 1) > max_length
            max_length = size(latents, 1)
            for e_3 = 1:(e_2-1)
                local length_i = size(ret[e_3], 1)
                if length_i < max_length
                    ret[e_3] .= vcat(ret[e_3], zeros(Float64, max_length - length_i))
                end
            end
        end
    end
    #hcat vectors
    return reduce(hcat, ret)::Matrix{Float64}
end

"""
    set_prior!(
        examinees::Vector{<:AbstractExaminee},
        prior::Union{Distributions.DiscreteUnivariateDistribution, Distributions.ContinuousUnivariateDistribution}
    )
"""
function  set_prior!(
    examinees::Vector{<:AbstractExaminee},
    prior::Union{Distributions.DiscreteUnivariateDistribution, Distributions.ContinuousUnivariateDistribution}
)
    map( e -> _set_prior!(e.latent, prior), examinees)
    return nothing
end


text:
export pathclip,pathdifference,pathintersection,pathunion,pathxor
function pathclip(path::Vector{NTuple{2,T}},clip::Vector{NTuple{2,R}},
                  cty::ClipType=ClipTypeIntersection,
                  pft::PolyFillType=PolyFillTypeNonZero,
                  cft::PolyFillType=PolyFillTypeNonZero;
                  pathclosed::Bool=true,
                  clipclosed::Bool=true) where {T,R}
    clipper=Clip()
    pathclip=to_intpoints([path,clip])
    add_path!(clipper,pathclip[1], PolyTypeSubject,pathclosed)
    add_path!(clipper,pathclip[2], PolyTypeClip,   clipclosed)
    (success,result)=execute(clipper,cty,pft,cft)
    #return (success,result)
    out = success ? from_intpoints([path;clip],result) : Vector{Vector{NTuple{2,T}}}()
    return (success,out)
end

function pathclip(path::Vector{Vector{NTuple{2,T}}},clip::Vector{NTuple{2,R}},
                  cty::ClipType=ClipTypeIntersection,
                  pft::PolyFillType=PolyFillTypeNonZero,
                  cft::PolyFillType=PolyFillTypeNonZero;
                  pathclosed::Bool=true,
                  clipclosed::Bool=true) where {T,R}
    clipper=Clip()
    pathclip=to_intpoints([path;[clip]])
    add_paths!(clipper,pathclip[1:end-1], PolyTypeSubject,pathclosed)
    add_path!( clipper,pathclip[end],     PolyTypeClip,   clipclosed)
    (success,result)=execute(clipper,cty,pft,cft)
    #return (success,result)
    out = success ? from_intpoints(vcat([path;[clip]]...),result) : Vector{Vector{NTuple{2,T}}}()
    return (success,out)
end

for (f,v) in zip( (:pathdifference,:pathintersection,:pathunion,:pathxor),
                  (ClipTypeDifference,ClipTypeIntersection,ClipTypeUnion,ClipTypeXor) )
  @eval $f(p1::Vector{NTuple{2,T}},p2::Vector{NTuple{2,R}},o...;k...) where {T,R} = pathclip(p1,p2,$v,o...;k...)
  @eval $f(p1::Vector{Vector{NTuple{2,T}}},p2::Vector{NTuple{2,R}},o...;k...) where {T,R} = pathclip(p1,p2,$v,o...;k...)
end
for f in (:pathdifference,:pathintersection,:pathunion,:pathxor)
  @eval begin
  function $f(p1::Vector{T},p2::Vector{R},o...;k...) where {T<:Vector,R<:Vector}
    @assert eltype(T)<:Number && eltype(R)<:Number "Paths should be constructed of Vector{Vector{Number}} or Vector{NTuple{2,Number}}"
    @assert all(length.(p1).==2)&&all(length.(p2).==2) "Clipper paths are restricted to two dimensions"
    S=promote_type(eltype(T),eltype(R))
    (success,result)=$f(vector2ntuple(2,S,p1),vector2ntuple(2,S,p2),o...;k...)
    out=success ? ntuple2vector(result) : Vector{Vector{Vector{S}}}()
    return (success,out)
  end
  function $f(p1::Vector{Vector{T}},p2::Vector{R},o...;k...) where {T<:Vector,R<:Vector}
    @assert eltype(T)<:Number && eltype(R)<:Number "Paths should be constructed of Vector{Vector{Vector{Number}}} or Vector{Vector{NTuple{2,Number}}}"
    @assert all(length.(vcat(p1...)).==2)&&all(length.(p2).==2) "Clipper paths are restricted to two dimensions"
    S=promote_type(eltype(T),eltype(R))
    (success,result)=$f(vector2ntuple(2,S,p1),vector2ntuple(2,S,p2),o...;k...)
    out=success ? ntuple2vector(result) : Vector{Vector{Vector{S}}}()
    return (success,out)
  end
  end
end

function area(path::Vector{NTuple{2,T}}) where T
  minx,width=min_width(path)
  areac=area(to_intpoints(path)) # the area in rescaled Cint^2
  areac/4/typemax(Cint)/typemax(Cint)*prod(width)
end
area(path::Vector{Vector{T}}) where T=area(vector2ntuple(2,T,path))


text:
using KernelDensity
import Statistics.std

function _defaultBandwidth(vs::Vector{Float64})
  return std(vs)*1.06*length(vs)^(-1/5)
end

function _defaultBandwidth(xs::Vector{Float64}, ys::Vector{Float64})
  @assert length(xs) == length(ys)
  len = length(xs)
  h1 = std(xs)*1.06*len^(-1/6)
  h2 = std(ys)*1.06*len^(-1/6)
  return [h1, h2]
end

## Compute adjusted kernel density estimates using an conservative approximation
## of the quantity proposed by:
## Sköld, M. and Roberts, G.O., 2003. Density estimation for the
## Metropolis–Hastings algorithm. Scandinavian Journal of Statistics, 30(4),
## pp.699-718.

## can call plot(kde(vs, [acceptanceRate]))
function kde(vs::Vector{Float64}, acceptanceRate::Float64, adjust::Float64)
  hDefault = _defaultBandwidth(vs)
  hAdjusted = hDefault * (2/acceptanceRate-1)^(0.2) * adjust
  left = minimum(vs) - 3 * hAdjusted
  right = maximum(vs) + 3 * hAdjusted
  xs = range(left, stop=right, length=512)
  ys = pdf(InterpKDE(KernelDensity.kde(vs; bandwidth = hAdjusted)), xs)
  return xs, ys
end

function kde(vs::Vector{Float64}, acceptanceRate::Float64 = 1.0)
  return kde(vs, acceptanceRate, 1.0)
end

## can call contour(kde(xs, ys, [acceptanceRate]))
function kde(xs::Vector{Float64}, ys::Vector{Float64},
  acceptanceRate::Float64, adjust::Float64)
  hAdjusted = _defaultBandwidth(xs, ys) * (2/acceptanceRate-1)^(1/6) * adjust

  left = minimum(xs) - 3 * hAdjusted[1]
  right = maximum(xs) + 3 * hAdjusted[1]
  bottom = minimum(ys) - 3 * hAdjusted[2]
  top = maximum(ys) + 3 * hAdjusted[2]
  xOut = range(left, stop=right, length=128)
  yOut = range(bottom, stop=top, length=128)
  ikde = InterpKDE(KernelDensity.kde((xs,ys); bandwidth = Tuple(hAdjusted)))
  function de(x, y)
    return pdf(ikde, x, y)
  end
  return xOut, yOut, de
end

function kde(xs::Vector{Float64}, ys::Vector{Float64},
  acceptanceRate::Float64 = 1.0)
  return kde(xs, ys, acceptanceRate, 1.0)
end


text:
# Fields should be self-explanatory. Position is an absolute time (in ticks) within the track.
# Please note that velocity cannot be higher than 127 (0x7F).
type Note
    value::UInt8
    duration::UInt64
    position::UInt64
    channel::UInt8
    velocity::UInt8

    Note(value, duration, position, channel, velocity=0x7F) =
        if channel > 0x7F
            error( "Channel must be less than 128" )
        elseif velocity > 0x7F
            error( "Velocity must be less than 128" )
        else
            new(value, duration, position, channel, velocity)
        end
end

const Bs = 0
const C = 0
const Cs = 1
const Db = 1
const D = 2
const Ds = 3
const Eb = 3
const E = 4
const Fb = 4
const Es = 5
const F = 5
const Fs = 6
const Gb = 6
const G = 7
const Gs = 8
const Ab = 8
const A = 9
const As = 10
const Bb = 10
const B = 11
const Cb = 11

import Base.+, Base.-, Base.==

+(n::Note, i::Integer) = Note(n.value + i, n.duration, n.position, n.channel, n.velocity)
+(i::Integer, n::Note) = n + i

-(n::Note, i::Integer) = Note(n.value - i, n.duration, n.position, n.channel, n.velocity)
-(i::Integer, n::Note) = n - i

==(n1::Note, n2::Note) =
    n1.value == n2.value &&
    n1.duration == n2.duration &&
    n1.position == n2.position &&
    n1.channel == n2.channel &&
    n1.velocity == n2.velocity

export Bs, C, Cs, Db, D, Ds, Eb, E, Fb, Es, F, Fs, Gb, G, Gs, Ab, A, As, Bb, B, Cb


text:
using Distributed
using Dates
using RoME
using RoMEPlotting
using Gadfly
@everywhere using RoME, RoMEPlotting, Gadfly

total_meas = parse(Int, ARGS[1])

# Let's load the Manhattan scenario using the g2o file.
file = (normpath(Base.find_package("RoME"), "../..", "examples", "manhattan_incremental.g2o"))
global instructions = importG2o(file)

# Make sure plots look a bit nicer.
latex_fonts = Theme(major_label_font="CMU Serif", major_label_font_size=16pt,
                    minor_label_font="CMU Serif", minor_label_font_size=14pt,
                    key_title_font="CMU Serif", key_title_font_size=12pt,
                    key_label_font="CMU Serif", key_label_font_size=10pt)
Gadfly.push_theme(latex_fonts)

function solve_batch(total_meas::Integer)
    # Choose where to save the data.
    manhattan_total_meas = total_meas
    data_logpath = ENV["HOME"]*"/Documents/wafr/manhattan-batch-$(total_meas)-$(now())"

    # Create initial factor graph with specified logging path.
    fg = LightDFG{SolverParams}(solverParams=SolverParams(logpath=data_logpath))

    # Add initial variable with a prior measurement to anchor the graph.
    addVariable!(fg, :x0, Pose2)
    initial_pose = MvNormal([0.0; 0.0; 0.0], Matrix(Diagonal([0.1;0.1;0.05].^2)))
    addFactor!(fg, [:x0], PriorPose2(initial_pose))

    # Add all variables and measurements.
    # manhattan_total_meas = 5453
    for i in 1:manhattan_total_meas
        parseG2oInstruction!(fg, instructions[i])
    end

    # Solve the graph, and save a copy of the tree.
    saveDFG(fg, "$(getLogPath(fg))/fg-before-solve")
    getSolverParams(fg).maxincidence =
    tree, smt, hist = solveTree!(fg)
    saveDFG(fg, "$(getLogPath(fg))/fg-after-solve")
    saveTree(tree, "$(getLogPath(fg))/tree$(manhattan_total_meas).jld2")
    drawTree(tree, show=false, filepath="$(getLogPath(fg))/bt.pdf")

    # Just store some quick plots.
    pl1 = drawPoses(fg, spscale=0.6, lbls=false)
    Gadfly.draw(PDF("$(getLogPath(fg))/poses$(manhattan_total_meas).pdf", 20cm, 10cm), pl1)

    # Run the garbage collector.
    GC.gc()
end

solve_batch(total_meas)


text:
using Test
using StaticArrays: SVector
using QuantumLattices.Essentials.FockPackage
using QuantumLattices.Essentials.Spatials: Bond,Point,PID,rcoord,azimuthd
using QuantumLattices.Essentials.DegreesOfFreedom: Table,IDFConfig,OID,Operators,oidtype,otype,script,latexformat
using QuantumLattices.Essentials.Terms: Couplings,@subscript,statistics,abbr
using QuantumLattices.Interfaces: dims,inds,⊗,⋅,expand,permute
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.VectorSpaces: IsMultiIndexable,MultiIndexOrderStyle
using QuantumLattices.Mathematics.AlgebraOverFields: ID
import QuantumLattices.FockPackage: fockcouplingnambus

@testset "FID" begin
    fid=FID(orbital=1,spin=1)
    @test fid|>typeof|>fieldnames==(:orbital,:spin,:nambu)
    @test fid'==FID(1,1,2)
    @test fid''==FID(1,1,1)

    fid=FID(1,1,0)
    @test fid'==fid
end

@testset "Fock" begin
    fock=Fock(norbital=1,nspin=2,nnambu=2)
    @test IsMultiIndexable(Fock)==IsMultiIndexable(true)
    @test MultiIndexOrderStyle(Fock)==MultiIndexOrderStyle('C')
    @test dims(fock)==(1,2,2)
    @test inds(FID(1,1,1),fock)==(1,1,1)
    @test FID((1,1,1),fock)==FID(1,1,1)
    @test fock|>collect==[FID(1,1,1),FID(1,1,2),FID(1,2,1),FID(1,2,2)]
end

@testset "FIndex" begin
    @test FIndex|>fieldnames==(:scope,:site,:orbital,:spin,:nambu)
    @test union(PID{Int},FID)==FIndex{Int}
    @test angle(OID(FIndex(1,1,1,1,1),[0.0,0.0],[1.0,2.0],1),[[1.0,0.0],[0.0,1.0]],[0.1,0.0])≈2pi*0.1
    @test angle(OID(FIndex(1,1,1,1,2),[0.0,0.0],[1.0,2.0],1),[[1.0,0.0],[0.0,1.0]],[0.0,0.2])≈-2pi*0.4
end

@testset "script" begin
    @test script(OID(FIndex('c',1,2,1,1)),Val(:site))==1
    @test script(OID(FIndex('c',1,2,1,1)),Val(:orbital))==2
    @test script(OID(FIndex('c',1,2,3,1)),Val(:spinint))==3
    @test script(OID(FIndex('c',1,2,2,1)),Val(:spinsym))=="↑"
    @test script(OID(FIndex('c',1,2,1,1)),Val(:spinsym))=="↓"
    @test script(OID(FIndex('c',1,2,3,1)),Val(:nambu))==""
    @test script(OID(FIndex('c',1,2,3,2)),Val(:nambu))=="\\dagger"
end

@testset "oidtype" begin
    @test oidtype(FID,Point{2,PID{Int}},Nothing,Val(true))==OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Nothing}
    @test oidtype(FID,Point{2,PID{Int}},Table,Val(true))==OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Int}
    @test oidtype(FID,Point{2,PID{Int}},Nothing,Val(false))==OID{FIndex{Int},Nothing,Nothing,Nothing}
    @test oidtype(FID,Point{2,PID{Int}},Table,Val(false))==OID{FIndex{Int},Nothing,Nothing,Int}
end

@testset "FockOperator" begin
    @test latexformat(FOperator)==foptdefaultlatex

    opt=FOperator(1.0,(FIndex(1,1,1,1,2),FIndex(1,1,1,1,1)))
    @test opt|>isnormalordered
    @test opt|>statistics==opt|>typeof|>statistics=='F'

    opt=FOperator(1.0,(FIndex(1,2,1,1,2),FIndex(1,2,1,1,1),FIndex(1,1,1,2,2),FIndex(1,1,1,2,1)))
    @test opt|>isnormalordered==false
    @test repr(opt)=="c^{\\dagger}_{2,1,↓}c^{}_{2,1,↓}c^{\\dagger}_{1,1,↑}c^{}_{1,1,↑}"

    opt1=FOperator(1.5,(FIndex(1,2,1,1,2),FIndex(1,2,1,1,1)))
    opt2=FOperator(2.0,(FIndex(1,2,1,1,1),FIndex(1,2,1,1,2)))
    @test opt1*opt2==nothing

    opt1=FOperator(1.5,(FIndex(1,2,1,1,2),FIndex(1,2,1,1,1)))
    opt2=FOperator(2.0,(FIndex(1,2,1,1,2),FIndex(1,2,1,1,1)))
    @test opt1*opt2==FOperator(3.0,(FIndex(1,2,1,1,2),FIndex(1,2,1,1,1),FIndex(1,2,1,1,2),FIndex(1,2,1,1,1)))

    @test latexformat(BOperator)==boptdefaultlatex

    opt=BOperator(1.0,(FIndex(1,1,1,1,2),FIndex(1,1,1,1,1)))
    @test opt|>statistics==opt|>typeof|>statistics=='B'
    @test repr(opt)=="b^{\\dagger}_{1,1,↓}b^{}_{1,1,↓}"
end

@testset "permute" begin
    i1,i2=OID(FIndex(1,1,1,1,2)),OID(FIndex(1,1,1,1,1))
    @test permute(FOperator,i1,i2)==(FOperator(1),FOperator(-1,ID(i2,i1)))
    @test permute(FOperator,i2,i1)==(FOperator(1),FOperator(-1,ID(i1,i2)))
    @test permute(BOperator,i1,i2)==(BOperator(1),BOperator(1,ID(i2,i1)))
    @test permute(BOperator,i2,i1)==(BOperator(-1),BOperator(1,ID(i1,i2)))

    i1,i2=OID(FIndex(1,1,1,2,2)),OID(FIndex(1,1,1,1,1))
    @test permute(FOperator,i1,i2)==(FOperator(-1,ID(i2,i1)),)
    @test permute(FOperator,i2,i1)==(FOperator(-1,ID(i1,i2)),)
    @test permute(BOperator,i1,i2)==(BOperator(1,ID(i2,i1)),)
    @test permute(BOperator,i2,i1)==(BOperator(1,ID(i1,i2)),)
end

@testset "FCID" begin
    @test FCID(atom=1,orbital=1,spin=1,nambu=1)|>typeof|>fieldnames==(:center,:atom,:orbital,:spin,:nambu,:obsub,:spsub)
end

@testset "FockCoupling" begin
    @test FockCoupling{2}(1.0)|>string=="FockCoupling{2}(value=1.0)"
    @test FockCoupling{2}(1.0,atoms=(1,1))|>string=="FockCoupling{2}(value=1.0,atoms=(1,1))"
    @test FockCoupling{2}(1.0,atoms=(1,1),spins=(1,2))|>string=="FockCoupling{2}(value=1.0,atoms=(1,1),spins=(1,2))"
    @test FockCoupling{2}(2.0)|>repr=="2.0 {2}"

    fc1=FockCoupling{2}(1.5,atoms=(2,1),spins=(@subscript (x,1)),centers=(1,2))
    fc2=FockCoupling{2}(2.0,atoms=(1,2),orbitals=(@subscript (x,y) with x<y),centers=(1,2))
    @test fc1|>repr=="1.5 sl(2,1)⊗sp(x,1)@(1,2) with (*,$(fc1.id[1].spsub)) && (*,$(fc1.id[2].spsub))"
    @test fc2|>repr=="2.0 sl(1,2)⊗ob(x,y)@(1,2) with ($(fc2.id[1].obsub),*) && ($(fc2.id[2].obsub),*)"
    fc=fc1*fc2
    @test fc|>repr=="3.0 sl(2,1,1,2)⊗ob(*,*,x,y)⊗sp(x,1,*,*)@(1,2,1,2) with (*,$(fc1.id[1].spsub)) && (*,$(fc1.id[2].spsub)) && ($(fc2.id[1].obsub),*) && ($(fc2.id[2].obsub),*)"

    fc1=FockCoupling{2}(1.5,spins=(@subscript (x,1)),centers=(1,2))
    fc2=FockCoupling{2}(2.0,orbitals=(@subscript (x,y) with x<y),centers=(1,2))
    fc=fc1⊗fc2
    @test fc|>repr=="3.0 ob(x,y)⊗sp(x,1)@(1,2) with ($(fc2.id[1].obsub),$(fc1.id[2].spsub)) && ($(fc2.id[1].obsub),$(fc1.id[2].spsub))"

    fc1=FockCoupling{2}(1.5,atoms=(2,1),centers=(1,2))
    fc2=FockCoupling{2}(2.0,atoms=(1,2),centers=(1,2))
    fc=fc1⋅fc2
    @test fc|>repr=="3.0 sl(2,2)@(1,2)"

    ex=expand(FockCoupling{2}(2.0,atoms=(1,1)),PID(1,1),Fock(atom=2,norbital=2,nspin=2,nnambu=2))
    @test collect(ex)==[]

    ex=expand(FockCoupling{2}(2.0,atoms=(1,2),orbitals=(1,2)),(PID(1,1),PID(1,2)),(Fock(atom=1,norbital=2,nspin=2,nnambu=2),Fock(atom=2,norbital=2,nspin=2,nnambu=2)))
    @test IsMultiIndexable(typeof(ex))==IsMultiIndexable(true)
    @test MultiIndexOrderStyle(typeof(ex))==MultiIndexOrderStyle('C')
    @test dims(ex)==(1,2)
    @test collect(ex)==[(2.0,(FIndex(1,1,1,1,2),FIndex(1,2,2,1,1))),
                        (2.0,(FIndex(1,1,1,2,2),FIndex(1,2,2,2,1)))
    ]

    ex=expand(FockCoupling{4}(2.0,centers=(1,1,1,1),spins=(2,2,1,1)),PID(1,1),Fock(atom=1,norbital=2,nspin=2,nnambu=2))
    @test dims(ex)==(2,1)
    @test collect(ex)==[(2.0,(FIndex(1,1,1,2,2),FIndex(1,1,1,2,1),FIndex(1,1,1,1,2),FIndex(1,1,1,1,1))),
                        (2.0,(FIndex(1,1,2,2,2),FIndex(1,1,2,2,1),FIndex(1,1,2,1,2),FIndex(1,1,2,1,1)))
    ]

    ex=expand(FockCoupling{4}(2.0,orbitals=(@subscript (α,α,β,β) with α<β),spins=(2,1,1,2),nambus=(2,2,1,1)),PID(1,1),Fock(atom=1,norbital=3,nspin=2,nnambu=2))
    @test dims(ex)==(3,1)
    @test collect(ex)==[(2.0,(FIndex(1,1,1,2,2),FIndex(1,1,1,1,2),FIndex(1,1,2,1,1),FIndex(1,1,2,2,1))),
                        (2.0,(FIndex(1,1,1,2,2),FIndex(1,1,1,1,2),FIndex(1,1,3,1,1),FIndex(1,1,3,2,1))),
                        (2.0,(FIndex(1,1,2,2,2),FIndex(1,1,2,1,2),FIndex(1,1,3,1,1),FIndex(1,1,3,2,1)))
    ]

    fc1=FockCoupling{2}(1.0,spins=(2,2))
    fc2=FockCoupling{2}(-1.0,spins=(1,1))
    ex=expand(fc1*fc2,PID(1,1),Fock(atom=1,norbital=2,nspin=2,nnambu=2))
    @test dims(ex)==(4,1)
    @test collect(ex)==[(-1.0,(FIndex(1,1,1,2,2),FIndex(1,1,1,2,1),FIndex(1,1,1,1,2),FIndex(1,1,1,1,1))),
                        (-1.0,(FIndex(1,1,1,2,2),FIndex(1,1,1,2,1),FIndex(1,1,2,1,2),FIndex(1,1,2,1,1))),
                        (-1.0,(FIndex(1,1,2,2,2),FIndex(1,1,2,2,1),FIndex(1,1,1,1,2),FIndex(1,1,1,1,1))),
                        (-1.0,(FIndex(1,1,2,2,2),FIndex(1,1,2,2,1),FIndex(1,1,2,1,2),FIndex(1,1,2,1,1)))
    ]
end

@testset "fockcouplingnambus" begin
    @test fockcouplingnambus(nothing,('*','*','*','*'),(1,1,2,2))==(0,0,2,1)
    @test fockcouplingnambus(nothing,(0,0,1,1),(1,1,2,2))==(0,0,1,1)
    @test fockcouplingnambus(Val(:Pairing),('*','*'),(2,2))==(1,1)
    @test fockcouplingnambus(Val(:Pairing),(1,2),(2,2))==(1,2)
end

@testset "σ⁰" begin
    @test σ⁰("sp")==FockCoupling{2}(1,spins=(1,1))+FockCoupling{2}(1,spins=(2,2))
    @test σ⁰("ob")==FockCoupling{2}(1,orbitals=(1,1))+FockCoupling{2}(1,orbitals=(2,2))
    @test σ⁰("sl")==FockCoupling{2}(1,atoms=(1,1))+FockCoupling{2}(1,atoms=(2,2))
    @test σ⁰("ph")==FockCoupling{2}(1,nambus=(1,2))+FockCoupling{2}(1,nambus=(2,1))
end

@testset "σˣ" begin
    @test σˣ("sp")==FockCoupling{2}(1,spins=(1,2))+FockCoupling{2}(1,spins=(2,1))
    @test σˣ("ob")==FockCoupling{2}(1,orbitals=(1,2))+FockCoupling{2}(1,orbitals=(2,1))
    @test σˣ("sl")==FockCoupling{2}(1,atoms=(1,2))+FockCoupling{2}(1,atoms=(2,1))
    @test σˣ("ph")==FockCoupling{2}(1,nambus=(1,1))+FockCoupling{2}(1,nambus=(2,2))
end

@testset "σʸ" begin
    @test σʸ("sp")==FockCoupling{2}(1im,spins=(1,2))+FockCoupling{2}(-1im,spins=(2,1))
    @test σʸ("ob")==FockCoupling{2}(1im,orbitals=(1,2))+FockCoupling{2}(-1im,orbitals=(2,1))
    @test σʸ("sl")==FockCoupling{2}(1im,atoms=(1,2))+FockCoupling{2}(-1im,atoms=(2,1))
    @test σʸ("ph")==FockCoupling{2}(1im,nambus=(1,1))+FockCoupling{2}(-1im,nambus=(2,2))
end

@testset "σᶻ" begin
    @test σᶻ("sp")==FockCoupling{2}(-1,spins=(1,1))+FockCoupling{2}(1,spins=(2,2))
    @test σᶻ("ob")==FockCoupling{2}(-1,orbitals=(1,1))+FockCoupling{2}(1,orbitals=(2,2))
    @test σᶻ("sl")==FockCoupling{2}(-1,atoms=(1,1))+FockCoupling{2}(1,atoms=(2,2))
    @test σᶻ("ph")==FockCoupling{2}(-1,nambus=(1,2))+FockCoupling{2}(1,nambus=(2,1))
end

@testset "σ⁺" begin
    @test σ⁺("sp")==Couplings(FockCoupling{2}(1,spins=(2,1)))
    @test σ⁺("ob")==Couplings(FockCoupling{2}(1,orbitals=(2,1)))
    @test σ⁺("sl")==Couplings(FockCoupling{2}(1,atoms=(2,1)))
    @test σ⁺("ph")==Couplings(FockCoupling{2}(1,nambus=(2,2)))
end

@testset "σ⁻" begin
    @test σ⁻("sp")==Couplings(FockCoupling{2}(1,spins=(1,2)))
    @test σ⁻("ob")==Couplings(FockCoupling{2}(1,orbitals=(1,2)))
    @test σ⁻("sl")==Couplings(FockCoupling{2}(1,atoms=(1,2)))
    @test σ⁻("ph")==Couplings(FockCoupling{2}(1,nambus=(1,1)))
end

@testset "fockcoupling" begin
    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1)@(1,1,2,2) with (α<β,σ≠γ)"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1)@(1,1,2,2) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1)@(1,1,2,2)"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1)@(1,1,2,2) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1)"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with (α<β,σ≠γ)"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with α<β,σ≠γ"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with (α<β,*)"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(σ,γ,σ,γ)⊗ph(2,1,2,1) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗ph(2,1,2,1) with α<β"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗ph(2,1,2,1) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(2,1,2,1) with α<β"
    ob,sp=fc.obsubscripts[1].identifier,fc.spsubscripts[1].identifier
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(α,α,β,β)⊗sp(2,1,2,1) with ($ob,$sp) && ($ob,$sp) && ($ob,$sp) && ($ob,$sp)"

    fc=fc"1.0 sl(1,1,1,1)⊗ob(1,1,1,1)⊗ph(2,1,2,1)"
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ob(1,1,1,1)⊗ph(2,1,2,1)"

    fc=fc"1.0 sl(1,1,1,1)⊗ph(2,1,2,1)"
    @test repr(fc)=="1.0 sl(1,1,1,1)⊗ph(2,1,2,1)"

    fc=fc"1.0 @(1,1,1,1)"
    @test repr(fc)=="1.0 @(1,1,1,1)"

    fc=fc"1.0im {2}"
    @test repr(fc)=="1.0im {2}"
end

@testset "fockcouplings" begin
    @test σ⁰"sp"==σ⁰("sp") && σ⁰"sp@(1,2)"==σ⁰("sp",centers=(1,2))
    @test σ⁰"ob"==σ⁰("ob") && σ⁰"ob@(1,2)"==σ⁰("ob",centers=(1,2))
    @test σ⁰"sl"==σ⁰("sl") && σ⁰"sl@(1,2)"==σ⁰("sl",centers=(1,2))
    @test σ⁰"ph"==σ⁰("ph") && σ⁰"ph@(1,2)"==σ⁰("ph",centers=(1,2))

    @test σˣ"sp"==σˣ("sp") && σˣ"sp@(1,2)"==σˣ("sp",centers=(1,2))
    @test σˣ"ob"==σˣ("ob") && σˣ"ob@(1,2)"==σˣ("ob",centers=(1,2))
    @test σˣ"sl"==σˣ("sl") && σˣ"sl@(1,2)"==σˣ("sl",centers=(1,2))
    @test σˣ"ph"==σˣ("ph") && σˣ"ph@(1,2)"==σˣ("ph",centers=(1,2))

    @test σʸ"sp"==σʸ("sp") && σʸ"sp@(1,2)"==σʸ("sp",centers=(1,2))
    @test σʸ"ob"==σʸ("ob") && σʸ"ob@(1,2)"==σʸ("ob",centers=(1,2))
    @test σʸ"sl"==σʸ("sl") && σʸ"sl@(1,2)"==σʸ("sl",centers=(1,2))
    @test σʸ"ph"==σʸ("ph") && σʸ"ph@(1,2)"==σʸ("ph",centers=(1,2))

    @test σᶻ"sp"==σᶻ("sp") && σᶻ"sp@(1,2)"==σᶻ("sp",centers=(1,2))
    @test σᶻ"ob"==σᶻ("ob") && σᶻ"ob@(1,2)"==σᶻ("ob",centers=(1,2))
    @test σᶻ"sl"==σᶻ("sl") && σᶻ"sl@(1,2)"==σᶻ("sl",centers=(1,2))
    @test σᶻ"ph"==σᶻ("ph") && σᶻ"ph@(1,2)"==σᶻ("ph",centers=(1,2))

    @test σ⁺"sp"==σ⁺("sp") && σ⁺"sp@(1,2)"==σ⁺("sp",centers=(1,2))
    @test σ⁺"ob"==σ⁺("ob") && σ⁺"ob@(1,2)"==σ⁺("ob",centers=(1,2))
    @test σ⁺"sl"==σ⁺("sl") && σ⁺"sl@(1,2)"==σ⁺("sl",centers=(1,2))
    @test σ⁺"ph"==σ⁺("ph") && σ⁺"ph@(1,2)"==σ⁺("ph",centers=(1,2))

    @test σ⁻"sp"==σ⁻("sp") && σ⁻"sp@(1,2)"==σ⁻("sp",centers=(1,2))
    @test σ⁻"ob"==σ⁻("ob") && σ⁻"ob@(1,2)"==σ⁻("ob",centers=(1,2))
    @test σ⁻"sl"==σ⁻("sl") && σ⁻"sl@(1,2)"==σ⁻("sl",centers=(1,2))
    @test σ⁻"ph"==σ⁻("ph") && σ⁻"ph@(1,2)"==σ⁻("ph",centers=(1,2))
end

@testset "Onsite" begin
    term=Onsite{'F'}(:mu,1.5)
    @test term|>abbr==:st
    @test otype(term|>typeof,OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Nothing})==FOperator{Float,ID{OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Nothing},2}}
    @test otype(term|>typeof,OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Int})==FOperator{Float,ID{OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Int},2}}

    term=Onsite{'B'}(:mu,1.5,couplings=σˣ("sp")⊗σᶻ("ob"),modulate=true)
    @test term|>abbr==:st
    @test otype(term|>typeof,OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Nothing})==BOperator{Float,ID{OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Nothing},2}}
    @test otype(term|>typeof,OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Int})==BOperator{Float,ID{OID{FIndex{Int},SVector{2,Float},SVector{2,Float},Int},2}}

    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[point.pid])
    table=Table(config,usualfockindextotuple)

    term=Onsite{'F'}(:mu,1.5,couplings=σˣ("sp")⊗σᶻ("ob"),modulate=true)
    operators=Operators(FOperator(+1.5,ID(OID(FIndex('a',1,2,2,2),[0.5,0.5],[0.0,0.0],4),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3))),
                        FOperator(-1.5,ID(OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2)))
    )
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators+operators'
    operators=Operators(FOperator(+1.5,ID(OID(FIndex('a',1,2,2,2),[0.5,0.5],[0.0,0.0],nothing),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],nothing))),
                        FOperator(-1.5,ID(OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],nothing),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],nothing)))
    )
    @test expand(term,point,config,nothing,true)==operators
    @test expand(term,point,config,nothing,false)==operators+operators'

    term=Onsite{'F'}(:mu,1.5,couplings=σᶻ("sp")⊗σᶻ("ob"),modulate=true)
    operators=Operators(FOperator(-0.75,ID(OID(FIndex('a',1,2,1,2),[0.5,0.5],[0.0,0.0],3),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3))),
                        FOperator(-0.75,ID(OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2))),
                        FOperator(+0.75,ID(OID(FIndex('a',1,2,2,2),[0.5,0.5],[0.0,0.0],4),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4))),
                        FOperator(+0.75,ID(OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1)))
    )
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators+operators'
end

@testset "Hopping" begin
    bond=Bond(1,Point(PID('a',1),(0.5,0.5),(0.0,0.0)),Point(PID('b',2),(0.0,0.0),(0.0,0.0)))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[bond.spoint.pid,bond.epoint.pid])
    table=Table(config,usualfockindextotuple)
    term=Hopping{'F'}(:t,1.5,1)
    operators=Operators(FOperator(1.5,ID(OID(FIndex('b',2,2,2,2),[0.0,0.0],[0.0,0.0],8),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4))),
                        FOperator(1.5,ID(OID(FIndex('b',2,2,1,2),[0.0,0.0],[0.0,0.0],7),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3))),
                        FOperator(1.5,ID(OID(FIndex('b',2,1,1,2),[0.0,0.0],[0.0,0.0],5),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1))),
                        FOperator(1.5,ID(OID(FIndex('b',2,1,2,2),[0.0,0.0],[0.0,0.0],6),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2)))
    )
    @test term|>abbr==:hp
    @test expand(term,bond,config,table,true)==operators
    @test expand(term,bond,config,table,false)==operators+operators'
end

@testset "Pairing" begin
    bond=Bond(1,Point(PID('a',1),(0.5,0.5),(0.0,0.0)),Point(PID('b',2),(0.0,0.0),(0.0,0.0)))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=1,nspin=2,nnambu=2),[bond.spoint.pid,bond.epoint.pid])
    table=Table(config,nambufockindextotuple)
    term=Pairing{'F'}(:Δ,1.5,1,couplings=FockCoupling{2}(spins=(2,2)),amplitude=bond->bond|>rcoord|>azimuthd≈45 ? 1 : -1)
    operators=Operators(FOperator(-1.5,ID(OID(FIndex('b',2,1,2,1),[0.0,0.0],[0.0,0.0],6),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2))),
                        FOperator(+1.5,ID(OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2),OID(FIndex('b',2,1,2,1),[0.0,0.0],[0.0,0.0],6)))
    )
    @test term|>abbr==:pr
    @test expand(term,bond,config,table,true)==operators
    @test expand(term,bond,config,table,false)==operators+operators'

    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=1,nspin=2,nnambu=2),[point.pid])
    table=Table(config,nambufockindextotuple)
    term=Pairing{'F'}(:Δ,1.5,0,couplings=FockCoupling{2}(spins=(2,1))-FockCoupling{2}(spins=(1,2)))
    operators=Operators(FOperator(+1.5,ID(OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1))),
                        FOperator(-1.5,ID(OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2)))
    )
    @test term|>abbr==:pr
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators+operators'
end

@testset "Hubbard" begin
    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[point.pid])
    table=Table(config,usualfockindextotuple)
    term=Hubbard{'F'}(:H,2.5)
    operators=Operators(FOperator(1.25,ID(  OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2),
                                            OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1))),
                        FOperator(1.25,ID(  OID(FIndex('a',1,2,2,2),[0.5,0.5],[0.0,0.0],4),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4),
                                            OID(FIndex('a',1,2,1,2),[0.5,0.5],[0.0,0.0],3),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3)))
    )
    @test term|>abbr==:hb
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators*2
end

@testset "InterOrbitalInterSpin" begin
    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[point.pid])
    table=Table(config,usualfockindextotuple)
    term=InterOrbitalInterSpin{'F'}(:H,2.5)
    operators=Operators(FOperator(1.25,ID(  OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2),
                                            OID(FIndex('a',1,2,1,2),[0.5,0.5],[0.0,0.0],3),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3))),
                        FOperator(1.25,ID(  OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1),
                                            OID(FIndex('a',1,2,2,2),[0.5,0.5],[0.0,0.0],4),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4)))
    )
    @test term|>abbr==:nons
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators*2
end

@testset "InterOrbitalIntraSpin" begin
    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[point.pid])
    table=Table(config,usualfockindextotuple)
    term=InterOrbitalIntraSpin{'F'}(:H,2.5)
    operators=Operators(FOperator(1.25,ID(  OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1),
                                            OID(FIndex('a',1,2,1,2),[0.5,0.5],[0.0,0.0],3),OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3))),
                        FOperator(1.25,ID(  OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2),
                                            OID(FIndex('a',1,2,2,2),[0.5,0.5],[0.0,0.0],4),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4)))
    )
    @test term|>abbr==:noes
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators*2
end

@testset "SpinFlip" begin
    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[point.pid])
    table=Table(config,usualfockindextotuple)
    term=SpinFlip{'F'}(:H,2.5)
    operators=Operators(FOperator(2.5,ID(   OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,2,1,2),[0.5,0.5],[0.0,0.0],3),
                                            OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4)))
    )
    @test term|>abbr==:sf
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators+operators'
end

@testset "PairHopping" begin
    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=2,nspin=2,nnambu=2),[point.pid])
    table=Table(config,usualfockindextotuple)
    term=PairHopping{'F'}(:H,2.5)
    operators=Operators(FOperator(2.5,ID(   OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),
                                            OID(FIndex('a',1,2,1,1),[0.5,0.5],[0.0,0.0],3),OID(FIndex('a',1,2,2,1),[0.5,0.5],[0.0,0.0],4)))
    )
    @test term|>abbr==:ph
    @test expand(term,point,config,table,true)==operators
    @test expand(term,point,config,table,false)==operators+operators'
end

@testset "Coulomb" begin
    bond=Bond(1,Point(PID('a',1),(0.5,0.5),(0.0,0.0)),Point(PID('b',2),(0.0,0.0),(0.0,0.0)))
    config=IDFConfig{Fock}(pid->Fock(atom=pid.site%2,norbital=1,nspin=2,nnambu=2),[bond.spoint.pid,bond.epoint.pid])
    table=Table(config,usualfockindextotuple)

    term=Coulomb{'F'}(:V,2.5,1,couplings=σᶻ("sp")*σᶻ("sp"))
    operators=Operators(FOperator(-1.25,ID( OID(FIndex('b',2,1,1,2),[0.0,0.0],[0.0,0.0],3),OID(FIndex('b',2,1,1,1),[0.0,0.0],[0.0,0.0],3),
                                            OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2))),
                        FOperator(+1.25,ID( OID(FIndex('b',2,1,1,2),[0.0,0.0],[0.0,0.0],3),OID(FIndex('b',2,1,1,1),[0.0,0.0],[0.0,0.0],3),
                                            OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1))),
                        FOperator(-1.25,ID( OID(FIndex('b',2,1,2,2),[0.0,0.0],[0.0,0.0],4),OID(FIndex('b',2,1,2,1),[0.0,0.0],[0.0,0.0],4),
                                            OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1))),
                        FOperator(+1.25,ID( OID(FIndex('b',2,1,2,2),[0.0,0.0],[0.0,0.0],4),OID(FIndex('b',2,1,2,1),[0.0,0.0],[0.0,0.0],4),
                                            OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2)))
    )
    @test term|>abbr==:cl
    @test expand(term,bond,config,table,true)==operators
    @test expand(term,bond,config,table,false)==operators+operators'

    term=Coulomb{'F'}(:V,2.5,1,couplings=σˣ("sp")*σᶻ("sp"))
    operators=Operators(FOperator(-2.5,ID(  OID(FIndex('b',2,1,2,2),[0.0,0.0],[0.0,0.0],4),OID(FIndex('b',2,1,1,1),[0.0,0.0],[0.0,0.0],3),
                                            OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1))),
                        FOperator(+2.5,ID(  OID(FIndex('b',2,1,1,2),[0.0,0.0],[0.0,0.0],3),OID(FIndex('b',2,1,2,1),[0.0,0.0],[0.0,0.0],4),
                                            OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2))),
                        FOperator(+2.5,ID(  OID(FIndex('b',2,1,2,2),[0.0,0.0],[0.0,0.0],4),OID(FIndex('b',2,1,1,1),[0.0,0.0],[0.0,0.0],3),
                                            OID(FIndex('a',1,1,2,2),[0.5,0.5],[0.0,0.0],2),OID(FIndex('a',1,1,2,1),[0.5,0.5],[0.0,0.0],2))),
                        FOperator(-2.5,ID(  OID(FIndex('b',2,1,1,2),[0.0,0.0],[0.0,0.0],3),OID(FIndex('b',2,1,2,1),[0.0,0.0],[0.0,0.0],4),
                                            OID(FIndex('a',1,1,1,2),[0.5,0.5],[0.0,0.0],1),OID(FIndex('a',1,1,1,1),[0.5,0.5],[0.0,0.0],1)))
    )
    @test term|>abbr==:cl
    @test expand(term,bond,config,table,true)==operators
    @test expand(term,bond,config,table,false)==operators+operators'
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -0.9705716506753843 <= q <= 1.6199082189427347)

                     add_NL_constraint(m, :(gelu(0.27812717705086953 + -0.4967460082986306*gelu(-0.6740257057257448 + -0.0007917455797805673*gelu(-0.08584163188401472 + 0.6683752749724832*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*gelu(-0.7769463282105753 + 0.25325339074048747*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + -0.883304587160544*gelu(0.3709718635886192 + -0.47514710764751555*gelu(-0.08584163188401472 + 0.6683752749724832*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*gelu(-0.7769463282105753 + 0.25325339074048747*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) + gelu(0.48450040123340576 + -0.32897580116872405*gelu(-0.6740257057257448 + -0.0007917455797805673*gelu(-0.08584163188401472 + 0.6683752749724832*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*gelu(-0.7769463282105753 + 0.25325339074048747*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + 0.43090905819531766*gelu(0.3709718635886192 + -0.47514710764751555*gelu(-0.08584163188401472 + 0.6683752749724832*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*gelu(-0.7769463282105753 + 0.25325339074048747*gelu(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*gelu(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
const Re       = 6370e3;                      # Earth radius, f64
const c        = 3e8;                         # speedo lite, f64
eta(L)         = Omegace0*L*Re/c;

B0 = 188;                     # nT, measured from THEMIS A at L of 5.2 (conjunction w/ ELFIN)
L = 5.2;
Ew = [21.358639,36.114532,25.813985,48.348956,32.235131,17.927961,23.999402]; # mV/m, measured from THEMIS E bursts

Omegace0 = (1e-9*1.6e-19*B0)/(9.11e-31);  # electron gyrofreq @ the equator
Omegape = 8; # from THEMIS A density measurements
omegam = 0.22; # from THEMIS A at L of 5.2

N = (Omegape/omegam)*((1/omegam)-1)^(-1/2)

Bw = (Ew .* N ./ 300 ) # nT

wave_amp_modifier = ( Bw ./ B0 )*eta(L)
avg_wave_amp_mod = sum(wave_amp_modifier)/(2*length(wave_amp_modifier))

using Plots, LaTeXStrings

u(lambda) = .5*(tanh(deg2rad(lambda)/deg2rad(1))+1)
u20(lambda) = 0.929745945424938^-1*tanh((deg2rad(lambda)/(deg2rad(2)))) * (exp(-(deg2rad(lambda)/(deg2rad(20)))^2));
u30(lambda) = 0.9597545683166291^-1*tanh((deg2rad(lambda)/(deg2rad(2)))) * (exp(-(deg2rad(lambda)/(deg2rad(30)))^2));
u40(lambda) = 0.9733653674717542^-1*tanh((deg2rad(lambda)/(deg2rad(2)))) * (exp(-(deg2rad(lambda)/(deg2rad(40)))^2));

plot(u, 0,90,label=L"0.5\tanh(\lambda)/1 + 1")
plot!(u20, 0,90,label=L"\delta\lambda_2 = 20\degree")
plot!(u30, 0,90,label=L"\delta\lambda_2 = 30\degree")
plot!(u40, 0,90,label=L"\delta\lambda_2 = 40\degree")
plot!(xlabel="Latitude (degrees)", ylabel="B_w modifier", title="Latitude distribution of chorus waves")
plot!(yscale=:log10, ylim = (0.1,1.))

maximum(u20.(0:.1:90)) # 0.929745945424938^-1*
maximum(u30.(0:.1:90)) # 0.9597545683166291^-1*
maximum(u40.(0:.1:90)) # 0.9733653674717542^-1*

max_stars_repo_name: ethantsai/nlwhistlers

max_issues_repo_name: ethantsai/nlwhistlers

max_forks_repo_name: ethantsai/nlwhistlers


text:
include("hlt.jl")
include("networking.jl")

myID, gameMap = getInit()
sendInit("RandomJuliaBot")

while true
  moves = Vector{Move}()
  gameMap = getFrame()
  for y in 0:gameMap.height-1
    for x in 0:gameMap.width-1
      if getSite(gameMap, Location(x, y)).owner == myID
        push!(moves, Move(Location(x, y), rand(0:4)))
      end
    end
  end
  sendFrame(moves)
end

max_stars_repo_name: yangle/HaliteIO

max_issues_repo_name: yangle/HaliteIO

max_forks_repo_name: yangle/HaliteIO


text: # code for implementing the metropolis algorithm (right now only for gauge fields)


text:
"""
# Stipple.Layout

Utilities for rendering the general layout of a Stipple app, such as of a data dashboard web page or supporting themes.
"""
module Layout

using Genie, Stipple

export layout
export page, row, cell

export theme
const THEMES = Function[]

"""
    `function layout(output::Union{String,Vector}; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                      head_content::String = "", channel::String = Genie.config.webchannels_default_route) :: String`

Utility for creating a basic web page structure, including doctype as well as <HTML>, <HEAD>, <TITLE>, <META viewport>,
  and <BODY> tags, together with the output content.

If `partial` is `true`, the page layout HTML elements are not returned.

### Examples

```julia
julia> layout([
        span("Hello", @text(:greeting))
        ])
"<!DOCTYPE html>\n<html><head><title></title><meta name=\"viewport\" content=\"width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui\" /></head><body class style><link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><span v-text='greeting'>Hello</span><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/underscore-min.js\"></script><script src=\"/js/vue.js\"></script><script src=\"/js/quasar.umd.min.js\"></script>\n<script src=\"/js/apexcharts.min.js\"></script><script src=\"/js/vue-apexcharts.min.js\"></script><script src=\"/js/stipplecore.js\" defer></script><script src=\"/js/vue_filters.js\" defer></script></body></html>"
```

```julia
julia> layout([
        span("Hello", @text(:greeting))
        ], partial = true)
"<link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><span v-text='greeting'>Hello</span><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/underscore-min.js\"></script><script src=\"/js/vue.js\"></script><script src=\"/js/quasar.umd.min.js\"></script>\n<script src=\"/js/apexcharts.min.js\"></script><script src=\"/js/vue-apexcharts.min.js\"></script><script src=\"/js/stipplecore.js\" defer></script><script src=\"/js/vue_filters.js\" defer></script>"
```
"""
function layout(output::Union{S,Vector}; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                head_content::String = "", channel::String = Genie.config.webchannels_default_route,
                core_theme::Bool = true)::ParsedHTMLString where {S<:AbstractString}

  isa(output, Vector) && (output = join(output, '\n'))

  content = string(
    output,
    theme(; core_theme = core_theme),
    Stipple.deps(channel; core_theme = core_theme)
  )

  partial && return content

  Genie.Renderer.Html.doc(
    Genie.Renderer.Html.html([
      Genie.Renderer.Html.head([
        Genie.Renderer.Html.title(title)
        Genie.Renderer.Html.meta(name="viewport", content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no")
        head_content
      ])
      Genie.Renderer.Html.body(content, class=class, style=style)
    ])
  )
end

"""
    `function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                    channel::String = Genie.config.webchannels_default_route , head_content::String = "", kwargs...)`

Generates the HTML code corresponding to an SPA (a single page application), defining the root element of the Vue app.

### Example

```julia
julia> page(:elemid, [
        span("Hello", @text(:greeting))
        ])
"<!DOCTYPE html>\n<html><head><title></title><meta name=\"viewport\" content=\"width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui\" /></head><body class style><link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><div id=elemid><span v-text='greeting'>Hello</span></div><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/underscore-min.js\"></script><script src=\"/js/vue.js\"></script><script src=\"/js/quasar.umd.min.js\"></script>\n<script src=\"/js/apexcharts.min.js\"></script><script src=\"/js/vue-apexcharts.min.js\"></script><script src=\"/js/stipplecore.js\" defer></script><script src=\"/js/vue_filters.js\" defer></script></body></html>"
```
"""
function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "container", style::String = "",
              channel::String = Genie.config.webchannels_default_route, head_content::String = "",
              prepend::Union{S,Vector} = "", append::Union{T,Vector} = [],
              core_theme::Bool = true, kwargs...)::ParsedHTMLString where {S<:AbstractString,T<:AbstractString}

  layout(
    [
      join(prepend)
      Genie.Renderer.Html.div(id = elemid, args...; class = class, kwargs...)
      join(append)
    ],
    partial = partial, title = title, style = style, head_content = head_content, channel = channel,
    core_theme = core_theme)
end
function page(model::T, args...; kwargs...)::ParsedHTMLString where {T<:Stipple.ReactiveModel}
  page(vm(model), args...; channel = getchannel(model), kwargs...)
end

"""
    `function row(args...; kwargs...)`

Creates a `div` HTML element with a CSS class named `row`. This works with Stipple's Twitter Bootstrap to create the
responsive CSS grid of the web page. The `row` function creates rows which should include `cell`s.

### Example

```julia
julia> row(span("Hello"))
"<div class=\"row\"><span>Hello</span></div>"
```
"""
function row(args...; kwargs...)
  kwargs = NamedTuple(Dict{Symbol,Any}(kwargs...), :class, "row")

  Genie.Renderer.Html.div(args...; kwargs...)
end

"""
    `function cell(args...; size::Int=0, kwargs...)`

Creates a `div` HTML element with CSS classes named `col col-12` and `col-sm-$size`.
This works with Stipple's Twitter Bootstrap to create the responsive CSS grid of the web page. The `cell`s should be
included within `row`s.

### Example

```julia
julia> row(cell(size=2, span("Hello")))
"<div class=\"row\"><div class=\"col col-12 col-sm-2\"><span>Hello</span></div></div>"
```
"""
function cell(args...; size::Union{Int,AbstractString} = 0, kwargs...)
  isa(size, AbstractString) && (size = parse(Int, size))

  kwargs = NamedTuple(Dict{Symbol,Any}(kwargs...), :class, "col col-12 col-sm$(size > 0 ? "-$size" : "")")

  Genie.Renderer.Html.div(args...; kwargs...)
end

"""
    `function theme() :: String

Provides theming support for Stipple apps and pages. It includes Stipple's default CSS files and additional elements,
  in the form of HTML tags, can be injected by pushing to the `Stipple.Layout.THEMES` collection.

### Example

```julia
julia> theme()
"<link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.css\" rel=\"stylesheet\" />"

julia> StippleUI.theme()
"<link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" />"

julia> push!(Stipple.Layout.THEMES, StippleUI.theme)
```
"""
function theme(; core_theme::Bool = true) :: String
  if core_theme
    if ! Genie.Assets.external_assets(Stipple.assets_config)
      Genie.Router.route(Genie.Assets.asset_path(Stipple.assets_config, :css, file="stipplecore")) do
        Genie.Renderer.WebRenderable(
          Genie.Assets.embedded(Genie.Assets.asset_file(cwd=dirname(@__DIR__), type="css", file="stipplecore")),
          :css) |> Genie.Renderer.respond
      end
    end

    output = string(
      stylesheet("https://fonts.googleapis.com/css?family=Material+Icons"),
      stylesheet(Genie.Assets.asset_path(Stipple.assets_config, :css, file="stipplecore"))
    )
  end

  string(output, join([f() for f in THEMES], "\n")) |> ParsedHTMLString
end

end