text:
function to_256_colors(crayon::Crayon)
    fg = crayon.fg
    bg = crayon.bg
    crayon.fg.style == COLORS_24BIT && (fg = to_256_colors(crayon.fg))
    crayon.bg.style == COLORS_24BIT && (bg = to_256_colors(crayon.bg))
    return Crayon(
        fg,
        bg,
        crayon.reset,
        crayon.bold,
        crayon.faint,
        crayon.italics,
        crayon.underline,
        crayon.blink,
        crayon.negative,
        crayon.conceal,
        crayon.strikethrough,
    )
end

function to_256_colors(color::ANSIColor)
    @assert color.style == COLORS_24BIT
    r, g, b = color.r, color.g, color.b
    r24, g24, b24 = map(c->round(Int, c * 23 / 256), (r, g, b))
    if r24 == g24 == b24
        return ANSIColor(UInt8(232 + r24), COLORS_256, color.active)
    else
        r6, g6, b6 = map(c->round(Int, c * 5  / 256), (r, g, b))
        return ANSIColor(UInt8(16 + 36 * r6 + 6 * g6 + b6), COLORS_256, color.active)
    end
end

# 24bit -> 16 system colors
function to_system_colors(crayon::Crayon)
    fg = crayon.fg
    bg = crayon.bg
    crayon.fg.style == COLORS_24BIT && (fg = to_system_colors(crayon.fg))
    crayon.bg.style == COLORS_24BIT && (bg = to_system_colors(crayon.bg))
    return Crayon(
        fg,
        bg,
        crayon.reset,
        crayon.bold,
        crayon.faint,
        crayon.italics,
        crayon.underline,
        crayon.blink,
        crayon.negative,
        crayon.conceal,
        crayon.strikethrough,
    )
end

function compute_value(r, g, b)
    r′, g′, b′ = (r, g, b) ./ 255
    Cmax = max(r′, g′, b′)
    return Cmax * 100
    #=
    # This is not needed
    Cmin = min(r′, g′, b′)
    Δ = Cmax - Cmin
    H = begin
        if Cmax == r′
            60 * (((g′ - b′) / Δ) % 6)
        elseif Cmax == g′
            60 * ((b′ - r′) / Δ + 2)
        else
            60 * ((r′ - g′) / Δ + 4)
        end
    end

    S = Cmax == 0 ? 0 : (Δ / Cmax)
    V = Cmax
    return H * 360, S * 100, V * 100
    =#
end

function to_system_colors(color::ANSIColor)
    @assert color.style == COLORS_24BIT
    r, g, b = color.r, color.g, color.b
    value = compute_value(r, g, b)

    value = round(Int, value / 50)

    if (value == 0)
        ansi = 0
    else
        ansi =
            ((round(Int, b / 255) << 2) |
             (round(Int, g / 255) << 1) |
              round(Int, r / 255))
        value == 2 && (ansi += 60)
    end
    return ANSIColor(UInt8(ansi), COLORS_16, color.active)
end


text:
using BasisFunctions
using Common, StaticArrays

@testset "BasisFunctions.jl" begin
    μ = Gaussian(42, [1.,2.,3.])

    @test μ.R ≈ [1,2,3]
    @test μ.ζ ≈ 42

    @test_throws TypeError Gaussian(1,[0,0,0])
    @test_throws DimensionMismatch Gaussian(1,[0.,0.])
end

ext: jl


text:
# stores all datafile edit information that pertains to all simulations

### UPDATE FUNCTION FOR DATA FILE ARRAYS

function UpdateArray(data_dictionary,string,alterations=[])
  # access array
  data_array = data_dictionary[string]
  # set modifications for all conditions here
  for i=1:size(alterations)
    updates = alterations[i]
    index = updates[1]
    value = updates[2]
    data_array[index]=value;
  data_dictionary[string] = data_array
  return data_dictionary
end


text:
import JET.VSCode:
    get_reports,
    gen_postprocess,
    vscode_diagnostics

@testset "sources" begin
    let res = @report_call sum("julia")
        @test res.source == "@report_call sum(::String)"
    end
    let res = @report_call sum(1:100)
        @test res.source == "@report_call sum(::$(typeof(1:100)))"
    end
    let
        local filename
        res = mktemp() do path, io
            filename = path
            report_file2(filename)
        end
        @test occursin(filename, res.source)
    end
end

@testset "diagnostics" begin
    hasfield′(obj::T, sym) where T = hasfield(T, sym)

    function check_basic_integration(diagnostics, reports)
        @test hasfield′(diagnostics, :source)
        @test hasfield′(diagnostics, :items)
        @test length(diagnostics.items) == length(reports)
    end
    function check_inference_integration(item, report)
        @test hasfield′(item, :msg)
        @test hasfield′(item, :path)
        @test hasfield′(item, :line)
        @test hasfield′(item, :severity)
        @test hasfield′(item, :relatedInformation)
        @test length(item.relatedInformation) == length(report.vst)
        @test !isempty(item.relatedInformation)
        ri = first(item.relatedInformation)
        @test hasfield′(ri, :msg)
        @test hasfield′(ri, :path)
        @test hasfield′(ri, :line)
    end
    function check_toplevel_integration(item, report)
        @test hasfield′(item, :msg)
        @test hasfield′(item, :path)
        @test hasfield′(item, :line)
        @test hasfield′(item, :severity)
    end

    # basic case
    let res = @report_call sum("julia")
        reports = get_reports(res)
        diagnostics = vscode_diagnostics(res.analyzer,
                                         reports,
                                         res.source,
                                         )
        check_basic_integration(diagnostics, reports)
        @test !isempty(diagnostics.items) && !isempty(reports)
        item = first(diagnostics.items)
        report = first(reports)
        check_inference_integration(item, report)
    end

    # no error case
    let res = @report_call sum(1:100)
        reports = get_reports(res)
        diagnostics = vscode_diagnostics(res.analyzer,
                                         reports,
                                         res.source,
                                         )
        check_basic_integration(diagnostics, reports)
        @test isempty(diagnostics.items) && isempty(reports)
    end

    # top-level integration (from top-level error)
    let res = mktemp() do path, io
            s = quote
                macro foo()
                    throw("foo")
                end
                @foo
            end |> string
            write(path, s)
            report_file2(path)
        end
        reports = get_reports(res)
        diagnostics = vscode_diagnostics(res.analyzer,
                                         reports,
                                         res.source,
                                         gen_postprocess(res.res.actual2virtual),
                                         )
        check_basic_integration(diagnostics, reports)
        @test !isempty(diagnostics.items) && !isempty(reports)
        item = first(diagnostics.items)
        report = first(reports)
        check_toplevel_integration(item, report)
    end

    # top-level integration (from inference error)
    let res = mktemp() do path, io
            s = quote
                undefvar
            end |> string
            write(path, s)
            report_file2(path)
        end
        reports = get_reports(res)
        diagnostics = vscode_diagnostics(res.analyzer,
                                         reports,
                                         res.source,
                                         gen_postprocess(res.res.actual2virtual),
                                         )
        check_basic_integration(diagnostics, reports)
        @test !isempty(diagnostics.items) && !isempty(reports)
        item = first(diagnostics.items)
        report = first(reports)
        check_inference_integration(item, report)
    end
end


text:
@testset "wigner" begin

end


text:
# using Revise
using POMDPStressTesting
using Distributions
using Parameters
using LinearAlgebra

@with_kw mutable struct Walk2DParams
    startx::Float64 = 0 # Starting x-position
    starty::Float64 = 0 # Starting y-position
    threshx::Float64 = 10 # +- x boundary threshold
    threshy::Float64 = 10 # +- y boundary threshold
    endtime::Int64 = 30 # Simulate end time
    corners::Matrix = [threshx  threshy;
                      -threshx  threshy;
                       threshx -threshy;
                      -threshx -threshy]
end

# Implement abstract GrayBox.Simulation
@with_kw mutable struct Walk2DSim <: GrayBox.Simulation
    params::Walk2DParams = Walk2DParams() # Parameters
    x::Float64 = 0 # Current x-position
    y::Float64 = 0 # Current x-position
    t::Int64 = 0 # Current time
    x_distribution::Distribution = Normal(0, 1) # x transition distribution
    y_distribution::Distribution = Normal(0, 1) # y transition distribution
end

# Override from GrayBox
GrayBox.environment(sim::Walk2DSim) = GrayBox.Environment(:x => sim.x_distribution, :y => sim.y_distribution)

# Override from GrayBox (NOTE: used with ASTSeedAction)
function GrayBox.transition!(sim::Walk2DSim)
    # We sample the environment and apply the transition
    environment::GrayBox.Environment = GrayBox.environment(sim) # Get the environment distributions
    sample::GrayBox.EnvironmentSample = rand(environment) # Sample from the environment
    sim.t += 1 # Keep track of time
    sim.x += sample[:x].value # Move agent using sampled x-value
    sim.y += sample[:y].value # Move agent using sampled y-value
    return logpdf(sample)::Real # Summation handled by `logpdf()`
end

# Override from GrayBox (NOTE: used with ASTSampleAction)
function GrayBox.transition!(sim::Walk2DSim, sample::GrayBox.EnvironmentSample)
    # The environment was sampled for us, and we just apply the transition
    sim.t += 1 # Keep track of time
    sim.x += sample[:x].value # Move agent using sampled value from input
    sim.y += sample[:y].value # Move agent using sampled value from input
    return logpdf(sample)::Real # Summation handled by `logpdf()`
end

# Override from BlackBox
function BlackBox.initialize!(sim::Walk2DSim)
    sim.t = 0
    sim.x = sim.params.startx
    sim.y = sim.params.starty
end

# Override from BlackBox
BlackBox.distance(sim::Walk2DSim) = minimum(mapslices(corner->norm([sim.x, sim.y] - corner), sim.params.corners, dims=2))

# Override from BlackBox
BlackBox.isevent(sim::Walk2DSim) = abs(sim.x) >= sim.params.threshx && abs(sim.y) >= sim.params.threshy

# Override from BlackBox
BlackBox.isterminal(sim::Walk2DSim) = BlackBox.isevent(sim) || sim.t >= sim.params.endtime

# Override from BlackBox (NOTE: used with ASTSeedAction)
function BlackBox.evaluate!(sim::Walk2DSim)
    logprob::Real = GrayBox.transition!(sim) # Step simulation
    d::Real = BlackBox.distance(sim) # Calculate miss distance
    event::Bool = BlackBox.isevent(sim) # Check event indication
    return (logprob::Real, d::Real, event::Bool)
end

# Override from BlackBox (NOTE: used with ASTSampleAction)
function BlackBox.evaluate!(sim::Walk2DSim, sample::GrayBox.EnvironmentSample)
    logprob::Real = GrayBox.transition!(sim, sample) # Step simulation given input sample
    d::Real = BlackBox.distance(sim) # Calculate miss distance
    event::Bool = BlackBox.isevent(sim) # Check event indication
    return (logprob::Real, d::Real, event::Bool)
end

function setup_ast(seed=AST.DEFAULT_SEED; solver=PPOSolver)
    # Create gray-box simulation object
    sim::GrayBox.Simulation = Walk2DSim()

    # AST MDP formulation object
    # NOTE: Use either {ASTSeedAction} or {ASTSampleAction}
    mdp::ASTMDP = ASTMDP{ASTSampleAction}(sim) # ASTSampleAction for use with DRL solvers (TRPO/PPO)
    mdp.params.debug = true # record metrics
    mdp.params.top_k = 10 # record top k best trajectories
    mdp.params.seed = seed # set RNG seed for determinism
    n_iterations = 1000 # number of algorithm iterations

    # Choose a solver
    if solver == RandomSearchSolver
        solver = RandomSearchSolver(n_iterations=n_iterations,
                                    episode_length=sim.params.endtime)
    elseif solver == MCTSPWSolver
        solver = MCTSPWSolver(n_iterations=n_iterations,
                              exploration_constant=1.0, # UCT exploration
                              k_action=1.0, # action widening
                              alpha_action=0.5, # action widening
                              depth=sim.params.endtime) # tree depth (i.e. episode length)
    elseif solver == CEMSolver
        solver = CEMSolver(n_iterations=n_iterations,
                           episode_length=sim.params.endtime)
    elseif solver == TRPOSolver
        solver = TRPOSolver(num_episodes=n_iterations,
                            episode_length=sim.params.endtime)
    elseif solver == PPOSolver
        solver = PPOSolver(num_episodes=n_iterations,
                           episode_length=sim.params.endtime)
    end

    # Get online planner (no work done, yet)
    planner = solve(solver, mdp)

    return planner
end

function run_ast(seed=AST.DEFAULT_SEED; kwargs...)
    planner = setup_ast(seed; kwargs...)

    action_trace::Vector{ASTAction} = search!(planner) # work done here.
    final_state::ASTState = playback(planner, action_trace, sim->(sim.x, sim.y))
    failure_rate::Float64 = print_metrics(planner)

    return planner, action_trace::Vector{ASTAction}, failure_rate::Float64
end

(planner, action_trace, failure_rate) = run_ast()

using PyPlot
function plot_trace(mdp, k=mdp.params.top_k)
    trace = playback(mdp, get_top_path(mdp, k), sim->[sim.x, sim.y]; return_trace=true)
    x = map(first, trace)
    y= map(last, trace)
    plot(x, y)
end

nothing # Suppress REPL


text:
@testset "System" begin
    mesh = buildSquare(0,1,0,1,2,2)

    sys = System{Float64}(mesh)

    @test sys.mesh.elements[1].id == 1

    dog = addVariable!(sys, "dog")
    cat = addVariable!(sys, "cat")

    @test dog.id == 1
    @test dog.name == "dog"

    @test cat.id == 2
    @test cat.name == "cat"

    initialize!(sys)

    @test sys.n_dofs == length(mesh.nodes) * length(sys.variables)

    @test mesh.nodes[1].dofs == [1,2]
    @test mesh.nodes[3].dofs == [9,10]

    # Test getting the dof indices
    dog_dofs = MOOSE.connectedDofIndices(mesh.elements[1], dog)
    @test dog_dofs == [1, 3, 5, 7]
    cat_dofs = MOOSE.connectedDofIndices(mesh.elements[1], cat)
    @test cat_dofs == [2, 4, 6, 8]

    dog_dofs = MOOSE.connectedDofIndices(mesh.elements[4], dog)
    @test dog_dofs == [5, 11, 17, 13]
    cat_dofs = MOOSE.connectedDofIndices(mesh.elements[4], cat)
    @test cat_dofs == [6, 12, 18, 14]

    # Test reinit
    MOOSE.reinit!(sys, mesh.elements[1], ones(18))

    @test_approx_eq dog.value [1,1,1,1]
    @test_approx_eq cat.value [1,1,1,1]
end

@testset "SystemFD" begin
    mesh = buildSquare(0,1,0,1,2,2)

    sys = System{Dual{8, Float64}}(mesh)

    @test sys.mesh.elements[1].id == 1

    dog = addVariable!(sys, "dog")
    cat = addVariable!(sys, "cat")

    @test dog.id == 1
    @test dog.name == "dog"

    @test cat.id == 2
    @test cat.name == "cat"

    initialize!(sys)

    @test sys.n_dofs == length(mesh.nodes) * length(sys.variables)

    @test mesh.nodes[1].dofs == [1,2]
    @test mesh.nodes[3].dofs == [9,10]

    # Test getting the dof indices
    dog_dofs = MOOSE.connectedDofIndices(mesh.elements[1], dog)
    @test dog_dofs == [1, 3, 5, 7]
    cat_dofs = MOOSE.connectedDofIndices(mesh.elements[1], cat)
    @test cat_dofs == [2, 4, 6, 8]

    dog_dofs = MOOSE.connectedDofIndices(mesh.elements[4], dog)
    @test dog_dofs == [5, 11, 17, 13]
    cat_dofs = MOOSE.connectedDofIndices(mesh.elements[4], cat)
    @test cat_dofs == [6, 12, 18, 14]

    # Test reinit
    MOOSE.reinit!(sys, mesh.elements[1], ones(18))

    @test_approx_eq dog.value [1,1,1,1]
    @test_approx_eq cat.value [1,1,1,1]

    @test partials(dog.value[1])[1] == dog.phi[1][1]
    @test partials(dog.value[1])[2] == dog.phi[1][2]
    @test partials(dog.value[1])[3] == dog.phi[1][3]
    @test partials(dog.value[1])[4] == dog.phi[1][4]

    @test partials(dog.value[1])[5] == 0.
    @test partials(dog.value[1])[6] == 0.
    @test partials(dog.value[1])[7] == 0.
    @test partials(dog.value[1])[8] == 0.

    @test partials(cat.value[1])[1] == 0.
    @test partials(cat.value[1])[2] == 0.
    @test partials(cat.value[1])[3] == 0.
    @test partials(cat.value[1])[4] == 0.

    @test partials(cat.value[1])[5] == cat.phi[1][1]
    @test partials(cat.value[1])[6] == cat.phi[1][2]
    @test partials(cat.value[1])[7] == cat.phi[1][3]
    @test partials(cat.value[1])[8] == cat.phi[1][4]
end


text:
export put_params!

"""
    P = put_params!(NL::NeuralNetLayer, Params::Array{Parameter,1})
 Inputs a initialized ResidualBlock and desired parameters values to load.
"""
function put_params!(RB::ResidualBlock, Params::Array{Any,1})
    RB.W1.data = Params[1].data
    RB.W2.data = Params[2].data
    RB.W3.data = Params[3].data
    RB.b1.data = Params[4].data
    RB.b2.data = Params[5].data
end

function put_params!(FB::FluxBlock, Params::Array{Any,1})
    Flux.loadparams!(FB, Params)
end

function put_params!(AL::AffineLayer, Params::Array{Any,1})
    AL.s.data = Params[1].data
    AL.b.data = Params[2].data
end

function put_params!(L::CouplingLayerIRIM, Params::Array{Any,1})
    put_params!(L.C, Params[1:3])
    put_params!(L.RB, Params[4:end])
end

function put_params!(HL::HyperbolicLayer, Params::Array{Any,1})
    HL.W.data = Params[1].data
    HL.b.data = Params[2].data
end

function put_params!(H::CouplingLayerHINT, Params::Array{Any,1})
    nlayers = length(H.CL)
    put_params!(H.CL[1], Params[1:5])
    if nlayers > 1
        for j = 2:nlayers
            put_params!(H.CL[j], Params[5*(j-1)+1:5*j])
        end
    end
    ~isnothing(H.C) && put_params!(H.C, Params[5*nlayers+1:end])
end

function put_params!(L::CouplingLayerGlow, Params::Array{Any,1})
    put_params!(L.C, Params[1:3])
    put_params!(L.RB, Params[4:end])
end

function put_params!(C::Conv1x1, Params::Array{Any,1})
    C.v1.data = Params[1].data
    C.v2.data = Params[2].data
    C.v3.data = Params[3].data
end

put_params!(L::CouplingLayerBasic, Params::Array{Any,1}) = put_params!(L.RB, Params)

function put_params!(AN::ActNorm, Params::Array{Any,1})
    AN.s.data = Params[1].data
    AN.b.data = Params[2].data
end

function put_params!(CH::ConditionalLayerHINT, Params::Array{Any,1})
    idx_s = 1
    counter = 0

    nlayers_x = length(CH.CL_X.CL)
    isnothing(CH.CL_X.C) ? counter += 5 * nlayers_x : counter += 5 * nlayers_x + 3
    put_params!(CH.CL_X, Params[idx_s:idx_s+counter-1])

    idx_s += counter
    counter = 0
    nlayers_y = length(CH.CL_Y.CL)
    isnothing(CH.CL_Y.C) ? counter += 5 * nlayers_y : counter += 5 * nlayers_y + 3
    put_params!(CH.CL_Y, Params[idx_s:idx_s+counter-1])

    idx_s += counter
    counter = 5
    put_params!(CH.CL_YX, Params[idx_s:idx_s+counter-1])

    idx_s += counter
    counter = 3
    ~isnothing(CH.C_X) && put_params!(CH.C_X, Params[idx_s:idx_s+counter-1])
    idx_s += counter
    ~isnothing(CH.C_Y) && put_params!(CH.C_Y, Params[idx_s:idx_s+counter-1])
    idx_s += counter
end

function put_params!(RB::ConditionalResidualBlock, Params::Array{Any,1})
    RB.W0.data = Params[1].data
    RB.W1.data = Params[2].data
    RB.W2.data = Params[3].data
    RB.W3.data = Params[4].data
    RB.b0.data = Params[5].data
    RB.b1.data = Params[6].data
    RB.b2.data = Params[7].data
end

function put_params!(CH::NetworkConditionalHINT, Params::Array{Any,1})
    depth = length(CH.CL)
    idx_s = 1
    counter = 0

    for j = 1:depth
        idx_s += counter
        counter = 2
        put_params!(CH.AN_X[j], Params[idx_s:idx_s+counter-1])

        idx_s += counter
        counter = 2
        put_params!(CH.AN_Y[j], Params[idx_s:idx_s+counter-1])

        idx_s += counter
        counter = 0
        nlayers_x = length(CH.CL[j].CL_X.CL)
        isnothing(CH.CL[j].CL_X.C) ? counter += 5 * nlayers_x : counter += 5 * nlayers_x + 3
        nlayers_y = length(CH.CL[j].CL_Y.CL)
        isnothing(CH.CL[j].CL_Y.C) ? counter += 5 * nlayers_y : counter += 5 * nlayers_y + 3
        counter += 5
        ~isnothing(CH.CL[j].C_X) && (counter += 3)
        ~isnothing(CH.CL[j].C_Y) && (counter += 3)

        put_params!(CH.CL[j], Params[idx_s:idx_s+counter-1])
    end
end

function put_params!(G::NetworkGlow, Params::Array{Any,1})
    L, K = size(G.AN)
    idx_s = 1
    counter = 0

    for i = 1:L
        for j = 1:K
            idx_s += counter
            counter = 2
            put_params!(G.AN[i, j], Params[idx_s:idx_s+counter-1])

            idx_s += counter
            counter = 8
            put_params!(G.CL[i, j], Params[idx_s:idx_s+counter-1])
        end
    end
end

function put_params!(H::NetworkHyperbolic, Params::Array{Any,1})
    depth = length(H.CL)
    put_params!(H.AL, Params[1:2])
    for j = 1:depth
        put_params!(H.HL[j], Params[2*j+1:2*(j+1)])
    end
end

function put_params!(UL::NetworkLoop, Params::Array{Any,1})
    maxiter = length(UL.L)

    if UL.L[1] isa CouplingLayerHINT
        idx_s = 1
        counter = 0
        nlayers_x = length(UL.L[1].CL)
        isnothing(UL.L[1].C) ? counter += 5 * nlayers_x : counter += 5 * nlayers_x + 3
        put_params!(UL.L[1], Params[idx_s:idx_s+counter-1])

        if maxiter > 1
            for j = 2:maxiter
                idx_s += counter
                counter = 0
                nlayers_x = length(UL.L[j].CL)
                isnothing(UL.L[j].C) ? counter += 5 * nlayers_x :
                counter += 5 * nlayers_x + 3
                put_params!(UL.L[j], Params[idx_s:idx_s+counter-1])
            end
        end

    elseif UL.L[1] isa CouplingLayerIRIM
        idx_s = 1
        counter = 8
        put_params!(UL.L[1], Params[idx_s:idx_s+counter-1])

        if maxiter > 1
            for j = 2:maxiter
                idx_s += counter
                put_params!(UL.L[j], Params[idx_s:idx_s+counter-1])
            end
        end
    end
end


text:
#
# 2D space, where all organisms are live. In reality this is
# just a peace of memory, where all organisms are located. It
# doesn't contain organisms codes, only rectangular with points.
# It's possible to run our application only in memory. In this
# case only this 2D world will be used (without visual
# presentation).
#
# Events:
#     dot Fires if new dot was paited.
#
# Usage:
#     import World
#     import Helper
#     ...
#     world = World.create()
#     World.setEnergy(world, Helper.Point(30, 40), UInt16(123456))
#     World.getEnergy(world, Helper.Point(30, 40)) # returns UInt16(123456))
#     ...
#     World.destroy(world)
#
# @author DeadbraiN
#
module World
  import Helper
  import Event
  import Dots
  import Creature

  export Plane

  export EVENT_DOT
  export EVENT_MOVE

  export create
  export destroy
  export setEnergy
  export getEnergy
  export moveEnergy
  export getFreePos
  export grabEnergy
  export getNearFreePos
  #
  # 2D plane for living
  #
  type Plane
    #
    # World width in pixels
    #
    width::Int
    #
    # World height in pixels
    #
    height::Int
    #
    # Array of pixels (RGB+8b)
    #
    data::Array{UInt16, 2}
    #
    # Observer for sending events
    #
    obs::Event.Observer
  end
  #
  # Event, which is fired if new dot was painted in a world
  #
  const EVENT_DOT  = "dot"
  #
  # Event, which is fired if some dot was moved from one
  # location to another.
  #
  const EVENT_MOVE = "move"
  #
  # Creates the world
  # @param width World width
  # @param height World height
  # @return {Plane} filled by zero values Plane
  #
  function create(width::Int, height::Int)
    Plane(width, height, fill(UInt16(Dots.INDEX_EMPTY), height, width), Event.create())
  end
  #
  # Opposite to create. Destroys the World
  # @param plane World instance
  #
  function destroy(plane::Plane)
    plane.data = Array{UInt16, 2}()
    Event.clear(plane.obs)
    plane.width  = 0
    plane.height = 0
  end
  #
  # Adds energy point by specified coordinates
  # @param plane World's plane
  # @param pos Position of the energy point
  # @param energy Amount of energy to add
  # @return {Bool} It's possible, that energy position is not empty.
  # false in this case. true - if dot was set.
  #
  function setEnergy(plane::Plane, pos::Helper.Point, energy::UInt16)
    if pos.x < 1 || pos.x > plane.width || pos.y < 1 || pos.y > plane.height return false end
    plane.data[pos.y, pos.x] = energy
    Event.fire(plane.obs, EVENT_DOT, pos, energy)

    true
  end
  #
  # Moves energy point from one location to another. Old coordinates
  # will be filled by empty(zero) color. Moving is only works with
  # one pixel (1px distance).
  # @param plane World's plane
  # @param pos Position of the energy point
  # @param energy Amount of energy to add
  # @return {Bool} It's possible, that energy position is not empty.
  # false in this case. true - if dot was set.
  #
  function moveEnergy(plane::Plane, oldPos::Helper.Point, newPos::Helper.Point, org::Creature.Organism, outOfBorder::Bool = false)
    if newPos.x < 1 || newPos.x > plane.width || newPos.y < 1 || newPos.y > plane.height return false end
    local dir::Int
    local energy::UInt16 = UInt16(min(org.energy, Creature.ORGANISM_MAX_ENERGY))

    if     newPos.x > oldPos.x dir = Dots.DIRECTION_RIGHT
    elseif newPos.x < oldPos.x dir = Dots.DIRECTION_LEFT
    elseif newPos.y > oldPos.y dir = Dots.DIRECTION_DOWN
    elseif newPos.y < oldPos.y dir = Dots.DIRECTION_UP
    else dir = Dots.DIRECTION_NO
    end
    plane.data[oldPos.y, oldPos.x] = Dots.INDEX_EMPTY
    plane.data[newPos.y, newPos.x] = energy
    #
    # infoBits: 7th bit - organism goes out of border
    #
    Event.fire(plane.obs, EVENT_MOVE, newPos, dir, org.color, org.id, outOfBorder)

    true
  end
  #
  # Returns amount of energy in a point with
  # specified coordinates
  # @param plane World's plane, where to check
  # @param pos Position to check
  # @return {UInt16} Amount of energy
  #
  function getEnergy(plane::Plane, pos::Helper.Point)
    if pos.x < 1 || pos.x > plane.width || pos.y < 1 || pos.y > plane.height return UInt16(0) end
    plane.data[pos.y, pos.x]
  end
  #
  # Returns free position in the world. Free means, the point
  # with zero energy. It randomly finds the point in a world
  # and checks it's energy. If this point already has an energy
  # it finds another point randomly.
  # @return {Helper.Point|Bool}
  #
  function getFreePos(plane::Plane)
    local pos::Helper.Point = Helper.Point(ceil(Int, plane.width / 2), ceil(Int, plane.height / 2))
    local width::Int = plane.width
    local height::Int = plane.height
    local i::Int = width * height < 1000 ? 100 : div(width * height, 10)

    while World.getEnergy(plane, pos) > UInt16(0) && i > 0
      pos.x = Helper.fastRand(width)
      pos.y = Helper.fastRand(height)
      i -= 1
    end
    i > 0 ? pos : false
  end
  #
  # Grabs energy in specified point in a world. If this point
  # has no energy, 0 value will be returned. If point has less
  # energy, then "amount", then all point energy will be returned.
  # @param plane Plane
  # @param pos Position to grab
  # @param amount Amount of energy we want to grab
  # @return {UInt} amount of grabbed energy
  #
  function grabEnergy(plane::Plane, pos::Helper.Point, amount::UInt16)
    local energy::UInt16 = min(getEnergy(plane, pos), amount)

    if energy > UInt16(0)
      Event.fire(plane.obs, EVENT_DOT, pos, (plane.data[pos.y, pos.x] -= energy))
    end

    energy
  end
  #
  # Returns position near specified. It checks up, down, left and right
  # sides of "pos". If every sides are full, then returns false
  # @param plane Plane
  # @param pos Start position
  # @return {Helper.Point|Bool}
  # TODO: this method should return Helper.Point() in any case. Point(0,0)
  # TODO: should mean no free place
  function getNearFreePos(plane::Plane, pos::Helper.Point)
    local i::Int
    local j::Int
    local x::Int = pos.x
    local y::Int = pos.y
    local positions::Array{Int, 1} = [
      x + 1, y,     # right
      x + 1, y + 1, # right down
      x    , y + 1, # down
      x - 1, y + 1, # down left
      x - 1, y,     # left
      x - 1, y - 1, # left up
      x    , y - 1, # up
      x + 1, y - 1  # up right
    ]
    #
    # We have to have a copy of position to prevent of organism moving
    #
    pos = Helper.Point(0, 0)

    j = 1
    for i = 1:8
      pos.x = positions[j]
      pos.y = positions[j + 1]
      if plane.width >= pos.x && plane.height >= pos.y && pos.x > 0 && pos.y > 0 && World.getEnergy(plane, pos) === UInt16(0)
        return pos
      end
      j += 2
    end

    false
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
using Documenter, DocumenterCitations, FastDEQ

bib = CitationBibliography(joinpath(@__DIR__, "ref.bib"), sorting = :nyt)

makedocs(
    bib,
    sitename = "Fast Deep Equilibrium Networks",
    authors = "Avik Pal et al.",
    clean = true,
    doctest = false,
    modules = [FastDEQ],
    format = Documenter.HTML(analytics = "UA-90474609-3",
                             #  assets = ["assets/favicon.ico"],
                             canonical="https://fastdeq.sciml.ai/stable/"),
    pages = [
        "FastDEQ: Fast Deep Equilibrium Networks" => "index.md",
        "API" => [
            "Dynamical Systems" => "api/solvers.md",
            "Non Linear Solvers" => "api/nlsolve.md",
            "General Purpose Layers" => "api/layers.md",
            "DEQ Layers" => "api/deqs.md",
            "Miscellaneous" => "api/misc.md",
        ],
        "References" => "references.md",
    ]
)

deploydocs(
   repo = "github.com/SciML/FastDEQ.jl.git";
   push_preview = true
)


text:
module Mixers

using MacroTools

export @mix, @premix, @pour

"""
Generates simple macros for inserting any lines anywhere.
"""
macro pour(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

"""
Generates a mixin macro that preserve parametric types.
Fields and parametric types are appended to the struct.
Identical parametric types are merged.
"""
macro mix(ex)
    defmacro(ex, false)
end

"""
Just like @mix but generated macro insert fields and types
at the *start* of the definition.
"""
macro premix(ex)
    defmacro(ex, true)
end

function defmacro(ex, prepend)
    # get chained macros
    macros = chained_macros(ex)

    # get name and parametric types
    @capture(firsthead(ex, :struct).args[2], mixname_{mixtypes__} | mixname_ )
    # get fields
    mixfields = firsthead(ex, :block).args

    # deal with empty values
    if mixtypes == nothing mixtypes = [] end
    if mixfields == nothing mixfields = [] end

    @esc macros mixname mixtypes mixfields
    return quote
        macro $mixname(ex)
            macros = $macros
            mixtypes = $mixtypes
            mixfields = $mixfields
            prepend = $prepend
            return mix(ex, macros, mixtypes, mixfields, prepend)
        end
    end
end

function mix(ex, macros, mixtypes, mixfields, prepend)
    # merge type parameters
    firsthead(ex, :curly) do x
        x.args = vcat(x.args[1], mergetypes(x.args[2:end], mixtypes, prepend))
    end
    # merge fields
    firsthead(ex, :block) do x
        x.args = vcat(x.args[1], mergefields(x.args[2:end], mixfields, prepend))
    end

    localmacros = chained_macros(ex)

    # get struct without macros
    firsthead(x -> ex = x, ex, :struct)

    # wrap local and mixed macros around the struct
    for mac in reverse(union(localmacros, macros))
        ex = Expr(:macrocall, mac, LineNumberNode(79, Symbol("Mixers.jl")), ex)
    end
    esc(ex)
end

chained_macros(ex) = chained_macros!(Symbol[], ex)
chained_macros!(macros, ex) = macros
chained_macros!(macros, ex::Expr) = begin
    if ex.head == :macrocall
        push!(macros, ex.args[1])
        length(ex.args) > 2 && chained_macros!(macros, ex.args[3])
    end
    macros
end

firsthead(ex, match) = firsthead(x->x, ex, match)
firsthead(f, ex::Expr, match) =
    if ex.head == match
        return f(ex)
    else
        for arg in ex.args
            x = firsthead(f, arg, match)
            x == nothing || return x
        end
        return nothing
    end
firsthead(f, ex, match) = nothing

mergetypes(t1, t2, prepend) = prepend ? union(t2, t1) : union(t1, t2)
mergefields(f1, f2, prepend) = prepend ? vcat(f2, f1) : vcat(f1, f2)

end # module


text:
@doc raw"""
    LatentFactorMOKernel(g, e::MOKernel, A::AbstractMatrix)

The kernel associated with the Semiparametric Latent Factor Model, introduced by
Seeger, Teh and Jordan (2005).

``k((x, p_x), (y, p_y)) = \Sum^{Q}_{q=1} A_{p_xq}g_q(x, y)A_{p_yq} + e((x, p_x), (y, p_y))``

# Arguments
- `g`: a collection of kernels, one for each latent process
- `e`: a [`MOKernel`](@ref) - multi-output kernel
- `A::AbstractMatrix`: a matrix of weights for the kernels of size `(out_dim, length(g))`

# Reference:
- [Seeger, Teh, and Jordan (2005)](https://infoscience.epfl.ch/record/161465/files/slfm-long.pdf)

"""
struct LatentFactorMOKernel{Tg, Te <: MOKernel, TA <: AbstractMatrix} <: MOKernel
    g::Tg
    e::Te
    A::TA
    function LatentFactorMOKernel(g, e::MOKernel, A::AbstractMatrix)
        all(gi isa Kernel for gi in g) || error("`g` should be an collection of kernels")
        length(g) == size(A, 2) ||
            error("Size of `A` not compatible with the given array of kernels `g`")
        return new{typeof(g), typeof(e), typeof(A)}(g, e, A)
    end
end

function (κ::LatentFactorMOKernel)((x, px)::Tuple{Any, Int}, (y, py)::Tuple{Any, Int})
    cov_f = sum(κ.A[px, q] * κ.g[q](x, y) * κ.A[py, q] for q in 1:length(κ.g))
    return cov_f + κ.e((x, px), (y, py))
end

function kernelmatrix(k::LatentFactorMOKernel, x::MOInput, y::MOInput)
    x.out_dim == y.out_dim || error("`x` and `y` should have the same output dimension")
    x.out_dim == size(k.A, 1) ||
        error("Kernel not compatible with the given multi-output inputs")

    # Weights matrix ((out_dim x out_dim) x length(k.g))
    W = [col * col' for col in eachcol(k.A)]

    # Latent kernel matrix ((N x N) x length(k.g))
    H = [gi.(x.x, permutedims(y.x)) for gi in k.g]

    # Weighted latent kernel matrix ((N*out_dim) x (N*out_dim))
    W_H = sum(kron(Wi, Hi) for (Wi, Hi) in zip(W, H))

    return W_H .+ kernelmatrix(k.e, x, y)
end

function Base.show(io::IO, k::LatentFactorMOKernel)
    print(io, "Semi-parametric Latent Factor Multi-Output Kernel")
end

function Base.show(io::IO, ::MIME"text/plain", k::LatentFactorMOKernel)
    print(io, "Semi-parametric Latent Factor Multi-Output Kernel\n\tgᵢ: ")
    join(io, k.g, "\n\t\t")
    print(io, "\n\teᵢ: ")
    join(io, k.e, "\n\t\t")
end


text:
module Err
export PatternUnsolvedException, InternalException, SyntaxError, UnknownExtension, @syntax_err

struct PatternUnsolvedException <: Exception
    msg :: String
    PatternUnsolvedException(arg) =
        if isa(arg, String)
            new(arg)
        else
            new("Non-exhaustive pattern found for target `$(string(arg))`.")
        end
end

struct InternalException <: Exception
    msg :: String
end

struct SyntaxError <: Exception
    msg :: String
end

struct UnknownExtension <: Exception
    ext :: Union{String, Symbol}
end

macro syntax_err(msg)
    esc(:($throw($SyntaxError($msg))))
end

end


text:
using LCPSim
using Base.Test
using RigidBodyDynamics
using RigidBodyDynamics: Bounds
using MeshCat
using MeshCatMechanisms
using Cbc: CbcSolver
using Gurobi: GurobiSolver
using Rotations: RotY
using JuMP: setvalue, getvalue

urdf = joinpath(@__DIR__, "urdf", "box.urdf")

function box_with_planar_base()
    mechanism = parse_urdf(Float64, urdf)
    core = findbody(mechanism, "core")
    fixed_joint = joint_to_parent(core, mechanism)
    floating_base = Joint(fixed_joint.name, frame_before(fixed_joint), frame_after(fixed_joint),
                          Planar([1., 0, 0], [0., 0, 1.]),
                          position_bounds=[Bounds(-5., 5), Bounds(0., 3), Bounds(-2π, 2π)],
                          velocity_bounds=[Bounds(-10., 10), Bounds(-10., 10), Bounds(-10., 10)],
                          effort_bounds=[Bounds(0., 0) for i in 1:3])
    replace_joint!(mechanism, fixed_joint, floating_base)
    x0 = MechanismState{Float64}(mechanism)
    set_configuration!(x0, floating_base, [-1, 0.3, -0.3])
    configuration_derivative_to_velocity!(velocity(x0, floating_base), floating_base, configuration(x0, floating_base), [1., 0, 0])
    mechanism, x0
end

function environment_with_floor(mechanism)
    world = root_body(mechanism)
    floor = planar_obstacle(default_frame(world), [0, 0, 1.], [0, 0, 0.], 0.5)

    core = findbody(mechanism, "core")
    env = Environment([
        (core, pt, floor) for pt in [
                    Point3D(default_frame(core), SVector(0.1, 0, 0.2)),
                    Point3D(default_frame(core), SVector(-0.1, 0, 0.2)),
                    Point3D(default_frame(core), SVector(0.1, 0, -0.2)),
                    Point3D(default_frame(core), SVector(-0.1, 0, -0.2)),
                     ]
            ])
    env
end

@testset "simulation warmstarts" begin
    mech1, x1 = box_with_planar_base()
    q0 = copy(configuration(x1))
    v0 = copy(velocity(x1))
    env1 = environment_with_floor(mech1)

    controller = LCPSim.passive_controller()
    Δt = 0.01
    N = 50
    eps = 2e-3

    mvis = MechanismVisualizer(mech1, URDFVisuals(urdf))
    if !haskey(ENV, "CI")
        open(mvis.visualizer)
        wait(mvis.visualizer)
    end

    @testset "linear dynamics" begin
        # First test the linear dynamics
        set_configuration!(x1, q0)
        set_velocity!(x1, v0)
        results1 = LCPSim.simulate(x1, controller, env1, Δt, N, GurobiSolver(Gurobi.Env(), OutputFlag=0); relinearize=false)

        for r in results1
            set_configuration!(mvis, configuration(r.state))
            sleep(Δt)
        end

        # Construct an optimization linearized around the initial state
        set_configuration!(x1, q0)
        set_velocity!(x1, v0)
        model2, results2 = LCPSim.optimize(x1, env1, Δt, N)
        setvalue.(results2, results1)
        ConditionalJuMP.warmstart!(model2, false)

        # Test that all the optimization constraints are satisfied by the
        # simulation result
        for c in model2.linconstr
            @test c.lb - eps <= getvalue(c.terms) <= c.ub + eps
        end
    end

    @testset "nonlinear dynamics" begin
        # Now simulate the nonlinear dynamics
        set_configuration!(x1, q0)
        set_velocity!(x1, v0)
        results3 = LCPSim.simulate(x1, controller, env1, Δt, N, GurobiSolver(Gurobi.Env(), OutputFlag=0); relinearize=true)

        for r in results3
            set_configuration!(mvis, configuration(r.state))
            sleep(Δt)
        end

        # Construct an optimization linearzed around the trajectory
        set_configuration!(x1, q0)
        set_velocity!(x1, v0)
        model4, results4 = LCPSim.optimize(x1, env1, Δt, results3)
        setvalue.(results4, results3)
        ConditionalJuMP.warmstart!(model4, false)

        # Test that all the optimization constraints are satisfied by the
        # simulation result
        for c in model4.linconstr
            @test c.lb - eps <= getvalue(c.terms) <= c.ub + eps
            if !(c.lb - eps <= getvalue(c.terms) <= c.ub + eps)
                @show c
            end
        end
    end
end


text:
import Base: iterate

mutable struct BatchIter
    N::Int
    batch_size::Int
end

function iterate(bi::BatchIter)

    l_idx = 1
    r_idx = min(bi.N, l_idx + bi.batch_size - 1)

    if r_idx < l_idx
        return nothing
    end

    return l_idx:r_idx, r_idx
end

function iterate(bi::BatchIter, prev_r_idx)

    l_idx = prev_r_idx + 1

    r_idx = min(bi.N, l_idx + bi.batch_size - 1)

    if r_idx < l_idx
        return nothing
    end

    return l_idx:r_idx, r_idx

end

max_stars_repo_name: dpmerrell/BatchMatFac

max_issues_repo_name: dpmerrell/BatchMatFac

max_forks_repo_name: dpmerrell/BatchMatFac


text:
module ESDLHelpers

using ESDL
using YAXArrays
using Statistics
using RecurrenceAnalysis
using DataStructures
using StatsBase:skewness, kurtosis, mad
using EmpiricalModeDecomposition; const EMD=EmpiricalModeDecomposition
using LombScargle

export timestats, lombscargle, tslength

include("tsanalysis.jl")
include("shapesampling.jl")

include("gdalcubes.jl")

end # module


text:
notes_pattern = r"((?:[^:]+: \d+-\d+ or \d+-\d+\n)+)\nyour ticket:\n((?:\d+,?)+)\n\nnearby tickets:\n((?:(?:\d+,?)+\n?)+)"
fields, my_ticket, tickets = open("ticket_notes.txt") do f
    fields, my_ticket, tickets = match(notes_pattern, read(f, String)).captures
    my_ticket = parse.(Int, split(my_ticket, ","))
    tickets = map(ticket -> parse.(Int, split(ticket, ",")), split(tickets, "\n"))
    fields = map(r -> (r.captures[1], [parse(Int, r.captures[2]):parse(Int, r.captures[3]); parse(Int, r.captures[4]):parse(Int, r.captures[5])]), match.(r"([^:]+): (\d+)-(\d+) or (\d+)-(\d+)", split(strip(fields), "\n")))
    fields, my_ticket, tickets
end

function valid_ticket(ticket, ranges)::Bool
    for val in ticket
        if !any(r -> val in r, ranges)
            return false
        end
    end
    return true
end

ranges = Iterators.flatten(map(x -> (x[2]), fields))
invalid = Vector{Int64}()
for ticket in tickets
    for val in ticket
        if !any(r -> val in r, ranges)
            # println("$(ticket) is invalid ticket ($(val))")
            push!(invalid, val)
        end
    end
end

println(sum(invalid))

candidates = map(r -> map(f -> f[1], r), map(r -> filter(f -> all(x -> x in f[2], r), fields), eachrow(hcat(filter(t -> valid_ticket(t, ranges), tickets)...))))

while any(x -> length(x) > 1, candidates)
    solved = Iterators.flatten(filter(x -> length(x) == 1, candidates))
    for i in 1:length(candidates)
        if length(candidates[i]) == 1
            continue
        end
        candidates[i] = filter(c -> c ∉ solved, candidates[i])
    end
end

println(prod([my_ticket[i] for i in findall(x -> startswith(x, "departure"), collect(Iterators.flatten(candidates)))]))


text:
open("out-elemental.txt", "w") do outf
  write(outf, "# nbprocs nbthreads distributed_dimension mintime (s) maxtime (s)\n")
  for fname in ARGS
    open(fname) do f
      nthreads = 0
      ntasks = 0
      timing = 0.0
      readnext = false
      for ln in eachline(f)
        if readnext
          timing = parse(Float64, split(ln)[2])
          break
        end
        if startswith(ln, "OPENBLAS_NUM_THREADS")
          nthreads = parse(Int,split(ln)[2])
        end
        if startswith(ln, "TASKS")
          ntasks = parse(Int,split(ln)[2])
        end
        if contains(ln, "Field=double")
          readnext = true
        end
      end
      write(outf, "$ntasks $nthreads 1 $timing $timing\n")
    end
  end
end


text:
function adaptive_simpsons_inner(f::Function, a::Real, b::Real,
                                 epsilon::Real, S::Real,
                                 fa::Real, fb::Real, fc::Real, bottom::Int)
    c = (a + b) / 2
    h = b - a
    d = (a + c) / 2
    g = (c + b) / 2
    fd = f(d)
    fe = f(g)
    Sleft = (h / 12) * (fa + 4 * fd + fc)
    Sright = (h / 12) * (fc + 4 * fe + fb)
    S2 = Sleft + Sright
    if bottom <= 0 || abs(S2 - S) <= 15 * epsilon
        return S2 + (S2 - S) / 15
    end
    return adaptive_simpsons_inner(f, a, c, epsilon / 2, Sleft,  fa, fc, fd, bottom - 1) +
           adaptive_simpsons_inner(f, c, b, epsilon / 2, Sright, fc, fb, fe, bottom - 1)
end

function adaptive_simpsons_outer(f::Function, a::Real, b::Real, accuracy::Real, max_iterations::Int)
    c = (a + b) / 2
    h = b - a
    fa = f(a)
    fb = f(b)
    fc = f(c)
    S = (h / 6) * (fa + 4 * fc + fb)
  return adaptive_simpsons_inner(f, a, b, accuracy, S, fa, fb, fc, max_iterations)
end
adaptive_simpsons_outer(f::Function, a::Real, b::Real) = adaptive_simpsons_outer(f, a, b, 10e-10, 50)

function monte_carlo(f::Function, a::Real, b::Real, iterations::Int)
    estimate = 0.0
    width = (b - a)
    for i in 1:iterations
        x = width * rand() + a
        estimate += f(x) * width
    end
    return estimate / iterations
end

function integrate(f::Function, a::Real, b::Real, method::Symbol)
    if method == :simpsons
        adaptive_simpsons_outer(f, a, b)
    elseif method == :monte_carlo
        monte_carlo(f, a, b, 10_000)
    else
        error("Unknown method of integration: $(method)")
    end
end
integrate(f::Function, a::Real, b::Real) = integrate(f, a, b, :simpsons)


text:
module SimplexPoly

using ComputedFieldTypes
using DifferentialForms
using LinearAlgebra
using SparseArrays
using StaticArrays

################################################################################

export Term
struct Term{D,T}
    powers::SVector{D,Int}
    coeff::T
    function Term{D,T}(powers::SVector{D,Int}, coeff::T) where {D,T}
        D::Int
        @assert D >= 0
        @assert all(>=(0), powers)
        return new{D,T}(powers, coeff)
    end
end

function Base.show(io::IO, ::MIME"text/plain", x::Term{D,T}) where {D,T}
    skiptype = get(io, :typeinfo, Any) <: Term{D,T}
    if !skiptype
        print(io, "Term(")
    end
    print(io, x.coeff, " * ", x.powers)
    if !skiptype
        print(io, ")")
    end
    return nothing
end

function Base.:(==)(x::Term{D}, y::Term{D}) where {D}
    return x.powers == y.powers && x.coeff == y.coeff
end

function Base.isequal(x::Term, y::Term)
    return isequal((x.powers, x.coeff), (y.powers, y.coeff))
end
Base.isless(x::Term, y::Term) = isless((x.powers, x.coeff), (y.powers, y.coeff))
Base.hash(x::Term, h::UInt) = hash((x.powers, x.coeff), hash(0x11325cd9, h))
export compare
function compare(x::Term, y::Term)
    # Note reversed signs; we want (1,0) < (0,1)
    x.powers > y.powers && return -1
    x.powers < y.powers && return +1
    return 0
end

function Base.map(f, x::Term{D,T}, ys::Term{D,T}...) where {D,T}
    @assert all(x.powers == y.powers for y in ys)
    return Term{D,T}(x.powers, map(f, x.coeff, (y.coeff for y in ys)...))
end

function Base.zero(::Type{Term{D,T}}) where {D,T}
    return Term{D,T}(zero(SVector{D,Int}), zero(T))
end
Base.zero(::Term{D,T}) where {D,T} = zero(Term{D,T})
Base.iszero(x::Term{D,T}) where {D,T} = iszero(x.coeff)

function Base.one(::Type{Term{D,T}}) where {D,T}
    return Term{D,T}(zero(SVector{D,Int}), one(T))
end
Base.one(::Term{D,T}) where {D,T} = one(Term{D,T})
Base.isone(x::Term{D,T}) where {D,T} = iszero(x.powers) && isone(x.coeff)

Base.:+(x::Term{D,T}) where {D,T} = map(+, x)
Base.:-(x::Term{D,T}) where {D,T} = map(-, x)

Base.:+(x::Term{D,T}, y::Term{D,T}) where {D,T} = map(+, x, y)
Base.:-(x::Term{D,T}, y::Term{D,T}) where {D,T} = map(-, x, y)
Base.:*(a, x::Term{D,T}) where {D,T} = map(c -> a * c, x)
Base.:*(x::Term{D,T}, a) where {D,T} = map(c -> c * a, x)
Base.:\(a, x::Term{D,T}) where {D,T} = map(c -> a \ c, x)
Base.:/(x::Term{D,T}, a) where {D,T} = map(c -> c / a, x)
Base.div(x::Term{D,T}, a) where {D,T} = map(c -> div(c, a), x)
Base.mod(x::Term{D,T}, a) where {D,T} = map(c -> mod(c, a), x)

function Base.:*(x::Term{D,T}, y::Term{D,T}) where {D,T}
    return Term{D,T}(x.powers + y.powers, x.coeff * y.coeff)
end
function Base.:^(x::Term{D,T}, n::Integer) where {D,T}
    @assert n >= 0
    r = one(x)
    while n > 0
        if n % 2 != 0
            r *= x
        end
        x *= x
        n ÷= 2
    end
    return r
end

export deriv
function deriv(term::Term{D,T}, dir::Int) where {D,T}
    D::Int
    @assert D >= 0
    @assert 1 <= dir <= D
    p = term.powers[dir]
    p == 0 && return zero(Term{D,T})
    return Term{D,T}(Base.setindex(term.powers, p - 1, dir), p * term.coeff)
end

"""
    koszul

The Koszul operator `κ` (kappa)

See e.g.: Douglas Arnold, Richard Falk, Ragnar Winther, "Finite
element exterior calculus, homological techniques, and applications",
Acta Numerica 15, 1-155 (2006), DOI:10.1017/S0962492906210018.
"""
function koszul end

export koszul
function koszul(term::Term{D,T}, dir::Int) where {D,T}
    D::Int
    @assert D >= 0
    @assert 1 <= dir <= D
    p = term.powers[dir]
    return Term{D,T}(Base.setindex(term.powers, p + 1, dir), term.coeff)
end

################################################################################

export Poly
struct Poly{D,T} <: Number
    terms::Vector{Term{D,T}}
    function Poly{D,T}(terms::Vector{Term{D,T}}) where {D,T}
        D::Int
        @assert D >= 0
        terms = combine(terms)
        return new{D,T}(terms)
    end
end

function combine(terms::Vector{Term{D,T}}) where {D,T}
    terms = sort(terms; by=(t -> t.powers), rev=true)
    isempty(terms) && return terms
    newterms = Term{D,T}[]
    i = 0
    term = terms[i += 1]
    while i < length(terms)
        term2 = terms[i += 1]
        if term2.powers == term.powers
            term += term2
        else
            !iszero(term) && push!(newterms, term)
            term = term2
        end
    end
    !iszero(term) && push!(newterms, term)
    return newterms
end

function Base.show(io::IO, mime::MIME"text/plain", p::Poly{D,T}) where {D,T}
    skiptype = get(io, :typeinfo, Any) <: Poly{D,T}
    if !skiptype
        print(io, "Poly")
    end
    # print(io, "[")
    print(io, "(")
    for (i, term) in enumerate(p.terms)
        # i > 1 && print(io, ", ")
        i > 1 && print(io, " + ")
        show(IOContext(io, :compact => true, :typeinfo => Term{D,T}), mime,
             term)
    end
    # print(io, "]")
    print(io, ")")
    return nothing
end

Base.:(==)(p::Poly{D}, q::Poly{D}) where {D} = p.terms == q.terms

Base.isequal(p::Poly, q::Poly) = isequal(p.terms, q.terms)
Base.isless(p::Poly, q::Poly) = isless(p.terms, q.terms)
Base.hash(p::Poly, h::UInt) = hash(p.terms, hash(0x0b9503a7, h))
export compare
function compare(p::Poly, q::Poly)
    for i in 1:min(length(p.terms), length(q.terms))
        c = compare(p.terms[i], q.terms[i])
        c != 0 && return c
    end
    # Polynomials with more terms come first since a non-zero term
    # comes before the zero term
    length(p.terms) > length(q.terms) && return -1
    length(p.terms) < length(q.terms) && return +1
    return 0
end

Base.map(f, p::Poly{D,T}) where {D,T} = Poly{D,T}(map(t -> map(f, t), p.terms))

function Base.map(f, p::Poly{D,T}, q::Poly{D,T}) where {D,T}
    terms = Term{D,T}[]
    i = j = 1
    ni = length(p.terms)
    nj = length(q.terms)
    while i <= ni || j <= nj
        usei = usej = false
        if i <= ni && j <= nj
            c = cmp(p.terms[i].powers, q.terms[j].powers)
            usei = c <= 0
            usej = c >= 0
        else
            usei = i <= ni
            usej = !usei
        end
        pi = usei ? p.terms[i] : map(zero, q.terms[j])
        qj = usej ? q.terms[j] : map(zero, p.terms[i])
        push!(terms, map(f, pi, qj))
        i += usei
        j += usej
    end
    return Poly{D,T}(terms)
end

function Base.zero(::Type{Poly{D,T}}) where {D,T}
    D::Int
    @assert D >= 0
    return Poly{D,T}(Vector{Term{D,T}}())
end
Base.zero(::Poly{D,T}) where {D,T} = zero(Poly{D,T})
Base.iszero(x::Poly{D,T}) where {D,T} = isempty(x.terms)

function Base.one(::Type{Poly{D,T}}) where {D,T}
    return Poly{D,T}([one(Term{D,T})])
end
Base.one(::Poly{D,T}) where {D,T} = one(Poly{D,T})
Base.isone(x::Poly{D,T}) where {D,T} = length(x.terms) == 1 && isone(x.terms[1])

export unit
function Forms.unit(::Type{Poly{D,T}}, dir::Int, coeff=one(T)) where {D,T}
    D::Int
    @assert D >= 0
    @assert 1 <= dir <= D
    term = Term{D,T}(SVector{D,T}(d == dir for d in 1:D), coeff)
    return Poly{D,T}([term])
end

Base.:+(p::Poly{D,T}) where {D,T} = map(+, p)
Base.:-(p::Poly{D,T}) where {D,T} = map(-, p)

Base.:+(p::Poly{D,T}, q::Poly{D,T}) where {D,T} = map(+, p, q)
Base.:-(p::Poly{D,T}, q::Poly{D,T}) where {D,T} = map(-, p, q)
Base.:*(a::Number, p::Poly{D,T}) where {D,T} = map(t -> a * t, p)
Base.:*(p::Poly{D,T}, a::Number) where {D,T} = map(t -> t * a, p)
Base.:\(a::Number, p::Poly{D,T}) where {D,T} = map(t -> a \ t, p)
Base.:/(p::Poly{D,T}, a::Number) where {D,T} = map(t -> t / a, p)
Base.div(p::Poly{D,T}, a::Number) where {D,T} = map(t -> div(t, a), p)
Base.mod(p::Poly{D,T}, a::Number) where {D,T} = map(t -> mod(t, a), p)

function Base.:*(p::Poly{D,T}, q::Poly{D,T}) where {D,T}
    return Poly{D,T}(Term{D,T}[t * u for t in p.terms for u in q.terms])
end
function Base.:^(x::Poly{D,T}, n::Integer) where {D,T}
    @assert n >= 0
    r = one(x)
    while n > 0
        if n % 2 != 0
            r *= x
        end
        x *= x
        n ÷= 2
    end
    return r
end

function deriv(poly::Poly{D,T}, dir::Int) where {D,T}
    return Poly{D,T}(map(t -> deriv(t, dir), poly.terms))
end

function koszul(poly::Poly{D,T}, dir::Int) where {D,T}
    return Poly{D,T}(map(t -> koszul(t, dir), poly.terms))
end

################################################################################

# Note: PolySpace is unused

export PolySpace
struct PolySpace{D,T}
    polys::Vector{Poly{D,T}}
    function PolySpace{D,T}(polys::Vector{Poly{D,T}}) where {D,T}
        D::Int
        @assert D >= 0
        polys = combine(polys)
        return new{D,T}(polys)
    end
end

function combine(polys::Vector{Poly{D,T}}) where {D,T}
    # TODO: This is wrong -- need to combine polynomials that are
    # multiples of each other
    polys = sort(polys; rev=true)
    unique!(polys)
    filter!(p -> !iszero(p), polys)
    return polys
end

PolySpace{D,T}() where {D,T} = PolySpace{D,T}(Poly{D,T}[])

function Base.show(io::IO, mime::MIME"text/plain",
                   ps::PolySpace{D,T}) where {D,T}
    skiptype = get(io, :typeinfo, Any) <: PolySpace{D,T}
    if !skiptype
        print(io, "PolySpace")
    end
    print(io, "[")
    for (i, poly) in enumerate(ps.polys)
        i > 1 && print(io, ", ")
        show(IOContext(io, :compact => true, :typeinfo => Poly{D,T}), mime,
             poly)
    end
    print(io, "]")
    return nothing
end

function Base.:(==)(ps::PolySpace{D}, qs::PolySpace{D}) where {D}
    return ps.polys == qs.polys
end

Base.isequal(ps::PolySpace, qs::PolySpace) = isequal(ps.polys, qs.polys)
Base.isless(ps::PolySpace, qs::PolySpace) = isless(ps.polys, qs.polys)
Base.hash(ps::PolySpace, h::UInt) = hash(ps.polys, hash(0xa3e4cfbf, h))

Base.map(f, ps::PolySpace{D,T}) where {D,T} = PolySpace{D,T}(map(f, ps.polys))

Base.isempty(ps::PolySpace) = isempty(ps.polys)
Base.length(ps::PolySpace) = length(ps.polys)
function Base.issubset(ps::PolySpace{D,T}, qs::PolySpace{D,T}) where {D,T}
    return issubset(Set(ps.polys), Set(qs.polys))
end
function Base.union(ps::PolySpace{D,T}, qs::PolySpace{D,T}) where {D,T}
    return PolySpace{D,T}([ps.polys; qs.polys])
end
function Base.setdiff(ps::PolySpace{D,T}, qs::PolySpace{D,T}) where {D,T}
    return PolySpace{D,T}(collect(setdiff(Set(ps.polys), Set(qs.polys))))
end

Base.zero(::Type{PolySpace{D,T}}) where {D,T} = PolySpace{D,T}()
Base.iszero(ps::PolySpace) = isempty(ps)
# Base.:+(ps::PolySpace) = ps
# Base.:+(ps::PolySpace{D,T}, qs::PolySpace{D,T}) where {D,T} = union(ps, qs)
# Base.:*(a, ps::PolySpace) = map(p -> a * p, ps)
# Base.:*(ps::PolySpace, a) = map(p -> p * a, ps)
# Base.:\(a, ps::PolySpace) = map(p -> a \ p, ps)
# Base.:/(ps::PolySpace, a) = map(p -> p / a, ps)

deriv(ps::PolySpace, dir::Int) = map(p -> deriv(p, dir), ps)
koszul(ps::PolySpace, dir::Int) = map(p -> koszul(p, dir), ps)

################################################################################

function deriv(f::Form{D,R,Poly{D,T}}) where {D,R,T}
    D::Int
    R::Int
    @assert 0 <= R < D
    r = zero(Form{D,R + 1,Poly{D,T}})
    N = length(f)
    for n in 1:N
        bits = Forms.lin2bit(Val(D), Val(R), n)
        for dir in 1:D
            if !bits[dir]
                parity = false
                for d in 1:(dir - 1)
                    parity ⊻= bits[d]
                end
                s = bitsign(parity)
                rbits = Base.setindex(bits, true, dir)
                rn = Forms.bit2lin(Val(D), Val(R + 1), rbits)
                r = Base.setindex(r, r[rn] + s * deriv(f[n], dir), rn)
            end
        end
    end
    return r::Form{D,R + 1,Poly{D,T}}
end

function koszul(f::Form{D,R,Poly{D,T}}) where {D,R,T}
    D::Int
    R::Int
    @assert 0 < R <= D
    r = zero(Form{D,R - 1,Poly{D,T}})
    N = length(f)
    for n in 1:N
        bits = Forms.lin2bit(Val(D), Val(R), n)
        for dir in 1:D
            if bits[dir]
                parity = false
                for d in 1:(dir - 1)
                    parity ⊻= bits[d]
                end
                s = bitsign(parity)
                rbits = Base.setindex(bits, false, dir)
                rn = Forms.bit2lin(Val(D), Val(R - 1), rbits)
                r = Base.setindex(r, r[rn] + s * koszul(f[n], dir), rn)
            end
        end
    end
    return r::Form{D,R - 1,Poly{D,T}}
end

################################################################################

function powers2ind(powers::SVector{D,Int}, maxp::Int) where {D}
    ind = 0
    for d in D:-1:1
        @assert 0 <= powers[d] <= maxp
        ind = (maxp + 1) * ind + powers[d]
    end
    return ind + 1
end

function poly2vec(poly::Poly{D,T}, maxp::Int) where {D,T}
    nrows = (maxp + 1)^D
    I = Int[]
    V = T[]
    for term in poly.terms
        push!(I, powers2ind(term.powers, maxp))
        push!(V, term.coeff)
    end
    return sparsevec(I, V, nrows)
end

export form2vec
function form2vec(form::Form{D,R,Poly{D,T}}, maxp::Int) where {D,R,T}
    stride = (maxp + 1)^D
    nrows = length(form) * stride
    I = Int[]
    V = T[]
    for (n, poly) in enumerate(form)
        offset = (n - 1) * stride
        vec = poly2vec(poly, maxp)
        Is, Vs = findnz(vec)
        append!(I, Is .+ offset)
        append!(V, Vs)
    end
    return sparsevec(I, V, nrows)
end

export forms2mat
function forms2mat(forms::Vector{<:Form{D,R,Poly{D,T}}},
                   maxp::Int) where {D,R,T}
    nrows = length(Form{D,R}) * (maxp + 1)^D
    ncols = length(forms)
    I = Int[]
    J = Int[]
    V = T[]
    for (j, form) in enumerate(forms)
        vec = form2vec(form, maxp)
        Is, Vs = findnz(vec)
        append!(I, Is)
        append!(J, fill(j, length(Is)))
        append!(V, Vs)
    end
    return sparse(I, J, V, nrows, ncols)
end

export maxpower
maximum0(xs) = isempty(xs) ? 0 : maximum(xs)
maxpower(term::Term) = maximum0(term.powers)
maxpower(poly::Poly) = maximum0(maxpower(term) for term in poly.terms)
maxpower(form::Form) = maximum0(maxpower(elt) for elt in form)
maxpower(forms::Vector{<:Form}) = maximum0(maxpower(form) for form in forms)

################################################################################

export Basis
@computed struct Basis{D,R,T}
    forms::Vector{fulltype(Form{D,R,Poly{D,T}})}
    function Basis{D,R,T}(forms::Vector{<:Form{D,R,Poly{D,T}}}) where {D,R,T}
        D::Int
        R::Int
        @assert 0 ≤ R ≤ D
        forms = combine(forms)
        return new{D,R,T}(forms)
    end
end

export compare
function compare(x::Form{D,R,<:Poly{D}}, y::Form{D,R,<:Poly{D}}) where {D,R}
    for (px, py) in zip(x.elts, y.elts)
        c = compare(px, py)
        c != 0 && return c
    end
    return 0
end

function combine(forms::Vector{<:Form{D,R,Poly{D,T}}}) where {D,R,T}
    forms = sort(forms; rev=true)
    unique!(forms)
    filter!(form -> !iszero(form), forms)
    map!(normalize, forms, forms)
    sort!(forms; lt=(x, y) -> compare(x, y) < 0)
    # TODO: filter out linear combinations
    return forms
end

function normalize(form::Form{D,R,Poly{D,T}}) where {D,R,T<:Integer}
    # coeffs = T[term.coeff for elt in form for poly in elt for term in poly]
    coeffs = T[]
    for poly in form
        poly::Poly{D,T}
        for term in poly.terms
            term::Term{D,T}
            push!(coeffs, term.coeff)
        end
    end
    q = gcd(coeffs)
    return map(x -> x ÷ q, form)
end

# function Basis{D,R,T}() where {D,R,T}
#     return Basis{D,R,T}([])
# end

function Base.show(io::IO, mime::MIME"text/plain",
                   b::Basis{D,R,T}) where {D,R,T}
    skiptype = get(io, :typeinfo, Any) <: Basis{D,R,T}
    if !skiptype
        print(io, "Basis{$D,$R,$T}")
    end
    print(io, "[")
    for (i, form) in enumerate(b.forms)
        i > 1 && print(io, ", ")
        show(IOContext(io, :compact => true, :typeinfo => Form{D,R,Poly{D,T}}),
             mime, form)
    end
    print(io, "]")
    return nothing
end

function Base.:(==)(b1::Basis{D,R}, b2::Basis{D,R}) where {D,R}
    # return b1.forms == b2.forms
    # Fast paths
    b1.forms == b2.forms && return true
    isempty(b1.forms) != isempty(b2.forms) && return false
    # Use subset relations
    return issubset(b1, b2) && issubset(b2, b1)
end

Base.isequal(b1::Basis, b2::Basis) = isequal(b1.forms, b2.forms)
Base.isless(b1::Basis, b2::Basis) = isless(b1.forms, b2.forms)
Base.hash(b::Basis, h::UInt) = hash(b.forms, hash(0xe8439c4b, h))

Base.zero(::Type{<:Basis{D,R,T}}) where {D,R,T} = Basis{D,R,T}([])
Base.zero(::Basis{D,R,T}) where {D,R,T} = zero(Basis{D,R,T})
Base.iszero(basis::Basis) = isempty(basis.forms)

function Forms.unit(form::Form{D,R,Poly{D,T}}) where {D,R,T}
    return Basis{D,R,T}([form])
end

export tensorsum, ⊕
function Forms.tensorsum(basis1::Basis{D,R,T},
                         basis2::Basis{D,R,T}) where {D,R,T}
    return Basis{D,R,T}([basis1.forms; basis2.forms])
end

function tensordiff end
const ⊖ = tensordiff
export tensordiff, ⊖
function tensordiff(basis1::Basis{D,R,T}, basis2::Basis{D,R,T}) where {D,R,T}
    forms2 = Set(basis2.forms)
    forms = fulltype(Form{D,R,Poly{D,T}})[]
    for form in basis1.forms
        if form ∉ forms2
            push!(forms, form)
        end
    end
    return Basis{D,R,T}(forms)
end

Base.isempty(basis::Basis) = isempty(basis.forms)
Base.length(basis::Basis) = length(basis.forms)

function Base.in(form::Form{D,R,Poly{D,T}}, basis::Basis{D,R,T}) where {D,R,T}
    # Fast paths
    isempty(basis.forms) && return false
    any(form == f for f in basis.forms) && return true
    # Convert representation to sparse vectors/matrices
    maxp = max(maxpower(form), maxpower(basis.forms))
    fvec = form2vec(form, maxp)
    bmat = forms2mat(basis.forms, maxp)
    # Sparse solver doesn't handle sparse vectors
    fvec = collect(fvec)
    # Is there a linear combination of the vectors in bmat that yields fvec?
    #     fⁱ = bⁱⱼ aʲ
    if T <: Integer
        # Use floating point numbers for integer polynomials
        avec = float.(bmat) \ float.(fvec)
        avec = rationalize.(avec; tol=sqrt(eps()))
    else
        @assert false
    end
    rvec = fvec - bmat * avec
    return iszero(rvec)
end

function Base.issubset(basis1::Basis{D,R}, basis2::Basis{D,R}) where {D,R}
    isempty(basis1.forms) && return true
    isempty(basis2.forms) && return false
    return all(f in basis2 for f in basis1.forms)
end

# Unused
function make_basis1(::Type{<:Basis{D,R,T}}, p::Int, cond) where {D,R,T}
    D::Int
    R::Int
    @assert 0 <= R <= D
    @assert p >= -1             # -1 is the empty basis
    polys = Poly{D,T}[]
    for i0 in
        CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> p, D))
        i = SVector{D,Int}(i0.I)
        if cond(sum(i))
            push!(polys, Poly{D,T}([Term{D,T}(i, one(T))]))
        end
    end
    N = length(Form{D,R})
    forms = fulltype(Form{D,R,Poly{D,T}})[zero(Form{D,R,Poly{D,T}})]
    for n in 1:N
        newforms = fulltype(Form{D,R,Poly{D,T}})[]
        for f in forms
            for p in polys
                push!(newforms, setindex(f, p, n))
            end
        end
        forms = newforms
    end

    return Basis{D,R,T}(forms)
end

# Unused
function make_basis2(::Type{<:Basis{D,R,T}}, p::Int, cond) where {D,R,T}
    D::Int
    R::Int
    @assert 0 <= R <= D
    @assert p >= -1             # -1 is the empty basis
    N = length(Form{D,R})
    forms = fulltype(Form{D,R,Poly{D,T}})[zero(Form{D,R,Poly{D,T}})]
    for i0 in
        CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> p, D))
        i = SVector{D,Int}(i0.I)
        if cond(sum(i))
            poly = Poly{D,T}([Term{D,T}(i, one(T))])
            form = Form{D,R,Poly{D,T}}(ntuple(n -> poly, N))
            push!(forms, form)
        end
    end
    return Basis{D,R,T}(forms)
end

function make_basis(::Type{<:Basis{D,R,T}}, p::Int, cond) where {D,R,T}
    D::Int
    R::Int
    @assert 0 <= R <= D
    @assert p >= -1             # -1 is the empty basis
    polys = Poly{D,T}[]
    for i0 in
        CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> p, D))
        i = SVector{D,Int}(i0.I)
        if cond(sum(i))
            push!(polys, Poly{D,T}([Term{D,T}(i, one(T))]))
        end
    end
    N = length(Form{D,R})
    forms = fulltype(Form{D,R,Poly{D,T}})[]
    for n in 1:N
        for poly in polys
            form = Form{D,R,Poly{D,T}}(ntuple(m -> m == n ? poly :
                                                   zero(Poly{D,T}), N))
            push!(forms, form)
        end
    end
    return Basis{D,R,T}(forms)
end

export full_basis
function full_basis(::Type{<:Basis{D,R,T}}, p::Int) where {D,R,T}
    return make_basis(Basis{D,R,T}, p, <=(p))
end

export homogeneous_basis
function homogeneous_basis(::Type{<:Basis{D,R,T}}, p::Int) where {D,R,T}
    return make_basis(Basis{D,R,T}, p, ==(p))
end

function deriv(basis::Basis{D,R,T}) where {D,R,T}
    D::Int
    R::Int
    @assert 0 <= R < D
    dforms = fulltype(Form{D,R + 1,Poly{D,T}})[]
    for form in basis.forms
        dform = deriv(form)
        !iszero(dform) && push!(dforms, dform)
    end
    return Basis{D,R + 1,T}(dforms)
end

function koszul(basis::Basis{D,R,T}) where {D,R,T}
    D::Int
    R::Int
    @assert 0 < R <= D
    κforms = fulltype(Form{D,R - 1,Poly{D,T}})[]
    for form in basis.forms
        κform = koszul(form)
        !iszero(κform) && push!(κforms, κform)
    end
    return Basis{D,R - 1,T}(κforms)
end

################################################################################

# function bernstein(::Type{<:Basis{D,0,T}}, p::Int) where {D,T}
#     D::Int
#     R = 0
#     @assert 0 <= R <= D
#     @assert p >= 0
#     polys = fulltype(Form{D,R,Poly{D,T}})[]
#     for i0 in
#         CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> p, D))
#         i = SVector{D,Int}(i.I)
#         if sum(i) == p
#             poly = Poly{D,T}([Term{D,T}(i, one(T))])
#             push!(polys, Form{D,R}((poly,)))
#         end
#     end
#     return Basis{D,T}(polys)
# end
#
# function bernstein(::Type{<:Basis{D,R,T}}, p::Int) where {D,R,T}
#     D::Int
#     R::Int
#     @assert 0 <= R <= D
#     @assert p >= 0
#     polys = Poly{D,T}[]
#     for i0 in
#         CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(ntuple(d -> p, D))
#         i = SVector{D,Int}(i.I)
#         if sum(i) == p
#             push!(polys, Poly{D,T}([Term{D,T}(i, one(T))]))
#         end
#     end
#     return Basis{D,T}(polys)
# end

################################################################################

export polynomial_complex
"""
polynomial_complex

Pᵣ Λᵏ
"""
function polynomial_complex(::Val{D}, ::Type{T}, p::Int) where {D,T}
    D::Int
    @assert 0 <= D
    @assert p >= -1

    cc = Dict{Int,Basis{D,R,T} where R}()

    for R in 0:D
        cc[R] = full_basis(Basis{D,R,T}, max(-1, p - R))::Basis{D,R,T}
    end

    return cc
end

export trimmed_polynomial_complex
"""
trimmed_polynomial_complex

Pᵣ⁻ Λᵏ
"""
function trimmed_polynomial_complex(::Val{D}, ::Type{T}, p::Int) where {D,T}
    D::Int
    @assert 0 <= D
    @assert p >= -1

    cc = Dict{Int,Basis{D,R,T} where R}()

    for R in 0:D
        if R == 0
            b = full_basis(Basis{D,R,T}, p)
        else
            b = full_basis(Basis{D,R,T}, max(-1, p - 1))
            if R + 1 <= D
                b = b ⊕ koszul(homogeneous_basis(Basis{D,R + 1,T}, p - 1))
            end
            # if R - 1 >= 0
            #     b = b ⊕ deriv(homogenous_basis(Basis{D,R - 1,T}, p + 1))
            # end
        end
        cc[R] = b::Basis{D,R,T}
    end

    return cc
end

export extended_trimmed_polynomial_complex
"""
extended_trimmed_polynomial_complex

This is experimental; note that this is not actually a complex.
"""
function extended_trimmed_polynomial_complex(::Val{D}, ::Type{T},
                                             p::Int) where {D,T}
    D::Int
    @assert 0 <= D
    @assert p >= -1

    cc = Dict{Int,Basis{D,R,T} where R}()

    for R in 0:D
        if R == 0
            b = full_basis(Basis{D,R,T}, p)
        else
            b = full_basis(Basis{D,R,T}, max(-1, p - 1))
            if R + 1 <= D
                b = b ⊕ koszul(homogeneous_basis(Basis{D,R + 1,T}, p - 1))
            end
            if R - 1 >= 0
                b = b ⊕ deriv(homogeneous_basis(Basis{D,R - 1,T}, p + 1))
            end
        end
        cc[R] = b::Basis{D,R,T}
    end

    return cc
end

################################################################################

export whitney
"""
whitney

Calculate Whitney forms
"""
function whitney(::Val{D1}, inds::SVector{N,Int}) where {D1,N}
    D1::Int
    N::Int
    @assert 1 <= D1
    @assert 1 <= N <= D1
    R = N - 1
    @assert all(1 <= inds[n] <= D1 for n in 1:N)
    @assert all(inds[n] < inds[n + 1] for n in 1:(N - 1))
    ϕ = zero(Form{D1,R,Poly{D1,Int}})
    q = factorial(N - 1)
    for n in 1:N
        inds′ = deleteat(inds, n)
        f = unit(Form{D1,R,Int}, inds′)
        p = q * bitsign(n - 1) * unit(Poly{D1,Int}, inds[n])
        ϕ += p * f
    end
    return ϕ::Form{D1,R,Poly{D1,Int}}
end
function whitney(::Val{D1}, inds::NTuple{N,Int}) where {D1,N}
    return whitney(Val(D1), SVector{N,Int}(inds))
end
function whitney(::Type{Basis{D1,R,Int}}) where {D1,R}
    D1::Int
    R::Int
    @assert 1 <= D1
    @assert 0 <= R <= D1 - 1
    nelts = binomial(D1, R + 1)
    forms = [whitney(Val(D1), Forms.lin2lst(Val(D1), Val(R + 1), n))
             for n in 1:nelts]
    return Basis{D1,R,Int}(forms)
end

export whitney_support
"""
whitney_support

Determine on what part of a simplex (vertex, edge, face, ...) a
particular Whitney forms should live.
"""
function whitney_support(::Val{D1}, inds::SVector{N,Int}) where {D1,N}
    return inds
end

################################################################################

export barycentric2cartesian
function barycentric2cartesian(λterm::Term{D1,T}) where {D1,T}
    D1::Int
    @assert 1 <= D1
    D = D1 - 1
    # This assumes a standard simplex:
    #     λᵢ = xᵢ
    #     λₙ = 1 - Σᵢ λᵢ
    λn = one(Poly{D,T})
    if D > 0
        λn -= sum(unit(Poly{D,T}, d) for d in 1:D)
    end
    xpoly = Poly{D,T}([Term{D,T}(deleteat(λterm.powers, D1), λterm.coeff)])
    xpoly *= λn^λterm.powers[D1]
    return xpoly::Poly{D,T}
end
function barycentric2cartesian(λpoly::Poly{D1,T}) where {D1,T}
    D1::Int
    @assert 1 <= D1
    D = D1 - 1
    isempty(λpoly.terms) && return zero(Poly{D,T})
    xpoly = sum(barycentric2cartesian(λterm) for λterm in λpoly.terms)
    return xpoly::Poly{D,T}
end
function barycentric2cartesian(λform::Form{D1,R,Poly{D1,T}}) where {D1,R,T}
    D1::Int
    R::Int
    @assert 1 <= D1
    D = D1 - 1
    @assert 0 <= R <= D
    # This assumes a standard simplex:
    #     λᵢ = xᵢ
    #     λₙ = 1 - Σᵢ λᵢ
    #     dλᵢ = dxᵢ
    #     dλₙ = -Σᵢ dλᵢ
    if D == 0
        λpoly = λform[]
        xpoly = barycentric2cartesian(λpoly)
        xform = Form{D,R,Poly{D,T}}((xpoly,))
    else
        dλn = -sum(unit(Form{D,1,Poly{D,T}}, d) for d in 1:D)
        xform = zero(Form{D,R,Poly{D,T}})
        for (λlin, λpoly) in enumerate(λform)
            λbits = Forms.lin2bit(Val(D1), Val(R), λlin)
            xbits = deleteat(λbits, D1)
            if λbits[D1]
                xlin = Forms.bit2lin(Val(D), Val(R - 1), xbits)
                xunit = unit(Form{D,R - 1,Poly{D,T}}, xlin) ∧ dλn
            else
                xlin = Forms.bit2lin(Val(D), Val(R), xbits)
                xunit = unit(Form{D,R,Poly{D,T}}, xlin)
            end
            xpoly = barycentric2cartesian(λpoly)
            xform += xpoly * xunit
        end
    end
    return xform::Form{D,R,Poly{D,T}}
end
function barycentric2cartesian(λbasis::Basis{D1,R,T}) where {D1,R,T}
    D1::Int
    R::Int
    @assert 1 <= D1
    D = D1 - 1
    @assert 0 <= R <= D
    return Basis{D,R,T}(barycentric2cartesian.(λbasis.forms))
end

end


text:
"""
    Descent(η = 1f-1)

Classic gradient descent optimiser with learning rate `η`.
For each parameter `p` and its gradient `dp`, this runs `p -= η*dp`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
"""
struct Descent{T}
  eta::T
end
Descent() = Descent(1f-1)

init(o::Descent, x::AbstractArray) = nothing

function apply(o::Descent, state, x, dx)
  η = convert(eltype(dx), o.eta)

  return state, dx .* η
end

(o::Descent)(state, m, dm) = update(o, state, m, dm)

"""
    Momentum(η = 1f-2, ρ = 9f-1)

Gradient descent optimizer with learning rate `η` and momentum `ρ`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Momentum (`ρ`): Controls the acceleration of gradient descent in the
                  prominent direction, in effect dampening oscillations.
"""
struct Momentum{T}
  eta::T
  rho::T
end
Momentum(η = 1f-2, ρ = 9f-1) = Momentum{typeof(η)}(η, ρ)

init(o::Momentum, x::AbstractArray) = zero(x)

function apply(o::Momentum, state, x, dx)
  η, ρ, v = o.eta, o.rho, state
  @. v = ρ * v - η * dx

  return v, -v
end

(o::Momentum)(state, m, dm) = update(o, state, m, dm)

"""
    Nesterov(η = 1f-3, ρ = 9f-1)

Gradient descent optimizer with learning rate `η` and Nesterov momentum `ρ`.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Nesterov momentum (`ρ`): Controls the acceleration of gradient descent in the
                           prominent direction, in effect dampening oscillations.
"""
struct Nesterov{T}
  eta::T
  rho::T
end
Nesterov(η = 1f-3, ρ = 9f-1) = Nesterov{typeof(η)}(η, ρ)

init(o::Nesterov, x::AbstractArray) = zero(x)

(o::Nesterov)(state, m, dm) = update(o, state, m, dm)

function apply(o::Nesterov, state, x, dx)
  η, ρ, v = o.eta, o.rho, state
  d = @. ρ^2 * v - (1+ρ) * η * dx
  @. v = ρ * v - η * dx

  return v, -d
end

"""
    RMSProp(η = 1f-3, ρ = 9f-1, ϵ = eps(typeof(η)))

Optimizer using the
[RMSProp](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
algorithm. Often a good choice for recurrent networks. Parameters other than learning rate
generally don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Momentum (`ρ`): Controls the acceleration of gradient descent in the
                  prominent direction, in effect dampening oscillations.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct RMSProp{T}
  eta::T
  rho::T
  epsilon::T
end
RMSProp(η = 1f-3, ρ = 9f-1, ϵ = eps(typeof(η))) = RMSProp{typeof(η)}(η, ρ, ϵ)

init(o::RMSProp, x::AbstractArray) = zero(x)

function apply(o::RMSProp, state, x, dx)
  η, ρ, ϵ, acc = o.eta, o.rho, o.epsilon, state
  @. acc = ρ * acc + (1 - ρ) * dx^2
  dx = @. dx * (η / (sqrt(acc) + ϵ))

  return acc, dx
end

(o::RMSProp)(state, m, dm) = update(o, state, m, dm)

"""
    ADAM(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η)))

[ADAM](https://arxiv.org/abs/1412.6980) optimiser.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct ADAM{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
ADAM(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η))) = ADAM{typeof(η)}(η, β, ϵ)

init(o::ADAM, x::AbstractArray) = (zero(x), zero(x), o.beta)

(o::ADAM)(state, m, dm) = update(o, state, m, dm)

function apply(o::ADAM{T}, state, x, dx) where T
  η, β, ϵ = o.eta, o.beta, o.epsilon
  mt, vt, βt = state

  @. mt = β[1] * mt + (one(T) - β[1]) * dx
  @. vt = β[2] * vt + (one(T) - β[2]) * dx ^ 2
  dx = @. mt / (one(T) - βt[1]) / (sqrt(vt / (one(T) - βt[2])) + ϵ) * η

  return (mt, vt, βt .* β), dx
end

"""
    RADAM(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η)))

[Rectified ADAM](https://arxiv.org/abs/1908.03265) optimizer.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct RADAM{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
RADAM(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η))) = RADAM{typeof(η)}(η, β, ϵ)

init(o::RADAM, x::AbstractArray) = (zero(x), zero(x), o.beta, 1)

(o::RADAM)(state, m, dm) = update(o, state, m, dm)

function apply(o::RADAM, state, x, dx)
  η, β, ϵ = o.eta, o.beta, o.epsilon
  ρ∞ = 2/(1-β[2])-1

  mt, vt, βt, t = state

  @. mt = β[1] * mt + (1 - β[1]) * dx
  @. vt = β[2] * vt + (1 - β[2]) * dx^2
  ρ = ρ∞ - 2*t * βt[2] / (1 - βt[2])
  if ρ > 4
    r = sqrt((ρ - 4) * (ρ - 2) * ρ∞/((ρ∞ - 4) * (ρ∞ - 2) * ρ))
    dx = @. mt / (1 - βt[1]) / (sqrt(vt / (1 - βt[2])) + ϵ) * η * r
  else
    dx = @. mt / (1 - βt[1]) * η
  end

  return (mt, vt, βt .* β, t + 1), dx
end

"""
    AdaMax(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η)))

[AdaMax](https://arxiv.org/abs/1412.6980) is a variant of ADAM based on the ∞-norm.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct AdaMax{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
AdaMax(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η))) = AdaMax{typeof(η)}(η, β, ϵ)

init(o::AdaMax, x::AbstractArray) = (zero(x), zero(x), o.beta)

(o::AdaMax)(state, m, dm) = update(o, state, m, dm)

function apply(o::AdaMax, state, x, dx)
  η, β, ϵ = o.eta, o.beta, o.epsilon

  mt, ut, βt = state

  @. mt = β[1] * mt + (1 - β[1]) * dx
  @. ut = max(β[2] * ut, abs(dx))
  dx = @. (η/(1 - βt[1])) * mt/(ut + ϵ)

  return (mt, ut, βt .* β), dx
end

"""
    OADAM(η = 1f-3, β = (5f-1, 9f-1), ϵ = eps(typeof(η)))

[OADAM](https://arxiv.org/abs/1711.00141) (Optimistic ADAM)
is a variant of ADAM adding an "optimistic" term suitable for adversarial training.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct OADAM{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
OADAM(η = 1f-3, β = (5f-1, 9f-1), ϵ = eps(typeof(η))) = OADAM{typeof(η)}(η, β, ϵ)

init(o::OADAM, x::AbstractArray) = (zero(x), zero(x), o.beta, zero(x))

(o::OADAM)(state, m, dm) = update(o, state, m, dm)

function apply(o::OADAM, state, x, dx)
  η, β, ϵ = o.eta, o.beta, o.epsilon

  mt, vt, βt, dx_ = state

  @. mt = β[1] * mt + (1 - β[1]) * dx
  @. vt = β[2] * vt + (1 - β[2]) * dx^2
  @. dx = -dx_
  @. dx_ = η * mt / (1 - βt[1]) / (sqrt(vt / (1 - βt[2])) + ϵ)
  dx = @. dx + 2*dx_

  return (mt, vt, βt .* β, dx_), dx
end

"""
    ADAGrad(η = 1f-1, ϵ = eps(typeof(η)))

[ADAGrad](http://www.jmlr.org/papers/volume12/duchi11a/duchi11a.pdf) optimizer. It has
parameter specific learning rates based on how frequently it is updated.
Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct ADAGrad{T}
  eta::T
  epsilon::T
end
ADAGrad(η = 1f-1, ϵ = eps(typeof(η))) = ADAGrad{typeof(η)}(η, ϵ)

init(o::ADAGrad, x::AbstractArray) = fill!(similar(x), o.epsilon)

(o::ADAGrad)(state, m, dm) = update(o, state, m, dm)

function apply(o::ADAGrad, state, x, dx)
  η, ϵ = o.eta, o.epsilon
  acc = state

  @. acc += dx^2
  dx = @. dx * η / (sqrt(acc) + ϵ)

  return acc, dx
end

"""
    ADADelta(ρ = 9f-1, ϵ = eps(typeof(ρ)))

[ADADelta](https://arxiv.org/abs/1212.5701) is a version of ADAGrad adapting its learning
rate based on a window of past gradient updates.
Parameters don't need tuning.

# Parameters
- Rho (`ρ`): Factor by which the gradient is decayed at each time step.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct ADADelta{T}
  rho::T
  epsilon::T
end
ADADelta(ρ = 9f-1, ϵ = eps(typeof(ρ))) = ADADelta{typeof(ρ)}(ρ, ϵ)

init(o::ADADelta, x::AbstractArray) = (zero(x), zero(x))

(o::ADADelta)(state, m, dm) = update(o, state, m, dm)

function apply(o::ADADelta, state, x, dx)
  ρ, ϵ = o.rho, o.epsilon
  acc, Δacc = state

  @. acc = ρ * acc + (1 - ρ) * dx^2
  # DON'T remove epsilon from numerator
  # or even out of the square roots
  dx = @. dx * sqrt(Δacc + ϵ) / sqrt(acc + ϵ)
  @. Δacc = ρ * Δacc + (1 - ρ) * dx^2

  return (acc, Δacc), dx
end

"""
    AMSGrad(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η)))

The [AMSGrad](https://openreview.net/forum?id=ryQu7f-RZ) version of the ADAM
optimiser. Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct AMSGrad{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
AMSGrad(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η))) = AMSGrad{typeof(η)}(η, β, ϵ)

init(o::AMSGrad, x::AbstractArray) =
  (fill!(similar(x), o.epsilon), fill!(similar(x), o.epsilon), fill!(similar(x), o.epsilon))

(o::AMSGrad)(state, m, dm) = update(o, state, m, dm)

function apply(o::AMSGrad, state, x, dx)
  η, β, ϵ = o.eta, o.beta, o.epsilon

  mt, vt, v̂t = state

  @. mt = β[1] * mt + (1 - β[1]) * dx
  @. vt = β[2] * vt + (1 - β[2]) * dx ^ 2
  @. v̂t = max(v̂t, vt)
  dx = @. η * mt / (sqrt(v̂t) + ϵ)

  return (mt, vt, v̂t), dx
end

"""
    NADAM(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η)))

[NADAM](https://openreview.net/forum?id=OM0jvwB8jIp57ZJjtNEZ) is a Nesterov variant of ADAM.
Parameters don't need tuning.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
struct NADAM{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
NADAM(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η))) = NADAM{typeof(η)}(η, β, ϵ)

init(o::NADAM, x::AbstractArray) = (zero(x), zero(x), o.beta)

(o::NADAM)(state, m, dm) = update(o, state, m, dm)

function apply(o::NADAM, state, x, dx)
  η, β, ϵ = o.eta, o.beta, o.epsilon

  mt, vt, βt = state

  @. mt = β[1] * mt + (1 - β[1]) * dx
  @. vt = β[2] * vt + (1 - β[2]) * dx^2
  dx = @. (β[1] * mt / (1 - β[1] * βt[1]) + (1 - β[1]) * dx / (1 - βt[1])) /
          (sqrt(vt * β[2] / (1 - βt[2])) + ϵ) * η

  return (mt, vt, βt .* β), dx
end

"""
    ADAMW(η = 1f-3, β = (9f-1, 9.99f-1), γ = 0, ϵ = eps(typeof(η)))

[ADAMW](https://arxiv.org/abs/1711.05101) is a variant of ADAM fixing (as in repairing) its
weight decay regularization.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Weight decay (`γ`): Decay applied to weights during optimisation.
- Machine epsilon (`ϵ`): Constant to prevent division by zero
                         (no need to change default)
"""
ADAMW(η = 1f-3, β = (9f-1, 9.99f-1), γ = 0, ϵ = eps(typeof(η))) =
  OptimiserChain(ADAM{typeof(η)}(η, β, ϵ), WeightDecay(γ))

"""
    AdaBelief(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η)))

The [AdaBelief](https://arxiv.org/abs/2010.07468) optimiser is a variant of the well-known
ADAM optimiser.

# Parameters
- Learning rate (`η`): Amount by which gradients are discounted before updating
                       the weights.
- Decay of momentums (`β::Tuple`): Exponential decay for the first (β1) and the
                                   second (β2) momentum estimate.
- Machine epsilon (`ϵ::Float32`): Constant to prevent division by zero
                                  (no need to change default)
"""
struct AdaBelief{T}
  eta::T
  beta::Tuple{T, T}
  epsilon::T
end
AdaBelief(η = 1f-3, β = (9f-1, 9.99f-1), ϵ = eps(typeof(η))) = AdaBelief{typeof(η)}(η, β, ϵ)

init(o::AdaBelief, x::AbstractArray) = (zero(x), zero(x))

(o::AdaBelief)(state, m, dm) = update(o, state, m, dm)

function apply(o::AdaBelief, state, x, dx)
  η, β, ϵ = o.eta, o.beta, o.epsilon
  mt, st = state

  @. mt = β[1] * mt + (1 - β[1]) * dx
  @. st = β[2] * st + (1 - β[2]) * (dx - mt)^2
  dx = @. η * mt / (sqrt(st) + ϵ)

  return (mt, st), dx
end

"""
    WeightDecay(γ = 5f-4)

Decay weights by `γ`.

# Parameters
- Weight decay (`γ`): Decay applied to weights during optimisation.
"""
struct WeightDecay{T}
  wd::T
end
WeightDecay() = WeightDecay(5f-4)

init(o::WeightDecay, x::AbstractArray) = nothing

(o::WeightDecay)(state, m, dm) = update(o, state, m, dm)

function apply(o::WeightDecay, state, x, dx)
  dx = @. dx + o.wd * x

  return state, dx
end

"""
    OptimiserChain(opts...)

Compose a chain (sequence) of optimisers so that each `opt` in `opts`
updates the gradient in the order specified.
"""
struct OptimiserChain{O}
  opts::O
end
OptimiserChain(opts...) = OptimiserChain(opts)

init(o::OptimiserChain, x::AbstractArray) = [init(opt, x) for opt in o.opts]

(o::OptimiserChain)(state, m, dms...) = update(o, state, m, dms...)

function apply(o::OptimiserChain, states, x, dx, dxs...)
  new_states = similar(states)
  for (i, (opt, state)) in enumerate(zip(o.opts, states))
    new_states[i], dx = apply(opt, state, x, dx, dxs...)
  end

  return new_states, dx
end

for Opt in (:Descent, :ADAM, :Momentum, :Nesterov, :RMSProp,
            :ADAGrad, :AdaMax, :ADADelta, :AMSGrad, :NADAM,
            :RADAM, :OADAM, :AdaBelief)
  @eval function $Opt(m::$Opt; kwargs...)
    fs = fieldnames($Opt)
    args = NamedTuple{fs}(f in keys(kwargs) ? getindex(kwargs, f) : getfield(m, f) for f in fs)
    $Opt(args...)
  end
end


text:
@testset "invalid row bounds" begin
    v = [1, 2, 3]
    @test_throws BoundsError RowProjectionVector(v, 0)
    @test_throws BoundsError RowProjectionVector(v, 4)

    M = [1 2 3; 4 5 6; 7 8 9]
    @test_throws BoundsError RowProjectionMatrix(M, 0)
    @test_throws BoundsError RowProjectionMatrix(M, 0, 1)
    @test_throws BoundsError RowProjectionMatrix(M, 1, 0)
    @test_throws BoundsError RowProjectionMatrix(M, 4)
    @test_throws BoundsError RowProjectionMatrix(M, 1, 4)
    @test_throws BoundsError RowProjectionMatrix(M, 4, 1)
end

@testset "single row" begin
    v = [1, 2, 3]
    for row ∈ 1:length(v)
        u = RowProjectionVector(v, row)
        @test u isa RowProjectionVector{Int}
        @test size(u) == (1,)
        @test length(u) == 1
        @test_throws BoundsError u[0]
        @test_throws BoundsError u[0] = 41
        @test u[] == v[row]
        @test u[1] == v[row]
        u[] = 10 + row
        @test u[] == 10 + row
        u[1] = 100 + row
        @test u[1] == 100 + row
        @test_throws BoundsError u[length(u) + 1]
        @test_throws BoundsError u[length(u) + 1] = 41
    end

    M = [1 2 3; 4 5 6; 7 8 9]
    for row ∈ 1:size(M, 1)
        A = RowProjectionMatrix(M, row)
        @test A isa RowProjectionMatrix{Int}
        @test size(A) == (1, size(M, 2))
        @test length(A) == size(M, 1)
        @test_throws BoundsError A[]
        @test_throws BoundsError A[] = 41
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0] = 41
        for col ∈ 1:size(M, 2)
            @test A[col] == M[row, col]
            @test A[1, col] == M[row, col]

            A[col] = row + col
            @test A[col] == row + col
            @test A[1, col] == row + col
            @test M[row, col] == row + col

            A[1, col] = 10 + row + col
            @test A[col] == 10 + row + col
            @test A[1, col] == 10 + row + col
            @test M[row, col] == 10 + row + col

            M[row, col] = 100 + row + col
            @test A[col] == 100 + row + col
            @test A[1, col] == 100 + row + col
        end
        @test_throws BoundsError A[length(A) + 1]
        @test_throws BoundsError A[length(A) + 1] = 41
    end
end

@testset "row subset" begin
    v = [1, 2, 3]
    u = RowProjectionVector(v, 3, 1)
    @test u isa RowProjectionVector{Int}
    @test size(u) == (2,)
    @test length(u) == 2
    @test_throws BoundsError u[0]
    @test_throws BoundsError u[0] = 41
    @test u[1] == v[3]
    u[1] = 11
    @test u[1] == v[3] == 11
    @test u[2] == v[1]
    u[2] = 12
    @test u[2] == v[1] == 12
    @test_throws BoundsError u[length(u) + 1]
    @test_throws BoundsError u[length(u) + 1] = 41

    M = [1 2 3; 4 5 6; 7 8 9]
    A = RowProjectionMatrix(M, 3, 1)
    @test A isa RowProjectionMatrix{Int}
    @test size(A) == (2, size(M, 2))
    @test length(A) == 2 * size(M, 2)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[] = 41
    @test_throws BoundsError A[0]
    @test_throws BoundsError A[0] = 41
    @test A[1, :] == M[3, :]
    A[1, :] .= [11, 12, 13]
    @test A[1, :] == M[3, :] == [11, 12, 13]
    @test A[2, :] == M[1, :]
    A[2, :] .= [21, 22, 23]
    @test A[2, :] == M[1, :] == [21, 22, 23]
    @test_throws BoundsError A[length(A) + 1]
    @test_throws BoundsError A[length(A) + 1] = 41
end

@testset "row bag" begin
    v = [1, 2, 3]
    u = RowProjectionVector(v, 2, 2)
    @test u isa RowProjectionVector{Int}
    @test size(u) == (2,)
    @test length(u) == 2
    @test_throws BoundsError u[0]
    @test_throws BoundsError u[0] = 41
    @test u[1] == u[2] == v[2]
    u[1] = 11
    @test u[1] == u[2] == v[2] == 11
    u[2] = 12
    @test u[1] == u[2] == v[2] == 12
    @test_throws BoundsError u[length(u) + 1]
    @test_throws BoundsError u[length(u) + 1] = 41

    M = [1 2 3; 4 5 6; 7 8 9]
    A = RowProjectionMatrix(M, 2, 2)
    @test A isa RowProjectionMatrix{Int}
    @test size(A) == (2, size(M, 2))
    @test length(A) == 2 * size(M, 2)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[] = 41
    @test_throws BoundsError A[0]
    @test_throws BoundsError A[0] = 41
    @test A[1, :] == A[2, :] == M[2, :]
    A[1, :] .= [11, 12, 13]
    @test A[1, :] == A[2, :] == M[2, :] == [11, 12, 13]
    A[2, :] .= [21, 22, 23]
    @test A[1, :] == A[2, :] == M[2, :] == [21, 22, 23]
    @test_throws BoundsError A[length(A) + 1]
    @test_throws BoundsError A[length(A) + 1] = 41
end

@testset "full projection" begin
    v = [1, 2, 3]
    u = RowProjectionVector(v, 1, 2, 3)
    @test u isa RowProjectionVector{Int}
    @test size(u) == size(v)
    @test length(u) == length(v)
    @test_throws BoundsError u[0]
    @test_throws BoundsError u[0] = 41
    @test u == v
    for row in 1:length(v)
        u[row] = 10 + row
        @test u[row] == 10 + row
        @test v[row] == 10 + row
    end
    @test u == v
    @test_throws BoundsError u[length(u) + 1]
    @test_throws BoundsError u[length(u) + 1] = 41

    M = [1 2 3; 4 5 6; 7 8 9]
    A = RowProjectionMatrix(M, 1, 2, 3)
    @test A isa RowProjectionMatrix{Int}
    @test size(A) == size(M)
    @test length(A) == length(M)
    @test_throws BoundsError A[]
    @test_throws BoundsError A[] = 41
    @test_throws BoundsError A[0]
    @test_throws BoundsError A[0] = 41
    @test A == M
    for row ∈ 1:size(M, 1)
        for col ∈ 1:size(M, 2)
            @test A[row, col] == M[row, col]

            A[row, col] = 10 + row + col
            @test A[row, col] == 10 + row + col
            @test M[row, col] == 10 + row + col
        end
    end
    @test_throws BoundsError A[length(A) + 1]
    @test_throws BoundsError A[length(A) + 1] = 41
end


text:
using Boltzmann
using Distributions
using Base.Test

X = rand(1000, 200)

function brbm_smoke_test()
    model = BernoulliRBM(1000, 500, (1000,1))
    fit(model, X)
end

function grbm_smoke_test()
    model = GRBM(1000, 500, (1000,1))
    fit(model, X)
end

function conf_smoke_test()
    model = RBM(Normal, Normal, 1000, 500,(1000,1))
    fit(model, X)
end

function init_smoke_test()
    model = RBM(Normal, Normal, 1000, 500, (1000,1);dataset=X)
    fit(model, X)
end

info("brbm_smoke_test()")
brbm_smoke_test()
info("grbm_smoke_test()")
grbm_smoke_test()
info("conf_smoke_test()")
conf_smoke_test()
info("init_smoke_test()")
init_smoke_test()


text:
function transform_regex_ast(ast::ASTNode)
  # flatten :or nodes with only child
  process_regex_or_nodes(ast)
  # flatten :and nodes with only child, and concatenate adjacent terminal children
  process_regex_and_nodes(ast)
  # add start node as reference to root to enable reachability check
  push!(ast.refs, ast.children[1])
end

function process_regex_or_nodes(parentnode::ASTNode)
  for idx in 1:length(parentnode.children)
    node = parentnode.children[idx]
    process_regex_or_nodes(node)
    if node.func == :or
      if length(node.children) == 1
        # flatten if a single child since no need for a choice
        append!(node.children[1].warnings, parentnode.children[idx].warnings)
        append!(node.children[1].source, parentnode.children[idx].source)
        parentnode.children[idx] = node.children[1]
      end
    end
  end
end

function process_regex_and_nodes(parentnode::ASTNode)
  for idx in 1:length(parentnode.children)
    node = parentnode.children[idx]
    process_regex_and_nodes(node)
    if node.func == :and
      # concatenate adjacent terminals
      i = 1
      while i < length(node.children)
        if (node.children[i].func == :terminal) && (node.children[i+1].func == :terminal)
          append!(node.children[i].warnings, node.children[i+1].warnings)
          append!(node.children[i].source, node.children[i+1].source)
          node.children[i].args[:value] *= node.children[i+1].args[:value]
          deleteat!(node.children, i+1)
        else
          i += 1
        end
      end
      # flatten if a single child since no need for and rule
      if length(node.children) == 1
        append!(node.children[1].warnings, parentnode.children[idx].warnings)
        append!(node.children[1].source, parentnode.children[idx].source)
        parentnode.children[idx] = node.children[1]
      end
    end
  end
end


text:
using MPAC
using Test

@testset "MPAC.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: syscop/MPAC

max_issues_repo_name: syscop/MPAC

max_forks_repo_name: syscop/MPAC