text:
"""
    TensorProduct(kernels...)

Create a tensor product kernel from kernels ``k_1, \\ldots, k_n``, i.e.,
a kernel ``k`` that is given by
```math
k(x, y) = \\prod_{i=1}^n k_i(x_i, y_i).
```

The `kernels` can be specified as individual arguments, a tuple, or an iterable data
structure such as an array. Using a tuple or individual arguments guarantees that
`TensorProduct` is concretely typed but might lead to large compilation times if the
number of kernels is large.
"""
struct TensorProduct{K} <: Kernel
    kernels::K
end

function TensorProduct(kernel::Kernel, kernels::Kernel...)
    return TensorProduct((kernel, kernels...))
end

@functor TensorProduct

Base.length(kernel::TensorProduct) = length(kernel.kernels)

function (kernel::TensorProduct)(x, y)
    if !(length(x) == length(y) == length(kernel))
        throw(DimensionMismatch("number of kernels and number of features
are not consistent"))
    end
    return prod(k(xi, yi) for (k, xi, yi) in zip(kernel.kernels, x, y))
end

function validate_domain(k::TensorProduct, x::AbstractVector)
    dim(x) == length(k) ||
        error("number of kernels and groups of features are not consistent")
end

# Utility for slicing up inputs.
slices(x::AbstractVector{<:Real}) = (x,)
slices(x::ColVecs) = eachrow(x.X)
slices(x::RowVecs) = eachcol(x.X)

function kernelmatrix!(K::AbstractMatrix, k::TensorProduct, x::AbstractVector)
    validate_inplace_dims(K, x)
    validate_domain(k, x)

    kernels_and_inputs = zip(k.kernels, slices(x))
    kernelmatrix!(K, first(kernels_and_inputs)...)
    for (k, xi) in Iterators.drop(kernels_and_inputs, 1)
        K .*= kernelmatrix(k, xi)
    end

    return K
end

function kernelmatrix!(
    K::AbstractMatrix, k::TensorProduct, x::AbstractVector, y::AbstractVector,
)
    validate_inplace_dims(K, x, y)
    validate_domain(k, x)

    kernels_and_inputs = zip(k.kernels, slices(x), slices(y))
    kernelmatrix!(K, first(kernels_and_inputs)...)
    for (k, xi, yi) in Iterators.drop(kernels_and_inputs, 1)
        K .*= kernelmatrix(k, xi, yi)
    end

    return K
end

function kerneldiagmatrix!(K::AbstractVector, k::TensorProduct, x::AbstractVector)
    validate_inplace_dims(K, x)
    validate_domain(k, x)

    kernels_and_inputs = zip(k.kernels, slices(x))
    kerneldiagmatrix!(K, first(kernels_and_inputs)...)
    for (k, xi) in Iterators.drop(kernels_and_inputs, 1)
        K .*= kerneldiagmatrix(k, xi)
    end

    return K
end

function kernelmatrix(k::TensorProduct, x::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kernelmatrix, hadamard, k.kernels, slices(x))
end

function kernelmatrix(k::TensorProduct, x::AbstractVector, y::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kernelmatrix, hadamard, k.kernels, slices(x), slices(y))
end

function kerneldiagmatrix(k::TensorProduct, x::AbstractVector)
    validate_domain(k, x)
    return mapreduce(kerneldiagmatrix, hadamard, k.kernels, slices(x))
end

Base.show(io::IO, kernel::TensorProduct) = printshifted(io, kernel, 0)

function Base.:(==)(x::TensorProduct, y::TensorProduct)
    return (
        length(x.kernels) == length(y.kernels) &&
        all(kx == ky for (kx, ky) in zip(x.kernels, y.kernels))
    )
end

function printshifted(io::IO, kernel::TensorProduct, shift::Int)
    print(io, "Tensor product of ", length(kernel), " kernels:")
    for k in kernel.kernels
        print(io, "\n")
        for _ in 1:(shift + 1)
            print(io, "\t")
        end
        print(io, "- ")
        printshifted(io, k, shift + 2)
    end
end


text:
#=##############################################################################
# DESCRIPTION
    Non-lifting paneled body types definition.
# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jun 2018
  * License   : MIT License
=###############################################################################

################################################################################
# NON-LIFTING BODY TYPE
################################################################################

"""
  `NonLiftingBody(grid::gt.GridTriangleSurface)`

Non-lifting paneled body that is solved using a constant source distribution.
`grid` is the grid surface (paneled geometry).

  **Properties**
  * `nnodes::Int64`                     : Number of nodes
  * `ncells::Int64`                     : Number of cells
  * `fields::Array{String, 1}`          : Available fields (solutions)
  * `Oaxis::Array{T,2} where {T<:Real}` : Coordinate system of original grid
  * `O::Array{T,1} where {T<:Real}`     : Position of CS of original grid

"""
struct NonLiftingBody <: AbstractBody

  # User inputs
  grid::gt.GridTriangleSurface              # Paneled geometry

  # Properties
  nnodes::Int64                             # Number of nodes
  ncells::Int64                             # Number of cells
  fields::Array{String, 1}                  # Available fields (solutions)
  Oaxis::Array{T1,2} where {T1<:RType}      # Coordinate system of original grid
  O::Array{T2,1} where {T2<:RType}          # Position of CS of original grid

  # Internal variables
  _G::Array{T3,2} where {T3<:RType}         # Geometric solution matrix

  NonLiftingBody( grid,
                  nnodes=grid.nnodes, ncells=grid.ncells,
                    fields=Array{String,1}(),
                    Oaxis=Array(1.0I, 3, 3), O=zeros(3),
                  _G=_calc_Gsource(grid)
         ) = new( grid,
                  nnodes, ncells,
                    fields,
                    Oaxis, O,
                  _G
         )
end

function solve(self::NonLiftingBody, Vinfs::Array{Array{T,1},1}) where {T<:RType}
  if size(Vinfs,1) != self.ncells
    error("Invalid Vinfs; expected size $(self.ncells), got $(size(Vinfs,1))")
  end

  lambda = [-dot(Vinfs[i], get_normal(self, i)) for i in 1:self.ncells]
  sigma = self._G\lambda

  add_field(self, "Vinf", Vinfs)
  add_field(self, "sigma", sigma)
  _solvedflag(self, true)
end

##### INTERNAL FUNCTIONS  ######################################################
function _calc_Gsource(grid::gt.GridTriangleSurface)
  return PanelSolver.G_constant_source(
                  grid.orggrid.nodes,                                  # Nodes
                  [gt.get_cell(grid, i) for i in 1:grid.ncells],       # Panels
                  [_get_controlpoint(grid, i) for i in 1:grid.ncells], # CPs
                  [gt.get_normal(grid, i) for i in 1:grid.ncells],     # Normals
                )
end

"""
Returns the velocity induced by the body on the targets `targets`. It adds the
velocity at the i-th target to out[i].
"""
function _Vind(self::NonLiftingBody, targets::Array{Array{T1,1},1},
                          out::Array{Array{T2,1},1}) where{T1<:RType, T2<:RType}
  # Iterates over panels
  for i in 1:self.ncells
    # Velocity of i-th  panel on every target
    PanelSolver.Vconstant_source(
                    gt.get_cellnodes(self.grid, i),    # Nodes in i-th panel
                    get_fieldval(self, "sigma", i; _check=false),  # Strength
                    targets,                           # Targets
                    out;                               # Outputs
                  )
  end
end
##### END OF NON-LIFTING BODY ##################################################

################################################################################
# NON-LIFTING DOUBLET BODY TYPE
################################################################################
"""
  `NonLiftingBodyDoublet(grid::gt.GridTriangleSurface)`

Non-lifting paneled body that is solved using a constant doublet distribution.
`grid` is the grid surface (paneled geometry).

  **Properties**
  * `nnodes::Int64`                     : Number of nodes
  * `ncells::Int64`                     : Number of cells
  * `fields::Array{String, 1}`          : Available fields (solutions)
  * `Oaxis::Array{T,2} where {T<:Real}` : Coordinate system of original grid
  * `O::Array{T,1} where {T<:Real}`     : Position of CS of original grid

"""
struct NonLiftingBodyDoublet <: AbstractBody

  # User inputs
  grid::gt.GridTriangleSurface              # Paneled geometry

  # Properties
  nnodes::Int64                             # Number of nodes
  ncells::Int64                             # Number of cells
  fields::Array{String, 1}                  # Available fields (solutions)
  Oaxis::Array{T1,2} where {T1<:RType}      # Coordinate system of original grid
  O::Array{T2,1} where {T2<:RType}          # Position of CS of original grid

  # Internal variables
  _G::Array{T3,2} where {T3<:RType}         # Geometric solution matrix

  NonLiftingBodyDoublet( grid,
                  nnodes=grid.nnodes, ncells=grid.ncells,
                    fields=Array{String,1}(),
                    Oaxis=Array(1.0I, 3, 3), O=zeros(3),
                  _G=_calc_Gdoublet(grid)
         ) = new( grid,
                  nnodes, ncells,
                    fields,
                    Oaxis, O,
                  _G
         )
end

function solve(self::NonLiftingBodyDoublet, Vinfs::Array{Array{T,1},1}
                                                              ) where {T<:RType}
  if size(Vinfs,1) != self.ncells
    error("Invalid Vinfs; expected size $(self.ncells), got $(size(Vinfs,1))")
  end

  lambda = [-dot(Vinfs[i], get_normal(self, i)) for i in 1:self.ncells]
  mu = self._G\lambda

  add_field(self, "Vinf", Vinfs)
  add_field(self, "mu", mu)
  _solvedflag(self, true)
end

##### INTERNAL FUNCTIONS  ######################################################
function _calc_Gdoublet(grid::gt.GridTriangleSurface)
  return PanelSolver.G_constant_doublet(
                  grid.orggrid.nodes,                                  # Nodes
                  [gt.get_cell(grid, i) for i in 1:grid.ncells],       # Panels
                  [_get_controlpoint(grid, i) for i in 1:grid.ncells], # CPs
                  [gt.get_normal(grid, i) for i in 1:grid.ncells],     # Normals
                )
end

"""
Returns the velocity induced by the body on the targets `targets`. It adds the
velocity at the i-th target to out[i].
"""
function _Vind(self::NonLiftingBodyDoublet, targets::Array{Array{T1,1},1},
                          out::Array{Array{T2,1},1}) where{T1<:RType, T2<:RType}
  # Iterates over panels
  for i in 1:self.ncells
    # Velocity of i-th  panel on every target
    PanelSolver.Vconstant_doublet(
                    gt.get_cellnodes(self.grid, i),    # Nodes in i-th panel
                    get_fieldval(self, "mu", i; _check=false),  # Strength
                    targets,                           # Targets
                    out;                               # Outputs
                  )
  end
end
##### END OF NON-LIFTING DOUBLET BODY ##########################################

################################################################################
# NON-LIFTING VORTEX-RING BODY TYPE
################################################################################
"""
  `NonLiftingBodyVRing(grid::gt.GridTriangleSurface)`

Non-lifting paneled body that is solved using vortex ring panels.
`grid` is the grid surface (paneled geometry).

  **Properties**
  * `nnodes::Int64`                     : Number of nodes
  * `ncells::Int64`                     : Number of cells
  * `fields::Array{String, 1}`          : Available fields (solutions)
  * `Oaxis::Array{T,2} where {T<:Real}` : Coordinate system of original grid
  * `O::Array{T,1} where {T<:Real}`     : Position of CS of original grid

"""
struct NonLiftingBodyVRing <: AbstractBody

  # User inputs
  grid::gt.GridTriangleSurface              # Paneled geometry

  # Properties
  nnodes::Int64                             # Number of nodes
  ncells::Int64                             # Number of cells
  fields::Array{String, 1}                  # Available fields (solutions)
  Oaxis::Array{T1,2} where {T1<:RType}      # Coordinate system of original grid
  O::Array{T2,1} where {T2<:RType}          # Position of CS of original grid

  # Internal variables
  _G::Array{T3,2} where {T3<:RType}         # Geometric solution matrix

  NonLiftingBodyVRing( grid,
                  nnodes=grid.nnodes, ncells=grid.ncells,
                    fields=Array{String,1}(),
                    Oaxis=Array(1.0I, 3, 3), O=zeros(3),
                  _G=_calc_Gvring(grid)
         ) = new( grid,
                  nnodes, ncells,
                    fields,
                    Oaxis, O,
                  _G
         )
end

function solve(self::NonLiftingBodyVRing, Vinfs::Array{Array{T,1},1}
                                                              ) where {T<:RType}
  if size(Vinfs,1) != self.ncells
    error("Invalid Vinfs; expected size $(self.ncells), got $(size(Vinfs,1))")
  end

  lambda = [-dot(Vinfs[i], get_normal(self, i)) for i in 1:self.ncells]
  Gamma = self._G\lambda

  add_field(self, "Vinf", Vinfs)
  add_field(self, "Gamma", Gamma)
  _solvedflag(self, true)
end

##### INTERNAL FUNCTIONS  ######################################################
function _calc_Gvring(grid::gt.GridTriangleSurface)
  return PanelSolver.G_vortexring(
                  grid.orggrid.nodes,                                  # Nodes
                  [gt.get_cell(grid, i) for i in 1:grid.ncells],       # Panels
                  [_get_controlpoint(grid, i) for i in 1:grid.ncells], # CPs
                  [gt.get_normal(grid, i) for i in 1:grid.ncells],     # Normals
                )
end

"""
Returns the velocity induced by the body on the targets `targets`. It adds the
velocity at the i-th target to out[i].
"""
function _Vind(self::NonLiftingBodyVRing, targets::Array{Array{T1,1},1},
                          out::Array{Array{T2,1},1}) where{T1<:RType, T2<:RType}
  # Iterates over panels
  for i in 1:self.ncells
    # Velocity of i-th  panel on every target
    PanelSolver.Vvortexring(
                    gt.get_cellnodes(self.grid, i),    # Nodes in i-th panel
                    get_fieldval(self, "Gamma", i; _check=false),  # Strength
                    targets,                           # Targets
                    out;                               # Outputs
                  )
  end
end
##### END OF NON-LIFTING VORTEX-RING BODY ######################################

################################################################################
# COMMON FUNCTIONS
################################################################################
"""
  `generate_loft_nonliftbody(args...; optargs...)`
Generates a lofted non-lifting body. See documentation of
`GeometricTools.generate_loft` for a description of the arguments of this
function.
"""
function generate_loft_nonliftbody(args...; vortexring=false, dimsplit::Int64=2,
                                                                    optargs...)
  # Lofts the surface geometry
  grid = gt.generate_loft(args...; optargs...)

  # Splits the quadrialateral panels into triangles
  # dimsplit = 2              # Dimension along which to split
  triang_grid = gt.GridTriangleSurface(grid, dimsplit)

  if vortexring
    return NonLiftingBodyVRing(triang_grid)
  else
    return NonLiftingBody(triang_grid)
  end
end

"""
  `generate_revolution_nonliftbody(args...; optargs...)`
Generates a non-lifting body of a body of revolution. See documentation of
`GeometricTools.surface_revolution` for a description of the arguments of this
function.
"""
function generate_revolution_nonliftbody(args...; vortexring=false,
                                                  dimsplit::Int64=2,
                                                  loop_dim::Int64=2, optargs...)
  # Revolves the geometry
  grid = gt.surface_revolution(args...; loop_dim=loop_dim, optargs...)

  # Splits the quadrialateral panels into triangles
  # dimsplit = 2              # Dimension along which to split
  triang_grid = gt.GridTriangleSurface(grid, dimsplit)

  if vortexring
    return NonLiftingBodyVRing(triang_grid)
  else
    return NonLiftingBody(triang_grid)
  end
end
##### END OF COMMON FUNTIONS ###################################################

ext: jl


text:
function lengthstring(folders)
    n = 0;
    for i = 1:length(folders)
        n = max(length(folders[i]), n)
    end
    return n
end

function replacestring(str1, str2)
     str3 = collect(str1)
     str4 = collect(str2)
     for i = 1:length(str4)
        str3[i] = str4[i]
     end
     return join(str3)
end

function num2string(a, n)
    s = collect(string(a))
    if length(s) < n
        s2 = repeat('0',(n-length(s)))
        s = join(s) * s2
    else
        for i = length(s):-1:(n+1)
            deleteat!(s, i)
        end
        s = join(s)
    end
    return s
end

function printerrors(folders, errors, errstr)

    n = max(lengthstring(folders), length(errstr))
    str1 = repeat(' ', n) * " |   ";

    num = 16;
    print("\n" * repeat('-', 88) * "\n")
    print(replacestring(str1 * "       MAE          |          RMSE          |          RSQ          |\n", errstr))
    print(repeat('-', 88) * "\n")
    e1 = errors[1,1]; e2 = errors[1,2]; e3 = errors[1,3];
    print(replacestring(str1 * num2string(e1, num) * "    |    " * num2string(e2, num) * "    |    " * num2string(e3, num) * "  " * " |\n", "ALL"))
    for i = 1:length(folders)
        a1 = errors[i+1,1]; a2 = errors[i+1,2]; a3 = errors[i+1,3];
        print(replacestring(str1 * num2string(a1, num) * "    |    " * num2string(a2, num) * "    |    " * num2string(a3, num) * "  " * " |\n", folders[i]))
    end
    print(repeat('-', 88) * "\n")

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
m = page_model()
load_parameters(m)

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
using ValidatedNumerics
using IntervalConstraintProgramming

S = @separator 1 <= x^2 + y^2 <= 3
X = IntervalBox(-10..10, -10..10)

@time inner, boundary = pave(S, X, 0.125)

@show length(inner), length(boundary)

include("draw_boxes.jl")

draw_boxes(inner, "green", 0.5, 1)
draw_boxes(boundary, "grey", 0.2)
axis("image")  # set aspect ratio


text:
using Manopt, Documenter, Literate

# generate examples using Literate
tutorialsInputPath = joinpath(@__DIR__, "..", "src/tutorials")
tutorialsRelativePath = "tutorials/"
tutorialsOutputPath = joinpath(@__DIR__,"src/"*tutorialsRelativePath)
tutorials  = [
    "MeanAndMedian",
    "GradientOfSecondOrderDifference",
    "JacobiFields",
    ]
menuEntries = [
    "Getting Started: Optimize!",
    "Gradient of \$d_2\$",
    "Jacobi Fields"]
TutorialMenu = Array{Pair{String,String},1}()
for (i,tutorial) in enumerate(tutorials)
    global TutorialMenu
    sourceFile = joinpath(tutorialsInputPath,tutorial*".jl")
    targetFile = joinpath(tutorialsOutputPath,tutorial*"md")
    Literate.markdown(sourceFile,tutorialsOutputPath; name=tutorial,
    # codefence = "```julia" => "```",
    credit=false)
    push!(TutorialMenu, menuEntries[i] => joinpath(tutorialsRelativePath,tutorial*".md") )
end
makedocs(
    # for development, we disable prettyurls
    # format = Documenter.HTML(prettyurls = false),
    modules = [Manopt],
    sitename = "Manopt.jl",
    pages = [
        "Home" => "index.md",
        "About" => "about.md",
        "Manifolds" => [
            "Introduction" => "manifolds/index.md",
            "Combinations of Manifolds" => "manifolds/combined.md",
            "The Circle \$\\mathbb S^1\$" => "manifolds/circle.md",
            "The Euclidean Space \$\\mathbb R^n\$" => "manifolds/euclidean.md",
            "The Grassmannian Manifold \$\\mathrm{Gr}(k,n)\$" => "manifolds/grassmannian.md",
            "The Hyperbolic Space \$\\mathbb H^n\$" => "manifolds/hyperbolic.md",
            "The Special Orthogonal Group \$\\mathrm{SO}(n)\$" => "manifolds/rotations.md",
            "The Sphere \$\\mathbb S^n\$" => "manifolds/sphere.md",
            "The Stiefel Manifold \$\\mathrm{St}(k,n)\$" => "manifolds/stiefel.md",
            "The Symmetric Matrices \$\\mathrm{Sym}(n)\$" => "manifolds/symmetric.md",
            "The Symmetric Positive Definite Matrices \$\\mathcal P(n)\$" => "manifolds/symmetricpositivedefinite.md",
        ],
        "Plans" => "plans/index.md",
        "Solvers" => [
            "Introduction" => "solvers/index.md",
            "Cyclic Proximal Point" => "solvers/cyclicProximalPoint.md",
            "Douglas–Rachford" => "solvers/DouglasRachford.md",
            "Gradient Descent" => "solvers/gradientDescent.md",
            "Subgradient Method" => "solvers/subGradientMethod.md",
         ],
        "Functions" => [
            "Introduction" => "functions/index.md",
            "cost functions" => "functions/costFunctions.md",
            "Differentials" => "functions/differentials.md",
            "Adjoint Differentials" => "functions/adjointDifferentials.md",
            "Gradients" => "functions/gradients.md",
            "JacobiFields" => "functions/jacobiFields.md",
            "Proximal Maps" => "functions/proximalMaps.md"
        ],
        "Helpers" => [
            "Data" => "helpers/data.md",
            "Error Measures" => "helpers/errorMeasures.md",
            "Exports" => "helpers/exports.md",
        ],
        "Tutorials" => TutorialMenu,
        "Function Index" => "list.md",
    ]
)
deploydocs(
    repo   = "github.com/kellertuer/Manopt.jl",
   # devbranch = "development"
)


text:
module Layer3Tests

using Test
using LinearAlgebra

using PicoQuant
include("test_utils.jl")

@testset "Layer 3 tests" begin
    @testset "Test loading and conversion of qasm" begin

        qasm_str = """OPENQASM 2.0;
                    include "qelib1.inc";
                    qreg q[3];
                    h q[0];
                    cx q[0],q[1];
                    cx q[1],q[2];"""

        circ = load_qasm_as_circuit(qasm_str)
        tng = convert_qiskit_circ_to_network(circ, InteractiveBackend())
        tng_json = to_json(tng)

        @test begin
            circ.n_qubits == 3
        end

        @test begin
            length(circ.data) == 3
        end

        @test begin
            # test conversion to json and back maintains graph
            tng_to_json_from_json = to_json(network_from_json(tng_json))
            tng_to_json_from_json == tng_json
        end
    end

    @testset "Test tensor network data structure" begin
        tn = TensorNetworkCircuit(3, InteractiveBackend())

        @test begin
            # Check if tn has the correct number of nodes
            length(tn.nodes) == 0
        end

        hadamard = [1 1; 1 -1]./sqrt(2)
        add_gate!(tn, hadamard, [1])
        @test begin
            length(tn.nodes) == 1
        end

        @test begin
            length(edges(tn)) == 4
        end
    end

    @testset "Test tensor decomposition" begin
        tn = TensorNetworkCircuit(2, InteractiveBackend())
        gate_data = rand(ComplexF64, 2, 2, 2, 2)

        @test begin
            gate_labels = add_gate!(tn, gate_data, [1, 2], decompose=true)
            result_label = contract_pair!(tn, gate_labels...)
            # should now only be a single node left
            # index order will have changed so permute back before comparing
            # data = reshape(collect(values(tn.backend.tensors))[1], 2, 2, 2, 2)
            data = tn.backend.tensors[result_label]
            data = permutedims(data, (1, 3, 2, 4))
            data ≈ gate_data
        end
    end

    @testset "Test gate addition with decomposition" begin
        qasm_str = """OPENQASM 2.0;
                    include "qelib1.inc";
                    qreg q[3];
                    h q[0];
                    cx q[0],q[1];
                    cx q[1],q[2];"""

        circ = load_qasm_as_circuit(qasm_str)
        tng = convert_qiskit_circ_to_network(circ, InteractiveBackend(), decompose=true)

        @test begin
            length(tng.nodes) == 5
        end
    end

    @testset "Test transpilation" begin
        qasm_str = """OPENQASM 2.0;
                    include "qelib1.inc";
                    qreg q[3];
                    h q[0];
                    cx q[0],q[2];"""

        circ = load_qasm_as_circuit(qasm_str)
        tng = convert_qiskit_circ_to_network(circ, InteractiveBackend(), transpile=true)

        @test begin
            tng.qubit_ordering == [2, 1, 3]
        end
    end

    @testset "Test add input/output functions" begin

        tng = TensorNetworkCircuit(3)

        # test that add_input function adds nodes
        add_input!(tng, "000")
        @test length(tng.nodes) == 3

        # test that add_output function adds nodes
        add_output!(tng, "000")
        @test length(tng.nodes) == 6

        # test that add_input function does not add nodes if already present
        add_input!(tng, "000")
        @test length(tng.nodes) == 6

        # test that add_output function does not add nodes if already present
        add_output!(tng, "000")
        @test length(tng.nodes) == 6
    end
end
end


text:
using GeometricIntegrators.BasisFunctions
using GeometricIntegrators.CommonFunctions
using Test

include("vandermonde_tests.jl")
include("lagrange_basis_tests.jl")
include("lagrange_polynomials_tests.jl")


text:
#=
mutable struct fmpr
  _man::Int
  _exp::Int

  function fmpr()
    z = new()
    ccall((:fmpr_init, libarb), Nothing, (Ref{fmpr}, ), z)
    finalizer(_fmpr_clear_fn, z)
    return z
  end

  function fmpr(x::Ref{arf_struct})
    z = new()
    ccall((:fmpr_init, libarb), Nothing, (Ref{fmpr}, ), z)
    ccall((:arf_get_fmpr, libarb), Nothing, (Ref{fmpr}, Ref{arf_struct}), z, x)
    finalizer(_fmpr_clear_fn, z)
    return z
  end

  function fmpr(x::arf_struct)
    z = new()
    ccall((:fmpr_init, libarb), Nothing, (Ref{fmpr}, ), z)
    ccall((:arf_get_fmpr, libarb), Nothing, (Ref{fmpr}, Ref{arf_struct}), z, x)
    finalizer(_fmpr_clear_fn, z)
    return z
  end
end

function _fmpr_clear_fn(x::fmpr)
  ccall((:fmpr_clear, libarb), Nothing, (Ref{fmpr}, ), x)
end

function fmpq(x::fmpr)
  z = fmpq()
  ccall((:fmpr_get_fmpq, libarb), Nothing, (Ref{fmpq}, Ref{fmpr}), z, x)
  return z
end
=#

mutable struct cfrac
  coeff::Ptr{fmpz}
  n::Int
  l::Int # real length

  function cfrac(x::Int)
    z = new()
    z.coeff = ccall((:_fmpz_vec_init, libflint), Ptr{fmpz}, (Int, ), x)
    z.n = x
    finalizer(_cfrac_clear_fn, z)
    return z
  end
end

function _cfrac_clear_fn(x::cfrac)
  ccall((:_fmpz_vec_clear, libflint), Nothing, (Ptr{fmpz}, Int), x.coeff, x.n)
end

function show(io::IO, x::cfrac)
  print(io, "[ ")
  for i in 1:x.l
    print(io, unsafe_load(x.coeff, i))
    print(io, " ,")
  end
  print(io, "]")
end

# THIS LEAKS MEMORY
function cfrac(x::fmpq, y::Int)
  r = fmpq()
  z = cfrac(y)
  l = ccall((:fmpq_get_cfrac, libflint), Int, (Ptr{fmpz}, Ref{fmpq}, Ref{fmpq}, Int), z.coeff, r, x, y)
  z.l = l
  return z, r
end

function fmpq(x::cfrac)
  z = fmpq()
  ccall((:fmpq_set_cfrac, libflint), Nothing, (Ref{fmpq}, Ptr{fmpz}, Int), z, x.coeff, x.l)
  return z
end

function fmpq(x::cfrac, y::Int)
  z = fmpq()
  ccall((:fmpq_set_cfrac, libflint), Nothing, (Ref{fmpq}, Ptr{fmpz}, Int), z, x.coeff, y)
  return z
end


text:
function get_data_type()

    data_type_prompt = "Step 1: Choose your input data type: "

    # Select between raster and network
    o = Observable("Raster")
    data_type = vbox(node(:div, data_type_prompt,
                          attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                     vskip(0.5em),
                     dropdown(["Raster", "Network"]))

    s = Scope()
    s.dom = data_type
    s["value"] = o

    s
end

function get_mod_mode_network()

    mod_mode_prompt = "Step 2: Choose a Modelling Mode: "
    o = Observable("Pairwise")
    mod_mode = vbox(node(:div, mod_mode_prompt,
                          attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                    vskip(0.5em),
                    dropdown(["Pairwise", "Advanced"]))

    s = Scope()
    s.dom = mod_mode
    s["value"] = o

    s
end

function get_mod_mode_raster()

    mod_mode_prompt = "Step 2: Choose a Modelling Mode: "
    o = Observable("Pairwise")
    mod_mode = vbox(node(:div, mod_mode_prompt,
                          attributes = Dict(:style => "margin-top: 12px")) |> class"b",
                    vskip(0.5em),
                    dropdown(["Pairwise",
                              "Advanced",
                              "One To All",
                              "All To One"], value = o))

    s = Scope()
    s.dom = mod_mode
    s["value"] = o

    s
end
function run_button()

    ob = Observable(rand())
    run = button("Run", value = ob)

	run, ob
end

const TESTPATH = joinpath(dirname(pathof(Circuitscape)), "..", "test")
function runtests_button(w)
    ob = Observable(rand())
    tests = button("Run tests!", value = ob)

    on(ob) do x
        p = pwd()
        cd(TESTPATH)
        try
            Circuitscape.runtests()
            Blink.@js_ w alert("Circuitscape tests passed!")
        catch
            Blink.@js_ w alert("Circuitscape tests failed!")
        end
        cd(p)
    end

    tests
end


text:
using Test
using st16SseqJuliaTools
using Suppressor

@testset "extract_properly_matching_rRNA_names" begin
    names_file, names_fo = mktemp()
    close(names_fo)
    @suppress_err begin
        @testset "Case: BOUNDARIES" begin
            hmm_str = "100366               -          NAME -              664     932     276       1     280       1     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
            @testset "Case: IN" begin
                input_file, input_fo = mktemp()
                println(input_fo, "#COMMENT")
                println(input_fo, "# COMMENT2")
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "0:500", "600:950", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:500", "600:950", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:280", "600:950", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:280", "664:950", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:280", "664:932", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:280", "664:932", 280, input_file)
            end

            @testset "Case: OUT" begin
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(String[]) == extract_rRNA_names(names_file, "2:500", "600:950", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(String[]) == extract_rRNA_names(names_file, "1:279", "600:950", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(String[]) == extract_rRNA_names(names_file, "1:280", "665:932", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(String[]) == extract_rRNA_names(names_file, "1:280", "664:931", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(String[]) == extract_rRNA_names(names_file, "1:280", "664:932", 281, input_file)
            end

            @testset "Case: MIX" begin
                hmm_str2 = "100365               -          NAME2 -              500     700     276       1     300       150     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100366", "100365"]) == extract_rRNA_names(names_file, "1:300", "500:932", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "1:300", "500:700", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "150:300", "500:1000", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:300", "500:932", 270, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:280", "600:932", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str)
                close(input_fo)
                @test Set(["100366"]) == extract_rRNA_names(names_file, "1:280", "600:932", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "150:300", "500:1000", 150, input_file)

                input_file, input_fo = mktemp()
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                println(input_fo, hmm_str)
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100366", "100365"]) == extract_rRNA_names(names_file, "1:300", "500:932", 150, input_file)
            end

            @testset "Case: POS MIX" begin
                hmm_str2 = "100365               -          NAME2 -              700     500     276       1     300       150     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "150:300", "500:700", 150, input_file)

                hmm_str2 = "100365               -          NAME2 -              500     700     276       1     150       300     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "150:300", "500:700", 150, input_file)

                hmm_str2 = "100365               -          NAME2 -              700     500     276       1     150       300     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "150:300", "500:700", 150, input_file)

                hmm_str2 = "100365               -          NAME2 -              700     500     276       1     150       300     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "300:150", "500:700", 150, input_file)

                hmm_str2 = "100365               -          NAME2 -              700     500     276       1     150       300     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "300:150", "700:500", 150, input_file)

                hmm_str2 = "100365               -          NAME2 -              700     500     276       1     150       300     280    -     1.1e-83  279.3   8.9  reference=Bacillus_subtilis_:2949495-2951495(+)/135-2000_1591 amplicon=1..726 position=1..300"
                input_file, input_fo = mktemp()
                println(input_fo, hmm_str2)
                close(input_fo)
                @test Set(["100365"]) == extract_rRNA_names(names_file, "150:300", "700:500", 150, input_file)
            end
        end
    end
end


text:
@generated function LinearAlgebra.mul!(d::RecursiveVector{T,N}, S::SymmetricMatrix{T,N,L}, x::RecursiveVector{T,N}) where {T,N,L}
    q, qa = create_quote()
    extract_linear!(qa, L, :S)
    extract_linear!(qa, N, :x)

    dnr = Symbol(:d_, 1)
    push!(qa, :( $dnr = S[1] * x[1] ) )
    for c ∈ 2:N
        push!(qa, :( $dnr += S[ $(small_triangle(c)+1) ] * x[$c] ) )
    end

    for r ∈ 2:N
        dnr = Symbol(:d_, r)
        push!(qa, :( $dnr = S[ $(small_triangle(r)+1) ] * x[1] ) )
        for c ∈ 2:r-1
            push!(qa, :( $dnr += S[ $(small_triangle(r)+c) ] * x[$c] ) )
        end
        # push!(qa, :() )
        for c ∈ r:N
            push!(qa, :( $dnr += S[ $(small_triangle(c)+r) ] * x[$c] ) )
        end
    end
    insert_linear!(qa, N, :d)
    push!(q.args, :(d))
    q
end


text:
for a = 1:500
    for b = a:500
        for c = b:500
            if a + b + c == 1000 && a^2 + b^2 == c^2
                println("a:$a,b:$b,c:$c")
                println("a×b×c = $(a*b*c)")
            end
        end
    end
end


text:
"""
    module HelmholtzOperators
"""
module HelmholtzOperators

export Helmholtz

files = (
    "1D/HelmholtzOperators1D.jl",
    "2D/Symmetric/HelmholtzOperators2D.jl",
    # "3D/HelmholtzOperators3D.jl"
    )

using ..VectorFields
using ..Simulations
using ..Dispersions
using SparseArrays

import ..Symmetric, ..Unsymmetric
import LinearAlgebra: I

"""
    maxwell(simulation{1}; ky=0, kz=0) -> m
"""
struct Helmholtz{N,CLASS,TΣ,TSIM}
    A::SparseMatrixCSC{ComplexF64,Int}
    D²::SparseMatrixCSC{ComplexF64,Int}
    αεpFχ::SparseMatrixCSC{Complex{Float64},Int}
    laplacian::SparseMatrixCSC{ComplexF64,Int}
    Σs::TΣ
    αε::SparseMatrixCSC{ComplexF64,Int}
    Fχ::SparseMatrixCSC{ComplexF64,Int}
    Fχs::Vector{SparseMatrixCSC{ComplexF64,Int}}
    dFχdψr::Array{ComplexF64,3}
    dFχdψi::Array{ComplexF64,3}
    dFχdω::Array{ComplexF64,3}
    dFχdϕ::Array{ComplexF64,3}
    simulation::TSIM
    ka::Float64
    kb::Float64
    kc::Float64
end

################################################################################
# susceptabilities and jacobians

# linear
@inline function helmholtz_susceptability!(h::Helmholtz,ω::Number)
    sim = h.sim
    @fastmath @inbounds @simd for i ∈ eachindex(h.Fχ.nzval)
        χ = susceptability(sim.χ[i], ω)
        h.Fχ.nzval[i] = sim.F[i]*χ
    end
    @fastmath @inbounds @simd for i ∈ eachindex(h.αε.nzval) h.αεpFχ.nzval[i] = h.αε.nzval[i] + h.Fχ.nzval[i] end
    return nothing
end

# nonlinear + local jacobian
@inline function helmholtz_susceptability!(h::Helmholtz, ω, ωs::Vector, ψs::ScalarField)
    sim = h.sim
    foreach(d->susceptability(d.χ,ω,ωs,ψs), sim.dispersive_domains)
    @inbounds for μ ∈ eachindex(h.Fχs)
        @inbounds for i ∈ eachindex(h.Fχ.nzval)
            if typeof(sim.χ[i])<:TwoLevelSystem
                χ::TwoLevelSystem = sim.χ[i]
                μ==1 ? h.Fχ.nzval[i] = sim.F[i]*χ.chi[i] : nothing
                h.Fχs[μ].nzval[i] = sim.F[i]*χ.chis[i,μ]
                @inbounds @simd for ν ∈ eachindex(h.Fχs)
                    h.dFχdψr[i,μ,ν] = sim.F[i]*χ.dχdψr[i,μ,ν]
                    h.dFχdψi[i,μ,ν] = sim.F[i]*χ.dχdψi[i,μ,ν]
                    h.dFχdω[i,μ,ν] = sim.F[i]*χ.dχdω[i,μ,ν]
                    h.dFχdϕ[i,μ,ν] = sim.F[i]*χ.dχdϕ[i,μ,ν]
                end
            end
        end
    end
    @fastmath @inbounds @simd for i ∈ eachindex(h.αε.nzval) h.αεpFχ.nzval[i] = h.αε.nzval[i] + h.Fχ.nzval[i] end
    return nothing
end

################################################################################
# Pretty Printing

import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_VARIABLE
import ..PRINTED_COLOR_INSTRUCTION

function Base.show(io::IO,::Helmholtz{N,CLASS}) where {N,CLASS}
    print(io,N,"D ")
    CLASS<:Symmetric ? print(io,"Symmetric ") : nothing
    CLASS<:Unsymmetric ? print(io,"Unsymmetric ") : nothing
    printstyled(io,"Helmholtz ",color=PRINTED_COLOR_DARK)
    print(io,"Operator ")
    printstyled(io,"(call w/ args ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ω",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,", [",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ωs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ψs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,"]) -> ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"A",color=PRINTED_COLOR_VARIABLE)
end

function Base.getproperty(h::Helmholtz,sym::Symbol)
    if sym==:sim
        return getfield(h,:simulation)
    elseif sym==:ky
        return getfield(h,:kb)
    elseif sym==:kz
        return getfield(h,:kc)
    else
        return getfield(h,sym)
    end
end

function Base.propertynames(::Helmholtz{1}, private=false)
    if private
        return fieldnames(Helmholtz)
    else
        (:simulation,:ky,:kz)
    end
end

# load 1D, 2D, 3D
foreach(include,files)

end #module


text:
using BenchmarkTools
using Printf

benchmark = true

exampleFilename = "src/advent/day12/example.txt"
inputFilename = "src/advent/day12/input.txt"

#### Direction

struct Direction
    name::Char
    heading::Int64
    dx::Int64
    dy::Int64
end

const north = Direction('N', 0, 0, 1)
const east = Direction('E', 90, 1, 0)
const south = Direction('S', 180, 0, -1)
const west = Direction('W', 270, -1, 0)

const directions = [north, east, south, west]
const directionByName = Dict('N' => north, 'E' => east, 'S' => south, 'W' => west)
const directionByHeading = Dict(0 => north, 90 => east, 180 => south, 270 => west)

function turnRight(direction::Direction, degrees::Int64)::Direction
    heading = direction.heading + degrees
    while heading < 0
        heading += 360
    end
    heading = heading % 360
    directionByHeading[heading]
end

turnLeft(direction::Direction, degrees::Int64)::Direction = turnRight(direction, -degrees)

#### Runner

function run(filename, entity; trace = false)
    f = open(filename)
    while !eof(f)
        line = readline(f)
        c, value = line[1], parse(Int64, line[2:length(line)])
        if c == 'F'
            doF(entity, value)
        elseif c == 'L'
            doL(entity, value)
        elseif c == 'R'
            doR(entity, value)
        else
            doDirection(entity, directionByName[c], value)
        end
        trace && @printf("%s -> %s\n", line, toString(entity))
    end
    close(f)
end

#### Point

struct Point
    x::Int64
    y::Int64
end

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

toString(point::Point) = @sprintf("(%s, %s)", point.x, point.y)

manhattanDistance(point::Point) = abs(point.x) + abs(point.y)

#### part 1

## Vehicle

mutable struct Vehicle
    position::Point
    heading::Direction
end

doDirection(vehicle::Vehicle, direction::Direction, distance::Int64) = vehicle.position += (direction.dx * distance, direction.dy * distance)

doF(vehicle::Vehicle, distance::Int64) = doDirection(vehicle, vehicle.heading, distance)

function doR(vehicle::Vehicle, degrees::Int64)
    vehicle.heading = turnRight(vehicle.heading, degrees)
end

doL(vehicle::Vehicle, degrees::Int64) = doR(vehicle, -degrees)

toString(vehicle::Vehicle) = @sprintf("%s %s", toString(vehicle.position), vehicle.heading.name)

##

function part1(filename; trace=false)
    vehicle = Vehicle(Point(0, 0), east)
    run(filename, vehicle; trace = trace)
    manhattanDistance(vehicle.position)
end

@assert part1(exampleFilename) == 25

@show part1(inputFilename)
benchmark && @btime part1(inputFilename);

#### part 2

## Vehicle2

mutable struct Vehicle2
    position::Point
    waypoint::Point
    heading::Direction
end

doDirection(vehicle::Vehicle2, direction::Direction, distance::Int64) = vehicle.waypoint += (direction.dx * distance, direction.dy * distance)

doF(vehicle::Vehicle2, value::Int64) = vehicle.position += (vehicle.waypoint.x * value, vehicle.waypoint.y * value)

function doR(vehicle::Vehicle2, degrees::Int64)
    radians = Base.Math.deg2rad(degrees)
    vehicle.waypoint = Point(round(vehicle.waypoint.x * cos(radians) + vehicle.waypoint.y * sin(radians)), round(vehicle.waypoint.y * cos(radians) - vehicle.waypoint.x * sin(radians)))
end

doL(vehicle::Vehicle2, degrees::Int64) = doR(vehicle, -degrees)

toString(vehicle::Vehicle2) = @sprintf("%s %s %s", toString(vehicle.position), toString(vehicle.waypoint), vehicle.heading.name)

##

function part2(filename; trace=false)
    vehicle = Vehicle2(Point(0, 0), Point(10, 1), east)
    run(filename, vehicle; trace = trace)
    manhattanDistance(vehicle.position)
end

@assert part2(exampleFilename) == 286

@show part2(inputFilename)
benchmark && @btime part2(inputFilename);


text:
int_rules_6_5_10 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.5.10*(c+d*x)^m*(a+b*sech)^n =#
    @apply_utils Antiderivative((~u) ^ ~(m') * sech(~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * sech(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ([~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
    @apply_utils Antiderivative((~u) ^ ~(m') * csch(~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * csch(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ([~m, ~n], ~x) && (LinearQ([~u, ~v], ~x) && Not(LinearMatchQ([~u, ~v], ~x)))
end

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.5 Hyperbolic secant/6.5.10 (c+d x)^m (a+b sech)^n.jl


text:
@testset "io" begin
    @testset "corpus" begin
        db = SQLite.DB()
        tablename = "test"
        colname = "text"
        schema = Tables.Schema([colname], [String])
        SQLite.createtable!(db, tablename, schema)
        DBInterface.execute(db, """
            INSERT INTO $tablename ($colname) VALUES ("a"), ("b");
        """)

        crps = corpus(db, tablename, colname)
        @test isa(crps, Corpus)
        @test text(crps[1]) == "a" && text(crps[2]) == "b"
    end
end


text:
@doc showText.function
###Функции для работы с текстом
"""
Отображение привественного текста
"""

function showText(name::String)
    t = "Hello, "*name*"! Welcome to this beautiful world!"
    return t
end

@doc showFirst.function

"""
Отображение первого символа в строке
"""

function showFirst(name::String)
    return name[1]

end


text:
# # Prepare Greenberger–Horne–Zeilinger state with Quantum Circuit

# First, you have to use this package in Julia.

using Yao, Yao.Blocks

# Then let's define the oracle, it is a function of the number of qubits.
# The circuit looks like this:

# ![ghz](../assets/figures/ghz4.png)

n = 4
circuit(n) = chain(
    n,
    put(1=>X),
    repeat(H, 2:n),
    control(2, 1=>X),
    control(4, 3=>X),
    control(3, 1=>X),
    control(4, 3=>X),
    repeat(H, 1:n),
)

# Let me explain what happens here. Firstly, we have an `X` gate which is applied to the first
# qubit. We need decide how we calculate this numerically, `Yao` offers serveral different approach
# to this. The simplest one is to use `put(n, ibit=>gate)` to apply a gate on the register.
# The first argument `n` means the number of qubits, it can be lazy evaluated.

put(n, 1=>X) == put(1=>X)(n)

# If you wanted to apply a two qubit gate,
put(n, (2,1)=>CNOT)

# However, this kind of general apply is not as efficient as the following statement
mat(put(n, (2,1)=>CNOT)) ≈ mat(control(n, 2, 1=>X))

# This means there is a `X` gate on the first qubit that is controled by the second qubit.
# `Yao.jl` providea a simple API `mat` to obtain the matrix representation of a block SUPER efficiently.
# This distinct feature helps users debug their quantum programs easily, and is equally useful in time evolution and ground state solving problems.

# For a multi-controlled gate like Toffoli gate, the construction is quite intuitive
control(n, (2, 1), 3=>X)
# Do you know how to construct a general multi-control, multi-qubit gate? Just have a guess and try it out!

# In the begin and end, we need to apply `H` gate to all lines, you can do it by `repeat`,
# For some specific types of gates such as `X`, `Y` and `Z`, applying multiple of them can be as efficient as applying single gate.

# The whole circuit is a chained structure of the above blocks. And we actually store a quantum
# circuit in a tree structure.

circuit

# After we have an circuit, we can construct a quantum register, and
# input it into the oracle. You will then receive this register after
# processing it.

r = apply!(register(bit"0000"), circuit(4))

# Let's check the output:

statevec(r)

# We have a GHZ state here, try to measure the first qubit

measure(r, nshot=1000)

# ![GHZ](../assets/figures/GHZ.png)

# GHZ state will collapse to ``|0000\rangle`` or ``|1111\rangle`` due to entanglement!


text:
@testset "CountCompare2Count: $(fct_type), $(array_size) items, $(T)" for fct_type in ["vector of variables", "vector affine function"], array_size in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(CountModel{T}())
    model = COIB.CountCompare2Count{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Count{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.CountCompare,
    )

    x_count, _ = MOI.add_constrained_variable(model, MOI.Integer())

    if T == Int
        x_array_1, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:array_size])
        x_array_2, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:array_size])
    elseif T == Float64
        x_array_1 = MOI.add_variables(model, array_size)
        x_array_2 = MOI.add_variables(model, array_size)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x_count, x_array_1..., x_array_2...])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x_count, x_array_1..., x_array_2...])
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct, CP.CountCompare(array_size))

    @test MOI.is_valid(model, x_count)
    for i in 1:array_size
        @test MOI.is_valid(model, x_array_1[i])
        @test MOI.is_valid(model, x_array_2[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.CountCompare}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Count{MOI.EqualTo{T}}) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == Tuple{Type}[]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 0
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}()) == 1

        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}()) == [bridge.con]
    end

    @testset "Sum" begin
        @test MOI.is_valid(model, bridge.con)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con)
        @test length(f.terms) == 2 * array_size + 1
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con) == CP.Count(array_size, MOI.EqualTo(zero(T)))

        t = f.terms[1]
        @test t.output_index == 1
        @test t.scalar_term.coefficient === one(T)
        @test t.scalar_term.variable == x_count

        for i in 1:array_size
            t1 = f.terms[2 + 2 * (i - 1)]
            @test t1.output_index == 1 + i
            @test t1.scalar_term.coefficient === one(T)
            @test t1.scalar_term.variable == x_array_1[i]

            t2 = f.terms[3 + 2 * (i - 1)]
            @test t2.output_index == 1 + i
            @test t2.scalar_term.coefficient === -one(T)
            @test t2.scalar_term.variable == x_array_2[i]
        end
    end
end


text:
module ExtendedWrapper

using JavaCall

const AbstractKalmanFilter = @jimport org.hipparchus.filtering.kalman.AbstractKalmanFilter
const Class = @jimport java.lang.Class
const ExtendedKalmanFilter = @jimport org.hipparchus.filtering.kalman.extended.ExtendedKalmanFilter
const JString = @jimport java.lang.String
const MatrixDecomposer = @jimport org.hipparchus.linear.MatrixDecomposer
const Measurement = @jimport org.hipparchus.filtering.kalman.Measurement
const NonLinearEvolution = @jimport org.hipparchus.filtering.kalman.extended.NonLinearEvolution
const NonLinearProcess = @jimport org.hipparchus.filtering.kalman.extended.NonLinearProcess
const Object = @jimport java.lang.Object
const ProcessEstimate = @jimport org.hipparchus.filtering.kalman.ProcessEstimate
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const RealVector = @jimport org.hipparchus.linear.RealVector

export ExtendedKalmanFilter
export NonLinearEvolution
export NonLinearProcess
export equals
export estimation_step
export get_class
export get_corrected
export get_current_state
export get_current_time
export get_evolution
export get_innovation
export get_measurement_jacobian
export get_predicted
export get_process_noise_matrix
export get_state_transition_matrix
export hash_code
export notify
export notify_all
export to_string
export wait

include("extended_kalman_filter.jl")
include("non_linear_evolution.jl")
include("non_linear_process.jl")

end


text:
include("jset.jl")
include("pgmfile.jl")

h = 400; w = 800; m = Array(Int64, h, w);
c0 = -0.8+0.16im;
pgm_name = "julia.pgm";

t0 = time();
for y=1:h, x=1:w
    c = complex((x-w/2)/(w/2), (y-h/2)/(w/2))
    m[y,x] = juliaset(c, c0, 256)
end
t1 = time();

create_pgmfile(m, pgm_name);
print("Written $pgm_name\nFinished in $(t1-t0) seconds.\n");

max_stars_repo_path: Module 3/Chapter02/jmain.jl

max_issues_repo_path: Module 3/Chapter02/jmain.jl

max_forks_repo_path: Module 3/Chapter02/jmain.jl


text:
# Tx dimensions
rTx = 10.4
# Low Pass filter corner frequencies
lowpassfcs = [300000.0, 450000.0]
# halt noise
LM_noise = vec([6.4087e+00 4.0463e+00 4.2526e+00 3.0585e+00 2.3754e+00 2.0971e+00 1.7919e+00 1.6001e+00 1.3299e+00 9.6707e-01 7.7027e-01 7.4084e-01 8.1148e-01 6.4811e-01 6.4059e-01 6.3400e-01 5.9044e-01])
HM_noise = vec([1.0652e-02 1.0436e-02 1.0185e-02 1.0033e-02 9.9405e-03 9.8320e-03 9.6964e-03 9.4945e-03 9.2693e-03 8.9339e-03 8.5681e-03 8.2023e-03     8.4102e-03 8.9711e-03 8.3685e-03 8.1345e-03 7.8504e-03 7.4964e-03 7.0504e-03 6.7447e-03 6.7160e-03 5.8913e-03 5.2938e-03 4.7887e-03 4.3874e-03])
# LM settings
LM_ramp = [
-0.0008000000   0.0000000000
-0.0007647300   0.0634310000
-0.0006281800   0.2497100000
-0.0004349700   0.4745300000
-0.0000921970   0.8904400000
-0.0000020929   1.0000000000
 0.0000000663   1.0000000000
 0.0000006956   0.9732500000
 0.0000021480   0.7986500000
 0.0000037941   0.5317200000
 0.0000056822   0.2765300000
 0.0000071829   0.1506200000
 0.0000085385   0.0750730000
 0.0000101360   0.0314230000
 0.0000119760   0.0079197000
 0.0000131380   0.0000000000
 0.0010181819   0.0000000000]

LM_times = vec(10 .^mean(log10.([
            0.0000186300    0.0000232000
            0.0000236300    0.0000292000
            0.0000296300    0.0000372000
            0.0000376300    0.0000472000
            0.0000476300    0.0000602000
            0.0000606300    0.0000762000
            0.0000766300    0.0000962000
            0.0000966300    0.0001212000
            0.0001216300    0.0001522000
            0.0001526300    0.0001922000
            0.0001926300    0.0002432000
            0.0002436300    0.0003062000
            0.0003066300    0.0003872000
            0.0003876300    0.0004882000
            0.0004886300    0.0006152000
            0.0006156300    0.0007762000
            0.0007766300    0.0009782000
]), dims=2))

# HM settings
HM_ramp = [
-0.005       0.0
-0.0048681   0.4938
-0.0047076   0.92737
-0.0046904   0.95949
 0.0         1.0
 7.2018e-6   0.9792
 9.9971e-5   0.69067
 0.00028748  0.065734
 0.00030229  0.015661
 0.00030504  0.0070146
 0.00030705  0.0028248
 0.00031197  0.0
 0.015       0.0]

HM_times = vec(10 .^mean(log10.([
            0.0003976300    0.0004102000
            0.0004106300    0.0004262000
            0.0004266300    0.0004462000
            0.0004466300    0.0004712000
            0.0004716300    0.0005022000
            0.0005026300    0.0005422000
            0.0005426300    0.0005932000
            0.0005936300    0.0006562000
            0.0006566300    0.0007372000
            0.0007376300    0.0008382000
            0.0008386300    0.0009652000
            0.0009656300    0.0011262000
            0.0011266300    0.0013282000
            0.0013286300    0.0015832000
            0.0015836300    0.0019052000
            0.0019056300    0.0023112000
            0.0023116300    0.0028222000
            0.0028226300    0.0034672000
            0.0034676300    0.0042802000
            0.0042806300    0.0053052000
            0.0053056300    0.0065972000
            0.0065976300    0.0082252000
            0.0082256300    0.0102792000
            0.0102796300    0.0125052000
            0.0125056300    0.0148152000
]), dims=2))

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
###
### Mer specific specializations of src/biosequence/predicates.jl
###

Base.cmp(seq1::T, seq2::T) where {T<:AbstractMer} = cmp(encoded_data(seq1), encoded_data(seq2))
Base.:(==)(x::T, y::T) where {T<:AbstractMer} = encoded_data(x) == encoded_data(y)
Base.isless(x::T, y::T) where {T<:AbstractMer} = isless(encoded_data(x), encoded_data(y))


text:
#########################################
### Overloads for `TensorBoardLogger` ###
#########################################
# `tb_name` is used by `preprocess` to decide how a given `arg` should look
"""
    tb_name(args...)

Returns a `string` representing the name for `arg` or `args` in TensorBoard.

If `length(args) > 1`, `args` are joined together by `"/"`.
"""
tb_name(arg) = string(arg)
tb_name(stat::OnlineStat) = string(nameof(typeof(stat)))
tb_name(o::Skip) = "Skip($(o.b))"
tb_name(o::Thin) = "Thin($(o.b))"
tb_name(o::WindowStat) = "WindowStat($(o.window.b))"
tb_name(o::AutoCov, b::Int) = "AutoCov(lag=$b)/corr"

# Recursive impl
tb_name(s1::String, s2::String) = s1 * "/" * s2
tb_name(arg1, arg2) = tb_name(arg1) * "/" * tb_name(arg2)
tb_name(arg, args...) = tb_name(arg) * "/" * tb_name(args...)

function TBL.preprocess(name, stat::OnlineStat, data)
    if nobs(stat) > 0
        TBL.preprocess(tb_name(name, stat), value(stat), data)
    end
end

function TBL.preprocess(name, stat::Skip, data)
    return TBL.preprocess(tb_name(name, stat), stat.stat, data)
end

function TBL.preprocess(name, stat::Thin, data)
    return TBL.preprocess(tb_name(name, stat), stat.stat, data)
end

function TBL.preprocess(name, stat::WindowStat, data)
    return TBL.preprocess(tb_name(name, stat), value(stat), data)
end

function TBL.preprocess(name, stat::AutoCov, data)
    autocor = OnlineStats.autocor(stat)
    for b = 1:(stat.lag.b - 1)
        # `autocor[i]` corresponds to the lag of size `i - 1` and `autocor[1] = 1.0`
        bname = tb_name(stat, b)
        TBL.preprocess(tb_name(name, bname), autocor[b + 1], data)
    end
end

function TBL.preprocess(name, stat::Series, data)
    # Iterate through the stats and process those independently
    for s in stat.stats
        TBL.preprocess(name, s, data)
    end
end

function TBL.preprocess(name, hist::KHist, data)
    if nobs(hist) > 0
        # Creates a NORMALIZED histogram
        edges = OnlineStats.edges(hist)
        cnts = OnlineStats.counts(hist)
        TBL.preprocess(
            name, (edges, cnts ./ sum(cnts)), data
        )
    end
end

# Unlike the `preprocess` overload, this allows us to specify if we want to normalize
function TBL.log_histogram(
    logger::TBLogger, name::AbstractString, hist::OnlineStats.HistogramStat;
    step=nothing, normalize=false
)
    edges = edges(hist)
    cnts = Float64.(OnlineStats.counts(hist))
    if normalize
        return TBL.log_histogram(logger, name, (edges, cnts ./ sum(cnts)); step=step)
    else
        return TBL.log_histogram(logger, name, (edges, cnts); step=step)
    end
end


text:
"""
    Reification{S <: MOI.AbstractSet}(set::S)

``\\{(y, x) \\in \\{0, 1\\} \\times \\mathbb{R}^n | y = 1 \\iff x \\in set, y = 0 otherwise\\}``.

This set serves to find out whether a given constraint is satisfied.

The only possible values are 0 and 1 for the first variable of the set.
"""
struct Reification{S <: MOI.AbstractSet} <: MOI.AbstractVectorSet
    set::S
end

MOI.dimension(set::Reification{S}) where {S} = 1 + MOI.dimension(set.set)
copy(set::Reification{S}) where {S} = Reification(copy(set.set))
Base.:(==)(x::Reification{S}, y::Reification{S}) where {S} = x.set == y.set

"""
    Equivalence{S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet}(set1::S1,
                                                              set2::S2)

The logical equivalence operator ≡ or ⇔.

``\\{(x, y) \\in \\mathbb{R}^{a+b} | x \\in S1 \\iff y \\in S2\\}``.

The two constraints must be either satisfied or not satisfied at the same time.
More explicitly, if the first one is satisfied, then the second one is implied
to be satisfied too; if the second one is satisfied, then the first one is
implied.
"""
struct Equivalence{S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet} <:
       MOI.AbstractVectorSet
    set1::S1
    set2::S2
end

function MOI.dimension(set::Equivalence{S, T}) where {S, T}
    return MOI.dimension(set.set1) + MOI.dimension(set.set2)
end

function copy(set::Equivalence{S, T}) where {S, T}
    return Equivalence(copy(set.set1), copy(set.set2))
end

function Base.:(==)(x::Equivalence{S, T}, y::Equivalence{S, T}) where {S, T}
    return x.set1 == y.set1 && x.set2 == y.set2
end

"""
EquivalenceNot{S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet}(set1::S1,
                                                                set2::S2)

The logical equivalence operator ≡ or ⇔, with the second argument negated.

``\\{(x, y) \\in \\mathbb{R}^{a+b} | x \\in S1 \\iff y \\not\\in S2\\}``.
"""
struct EquivalenceNot{S1 <: MOI.AbstractSet, S2 <: MOI.AbstractSet} <:
       MOI.AbstractVectorSet
    set1::S1
    set2::S2
end

function MOI.dimension(set::EquivalenceNot{S, T}) where {S, T}
    return MOI.dimension(set.set1) + MOI.dimension(set.set2)
end

function copy(set::EquivalenceNot{S, T}) where {S, T}
    return EquivalenceNot(copy(set.set1), copy(set.set2))
end

function Base.:(==)(x::EquivalenceNot{S, T}, y::EquivalenceNot{S, T}) where {S, T}
    return x.set1 == y.set1 && x.set2 == y.set2
end

"""
    IfThenElse{
        Condition <: MOI.AbstractSet,
        TrueConstraint <: MOI.AbstractSet,
        FalseConstraint <: MOI.AbstractSet
    }(condition::Condition, true_constraint::TrueConstraint,
      false_constraint::FalseConstraint)

The ternary operator.

If the `condition` is satisfied, then the first constraint (of type
`TrueConstraint`) will be implied. Otherwise, the second constraint
(of type `FalseConstraint`) will be implied.

``\\{(x, y, z) \\in \\mathbb{R}^(a+b+c) | y \\in TrueConstraint \\iff x \\in set, z \\in FalseConstraint otherwise\\}``.
"""
struct IfThenElse{
    Condition <: MOI.AbstractSet,
    TrueConstraint <: MOI.AbstractSet,
    FalseConstraint <: MOI.AbstractSet,
} <: MOI.AbstractVectorSet
    condition::Condition
    true_constraint::TrueConstraint
    false_constraint::FalseConstraint
end

function MOI.dimension(set::IfThenElse{S, T, U}) where {S, T, U}
    return MOI.dimension(set.condition) +
           MOI.dimension(set.true_constraint) +
           MOI.dimension(set.false_constraint)
end

function copy(set::IfThenElse{S, T, U}) where {S, T, U}
    return IfThenElse(
        copy(set.condition),
        copy(set.true_constraint),
        copy(set.false_constraint),
    )
end

function Base.:(==)(x::IfThenElse{S, T, U}, y::IfThenElse{S, T, U}) where {S, T, U}
    return x.condition == y.condition && x.true_constraint == y.true_constraint && x.false_constraint == y.false_constraint
end

"""
    Implication{
        Antecedent <: MOI.AbstractSet,
        Consequent <: MOI.AbstractSet
    }(antecedent::Antecedent, consequent::Consequent)

The logical implication operator ⇒.

If the `antecedent` is satisfied, then the `consequent` will be implied to be
satisfied. Otherwise, nothing is implied on the truth value of `consequent`.

``\\{(x, y) \\in \\mathbb{R}^a \\times \\mathbb{R}^b | y \\in Consequent if x \\in Antecedent\\}``.

Also called `if_then`, material implication, or material conditional.
"""
struct Implication{Antecedent <: MOI.AbstractSet, Consequent <: MOI.AbstractSet} <:
       MOI.AbstractVectorSet
    antecedent::Antecedent
    consequent::Consequent
end

function MOI.dimension(set::Implication{S, T}) where {S, T}
    return MOI.dimension(set.antecedent) + MOI.dimension(set.consequent)
end

function copy(set::Implication{S, T}) where {S, T}
    return Implication(copy(set.antecedent), copy(set.consequent))
end

function Base.:(==)(x::Implication{S, T}, y::Implication{S, T}) where {S, T}
    return x.antecedent == y.antecedent && x.consequent == y.consequent
end

"""
    Conjunction{Ts}(constraints::Ts)

The logical conjunction operator ∧ (AND): all the constraints in the
conjunction must be satisfied.

``\\{(x, y\\dots) \\in \\mathbb{R}^a \\times \\mathbb{R}^b\\dots | x \\in \\mathbb{S_1} \\land y \\in \\mathbb{S_2} \\dots \\}``.
"""
struct Conjunction{Ts} <: MOI.AbstractVectorSet where {Ts <: Tuple}
    constraints::Ts
end

# Currently, no varargs for parametric types... For instance, see:
# https://discourse.julialang.org/t/user-defined-variadic-parametric-types/25487/4

function MOI.dimension(set::Conjunction{Ts}) where {Ts}
    return sum(MOI.dimension(s) for s in set.constraints)
end

function copy(set::Conjunction{Ts}) where {Ts}
    return Conjunction(deepcopy(set.constraints))
end

function Base.:(==)(x::Conjunction{Ts}, y::Conjunction{Ts}) where {Ts}
    return x.constraints == y.constraints
end

"""
    Disjunction{Ts}(constraints::Ts)

The logical disjunction operator ∨ (OR): at least one of the constraints in
the disjunction must be satisfied.

``\\{(x, y\\dots) \\in \\mathbb{R}^a \\times \\mathbb{R}^b\\dots | x \\in \\mathbb{S_1} \\lor y \\in \\mathbb{S_2} \\dots \\}``.
"""
struct Disjunction{Ts} <: MOI.AbstractVectorSet where {Ts <: Tuple}
    constraints::Ts
end

# Currently, no varargs for parametric types... For instance, see:
# https://discourse.julialang.org/t/user-defined-variadic-parametric-types/25487/4

function MOI.dimension(set::Disjunction{Ts}) where {Ts}
    return sum(MOI.dimension(s) for s in set.constraints)
end

function copy(set::Disjunction{Ts}) where {Ts}
    return Disjunction(deepcopy(set.constraints))
end

function Base.:(==)(x::Disjunction{Ts}, y::Disjunction{Ts}) where {Ts}
    return x.constraints == y.constraints
end

"""
    ExclusiveDisjunction{Ts}(constraints::Ts)

The logical exclusive disjunction operator ⊻ (XOR): an odd number of
constraints in the exclusive disjunction must be satisfied.

``\\{(x, y\\dots) \\in \\mathbb{R}^a \\times \\mathbb{R}^b\\dots | x \\in \\mathbb{S_1} \\veebar y \\in \\mathbb{S_2} \\dots \\}``.
"""
struct ExclusiveDisjunction{Ts} <: MOI.AbstractVectorSet where {Ts <: Tuple}
    constraints::Ts
end

# Currently, no varargs for parametric types... For instance, see:
# https://discourse.julialang.org/t/user-defined-variadic-parametric-types/25487/4

function MOI.dimension(set::ExclusiveDisjunction{Ts}) where {Ts}
    return sum(MOI.dimension(s) for s in set.constraints)
end

function copy(set::ExclusiveDisjunction{Ts}) where {Ts}
    return ExclusiveDisjunction(deepcopy(set.constraints))
end

function Base.:(==)(x::ExclusiveDisjunction{Ts}, y::ExclusiveDisjunction{Ts}) where {Ts}
    return x.constraints == y.constraints
end

"""
    Negation{S <: MOI.AbstractSet}(set::S)

The logical negation operator ¬ (NOT).

``\\{x \\in \\times \\mathbb{R}^n | x \\not\\in set\\}``.
"""
struct Negation{S <: MOI.AbstractSet} <: MOI.AbstractVectorSet
    set::S
end

MOI.dimension(set::Negation{S}) where {S} = MOI.dimension(set.set)
copy(set::Negation{S}) where {S} = Negation(copy(set.set))
Base.:(==)(x::Negation{S}, y::Negation{S}) where {S} = x.set == y.set

"""
    True()

A constraint that is always true.

It is only useful with reification-like constraints.
"""
struct True <: MOI.AbstractVectorSet end

MOI.dimension(set::True) = 0
Base.:(==)(::True, ::True) = true

"""
    False()

A constraint that is always false.

It is only useful with reification-like constraints.
"""
struct False <: MOI.AbstractVectorSet end

MOI.dimension(set::False) = 0
Base.:(==)(::False, ::False) = true

# isbits types, nothing to copy
function copy(set::Union{True, False})
    return set
end


text:
module TestDoctest

using Documenter
using Test
using LiterateTest

@testset "/docs" begin
    doctest(LiterateTest; manual = true)
end

@testset "/test/doctests" begin
    doctest(joinpath((@__DIR__), "doctests"), Module[])
end

end  # module


text:
# add Interact   # pkg mode

using Interact, Blink

ddn = dropdown(["Getting groceries",
                "Visiting my therapist",
                "Getting a haircut",
                "Paying the energy bill",
                "Blog about workspace management"])

w = Window()
ui = dom"div"(ddn)
body!(w, ui)


text:
# Procedure derived from James Slagle's MIT Ph.D. Thesis
#
# ------------------------------------------------------------------------------
# Integral struct
# ------------------------------------------------------------------------------

struct Integral{I <: Union{SymbolicUtils.Symbolic, Number}, T} <: SymbolicUtils.Symbolic{T}
# struct Integral{I <: SymbolicUtils.Symbolic, T} <: SymbolicUtils.Symbolic{T}
    integrand::I
    variable::SymbolicUtils.Sym
    antiderivative::Ref{Union{SymbolicUtils.Symbolic, Nothing}}

    # Not currently used for anything, but will eventually support definite integrals.
    limits::Union{NTuple{2, T}, Nothing}
    value::Ref{Union{SymbolicUtils.Symbolic, Nothing}}

    equivalent_forms::Vector{SymbolicUtils.Symbolic}
    ef_integrals::Vector{Vector{Integral}}
    parents::Vector{Integral}

    # These define the "character" of the integral. They are not evaluated
    # immediately for performance, but the evaluation must be done before the
    # goal is placed into the heuristic list.
    depth::Int
    size::Int
    rational::Bool
    algebraic::Bool
    exponent_base::Union{Number, SymbolicUtils.Symbolic, Nothing}
    rat_sin_cos::Bool
    elf_sin_cos::Bool
    elf_sec_tan::Bool
    elf_csc_cot::Bool
    elf_trig::Bool
    partial::Bool

    function Integral(integrand, variable)

        return new{typeof(integrand), Float64}(
            integrand,
            variable,
            Ref{Union{SymbolicUtils.Symbolic, Nothing}}(nothing),
            nothing,
            Ref{Union{SymbolicUtils.Symbolic, Nothing}}(nothing),
            SymbolicUtils.Symbolic[],
            Vector{Integral}[],
            Integral[])
    end
end

function Base.isless(a::Integral, b::Integral)
    return SymbolicUtils.:<ₑ(a.integrand, b.integrand)
end

function Base.show(io::IO, i::Integral)
    print(io, "Integral(", i.integrand, ", ", i.variable, ")")
end

isintegral(i::Integral) = true
isintegral(x) = false

issolved(i::Integral) = i.antiderivative[] !== nothing

function isneeded(i::Integral)
    issolved(i) && return false
    length(i.parents) == 0 && return true
    return all(isneeded.(i.parents))
end

function check_solution!(i::Integral)
    for (n, form_integrals) in enumerate(i.ef_integrals)
        all(issolved, form_integrals) || continue

        r = @rule(~i::isintegral => (~i).antiderivative[])
        sub = SymbolicUtils.Postwalk(SymbolicUtils.PassThrough(r))

        antiderivative = sub(i.equivalent_forms[n])
        set_solution!(i, antiderivative)
    end
end

function set_solution!(i::Integral, solution)
    i.antiderivative[] = solution

    for parent in i.parents
        check_solution!(parent)
    end

    return solution
end

function compute_character(i::Integral)
    return
end

# ------------------------------------------------------------------------------
# Standard forms
# ------------------------------------------------------------------------------

function _apply_standard_forms(ex, v)
    is_const_wrt_v(ex) = is_const(ex, v)

    standard_forms = [
        @rule(~c::is_const_wrt_v => ~c*v)

        @rule(exp(v) => exp(v))
        @rule((~c::is_const_wrt_v)^v => ~c^v / log(~c))

        @rule(log(v) => v*log(v) + -1*v)
        @rule(log(~c::is_const_wrt_v, v) => v*log(~c, v) - v / log(~c))

        @rule(sin(v) => -1*cos(v))
        @rule(cos(v) => sin(v))
        @rule(tan(v) => log(sec(v)))
        @rule(cot(v) => log(sin(v)))
        @rule(sec(v) => log(sec(v) + tan(v)))
        @rule(csc(v) => log(csc(v) - cot(v)))

        @rule(asin(v) => v*asin(v) + sqrt(1 - v^2))
        @rule(acos(v) => v*acos(v) - sqrt(1 - v^2))
        @rule(atan(v) => v*atan(v) - 1//2*log(1 + v^2))
        @rule(acot(v) => v*acot(v) + 1//2*log(1 + v^2))
        @rule(asec(v) => v*asec(v) - log(v + sqrt(v^2 - 1)))
        @rule(acsc(v) => v*acsc(v) + log(v + sqrt(v^2 - 1)))

        @rule(sec(v)^2 => tan(v))
        @rule(csc(v)^2 => -1*cot(v))

        @rule(v^-1 => log(v))
        # Need this extra rule because the next one won't match a linear v
        @rule(v => 1//2 * v^2)
        @rule(v^(~c::is_const_wrt_v) => 1//(~c + 1)*v^(~c + 1))

        @acrule(sec(v)*tan(v) => sec(v))
        @acrule(csc(v)*cot(v) => -1*csc(v))

        # TODO: For the following three rules, I need to ensure that m != +/- n.
        # I don't believe that there is a way to do this using the current rule
        # system.
        # Although, Mathematica doesn't impose this restriction either...
        @acrule(sin(~m*v) * cos(~n*v) => -1//(2*(~m - ~n))*cos((~m - ~n)*v) -
                                          1//(2*(~m + ~n))*cos((~m + ~n)*v))
        @rule(sin(~m*v) * sin(~n*v) => 1//(2*(~m - ~n))*sin((~m - ~n)*v) -
                                       1//(2*(~m + ~n))*sin((~m + ~n)*v))
        @rule(cos(~m*v) * cos(~n*v) => 1//(2*(~m - ~n))*sin((~m - ~n)*v) +
                                       1//(2*(~m + ~n))*sin((~m + ~n)*v))
    ]

    for form in standard_forms
        eval = form(ex)
        eval != nothing && return eval
    end

    return nothing
end
#
# ------------------------------------------------------------------------------
# Algorithmic transformations
# ------------------------------------------------------------------------------

function factor_constant(i)
    typeof(i.integrand) <: SymbolicUtils.Mul || return []

    const_term = i.integrand.coeff
    variable_term = one(i.integrand.coeff)

    for (expr, pow) in pairs(i.integrand.dict)
        if is_const(expr, i.variable)
            const_term *= expr^pow
        else
            variable_term *= expr^pow
        end
    end

    const_term === one(i.integrand.coeff) && return []

    new_integral = Integral(variable_term, i.variable)
    push!(i.equivalent_forms, const_term * new_integral)
    push!(i.ef_integrals, [new_integral])
    push!(new_integral.parents, i)

    return [new_integral]
end

function decompose(i)
    typeof(i.integrand) <: SymbolicUtils.Add || return []
    length(i.integrand.dict) + (i.integrand.coeff == 0 ? 0 : 1) > 1 || return []

    integrals = Integral[]
    integral_sum = zero(i.integrand.coeff)
    if i.integrand.coeff != 0
        const_int = Integral(i.integrand.coeff, i.variable)
        push!(const_int.parents, i)
        push!(integrals, const_int)
        integral_sum += const_int
    end

    for (expr, coeff) in pairs(i.integrand.dict)
        int = Integral(coeff * expr, i.variable)
        push!(int.parents, i)
        push!(integrals, int)
        integral_sum += int
    end

    push!(i.equivalent_forms, integral_sum)
    push!(i.ef_integrals, integrals)

    return integrals
end

function linear_sub(p)
    # TODO: I expect that this section will benefit from CSE. I'll leave it
    # alone for now.
end

# This doesn't quite have the same function as in the paper because it
# recursively expands polynomials, instead of only at the top level. Not sure
# if that's a problem yet...
function expand_transform(i::Integral)
    new_integrand = expand(i.integrand)
    isequal(new_integrand, i.integrand) && return []

    new_integral = Integral(new_integrand, i.variable)
    push!(new_integral.parents, i)
    push!(i.equivalent_forms, new_integral)
    push!(i.ef_integrals, [new_integral])
    return [new_integral]
end

# Again, I am not strictly following the paper here because simplify can do
# much more than just eliminate common factors between a numerator and
# denominator.
function combine_factors(i)
    new_integrand = simplify(i.integrand)
    isequal(new_integrand, i.integrand) && return []

    new_integral = Integral(new_integrand, i.variable)
    push!(new_integral.parents, i)
    push!(i.equivalent_forms, new_integral)
    push!(i.ef_integrals, [new_integral])
    return [new_integral]
end

function divide_polynomials(p)
    # TODO
end

function half_angle_identities(p)
    # TODO
end

algo_transforms = [
    factor_constant
    decompose
    # linear_sub
    expand_transform
    combine_factors
    # divide_polynomials
    # half_angle_identities
]

# ------------------------------------------------------------------------------
# Heuristic transformations
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# Integral solver
# ------------------------------------------------------------------------------

function solve_or_add_to_list!(list, integral)
    # Try for a direction solution before adding to the goal list
    new_ex = _apply_standard_forms(integral.integrand, integral.variable)
    if new_ex !== nothing
        set_solution!(integral, new_ex)
        return
    end

    # I am deviating from SAINT by adding the integral to the list without
    # checking if an algorithmic tranformation is applicable. I expect that
    # this will be fine because the algorithmic transformations will probably
    # reduce the complexity of the integral, and so they will be considered
    # first.
    enqueue!(list, integral)

    # Check for algorithmic transformations
    for algo in algo_transforms
        new_integrals = algo(integral)
        solve_or_add_to_list!.(Ref(list), new_integrals)
    end
end

function solve!(initial_integral::Integral)
    temp_list = Queue{Integral}()
    heuristic_list = PriorityQueue{Integral, Int}()

    solve_or_add_to_list!(temp_list, initial_integral)

    # TODO: I should track resource usage here too.
    while true
        issolved(initial_integral) && return true

        # Move problems from temp_list to heuristic_list

        # If no heuristic problems remain, then we failed to solve the integral
        isempty(heuristic_list) && return false

        # Take the cheapest problem from the heuristic_list
        current_goal = dequeue!(heuristic_list)
    end

    return issolved(initial_integral)
end

function integrate(ex, v)
    integral = Integral(ex, v)
    solve!(integral)
    return integral.antiderivative[]
end


text:
using DynamicPolynomials
using MultivariateMoments
using SemialgebraicSets
using BeurlingLasso
using Test

const MM = MultivariateMoments

@testset "poly" begin
    @polyvar x[1:2]

    # only valid for odd numbers
    function doublefactorial(n::Integer)
        if n <= 0
            1
        else
            prod(2*k+1 for k=0:div(n-1, 2))
        end
    end

    nmoments = 15

    φ = [sum(binomial(n, 2*k) * x[1]^(n-2*k) * x[2]^(2*k) * doublefactorial(2*k-1)
         for k=0:div(n, 2)) for n=0:nmoments-1]

    params = [[-0.5, 0.5], [0.6, 0.3]]
    weights = [0.7, 0.3]

    natoms = length(params)

    μ0 = AtomicMeasure(x, [WeightedDiracMeasure(params[i], weights[i]) for i=1:natoms])
    y = [MM.expectation(μ0, φ[n]) for n=1:nmoments]

    λ = 1e-5

    φ_func(u, n) = φ[n](x=>u)
    bounds = [[-1 1]; [-1 1]]

    prob = blasso(y, φ_func, λ, bounds)
    solver = sfwsolver(2, 1000)
    solve!(prob, solver)
end


text:
#=
  8 schools problem

  Cf ~/stan/eight_schools.stan
"""
           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu         7.88    0.21   5.31  -2.91    4.6   7.91  11.12  18.27  658.0    1.0
tau        6.83    0.22   5.63   0.19   2.64   5.45   9.73  21.12  647.0   1.01
eta[1]     0.42    0.02   0.91  -1.52  -0.17   0.44   1.06   2.19 1843.0   1.01
eta[2]     0.01    0.02   0.89  -1.73  -0.57   0.01   0.59   1.78 2078.0    1.0
eta[3]    -0.19    0.02   0.89  -1.89   -0.8   -0.2   0.38   1.56 1837.0    1.0
eta[4]    -0.02    0.02   0.88  -1.78   -0.6  -0.02   0.54   1.66 1816.0    1.0
eta[5]    -0.37    0.02   0.88  -2.03  -0.97  -0.39   0.21   1.45 1511.0    1.0
eta[6]    -0.21    0.02   0.87  -1.96   -0.8  -0.24   0.38   1.57 1846.0    1.0
eta[7]     0.33    0.02   0.89  -1.48  -0.24   0.34   0.94    2.0 2000.0    1.0
eta[8]     0.08    0.02   0.91   -1.7  -0.53   0.08   0.69   1.85 2429.0    1.0
theta[1]   11.6    0.21   8.15  -1.93   6.18  10.61  15.69  30.96 1506.0    1.0
theta[2]   7.99    0.14   6.43  -4.29   3.86   7.86  11.95  20.89 2090.0    1.0
theta[3]   6.08    0.19   7.82 -12.72   2.22   6.59  10.81  20.56 1613.0    1.0
theta[4]   7.55    0.14   6.78  -6.74   3.51   7.69  11.79  21.31 2186.0    1.0
theta[5]   4.69    0.16   6.43  -9.19   0.82   5.34   9.12   16.0 1647.0    1.0
theta[6]   6.15    0.15   6.61  -8.13   2.34   6.55  10.63   18.5 1873.0    1.0
theta[7]  10.77    0.18   6.87  -1.07    6.2   9.98  14.85  25.87 1500.0    1.0
theta[8]   8.58     0.2   8.17  -7.99   3.87   8.24  12.88  26.77 1709.0    1.0
lp__      -4.74    0.11   2.63 -10.66  -6.35  -4.48  -2.85  -0.26  552.0    1.0
"""

  See ~/blog/8_shools.blog
      ~/webppl/8_shools.wppl

=#

# Fixing plot warnings GKS: character ignored due to unicode error
ENV["GKS_ENCODING"] = "utf-8"
using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function eight_schools(ys, sigmas)

    n = length(ys)

    # In contrast to stan, Turing requires a prior distribution.
    mu ~ Normal(10.0,10.0);
    tau ~ truncated(Normal(10.0,10.0),0.0,120.0)

    eta   =  tzeros(n)
    theta =  tzeros(n)
    sig   =  tzeros(n)
    for i in 1:n
        eta[i] ~ Normal(0.0,1.0)
        theta[i] = mu + tau * eta[i]
        # sig[i] ~ truncated(Normal(10,10),0,100) # s.e. of effect estimates
        sig[i] = sigmas[i]

        ys[i] ~ Normal(theta[i], sig[i]) # avg_likelihood
    end

end

ys     = [28.0,  8.0, -3.0,  7.0, -1.0,  1.0, 18.0, 12.0]
sigmas = [15.0, 10.0, 16.0, 11.0,  9.0, 11.0, 10.0, 18.0]
model = eight_schools(ys,sigmas)

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, MH(), 10_000)
chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(20), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, PG(20), 1_000)

# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), 10_000)

# chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 10_000)

# chains = sample(model, NUTS(1000,0.65), 1_000)
# chains = sample(model, HMC(0.1,5), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(1000,0.65,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(10,0.65,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),HMC(0.1,5,:height)), 1_000)
# chains = sample(model, Gibbs(PG(10,:gender),HMC(0.1,5,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(1_000,0.65,:height)), 1_000)

display(chains)
display(plot(chains))

# gen = generated_quantities(model, chains)
# show_var_dist_pct(gen, 20)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
function elapsed_time_criterion(duration::Number)
    start = time()
    while true
        produce((time() - start) >= duration)
    end
end

function iterations_criterion(duration::Number)
    start = 0
    while true
        produce(start >= duration)
        start += 1
    end
end


text:
"""
Struct holding parameters for ODE solver
"""
mutable struct DiffEqParams{T}
    problem::ODEProblem             # ODEProblem which is simulated
    dt::T                           # time step size
    solver::AbstractODEAlgorithm    # solver to use to solve ODE
    solve_args::Dict                # dictionary holding kwargs for ODE solve command
end

mutable struct DiffEqEnv{T} <: AbstractEnv
    ode_params::DiffEqParams{T}
    ic_sampler::ICSampler
    # Parameters for ReinforcementLearning
    observation_space::Union{ContinuousSpace{T},MultiContinuousSpace{Vector{T}}}
    action_space::Union{ContinuousSpace{T},MultiContinuousSpace{Vector{T}}}
    # Observation & reward functions
    observation_fn::ObservationFunction     # obvervation o=f(s)
    reward_fn::RewardFunction               # reward function r(s,a,s')
    # Buffer for previous transition
    state::Union{T,Vector{T}}               # current state
    observation::Union{T,Vector{T}}         # current observation
    action::Union{Nothing,T,Vector{T}}      # last action
    reward::Union{Nothing,T}                # last scalar reward
    done::Bool                              # true if in terminal state
    steps::Int                              # counter for steps taken in episode
    t::T                                    # current time
end

"""
User-facing constructor for DiffEqEnv with IC sampling
"""
function DiffEqEnv(
    problem::ODEProblem,
    reward_fn::RewardFunction,
    n_actions::Int,
    dt::Real,
    s0_lb::Union{Real,Vector{<:Real}},
    s0_ub::Union{Real,Vector{<:Real}};
    #= Keyword arguments =#
    observation_fn::ObservationFunction=FullStateObservation(),
    o_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for observation space
    o_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for observation space
    a_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for action space
    a_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for action space
    solver::DiffEqBase.AbstractODEAlgorithm=Euler(),
    reltol::Real=1e-8,
    abstol::Real=1e-8, # TODO: add more kwargs for integrator
    T=Float64
    )
    # Save solver arguments
    solve_args = Dict{Symbol,Any}(
        :reltol => reltol,
        :abstol => abstol,
        :save_everystep => false,
        :save_start => false) # only output terminal value

    ode_params = DiffEqParams{T}(problem, dt, solver, solve_args)

    # Set inital state and observation
    ic_sampler = UniformSampler(s0_lb, s0_ub)
    s0 = T.(ic_sampler()) # Convert type
    o0 = T.(observation_fn(s0, nothing)) # initial observations

    n_states = length(s0)
    n_observations = length(o0)

    # Helper function to set action and observation spaces
    function set_space(lb, ub, dim)
        # Check if types match
        typeof(lb) == typeof(ub) ||
            throw(ArgumentError("$(typeof(lb)) != $(typeof(ub)), types must match"))

        # Set defaults if no bounds are provided
        if isnothing(lb)
            if dim == 1
                ub = T(1e38)
                lb = -ub
            else
                ub = T(1e38) * ones(T, dim)
                lb = - ub
            end
        end

        # Check if sizes match
        length(lb) == length(ub) ||
            throw(ArgumentError("$(size(lb)) != $(size(ub)), size must match"))
        length(lb) == dim ||
            throw(ArgumentError("$(size(lb)) != $(dim), size must match"))

        # Return ContinuousSpace or MultiContinuousSpace
        if lb isa Real
            return ContinuousSpace(T(lb), T(ub))
        elseif lb isa Vector{<:Real}
            return MultiContinuousSpace(T.(lb), T.(ub))
        end
    end

    # Set observation and action spaces
    observation_space = set_space(o_lb, o_ub, n_observations)
    action_space = set_space(a_lb, a_ub, n_actions)

    # Initialize buffer holding previous step
    state = s0
    observation = o0
    action = nothing
    reward = nothing
    done = false
    steps = 0
    t = problem.tspan[1]

    env = DiffEqEnv{T}(ode_params, ic_sampler,
            observation_space, action_space,
            observation_fn, reward_fn,
            state, observation, action, reward, done, steps, t)
    return env
end

"""
User-facing constructor for DiffEqEnv with constant ICs
"""
function DiffEqEnv(
    problem::ODEProblem,
    reward_fn::RewardFunction,
    n_actions::Int,
    dt::Real;
    #= Keyword arguments =#
    observation_fn::ObservationFunction=FullStateObservation(),
    o_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for observation space
    o_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for observation space
    a_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for action space
    a_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for action space
    solver::DiffEqBase.AbstractODEAlgorithm=Euler(),
    reltol::Real=1e-8,
    abstol::Real=1e-8, # TODO: add more kwargs for integrator
    T=Float64
    )

    #= Set bounds for ICs to same values
    This ensures that the same IC gets sampled by env.ic_sampler =#
    s0_lb = T.(problem.u0)
    s0_ub = s0_lb

    return DiffEqEnv(problem, reward_fn, n_actions, dt, s0_lb, s0_ub;
        observation_fn=observation_fn,
        o_ub=o_ub, o_lb=o_lb, a_ub=a_ub, a_lb=a_lb,
        solver=solver, reltol=reltol, abstol=abstol, T=T)

end

"""
RLBase interface for use with ReinforcementLearning.jl
"""
RLBase.get_actions(env::DiffEqEnv) = env.action_space
RLBase.get_state(env::DiffEqEnv) = env.observation # returns observed state, not markov state!
RLBase.get_reward(env::DiffEqEnv) = env.reward
RLBase.get_terminal(env::DiffEqEnv) = env.done

function RLBase.reset!(env::DiffEqEnv{T}) where {T <: Real}
    # Reset environment
    env.state = T.(env.ic_sampler())
    env.action = nothing
    env.observation = T.(env.observation_fn(env.state, env.action))
    env.reward = nothing
    env.done = false
    env.steps = 0
    env.t = env.ode_params.problem.tspan[1]

    return nothing
end

function (env::DiffEqEnv{T})(action) where {T <: Real}
    # type of action must fit with type of env.action_space
    action = T.(action)

    # unpack Array{T,1} with single value for use in ContinuousSpace{T}
    if length(action) == 1
        action = action[]
    end

    @assert action in env.action_space

    # Remake ODEProblem over new tspan
    t_end =  env.t + env.ode_params.dt
    tspan = (env.t, t_end)
    prob = remake(env.ode_params.problem;
                u0=env.state, tspan=tspan, p=action)

    # Integrate ODE
    if isadaptive(env.ode_params.solver)
        sol = solve(prob, env.ode_params.solver; env.ode_params.solve_args...)
    else # add timestep argument dt for fixed step-width solvers
        sol = solve(prob, env.ode_params.solver; env.ode_params.solve_args..., dt=t_end)
    end
    state_next = T.(sol.u[1]) # unpack Array{Array{T,1},1}

    # Update environment buffer
    env.observation = T.(env.observation_fn(state_next, action))
    env.action = T.(action)
    env.reward = env.reward_fn(env.state, action, state_next) # update before env.state!
    env.state = state_next
    env.t = t_end # update before env.done!
    env.done = env.t >= env.ode_params.problem.tspan[2]
    env.steps += 1

    return nothing
 end


text:
abstract type AbstractClique end  # Defines a clique object

struct Clique <: AbstractClique
    head::BitVector
    member::BitVector

    function Clique(head::AbstractVector{Bool}, member::AbstractVector{Bool})
        size(head) == size(member) || throw(DimensionMismatch("size mismatch"))
        return new(BitVector(head), BitVector(member))
    end
end

function Clique(member::AbstractVector{Bool}; allowsingle=false)
    head = falses(length(member))
    head[findfirst(member)] = true
    @assert sum(member) > 1
    return Clique(head, member)
end

head(c::Clique) = c.head
member(c::Clique) = c.member
sizeindex(c::Clique) = length(head(c))