text:
function synthetic_spiral_data(N, degrees, start, σ)
    s = deg2rad(start)
    N₀ = floor(Int, N/2)

    X = Matrix{Float64}(undef, 2, N)
    t = [ones(Int, N₀); -ones(Int, N-N₀)]

    for n in 1:N₀
        θ = s + sqrt(first(rand(1))) * deg2rad(degrees)
        X[:,n] = [-cos(θ) * θ + σ * first(rand(1)), sin(θ) * θ + σ * first(rand(1))]
    end

    for n in N₀+1:N
        θ = s + sqrt(first(rand(1))) * deg2rad(degrees)
        X[:,n] = [cos(θ) * θ + σ * first(rand(1)), -sin(θ) * θ + σ * first(rand(1))]
    end

    X, t
end

function synthetic_circular_data(N, radius₁, radius₂, σ)
    N₀ = floor(Int, N/2)

    X = Matrix{Float64}(undef, 2, N)
    t = [ones(Int, N₀); -ones(Int, N-N₀)]

    for n in 1:N₀
        θ = n / (2 * π)
        X[1,n] = radius₁ * cos(θ) + σ * randn()
        X[2,n] = radius₁ * sin(θ) + σ * randn()
    end

    for n in N₀+1:N
        θ = n / (2 * π)
        X[1,n] = radius₂ * cos(θ) + σ * randn()
        X[2,n] = radius₂ * sin(θ) + σ * randn()
    end

    X, t
end


text:
const DIRS = [:left, :right, :up, :down]
const VALS = vcat(0,[(2).^(1:13)...])
using StatsBase, DataFrames, Flux, Plots, JLD2, FileIO, BSON, Random

function grid_to_feature(grid)
    vec(Flux.onehotbatch(vec(grid),VALS))
end

policy = Flux.Chain(
  grid_to_feature,
  Dense(16*14, 128, relu),
  Dense(128, 4, relu),
  softmax
  )

loss(x, (move, reward)) = -reward*policy(x)[move .== DIRS][1]
opt=ADAM(Flux.params(policy))

@load "res_vec.jld2" res_vec

Random.seed!(0)
for i =1:length(res_vec)
  println(i)
  Flux.train!(loss, res_vec[i:i], opt)
end

polires_vec[800]

loss(res_vec[870]...)

Flux.crossentropy(policy(res_vec[869][1]), Flux.onehot(:up,DIRS))

policy(res_vec[869][1])

Flux.onehot(:up,DIRS)

Flux.crossentropy([false, true, false, false], [1.0, 0.0, 0.0, 0.0] )

Flux.crossentropy([1.0, 0.0,0.0,0.0], [0.0, 1.0, 0.0, 0.0])

-reward*policy(res_vec[1][1])[:left .== DIRS]

ext: jl


text:
import   TexTables: tuplefy

@testset "tuplefy" begin
    vals = [1, 1.0, "1", :1]
    for x in vals
        @test tuplefy(x) == (x,)
    end
    for x in [(val,) for val in vals]
        @test tuplefy(x) == x
    end
    for x in vals, y in vals
        @test tuplefy((x, y)) == (x,y)
    end
end

@testset "summarize" begin
    df = DataFrame(A=[1,2,3],
                   B=[1.0, 2.0, 3.0],
                   C=[true, true, false],
                   D=BitArray([true, true, false]),
                   E=[1, 2, missing],
                   F=[1.0, 2.0, missing],
                   G=["test1", "test2", "test3"])
    t = summarize(df)
    compare_file(34, to_ascii(t))
    compare_file(35, to_tex(t))
end

@testset "twoway tabulate" begin

    # Check a generic one
    Random.seed!(1234)
    iris            = dataset("datasets", "iris")
    iris[:TestVar]  = rand('A':'B', 150)
    t = tabulate(iris, :Species, :TestVar)

    compare_file(36, to_ascii(t))
    compare_file(37, to_tex(t))

end


text:
using Revise
using AdFem

mesh = Mesh(8,8,1/8)
ρ = ones(get_ngauss(mesh))
B = compute_fem_advection_matrix1(2*ρ, 3*ρ, mesh)

ρ = constant(ones(get_ngauss(mesh)))
B2 = compute_fem_advection_matrix1(2*ρ, 3*ρ, mesh)

sess = Session(); init(sess)
maximum(abs.(run(sess, B2) - B))


text:
module Vecchia

  using LinearAlgebra, NearestNeighbors, StaticArrays, SparseArrays
  using FLoops, BangBang, MicroCollections
  using LoopVectorization

  export nll

  include("structstypes.jl")

  include("methods.jl")

  include("utils.jl")

  include("nll.jl")

  include("interfaces.jl")

  include("precision.jl")

  include("rcholesky.jl")

end


text:
## circulant_embedding.jl : Gaussian random field generator using fft; only for uniformly spaced GRFs

"""
    CirculantEmbedding()

Returns a [`GaussianRandomFieldGenerator`](@ref) that uses FFTs to compute samples of the Gaussian random field.

!!! warning
    Circulant embedding can only be applied if the points are specified on a structured grid.

# Optional Arguments for `GaussianRandomField`
- `minnpadding::Integer`: minimum amount of padding.
- `measure::Bool`: optimize the FFT to increase the efficiency of the [`sample`](@ref) method. Default is `true`.
- `primes::Bool`: the size of the minimum circulant embedding of the covariance matrix can be written as a product of small primes (2, 3, 5 and 7). Default is `false`.

# Examples
```jldoctest
julia> cov = CovarianceFunction(2, Matern(.1, 1))
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, CirculantEmbedding(), pts, pts)
Gaussian random field with 2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a circulant embedding

julia> contourf(grf)
[...]

julia> plot_eigenvalues(grf)
[...]

```
!!! note
    With appropriate ordering, the covariance matrix of a Gaussian random field is a (nested block) Toeplitz matrix. This matrix can be embedded into a larger (nested block) circulant matrix, whose eigenvalues can be rapidly computed using FFT. A difficulty here is that although the covariance matrix is positive semi-definite, its circulant extension in general is not. As a remedy, one can add so-called *ghost points* outside the domain of interest using the optional flag `minpadding`.
```jldoctest
julia> cov = CovarianceFunction(2, Matern(.3, 1))
2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0)

julia> pts = range(0, stop=1, length=51)
0.0:0.02:1.0

julia> grf = GaussianRandomField(cov, CirculantEmbedding(), pts, pts)
┌ Warning: 318 negative eigenvalues ≥ -0.5828339433508111 detected, Gaussian random field will be approximated (ignoring all negative eigenvalues); increase padding if possible
└ @ GaussianRandomFields ~/.julia/dev/GaussianRandomFields/src/generators/circulant_embedding.jl:94
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a circulant embedding

julia> grf = GaussianRandomField(cov, CirculantEmbedding(), pts, pts, minpadding=79)
Gaussian random field with 2d Matérn covariance function (λ=0.3, ν=1.0, σ=1.0, p=2.0) on a 51×51 structured grid, using a circulant embedding

```

See also: [`Cholesky`](@ref), [`Spectral`](@ref), [`KarhunenLoeve`](@ref)
"""
struct CirculantEmbedding <: GaussianRandomFieldGenerator end

# whether a covariance function is even in every dimension
# this simplifies the circulant embedding algorithm
Base.iseven(::CovarianceStructure) = false
Base.iseven(::IsotropicCovarianceStructure) = true

function _GaussianRandomField(mean, cov::CovarianceFunction{d}, method::CirculantEmbedding,
                              pts::Vararg{AbstractRange,d};
                              minpadding::Union{Int,Dims{d}} = 0,
                              measure::Bool = true,
                              primes::Bool = false) where {d}
    # normalize points
    normedpts = map(x -> x .- first(x), pts)

    # compute size of minimum circulant embedding
    ρ = cov.cov
    factors = primes ? [2, 3, 5, 7] : [2]
    dims = circulant_minsize.(Ref(ρ), normedpts, minpadding, Ref{Vector{Int}}(factors))

    # compute eigenvalues of the circulant embedding
    Λ = circulant_eigvals(ρ, normedpts, dims)

    # store sqrt of eigenvalues for more efficient sampling
    n₋, λ₋ = 0, 0.0
    M = length(Λ)
    @inbounds for i in eachindex(Λ)
        λ = Λ[i]
        if λ < 0
            # reset negative eigenvalues to zero
            Λ[i] = 0

            # update statistics about negative eigenvalues
            n₋ += 1
            λ < λ₋ && (λ₋ = λ)
        else
            Λ[i] = sqrt(λ / M)
        end
    end
    λ₋ < 0 && @warn "$n₋ negative eigenvalues ≥ $λ₋ detected, Gaussian random field will be
        approximated (ignoring all negative eigenvalues); increase padding if possible"

    # optimize
    copy_of_Λ = Complex.(copy(Λ))
    P = measure ? plan_fft!(copy_of_Λ, flags=FFTW.MEASURE) : plan_fft!(copy_of_Λ)
    data = (Λ, P)

    GaussianRandomField{CirculantEmbedding,typeof(cov),typeof(pts),typeof(mean),typeof(data)}(mean, cov, pts, data)
end

"""
    circulant_minsize(cov::CovarianceStructure, pts::AbstractRange,
                      minpadding::Int, factors::Vector{Int})

Return size of minimum circulant embedding of covariance matrix of points `pts` with
covariance function `cov` and minimum padding `minpadding` that can be written as product of
integers in `factors`.

Typically `factors` is chosen to be `[2]` or `[2, 3, 5, 7]` to speed up FFT computations.
"""
function circulant_minsize(cov::CovarianceStructure, pts::AbstractRange,
                           minpadding::Int, factors::Vector{Int})
    if iseven(cov)
        2 * nextprod(factors, length(pts) + minpadding - 1)
    else
        2 * nextprod(factors, length(pts) + minpadding)
    end
end

"""
    circulant_eigvals(cov::CovarianceStructure, pts::NTuple{N,AbstractRange},
                      dims::Dims{N}) where {N}

Compute eigenvalues of circulant embedding with dimensions `dims` of covariance matrix of
points `pts` with covariance function `cov`.
"""
@generated function circulant_eigvals(cov::CovarianceStructure{T},
                                      pts::NTuple{N,AbstractRange},
                                      dims::Dims{N}) where {T,N}
    quote
        if iseven(cov)
            # compute first row of circulant embedding
            # since the covariance function is even in every dimension,
            # instead of 2n only n+1 values have to be computed in every dimension
            C = Array{$T}(undef, div.(dims, 2) .+ 1)
            x = Array{$T}(undef, $N)

            @nloops $N i C d -> begin
                @inbounds x[d] = extrapolate(pts[d], i_d)
            end begin
                @inbounds (@nref $N C i) = apply(cov, x)
            end

            # compute eigenvalues of circulant embedding
            # since the covariance function is even in every dimension,
            # we can perform a real even FFT and extend the resulting array of eigenvalues
            FFTW.r2r!(C, FFTW.REDFT00)
            Λ = Array{eltype(C)}(undef, dims)
            Imax = CartesianIndex(size(Λ) .+ 2)
            @inbounds for i in CartesianIndices(Λ)
                Λ[i] = C[min(i, Imax - i)]
            end
            Λ
        else
            # compute first row of circulant embedding
            C = zeros($T, dims)
            x = Array{$T}(undef, $N)

            dims2 = dims .+ 2
            mids = div.(dims2, 2)
            @nloops $N i C d -> begin
                @inbounds begin
                    m = mids[d]
                    i_d == m && continue

                    if i_d < m
                        x[d] = extrapolate(pts[d], i_d)
                    else
                        x[d] = -extrapolate(pts[d], dims2[d] - i_d)
                    end
                end
            end begin
                @inbounds (@nref $N C i) = apply(cov, x)
            end

            # compute eigenvalues of circulant embedding
            real.(fft(C))
        end
    end
end

# extrapolate ranges
extrapolate(r::Union{StepRangeLen,LinRange}, i::Integer) = Base.unsafe_getindex(r, i)
extrapolate(r::UnitRange{T}, i::Integer) where T = convert(T, r.start + (i - 1))
extrapolate(r::Base.OneTo{T}, i::Integer) where T = convert(T, i)
extrapolate(r::AbstractRange{T}, i::Integer) where T =
    convert(T, first(r) + (i - 1)*step(r))

# returns the required dimension of the random points
randdim(grf::GaussianRandomField{CirculantEmbedding}) = length(grf.data[1])

# sample function
function _sample(grf::GaussianRandomField{CirculantEmbedding}, xi::AbstractArray{X}) where X
    v = grf.data[1]
    w = Array{complex(promote_type(eltype(v), X))}(undef, size(v))
    z = Array{eltype(grf.cov)}(undef, length.(grf.pts))
    _sample!(w, z, grf, xi)
end

# in-place sample function (#22)
function _sample!(w, z, grf::GaussianRandomField{CirculantEmbedding}, xi)
    v, P = grf.data

    # compute multiplication with square root of circulant embedding via FFT
    w .= complex.(v .* reshape(xi, size(v)))
    mul!(w, P, w)

    # extract realization of random field
    μ, σ = grf.mean, std(grf.cov)
    @inbounds for i in CartesianIndices(z)
        wi = w[i]
        z[i] = μ[i] + σ * (real(wi) + imag(wi))
    end
    z
end

Base.show(io::IO, ::CirculantEmbedding) = print(io, "circulant embedding")


text:
module ConvolutionalOperatorLearning

export CAOL, generatefilters

include("util.jl")          # utility functions
include("core.jl")          # core implementation
include("signatures.jl")    # CAOL function signatures
include("filters.jl")       # filters for initialization

"""
    CAOL(x, λ, h0; p=0, maxiters=2000, tol=1e-13, trace=false)

Learn convolutional analysis operators, i.e., sparsifying filters,
for signals `x` with sparsity regularization `λ` and initial filters `h0`.
Preserve the first `p` filters in `h0`, i.e., consider them handcrafted.

`x` can be either:
+ a vector of training signals each a D-dimensional array
+ a (D+1)-dimnensional array with training sample `i` being slice `x[:,...,:,i]`

`h0` can be either:
+ a vector of filters each a D-dimensional array
+ a tuple `(H0,R)` where each column of the matrix `H0` is a vectorized filter
  and where `R` gives the size/shape of the filter.
The filters must be orthogonal to one another
and each filter must be normalized to have norm `1/sqrt(filter length)`,
i.e., `H0'H0 ≈ (1/size(H0,1))*I`.

When `trace=false`, CAOL returns just the learned filters (in the same form as `h0`).
When `trace=true`, CAOL also returns
+ the iterates (only the learned filters since handcrafted filters do not change)
+ the objective function values (evaluated on only the learned filters)
+ the iterate differences `norm(H[t]-H[t-1])/norm(H[t])` used for the stopping criterion

# Examples
Passing in a vector of ten 100x50 training signals and a vector of 3x3 DCT filters
```julia-repl
julia> x = [randn(100,50) for _ in 1:10];
julia> h0 = generatefilters(:DCT,(3,3));
julia> h = CAOL(x,1e-3,h0,maxiters=10);
```
Passing the ten 100x50 training signals as an array
```julia-repl
julia> X = randn(100,50,10);
julia> h0 = generatefilters(:DCT,(3,3));
julia> h = CAOL(X,1e-3,h0,maxiters=10);
```
Passing the DCT filters in matrix form
```julia-repl
julia> x = [randn(100,50) for _ in 1:10];
julia> H0 = generatefilters(:DCT,(3,3),form=:matrix);
julia> H = CAOL(x,1e-3,(H0,(3,3)),maxiters=10);
```
Getting the trace
```julia-repl
julia> x = [randn(100,50) for _ in 1:10];
julia> H0 = generatefilters(:DCT,(3,3),form=:matrix);
julia> H, Htrace, objtrace, Hdifftrace = CAOL(x,1e-3,(H0,(3,3)),maxiters=10,trace=true);
```
"""
function CAOL end

end # module


text:
using Symbolics
using LinearAlgebra
using Test

@variables t p(t) x y(t)
D = Differential(t)
expr = x * p + y * (1 - p) ~ 0
sol = Symbolics.solve_for(expr, p)
@test sol isa Num
@test isequal(sol, y/(y - x))
a, b, islinear = Symbolics.linear_expansion(expr, p)
@test eltype((a, b)) <: Num
@test isequal((a, b, islinear), (-(x - y), -y, true))
@test isequal(Symbolics.solve_for(x * p ~ 0, p), 0)
@test_throws Any Symbolics.solve_for(1/x + p * p/x ~ 0, p)
@test isequal(Symbolics.solve_for(x * y ~ p, x), p / y)
@test isequal(Symbolics.solve_for(x * -y ~ p, y), -p / x)
@test isequal(Symbolics.solve_for(x * y ~ p, p), x * y)
@test isequal(Symbolics.solve_for(x^2 * y ~ p, y), p / x^2)
@test isequal(Symbolics.solve_for(x^2 * y ~ p, p), x^2 * y)
@test_throws Any Symbolics.solve_for(x^2 * y - sin(p) ~ p, p)
@test Symbolics.solve_for(x^2 * y - sin(p) ~ p, p, check=false) === nothing
@test_throws Any Symbolics.solve_for(t*D(x) ~ y, t)
@test isequal(Symbolics.solve_for(t*D(x) ~ y, D(x)), y/t)
@test isequal(Symbolics.solve_for([t*D(x) ~ y], [D(x)]), [y/t])
@test isequal(Symbolics.solve_for(t*D(x) ~ y, y), t*D(x))
Dx = D(x)
expr = Dx * x + Dx*t - 2//3*x + y*Dx
a, b, islinear = Symbolics.linear_expansion(expr, x)
@test iszero(expand(a * x + b - expr))
@test isequal(Symbolics.solve_for(expr ~ Dx, Dx), (-2//3*x)/(1 - t - x - y))
@test isequal(Symbolics.solve_for(expr ~ Dx, x), (t*Dx + y*Dx - Dx) / ((2//3) - Dx))

exprs = [
 3//2*x + 2y + 10
 7x + 3y - 8
]
xs = [x, y]
A, b, islinear = Symbolics.linear_expansion(exprs, xs)
@test islinear
@test isequal(A, [3//2 2; 7 3])
@test isequal(b, [10; -8])

@variables x y z
eqs = [
        2//1 * x + y - z ~ 2//1
        2//1 + y - z ~ 3//1*x
        2//1 + y - 2z ~ 3//1*z
      ]
@test [2 1 -1; -3 1 -1; 0 1 -5] * Symbolics.solve_for(eqs, [x, y, z]) == [2; -2; -2]
@test isequal(Symbolics.solve_for(2//1*x + y - 2//1*z ~ 9//1*x, 1//1*x), 1//7*y - 2//7*z)


text:
#!/usr/bin/env julia --project
using Test
using ClimateMachine
ClimateMachine.init()
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Grids: polynomialorder
using ClimateMachine.Ocean.HydrostaticBoussinesq
using ClimateMachine.Ocean.ShallowWater
using ClimateMachine.Ocean.SplitExplicit: VerticalIntegralModel
using ClimateMachine.Ocean.OceanProblems

using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.BalanceLaws
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.VTK

using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function test_vertical_integral_model(time; refDat = ())
    mpicomm = MPI.COMM_WORLD
    ArrayType = ClimateMachine.array_type()

    brickrange_2D = (xrange, yrange)
    topl_2D = BrickTopology(
        mpicomm,
        brickrange_2D,
        periodicity = (true, true),
        boundary = ((0, 0), (0, 0)),
    )
    grid_2D = DiscontinuousSpectralElementGrid(
        topl_2D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    brickrange_3D = (xrange, yrange, zrange)
    topl_3D = StackedBrickTopology(
        mpicomm,
        brickrange_3D;
        periodicity = (true, true, false),
        boundary = ((0, 0), (0, 0), (1, 2)),
    )
    grid_3D = DiscontinuousSpectralElementGrid(
        topl_3D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    problem = SimpleBox{FT}(Lˣ, Lʸ, H)

    model_3D = HydrostaticBoussinesqModel{FT}(
        param_set,
        problem;
        cʰ = FT(1),
        αᵀ = FT(0),
        κʰ = FT(0),
        κᶻ = FT(0),
        fₒ = FT(0),
        β = FT(0),
    )

    model_2D = ShallowWaterModel{FT}(
        param_set,
        problem,
        ShallowWater.ConstantViscosity{FT}(model_3D.νʰ),
        nothing;
        c = FT(1),
        fₒ = FT(0),
        β = FT(0),
    )

    integral_bl = VerticalIntegralModel(model_3D)

    integral_model = DGModel(
        integral_bl,
        grid_3D,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    dg_3D = DGModel(
        model_3D,
        grid_3D,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    dg_2D = DGModel(
        model_2D,
        grid_2D,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q_3D = init_ode_state(dg_3D, FT(time); init_on_cpu = true)
    Q_2D = init_ode_state(dg_2D, FT(time); init_on_cpu = true)
    Q_int = integral_model.state_auxiliary

    state_check = ClimateMachine.StateCheck.sccreate(
        [
            (Q_int, "∫u")
            (Q_2D, "U")
        ],
        1;
        prec = 12,
    )

    update_auxiliary_state!(integral_model, integral_bl, Q_3D, time)
    GenericCallbacks.call!(state_check, nothing, nothing, nothing, nothing)

    ClimateMachine.StateCheck.scprintref(state_check)
    if length(refDat) > 0
        @test ClimateMachine.StateCheck.scdocheck(state_check, refDat)
    end

    return nothing
end

#################
# RUN THE TESTS #
#################
FT = Float64

const N = 4
const Nˣ = 5
const Nʸ = 5
const Nᶻ = 8
const Lˣ = 1e6  # m
const Lʸ = 1e6  # m
const H = 400  # m

xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
zrange = range(FT(-H); length = Nᶻ + 1, stop = 0)

const cʰ = 1  # typical of ocean internal-wave speed
const cᶻ = 0

@testset "$(@__FILE__)" begin
    include("../refvals/test_vertical_integral_model_refvals.jl")

    times =
        [0, 86400, 30 * 86400, 365 * 86400, 10 * 365 * 86400, 100 * 365 * 86400]
    for (index, time) in enumerate(times)
        @testset "$(time)" begin
            test_vertical_integral_model(time, refDat = refVals[index])
        end
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# Some GUI is sending a command from the user clicking on a menu
# e.g. "open", "close", "exit", "help"

# Depending on the command, the program does different things.
function process_command(command::String, args...)
    if command == "open"
        # open a file
    elseif command == "close"
        # close current file
    elseif command == "exit"
        # exit program
    elseif command == "help"
        # pops up a help dialog
    else
        error("bug - this should have never happened.")
    end
end

# problem - this function can become very long
# problem - this function must be changed every time a new command is added

struct OpenCommand end

oc1 = OpenCommand()
oc2 = OpenCommand()
oc1 === oc2

# what is Val?
Val(1)
Val(1) |> typeof

#=
julia> Val(1)
Val{1}()

julia> Val(1) |> typeof
Val{1}
=#

#= as defined in Julia, it's nothing but:
struct Val{x}
end
=#

# is it a singleton?
Val(1) === Val(1)
Val(:foo) === Val(:foo)

#=
julia> Val(1) === Val(1)
true

julia> Val(:foo) === Val(:foo)
true
=#

# only bits-type is allowed, otherwise...
#=
julia> Val("julia")
ERROR: TypeError: in Type, in parameter, expected Type, got String
=#

Val(:open)
#=
julia> Val(:open)
Val{:open}()
=#

function process_command(::Val{:open}, filename)
    println("opening file $filename")
end

function process_command(::Val{:close}, filename)
    println("closing file $filename")
end

process_command(Val(:open), "julia.pdf")
#=
julia> process_command(Val(:open), "julia.pdf")
opening file julia.pdf
=#

# look ma, there's no if-then-else statements!
function process_command(command::String, args...)
    process_command(Val(Symbol(command)), args...)
end

process_command("open", "julia.pdf")
process_command("close", "julia.pdf")
#=
julia> process_command("open", "julia.pdf")
opening file julia.pdf

julia> process_command("close", "julia.pdf")
closing file julia.pdf
=#

# Let's build your own parametric type
# A parametric type that represents a specific command
struct Command{T} end

# Constructor function to create a new Command instance from a string
Command(s::AbstractString) = Command{Symbol(s)}()

# Dispatcher function
function process_command(command::String, args...)
    process_command(Command(command), args...)
end

# Actions
function process_command(::Command{:open}, filename)
    println("opening file $filename")
end

function process_command(::Command{:close}, filename)
    println("closing file $filename")
end

# Test
process_command("open", "julia.pdf")
process_command("close", "julia.pdf")

# -----------------------------------------------

ntuple(i -> 2i, 10)

using BenchmarkTools
@btime ntuple(i->2i, 10);
@btime ntuple(i->2i, 11);
#=
julia> ntuple(i -> 2i, 10)
(2, 4, 6, 8, 10, 12, 14, 16, 18, 20)

julia> @btime ntuple(i->2i, 10);
  0.032 ns (0 allocations: 0 bytes)

julia> @btime ntuple(i->2i, 11);
  812.209 ns (4 allocations: 336 bytes)
=#

@btime ntuple(i->2i, Val(10));
@btime ntuple(i->2i, Val(11));
@btime ntuple(i->2i, Val(100));
#= singleton type
julia> @btime ntuple(i->2i, Val(10));
  0.032 ns (0 allocations: 0 bytes)

julia> @btime ntuple(i->2i, Val(11));
  0.031 ns (0 allocations: 0 bytes)

julia> @btime ntuple(i->2i, Val(100));
  17.270 ns (0 allocations: 0 bytes)
=#

# An extreme example from Julia Base
# This is a huge if/then/else block of code for performance reason
function ntuple(f::F, n::Integer) where F
    t = n == 0  ? () :
        n == 1  ? (f(1),) :
        n == 2  ? (f(1), f(2)) :
        n == 3  ? (f(1), f(2), f(3)) :
        n == 4  ? (f(1), f(2), f(3), f(4)) :
        n == 5  ? (f(1), f(2), f(3), f(4), f(5)) :
        n == 6  ? (f(1), f(2), f(3), f(4), f(5), f(6)) :
        n == 7  ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7)) :
        n == 8  ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8)) :
        n == 9  ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8), f(9)) :
        n == 10 ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8), f(9), f(10)) :
        _ntuple(f, n)
    return t
end

function _ntuple(f, n)
    @_noinline_meta
    (n >= 0) || throw(ArgumentError(string("tuple length should be ≥ 0, got ", n)))
    ([f(i) for i = 1:n]...,)
end

# Using singleton type dynamic dispatch
# inferrable ntuple (enough for bootstrapping)
ntuple(f, ::Val{0}) = ()
ntuple(f, ::Val{1}) = (@_inline_meta; (f(1),))
ntuple(f, ::Val{2}) = (@_inline_meta; (f(1), f(2)))
ntuple(f, ::Val{3}) = (@_inline_meta; (f(1), f(2), f(3)))

@inline function ntuple(f::F, ::Val{N}) where {F,N}
    N::Int
    (N >= 0) || throw(ArgumentError(string("tuple length should be ≥ 0, got ", N)))
    if @generated
        quote
            @nexprs $N i -> t_i = f(i)
            @ncall $N tuple t
        end
    else
        Tuple(f(i) for i = 1:N)
    end
end

# Optional: let's understand what's going on here

using Base.Cartesian
@nexprs 3 i -> t_i = i

# It's the same thing as:
# t_1 = 1
# t_2 = 2
# t_3 = 3

#=
julia> @macroexpand(@nexprs 3 i -> t_i = i)
quote
    t_1 = 1
    t_2 = 2
    t_3 = 3
end
=#

#=
julia> @macroexpand(@ncall 3 tuple t)
:(tuple(t_1, t_2, t_3))
=#

# of course, we cannot hard code the number 3.
# but the macro can only work wth an integer not a variable...
# that's why the implementation wraps it in a @generated block
# and interpolate the value of N into the @ncall and @nexprs calls.

#=
julia> @macroexpand(@ncall N tuple t)
ERROR: LoadError: MethodError: no method matching @ncall(::LineNumberNode, ::Module, ::Symbol, ::Symbol, ::Symbol)
Closest candidates are:
  @ncall(::LineNumberNode, ::Module, ::Int64, ::Any, ::Any...) at cartesian.jl:105
=#


text:
using ImageFiltering, ImageCore
using Test
import StaticArrays
using Random

# Ambiguity test
if Base.VERSION >= v"1.6.0-DEV.1005"   # julia #37616
    @test isempty(detect_ambiguities(ImageFiltering, Kernel, KernelFactors))
else
    # Because StaticArrays may have ambiguities with Base, we have to "subtract" these
    aif = detect_ambiguities(ImageFiltering, Kernel, KernelFactors, Base)
    asa = detect_ambiguities(StaticArrays, Base)
    @test isempty(setdiff(aif, asa))
end

function typestring(::Type{T}) where T   # from https://github.com/JuliaImages/ImageCore.jl/pull/133
    buf = IOBuffer()
    show(buf, T)
    String(take!(buf))
end

include("compat.jl")
include("border.jl")
include("nd.jl")
include("2d.jl")
include("triggs.jl")
include("cascade.jl")
include("specialty.jl")
include("gradient.jl")
include("mapwindow.jl")
include("extrema.jl")
include("basic.jl")
include("gabor.jl")

nothing


text:
@testset "GenusRep" begin

  #
  # Some definite examples
  #

  # Lattice 214 from the database: `a(P) != P` in `_neighbours`

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2 - x - 1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2 - a + 3
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 3, 3, [1, 0, 0, 0, 1, 0, 0, 0, 1])
  gene = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 0, 0]), map(E, [0, 1, -1]), map(E, [0, -1//2*a*b - 1//2*a + 3//2, 0])]

  L = hermitian_lattice(E, gene, gram = D)
  gens, def, P0 = @inferred Hecke.genus_generators(L)
  @test isempty(gens)
  @test def

  gen_rep = @inferred genus_representatives(L)
  @test length(gen_rep) == 2
  @test L in gen_rep
  @test !isisometric(gen_rep[1], gen_rep[2])[1]

  # Lattice 324 from the database: `special == true` in `_neighbours`

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^2 - x - 1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2 - a + 3
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [2, 1, 0, 0]), map(E, [(a + 3//2)*b - 13//2*a + 2, 0, (1//2*a + 1//2)*b + 3*a - 33//2, 0]), map(E, [0, 0, 2, 1]), map(E, [(-92011//2*a - 62822)*b - 15049//2*a + 178467//2, 0, (-126700*a + 60200)*b + 160300*a - 109900, 0])]

  L = hermitian_lattice(E, gens, gram = D)
  gens, def, P0 = @inferred Hecke.genus_generators(L)
  @test isempty(gens)
  a = involution(L)
  @test a(P0) == P0
  ok, scale = ismodular(L,P0)
  if scale != 0 && isramified(base_ring(L), minimum(P0))
    special = isodd(scale)
  end
  @test special

  @test genus_representatives(L) == [L]

  # Lattice 441 from the database: `special == false` in `_neighbours`; calls `_all_row_span`

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x - 1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2 + 11
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 4, 4, [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [5//2*b + 121//2, 1//2*b + 55//2, 0, 0]), map(E, [-1, -3, -1, 0]), map(E, [987//2*b - 2407//2, 475//2*b - 679//2, 0, 2*b + 62]), map(E, [3906*b - 16305, 2074*b - 5477, 0, 70*b + 595])]

  L = hermitian_lattice(E, gens, gram = D)
  gens, def, P0 = @inferred Hecke.genus_generators(L)
  a = involution(L)
  @test a(P0) == P0
  ok, scale = ismodular(L,P0)
  @test scale == 0
  gen_rep = @inferred representatives(genus(L))
  @test any(LL -> isisometric(LL,L)[1], gen_rep)
  @test !all(LL -> isisometric(LL,L)[1], gen_rep)

  #
  # An indefinite example (see [Kir19, Page 9])
  #

  K, a = rationals_as_number_field()
  Kt, t = PolynomialRing(K, "t")
  g = t^2 + 17
  E, b = NumberField(g, "b", cached = false)
  D = matrix(E, 2, 2, [102, b, -b, 0])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 0]), map(E, [b, 0]), map(E, [0, 1]), map(E, [0, b])]

  L = hermitian_lattice(E, gens, gram = D)
  gens, def, P0 = @inferred Hecke.genus_generators(L)
  @test !def
  @test length(gens) == 1
  P, e = gens[1]
  @test minimum(P) == 3*maximal_order(K)
  @test e == 4
  @test P0 == 1*maximal_order(E)

  gen_rep = @inferred genus_representatives(L)
  @test length(gen_rep) == 4
  @test L in gen_rep
  a = involution(L)
  @test a(P) != P
  PB = [pseudo_basis(LL) for LL in gen_rep]
  @test all(i -> PB[i][1][2] == P^(i-1) && PB[i][2][2]^-1 == a(P)^(i-1), 1:length(PB))

  #
  # Other indefinite examples
  #

  Qx, x = FlintQQ["x"]
  K, a = NumberField(x - 1, "a")
  Kt, t = K["t"]
  E, b = NumberField(t^2 + 1, "b")
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  G = genus(HermLat, E, p, [(0, 3, -1, 0)])
  L = @inferred representative(G)
  @test length(@inferred Hecke.genus_representatives(L)) == 1

  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^6 + x^5 + x^4 + x^3 + x^2 + x + 1
  K, a = NumberField(f, "a", cached = false)
  Kt, t = PolynomialRing(K, "t")
  g = t^2 + 47
  E, z_7 = NumberField(g, "b", cached = false)
  b = z_7
  D = matrix(E, 6, 6, [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1])
  gens = Vector{Hecke.NfRelElem{nf_elem}}[map(E, [1, 0, 0, 0, 0, 0]), map(E, [z_7, 0, 0, 0, 0, 0]), map(E, [z_7^2, 0, 0, 0, 0, 0]), map(E, [z_7^3, 0, 0, 0, 0, 0]), map(E, [z_7^4, 0, 0, 0, 0, 0]), map(E, [z_7^5, 0, 0, 0, 0, 0]), map(E, [1//2*b + 1//2, 0, 0, 0, 0, 0]), map(E, [1//2*z_7*b + 1//2*z_7, 0, 0, 0, 0, 0]), map(E, [1//2*z_7^2*b + 1//2*z_7^2, 0, 0, 0, 0, 0]), map(E, [1//2*z_7^3*b + 1//2*z_7^3, 0, 0, 0, 0, 0]), map(E, [1//2*z_7^4*b + 1//2*z_7^4, 0, 0, 0, 0, 0]), map(E, [1//2*z_7^5*b + 1//2*z_7^5, 0, 0, 0, 0, 0]), map(E, [0, 1, 0, 0, 0, 0]), map(E, [0, z_7, 0, 0, 0, 0]), map(E, [0, z_7^2, 0, 0, 0, 0]), map(E, [0, z_7^3, 0, 0, 0, 0]), map(E, [0, z_7^4, 0, 0, 0, 0]), map(E, [0, z_7^5, 0, 0, 0, 0]), map(E, [0, 1//2*b + 1//2, 0, 0, 0, 0]), map(E, [0, 1//2*z_7*b + 1//2*z_7, 0, 0, 0, 0]), map(E, [0, 1//2*z_7^2*b + 1//2*z_7^2, 0, 0, 0, 0]), map(E, [0, 1//2*z_7^3*b + 1//2*z_7^3, 0, 0, 0, 0]), map(E, [0, 1//2*z_7^4*b + 1//2*z_7^4, 0, 0, 0, 0]), map(E, [0, 1//2*z_7^5*b + 1//2*z_7^5, 0, 0, 0, 0]), map(E, [0, 0, 1, 0, 0, 0]), map(E, [0, 0, z_7, 0, 0, 0]), map(E, [0, 0, z_7^2, 0, 0, 0]), map(E, [0, 0, z_7^3, 0, 0, 0]), map(E, [0, 0, z_7^4, 0, 0, 0]), map(E, [0, 0, z_7^5, 0, 0, 0]), map(E, [0, 0, 1//2*b + 1//2, 0, 0, 0]), map(E, [0, 0, 1//2*z_7*b + 1//2*z_7, 0, 0, 0]), map(E, [0, 0, 1//2*z_7^2*b + 1//2*z_7^2, 0, 0, 0]), map(E, [0, 0, 1//2*z_7^3*b + 1//2*z_7^3, 0, 0, 0]), map(E, [0, 0, 1//2*z_7^4*b + 1//2*z_7^4, 0, 0, 0]), map(E, [0, 0, 1//2*z_7^5*b + 1//2*z_7^5, 0, 0, 0]), map(E, [0, 0, 0, 1, 0, 0]), map(E, [0, 0, 0, z_7, 0, 0]), map(E, [0, 0, 0, z_7^2, 0, 0]), map(E, [0, 0, 0, z_7^3, 0, 0]), map(E, [0, 0, 0, z_7^4, 0, 0]), map(E, [0, 0, 0, z_7^5, 0, 0]), map(E, [0, 0, 0, 1//2*b + 1//2, 0, 0]), map(E, [0, 0, 0, 1//2*z_7*b + 1//2*z_7, 0, 0]), map(E, [0, 0, 0, 1//2*z_7^2*b + 1//2*z_7^2, 0, 0]), map(E, [0, 0, 0, 1//2*z_7^3*b + 1//2*z_7^3, 0, 0]), map(E, [0, 0, 0, 1//2*z_7^4*b + 1//2*z_7^4, 0, 0]), map(E, [0, 0, 0, 1//2*z_7^5*b + 1//2*z_7^5, 0, 0]), map(E, [0, 0, 0, 0, 1, 0]), map(E, [0, 0, 0, 0, z_7, 0]), map(E, [0, 0, 0, 0, z_7^2, 0]), map(E, [0, 0, 0, 0, z_7^3, 0]), map(E, [0, 0, 0, 0, z_7^4, 0]), map(E, [0, 0, 0, 0, z_7^5, 0]), map(E, [0, 0, 0, 0, 1//2*b + 1//2, 0]), map(E, [0, 0, 0, 0, 1//2*z_7*b + 1//2*z_7, 0]), map(E, [0, 0, 0, 0, 1//2*z_7^2*b + 1//2*z_7^2, 0]), map(E, [0, 0, 0, 0, 1//2*z_7^3*b + 1//2*z_7^3, 0]), map(E, [0, 0, 0, 0, 1//2*z_7^4*b + 1//2*z_7^4, 0]), map(E, [0, 0, 0, 0, 1//2*z_7^5*b + 1//2*z_7^5, 0]), map(E, [0, 0, 0, 0, 0, 1]), map(E, [0, 0, 0, 0, 0, z_7]), map(E, [0, 0, 0, 0, 0, z_7^2]), map(E, [0, 0, 0, 0, 0, z_7^3]), map(E, [0, 0, 0, 0, 0, z_7^4]), map(E, [0, 0, 0, 0, 0, z_7^5]), map(E, [0, 0, 0, 0, 0, 1//2*b + 1//2]), map(E, [0, 0, 0, 0, 0, 1//2*z_7*b + 1//2*z_7]), map(E, [0, 0, 0, 0, 0, 1//2*z_7^2*b + 1//2*z_7^2]), map(E, [0, 0, 0, 0, 0, 1//2*z_7^3*b + 1//2*z_7^3]), map(E, [0, 0, 0, 0, 0, 1//2*z_7^4*b + 1//2*z_7^4]), map(E, [0, 0, 0, 0, 0, 1//2*z_7^5*b + 1//2*z_7^5])]
  L = hermitian_lattice(E, gens, gram = D)
  @test length(Hecke.genus_representatives(L)) == 15

end


text:
module Example

export func

"""
    func(x)

Returns double the number `x` plus `1`.

Here's an equation:

``\\frac{n!}{k!(n - k)!} = \\binom{n}{k}``

"""
func(x) = 2x + 1

end


text:
using DensityEstimationML
using Base.Test

function testset_from_file(filename)
    @testset "$filename" begin
        include(filename)
    end
end

testset_from_file.("datagen.jl")


text:
xticks = (1:2, ["Load Shedding", "Load Shedding"])
p = plot(xticks = xticks, ylabel = "Energy [MWh]", aspect_ratio = 1.77777, xlims=(0.8, 1.2))
boxplot!(load_shedding_uc./12, color = :gray, label = "UC")
boxplot!(load_shedding_suc./12, color = :white, label = "SUC")
Plots.savefig("first_measures.pdf")

xticks = (1:4, ["Fuel Cost", "Fuel Cost"])
p = plot(xticks = xticks, ylabel = "Fuel Cost [1000-\$]", aspect_ratio =3)
boxplot!(fuel_cost_uc./12e3, color = :gray, label = "UC")
boxplot!(fuel_cost_suc./12e3, color = :white, label = "SUC")
Plots.savefig("Fuel.pdf")

boxplot!(curtailment_uc./12, color = :gray, label = "UC")
boxplot!(curtailment_suc./12, color = :white, label = "SUC")


text:
using Test, Random
using Knet.Ops20: dropout
using AutoGrad: gradcheck
using CUDA: CUDA, functional
using Knet.KnetArrays: KnetArray

@testset "dropout" begin
    dropout1(x,p)=dropout(x,p;seed=1,drop=true)
    a = rand(Random.MersenneTwister(2),100,100)
    @test gradcheck(dropout1,a,0.5; args=1)
    if CUDA.functional()
        k = KnetArray(a)
        @test gradcheck(dropout1,k,0.5; args=1)
        # This fails because seeds work differently on cpu vs gpu
        # @test isapprox(dropout1(k,0.5),dropout1(a,0.5))
        @test isapprox(sum(abs2,dropout1(k,0.5)), sum(abs2,dropout1(a,0.5)), rtol=0.1)
    end
end


text:
abstract type AbstractFiniteDifferences{T,I<:Integer} <: Basis{T} end
const RestrictedFiniteDifferences{T,B<:AbstractFiniteDifferences{T}} =
    RestrictedQuasiArray{T,2,B}
const FiniteDifferencesOrRestricted = BasisOrRestricted{<:AbstractFiniteDifferences}

eltype(::AbstractFiniteDifferences{T}) where T = T

axes(B::AbstractFiniteDifferences{T}) where T = (Inclusion(leftendpoint(B)..rightendpoint(B)), Base.OneTo(length(locs(B))))
size(B::AbstractFiniteDifferences) = (ℵ₁, length(locs(B)))
local_step(B::AbstractFiniteDifferences, _) = step(B)
weight(B::AbstractFiniteDifferences{T}, _) where T = one(T)

locs(B::RestrictedFiniteDifferences) = locs(parent(B))[indices(B,2)]

==(A::AbstractFiniteDifferences,B::AbstractFiniteDifferences) = locs(A) == locs(B)

tent(x, x₀::T, δx::T) where T =
    one(T) - abs(clamp((x-x₀)/δx, -one(T), one(T)))

getindex(B::AbstractFiniteDifferences, x::Real, i::Integer) =
    weight(B,i)*tent(x, locs(B)[i], local_step(B, i))

step(B::RestrictedFiniteDifferences) = step(parent(B))
local_step(B::RestrictedFiniteDifferences, i) = local_step(parent(B), j+first(indices(B,2))-1)

"""
    within_interval(x, interval)

Return the indices of the elements of `x` that lie within the given
closed `interval`.
"""
function within_interval(x::AbstractRange, interval::ClosedInterval)
    a = leftendpoint(interval)
    b = rightendpoint(interval)
    δx = step(x)
    max(ceil(Int, (a-x[1])/δx),1):min(floor(Int, (b-x[1])/δx)+1,length(x))
end

function getindex(B::AbstractFiniteDifferences{T}, x::AbstractRange, sel::AbstractVector) where T
    l = locs(B)
    χ = spzeros(T, length(x), length(sel))
    o = sel[1] - 1
    for j in sel
        x₀ = l[j]
        J = j - o
        δx = local_step(B, j)
        w = weight(B, j)
        for i ∈ within_interval(x, x₀ - δx..x₀ + δx)
            χ[i,J] = w*tent(x[i], x₀, δx)
        end
    end
    χ
end

getindex(B::FiniteDifferencesOrRestricted, x::AbstractRange, ::Colon) =
    getindex(parent(B), x, indices(B,2))

getindex(B::RestrictedFiniteDifferences, x::AbstractRange, sel::AbstractVector) =
    getindex(parent(B), x, indices(B,2)[sel])

# * Types

const FDArray{T,N,B<:FiniteDifferencesOrRestricted} = MulQuasiArray{T,N,<:Tuple{<:B,<:AbstractArray{T,N}}}
const FDVector{T,B} = FDArray{T,1,B}
const FDMatrix{T,B} = FDArray{T,2,B}
const FDVecOrMat{T,B} = Union{FDVector{T,B},FDMatrix{T,B}}

const AdjointFDArray{T,N,B<:FiniteDifferencesOrRestricted} = MulQuasiArray{T,<:Any,<:Tuple{<:Adjoint{T,<:AbstractArray{T,N}},
                                                                                <:QuasiAdjoint{T,<:B}}}
const AdjointFDVector{T,B} = AdjointFDArray{T,1,B}
const AdjointFDMatrix{T,B} = AdjointFDArray{T,2,B}
const AdjointFDVecOrMat{T,B} = Union{AdjointFDVector{T,B},AdjointFDMatrix{T,B}}

# This is an operator on the form O = B*M*B⁻¹, where M is a matrix
# acting on the expansion coefficients. When applied to e.g. a
# FDVector, v = B*c, we get O*v = B*M*B⁻¹*B*c = B*M*c. Acting on an
# AdjointFDVector is also trivial: v'O' = (B*c)'*(B*M*B⁻¹)' =
# B'c'M'. It is less clear if it is possible to form Hermitian
# operator that does _not_ need be transposed before application to
# adjoint vectors.
const FDOperator{T,B<:FiniteDifferencesOrRestricted,M<:AbstractMatrix} = MulQuasiArray{T,2,<:Tuple{<:B,<:M,
                                                                                <:PInvQuasiMatrix{<:Any,<:Tuple{<:B}}}}

const FDMatrixElement{T,B<:FiniteDifferencesOrRestricted,M<:AbstractMatrix,V<:AbstractVector} =
    MulQuasiArray{T,0,<:Mul{<:Any,
                            <:Tuple{<:Adjoint{<:Any,<:V},<:QuasiAdjoint{<:Any,<:B},
                                    <:B,<:M,<:QuasiAdjoint{<:Any,<:B},
                                    <:B,<:V}}}

const FDInnerProduct{T,U,B<:FiniteDifferencesOrRestricted{U},V1<:AbstractVector{T},V2<:AbstractVector{T}} =
    Mul{<:Any, <:Tuple{<:Adjoint{<:Any,<:V1},<:QuasiAdjoint{<:Any,<:B},<:B,<:V2}}

const LazyFDInnerProduct{FD<:FiniteDifferencesOrRestricted} = Mul{<:Any,<:Tuple{
    <:Mul{<:Any, <:Tuple{
        <:Adjoint{<:Any,<:AbstractVector},
        <:QuasiAdjoint{<:Any,<:FD}}},
    <:Mul{<:Any, <:Tuple{
        <:FD,
        <:AbstractVector}}}}

# * Mass matrix

@simplify function *(Ac::QuasiAdjoint{<:Any,<:AbstractFiniteDifferences}, B::AbstractFiniteDifferences)
    A = parent(Ac)
    A == B || throw(ArgumentError("Cannot multiply functions on different grids"))
    I*step(B)
end

# * Basis inverses

@simplify function *(A⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:AbstractFiniteDifferences}}, B::AbstractFiniteDifferences)
    A = parent(A⁻¹)
    A == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    one(eltype(A))*I
end

@simplify function *(A::AbstractFiniteDifferences, B⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:AbstractFiniteDifferences}})
    B = parent(B⁻¹)
    A == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    one(eltype(A))*I
end

@simplify function *(A⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:AbstractFiniteDifferences}}, v::FDArray)
    A = parent(A⁻¹)
    B,c = v.args
    A == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    c
end

@simplify function *(v::AdjointFDArray, B⁻¹::PInvQuasiMatrix{<:Any,<:Tuple{<:AbstractFiniteDifferences}})
    c,Ac = v.args
    B = parent(B⁻¹)
    parent(Ac) == B || throw(ArgumentError("Cannot multiply basis with inverse of other basis"))
    c
end

# * Various finite differences
# ** Cartesian finite differences
struct FiniteDifferences{T,I} <: AbstractFiniteDifferences{T,I}
    j::UnitRange{I}
    Δx::T
end
FiniteDifferences(j::UnitRange{I}, Δx::T) where {T<:Real,I<:Integer} =
    FiniteDifferences{T,I}(j, Δx)

FiniteDifferences(n::I, Δx::T) where {T<:Real,I<:Integer} =
    FiniteDifferences(1:n, Δx)

locs(B::FiniteDifferences) = B.j*B.Δx
step(B::FiniteDifferences{T}) where {T} = B.Δx

IntervalSets.leftendpoint(B::FiniteDifferences) = (B.j[1]-1)*step(B)
IntervalSets.rightendpoint(B::FiniteDifferences) = (B.j[end]+1)*step(B)

show(io::IO, B::FiniteDifferences{T}) where {T} =
    write(io, "Finite differences basis {$T} on $(axes(B,1).domain) with $(size(B,2)) points spaced by Δx = $(B.Δx)")

# ** Radial finite differences
# Specialized finite differences for the case where there is
# Dirichlet0 boundary condition at r = 0.
struct RadialDifferences{T,I} <: AbstractFiniteDifferences{T,I}
    j::Base.OneTo{I}
    ρ::T
    Z::T
    δβ₁::T # Correction used for bare Coulomb potentials, Eq. (22) Schafer2009

    RadialDifferences(n::I, ρ::T, Z=one(T),
                      δβ₁=Z*ρ/8 * (one(T) + Z*ρ)) where {I<:Integer, T} =
                          new{T,I}(Base.OneTo(n), ρ, convert(T,Z), convert(T,δβ₁))
end

"""
    RadialDifferences(rₘₐₓ::T, n::I, args...)

Convenience constructor for [`RadialDifferences`](@ref) covering the
open interval `(0,rₘₐₓ)` with `n` grid points.
"""
RadialDifferences(rₘₐₓ::T, n::I, args...) where {I<:Integer, T} =
    RadialDifferences(n, rₘₐₓ/(n+one(T)/2), args...)

locs(B::RadialDifferences{T}) where T = (B.j .- one(T)/2)*B.ρ
step(B::RadialDifferences{T}) where {T} = B.ρ

IntervalSets.leftendpoint(B::RadialDifferences{T}) where T = zero(T)
IntervalSets.rightendpoint(B::RadialDifferences{T}) where T = (B.j[end]+one(T)/2)*step(B)

==(A::RadialDifferences,B::RadialDifferences) = locs(A) == locs(B) && A.Z == B.Z && A.δβ₁ == B.δβ₁

show(io::IO, B::RadialDifferences{T}) where T =
    write(io, "Radial finite differences basis {$T} on $(axes(B,1).domain) with $(size(B,2)) points spaced by ρ = $(B.ρ)")

# ** Staggered finite differences

# Generalization of RadialDifferences to non-uniform grid, c.f.
#
# - Krause, J. L., & Schafer, K. J. (1999). Control of THz Emission from
#   Stark Wave Packets. The Journal of Physical Chemistry A, 103(49),
#   10118–10125. http://dx.doi.org/10.1021/jp992144

local_step(r::AbstractRange, _) = step(r)
function local_step(r::AbstractVector, j)
    if j == 1
        (r[2] - r[1])
    elseif j == length(r)
        r[j] - r[j-1]
    else
        (r[j+1]-r[j-1])/2
    end
end

schafer_corner_fix(ρ, Z) = 1/ρ^2*Z*ρ/8*(1 + Z*ρ)

function log_lin_grid!(r, ρₘᵢₙ, ρₘₐₓ, α, js)
    δρ = ρₘₐₓ-ρₘᵢₙ
    for j = js
        r[j] = r[j-1] + ρₘᵢₙ + (1-exp(-α*r[j-1]))*δρ
    end
end

function log_lin_grid(ρₘᵢₙ, ρₘₐₓ, α, n::Integer)
    r = zeros(n)
    r[1] = ρₘᵢₙ/2
    log_lin_grid!(r, ρₘᵢₙ, ρₘₐₓ, α, 2:n)
    r
end

function log_lin_grid(ρₘᵢₙ, ρₘₐₓ, α, rₘₐₓ)
    # Guess n
    n = ceil(Int, rₘₐₓ/ρₘₐₓ)
    r = zeros(n)
    r[1] = ρₘᵢₙ/2
    log_lin_grid!(r, ρₘᵢₙ, ρₘₐₓ, α, 2:n)
    nprev = n
    while r[end] < rₘₐₓ
        n *= 2
        resize!(r, n)
        log_lin_grid!(r, ρₘᵢₙ, ρₘₐₓ, α, nprev+1:n)
        nprev = n
    end
    j = findlast(<(rₘₐₓ), r)
    r[1:j+1]
end

struct StaggeredFiniteDifferences{T,V<:AbstractVector{T},A,B} <: AbstractFiniteDifferences{T,Int}
    r::V
    Z::T
    δβ₁::T # Correction used for bare Coulomb potentials, Eq. (22) Schafer2009
    α::A
    β::B
    function StaggeredFiniteDifferences(r::V, Z=one(T),
                                        δβ₁=schafer_corner_fix(local_step(r,1), Z)) where {T,V<:AbstractVector{T}}
        N = length(r)
        α = zeros(T, N-1)
        β = zeros(T, N)

        r̃ = vcat(r, 2r[end]-r[end-1])
        a = 2r[1]-r[2]
        b,c,d = r[1],r[2],r[3]

        β[1] = δβ₁

        for j = 1:N
            if j < N
                d = r̃[j+2]
                α[j] = 2/((c-b)*√((d-b)*(c-a)))*((b+c)/2)^2/(b*c)
            end

            β[j] += 1/(c-a)*(1/(c-b)*((c+b)/2b)^2 +
                             1/(b-a)*((b+a)/2b)^2)

            a,b,c = b,c,d
        end

        new{T,V,typeof(α),typeof(β)}(r, T(Z), T(δβ₁), α, β)
    end

    # Constructor of uniform grid
    StaggeredFiniteDifferences(n::Integer, ρ, args...) =
        StaggeredFiniteDifferences(((1:n) .- 0.5)*ρ, args...)

    # Constructor of log-linear grid
    StaggeredFiniteDifferences(ρₘᵢₙ, ρₘₐₓ, α, n::Integer, args...) =
        StaggeredFiniteDifferences(log_lin_grid(ρₘᵢₙ, ρₘₐₓ, α, n::Integer), args...)

    StaggeredFiniteDifferences(ρₘᵢₙ::T, ρₘₐₓ::T, α::T, rₘₐₓ::T, args...) where {T<:Real} =
        StaggeredFiniteDifferences(log_lin_grid(ρₘᵢₙ, ρₘₐₓ, α, rₘₐₓ), args...)
end

locs(B::StaggeredFiniteDifferences) = B.r
# The (locally varying) step is baked into the expansion coefficients
step(B::StaggeredFiniteDifferences{T}) where T = one(T)
local_step(B::StaggeredFiniteDifferences, j) = local_step(B.r, j)
weight(B::StaggeredFiniteDifferences, j) = 1/√(local_step(B, j))

IntervalSets.leftendpoint(B::StaggeredFiniteDifferences{T}) where T = zero(T)
IntervalSets.rightendpoint(B::StaggeredFiniteDifferences{T}) where T = 2B.r[end]-B.r[end-1]

==(A::StaggeredFiniteDifferences,B::StaggeredFiniteDifferences) = locs(A) == locs(B) && A.Z == B.Z && A.δβ₁ == B.δβ₁

show(io::IO, B::StaggeredFiniteDifferences{T}) where T =
    write(io, "Staggered finite differences basis {$T} on $(axes(B,1).domain) with $(size(B,2)) points")

# ** Numerov finite differences
#=
This is an implementation of finite difference scheme described in

- Muller, H. G. (1999). An Efficient Propagation Scheme for the
  Time-Dependent Schrödinger equation in the Velocity Gauge. Laser
  Physics, 9(1), 138–148.

where the first derivative is approximated as

\[\partial f =
\left(1+\frac{h^2}{6}\Delta_2\right)^{-1}
\Delta_1 f
\equiv
M_1^{-1}\tilde{\Delta}_1 f,\]
where
\[M_1 \equiv
\frac{1}{6}
\begin{bmatrix}
4+\lambda' & 1 &\\
1 & 4 & 1\\
& 1 & 4 & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]
and
\[\tilde{\Delta}_1 \equiv
\frac{1}{2h}
\begin{bmatrix}
\lambda & 1 &\\
-1 &  & 1\\
& -1 &  & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]

where \(\lambda=\lambda'=\sqrt{3}-2\) for problems with a singularity
at the boundary \(r=0\) and zero otherwise; and the second derivative
as

\[\partial^2 f =
\left(1+\frac{h^2}{12}\Delta_2\right)^{-1}
\Delta_2 f
\equiv
-2M_2^{-1}\Delta_2 f,\]
where
\[M_2 \equiv
-\frac{1}{6}
\begin{bmatrix}
10-2\delta\beta_1 & 1 &\\
1 & 10 & 1\\
& 1 & 10 & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]
and
\[\Delta_2 \equiv
\frac{1}{h^2}
\begin{bmatrix}
-2(1+\delta\beta_1) & 1 &\\
1 & -2 & 1\\
& 1 & -2 & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]

where, again, \(\delta\beta_1 = -Zh[12-10Zh]^{-1}\) is a correction
introduced for problems singular at the origin.

=#

struct NumerovFiniteDifferences{T,I} <: AbstractFiniteDifferences{T,I}
    j::UnitRange{I}
    Δx::T
    # Used for radial problems with a singularity at r = 0.
    λ::T
    δβ₁::T
end

function NumerovFiniteDifferences(j::UnitRange{I}, Δx::T, singular_origin::Bool=false, Z=zero(T)) where {T<:Real,I<:Integer}
    λ,δβ₁ = if singular_origin
        first(j) == 1 ||
            throw(ArgumentError("Singular origin correction only valid when grid starts at Δx (i.e. `j[1] == 1`)"))
        # Eqs. (20,17), Muller (1999)
        (√3 - 2),(-Z*Δx/(12 - 10Z*Δx))
    else
        zero(T),zero(T)
    end
    NumerovFiniteDifferences{T,I}(j, Δx, λ, δβ₁)
end

NumerovFiniteDifferences(n::I, Δx::T, args...) where {T<:Real,I<:Integer} =
    NumerovFiniteDifferences(1:n, Δx, args...)

locs(B::NumerovFiniteDifferences) = B.j*B.Δx
step(B::NumerovFiniteDifferences{T}) where {T} = B.Δx

IntervalSets.leftendpoint(B::NumerovFiniteDifferences) = (B.j[1]-1)*step(B)
IntervalSets.rightendpoint(B::NumerovFiniteDifferences) = (B.j[end]+1)*step(B)

show(io::IO, B::NumerovFiniteDifferences{T}) where {T} =
    write(io, "Numerov finite differences basis {$T} on $(axes(B,1).domain) with $(size(B,2)) points spaced by Δx = $(B.Δx)")

mutable struct NumerovDerivative{T,Tri,Mat,MatFact} <: AbstractMatrix{T}
    Δ::Tri
    M::Mat
    M⁻¹::MatFact
    c::T
end
NumerovDerivative(Δ::Tri, M::Mat, c::T) where {T,Tri,Mat} =
    NumerovDerivative(Δ, M, factorize(M), c)

Base.show(io::IO, ∂::ND) where {ND<:NumerovDerivative} =
    write(io, "$(size(∂,1))×$(size(∂,2)) $(ND)")

Base.show(io::IO, ::MIME"text/plain", ∂::ND) where {ND<:NumerovDerivative} =
    show(io, ∂)

Base.size(∂::ND, args...) where {ND<:NumerovDerivative} = size(∂.Δ, args...)
Base.eltype(∂::ND) where {ND<:NumerovDerivative} = eltype(∂.Δ)
Base.axes(∂::ND, args...) where {ND<:NumerovDerivative} = axes(∂.Δ, args...)

Base.:(*)(∂::ND,a::T) where {T<:Number,ND<:NumerovDerivative} =
    NumerovDerivative(∂.Δ, ∂.M, ∂.M⁻¹, ∂.c*a)

Base.:(*)(a::T,∂::ND) where {T<:Number,ND<:NumerovDerivative} =
    ∂ * a

Base.:(/)(∂::ND,a::T) where {T<:Number,ND<:NumerovDerivative} =
    ∂ * inv(a)

function LinearAlgebra.mul!(y::Y, ∂::ND, x::X) where {Y<:AbstractVector,
                                                      ND<:NumerovDerivative,
                                                      X<:AbstractVector}
    mul!(y, ∂.Δ, x)
    ldiv!(∂.M⁻¹, y)
    lmul!(∂.c, y)
    y
end

function Base.copyto!(y::Y, ∂::Mul{<:Any,Tuple{<:NumerovDerivative, X}}) where {X<:AbstractVector,Y<:AbstractVector}
    C = ∂.C
    mul!(C, ∂.A, ∂.B)
    lmul!(∂.α, C)
    C
end

for op in [:(+), :(-)]
    for Mat in [:Diagonal, :Tridiagonal, :SymTridiagonal, :UniformScaling]
        @eval begin
            function Base.$op(∂::ND, B::$Mat) where {ND<:NumerovDerivative}
                B̃ = inv(∂.c)*∂.M*B
                NumerovDerivative($op(∂.Δ, B̃), ∂.M, ∂.M⁻¹, ∂.c)
            end
        end
    end
end

struct NumerovFactorization{TriFact,Mat}
    Δ⁻¹::TriFact
    M::Mat
end

Base.size(∂⁻¹::NF, args...) where {NF<:NumerovFactorization} = size(∂⁻¹.M, args...)
Base.eltype(∂⁻¹::NF) where {NF<:NumerovFactorization} = eltype(∂⁻¹.M)

LinearAlgebra.factorize(∂::ND) where {ND<:NumerovDerivative} =
    NumerovFactorization(factorize(∂.c*∂.Δ), ∂.M)

function LinearAlgebra.ldiv!(y, ∂⁻¹::NF, x) where {NF<:NumerovFactorization}
    mul!(y, ∂⁻¹.M, x)
    ldiv!(∂⁻¹.Δ⁻¹, y)
    y
end

# * Projections

# Vandermonde interpolation for finite differences is equivalent to
# evaluating the function on the grid points, since the basis
# functions are orthogonal (in the sense of a quadrature) and there is
# no overlap between adjacent basis functions.
function Base.:(\ )(B::FD, f::BroadcastQuasiArray) where {T,FD<:AbstractFiniteDifferences{T}}
    axes(f,1) == axes(B,1) ||
        throw(DimensionMismatch("Function on $(axes(f,1).domain) cannot be interpolated over basis on $(axes(B,1).domain)"))

    x = locs(B)
    getindex.(Ref(f), x) ./ weight.(Ref(B), eachindex(x))
end


text:
#----------------------------#
# see page 57 in Sarkka book #
#----------------------------#

function filteringrecursion(y; A = A, H = H, Q = Q, R = R, m₀ = m₀, P₀ = P₀)

    N = length(y)

    # Store here state distribution

    μstate, Σstate = Array{Vector{Float64}}(undef, N), Array{Matrix{Float64}}(undef, N)

    # Start recursion with prior mean and prior covariance

    mₖ₋₁, Pₖ₋₁ = m₀, P₀

    for k in 1:N

        # prediction step

        mₖ⁻ = A * mₖ₋₁

        Pₖ⁻ = A * Pₖ₋₁ * A' + Q     # (4.20)

        # update step

        𝐯ₖ = y[k] - H * mₖ⁻

        Sₖ = H * Pₖ⁻ * H' + R

        Kₖ = (Pₖ⁻ * H') / Sₖ

        mₖ = mₖ⁻ + Kₖ * 𝐯ₖ

        Pₖ = Pₖ⁻ - Kₖ * Sₖ * (Kₖ)'  # (4.21)

        # store mean and covariance of Gaussian distribution N(xₖ|mₖ, Pₖ)

        Pₖ = (Pₖ+(Pₖ)') / 2

        μstate[k], Σstate[k] = mₖ, Pₖ

        # next state

        mₖ₋₁, Pₖ₋₁ = mₖ, Pₖ

    end

    return μstate, Σstate

end

function testfiltering() # Works ✅

    yclean, y, x, A, H, Q, R, m₀, P₀ = simulatedata()

    μ, Σ = filteringrecursion(y; A = A, H = H, Q = Q, R = R, m₀ = m₀, P₀ = P₀)

    figure()

    subplot(211)
    plot([yᵢ[1]    for yᵢ in yclean], label = "simulated 1 clean")
    plot([yᵢ[1]    for yᵢ in y], label = "simulated 1")
    plot([(H*m)[1] for m  in μ], label = "mean 1")
    legend()

    subplot(212)
    plot([yᵢ[2]    for yᵢ in y], label = "simulated 2")
    plot([(H*m)[2] for m  in μ], label = "mean 2")

    legend()

end


text:
#note: may have to remove the 'data' folder from this directory and run 'build PowerSystems'
base_dir = dirname(dirname(pathof(PowerSystems)))

ps_dict = PowerSystems.parsestandardfiles(joinpath(MATPOWER_DIR, "case5_re.m"))

buses, generators, storage, branches, loads, loadZones, shunts, services =
    PowerSystems.ps_dict2ps_struct(ps_dict);
sys = PowerSystems.System(buses, generators, loads, branches, storage,
                               ps_dict["baseMVA"]);

@testset "Check bus index" begin
    @test sort([b.number for b in sys.buses]) == [1, 2, 3, 4, 5]
    @test sort(collect(Set([b.connectionpoints.from.number for b in sys.branches]))) ==
        [1, 2, 3, 4]
    @test sort(collect(Set([b.connectionpoints.to.number for b in sys.branches]))) ==
        [2, 3, 4, 5]

    # TODO: add test for loadzones testing MAPPING_BUSNUMBER2INDEX

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/safe_mode.jl
# UNDER DEVELOPMENT. CURRENTLY NOT USED. Implements a correctly-rounded
# version of standard McCormick arithmetic.
#############################################################################

#=
# Need to add this to Project.toml... requires Julia 1.3+
using CRLibm
using RoundingEmulator

struct NSSafe <: RelaxTag end

# Safe addition
@inline function plus_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
	cv = add_down(x.cv, y.cv)
	cc = add_up(x.cc, y.cc)
	MC{N,T}(cv, cc, z, x.cv_grad + y.cv_grad, x.cc_grad + y.cc_grad, (x.cnst && y.cnst))
end

# Safe scalar addition
@inline function plus_kernel(y::Float64, x::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
	cv = add_down(x.cv, y)
	cc = add_up(x.cc, y)
	MC{N,T}(cv, cc, z, x.cv_grad, x.cc_grad, x.cnst)
end

# Safe subtraction
@inline function minus_kernel(x::MC{N,T}, y::MC{N,T}, z::Interval{Float64}) where {N, T <:  Union{NSSafe}}
	cv = sub_down(x.cv, y.cc)
	cc = sub_up(x.cc, y.cv)
	MC{N,T}(cv, cc, z, x.cv_grad - y.cc_grad, x.cc_grad - y.cv_grad, (x.cnst && y.cnst))
end

# Safe scalar subtraction
@inline function minus_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
	cv = sub_down(x.cv, c)
	cc = sub_up(x.cc, c)
	MC{N,T}(cv, cc, z, x.cv_grad, x.cc_grad, x.cnst)
end
@inline function minus_kernel(c::Float64, x::MC{N,T}, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
	cv = sub_down(c, x.cc)
	cc = sub_up(c, x.cv)
	MC{N,T}(cv, cc, z, -x.cc_grad, -x.cv_grad, x.cnst)
end

# Safe Scalar Multiplication
@inline function mult_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
	if c >= 0.0
		cv = mul_down(c, x.cv)
		cc = mul_up(c, x.cc)
		zMC = MC{N,T}(cv, cc, z, c*x.cv_grad, c*x.cc_grad, x.cnst)
	else
		cv = mul_down(c, x.cc)
		cc = mul_up(c, x.cv)
		zMC = MC{N,T}(cc, cv, z, c*x.cc_grad, c*x.cv_grad, x.cnst)
	end
	return zMC
end

# Scalar Safe Kernel
@inline function div_kernel(x::MC{N,T}, c::Float64, z::Interval{Float64}) where {N, T <: Union{NSSafe}}
	if c >= 0.0
		cv = div_down(x.cv, c)
		cc = div_up(x.cc, c)
	else
		cc = div_down(x.cv, c)
		cv = div_up(x.cc, c)
	end
	return MC{N,T}(cc, cv, z, c*x.cc_grad, c*x.cv_grad, x.cnst)
end
=#


text:
#=
    Implements mpohamiltonian algebra (addition,subtraction,multiplication)
    other actions on mpohamiltonian objects
    the mpohamiltonian transfers
=#

#addition / substraction
function Base.:+(a::MPOHamiltonian{S,T,E},e::AbstractArray{V,1}) where {S,T,E,V}

    @assert length(e) == a.period

    nOs = deepcopy(a.Os) # we don't want our addition to change different copies of the original hamiltonian

    for c = 1:a.period
        if nOs[c,1,end] isa E
            nOs[c,1,end] += e[c]
        else
            @tensor nOs[c,1,a.odim][-1 -2;-3 -4]:=a[c,1,a.odim][-1,-2,-3,-4]+(e[c]*one(a[c,1,a.odim]))[-1,-2,-3,-4]
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,a.domspaces,a.pspaces)
end
Base.:-(e::Array{T,1},a::MPOHamiltonian) where T = -1.0*a + e
Base.:+(e::Array{T,1},a::MPOHamiltonian) where T = a + e
Base.:-(a::MPOHamiltonian,e::AbstractArray{T,1}) where T = a + (-e)

function Base.:+(a::MPOHamiltonian{S,T,E},b::MPOHamiltonian{S,T,E}) where {S,T,E}
    @assert a.period == b.period
    @assert sanitycheck(a)
    @assert sanitycheck(b)

    #new odim;domspaces,imspaces (pspaces necessarily stays the same)
    nodim = a.odim+b.odim-2
    ndomspaces = PeriodicArray{S,2}(undef,a.period,nodim)

    nOs = PeriodicArray{Union{E,T},3}(fill(zero(E),a.period,nodim,nodim))

    for pos in 1:a.period
        ndomspace = Vector{S}(undef,nodim)

        #not entirely sure :3
        ndomspace[1:(a.odim-1)]=a.domspaces[pos,1:a.odim-1]
        ndomspace[a.odim:a.odim+b.odim-2]=b.domspaces[pos,2:b.odim]
        ndomspaces[pos,:] = ndomspace[:]

        for (i,j) in keys(a,pos)
            #A block
            if(i<a.odim && j<a.odim)
                nOs[pos,i,j]=a.Os[pos,i,j]
            end

            #right side
            if(i<a.odim && j==a.odim)
                nOs[pos,i,nodim]=a.Os[pos,i,j]
            end
        end

        for (i,j) in keys(b,pos)

            #upper Bs
            if(i==1 && j>1)
                if nOs[pos,1,a.odim+j-2] isa T
                    nOs[pos,1,a.odim+j-2] += b[pos,i,j]
                else
                    if b.Os[pos,i,j] isa T
                        if nOs[pos,1,a.odim+j-2] == zero(E)
                            nOs[pos,1,a.odim+j-2] = b[pos,i,j] #one doesn't exist when dom != codom; hacky workaround
                        else
                            nOs[pos,1,a.odim+j-2] = nOs[pos,1,a.odim+j-2]*one(b[pos,i,j])+b[pos,i,j]
                        end
                    else
                        nOs[pos,1,a.odim+j-2] += b.Os[pos,i,j]
                    end
                end
            end

            #B block
            if(i>1 && j>1)
                nOs[pos,a.odim+i-2,a.odim+j-2]=b.Os[pos,i,j]
            end
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,ndomspaces,a.pspaces)
end
Base.:-(a::MPOHamiltonian,b::MPOHamiltonian) = a+(-1.0*b)

#multiplication
Base.:*(b::Number,a::MPOHamiltonian)=a*b
function Base.:*(a::MPOHamiltonian{S,T,E},b::Number) where {S,T,E}
    nOs=deepcopy(a.Os)

    for i=1:a.period
        for j=1:(a.odim-1)
            if(contains(a,i,j,a.odim))
                nOs[i,j,a.odim]*=b
            end
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,a.domspaces,a.pspaces)
end

#this is the index-map used in the ham x ham multiplication function (also needed somewhere else)
#i think julia has a build in for this, but it got renamed somewhere (linearindices?)
multmap(a::MPOHamiltonian,b::MPOHamiltonian) = (i,j)->(i-1)*b.odim+j
function Base.:*(b::MPOHamiltonian{S,T,E},a::MPOHamiltonian{S,T,E}) where {S,T,E}
    nodim=a.odim*b.odim

    indmap=multmap(a,b)

    nOs = PeriodicArray{Union{E,T},3}(fill(zero(E),a.period,nodim,nodim))

    ndomspaces = PeriodicArray{S,2}(undef,a.period,nodim)

    for pos=1:a.period
        ndomspace=Array{eltype(a.domspaces[pos]),1}(undef,nodim)

        for i in 1:a.odim
            for j in 1:b.odim
                ndomspace[indmap(i,j)]=fuse(a.domspaces[pos,i]*b.domspaces[pos,j])
            end
        end
        ndomspaces[pos,:]=ndomspace[:]
    end

    for pos=1:a.period
        for (i,j) in keys(a,pos)
            for (k,l) in keys(b,pos)
                if isscal(a,pos,i,j) && isscal(b,pos,k,l)
                    nOs[pos,indmap(i,k),indmap(j,l)] = a.Os[pos,i,j]*b.Os[pos,k,l]
                else
                    @tensor newopp[-1 -2;-3 -4 -5 -6]:=a[pos,i,j][-1,1,-4,-6]*b[pos,k,l][-2,-3,-5,1]
                    newopp=TensorMap(newopp.data,ndomspaces[pos,indmap(i,k)],domain(newopp))
                    newopp=permute(newopp,(1,2,5),(3,4))
                    newopp=TensorMap(newopp.data,codomain(newopp),ndomspaces[pos+1,indmap(j,l)])
                    newopp=permute(newopp,(1,2),(4,3))

                    nOs[pos,indmap(i,k),indmap(j,l)]=newopp
                end
            end
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,ndomspaces,a.pspaces)
end

#without the copy, we get side effects when repeating + setindex
Base.repeat(x::MPOHamiltonian{S,T,E},n::Int) where {S,T,E} =
    MPOHamiltonian{S,T,E}(
                                            repeat(x.Os,n,1,1),
                                            repeat(x.domspaces,n,1),
                                            repeat(x.pspaces,n))

#transpo = false => inplace conjugate
#transpo = true => flip physical legs
function Base.conj(a::MPOHamiltonian;transpo=false)
    b = deepcopy(a.Os)

    for (i,j,k) in keys(a)
        b[i,j,k] = @tensor temp[-1 -2;-3 -4]:=conj(a[i,j,k][-1,-2,-3,-4])
        if transpo
            b[i,j,k]=permute(b[i,j,k],(1,4),(3,2))
        end
    end

    MPOHamiltonian(b)
end

#needed this; perhaps move to tensorkit?
TensorKit.fuse(f::T) where T<: VectorSpace = f

#the usual mpoham transfer
function transfer_left(vec::Array{V,1},ham::MPOHamiltonian,pos::Int,A::V,Ab::V=A) where V<:MPSTensor
    toreturn = Array{V,1}(undef,length(vec));
    assigned = [false for i in 1:ham.odim]

    for (j,k) in keys(ham,pos)
        if assigned[k]
            if isscal(ham,pos,j,k)
                toreturn[k]+=ham.Os[pos,j,k]*transfer_left(vec[j],A,Ab)
            else
                toreturn[k]+=transfer_left(vec[j],ham[pos,j,k],A,Ab)
            end
        else
            if isscal(ham,pos,j,k)
                toreturn[k]=ham.Os[pos,j,k]*transfer_left(vec[j],A,Ab)
            else
                toreturn[k]=transfer_left(vec[j],ham[pos,j,k],A,Ab)
            end
            assigned[k]=true
        end
    end

    for k in 1:ham.odim
        if !assigned[k]
            #prefereably this never happens, because it's a wasted step
            #it's also avoideable with a little bit more code
            toreturn[k]=transfer_left(vec[1],ham[pos,1,k],A,Ab)
        end
    end

    return toreturn
end
function transfer_right(vec::Array{V,1},ham::MPOHamiltonian,pos::Int,A::V,Ab::V=A) where V<:MPSTensor
    toreturn = Array{V,1}(undef,length(vec));
    assigned = [false for i in 1:ham.odim]

    for (j,k) in keys(ham,pos)
        if assigned[j]
            if isscal(ham,pos,j,k)
                toreturn[j]+=ham.Os[pos,j,k]*transfer_right(vec[k],A,Ab)
            else
                toreturn[j]+=transfer_right(vec[k],ham[pos,j,k],A,Ab)
            end

        else
            if isscal(ham,pos,j,k)
                toreturn[j]=ham.Os[pos,j,k]*transfer_right(vec[k],A,Ab)
            else
                toreturn[j]=transfer_right(vec[k],ham[pos,j,k],A,Ab)
            end
            assigned[j]=true
        end
    end

    for j in 1:ham.odim
        if !assigned[j]
            toreturn[j]=transfer_right(vec[1],ham[pos,j,1],A,Ab)
        end
    end

    return toreturn
end


text:
using Morbit
using Test

const m = Morbit;

# TODO test for different configurations (RbfConfig, LagrangeConfig etc.)

@testset "Adding Objectives" begin
    global mop, lb, ub;
    global f1, f2, f3, F56, F78

    mop_unconstrained = MixedMOP()
    @test !mop_unconstrained.is_constrained

    ub = 4 .* ones(2);
    lb = -ub;
    mop = MixedMOP( lb = lb, ub = ub )
    @test mop.is_constrained

    # 1) add a scalar expensive objective
    f1(x) = (x[1]-1)^2 + (x[2] -1)^2
    add_objective!( mop , f1, :expensive )
    @test mop.n_objfs == 1

    # 2a) add a scalar cheap objective with gradient provided
    function f2(x)
        (x[1] + 1)^2 + (x[2]+1)^2;
    end
    ∇f2(x) = [ 2*(x[1]+1); 2*(x[2]+1) ];
    add_objective!( mop, f2, ∇f2 );
    @test mop.n_objfs == 2

    # 2b) add scalar cheap objective with autodiff enabled
    f3 = x -> 2 * x[1] + 2 * x[2];
    add_objective!( mop, f3, :cheap );
    @test mop.n_objfs == 3

    # test internal indexing and sorting
    m.set_non_exact_indices!(mop)
    @test mop.non_exact_indices == [1,]

    m.set_sorting!(mop)
    @test mop.internal_sorting == [1, 2, 3]
    add_objective!( mop , f1, :expensive )
    m.set_sorting!(mop)
    @test mop.internal_sorting == [2,3,1,4]
    @test mop.reverse_sorting == [3,1,2,4]
    @test m.apply_internal_sorting( mop, [1,2,3,4] ) == mop.internal_sorting

    # 3) add vector cheap objective ( not super effective )
    function F56(x)
        y1 = exp( sum(x) );
        y2 = sin( x[1] );
        return [y1;y2]
    end

    add_objective!( mop, F56, :cheap, 2 )
    @test mop.n_objfs == 6

    # 4) add vector expensive objective ( not super effective )
    function F78(x)
        y1 = sum(x);
        y2 = 1.0;
        return [y1;y2]
    end

    add_objective!( mop, F78, :expensive, 2 )
    @test mop.n_objfs == 8
    m.set_sorting!(mop)
    m.set_non_exact_indices!(mop)
    # now 1,4,7,8 (and 5,6) should have been grouped together
    @test mop.internal_sorting == [2, 3, 5, 6, 1, 4, 7, 8]
    @test mop.non_exact_indices == [5,6,7,8]
end

@testset "Objective Evaluations" begin
    global mop, lb, ub;
    global f1, f2, f3, F56, F78 # f1 & F78 are expensive

    # test evaluation for constrained problem

    X = lb .+ (ub .- lb) .* rand(2);
    @test m.unscale(mop, m.scale(mop, X)) ≈ X

    # internally all vectors are in [0,1]^n
    x = m.scale(mop, X);

    all_eval = m.eval_all_objectives( mop, x);
    exp_eval = all_eval[mop.non_exact_indices];

    exp_comparison = [f1(X); f1(X); F78(X)][:];
    @test exp_eval ≈ exp_comparison

    Y = [ f1(X); f2(X); f3(X); f1(X); F56(X); F78(X) ][:];
    @test all_eval ≈ m.apply_internal_sorting(mop, Y)
    @test m.reverse_internal_sorting(mop, all_eval) ≈ Y
end

@testset "Internal Evaluations" begin
    lb = zeros(2);
    ub = 4 .* ones(2);

    s = [ 4 .* rand(2) for i=1:10 ];

    f1(x) = 1.0;
    f2(x) = (x[1]+1)^2 + (x[2] +1)^2;
    f(x) = [ f1(x); f2(x) ];

    p = MixedMOP(lb = lb, ub = ub);
    add_objective!(p,f1,:cheap,1);
    add_objective!(p,f2,:expensive,1);

    eval1 = Morbit.eval_all_objectives(p, s[1])
    eval2 = Morbit.eval_all_objectives.(p, s)

    p = MixedMOP(lb=lb, ub = ub)
    add_objective!(p,f,:expensive, 2)
    eval3 = m.eval_all_objectives.(p,s);

    global counter = 0;
    function batch_fun(x)
        global counter
        counter += 1
        if isa(x, Vector{T} where{T<:Real})
            return f(x)
        else
            return f.(x)
        end
    end

    p = MixedMOP(lb=lb, ub = ub);
    add_objective!(p, batch_fun, :expensive, 2, true);
    eval4 = Morbit.eval_all_objectives.(p,s);

    @test isa( eval1, Vector{Float64})
    @test length(eval1) == 2
    @test all( isa(e, Vector{Vector{Float64}} ) for e in [eval2, eval3, eval4])
    @test eval2[1] == eval3[1]
    @test eval3 == eval4
    @test counter == 1
end


text:
""
function expression_bus_generation(pm::AbstractPowerModel, n::Int, i::Int, bus_gens)
    pg = get(var(pm, n), :pg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "active power", "generator")
    qg = get(var(pm, n), :qg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "reactive power", "generator")

    pg_total = 0.0
    qg_total = 0.0
    if length(bus_gens) > 0
        pg_total = sum(pg[g] for g in bus_gens)
        qg_total = sum(qg[g] for g in bus_gens)
    end

    var(pm, n, :bus_pg)[i] = pg_total
    var(pm, n, :bus_qg)[i] = qg_total
end


text:
using LaTeXChanges
using Test

test = joinpath(@__DIR__, "test.tex")
test_clean = joinpath(@__DIR__, "test_clean.tex")
test_applied = joinpath(@__DIR__, "test_applied.tex")

apply_changes(test, test_applied)

@testset "LaTeXChanges.jl" begin
    @test readlines(test_clean) == readlines(test_applied)
end


text:
# indicator of a simplex

export IndSimplex

"""
**Indicator of a simplex**

    IndSimplex(a=1.0)

Returns the indicator of the set
```math
S = \\left\\{ x : x \\geq 0, \\sum_i x_i = a \\right\\}.
```
By default `a=1.0`, therefore ``S`` is the probability simplex.
"""
struct IndSimplex{R <: Real} <: ProximableFunction
    a::R
    function IndSimplex{R}(a::R) where {R <: Real}
        if a <= 0
            error("parameter a must be positive")
        else
            new(a)
        end
    end
end

is_convex(f::IndSimplex) = true
is_set(f::IndSimplex) = true

IndSimplex(a::R=1.0) where {R <: Real} = IndSimplex{R}(a)

function (f::IndSimplex)(x::AbstractArray{R}) where R <: Real
    if all(x .>= 0) && sum(x) ≈ f.a
        return R(0)
    end
    return R(Inf)
end

function simplex_proj_condat!(y::AbstractArray{R}, a, x::AbstractArray{R}) where R
    # Implements algorithm proposed in:
    # Condat, L. "Fast projection onto the simplex and the l1 ball",
    # Mathematical Programming, 158:575–585, 2016.
    v = [x[1]]
    v_tilde = R[]
    rho = x[1] - a
    N = length(x)
    for k in 2:N
        if x[k] > rho
            rho += (x[k] - rho) / (length(v) + 1)
            if rho > x[k] - a
                push!(v, x[k])
            else
                append!(v_tilde, v)
                v = [x[k]]
                rho = x[k] - a
            end
        end
    end
    for z in v_tilde
        if z > rho
            push!(v, z)
            rho += (z - rho) / length(v)
        end
    end
    v_changed = true
    while v_changed == true
        v_changed = false
        k = 1
        while k <= length(v)
            z = v[k]
            if z <= rho
                deleteat!(v, k)
                v_changed = true
                rho += (rho - z) / length(v)
            else
                k = k + 1
            end
        end
    end
    y .= max.(x .- rho, R(0))
end

function prox!(y::AbstractArray{R}, f::IndSimplex, x::AbstractArray{R}, _::R=R(1)) where R <: Real
    simplex_proj_condat!(y, f.a, x)
    return R(0)
end

fun_name(f::IndSimplex) = "indicator of the probability simplex"
fun_dom(f::IndSimplex) = "AbstractArray{Real}"
fun_expr(f::IndSimplex) = "x ↦ 0 if x ⩾ 0 and sum(x) = a, +∞ otherwise"
fun_params(f::IndSimplex) = "a = $(f.a)"

function prox_naive(f::IndSimplex, x::AbstractArray{R}, _::R=R(1)) where R <: Real
    low = minimum(x)
    upp = maximum(x)
    v = x
    s = Inf
    for i = 1:100
        if abs(s)/f.a ≈ 0
            break
        end
        alpha = (low+upp)/2
        v = max.(x .- alpha, R(0))
        s = sum(v) - f.a
        if s <= 0
            upp = alpha
        else
            low = alpha
        end
    end
    return v, R(0)
end


text:
const DefaultPointsPerArc = 64

"""
    catmullrom(points, pointsperarc; extend=true)
    catmullrom(xs, ys, pointsperarc; extend=true)
    catmullrom(xs, ys, zs, pointsperarc; extend=true)

Given abcissa-sequenced path of points (as points or as xs and ys), and
the number of subdivisions to be fit inbetween pair of points
(all neighboring points except the first & last),
obtain the centripetal Catmull-Rom splines that cover each segment.

By default, extrapolated first and final points
are affixed to the input point sequence so that
all of the input points are present in the resulting sequence
(CatmullRom fitting does not include the first and final points).
If you prefer this not occur, call the function with `extend=false`.

If you prefer to specify the scale factor used in that extrapolation,
use `extendbounds(points, scale=scalefactor)`, and then pass the result
to this function with `extend=false`.
"""
function catmullrom(points::P, pointsperarc::Integer=DefaultPointsPerArc; extend::Bool=true) where P
    catmullrom_requirement(npoints(points))
    if eltype(points) <: Tuple
        apoints = map(x->[x...,], points)
        return catmullrom(apoints, pointsperarc, extend=extend)
    end
    pointsperarc += isodd(pointsperarc)     # force even

    # ensure that the 'x' values are not coinciding
    changes = norm.(diff(points))[1:end]
    relchanges = changes ./ sum(changes)
    cumrelchanges = cumsum(relchanges)
    pushfirst!(cumrelchanges, 0.0)
    m = vcat(cumrelchanges',reduce(hcat,points))
    m = permutedims(m)
    crpoints = [m[i,:] for i=1:size(m)[1]]

    crpoints = catmullrom_prep(crpoints, pointsperarc, extend=extend)

    # ensure the spline passes through the original values
    return catmullrom_splines(crpoints, pointsperarc)[2:end]
end

function catmullrom_prep(points::P, pointsperarc::Integer=DefaultPointsPerArc; extend::Bool=true) where P
    if isa(points[1], Tuple)
        cr_points = [Float64.([x...,]) for x in points]
    else
        cr_points = [Float64.(x) for x in points]
    end
    if extend
        cr_points = extend_seq(cr_points)
    end
    return cr_points
end

function catmullrom_splines(points::P, pointsperarc::Integer) where P
    coord_type = coordtype(points)

    vals_along_each_coord = catmullrom_core(points, pointsperarc)

    if coord_type === Float64
        return vals_along_each_coord
    else
        return [map(coord_type, vals) for vals in vals_along_each_coord]
    end
end

function catmullrom_splines(points::P, pointsperarc::Vector{Integer}) where P
    n_points = npoints(points)
    n_ppa    = length(pointsperarc)
    n_points != n_ppa+1 && throw(ErrorException("length(points) != length(pointsperarc)+1 ($(n_points) != $(n_ppa+1))"))

    vals_along_each_coord = []
    for idx = 1:n_ppa
        pts = points[idx:idx+3]
        ppa = pointsperarc[idx]
        vals_each_coord = catmullrom_core(points, ppa)
        push!(vals_along_each_coord, vals_each_coord)
    end
    return vals_along_each_coord

    coord_type = coordtype(points)
    if coord_type === Float64
        return vals_along_each_coord
    else
        return [map(coord_type, vals) for vals in vals_along_each_coord]
    end
end

function catmullrom_core(points::P, pointsperarc::Integer) where P
    catmullrom_requirement(npoints(points))
    n_coords = ncoords(points)

    # include the given points (knots) for poly generation
    n_through_points = pointsperarc + 2  # include both endpoints

    #=
        Separate sequences for each coordinate dimension, ordered first to last interpoint span.
        Each poly within a sequence fits an interpoint span for its respective coordinate axis.
        _polys_ is an `(n_points - 3) x (n_coords) 2D array`
           columns are of one coordinate axis
           each row holds polys that interpolate one interpoint span across all coordinates.
           (n_points - 2 - 1): drop two extremal points, spans are counted between fenceposts
    =#
    polys = catmullrom_polys(points)

    # interpolate using span-relative abcissae [0.0..1.0]
    abcissae01 = collect(range(0.0, stop=1.0, length=n_through_points))

    # vals is an `(n_points - 3) x (n_dims)` array of groups of successive point_along_coord placements

    xvals = flat(map(ply, abcissae01[1:end-1]) for ply in polys[:,1])
    yvals = flat(map(ply, abcissae01[1:end-1]) for ply in polys[:,2])
    coord_vals = collect(zip(xvals,yvals))

    # vals = polys.(Ref(abcissae01[1:end-1]))

    #coord_vals = [collect(Iterators.flatten( vals[:,i] )) for i=1:n_coords]
    #endvals = (polys[end,:]).(fill(1.0, n_coords))
    #@. push!(coord_vals, endvals)

    return coord_vals
end

flat(x) = collect(Base.Iterators.flatten(x))

"""
    catmullrom_polys(points)

Traverse a sequence of points (coordinates in 2,3,..n dims),
obtaining interpolatory centripetal Catmull-Rom polynomials.

Yields separate sequences for each coordinate dimension,
    each sequence is ordered from first to last interpoint span.

Each poly within a sequence fits an interpoint span with respect to
    the sequence's coordinate axis.

_polys_ is an `(n_points - 3) x (n_coords) 2D array`
   columns are of one coordinate axis
   each row holds polys that interpolate one interpoint span across all coordinates.
"""
function catmullrom_polys(points::P) where P
    n_points = npoints(points)
    catmullrom_requirement(n_points)

    T = coordtype(points)
                                    # omit the extremal points (-2)
    n_spans  = n_points - 3         # count between fenceposts (-1)
    n_coords = ncoords(points)      # each point has ncoordinates

    polys  = Array{Polynomial{T}, 2}(undef, n_spans, n_coords)

    for idx = 1:n_spans
        pt₋, pt₀, pt₁, pt₊ = points[idx:idx+3]
        plys = centripetal_catmullrom(pt₋, pt₀, pt₁, pt₊)
        polys[idx,:] = plys
    end

    return polys
end

"""
    centripetal_catmullrom(pt1, pt2, pt3, pt4)

Given four points on a path in their given sequence,
obtain the centripetal Catmull-Rom cubic polynomial
that interpolates the middle two points of the path.
"""
function centripetal_catmullrom(pt₋::T,
                                pt₀::T,
                                pt₁::T,
                                pt₊::T) where T

    tg₀,tg₁     = centripetal_tangents(pt₋, pt₀, pt₁, pt₊)
    c₀,c₁,c₂,c₃ = centripetal_hermite(pt₀, pt₁, tg₀, tg₁)

    coeffs = [[x...] for x in zip(c₀,c₁,c₂,c₃)]
    return Polynomial.(coeffs)
end

"""
    centripetal_tangents(pt1, pt2, pt3, pt4)

Given four points on a path in their given sequence,
obtain the tangents at the middle two points of the path,
pt2, pt3 delimiting the centripetal curve segment,
each as a unit vector.
"""
function centripetal_tangents(pt₋::T,
                              pt₀::T,
                              pt₁::T,
                              pt₊::T) where T

    # intrasegment centripetal speeds
    dt₋₀ = speed(pt₋, pt₀)
    dt₀₁ = speed(pt₀, pt₁)
    dt₁₊ = speed(pt₁, pt₊)

    dt₋₀, dt₀₁, dt₁₊ = prevent_overlap(T, dt₋₀, dt₀₁, dt₁₊)

    # Compute the tangents at the two boundary points
    #   that delimit the interpolatory curve segment.
    #   - tg₀ is the tangent at pt₀
    #   - tg₁ is the tangent at pt₁

    tg₀ = @. (pt₀ - pt₋) /  dt₋₀ -
             (pt₁ - pt₋) / (dt₋₀ + dt₀₁) +
             (pt₁ - pt₀) /  dt₀₁

    tg₁ = @. (pt₁ - pt₀) /  dt₀₁ -
             (pt₊ - pt₀) / (dt₀₁ + dt₁₊) +
             (pt₊ - pt₁) /  dt₁₊

    # Normalize the tangent vectors. These unit vector tangents
    #   provide (guide) centripetal interpolation in [0..1].
    tg₀ = tg₀ .* dt₀₁
    tg₁ = tg₁ .* dt₀₁

    return (tg₀, tg₁)
end

@inline function prevent_overlap(::Type{P}, dt₋₀::T, dt₀₁::T, dt₁₊::T) where {P,T}
    # check if any coordinates coincidee
    ε = coordeps(P)
    # correct for repeated coordinates
    dt₀₁ = @. ifelse(dt₀₁ <= ε, one(coordtype(P)), dt₀₁)
    dt₋₀ = @. ifelse(dt₋₀ <= ε, dt₀₁, dt₋₀)
    dt₁₊ = @. ifelse(dt₁₊ <= ε, dt₀₁, dt₁₊)
    return dt₋₀, dt₀₁, dt₁₊
 end

"""
    centripetal_hermite(pt₀, pt₁, tg₀, tg₁)

- ply(t)  = c0 + c1*t   +   c2*t^2 + c3*t^3
- dply(t) = c1 + 2*c2*t + 3*c3*t^2

Given the bounding points `(pt₀, pt₁)` of a  curve segment
and the unit tangent vectors `(tg₀, tg₁)` at those points,
determine the coeffs `(c0, c1, c2, c3)` of cubic `ply(t)`
and of associated quadratic `dply(t)` such that

- ply(0) == pt₀,  dply(0) == tg₀
- ply(1) == pt₁,  dply(1) == tg₁
"""
function centripetal_hermite(pt₀::T, pt₁::T, tg₀::G, tg₁::G) where {T, G}
    c0 = pt₀
    c1 = tg₀
    c2 = @. -3*pt₀ + 3*pt₁ - 2*tg₀ - tg₁
    c3 = @.  2*pt₀ - 2*pt₁ +   tg₀ + tg₁

    return (c0, c1, c2, c3)
end

"""
    speed(pt₀, pt₁)

Obtain the tangent vector along the path
at pt₀ moving to pt₁. The sqrt of that
magnitude is the centripetal speed.
"""
function speed(pt₀::T, pt₁::T) where T
    delta = pt₁ .- pt₀
    return sqrt(sqrtdot(delta, delta))
end

@inline function coordeps(::Type{T}) where T
    coord_type = coordtype(T)
    if coord_type <: AbstractFloat
        ε = eps(coord_type)
    elseif coord_type <: Integer
        ε = zero(coord_type)
    else
        ε = eps(float(dt₀₁))
    end
    return ε
end

sqrtdot(a::T,b::T) where {T} = sqrt(sum(a .* b))

function catmullrom_requirement(n_points::Int)
    n_points >= 4 || throw(DomainError("At least four points are required ($n_points points found)"))
end


text:
module GANs

using Flux
using Distributions
import Random
import Statistics
using LinearAlgebra

struct GAN
   g # generator
   gg # non-trainable generator copy
   d # discriminator
   dd # non-trainable discriminator copy
   pz # code distribution
end

Flux.@treelike GAN
freeze(m) = Flux.mapleaves(Flux.Tracker.data,m)

logpdf(b::Bernoulli, y::Bool) = y * log(b.p + eps()) + (1 - y) * log(1 - b.p + eps())

GAN(G::Flux.Chain, D::Flux.Chain; pz=randn) = GAN(G, freeze(G), D, freeze(D), pz)

getcode(gan::GAN) = Float64.(gan.pz(size(Flux.params(gan.g).order[1],2)))
getcode(gan::GAN, n::Int) = Float64.(gan.pz(size(Flux.params(gan.g).order[1],2), n))
generate(gan::GAN) = ((gan.g(getcode(gan)).data) .>= 0.5) .* 1.0
generate(gan::GAN, n::Int) = ((gan.g(getcode(gan, n)).data) .>= 0.5) .* 1.0

Dloss(gan::GAN, X, Z) = - Float64(0.5)*(Statistics.mean(log.(gan.d(X) .+ eps(Float64))) + Statistics.mean(log.(1 .- gan.d(gan.gg(Z)) .+ eps(Float64))))
Gloss(gan::GAN, Z) = - Statistics.mean(log.(gan.dd(gan.g(Z)) .+ eps(Float64)))

function initgan(p)
   generator = Chain(Dense(10, p, NNlib.leakyrelu), Dense(p, p, NNlib.σ))
   discriminator = Chain(Dense(p, 25, NNlib.leakyrelu), Dense(25, 1, NNlib.σ))
   gan = GAN(generator, discriminator, pz=randn)
end

function random_batch_index(x::AbstractArray, batch_size=1; dims=1)
   n = size(x,dims)
   Iterators.partition(Random.shuffle(1:n), batch_size)
end

function traingan!(gan::GAN, x; int_epochs::Int64=10,
      batch_size::Int64=10,etagen = 0.00001,etadisc=0.00001)

   Dopt=ADAM(etadisc)
   Gopt=ADAM(etagen)

   for epoch=1:int_epochs
      loss_epoch = [0.0, 0.0, 0]
      for I in random_batch_index(x, batch_size)
         m = length(I)
         # sample data and generate codes
         z = getcode(gan,m)

         # discriminator training
         Dl = Dloss(gan, x[I,:]',z)
         Flux.Tracker.back!(Dl)
         for param in Flux.params(gan.d)
            Flux.Tracker.update!(Dopt, param, Flux.Tracker.grad(param))
         end

         # generator training
         Gl = Gloss(gan, z)
         Flux.Tracker.back!(Gl)
         for param in Flux.params(gan.g)
            Flux.Tracker.update!(Gopt, param, Flux.Tracker.grad(param))
         end
      end
   end
end

function samples(gan, nsamples=1000)
   return convert(Array{Float64,2},Transpose(generate(gan,nsamples)))
end

end # module

max_stars_repo_name: stefan-m-lenz/dist-gen-comp

max_issues_repo_name: stefan-m-lenz/dist-gen-comp

max_forks_repo_name: stefan-m-lenz/dist-gen-comp


text:
"""
    ESS

Elliptical slice sampling algorithm.

# Examples
```jldoctest; setup = :(Random.seed!(1))
julia> @model gdemo(x) = begin
           m ~ Normal()
           x ~ Normal(m, 0.5)
       end
gdemo (generic function with 2 methods)

julia> sample(gdemo(1.0), ESS(), 1_000) |> mean
Mean

│ Row │ parameters │ mean     │
│     │ Symbol     │ Float64  │
├─────┼────────────┼──────────┤
│ 1   │ m          │ 0.824853 │
```
"""
struct ESS{space} <: InferenceAlgorithm end

ESS() = ESS{()}()
ESS(space::Symbol) = ESS{(space,)}()

mutable struct ESSState{V<:VarInfo} <: AbstractSamplerState
    vi::V
end

function Sampler(alg::ESS, model::Model, s::Selector)
    # sanity check
    vi = VarInfo(model)
    space = getspace(alg)
    vns = _getvns(vi, s, Val(space))
    length(vns) == 1 ||
        error("[ESS] does only support one variable ($(length(vns)) variables specified)")
    for vn in vns[1]
        dist = getdist(vi, vn)
        EllipticalSliceSampling.isgaussian(typeof(dist)) ||
            error("[ESS] only supports Gaussian prior distributions")
    end

    state = ESSState(vi)
    info = Dict{Symbol, Any}()

    return Sampler(alg, info, s, state)
end

# always accept in the first step
function AbstractMCMC.step!(
    ::AbstractRNG,
    model::Model,
    spl::Sampler{<:ESS},
    ::Integer,
    ::Nothing;
    kwargs...
)
    return Transition(spl)
end

function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:ESS},
    ::Integer,
    transition;
    kwargs...
)
    # obtain random variable that should be sampled
    vi = spl.state.vi

    # obtain previous sample
    f = vi[spl]

    # recompute log-likelihood in logp
    if spl.selector.tag !== :default
        runmodel!(model, vi, spl)
    end

    # define previous sampler state
    oldstate = EllipticalSliceSampling.EllipticalSliceSamplerState(f, getlogp(vi))

    # compute next state
    state = AbstractMCMC.step!(rng, ESSModel(model, spl),
                               EllipticalSliceSampling.EllipticalSliceSampler(), 1, oldstate;
                               progress = false)

    # update sample and log-likelihood
    vi[spl] = state.sample
    setlogp!(vi, state.loglikelihood)

    return Transition(spl)
end

struct ESSModel{M<:Model,S<:Sampler{<:ESS},T} <: AbstractMCMC.AbstractModel
    model::M
    spl::S
    μ::T
end

function ESSModel(model::Model, spl::Sampler{<:ESS})
    vi = spl.state.vi
    vns = _getvns(vi, spl)
    μ = mapreduce(vcat, vns[1]) do vn
        dist = getdist(vi, vn)
        vectorize(dist, mean(dist))
    end

    ESSModel(model, spl, μ)
end

# sample from the prior
function EllipticalSliceSampling.sample_prior(rng::Random.AbstractRNG, model::ESSModel)
    spl = model.spl
    vi = spl.state.vi
    vns = _getvns(vi, spl)
    set_flag!(vi, vns[1][1], "del")
    runmodel!(model.model, vi, spl)
    return vi[spl]
end

# compute proposal and apply correction for distributions with nonzero mean
function EllipticalSliceSampling.proposal(model::ESSModel, f, ν, θ)
    sinθ, cosθ = sincos(θ)
    a = 1 - (sinθ + cosθ)
    return @. cosθ * f + sinθ * ν + a * model.μ
end

function EllipticalSliceSampling.proposal!(out, model::ESSModel, f, ν, θ)
    sinθ, cosθ = sincos(θ)
    a = 1 - (sinθ + cosθ)
    @. out = cosθ * f + sinθ * ν + a * model.μ
    return out
end

# evaluate log-likelihood
function Distributions.loglikelihood(model::ESSModel, f)
    spl = model.spl
    vi = spl.state.vi
    vi[spl] = f
    runmodel!(model.model, vi, spl)
    getlogp(vi)
end

function DynamicPPL.tilde(ctx::DefaultContext, sampler::Sampler{<:ESS}, right, vn::VarName, inds, vi)
    if vn in getspace(sampler)
        return DynamicPPL.tilde(LikelihoodContext(), SampleFromPrior(), right, vn, inds, vi)
    else
        return DynamicPPL.tilde(ctx, SampleFromPrior(), right, vn, inds, vi)
    end
end

function DynamicPPL.tilde(ctx::DefaultContext, sampler::Sampler{<:ESS}, right, left, vi)
    return DynamicPPL.tilde(ctx, SampleFromPrior(), right, left, vi)
end

function DynamicPPL.dot_tilde(ctx::DefaultContext, sampler::Sampler{<:ESS}, right, left, vn::VarName, inds, vi)
    if vn in getspace(sampler)
        return DynamicPPL.dot_tilde(LikelihoodContext(), SampleFromPrior(), right, left, vn, inds, vi)
    else
        return DynamicPPL.dot_tilde(ctx, SampleFromPrior(), right, left, vn, inds, vi)
    end
end

function DynamicPPL.dot_tilde(ctx::DefaultContext, sampler::Sampler{<:ESS}, right, left, vi)
    return DynamicPPL.dot_tilde(ctx, SampleFromPrior(), right, left, vi)
end


text:
using CombinedParsers
using CombinedParsers.BNF
using CombinedParsers.Regexp
using Test

@testset "BNF" begin
    @test CombinedParsers.BNF.meta_identifier("abc") == :abc
    @test CombinedParsers.BNF.special_sequence("? ASCII char 32 ?") == "ASCII char 32"
    @test CombinedParsers.BNF.syntactic_primary("\"1b\"") == parser("1b")
    @test CombinedParsers.BNF.syntactic_factor("\"1b\"") == parser("1b")
    @test CombinedParsers.BNF.definitions_list("3 * \"aa\", \"B\"",log=true) == Repeat(3,3,"aa") * "B"
    @test parse(CombinedParsers.BNF.syntax_rule,"twelve = \"1\", \"2\" ;") == with_name(:twelve, Sequence("1","2"))
    @test CombinedParsers.BNF.definitions_list("digit excluding zero, { digit }") == Sequence(substitute(Symbol("digit excluding zero")),
                                                                                Repeat(substitute(:digit)))

    @test CombinedParsers.BNF.definitions_list("\"0\" | [ \"-\" ], natural number") == Either("0", Optional("-") * substitute("natural number"))

    bnf = parse(ebnf,
        """
        digit excluding zero = "1" | "2" | "3" | "4" | "5" | "6" | "7" | "8" | "9" ;
        digit                = "0" | digit excluding zero ;
        natural number       = digit excluding zero, { digit } ;
        integer              = "0" | [ "-" ], natural number ;
        """)

    @test bnf[:integer]("120") == ("", ("1", [ "2", "0" ]))
    @test deepmap(MatchedSubSequence, bnf, :integer)[:integer]("120") == "120"

    @test ebnf("""
    twelve a                          = "1", "2" ;
    """, log=true) isa NamedParser

    @test ebnf("""
    twelve                          = "1", "2" ; one hundred twelve = "1", twelve;
    """, log=true) isa Either

    @test ebnf("""
    twelve                          = "1", "2" ;    two hundredone                 = "2", "0", "1" ;
    """, log=true) isa Either

end


text:
# module Bukdu

function render(::Type{Text}, args...)::Conn
    filtering(render,Text,args...) do
        obj = isempty(args) ? "": first(args)
        Conn(200, Dict("Content-Type"=>"text/plain"), isa(obj, String) ? obj : string(obj))
    end
end


text:
abstract type AbstractAxis{IdxMap} end

@inline indexmap(::AbstractAxis{IdxMap}) where IdxMap = IdxMap
@inline indexmap(::Type{<:AbstractAxis{IdxMap}}) where IdxMap = IdxMap

# struct FlatAxis <: AbstractAxis{NamedTuple()} end

struct NullAxis <: AbstractAxis{nothing} end
const VarAxes = Tuple{Vararg{<:AbstractAxis}}

"""
    ax = Axis(nt::NamedTuple)

Gives named component access for `ComponentArray`s.
# Examples

```
julia> using ComponentArrays

julia> ax = Axis((a = 1, b = ViewAxis(2:7, PartitionedAxis(2, (a = 1, b = 2))), c = ViewAxis(8:10, (a = 1, b = 2:3))));

julia> A = [100, 4, 1.3, 1, 1, 4.4, 0.4, 2, 1, 45];

julia> ca = ComponentArray(A, ax)
ComponentVector{Float64}(a = 100.0, b = [(a = 4.0, b = 1.3), (a = 1.0, b = 1.0), (a = 4.4, b = 0.4)], c = (a = 2.0, b = [1.0, 45.0]))

julia> ca.a
100.0

julia> ca.b
3-element LazyArray{ComponentVector{Float64,SubArray...}}:
 ComponentVector{Float64,SubArray...}(a = 4.0, b = 1.3)
 ComponentVector{Float64,SubArray...}(a = 1.0, b = 1.0)
 ComponentVector{Float64,SubArray...}(a = 4.4, b = 0.4)

julia> ca.c
ComponentVector{Float64,SubArray...}(a = 2.0, b = [1.0, 45.0])

julia> ca.c.b
2-element view(::Vector{Float64}, 9:10) with eltype Float64:
  1.0
 45.0
```
"""
struct Axis{IdxMap} <: AbstractAxis{IdxMap} end
@inline Axis(IdxMap::NamedTuple) = Axis{IdxMap}()
Axis(;kwargs...) = Axis((;kwargs...))
function Axis(symbols::Union{AbstractVector{Symbol}, NTuple{N,Symbol}}) where {N}
    return Axis((;(symbols .=> eachindex(symbols))...))
end
Axis(symbols::Symbol...) = Axis(symbols)

const FlatAxis = Axis{NamedTuple()}
const NullorFlatAxis = Union{NullAxis, FlatAxis}

"""
    sa = ShapedAxis(shape, index_map)

Preserves higher-dimensional array components in `ComponentArray`s (matrix components, for
example)
"""
struct ShapedAxis{Shape, IdxMap} <: AbstractAxis{IdxMap} end
@inline ShapedAxis(Shape, IdxMap) = ShapedAxis{Shape, IdxMap}()
ShapedAxis(Shape) = ShapedAxis(Shape, NamedTuple())
ShapedAxis(::Tuple{<:Int}) = FlatAxis()

const Shape = ShapedAxis

unshape(ax) = ax
unshape(ax::ShapedAxis) = Axis(indexmap(ax))

Base.size(::ShapedAxis{Shape, IdxMap}) where {Shape, IdxMap} = Shape

"""
    pa = PartitionedAxis(partition_size, index_map)

Axis for creating arrays of `ComponentArray`s
"""
struct PartitionedAxis{PartSz, IdxMap, Ax<:AbstractAxis{IdxMap}} <: AbstractAxis{IdxMap}
    ax::Ax

    function PartitionedAxis(PartSz, ax::AbstractAxis{IdxMap}) where IdxMap
        return new{PartSz,IdxMap,typeof(ax)}(ax)
    end
end
PartitionedAxis{PartSz,IdxMap,Ax}() where {PartSz,IdxMap,Ax} = PartitionedAxis(PartSz, Ax())
PartitionedAxis(PartSz, IdxMap) = PartitionedAxis(PartSz, Axis(IdxMap))

const Partition = PartitionedAxis

Base.size(::PartitionedAxis{PartSz,IdxMap}) where {PartSz,IdxMap} = PartSz
Base.size(::Type{PartitionedAxis{PartSz,IdxMap}}) where {PartSz,IdxMap} = PartSz

"""
    va = ViewAxis(parent_index, index_map)

Axis for creating arrays of `ComponentArray`s
"""
struct ViewAxis{Inds, IdxMap, Ax<:AbstractAxis{IdxMap}} <: AbstractAxis{IdxMap}
    ax::Ax
    function ViewAxis(Inds, ax::AbstractAxis{IdxMap}) where IdxMap
        return new{Inds,IdxMap,typeof(ax)}(ax)
    end
    ViewAxis(Inds, ::NullorFlatAxis) = Inds
end
# ViewAxis{Inds,IdxMap,Ax}() where {Inds,IdxMap,Ax} = PartitionedAxis(Inds, Ax())
ViewAxis{Inds,IdxMap,Ax}() where {Inds,IdxMap,Ax} = ViewAxis(Inds, Ax())
ViewAxis(Inds, IdxMap) = ViewAxis(Inds, Axis(IdxMap))
ViewAxis(Inds) = Inds

const View = ViewAxis
const NullOrFlatView{Inds,IdxMap} = ViewAxis{Inds,IdxMap,<:NullorFlatAxis}

viewindex(::ViewAxis{Inds,IdxMap}) where {Inds,IdxMap} = Inds
viewindex(::Type{<:ViewAxis{Inds,IdxMap}}) where {Inds,IdxMap} = Inds
viewindex(i) = i

Axis(ax::AbstractAxis) = ax
Axis(ax::PartitionedAxis) = ax.ax
Axis(ax::ViewAxis) = ax.ax

# Get rid of this
Axis(::Number) = NullAxis()
Axis(::NamedTuple{()}) = FlatAxis()
Axis(x) = FlatAxis()

const NotShapedAxis = Union{Axis{IdxMap}, FlatAxis, NullAxis} where {IdxMap}
const NotPartitionedAxis = Union{Axis{IdxMap}, FlatAxis, NullAxis, ShapedAxis{Shape, IdxMap}} where {Shape, IdxMap}
const NotShapedOrPartitionedAxis = Union{Axis{IdxMap}, FlatAxis, NullAxis} where {IdxMap}

Base.merge(axs::Axis...) = Axis(merge(indexmap.(axs)...))

Base.lastindex(ax::AbstractAxis) = last(viewindex(last(indexmap(ax))))

Base.keys(ax::AbstractAxis) = keys(indexmap(ax))

reindex(i, offset) = i .+ offset
reindex(ax::FlatAxis, _) = ax
reindex(ax::Axis, offset) = Axis(map(x->reindex(x, offset), indexmap(ax)))
reindex(ax::ViewAxis, offset) = ViewAxis(viewindex(ax) .+ offset, indexmap(ax))

# Get AbstractAxis index
@inline Base.getindex(::AbstractAxis, idx::FlatIdx) = ComponentIndex(idx)
@inline Base.getindex(ax::AbstractAxis, ::Colon) = ComponentIndex(:, ax)
@inline Base.getindex(::AbstractAxis{IdxMap}, s::Symbol) where IdxMap =
    ComponentIndex(getproperty(IdxMap, s))