text:
"""
    RectilinearGrid{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch} <: AbstractRectilinearGrid{FT, TX, TY, TZ}

A rectilinear grid with with either constant or varying grid spacings between cell centers and cell faces
in all directions. Grid elements of type `FT`, topology `{TX, TY, TZ}`, grid spacings of type `{FX, FY, FZ}`
and coordinates in each direction of type `{VX, VY, VZ}`.
"""

struct RectilinearGrid{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch} <: AbstractRectilinearGrid{FT, TX, TY, TZ, Arch}
        architecture::Arch
        Nx :: Int
        Ny :: Int
        Nz :: Int
        Hx :: Int
        Hy :: Int
        Hz :: Int
        Lx :: FT
        Ly :: FT
        Lz :: FT
      # All directions can be either regular (FX, FY, FZ) <: Number or stretched (FX, FY, FZ) <: AbstractVector
      Δxᶠᵃᵃ :: FX
      Δxᶜᵃᵃ :: FX
      xᶠᵃᵃ  :: VX
      xᶜᵃᵃ  :: VX
      Δyᵃᶠᵃ :: FY
      Δyᵃᶜᵃ :: FY
      yᵃᶠᵃ  :: VY
      yᵃᶜᵃ  :: VY
      Δzᵃᵃᶠ :: FZ
      Δzᵃᵃᶜ :: FZ
      zᵃᵃᶠ  :: VZ
      zᵃᵃᶜ  :: VZ
end

const XRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number}
const YRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Number}
const ZRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Number}
const HRegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number}
const  RegRectilinearGrid = RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number, <:Number}

"""
    RectilinearGrid([architecture = CPU(), FT = Float64];
                    size,
                    x = nothing,
                    y = nothing,
                    z = nothing,
                    halo = nothing,
                    extent = nothing,
                    topology = (Periodic, Periodic, Bounded))

Creates a `RectilinearGrid` with `size = (Nx, Ny, Nz)` grid points.

Positional arguments
=================

- `architecture`: Specifies whether the array of coordinates, interfaces, and spacings are stored
                  on the CPU or GPU. Default: `architecture = CPU()`.

- `FT` : Floating point precision, might be Float32 or Float64

Keyword arguments
=================

- `size` (required): A tuple prescribing the number of grid points in non-`Flat` directions.
                     `size` is a 3-tuple for 3D models, a 2-tuple for 2D models, and either a
                     scalar or 1-tuple for 1D models.

- `topology`: A 3-tuple `(TX, TY, TZ)` specifying the topology of the domain.
              `TX`, `TY`, and `TZ` specify whether the `x`-, `y`-, and `z` directions are
              `Periodic`, `Bounded`, or `Flat`. The topology `Flat` indicates that a model does
              not vary in those directions so that derivatives and interpolation are zero.
              The default is `topology = (Periodic, Periodic, Bounded)`.

- `extent`: A tuple prescribing the physical extent of the grid in non-`Flat` directions.
            All directions are contructed with regular grid spacing and the domain (in the
            case that no direction is `Flat`) is x ∈ (0, Lx), y ∈ (0, Ly), and z ∈ (-Lz, 0),
            which is most appropriate for oceanic applications with z = 0 usually being the
            ocean's surface.

- `x`, `y`, and `z`: Each of `x, y, z` are either (i) 2-tuples that specify the end points of the domain
                     in their respect directions (in which case scalar values may be used in `Flat`
                     directions), or (ii) arrays or functions of the corresponding indices `i`, `j`, or `k`
                     that specify the locations of cell faces in the `x`-, `y`-, or `z`-direction, respectively.
                     For example, to prescribe the cell faces in `z` we need to provide a function that takes
                     `k` as argument and retuns the location of the faces for indices `k = 1` through `k = Nz + 1`,
                     where `Nz` is the `size` of the stretched `z` dimension.

*Note*: _Either_ `extent`, or all of `x`, `y`, and `z` must be specified.

- `halo`: A tuple of integers that specifies the size of the halo region of cells surrounding
          the physical interior for each non-`Flat` direction.

The physical extent of the domain can be specified via `x`, `y`, and `z` keyword arguments
indicating the left and right endpoints of each dimensions, e.g. `x = (-π, π)` or via
the `extent` argument, e.g. `extent = (Lx, Ly, Lz)`, which specifies the extent of each dimension
in which case 0 ≤ x ≤ Lx, 0 ≤ y ≤ Ly, and -Lz ≤ z ≤ 0.

A grid topology may be specified via a tuple assigning one of `Periodic`, `Bounded`, and `Flat`
to each dimension. By default, a horizontally periodic grid topology `(Periodic, Periodic, Bounded)`
is assumed.

Constants are stored using floating point values of type `FT`. By default this is `Float64`.
Make sure to specify the desired `FT` if not using `Float64`.

Grid properties
===============

- `(Nx, Ny, Nz) :: Int`: Number of physical points in the ``(x, y, z)``-direction.

- `(Hx, Hy, Hz) :: Int`: Number of halo points in the ``(x, y, z)``-direction.

- `(Lx, Ly, Lz) :: FT`: Physical extent of the grid in the ``(x, y, z)``-direction.

- `(Δxᶜᵃᵃ, Δyᵃᶜᵃ, Δzᵃᵃᶜ)`: Grid spacing in the ``(x, y, z)``-direction between cell centers.
                           Defined at cell centers in ``x``, ``y``, and ``z``.

- `(Δxᶠᵃᵃ, Δyᵃᶠᵃ, Δzᵃᵃᶠ)`: Grid spacing in the ``(x, y, z)``-direction between cell faces.
                           Defined at cell faces in ``x``, ``y``, and ``z``.

- `(xᶜᵃᵃ, yᵃᶜᵃ, zᵃᵃᶜ)`: ``(x, y, z)`` coordinates of cell centers.

- `(xᶠᵃᵃ, yᵃᶠᵃ, zᵃᵃᶠ)`: ``(x, y, z)`` coordinates of cell faces.

Examples
========

* A default grid with `Float64` type:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(32, 32, 32), extent=(1, 2, 3))
RectilinearGrid{Float64, Periodic, Periodic, Bounded}
             architecture: CPU()
                   domain: x ∈ [0.0, 1.0], y ∈ [0.0, 2.0], z ∈ [-3.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
        size (Nx, Ny, Nz): (32, 32, 32)
        halo (Hx, Hy, Hz): (1, 1, 1)
             spacing in x: Regular, with spacing 0.03125
             spacing in y: Regular, with spacing 0.0625
             spacing in z: Regular, with spacing 0.09375
```

* A default grid with `Float32` type:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(CPU(), Float32; size=(32, 32, 16), x=(0, 8), y=(-10, 10), z=(-π, π))
RectilinearGrid{Float32, Periodic, Periodic, Bounded}
             architecture: CPU()
                   domain: x ∈ [0.0, 8.0], y ∈ [-10.0, 10.0], z ∈ [-3.1415927, 3.1415927]
                 topology: (Periodic, Periodic, Bounded)
        size (Nx, Ny, Nz): (32, 32, 16)
        halo (Hx, Hy, Hz): (1, 1, 1)
             spacing in x: Regular, with spacing 0.25
             spacing in y: Regular, with spacing 0.625
             spacing in z: Regular, with spacing 0.3926991
```

* A two-dimenisional, horizontally-periodic grid:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=(32, 32), extent=(2π, 4π), topology=(Periodic, Periodic, Flat))
RectilinearGrid{Float64, Periodic, Periodic, Flat}
             architecture: CPU()
                   domain: x ∈ [0.0, 6.283185307179586], y ∈ [0.0, 12.566370614359172], z ∈ [1.0, 1.0]
                 topology: (Periodic, Periodic, Flat)
        size (Nx, Ny, Nz): (32, 32, 1)
        halo (Hx, Hy, Hz): (1, 1, 0)
             spacing in x: Regular, with spacing 0.19634954084936207
             spacing in y: Regular, with spacing 0.39269908169872414
             spacing in z: Flattened
```

* A one-dimensional "column" grid:

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size=256, z=(-128, 0), topology=(Flat, Flat, Bounded))
RectilinearGrid{Float64, Flat, Flat, Bounded}
             architecture: CPU()
                   domain: x ∈ [1.0, 1.0], y ∈ [1.0, 1.0], z ∈ [-128.0, 0.0]
                 topology: (Flat, Flat, Bounded)
        size (Nx, Ny, Nz): (1, 1, 256)
        halo (Hx, Hy, Hz): (0, 0, 1)
             spacing in x: Flattened
             spacing in y: Flattened
             spacing in z: Regular, with spacing 0.5
```

* A horizontally-periodic regular grid with cell interfaces stretched hyperbolically near the top:

```jldoctest
julia> using Oceananigans

julia> σ = 1.1; # stretching factor

julia> Nz = 24; # vertical resolution

julia> Lz = 32; # depth (m)

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = RectilinearGrid(size = (32, 32, Nz),
                       x = (0, 64),
                       y = (0, 64),
                       z = hyperbolically_spaced_faces)
RectilinearGrid{Float64, Periodic, Periodic, Bounded}
             architecture: CPU()
                   domain: x ∈ [0.0, 64.0], y ∈ [0.0, 64.0], z ∈ [-32.0, -0.0]
                 topology: (Periodic, Periodic, Bounded)
        size (Nx, Ny, Nz): (32, 32, 24)
        halo (Hx, Hy, Hz): (1, 1, 1)
             spacing in x: Regular, with spacing 2.0
             spacing in y: Regular, with spacing 2.0
             spacing in z: Stretched, with spacing min=0.6826950100338962, max=1.8309085743885056
```

* A three-dimensional grid with regular spacing in x, cell interfaces that are closely spaced
  close to the boundaries in y (closely mimicing the Chebychev nodes) and cell interfaces
  stretched in z hyperbolically near the top:

```jldoctest
julia> using Oceananigans

julia> Nx, Ny, Nz = 32, 30, 24;

julia> Lx, Ly, Lz = 200, 100, 32; # (m)

julia> chebychev_like_spaced_faces(j) = - Ly/2 * cos(π * (j - 1) / Ny);

julia> σ = 1.1; # stretching factor

julia> hyperbolically_spaced_faces(k) = - Lz * (1 - tanh(σ * (k - 1) / Nz) / tanh(σ));

julia> grid = RectilinearGrid(size = (Nx, Ny, Nz),
                              topology=(Periodic, Bounded, Bounded),
                              x = (0, Lx),
                              y = chebychev_like_spaced_faces,
                              z = hyperbolically_spaced_faces)
RectilinearGrid{Float64, Periodic, Bounded, Bounded}
             architecture: CPU()
                   domain: x ∈ [0.0, 200.0], y ∈ [-50.0, 50.0], z ∈ [-32.0, -0.0]
                 topology: (Periodic, Bounded, Bounded)
        size (Nx, Ny, Nz): (32, 30, 24)
        halo (Hx, Hy, Hz): (1, 1, 1)
             spacing in x: Regular, with spacing 6.25
             spacing in y: Stretched, with spacing min=0.2739052315863333, max=5.226423163382681
             spacing in z: Stretched, with spacing min=0.6826950100338962, max=1.8309085743885056
```
"""
function RectilinearGrid(architecture::AbstractArchitecture = CPU(),
                         FT = Float64;
                         size,
                         x = nothing,
                         y = nothing,
                         z = nothing,
                         halo = nothing,
                         extent = nothing,
                         topology = (Periodic, Periodic, Bounded))

    TX, TY, TZ = validate_topology(topology)
    size = validate_size(TX, TY, TZ, size)
    halo = validate_halo(TX, TY, TZ, halo)

    # Validate the rectilinear domain
    x, y, z = validate_rectilinear_domain(TX, TY, TZ, FT, extent, x, y, z)

    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo

    Lx, xᶠᵃᵃ, xᶜᵃᵃ, Δxᶠᵃᵃ, Δxᶜᵃᵃ = generate_coordinate(FT, topology[1], Nx, Hx, x, architecture)
    Ly, yᵃᶠᵃ, yᵃᶜᵃ, Δyᵃᶠᵃ, Δyᵃᶜᵃ = generate_coordinate(FT, topology[2], Ny, Hy, y, architecture)
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶠ, Δzᵃᵃᶜ = generate_coordinate(FT, topology[3], Nz, Hz, z, architecture)

    FX   = typeof(Δxᶠᵃᵃ)
    FY   = typeof(Δyᵃᶠᵃ)
    FZ   = typeof(Δzᵃᵃᶠ)
    VX   = typeof(xᶠᵃᵃ)
    VY   = typeof(yᵃᶠᵃ)
    VZ   = typeof(zᵃᵃᶠ)
    Arch = typeof(architecture)

    return RectilinearGrid{FT, TX, TY, TZ, FX, FY, FZ, VX, VY, VZ, Arch}(architecture,
        Nx, Ny, Nz, Hx, Hy, Hz, Lx, Ly, Lz, Δxᶠᵃᵃ, Δxᶜᵃᵃ, xᶠᵃᵃ, xᶜᵃᵃ, Δyᵃᶠᵃ, Δyᵃᶜᵃ, yᵃᶠᵃ, yᵃᶜᵃ, Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
end

# architecture = CPU() default, assuming that a DataType positional arg
# is specifying the floating point type.
RectilinearGrid(FT::DataType; kwargs...) = RectilinearGrid(CPU(), FT; kwargs...)

@inline x_domain(grid::RectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TX, grid.Nx, grid.xᶠᵃᵃ)
@inline y_domain(grid::RectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TY, grid.Ny, grid.yᵃᶠᵃ)
@inline z_domain(grid::RectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} = domain(TZ, grid.Nz, grid.zᵃᵃᶠ)

short_show(grid::RectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
    "RectilinearGrid{$FT, $TX, $TY, $TZ}(Nx=$(grid.Nx), Ny=$(grid.Ny), Nz=$(grid.Nz))"

function domain_string(grid::RectilinearGrid)
    x₁, x₂ = domain(topology(grid, 1), grid.Nx, grid.xᶠᵃᵃ)
    y₁, y₂ = domain(topology(grid, 2), grid.Ny, grid.yᵃᶠᵃ)
    z₁, z₂ = domain(topology(grid, 3), grid.Nz, grid.zᵃᵃᶠ)
    return "x ∈ [$x₁, $x₂], y ∈ [$y₁, $y₂], z ∈ [$z₁, $z₂]"
end

function show(io::IO, g::RectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ}
    print(io, "RectilinearGrid{$FT, $TX, $TY, $TZ}\n",
              "             architecture: $(g.architecture)\n",
              "                   domain: $(domain_string(g))\n",
              "                 topology: ", (TX, TY, TZ), '\n',
              "        size (Nx, Ny, Nz): ", (g.Nx, g.Ny, g.Nz), '\n',
              "        halo (Hx, Hy, Hz): ", (g.Hx, g.Hy, g.Hz), '\n',
              "             spacing in x: ", show_coordinate(g.Δxᶜᵃᵃ, TX), '\n',
              "             spacing in y: ", show_coordinate(g.Δyᵃᶜᵃ, TY), '\n',
              "             spacing in z: ", show_coordinate(g.Δzᵃᵃᶜ, TZ))
end

Adapt.adapt_structure(to, grid::RectilinearGrid{FT, TX, TY, TZ}) where {FT, TX, TY, TZ} =
             RectilinearGrid{FT, TX, TY, TZ,
                          typeof(Adapt.adapt(to, grid.Δxᶜᵃᵃ)),
                          typeof(Adapt.adapt(to, grid.Δyᵃᶠᵃ)),
                          typeof(Adapt.adapt(to, grid.Δzᵃᵃᶠ)),
                          typeof(Adapt.adapt(to, grid.xᶠᵃᵃ)),
                          typeof(Adapt.adapt(to, grid.yᵃᶠᵃ)),
                          typeof(Adapt.adapt(to, grid.zᵃᵃᶠ)),
                          Nothing}(
        nothing,
        grid.Nx, grid.Ny, grid.Nz,
        grid.Hx, grid.Hy, grid.Hz,
        grid.Lx, grid.Ly, grid.Lz,
        Adapt.adapt(to, grid.Δxᶠᵃᵃ),
        Adapt.adapt(to, grid.Δxᶜᵃᵃ),
        Adapt.adapt(to, grid.xᶠᵃᵃ),
        Adapt.adapt(to, grid.xᶜᵃᵃ),
        Adapt.adapt(to, grid.Δyᵃᶠᵃ),
        Adapt.adapt(to, grid.Δyᵃᶜᵃ),
        Adapt.adapt(to, grid.yᵃᶠᵃ),
        Adapt.adapt(to, grid.yᵃᶜᵃ),
        Adapt.adapt(to, grid.Δzᵃᵃᶠ),
        Adapt.adapt(to, grid.Δzᵃᵃᶜ),
        Adapt.adapt(to, grid.zᵃᵃᶠ),
        Adapt.adapt(to, grid.zᵃᵃᶜ))

@inline xnode(::Center, i, grid::RectilinearGrid) = @inbounds grid.xᶜᵃᵃ[i]
@inline xnode(::Face  , i, grid::RectilinearGrid) = @inbounds grid.xᶠᵃᵃ[i]

@inline ynode(::Center, j, grid::RectilinearGrid) = @inbounds grid.yᵃᶜᵃ[j]
@inline ynode(::Face  , j, grid::RectilinearGrid) = @inbounds grid.yᵃᶠᵃ[j]

@inline znode(::Center, k, grid::RectilinearGrid) = @inbounds grid.zᵃᵃᶜ[k]
@inline znode(::Face  , k, grid::RectilinearGrid) = @inbounds grid.zᵃᵃᶠ[k]

all_x_nodes(::Type{Center}, grid::RectilinearGrid) = grid.xᶜᵃᵃ
all_x_nodes(::Type{Face}  , grid::RectilinearGrid) = grid.xᶠᵃᵃ
all_y_nodes(::Type{Center}, grid::RectilinearGrid) = grid.yᵃᶜᵃ
all_y_nodes(::Type{Face}  , grid::RectilinearGrid) = grid.yᵃᶠᵃ
all_z_nodes(::Type{Center}, grid::RectilinearGrid) = grid.zᵃᵃᶜ
all_z_nodes(::Type{Face}  , grid::RectilinearGrid) = grid.zᵃᵃᶠ

@inline cpu_face_constructor_x(grid::XRegRectilinearGrid) = x_domain(grid)
@inline cpu_face_constructor_y(grid::YRegRectilinearGrid) = y_domain(grid)
@inline cpu_face_constructor_z(grid::ZRegRectilinearGrid) = z_domain(grid)

function with_halo(new_halo, old_grid::RectilinearGrid)

    size = (old_grid.Nx, old_grid.Ny, old_grid.Nz)
    topo = topology(old_grid)

    x = cpu_face_constructor_x(old_grid)
    y = cpu_face_constructor_y(old_grid)
    z = cpu_face_constructor_z(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    size     = pop_flat_elements(size, topo)
    new_halo = pop_flat_elements(new_halo, topo)

    new_grid = RectilinearGrid(architecture(old_grid), eltype(old_grid);
               size = size,
               x = x, y = y, z = z,
               topology = topo,
               halo = new_halo)

    return new_grid
end

function with_arch(new_arch, old_grid::RectilinearGrid)

    size = (old_grid.Nx, old_grid.Ny, old_grid.Nz)
    topo = topology(old_grid)

    x = cpu_face_constructor_x(old_grid)
    y = cpu_face_constructor_y(old_grid)
    z = cpu_face_constructor_z(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    size = pop_flat_elements(size, topo)
    halo = pop_flat_elements(halo_size(old_grid), topo)

    new_grid = RectilinearGrid(new_arch, eltype(old_grid);
               size = size,
               x = x, y = y, z = z,
               topology = topo,
               halo = halo)

    return new_grid
end

#
# Get minima of grid
#

function min_Δx(grid::RectilinearGrid)
    topo = topology(grid)
    if topo[1] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δxᶜᵃᵃ)
    end
end

function min_Δy(grid::RectilinearGrid)
    topo = topology(grid)
    if topo[2] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δyᵃᶜᵃ)
    end
end

function min_Δz(grid::RectilinearGrid)
    topo = topology(grid)
    if topo[3] == Flat
        return Inf
    else
        return min_number_or_array(grid.Δzᵃᵃᶜ)
    end
end

@inline min_number_or_array(var) = var
@inline min_number_or_array(var::AbstractVector) = minimum(parent(var))


text:
using Flux: @treelike
import Flux
import Base

include("./vocab.jl")
include("./onehot.jl")
include("./gather.jl")
include("./etype.jl")
include("./embed.jl")
include("./position_embed.jl")

ext: jl


text:
export oneplus_populate, ga_populate, nsga2_populate

"""
    nsga2_populate(e::NSGA2Evo)

Complete the population vector of the input `NSGA2Evo` object by mutating
elite/parents individuals that are already contained in this population vector.
"""
function nsga2_populate(e::NSGA2Evo{T}) where T
    parent_index = 1
    while length(e.population) < e.config.n_population
        push!(e.population, mutate(e.population[parent_index]))
        parent_index = mod(parent_index, e.config.n_elite) + 1
    end
end

"""
    function oneplus_populate(e::DualCGPEvolution)

Implements a 1+λ evolution scheme for both population of a Dual CGP Coevolution.
"""
function oneplus_populate(e::DualCGPEvolution)
    enc = Cambrian.max_selection(e.encoder_population)
    ctr = Cambrian.max_selection(e.controller_population)
    e.encoder_population[1] = enc
    e.controller_population[1] = ctr
    for i in 2:length(e.encoder_population)
        e.encoder_population[i] = mutate(enc, "encoder")
    end
    for i in 2:length(e.controller_population)
        e.controller_population[i] = mutate(ctr, "controller")
    end
end

function tournament_selection(e::DualCGPGAEvo, ci::CartesianIndices)
    n_pairs = length(e.fitness_matrix)
    selected = randperm(n_pairs)[1:e.tournament_size]
    fitnesses = e.fitness_matrix[selected]
    best_pair_index = sortperm(fitnesses)[end]
    best_pair_carte = ci[best_pair_index]
    echr = copy(e.encoder_sympop[best_pair_carte[1]].chromosome)
    cchr = copy(e.controller_sympop[best_pair_carte[2]].chromosome)
    return echr, cchr
end

function isin(v::Vector{SymInd}, candidate::SymInd)
    for vi in v
        if candidate.chromosome == vi.chromosome
            return true
        end
    end
    false
end

function find_symind_index(ind::SymInd, v::Vector{SymInd})
    for i in eachindex(v)
        if v[i].chromosome == ind.chromosome
            return i
        end
    end
end

function ga_populate(e::DualCGPGAEvo)
    enew = Vector{SymInd}()
    cnew = Vector{SymInd}()
    n_e = e.encoder_config.n_population
    n_c = e.controller_config.n_population
    # 1. Add elites
    new_elite_indexes = Vector{Tuple{Int64, Int64}}()
    for i in CartesianIndices(e.elites_matrix)
        if e.elites_matrix[i]
            ecandidate = SymIndCopy(e.encoder_sympop[i[1]])
            ccandidate = SymIndCopy(e.controller_sympop[i[2]])
            if !isin(enew, ecandidate)
                push!(enew, ecandidate)
            end
            if !isin(cnew, ccandidate)
                push!(cnew, ccandidate)
            end
            row = find_symind_index(ecandidate, enew)
            col = find_symind_index(ccandidate, cnew)
            push!(new_elite_indexes, (row, col))
        end
    end
    # 2. Run tournaments until population is full
    new_mutant_indexes = Vector{Tuple{Int64, Int64}}()
    ci = CartesianIndices(size(e.fitness_matrix))
    while length(enew) < n_e || length(cnew) < n_c
        echr, cchr = tournament_selection(e, ci)
        if length(enew) < n_e
            type = e.encoder_config.type
        	enco = IPCGPInd(e.encoder_config, echr)
            enco_m = mutate(enco, type)
            push!(enew, SymInd(enco_m.chromosome, 0, type))
        end
        if length(cnew) < n_c
            type = e.controller_config.type
        	cont = CGPInd(e.controller_config, cchr)
            cont_m = mutate(cont, type)
            push!(cnew, SymInd(cont_m.chromosome, 0, type))
        end
        row, col = length(enew), length(cnew)
        push!(new_mutant_indexes, (row, col))
    end
    # 3. Set new population and reset matrices
    e.encoder_sympop = enew
    e.controller_sympop = cnew
    mat_size = size(e.fitness_matrix)
    e.fitness_matrix = -Inf * ones(mat_size...)
    e.elites_matrix = falses(mat_size...)
    e.eval_matrix = falses(mat_size...)
    for index in new_elite_indexes
        e.elites_matrix[index...] = true
        e.eval_matrix[index...] = true
    end
    if e.eval_mutant
        for index in new_mutant_indexes
            e.eval_matrix[index...] = true
        end
    end
end


text:
# importall Base

immutable Dual{T <: Real} <: Number
    re::T
    du::T
end

Dual(x::Real, y::Real) = Dual(promote(x, y)...)

Dual(x::Real) = Dual(x, zero(x))

typealias Dual128 Dual{Float64}
typealias Dual64 Dual{Float32}
typealias DualPair Dual

Base.real(z::Dual) = z.re
epsilon(z::Dual) = z.du

Base.eps(z::Dual) = eps(real(z))
Base.eps{T}(::Type{Dual{T}}) = eps(T)
Base.one(z::Dual) = dual(one(real(z)))
Base.one{T}(::Type{Dual{T}}) = dual(one(T))
Base.inf{T}(::Type{Dual{T}}) = dual(inf(T))
Base.nan{T}(::Type{Dual{T}}) = nan(T)
Base.isnan(z::Dual) = isnan(real(z))

function Base.convert{T<:Real}(::Type{Dual{T}}, x::Real)
    return Dual{T}(convert(T, x), convert(T, 0))
end
Base.convert{T<:Real}(::Type{Dual{T}}, z::Dual{T}) = z
function Base.convert{T <: Real}(::Type{Dual{T}}, z::Dual)
    return Dual{T}(convert(T, real(z)), convert(T, epsilon(z)))
end

function Base.convert{T <: Real}(::Type{T}, z::Dual)
    if epsilon(z) == zero(T)
        return convert(T, real(z))
    else
        throw(InexactError())
    end
end

function Base.promote_rule{T<:Real, S<:Real}(::Type{Dual{T}}, ::Type{Dual{S}})
    return Dual{promote_type(T, S)}
end

# these promotion rules shouldn't be used for scalar operations -- they're slow
Base.promote_rule{T<:Real}(::Type{Dual{T}}, ::Type{T}) = Dual{T}
function Base.promote_rule{T<:Real, S<:Real}(::Type{Dual{T}}, ::Type{S})
  Dual{promote_type(T, S)}
end

dual(x, y) = Dual(x, y)
dual(x) = Dual(x)

Base.@vectorize_1arg Real dual

dual128(x::Float64, y::Float64) = Dual{Float64}(x, y)
dual128(x::Real, y::Real) = dual128(float64(x), float64(y))
dual128(z) = dual128(real(z), epsilon(z))
dual64(x::Float32, y::Float32) = Dual{Float32}(x, y)
dual64(x::Real, y::Real) = dual64(float32(x), float32(y))
dual64(z) = dual64(real(z), epsilon(z))

isdual(x::Dual) = true
isdual(x::Number) = false

Base.real_valued{T<:Real}(z::Dual{T}) = epsilon(z) == 0
Base.integer_valued(z::Dual) = real_valued(z) && integer_valued(real(z))

Base.isfinite(z::Dual) = isfinite(real(z))
Base.reim(z::Dual) = (real(z), epsilon(z))

function dual_show(io::IO, z::Dual, compact::Bool)
    x, y = reim(z)
    if isnan(x) || isfinite(y)
        compact ? showcompact(io,x) : show(io,x)
        if signbit(y)==1 && !isnan(y)
            y = -y
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        compact ? showcompact(io, y) : show(io, y)
        if !(isa(y,Integer) || isa(y,Rational) ||
             isa(y,FloatingPoint) && isfinite(y))
            print(io, "*")
        end
        print(io, "du")
    else
        print(io, "dual(", x, ",", y, ")")
    end
end
Base.show(io::IO, z::Dual) = dual_show(io, z, false)
Base.showcompact(io::IO, z::Dual) = dual_show(io, z, true)

function Base.read{T<:Real}(s::IO, ::Type{Dual{T}})
    x = read(s, T)
    y = read(s, T)
    Dual{T}(x, y)
end
function Base.write(s::IO, z::Dual)
    write(s, real(z))
    write(s, epsilon(z))
end

## Generic functions of dual numbers ##

Base.convert(::Type{Dual}, z::Dual) = z
Base.convert(::Type{Dual}, x::Real) = dual(x)

Base.(:(==))(z::Dual, w::Dual) = real(z) == real(w) && epsilon(z) == epsilon(w)
# ==(z::Dual, x::Real) = real_valued(z) && real(z) == x
# ==(x::Real, z::Dual) = real_valued(z) && real(z) == x

function Base.isequal(z::Dual, w::Dual)
    return isequal(real(z),real(w)) && isequal(epsilon(z), epsilon(w))
end
Base.isequal(z::Dual, x::Real) = real_valued(z) && isequal(real(z), x)
Base.isequal(x::Real, z::Dual) = real_valued(z) && isequal(real(z), x)

Base.isless(z::Dual,w::Dual) = real(z) < real(w)
Base.isless(z::Number,w::Dual) = z < real(w)
Base.isless(z::Dual,w::Number) = real(z) < w

function Base.hash(z::Dual)
    x = hash(real(z));
    return real_valued(z) ? x : bitmix(x, hash(epsilon(z)))
end

# we don't support Dual{Complex}, so conj is a noop
Base.conj(z::Dual) = z
Base.abs(z::Dual)  = (real(z) >= 0) ? z : -z
Base.abs2(z::Dual) = z*z

# algebraic definitions
conjdual(z::Dual) = Dual(real(z),-epsilon(z))
absdual(z::Dual) = abs(real(z))
abs2dual(z::Dual) = abs2(real(z))

Base.(:(+))(z::Dual, w::Dual) = dual(real(z)+real(w), epsilon(z)+epsilon(w))
Base.(:(+))(z::Number, w::Dual) = dual(z+real(w), epsilon(w))
Base.(:(+))(z::Dual, w::Number) = dual(real(z)+w, epsilon(z))

Base.(:(-))(z::Dual) = dual(-real(z), -epsilon(z))
Base.(:(-))(z::Dual, w::Dual) = dual(real(z)-real(w), epsilon(z)-epsilon(w))
Base.(:(-))(z::Number, w::Dual) = dual(z-real(w), -epsilon(w))
Base.(:(-))(z::Dual, w::Number) = dual(real(z)-w, epsilon(z))

# avoid ambiguous definition with Bool * Number
function Base.(:(*))(x::Bool, z::Dual)
    return ifelse(x, z, ifelse(signbit(real(z)) == 0, zero(z), -zero(z)))
end
Base.(:(*))(x::Dual, z::Bool) = z * x

function Base.(:(*))(z::Dual, w::Dual)
    return dual(real(z) * real(w), epsilon(z) * real(w) + real(z) * epsilon(w))
end
Base.(:(*))(x::Real, z::Dual) = dual(x * real(z), x * epsilon(z))
Base.(:(*))(z::Dual, x::Real) = dual(x * real(z), x * epsilon(z))

Base.(:(/))(z::Real, w::Dual) = dual(z / real(w), -z * epsilon(w) / real(w)^2)
Base.(:(/))(z::Dual, x::Real) = dual(real(z) / x, epsilon(z) / x)
function Base.(:(/))(z::Dual, w::Dual)
    return dual(
        real(z) / real(w),
        (epsilon(z) * real(w) - real(z) * epsilon(w)) / (real(w) * real(w))
    )
end

function Base.(:(^))(z::Dual, w::Dual)
    re = real(z)^real(w)
    du = epsilon(z) * real(w) * (real(z)^(real(w) - 1))
    du += epsilon(w) * (real(z)^real(w)) * log(real(z))
    return dual(re, du)
end

# these two definitions are needed to fix ambiguity warnings
function Base.(:(^))(z::Dual, n::Integer)
    return dual(real(z)^n, epsilon(z) * n * real(z)^(n - 1))
end
function Base.(:(^))(z::Dual, n::Rational)
    return dual(real(z)^n, epsilon(z) * n * real(z)^(n - 1))
end
function Base.(:(^))(z::Dual, n::Real)
    return dual(real(z)^n, epsilon(z) * n * real(z)^(n - 1))
end

for (funsym, exp) in derivative_rules
    @eval function (Base.$funsym)(z::Dual)
        xp = epsilon(z)
        x = real(z)
        Dual($(funsym)(x),$exp)
    end
end


text:
function get_afac(n::Float64)
    if n > 1.0
        A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
        5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
    elseif n < 1.0
        A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
    else
        A = 1.0
    end
    return A
end

### Semi-infinite models ###

function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{refl_DA_semiinf_TD})
    Rt = refl_DA_semiinf_TD(data.t, β, data.ndet, data.nmed)
end

# fit DTOF to fluence with EBPC conditions
function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{fluence_DA_semiinf_TD})
    Rt = fluence_DA_semiinf_TD(data.t, β, data.ndet, data.nmed, 0.0)
    return Rt./(2*get_afac(data.nmed/data.ndet))
end

### Slab models ###

function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{refl_DA_slab_TD})
    Rt = refl_DA_slab_TD(data.t, β, data.ndet, data.nmed, data.s)
end

### Parallelepiped models ###

function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{refl_DA_paralpip_TD})
    Rt = refl_DA_paralpip_TD(data.t, β, data.ndet, data.nmed, data.rd, data.rs, data.L)
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
# PGLib exports everything except internal symbols, which are defined as
# those whose name starts with an underscore. If you don't want all of these
# symbols in your environment, then use `import PGLib` instead of
# `using PGLib`.

const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
         Base.isidentifier(sym_string[2:end])))
       continue
    end
    #println("$(sym)")
    @eval export $sym
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Pre-defined highlighters for the text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export hl_cell, hl_col, hl_row, hl_lt, hl_leq, hl_gt, hl_geq, hl_value

"""
    hl_cell(i::Number, j::Number, crayon::Crayon)

Highlight the cell `(i,j)` with the crayon `crayon`.

    hl_cell(cells::AbstractVector{NTuple(2,Int)}, crayon::Crayon)

Highlights all the cells in `cells` with the crayon `crayon`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_cell(i::Number, j::Number, crayon::Crayon) = Highlighter(
    f = (data,x,y)->begin
        return (x == i) && (y == j)
    end,
    crayon = crayon
)

hl_cell(cells::AbstractVector{NTuple{2,Int}}, crayon::Crayon) = Highlighter(
    f = (data,x,y)->begin
        return (x,y) ∈ cells
    end,
    crayon = crayon
)

"""
    hl_col(i::Number, crayon::Crayon)

Highlight the entire column `i` with the crayon `crayon`.

    hl_col(cols::AbstractVector{Int}, crayon::Crayon)

Highlights all the columns in `cols` with the crayon `crayon`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_col(j::Number, crayon::Crayon) = Highlighter(
    f = (data,x,y)->begin
        return y == j
    end,
    crayon = crayon
)

hl_col(cols::AbstractVector{Int}, crayon::Crayon) = Highlighter(
    f = (data,x,y)->begin
        return y ∈ cols
    end,
    crayon = crayon
)

"""
    hl_row(i::Number, crayon::Crayon)

Highlight the entire row `i` with the crayon `crayon`.

    hl_row(rows::AbstractVector{Int}, crayon::Crayon)

Highlights all the rows in `rows` with the crayon `crayon`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_row(i::Number, crayon::Crayon) = Highlighter(
    f = (data,x,y)->begin
        return x == i
    end,
    crayon = crayon
)

hl_row(rows::AbstractVector{Int}, crayon::Crayon) = Highlighter(
    f = (data,x,y)->begin
        return x ∈ rows
    end,
    crayon = crayon
)

"""
    hl_lt(n::Number)

Highlight all elements that < `n`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_lt(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(<,data[i,j],n) && data[i,j] < n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :red)
)

"""
    hl_leq(n::Number)

Highlight all elements that ≤ `n`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_leq(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(≤,data[i,j],n) && data[i,j] ≤ n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :red)
)

"""
    hl_gt(n::Number)

Highlight all elements that > `n`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_gt(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(>,data[i,j],n) && data[i,j] > n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :blue)
)

"""
    hl_geq(n::Number)

Highlight all elements that ≥ `n`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_geq(n::Number) = Highlighter(
    f = (data,i,j)->begin
        if applicable(≥,data[i,j],n) && data[i,j] ≥ n
            return true
        else
            return false
        end
    end,
    crayon = Crayon(bold       = true,
                    foreground = :blue)
)

"""
    hl_value(v::Any)

Highlight all the values that matches `data[i,j] == v`.

# Remarks

Those functions return a `Highlighter` to be used with the text backend.

"""
hl_value(v) = Highlighter(
    f = (data,i,j)->data[i,j] == v,
    crayon = Crayon(bold       = true,
                    foreground = :yellow)
)


text:
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    return sense(model) == 1 ? MOI.MIN_SENSE : MOI.MAX_SENSE
end
function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    sense!(model, Val(sense == MOI.MAX_SENSE ? :max : :min))
end

"""
    ScalarFunction{F <: Function, V <: Union{Nothing, VOV}} <: MOI.AbstractScalarFunction

A container to express any function with real value in JuMP syntax. Used with the `@objective` macro.

# Arguments:
- `f::F`: function to be applied to `X`
- `X::V`: a subset of the variables of the model.

Given a `model`, and some (collection of) variables `X` to optimize. an objective function `f` can be added as follows. Note that only `Min` for minimization us currently defined. `Max` will come soon.

```julia
# Applies to all variables in order of insertion.
# Recommended only when the function argument order does not matter.
@objective(model, ScalarFunction(f))

# Generic use
@objective(model, ScalarFunction(f, X))
```
"""
struct ScalarFunction{F <: Function, V <: Union{Nothing, VOV}} <: MOI.AbstractScalarFunction
    f::F
    X::V

    ScalarFunction(f, X::Union{Nothing, VOV} = nothing) = new{typeof(f), typeof(X)}(f, X)
end

# external constructors

function ScalarFunction(f, X::A) where {A <: AbstractArray{VariableRef}}
    return ScalarFunction(f, VOV(vec(map(index, X))))
end
ScalarFunction(f, x::VariableRef) = ScalarFunction(f, [x])
ScalarFunction(f, x::VI) = ScalarFunction(f, VOV([x]))

# copy
Base.copy(func::ScalarFunction) = ScalarFunction(func.f, func.X)

# supports
MOI.supports(::Optimizer, ::OF{ScalarFunction{F, V}}) where {F <: Function, V <: Union{Nothing, SVF,VOV}} = true

# set
function MOI.set(optimizer::Optimizer, ::OF, func::ScalarFunction{F, Nothing}
) where {F <: Function} # for direct use of MOI
    return objective!(optimizer, func.f)
end

function MOI.set(optimizer::Optimizer, ::OF, func::ScalarFunction{F, VOV}
) where {F <: Function} # VOV, mainly for JuMP
        objective_func = _ -> func.f(map(y -> get_value(optimizer,y.value), func.X.variables))
        return objective!(optimizer, objective_func)
 end

#  @autodoc
 function MOIU.map_indices(index_map::Function, sf::ScalarFunction{F,VOV}
) where {F <: Function}
    return ScalarFunction(sf.f, MOIU.map_indices(index_map, sf.X))
 end

#  @autodoc
 function MOIU.map_indices(::Function, sf::ScalarFunction{F,Nothing}) where {F <: Function}
        return ScalarFunction(sf.f, nothing)
end


text:
# This file includes ODE networks.

"""
    $TYPEDSIGNATURES

Returns the differential equaiton problem of `net`. `tspan` is the time interval of the problem.
"""
function netproblem(net::ODENetwork, tspan = (0.0, 10.0), args...; kwargs...)
    f = function rhs(dx, x, (net, args...), t)
        @views for (i, idx) in Iterators.partition(1:length(x), size(net.coupling, 1)) |> enumerate
            net.nodes[i](dx[idx], x[idx], nothing, t)
        end
        W = laplacian_matrix(net.graph)
        C = W ./ diag(W)
        dx .+= -net.ϵ * kron(C, net.coupling) * x
    end
    x0 = vcat([node.x for node in net.nodes]...)
    ODEProblem(f, x0, tspan, (net, args...); kwargs...)
end

"""
    $TYPEDSIGNATURES

Returns the solution of the `net` problem.
"""
function solvenet!(net::AbstractNetwork, tspan = (0.0, 10.0), probargs = (), solverargs = (); alg = getsolver(net), probkwargs = (;), solverkwargs = (;))
    if typeof(net) <: ODENetwork
        solvenet!(netproblem(net, tspan, probargs...; probkwargs...), alg, solverargs...; solverkwargs...)
    else
        d = size(net.coupling, 1)
        n = length(net.nodes)
        l = size(net.noisemat, 2)
        noise_rate_prototype = zeros(n * d, l * d)
        solvenet!(
            netproblem(net, tspan, probargs...; noise_rate_prototype = noise_rate_prototype, probkwargs...),
            alg, solverargs...; solverkwargs...
        )
    end
end

"""
    $TYPEDSIGNATURES

Returns the solution of network problem `prob` with `alg`.
"""
function solvenet!(prob::DEProblem, alg = getsolver(prob), args...; kwargs...)
    integ = initnet(prob, alg, args...; kwargs...)
    _, tf = prob.tspan
    callback = get(prob.kwargs, :callback, nothing)
    callback === nothing || callback.affect!(integ)
    while integ.t < tf
        step!(integ)
    end
    updatenodes!(integ.sol)
end

"""
    $TYPEDSIGNATURES

Returns and differential equation integrator corresponding to `prob` and `alg`.
"""
initnet(prob::DEProblem, alg = getsolver(prob), args...; kwargs...) =
    init(prob, alg, args...; maxiters = typemax(UInt), kwargs...)

# Updates the nodes of the network with solution `sol`
function updatenodes!(sol)
    t, x = last(sol.t), last(sol.u)
    net = first(sol.prob.p)
    d = size(net.coupling, 1)
    @views for (i, idx) in Iterators.partition(1:length(x), d) |> enumerate
        node = net.nodes[i]
        node.x .= x[idx]
        node.t = t
    end
    sol
end


text:
"""
    FastLip(maxIter::Int64, ϵ0::Float64, accuracy::Float64)

FastLip adds Lipschitz estimation on top of FastLin.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hypercube
3. Output: halfspace

# Return
`AdversarialResult`

# Method
Lipschitz estimation + FastLin. All arguments are for FastLin.
- `max_iter` is the maximum iteration in search, default `10`;
- `ϵ0` is the initial search radius, default `100.0`;
- `accuracy` is the stopping criteria, default `0.1`;

# Property
Sound but not complete.

# Reference
[T.-W. Weng, H. Zhang, H. Chen, Z. Song, C.-J. Hsieh, D. Boning, I. S. Dhillon, and L. Daniel,
"Towards Fast Computation of Certified Robustness for ReLU Networks,"
*ArXiv Preprint ArXiv:1804.09699*, 2018.](https://arxiv.org/abs/1804.09699)
"""
@with_kw struct FastLip
    maxIter::Int64    = 10
    ϵ0::Float64       = 100.0
    accuracy::Float64 = 0.1
end

convert(::Type{FastLin}, S::FastLip) = FastLin(S.maxIter, S.ϵ0, S.accuracy)

function solve(solver::FastLip, problem::Problem)
    c, d = tosimplehrep(HPolytope(problem.output))
    y = compute_output(problem.network, problem.input.center)
    o = (c * y - d)[1]
    if o > 0
        return AdversarialResult(:UNSAT, -o)
    end
    result = solve(FastLin(solver), problem)
    result.status == :UNSAT && return result
    ϵ_fastLin = result.max_disturbance
    LG, UG = get_gradient(problem.network, problem.input)

    # C = problem.network.layers[1].weights
    # L = zeros(size(C))
    # U = zeros(size(C))
    # for l in 2:length(problem.network.layers)
    #   C, L, U = bound_layer_grad(C, L, U, problem.network.layers[l].weights, act_pattern[l])
    # end
    # v = max.(abs.(C+L), abs.(C+U))

    a, b = interval_map(c, LG, UG)
    v = max.(abs.(a), abs.(b))
    ϵ = min(-o/sum(v), ϵ_fastLin)

    if ϵ > minimum(problem.input.radius)
        return AdversarialResult(:SAT, ϵ)
    else
        return AdversarialResult(:UNSAT, ϵ)
    end
end

function bound_layer_grad(C::Matrix, L::Matrix, U::Matrix, W::Matrix, D::Vector{Float64})
    n_input = size(C)
    rows, cols = size(W)
    new_C = zeros(rows, n_input)
    new_L = zeros(rows, n_input)
    new_U = zeros(rows, n_input)
    for k in 1:n_input         # NOTE n_input is 2-D
        for j in 1:rows, i in 1:cols

            u = U[i, k]
            l = L[i, k]
            c = C[i, k]
            w = W[j, i]

            if D[i] == 1

                new_C[j,k] += w*c
                new_U[j,k] += (w > 0) ? u : l
                new_L[j,k] += (w > 0) ? l : u

            elseif D[i] == 0 && w*(c+u)>0

                new_U[j,k] += w*(c+u)
                new_L[j,k] += w*(c+l)

            end
        end
    end
    return (new_C, new_L, new_U)
end


text:
module LHEF

import Base.dot

using LightXML
using Compat

immutable EventHeader
    nup::Uint8          # Number of particles
    ldprup::Uint8       # Process type?
    xwgtup::Float64     # Event Wight
    scalup::Float64     # Scale
    αem::Float64        # AQEDUP
    αs::Float64         # AQCDUP
end

immutable FourVector
    data::NTuple{4,Float64}
end
getindex(x::FourVector,i) = x.data[i+1]

# Mostly negative convention
dot(x::FourVector,y::FourVector) = (x[0]*y[0]-x[1]*y[1]-x[2]*y[2]-x[3]*y[3])

immutable Particle
    particle::Int8
    status::Int8
    mothup::NTuple{2,Uint8}
    color::NTuple{2,Uint16}
    pμ::FourVector
    m::Float64
    vtimup::Float64
    spinup::Float64
end

immutable Event
    header::EventHeader
    data::Vector{Particle}
end

function parse_lhe(filename; format = nothing)
    if format === nothing
        # Format not declared, inferring from extension
        fparts = split(basename(filename),".")
        if fparts[end] == "lhe"
            format = :lhe
        elseif fparts[end] == "gz" && fparts[end-1] == "lhe"
            format = :lhegz
        end
    end

    @assert format == :lhe
    lhefile = parse_file(filename)
    lhenode = root(lhefile)

    (name(lhenode) == "LesHouchesEvents") || error("Invalid root node")
    (attributes_dict(lhenode)["version"] == "3.0") || error("Unsupported Version")

    events = get_elements_by_tagname(lhenode,"event")

    [begin
        data = content(event)
        lines = @compat split(data,'\n',keep=false)
        headerdata = @compat split(lines[1],' ',keep=false)
        header = EventHeader(parseint(headerdata[1]), parseint(headerdata[2]),
            parsefloat(Float64, headerdata[3]), parsefloat(Float64, headerdata[4]),
            parsefloat(Float64, headerdata[5]), parsefloat(Float64, headerdata[6]))
        data = [begin
            fields = @compat split(line,' ',keep=false)
            p = Particle(parseint(Int8,fields[1]),parseint(Int8,fields[2]),
                (parseint(Uint8,fields[3]),parseint(Uint8,fields[4])),
                (parseint(Uint16,fields[5]),parseint(Uint16,fields[6])),
                FourVector((parsefloat(Float64,fields[10]),parsefloat(Float64,fields[7]),
                    parsefloat(Float64,fields[8]),parsefloat(Float64,fields[9]))),
                parsefloat(Float64, fields[11]), parsefloat(Float64, fields[12]),
                parsefloat(Float64, fields[13]))
            p
        end for line in lines[2:end]]
        Event(header,data)
    end for event in events]
end

# package code goes here

end # module


text:
immutable FunctionSignature
    name::String
    argtypes::Vector{DataType}
    returntype::DataType
end

immutable Backend{Identifier}
    typemap::Dict{DataType, String} #
    funcnames::Dict{Function, Dict{Pair, String}}
end

function Backend(sym::Symbol)
    Backend{sym}(
        Dict{DataType, String}(),
        Dict{Function, Dict{Pair, String}}()
    )
end

"""
Most backends don't support type parameters, so we need to embedd them in the name
"""
function to_type_name(T)
    join((a.name.name, a.parameters...), "_")
end

function to_type_name{T}(backend, ::Type{T})
    if haskey(backend.typemap, T)
        return backend.typemap[T]
    else
        _declare_composite_type(backend, T)
    end
end

function FunctionSignature(backend, f, types::Tuple{Vararg{DataType}})
    signature = get!(backend.funcnames, (fun2name(f), types)) do

    end
    get!(backend.declared_functions, signature) do

        return declared_functions.typemap[T]
    else
        _declare_function(backend, f, types::Tuple)
    end
end
function to_type_assert(backend, slot, typ)
    error("not implemented for $backend")
end

function _declare_function(backend, f, types)
    error("not implemented for $backend")
end
function declare_function(backend, name, args, returntype, body)
    error("not implemented for $backend")
end
function _declare_composite_type{T}(backend, ::Type{T})
    name = declare_composite_type(backend, T)
    backend.declared_types[T] = name
    name
end
function declare_composite_type(backend, ::Type{T})
    error("not implemented for $backend")
end


text:
using covid19abm
const cv = covid19abm
using Test

const infectiousstates = (cv.LAT, cv.MILD, cv.MISO, cv.INF, cv.IISO, cv.HOS, cv.ICU)

@testset "parameters" begin
    ip = ModelParameters()
    mod_ip = cv.p ## get the module parameters
    ## set random parameters
    ip.β = 99
    ip.prov = :noknownlocation
    ip.τmild = 1 ## days before they self-isolate for mild cases
    ip.fmild = 0.05  ## percent of people practice self-isolation
    ip.fsevere = 0.80 # fixed at 0.80, always within 1 day.
    ip.eldq = 0.0 ## percentage quarantine of 60+ individuals, removes from chain of transmission
    ip.calibration = false
    ip.modeltime = 500
    reset_params(ip)
    for x in propertynames(cv.p)
        @test getfield(ip, x) == getfield(mod_ip, x)
    end
end

@testset "init" begin
    cv.reset_params_default()
    initialize()

    ## age groups
    inmodel = (1, 2, 3, 4, 5) ## possible age groups in model
    ags = []
    for x in cv.humans
        push!(ags, x.ag)
        @test x.ag in inmodel
        @test x.exp == 999
        @test x.health == cv.SUS
        @test x.swap == cv.UNDEF
        @test x.iso == false
        @test x.isovia == :null
        @test x.tis == 0
        @test x.exp == 999
        @test x.sickfrom == cv.UNDEF
    end
    @test length(unique(ags)) == length(inmodel) # check if all age groups are sampled

    ## insert_infected check if infected person is added in the correct age group
    for ag in inmodel
        initialize() # reset population
        insert_infected(cv.INF, 1, ag) # 1 infected in age group ag
        @test length(findall(x -> x.health == cv.INF && x.ag == ag, cv.humans)) == 1
    end

    ## check if the initial infected person is NOT IISO
    cv.p.fsevere = 0.0
    initialize() # reset population
    insert_infected(cv.INF, 1, 1) # 1 infected in age group ag
    @test length(findall(x -> x.health == cv.INF && x.swap == cv.REC, cv.humans)) == 1

    ## check if the initial infected person is REC (since simpleinf function only puts them to REC)
    cv.p.fsevere = 1.0
    initialize() # reset population
    insert_infected(cv.INF, 1, 1) # 1 infected in age group ag
    @test length(findall(x -> x.health == cv.INF && x.swap == cv.REC, cv.humans)) == 1

    ## check if durations are properly set
    initialize()
    randhumans = rand(1:10000, 100) # select 100 random humans, no need to test all 10000
    for i in randhumans
        @test minimum(humans[i].dur) != 0
    end
end

@testset "transitions" begin
    cv.reset_params_default()
    initialize()

    ## check if time in state is up by one
    cv.time_update()
    for x in cv.humans
        @test x.tis == 1
        @test x.exp == 999
        @test x.health == cv.SUS
        @test x.swap == cv.UNDEF ## shouldn't set a swap
    end

    for x in cv.humans
        x.exp = 0 ## this will trigger a swap
    end
    @test_throws ErrorException("swap expired, but no swap set.") cv.time_update()

    #check if it goes through all the move compartments and see if health/swap changes

    for h in 2:9  ## latent to ded, ignore susceptible
        initialize()
        rh = cv.HEALTH(h)
        for x in humans
            x.swap = rh
            x.exp = 0 ## to force the swap
        end
        time_update() ## since tis > exp
        for x in humans[1:5]
            @test x.health == rh
            if rh ∈ infectiousstates  ## for all states, except rec/ded there should be a swap
                @test x.swap != cv.UNDEF
            end
        end
    end

    # CHECKING LATENT
    initialize()
    x = humans[1]
    x.ag = 1 ## move to the first age group manually.
    myp = cv.ModelParameters()

    # check basic variables
    cv.reset_params(myp)
    move_to_latent(x)
    @test x.swap ∈ (cv.ASYMP, cv.PRE)
    @test x.iso == false
    @test x.doi == 0
    @test x.tis == 0
    @test x.exp == x.dur[1]

    # check split between asymp/pre
    myp.fasymp = 0.0
    cv.reset_params(myp)
    move_to_latent(x)
    @test x.swap == cv.PRE

    myp.fasymp = 1.0
    cv.reset_params(myp)
    move_to_latent(x)
    @test x.swap == cv.ASYMP

    # CHECKING PRE
    initialize()
    x = humans[1]
    cv.reset_params(myp)
    move_to_pre(x)
    @test x.health == cv.PRE
    @test x.swap ∈ (cv.MILD, cv.INF)

    ## check if individual moves through mild, miso through fmild, tmild parameters
    initialize()
    x = humans[1]
    x.ag = 1 ## move to the first age group manually.
    x.iso = false ## turn this off so we can test the effect of fmild, tmild
    myp.fmild = 0.0
    cv.reset_params(myp)
    move_to_mild(x)
    @test x.swap == cv.REC

    myp.fmild = 1.0
    myp.τmild = 1
    cv.reset_params(myp)
    move_to_mild(x)
    @test x.swap == cv.MISO
    @test x.exp == 1

    time_update()
    @test x.health == cv.MISO
    @test x.swap == cv.REC

    ## todo: check if x.iso and x.isovia property are set correctly.
    # 1) check if they are isolated through quarantine
    myp = cv.ModelParameters()
    myp.eldq = 1.0  ## turn on eldq
    myp.eldqag = 4
    cv.reset_params(myp)
    cv.initialize()
    for x in cv.humans
        if x.ag  == myp.eldqag
            @test x.iso == true
            @test x.isovia == :qu
        else
            @test x.iso == false
            @test x.isovia == :null
        end
    end
    # 2) check if they are isolated through presymptomatic capture
    # no need to test for mild/inf movement since it's a simple assignment in code
    # more important to check through `main` when all the function dynamics are happening.
    myp = cv.ModelParameters()
    myp.eldq = 0.0  ## turn off eldq
    myp.fsevere = 0.0 ## turn on fsevere
    myp.fmild = 0.0
    myp.fpreiso = 1.0
    cv.reset_params(myp)
    cv.initialize()
    for x in cv.humans
        cv.move_to_pre(x)
        @test x.iso == true && x.isovia == :pi
    end
end

@testset "tranmission" begin
    cv.reset_params_default()
    cv.initialize()
    grps = cv.get_ag_dist()

    # since beta = 0 default, and everyone sus
    totalinf = cv.dyntrans(1, grps)
    @test totalinf == 0

    # check with only a single infected person
    insert_infected(cv.INF, 1, 1) # 1 infected in age group ag
    totalinf = cv.dyntrans(1, grps)
    @test totalinf == 0 # still zero cuz beta = 0

    # now change beta
    cv.reset_params(ModelParameters(β = 1.0))
    totalinf = cv.dyntrans(1, grps)
    @test totalinf > 0 ## actually may still be zero because of stochasticity but very unliekly

    ## somehow check the transmission reduction and number of contacts
end

@testset "calibration" begin
    # to do, run the model and test total number of infections
    myp = ModelParameters()
    myp.β = 1.0
    myp.prov = :ontario
    myp.calibration = true
    myp.fmild = 0.0
    myp.fsevere = 0.0
    myp.fpreiso = 0.0
    myp.fasymp = 0.5
    myp.initialinf = 1
    cv.reset_params(myp)
    cv.initialize()
    grps = cv.get_ag_dist()
    cv.insert_infected(cv.PRE, 1, 4)
    # find the single insert_presymptomatic person
    h = findall(x -> x.health == cv.PRE, cv.humans)
    x = humans[h[1]]
    @test length(h) == 1
    @test x.ag == 4
    @test x.swap ∈ (cv.ASYMP, cv.MILD, cv.INF) ## always true for calibration
    for i = 1:20 ## run for 20 days
        cv.dyntrans(i, grps)
        cv.time_update()
    end
    @test x.health == cv.REC  ## make sure the initial guy recovered
    @test x.exp == 999
    ## everyone should really be in latent (or susceptible) except the recovered guy
    all = findall(x -> x.health ∈ (cv.PRE, cv.ASYMP, cv.MILD, cv.MISO, cv.INF, cv.IISO, cv.HOS, cv.ICU, cv.DED), cv.humans)
    @test length(all) == 0
    ## to do, make sure everyone stays latent

end

@testset "contact trace" begin
    # todo: check interaction between doi and tracestart, default values
    # myp = ModelParameters()
    # myp.β = 1.0
    # cv.reset_params(myp)
    # cv.initialize()
    # hdx = rand(1:10000, 20) ## sample 20 humans instead of all 10000
    # for i in hdx
    #     x = cv.humans[i]
    #     cv.ct_dynamics(x)
    #     @test x.tracing == false ## since fctcapture = 0.0
    #     @test x.tracinguntil == -1
    #     @test x.tracedby == 0
    #     @test x.tracedxp == 0
    # end
    # ## test the contact_tracing() function
    # cv.initialize()
    # grps = cv.get_ag_dist()
    # myp.fctcapture = 1.0
    # myp.fasymp = 0.0  # to force no asymp
    # cv.reset_params(myp)
    # cv.initialize()
    # tracer = cv.humans[1]
    # cv.move_to_pre(tracer) # newly presymptomatic
    # cv.ct_dynamics(tracer) # will turn on tracing
    # @test tracer.tracing == true
    # @test tracer.tracinguntil == 3

    # cv.dyntrans(1, grps) ## go through a single tranmission cycle
    # alltraced = findall(x -> x.tracedby > 0, cv.humans) ## since we havn't turned fctcapture > 0
    # @test length(alltraced) > 0
    # for i in alltraced
    #     y = cv.humans[i]
    #     @test y.tracedby == 1 ## since we used the first human as the tracing contact
    #     @test y.tracedxp == 14
    #     @test y.iso == false ## the first human is not in INF stage yet.. still in presymp
    #     @test y.isovia == :null
    # end
    # ## time update calls the ctdynamics function.. have to use this because otherwise tracer is stuck on PRE day 1.
    # cv.time_update()
    # #cv.ct_dynamics(tracer) # will turn on tracing
    # @test tracer.tracing == true
    # @test tracer.tracinguntil == 2
    # for i in alltraced
    #     y = cv.humans[i]
    #     @test y.tracedby == 1 ## since we used the first human as the tracing contact
    #     @test y.tracedxp == 14
    #     @test y.iso == false ## the first human is not in INF stage yet.. still in presymp
    #     @test y.isovia == :null
    # end
    # cv.time_update()
    # @test tracer.tracing == true
    # @test tracer.tracinguntil == 1
    # for i in alltraced
    #     y = cv.humans[i]
    #     @test y.tracedby == 1 ## since we used the first human as the tracing contact
    #     @test y.tracedxp == 14
    #     @test y.iso == false ## the first human is not in INF stage yet.. still in presymp
    #     @test y.isovia == :null
    # end
    # cv.time_update()
    # @test tracer.tracing == false
    # @test tracer.tracinguntil == 0
    # # for i in alltraced
    # #     y = cv.humans[i]
    # #     @test y.tracedby == 1
    # #     @test y.tracedxp == 13 ## one less day from 14
    # #     @test y.iso == true ## the first human is not in INF stage yet.. still in presymp
    # #     @test y.isovia == :ct
    # # end
    myp = cv.ModelParameters()
    myp.fctcapture = 1.0 # force contact tracing
    myp.fasymp = 0.0 # force to pre
    cv.reset_params(myp)
    cv.initialize()
    x = cv.humans[1]
    @test x.tracestart == -1
    @test x.traceend == -1

    cv.move_to_latent(x)
    cv.ct_dynamics(x) # since ct_dynamics will run when move_to_latent will run in time update func
    @test x.swap == cv.PRE
    @test x.doi == 0
    @test x.tis == 0
    @test x.exp == x.dur[1]
    @test x.tracestart > 0
    @test x.traceend > 0

    # go through entire inf period
    totalinfperiod = x.dur[1] + x.dur[3] + x.dur[4]
    totaltracedays = 0
    for i = 1:totalinfperiod
        # dyntrans
        cv.time_update()
        @test x.doi == i # should increase with the counter.
        if x.doi >= x.tracestart && x.doi < x.traceend
            @test x.tracing == true
            totaltracedays += 1
            println("trace true in health: $(x.HEALTH), tis: $(x.tis)")
        end
    end
    @test totaltracedays == cv.p.cdaysback
end

@testset "main run" begin
    ## run model with high beta
    myp = cv.ModelParameters()
    myp.β = 0.0525
    myp.prov = :newyork
    ## run empty scenario
    ## this wont return calibrated scenario since fasymp = 0
    myp.τmild = 0
    myp.fmild = 0.0
    myp.fsevere = 0.0
    myp.eldq = 0.0
    myp.fasymp = 0.5
    myp.fpre = 1.0
    myp.fpreiso = 0.0
    myp.tpreiso = 0
    myp.fctcapture = 0.0
    cv.runsim(1, myp) # warm up the functions
    println("time with contact tracing off:")
    @time results = cv.main(myp)

    myp.fctcapture = 1.0
    cv.runsim(1, myp)
    println("time with contact tracing on:")
    @time results = cv.main(myp)

#     function string_or_char(s::SubString)
#         length(s) == 1 && return first(s)
#         return String(s)
#     end
#     arr = Union{Char,String}[]
# for i in split_string
#     push!(arr, string_or_char(i))
# end

    # function addday(n)
    #     strtorep = "---" * "\u00B7"
    #     str = ""
    #     for i = 1:n
    #         str =str * strtorep
    #     end
    #     str = str[1:end-1]
    #     str = " " * str * " "
    #     return str
    # end

    # function printascii_epi(x)
    #     fs = ""
    #     fs = fs * "L" * addday(x.dur[1])

    #     fs = fs * "P" * addday(x.dur[3])
    #     fs = fs * "S" * addday(x.dur[4])
    #     return fs
    # end

end


text:
# ---
# title: 919. Complete Binary Tree Inserter
# id: problem919
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/complete-binary-tree-inserter/description/>
# hidden: true
# ---
#
# A _complete_ binary tree is a binary tree in which every level, except
# possibly the last, is completely filled, and all nodes are as far left as
# possible.
#
# Write a data structure `CBTInserter` that is initialized with a complete
# binary tree and supports the following operations:
#
#   * `CBTInserter(TreeNode root)` initializes the data structure on a given tree with head node `root`;
#   * `CBTInserter.insert(int v)` will insert a `TreeNode` into the tree with value `node.val = v` so that the tree remains complete, **and returns the value of the parent of the inserted`TreeNode`**;
#   * `CBTInserter.get_root()` will return the head node of the tree.
#
#
#
# **Example 1:**
#
#
#
#     Input: inputs = ["CBTInserter","insert","get_root"], inputs = [[[1]],[2],[]]
#     Output: [null,1,[1,2]]
#
#
# **Example 2:**
#
#
#
#     Input: inputs = ["CBTInserter","insert","insert","get_root"], inputs = [[[1,2,3,4,5,6]],[7],[8],[]]
#     Output: [null,3,4,[1,2,3,4,5,6,7,8]]
#
#
#
# **Note:**
#
#   1. The initial given tree is complete and contains between `1` and `1000` nodes.
#   2. `CBTInserter.insert` is called at most `10000` times per test case.
#   3. Every value of a given or inserted node is between `0` and `5000`.
#
#
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    create_application_command(c::Client; kwargs...) -> ApplicationCommand

Creates a global [`ApplicationCommand`](@ref).
"""
function create_application_command(c::Client; kwargs...)
    return Response{ApplicationCommand}(c, :POST, "/applications/$appid/commands"; body=kwargs)
end
"""
    create_application_command(c::Client, guild::Snowflake; kwargs...) -> ApplicationCommand

Creates a guild [`ApplicationCommand`](@ref).
"""
function create_application_command(c::Client, guild::Snowflake; kwargs...)
    appid = c.application_id
    return Response{ApplicationCommand}(c, :POST, "/applications/$appid/guilds/$guild/commands"; body=kwargs)
end

"""
    get_application_commands(c::Client) -> Vector{ApplicationCommand}

Gets all global [`ApplicationCommand`](@ref)s for the logged in client.
"""
function get_application_commands(c::Client)
    appid = c.application_id
    return Response{Vector{ApplicationCommand}}(c, :GET, "/applications/$appid/commands")
end
"""
    get_application_commands(c::Client, guild::Snowflake) -> Vector{ApplicationCommand}

Gets all guild [`ApplicationCommand`](@ref)s for the logged in client.
"""
function get_application_commands(c::Client, guild::Snowflake)
    appid = c.application_id
    return Response{Vector{ApplicationCommand}}(c, :GET, "/applications/$appid/guilds/$guild/commands")
end

"""
    respond_to_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...) -> Message

Respond to an interaction with code 4.
"""
function respond_to_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...)
    dict = Dict{Symbol, Any}(
        :data => kwargs,
        :type => 4,
    )
    return Response{Message}(c, :POST, "/interactions/$int_id/$int_token/callback"; body=dict)
end

function respond_to_interaction_with_a_modal(c::Client, int_id::Snowflake, int_token::String; kwargs...)
    dict = Dict{Symbol, Any}(
        :data => kwargs,
        :type => 9,
    )
    return Response{Message}(c, :POST, "/interactions/$int_id/$int_token/callback"; body=dict)
end

"""
    ack_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...)

Respond to an interaction with code 5.
"""
function ack_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...)
    dict = Dict{Symbol, Any}(
        :type => 5,
    )
    return Response(c, :POST, "/interactions/$int_id/$int_token/callback"; body=dict)
end

"""
    update_ack_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...) -> Message

Respond to an interaction with code 6.
"""
function update_ack_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...)
    dict = Dict{Symbol, Any}(
        :type => 6,
    )
    return Response{Message}(c, :POST, "/interactions/$int_id/$int_token/callback"; body=dict)
end

"""
    update_message_int(c::Client, int_id::Snowflake, int_token::String; kwargs...)

Respond to an interaction with code 7.
"""
function update_message_int(c::Client, int_id::Snowflake, int_token::String; kwargs...)
    dict = Dict{Symbol, Any}(
        :data => kwargs,
        :type => 7,
    )
    return Response(c, :POST, "/interactions/$int_id/$int_token/callback"; body=dict)
end

"""
    create_followup_message(c::Client, int_id::Snowflake, int_token::String; kwargs...) -> Message

Creates a followup message for an interaction.
"""
function create_followup_message(c::Client, int_id::Snowflake, int_token::String; kwargs...)
    appid = c.application_id
    return Response{Message}(c, :POST, "/webhooks/$appid/$int_token)"; body=kwargs)
end

"""
    edit_interaction(c::Client, int_id::Snowflake, int_token::String; kwargs...)

Edit a followup message for an interaction.
"""
function edit_interaction(c::Client, int_token::String, mid::Snowflake; kwargs...)
    appid = c.application_id
    return Response(c, :PATCH, "/webhooks/$appid/$int_token/messages/$mid"; body=kwargs)
end

"""
    bulk_overwrite_application_commands(c::Client, guild::Snowflake, cmds::Vector{ApplicationCommand}) -> Vector{ApplicationCommand}

Overwrites global [`ApplicationCommand`](@ref)s with the given cmds vector.
"""
function bulk_overwrite_application_commands(c::Client, guild::Snowflake, cmds::Vector{ApplicationCommand})
    appid = c.application_id
    return Response{Vector{ApplicationCommand}}(c, :PUT, "/applications/$appid/guilds/$guild/commands"; body=cmds)
end
"""
    bulk_overwrite_application_commands(c::Client, guild::Snowflake, cmds::Vector{ApplicationCommand}) -> Vector{ApplicationCommand}

Overwrites guild [`ApplicationCommand`](@ref)s with the given cmds vector.
"""
function bulk_overwrite_application_commands(c::Client, cmds::Vector{ApplicationCommand})
    appid = c.application_id
    return Response{Vector{ApplicationCommand}}(c, :PUT, "/applications/$appid/commands"; body=cmds)
end


text:
# Implementation of
#
# Áfra, Attila T., & Szirmay-Kalos, László (2013). Stackless
# multi-bvh traversal for CPU, mic and GPU ray tracing. Computer
# Graphics Forum. http://dx.doi.org/10.1111/cgf.12259
#

const AfraBVHAccelerator{T,O,B<:BVH{T,O,IntTree{T,<:Unsigned}}} = BVHAccelerator{T,O,B,:afra}

for (fun,R,getray) in [(:intersect,:Ray,:identity), (:intersect!,:Intersection,:(r->r.ray))]
    @eval function Base.$fun(ray::$R{T}, acc::AfraBVHAccelerator{T,O,B}) where {T,O,UI,B<:BVH{T,O,IntTree{T,O,UI}}}
        tree = acc.bvh.tree
        node = acc.bvh.root
        bitstack = zero(UI)
        while true
            if is_inner(node)
                lnode = tree[node.left]
                rnode = tree[node.right]
                left = intersect($getray(ray), aabb(lnode))
                right = intersect($getray(ray), aabb(rnode))
                if left || right
                    bitstack <<= 1
                    if left ⊻ right
                        node = left ? lnode : rnode
                    else
                        # node = nearest child
                        node = lnode
                        bitstack |= 1
                    end
                    continue
                end
            else # Leaf
                $fun(ray, node) && return true
            end
            # Backtrack
            while bitstack & 1 == 0
                bitstack == 0 && return false
                node = parent(node, tree)
                bitstack >>= 1
            end
            node = sibling(node, tree)
            bitstack ⊻= 1
        end
        false
    end
end

export AfraBVHAccelerator


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SoapyPlutoSDR_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SoapyPlutoSDR")
JLLWrappers.@generate_main_file("SoapyPlutoSDR", UUID("96b01697-91ed-557a-a63b-0996cb7d9743"))
end  # module SoapyPlutoSDR_jll


text:
# obtain the encoding of the page
using Thrift: isfilled

function page_encoding(page::Page)
    if isfilled(page.hdr, :data_page_header)
        return page.hdr.data_page_header.encoding
    elseif isfilled(page.hdr, :data_page_header_v2)
        return page.hdr.data_page_header_v2.encoding
    elseif isfilled(page.hdr, :dictionary_page_header)
        return page.hdr.dictionary_page_header.encoding
    else
        error("not supported page")
    end
end


text:
module PrincipalMomentAnalysisApp

export pmaapp

using PrincipalMomentAnalysis
using LinearAlgebra
using DataFrames
using CSV
using Statistics

using Blink
using JSExpr

using Colors
using PlotlyJS
import IterTools
using SparseArrays

include("Schedulers.jl")
using .Schedulers

include("pmaplots.jl")
include("pca.jl")
include("app.jl")

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/concave_increasing.jl
# Contains definitions of functions used to enforce bounds on intermediate
# terms in the computation of relaxations:
# positive, negative, lower_bnd, upper_bnd, bnd
#############################################################################

function cut_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    zL = z.lo
    zU = z.hi
    cvv, ccv, cv_grad, cc_grad = cut(zL, zU, x.cv, x.cc, x.cv_grad, x.cc_grad)
    return MC{N, T}(cvv, ccv, z, cv_grad, cc_grad, x.cnst)
end

"""
positive(x::MC)

Sets the lower interval bound and the convex relaxation of `x` to a value of
at least `McCormick.MC_DOMAIN_TOL`. (Sub)gradients are adjusted appropriately.
"""
positive(x) = x
function positive(x::Interval)
     x ∩ Interval(MC_DOMAIN_TOL, Inf)
end
function positive_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end
function positive(x::MC{N,T}) where {N, T<:Union{NS,MV}}
    positive_kernel(x, positive(x.Intv))
end

negative(x) = x
function negative(x::Interval)
     x ∩ Interval(-Inf, -MC_DOMAIN_TOL)
end
function negative_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
negative(x::MC)

Sets the upper interval bound and the concave relaxation of `x` to a value of
at most `-McCormick.MC_DOMAIN_TOL`. (Sub)gradients are adjusted appropriately.
"""
function negative(x::MC{N,T}) where {N, T<:Union{NS,MV}}
    negative_kernel(x, negative(x.Intv))
end

lower_bnd(x, lb) = x
function lower_bnd(x::Interval{Float64}, lb::Float64)
     x ∩ Interval(lb, Inf)
end
function lower_bnd_kernel(x::MC{N,T}, lb::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
lower_bnd(x::MC, lb::Float64)

Sets the lower interval bound and the convex relaxation of `x` to a value of
at least `lb`. (Sub)gradients are adjusted appropriately.
"""
function lower_bnd(x::MC{N,T}, lb::Float64) where {N, T<:Union{NS,MV}}
    lower_bnd_kernel(x, lb, lower_bnd(x.Intv, lb))
end

upper_bnd(x, lb) = x
function upper_bnd(x::Interval{Float64}, ub::Float64)
     x ∩ Interval(-Inf, ub)
end
function upper_bnd_kernel(x::MC{N,T}, ub::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
upper_bnd(x::MC, ub)

Sets the upper interval bound and the concave relaxation of `x` to a value of
at most `ub`. (Sub)gradients are adjusted appropriately.
"""
function upper_bnd(x::MC{N,T}, lb::Float64) where {N, T<:Union{NS,MV}}
    upper_bnd_kernel(x, lb, upper_bnd(x.Intv, lb))
end

bnd(x, lb, ub) = x
function bnd(x::Interval{Float64}, lb::Float64, ub::Float64)
     x ∩ Interval(lb, ub)
end
function bnd_kernel(x::MC{N,T}, lb::Float64, ub::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    cut_kernel(x, z)
end

"""
bnd(x::MC, lb, ub)

Sets the lower interval bound and the convex relaxation of `x` to a value of
at least `lb`. Sets the upper interval bound and the concave relaxation of `x`
to a value of at most `ub`. (Sub)gradients are adjusted appropriately.
"""
function bnd(x::MC{N,T}, lb::Float64, ub::Float64) where {N, T<:Union{NS,MV}}
    bnd_kernel(x, lb, ub, bnd(x.Intv, lb, ub))
end

function d_lower_bnd_grad(g, x::T, y::T) where T<:Number
    g[1] = one(T)
    g[2] = zero(T)
    return
end
function d_upper_bnd_grad(g, x::T, y::T) where T<:Number
    g[1] = one(T)
    g[2] = zero(T)
    return
end
function d_bnd_grad(g, x::T, y::T, z::T) where T<:Number
    g[1] = one(T)
    g[2] = zero(T)
    g[3] = zero(T)
    return
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
struct DataDistribution
    dd::Dict{Float64, Distributions.Distribution}
end

# assuming normally distributed data
function createDataDistribution(dfs::Vector{DataFrame}, obs::Vector{Symbol})
    nobs = length(obs)

    # extract unique timepoints
    times = Float64[]
    for df in dfs
        append!(times, convert(Vector{Float64}, df[:Time]))
    end
    times = sort!(unique(times))

    # create vector of data distributions over time
    npoints = length(times)
    dd = Dict{Float64, Distributions.Distribution}()
    for t in times
        ar = Float64[]
        for df in dfs
            if any(df[:Time] .== t)
                ar = vcat(ar, convert(Array{Float64}, df[df[:Time] .== t, obs]))
            end
        end
        m = vec(mean(ar, 1))
        s = vec(std(ar, 1))

        # if std = NaN? Implement 10% std as a rough approximation
        if any(s .== NaN)
            s[s .== NaN] = m[s .== NaN] ./ 10
        end
        if any(s .== 0.0)
            s[s .== 0.0] = m[s .== 0.0] ./ 10
        end

        dd[Float64(t)] = MvNormal(m, s)
    end
    return DataDistribution(dd), convert(Vector{Float64},times)
end

function createDataDistributionM(dfs::Vector{DataFrame}, obs::Vector{Symbol})
    nobs = length(obs)

    # extract unique timepoints
    times = Float64[]
    for df in dfs
        append!(times, convert(Vector{Float64}, df[!, :Time]))
    end
    times2 = sort!(unique(times))

    # create vector of data distributions over time
    npoints = length(times2)
    dd = Dict{Float64, Distributions.Distribution}()
    for t in times2
        ar = MvNormal[]
        for df in dfs
            if any(df[!, :Time] .== t)
                for aa in 1:size(df[df[!, :Time] .== t, obs])[1]
                    vv = vec(convert(Array{Float64},df[df[!, :Time] .== t, obs][aa,:]))
                    push!(ar,MvNormal(vv, vv./10))
                end
            end
        end
        dd[Float64(t)] = MixtureModel(vec(ar))
    end
    return DataDistribution(dd), times2
end

function genstate(d::DataDistribution, t::Float64)
    xa = round.(Int,rand(d.dd[t], 1))
    while any(xa.<0)
        xa = round.(Int,rand(d.dd[t], 1))
    end
    vec(xa)
end

function loadData(path::String)
    datasets = DataFrame[] # Init empty vector
    if ispath(path)
        if isfile(path)
            push!(datasets, CSV.read(path))
            return datasets
        else
            files = readdir(path) # grab all files in the given folder
            for f in files
                push!(datasets, CSV.read(path * "/" * f))
            end
            return datasets
        end
    else
        error("Enter a valid file or folder")
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
using MAT
using PyPlot

SMALL_SIZE = 22
MEDIUM_SIZE = 24
BIGGER_SIZE = 24

plt.rc("font", size=SMALL_SIZE)          # controls default text sizes
plt.rc("axes", titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc("axes", labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc("xtick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc("ytick", labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc("legend", fontsize=SMALL_SIZE)    # legend fontsize
plt.rc("figure", titlesize=BIGGER_SIZE)  # fontsize of the figure title

loss = matread("sess_figures/ex1/loss_nn/loss1.mat")["L"]

for k = 2:50
    loss_ = matread("sess_figures/ex1/loss_nn/loss$k.mat")["L"]
    global loss = [loss; loss_]
end

close("all"); figure();semilogy(loss); tight_layout()
savefig("sess_figures/ex1/nn/loss.png")

loss2 = matread("sess_figures/ex1/loss_ones/loss1.mat")["L"]

for k = 2:50
    loss_ = matread("sess_figures/ex1/loss_ones/loss$k.mat")["L"]
    global loss2 = [loss2; loss_]
end

close("all"); figure();semilogy(loss2); tight_layout()
savefig("sess_figures/ex1/ones/loss.png")

close("all"); figure();
semilogy(loss, lw=2, label="DNN");
semilogy(loss2, lw=2, label="pointwise"); legend();
xlabel("number of iterations");
ylabel("loss")
tight_layout()
savefig("sess_figures/ex1/loss.png")


text:
using Blosc
using Base.Test

@test_throws ArgumentError Blosc.set_num_threads(0)
@test_throws ArgumentError Blosc.set_num_threads(Blosc.MAX_THREADS + 1)
@test Blosc.set_num_threads(1) == 1

@test_throws ArgumentError Blosc.set_compressor("does_not_exist")
for comp in Blosc.compressors()
    Blosc.set_compressor(comp)
end
@test "blosclz" in Blosc.compressors()
@test Blosc.set_compressor("blosclz") != -1

Blosc.set_blocksize(16)
Blosc.set_blocksize(256)
Blosc.set_blocksize(0)
Blosc.set_blocksize()
@test_throws ArgumentError Blosc.set_blocksize(-1)

s = convert(String, rand('0':'z', 10000))
@test String(decompress(UInt8, compress(s))) == s
@test isempty(decompress(UInt8, compress("")))

x = rand(100)
@test decompress(eltype(x), compress(x)) == x
@test isempty(decompress(Int, compress(Int[])))

# round trip test
roundtrip(orig) = Blosc.decompress(eltype(orig), Blosc.compress(orig)) == orig
for ty in [Float16, Float32, Float64,
           Int8, Int16, Int32, Int64, Int128,
           UInt8, UInt16, UInt32, UInt64, UInt128]
    for i=1:2048
        a = rand(ty, i)
        @test roundtrip(a)
    end
end
# cannot compress element types that are not isbits
@test_throws ArgumentError Blosc.compress([BigInt(1)])

# test that we actually are compressing
a = ones(Float64, 1000)
ac = Blosc.compress(a)
@test sizeof(ac) < sizeof(a)
@test Blosc.decompress(Float64, ac) == a

# test all compressors
for (comp, name, _) in Blosc.compressors_info()
    Blosc.set_compressor(comp)
    for level=0:9
        for shuffle in (true, false)
            for i=1:2048
                a = rand(UInt8, i)
                ac = Blosc.compress(a, level=level, shuffle=shuffle)
                info = Blosc.compressor_info(ac)
                @test info.library == name
                @test info.typesize == sizeof(UInt8)
                @test info.shuffled == shuffle
                @test Blosc.decompress(UInt8, ac) == a
            end
        end
    end
end

# test compress invalid args
@test_throws ArgumentError Blosc.compress([BigInt(1)])
@test_throws ArgumentError Blosc.compress(ones(UInt8, 256), level=-1)
@test_throws ArgumentError Blosc.compress(ones(UInt8, 256), level=11)

@test Blosc.free_resources!()


text:
module TestTlsiid

using Distributions
using TimeData
using Base.Test

import Econometrics

## reload(joinpath(Pkg.dir(), "Econometrics", "src", "Econometrics.jl"))

############
## TlsIID ##
############

## constructors / display
##-----------------------

issubtype(Econometrics.TlsIID, Econometrics.AbstrUnivarModel)

@test_throws Exception Econometrics.TlsIID(-0.3, -0.2, 0.1)
@test_throws Exception Econometrics.TlsIID(0.3, -0.2, -0.1)

mod = Econometrics.TlsIID(0.3, -0.2, 0.1)

println("\nTesting display output:")
display(mod)
Econometrics.description(mod)

## getParams
##----------

mod = Econometrics.TlsIID(2.3, 0.2, 1.4)
@test Econometrics.getParams(mod) == (2.3, 0.2, 1.4)

## simulate
##---------

nSim = 5000
mod = Econometrics.TlsIID(2.3, 0.2, 1.4)
simVals = Econometrics.simulate(mod, nSim)

## estimate
##---------

mod_fit = Econometrics.estimate(Econometrics.TlsIID, simVals)

## compare parameters
[[Econometrics.getParams(mod)...]'
[Econometrics.getParams(mod_fit)...]']

end


text:
module SuffixArrays

export suffixsort

immutable SuffixArray{S<:AbstractString,N<:Signed}
    string::S
    n::Int
    index::Array{N,1}
end

function SuffixArray{S}(s::S)
    n = length(s)
    index = zeros(n < 256 ? Int8 :
               n < 65536 ? Int16 :
               n < 4294967296 ? Int32 : Int64,n)
    return SuffixArray(s,n,index)
end

include("sais.jl")

function suffixsort(s)
    isempty(s) && return SA
    SA = SuffixArray(s)
    SA.n <= 1 && return SA
    SuffixArrays.sais(s, SA.index, 0, SA.n, isascii(s) ? 256 : 65536, false)
    return SA
end

#=contains(haystack, needle)

matchall(substring, s::String)=#
const MAXCHAR = Char(255)

function lcp2(SA,s)
    inv = similar(SA)
    lcparr = similar(SA)
    n = length(SA)
    for i = 1:n
        inv[SA[i]+1] = i-1
    end
    m = 0
    for i = 1:n
        if inv[i] > 0
            j = SA[inv[i]]
            while s[m+i] == s[m+j+1]
                m += 1
            end
            lcparr[inv[i]+1] = m
            m > 0 && (m-=1)
        end
    end
    lcparr[1] = -1
    return lcparr
end

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
include("../common.jl")

gap_version = v"400.1192.000"
gap_lib_version = v"400.1192.000"
name = "Browse"
upstream_version = v"1.8.13" # when you increment this, reset offset to v"0.0.0"
offset = v"0.0.1" # increment this when rebuilding with unchanged upstream_version, e.g. gap_version changes
version = offset_version(upstream_version, offset)

# Collection of sources required to build libsingular-julia
sources = [
    ArchiveSource("https://www.math.rwth-aachen.de/~Browse/$(name)-$(upstream_version).tar.bz2",
                  "118d1f7fb67044c103cf44673ade8de57e4289a906dd01e15d7562dd0aface5b"),
]

# Bash recipe for building across all platforms
script = raw"""
cd Browse*

./configure ${prefix}/share/gap
make -j${nproc} CFLAGS="-I$includedir -I$includedir/ncurses"

# copy the loadable module
mkdir -p ${prefix}/lib/gap
cp bin/*/*.so ${prefix}/lib/gap/

install_license /usr/share/licenses/GPL3
"""

name = gap_pkg_name(name)
platforms, dependencies = setup_gap_package(gap_version, gap_lib_version)

# The products that we will ensure are always built
products = [
    FileProduct("lib/gap/ncurses.so", :ncurses),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: stemann/Yggdrasil

max_issues_repo_name: stemann/Yggdrasil

max_forks_repo_name: stemann/Yggdrasil


text:
using NESSie.Format

@testset "readstl" begin
    @test_skip readstl
end

@testset "writestl" begin
    for T in testtypes
        # empty system
        model = Model{T, Triangle{T}}()
        (fname, fh) = mktemp()
        try
            writestl(fh, model)
            seek(fh, 0)
            @test [read(fh, UInt8) for _ in 1:80] ≈ zeros(Float32, 80)
            @test read(fh, UInt32) == 0
            @test_throws EOFError read(fh, UInt8)
        finally
            close(fh)
            rm(fname)
        end

        # small system
        nodes = Vector{T}[[0, 0, 0], [0, 0, 3], [0, 3, 0], [1, -3, 3]]
        model = Model(
            nodes,
            Triangle{T}[
                Triangle(nodes[1], nodes[2], nodes[3]),
                Triangle(nodes[1], nodes[4], nodes[2])
            ]
        )
        (fname, fh) = mktemp()
        try
            writestl(fh, model)
            seek(fh, 0)
            @test [read(fh, UInt8) for _ in 1:80] ≈ zeros(Float32, 80)
            @test read(fh, UInt32) == 2
            for elem in model.elements
                @test [read(fh, Float32) for _ in 1:3] ≈ elem.normal
                @test [read(fh, Float32) for _ in 1:3] ≈ elem.v1
                @test [read(fh, Float32) for _ in 1:3] ≈ elem.v2
                @test [read(fh, Float32) for _ in 1:3] ≈ elem.v3
                @test read(fh, UInt16) == zero(UInt16)
            end
            @test_throws EOFError read(fh, UInt8)
        finally
            close(fh)
            rm(fname)
        end
    end
end


text:
#=
Neighbor is an element of a patch.

Element of array of points of the target image (both selection and context).
But kept as offset (arrow vectors) i.e. that work in the target and corpus frames.
That is, the coordinate system of an offset has as origin the distinguished point of the patch.

!!! Mutable.  A single vector of Neighbors is mutated.
(When created dynamically, need not be mutable.)

!!! All fields are NOT abstract types, since that lowers performance.
Fields are parameterized concrete types.
=#

include("null.jl")  # leastIndex
include("image.jl")

mutable struct Neighbor{ValueType, DimensionCount}

    #=
    The fundamental attribute: a vector from patch's distinguished point to this neighbor point.

    Not a PointInMaskedImage, i.e. doesn't carry a frame.
    =#
    offset::CartesianIndex{DimensionCount}

    #=
    Attributes for performance.
    1) attributes used many times computed once.
       Same targetPoint is used for many search probes over the corpus
    2) CPU cache coherence is improved if we precompute these attributes into one compact locality

    For example, the cache can hold the entire ScatterPatch,
    but might not hold the entire target image.
    Since on the first pass the target points might be scattered widely,
    and even on further passes, since the image is not tiled, nearby points
    can be far away in the cache memory.

    TODO the original code kept: targetPoint, color, sourceOf
    =#

    #=
    The corresponding point in the target image.
    Assert is inbounds of the target image.
    =#
    targetPoint::CartesianIndex{DimensionCount}

    # The value (e.g. color) of the point in the target
    value::ValueType

    # TODO
    # Coords of corpus point this target synthed from, or -1 if this neighbor is context
    # Coordinates sourceOf;
end

#=
Outer constructor taking only targetImage.
Creates more or less undefined Neighbor, subsequently mutated.

The types must be correct, but the values are don't care: we write but never read.
So we use the coordinates and values of the first index.
=#

function Neighbor(targetImage::MaskedImage{ValueType, DimensionCount}) where {ValueType, DimensionCount}
    firstIndex = leastIndex(targetImage.image)
    # WAS: CartesianIndex(1,1) but that is specialized to 2D
    return Neighbor(
        firstIndex, # point in target frame
        firstIndex, # point in corpus frame ??? is leastIndex same for corpus?
        targetImage.image[firstIndex]   # value
        )
end

#=
Outer constructor taking only the offset, targetImage, and targetPoint (and not the value.)
The value is taken from the targetImage.
=#
function Neighbor(
        offset::CartesianIndex{DimensionCount},
        targetMaskedImage::MaskedImage{ValueType, DimensionCount},
        targetPoint::CartesianIndex{DimensionCount}    # framed point
        ) where {ValueType, DimensionCount}
    # call default constructor
    return Neighbor(
        offset,
        targetPoint,
        targetMaskedImage.image[targetPoint]    # value
        )
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
@testset "bigfloat initial convergence #720" begin
    f(x) = x[1]^2
    x0 = BigFloat[0]
    obj = OnceDifferentiable(f, x0; autodiff=:forward)
    for method in (GradientDescent, BFGS, LBFGS, AcceleratedGradientDescent, MomentumGradientDescent, ConjugateGradient, Newton, NewtonTrustRegion, SimulatedAnnealing)
        result = Optim.optimize(f, x0, method())
    end
end


text:
"""
    promote_tspan(tspan)

Convert the `tspan` field of a `DEProblem` to a `(tmin, tmax)` tuple, where both
elements are of the same type. If `tspan` is a function, returns it as-is.
"""
promote_tspan((t1,t2)::Tuple{T,S}) where {T,S} = promote(t1, t2)
promote_tspan(tspan::Number) = (zero(tspan),tspan)
promote_tspan(tspan::Nothing) = (nothing,nothing)
promote_tspan(tspan::Function) = tspan
promote_tspan(tspan::AbstractArray) = length(tspan) == 2 ? promote_tspan((first(tspan),last(tspan))) : throw(error("The length of tspan must be two (and preferably, tspan should be a tuple, i.e. (0.0,1.0)). If you are trying to include other values for saving reasons, note see the [common solver arguments page](https://docs.juliadiffeq.org/latest/basics/common_solver_opts/) for information on the saving command saveat."))

### Displays

Base.summary(prob::DEProblem) = string(TYPE_COLOR, nameof(typeof(prob)),
                                       NO_COLOR, " with uType ",
                                       TYPE_COLOR, typeof(prob.u0),
                                       NO_COLOR, " and tType ",
                                       TYPE_COLOR,
                                       typeof(prob.tspan) <: Function ?
                                       "Unknown" : (prob.tspan === nothing ?
                                       "Nothing" : typeof(prob.tspan[1])),
                                       NO_COLOR, ". In-place: ",
                                       TYPE_COLOR, isinplace(prob),
                                       NO_COLOR)

Base.summary(prob::AbstractLinearProblem) = string(TYPE_COLOR, nameof(typeof(prob)),
                                                   NO_COLOR, ". In-place: ",
                                                   TYPE_COLOR, isinplace(prob),
                                                   NO_COLOR)
function Base.show(io::IO, A::AbstractLinearProblem)
  show(io,summary(A))
  print(io,"b: ")
  show(io, A.b)
end

Base.summary(prob::AbstractNonlinearProblem{uType,iip}) where {uType,iip} = string(
                                       TYPE_COLOR, nameof(typeof(prob)),
                                       NO_COLOR, " with uType ",
                                       TYPE_COLOR, uType,
                                       NO_COLOR, ". In-place: ",
                                       TYPE_COLOR, isinplace(prob),
                                       NO_COLOR)
function Base.show(io::IO, A::AbstractNonlinearProblem)
  println(io,summary(A))
  print(io,"u0: ")
  show(io, A.u0)
end

Base.summary(prob::AbstractOptimizationProblem) = string(
                                       TYPE_COLOR, nameof(typeof(prob)),
                                       NO_COLOR, ". In-place: ",
                                       TYPE_COLOR, isinplace(prob),
                                       NO_COLOR)
function Base.show(io::IO, A::AbstractOptimizationProblem)
  println(io,summary(A))
  print(io,"u0: ")
  show(io, A.u0)
end

Base.summary(prob::AbstractQuadratureProblem) = string(
                                                       TYPE_COLOR, nameof(typeof(prob)),
                                                       NO_COLOR, ". In-place: ",
                                                       TYPE_COLOR, isinplace(prob),
                                                       NO_COLOR)
function Base.show(io::IO, A::AbstractQuadratureProblem)
  println(io,summary(A))
end

Base.summary(prob::AbstractNoiseProblem) = string(nameof(typeof(prob))," with WType ",typeof(prob.noise.W[1])," and tType ",typeof(prob.tspan[1]),". In-place: ",isinplace(prob))
function Base.show(io::IO, A::DEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  println(io)
  print(io,"u0: ")
  show(io, A.u0)
end
function Base.show(io::IO, A::AbstractNoiseProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  println(io)
end
function Base.show(io::IO, A::AbstractDAEProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  println(io)
  print(io,"u0: ")
  show(io, A.u0)
  println(io)
  print(io,"du0: ")
  show(io, A.du0)
end

Base.summary(prob::AbstractEnsembleProblem) = string(
nameof(typeof(prob))," with problem ",nameof(typeof(prob.prob)))
Base.show(io::IO, A::AbstractEnsembleProblem) = print(io,summary(A))
TreeViews.hastreeview(x::DEProblem) = true
function TreeViews.treelabel(io::IO,x::DEProblem,
                             mime::MIME"text/plain" = MIME"text/plain"())
  show(io,mime,Base.Text(Base.summary(x)))
end

struct NullParameters end
Base.getindex(::NullParameters,i...) = error("Parameters were indexed but the parameters are `nothing`. You likely forgot to pass in parameters to the DEProblem!")
Base.iterate(::NullParameters) = error("Parameters were indexed but the parameters are `nothing`. You likely forgot to pass in parameters to the DEProblem!")

function Base.show(io::IO, A::AbstractPDEProblem)
  println(io,summary(A.prob))
  println(io)
end
Base.summary(prob::AbstractPDEProblem) = string(TYPE_COLOR,
                                                nameof(typeof(prob)),
                                                NO_COLOR)


text:
@testset "Wagner + Section" begin

    coupling = Coupling(Wagner(), Section())

    # test provided jacobians
    coupling_jacobian_tests(coupling)

    # test consistency of input/output functions
    coupling_io_tests(coupling)
end


text:
using OneDmin
using Test
using Printf
using Logging

fh(x) = - ( x + sin(x)) * exp(-x^2)
x0 = big(0.2); x1 = big(1.0)

using ForwardDiff

gh(x) = ForwardDiff.derivative(fh, x)  # first derivative
Hh(x) = ForwardDiff.derivative(gh, x)  # second derivative
Ch(x) = ForwardDiff.derivative(Hh, x)  # third derivative

setprecision(4096)

function test(algo::Function, x, xm)
    h = fh(x)
    dh = gh(x)
    d2h = Hh(x)
    d3h = Ch(x)

    hm = fh(xm)
    dhm = gh(xm)
    d2hm = Hh(xm)

    it = 0
    maxit = 25
    with_logger(Logging.NullLogger()) do
        @info "it  x  |h'|"
        while abs(dh) > BigFloat("1.0e-1200")
            xp = algo(x,xm,h,hm,dh,dhm,d2h,d2hm,d3h,0)
            x, xm = xp, x

            hm = h
            dhm = dh
            d2hm = d2h

            h = fh(x)
            dh = gh(x)
            d2h = Hh(x)
            d3h = Ch(x)

            it += 1
            if (it > maxit)  break  end
            @info @sprintf("%2d",it), @sprintf("%f",x), @sprintf("%.2e",abs(dh))
        end
    end
    return it

end

iter1 = Dict(base .=> [16,12,8])
iter2 = Dict(quintic .=> [8,7])
iter3 = Dict(cubic1 .=> [8,8,11])
iter4 = Dict(cubic0 .=> [11,11,15])

iter = merge(iter1, iter2, iter3, iter4)

@testset "OneDmin.jl" begin
    xm = x1; x = x0

    @info "Testing algos in base."
    for algo in base
        @info "testing" algo
        it = test(algo, x, xm)
        @test it == iter[algo]
    end

    @info "Testing algos in quintic."
    for algo in quintic
        @info "testing" algo
        it = test(algo, x, xm)
        @test it == iter[algo]
    end

    @info "Testing algos in cubic1."
    for algo in cubic1
        @info "testing" algo
        it = test(algo, x, xm)
        @test it == iter[algo]
    end

    @info "Testing algos in cubic0."
    for algo in cubic0
        @info "testing" algo
        it = test(algo, x, xm)
        @test it == iter[algo]
    end

    include("test1D.jl")

end


text:
# https://discourse.julialang.org/t/write-a-rest-interface-like-flask/18538

# Bukdu v0.4.2
using Bukdu

struct RestController <: ApplicationController
    conn::Conn
end

function init(c::RestController)
    render(JSON, (:init, c.params.region, c.params.site_id, c.params.channel_id))
end

function update(c::RestController)
    render(JSON, (:update, c.params.region, c.params.site_id, c.params.channel_id))
end

routes() do
    get("/init/region/:region/site/:site_id/channel/:channel_id/", RestController, init, :site_id=>Int, :channel_id=>Int)
    get("/update/region/:region/site/:site_id/channel/:channel_id/", RestController, update, :site_id=>Int, :channel_id=>Int)
end

Bukdu.start(8080)

#=
curl localhost:8080/init/region/west/site/1/channel/2/
curl localhost:8080/update/region/west/site/1/channel/2/
=#


text:
"""
    elimGaussNonPivoting(A,B)
adalah fungsi untuk menyelesaikan SPL `Ax=B` dengan `A` adalah matriks sembarang
dan `B` adalah vektor ruas kanan menggunakan metode eliminasi Gauss tanpa pivoting.

# Examples
```jl
julia> A = [ 2 5 4 2 3
            4 1 0.4 6 5
            0.9 4 8 7 10
            5 7 0.5 2 7
            2 6 8 9 4];

julia> B = [26; 24.4; 32.8; 30.5; 53];

julia> Xn = elimGaussNonPivoting(A,B)
5×1 Array{Float64,2}:
  2.0000000000000018
  2.999999999999999
  1.0000000000000013
  2.999999999999999
 -1.0000000000000009
```
return solusi SPL `X`
"""
function elimGaussNonPivoting(A,b)
    # Hitung ukuran matriks dan inisiasi solusi X
    n = length(b);
    X = zeros(n,1);
    # Buatlah matriks Gandeng [A|b]
    Aug = [A b];
    # Lakukan operasi baris dasar terhadap matriks gandeng
    for p = 1:n-1 # p menunjukkan kolom
        # Jika pivot bernilai nol, maka gagal.
        if Aug[p,p]==0
            error("Harus Pakai Pivoting")
        end
        # Lakukan eliminasi menggunakan operasi Eij(k)
        for i = p+1:n # i menunjukkan baris
            k = Aug[i,p]/Aug[p,p];
            Aug[i,:] = Aug[i,:] - k*Aug[p,:];
        end
    end
    # Pisahkan matriks gandeng, kemudian substitusi mundur.
    A = Aug[:,1:n];
    b = Aug[:,1+n];
    X = backsub(A,b);
end


text:
get_pre(::Synapse, sys::ODESystem) = sys.V
get_post(::Synapse, sys::ODESystem) = sys.V

my_pre(::Synapse, sys::ODESystem) = sys.Vpre
my_post(::Synapse, sys::ODESystem) = sys.Vpost
post_connector(::Synapse) = :Isyn

mutable struct Chol{T} <: Synapse
    ḡChol::T
    s::T
    Eₛ::T # mV
    k₋::T # ms
    Vth::T # mV
    δ::T   # mV
end
Chol(x) = Chol(x,0.,-80., 0.01, -35., 5.)
s̄(syn::Chol, Vpre)  = 1. / (1. + exp((syn.Vth - Vpre)/syn.δ))
τs(syn::Chol, Vpre) = (1. - s̄(syn, Vpre))/syn.k₋
ionic_current(::Chol, sys::ODESystem) = sys.IChol

mutable struct Glut{T} <: Synapse
    ḡGlut::T
    s::T
    Eₛ::T
    k₋::T
    Vth::T
    δ::T
end
Glut(x) = Glut(x,0.,-70., 0.025, -35., 5.)
s̄(syn::Glut, Vpre)  = 1. / (1. + exp((syn.Vth - Vpre)/syn.δ))
τs(syn::Glut, Vpre) = (1. - s̄(syn, Vpre))/syn.k₋
ionic_current(::Glut, sys::ODESystem) = sys.IGlut

function channel_dynamics(ch::Chol, Vpre, Vpost)
    states = @variables s(t) IChol(t)
    parameters = @parameters ḡChol
    eqs = [ D(s)  ~     (1/τs(ch, Vpre))*(s̄(ch, Vpre) - s),
            IChol ~     -ḡChol*s*(Vpost - ch.Eₛ)]
    current = [eqs[2]]
    defaultmap = [s => ch.s, ḡChol => ch.ḡChol]
    return eqs, states, parameters, current, defaultmap
end

function channel_dynamics(ch::Glut, Vpre, Vpost)
    states = @variables s(t) IGlut(t)
    parameters = @parameters ḡGlut
    eqs = [ D(s) ~ (1/τs(ch, Vpre))*(s̄(ch, Vpre) - s),
            IGlut ~ -ḡGlut*s*(Vpost - ch.Eₛ)]
    current = [eqs[2]]
    defaultmap = [s => ch.s, ḡGlut => ch.ḡGlut]
    return eqs, states, parameters, current, defaultmap
end


text:
include("preliminaries.jl")

const TIMEROUTPUT = TimerOutput()

@timed_testset "context" begin
    include("context.jl")
end

@timed_testset "cache" begin
    include("cache.jl")
end

@timed_testset "mimeoverride" begin
    include("mimeoverride.jl")
end

@timed_testset "with_terminal" begin
    include("with_terminal.jl")
end

@timed_testset "html" begin
    include("html.jl")
end

@timed_testset "style" begin
    include("style.jl")
end

@timed_testset "build" begin
    include("build.jl")
end

show(TIMEROUTPUT; compact=true, sortby=:firstexec)