text:
module DynamicsUtils

using NQCDynamics:
    AbstractSimulation,
    Simulation,
    RingPolymerSimulation,
    Calculators,
    masses

divide_by_mass!(dv, masses) = dv ./= masses'
velocity!(dr, v, r, sim, t) = dr .= v

"""
    apply_interbead_coupling!(du::DynamicalVariables, u::DynamicalVariables,
                              sim::RingPolymerSimulation)

Applies the force that arises from the harmonic springs between adjacent beads.

Only applies the force for atoms labelled as quantum within the `RingPolymerParameters`.
"""
function apply_interbead_coupling!(dr::AbstractArray{T,3}, r::AbstractArray{T,3}, sim::RingPolymerSimulation) where {T}
    for i in axes(dr, 3)
        iplus = mod1(i+1, sim.beads.n_beads)
        iminus = mod1(i-1, sim.beads.n_beads)
        for j in sim.beads.quantum_atoms
            for k in axes(dr, 1)
                dr[k,j,i] -= sim.beads.ω_n² * (2r[k,j,i] - r[k,j,iplus] - r[k,j,iminus])
            end
        end
    end
    return nothing
end

function classical_hamiltonian end

function classical_kinetic_energy(sim::Simulation, v::AbstractMatrix)
    kinetic = zero(eltype(v))
    for i in axes(v, 2)
        for j in axes(v, 1)
            kinetic += masses(sim, i) * v[j,i]^2
        end
    end
    return kinetic / 2
end

function classical_kinetic_energy(sim::RingPolymerSimulation, v::AbstractArray{T,3}) where {T}
    kinetic = zero(eltype(v))
    for k in axes(v, 3)
        for i in axes(v, 2)
            for j in axes(v, 1)
                kinetic += masses(sim, i) * v[j,i,k]^2
            end
        end
    end
    return kinetic / 2
end

function classical_potential_energy(sim::Simulation, r::AbstractMatrix)
    Calculators.evaluate_potential!(sim.calculator, r)
    sim.calculator.potential
end

function classical_potential_energy(sim::RingPolymerSimulation, r::AbstractArray{T,3}) where {T}
    Calculators.evaluate_potential!(sim.calculator, r)
    sum(sim.calculator.potential) + RingPolymers.get_spring_energy(sim.beads, masses(sim), r)
end

include("dynamics_variables.jl")
include("callbacks.jl")
export CellBoundaryCallback
export TerminatingCallback

include("density_matrix_dynamics.jl")
include("plot.jl")

end # module


text:
using Oceananigans

include("ConvergenceTests/ConvergenceTests.jl")

using .ConvergenceTests.DoublyPeriodicTaylorGreen: setup_and_run

Nx = [8, 16, 32, 64, 128]
stop_time = 0.25

# Calculate time step
max_h = 2π / maximum(Nx)
Δt = 0.01 * max_h^2 # satisfy diffusive constraint for finest resolution
Nt = round(Int, stop_time / Δt)
Δt = stop_time / Nt

for N in Nx
    setup_and_run(Nx=N, Δt=Δt, stop_iteration=Nt)
end

ext: jl


text:
"""
    function innercorners(length::Int, width::Int)

return innercorners in a checkerboard of size length * width
"""
function innercorners(length::Int, width::Int)
    (length - 1) * (width - 1)
end

"""
    allboardcorners(length::Int, width::Int)

returns allboardcorners in a checkerboard of size length * width

"""
function allcorners(length::Int, width::Int)
    (length + 1) * (width + 1)
end

function drawdots!(img, res, color; size = 5)
    for i in res
        img[i[1]-size:i[1]+size i[2]-size:i[2]+size] .= color
    end
end

function draw_rect(img, dots, color = Gray{N0f8}(1); size = 5)
    for i in dots
        # img[i[1]-5:i[1]+5, i[2]-5:i[2]+5] .= color
        draw!(
            img,
            Polygon(
                RectanglePoints(
                    Point(i[2] - size, i[1] - size),
                    Point(i[2] + size, i[1] + size),
                ),
            ),
            color,
        )
    end
end

"""

"""
function kxkneighboardhood(
    chessboard,
    refined1;
    stdatol = 0.1,
    cortol = 0.6,
    n = 25,
    m = 11,
    k = 13,
)
    reut = zeros(Bool, length(refined1))
    board = Float32.(chessboard)
    for (idx, i) in enumerate(refined1)
        x, y = i[1], i[2]
        std1 = !isapprox(std(p1), std(p2); atol = stdatol)
        std2 = !isapprox(std(p3), std(p4); atol = stdatol)
        # cor1 = 10^((cor(vec(p1), vec(reverse(p2))) / 0.8) -1) > cortol
        # cor2 = 10^((cor(vec(p3), vec(reverse(p4))) / 0.8) -1) > cortol
        # if  std1 || std2 || !cor1 || !cor2
        #     continue
        # end
        imgtest = board[x-n:x+n, y-n:y+n]
        res = cor(vec(imgtest), vec(reverse(imgtest)))
        if std1 || std2 || res < 0.7
            continue
        end
        reut[idx] = true
    end
    refined2 = map((x, y) -> y ? x : nothing, refined1, reut)
    refined3 = filter(x -> x !== nothing, refined2)
end

"""
    nonmaxsuppresion(refined3)

returns checkboard refined points after non maximal suppresion, currently uses mean
"""
function nonmaxsuppresion(refined3)
    checked = Set([])
    final1 = []
    function dista(i, j)
        sqrt((i[1] - j[1])^2 + (i[2] - j[2])^2)
    end
    for (idxi, i) in enumerate(refined3)
        if idxi ∉ checked
            dist = []
            for (idxj, j) in enumerate(refined3)
                if idxj ∉ checked
                    dist1 = dista(i, j)
                    if dist1 < 10
                        push!(dist, idxj)
                    end
                end
            end
            local mat = zeros(1, 2)
            for i in dist
                i = refined3[i]
                mat = vcat([i[1] i[2]], mat)
            end
            res = Int64.(floor.(mean(mat[1:end-1, :]; dims = 1)))
            map(d -> push!(checked, d), dist)
            value = CartesianIndex(res[1], res[2])
            push!(final1, value)
        end
    end
    final1
end

"""
    markcorners(img::AbstractArray; method = harris, crthres = Percentile(99), LoGparams = 2.0.^[3.0], filter = (5,5), returnimg = true)

returns corners of checkerboard in an image with size length * width

### Arguments
- `img`: image to be processed
- `method`: method to be used for corner detection
- `crthres`: threshold for corner imcorner method
- `LoGparams`: parameters for LoG filter
- `filter`: size of filter for mapwindow
- `returnimg`: if true, returns image with corners marked

### Example
```jl

using CameraCalibrations
```
"""
function markcorners(
    img::AbstractArray;
    method = harris,
    crthres = Percentile(99),
    LoGparams = 2.0 .^ [3.0],
    filter = (5, 5),
    returnimg = false,
)
    imagecorners = imcorner(img, crthres; method = harris)
    img_cleaned = dilate(mapwindow(median!, (Gray.(imagecorners)), filter))
    results = blob_LoG(Int64.(img_cleaned), LoGparams)
    if returnimg == true
        resultantimage = zeros(size(img))
        map(x -> resultantimage[x.location] = 1, results)
        return map(x -> x.location, results), Gray.(resultantimage)
    else
        return map(x -> x.location, results)
    end
end

"""
    segboundariescheck(imgs; numcheck = 4)

returns if the boundaries of the images has different segments to detect different segments:

### Arguments

- `imgs`: images to be processed
- `numcheck`: number of segments in boundary to be detected

Retuns a bool array for of the images and if they satisfy the condition, we return true else false
### Example

We have a corner that looks something like below:

000111
000111
111000
111000

We want to check if the boundary is made up of 4 segments.
Boundary in this case would be : 000111100001111 starting from top edges and going clockwise.
In this we can say we have 4 changes in the boundary.
Numcheck is the number of changes we want to check.

"""
function segboundariescheck(imgs; numcheck = 4)
    check = zeros(Bool, length(imgs))
    for (idx, i) in enumerate(imgs)
        a = vcat(i[1, :], i[:, end], reverse(i[end, :]), reverse(i[:, 1]))
        numchange = 0
        for num = 2:length(a)
            if a[num] == 1 && a[num-1] == 0 || a[num] == 0 && a[num-1] == 1
                numchange = numchange + 1
            end
        end
        if numchange == numcheck
            check[idx] = true
        end
    end
    check
end

"""
    checkboundaries(checkerboard, cords; pixels = [11,23,35])

returns true if boundaries satisfy segboundariescheck for a range of pixels regions

### Arguments
- 'checkerboard': checkerboard img to be processed
- 'cords' : array of cartesian indices which indicaes corners in a image
- `pixels`: array of pixels region to be checked centered at cords
"""
function checkboundaries(checkerboard, cords; pixels = [11, 23, 35])
    currentstate = zeros(Bool, length(cords))
    # assumes that checkboard is gray
    checkerboard = checkerboard .> 0.4
    for n in pixels
        n = Int(floor((n - 1) / 2)) - 1
        res = map(x -> checkerboard[x[1]-n:x[1]+n, x[2]-n:x[2]+n], cords)
        # # res = map(x-> Gray.(x .> meanfinite(x)), corners)
        # res = map(x-> x .> 0.4, corners)
        check = segboundariescheck(res)
        currentstate = map(x -> (x > 0) ? true : false, currentstate .+ check)
    end
    refined = map((x, y) -> y ? x : nothing, cords, currentstate)
    refined1 = filter(x -> x !== nothing, refined)
end

"""

processes the checkerboard
"""
function process_image(chessboard)
    # we need a algorithm to check if there is a checkerboard or not in image
    # still need to study how filters from ImageFiltering.jl can improve results
    imagecorners = imcorner(chessboard, Percentile(99); method = Images.harris)
    # imagecorners = fastcorners(chessboard, 11, 0.20) # still gotta check if this is worth it
    imagecorners = clearborder(imagecorners, 35) # 35 is the boundary width we change
    results =
        map(x -> imagecorners[x] == true ? x : nothing, CartesianIndices(imagecorners))
    results = filter(x -> x !== nothing, results)
    correlationcheck = kxkneighboardhood(chessboard, results;)
    bounds = checkboundaries(chessboard, correlationcheck; pixels = [11, 23, 35])
    # also we need algorithm for checking if we have a board now, with connected components still
    # also we need algorithm for checking if we have outliers and remove them  if they exist
    finalcorners = nonmaxsuppresion(bounds) # return checkboard points
end

"""
    videotrack()

To work with realtime data after corners have been detected, use videotrack.
"""
# function videotrack()
#     if :GLMakie ∉ names(Main,imported=true)
#         throw(error("GLMakie needs to be imported first"))
#     end
#     try

#         img = read(cam)
#         fig = GLMakie.Figure(size = (1000, 700), title = "Checkerboard detection")
#         ax = GLMakie.Axis(
#             fig[1, 1],
#             aspect = DataAspect(),
#             xlabel = "x",
#             xlabelcolor = :black,
#             ylabel = "y label",
#             ylabelcolor = :white,
#             title = "Image",
#             backgroundcolor = :black,
#             labelcolor = :white,
#         )
#         node = Node(rotr90(img))
#         makieimg = image!(ax, node)
#         while isopen(cam)
#             read!(cam, img)
#             res, img = markcorners(img; returnimg = true)
#             node[] = rotr90(img)
#             if ispressed(scene, Keyboard.q) == true
#                 close(cam)
#                 return
#             end
#             sleep(1 / VideoIO.framerate(cam))
#         end
#         close(cam)
#     catch e
#         try
#             close(cam)
#         catch e
#             throw(error("Could not close camera camera"))
#         end
#         throw(error("Unable to open webcam"))
#     end
# end


text:
module Load

using LinearAlgebra
using SparseArrays
using HCubature

export NodalForce,NodalDisp,
BeamForce,QuadForce

abstract type AbstractElementForce end

struct NodalForce
    id::String
    val::Vector{Float64}
    csys::String #This is a string to indicate whether the value is given in global or local csys.
    NodalForce(id,force,csys="global")=new(id,force,csys)
end

struct NodalDisp
    id::String
    val::Vector{Float64}
    csys::String
    NodalDisp(id,disp,csys="global")=new(id,disp,csys)
end

mutable struct BeamForce  <: AbstractElementForce
    id::String
    f::Vector{Float64}
    s::Vector{Float64}
    σ₀::Vector{Float64}
    ϵ₀::Vector{Float64}
    BeamForce(id)=new(id,zeros(12),zeros(12),zeros(2),zeros(2))
end

mutable struct QuadForce  <: AbstractElementForce
    id::String
    f::Vector{Float64}
    s::Vector{Float64}
    σ₀::Vector{Float64}
    ϵ₀::Vector{Float64}
    QuadForce(id)=new(id,zeros(24,1),zeros(24,1),zeros(4,1),zeros(4,1))
end

end


text:
module Flatten

# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    # Use [`XX`](@ref) in the docs but not the readme
    replace(read(path, String), r"`(\w+\w)`" => s"[`\1`](@ref)")
    # Use doctests
    replace(read(path, String), "```julia" => "```jldoctest")
end Flatten

using ConstructionBase, FieldMetadata

import FieldMetadata: @flattenable, flattenable

export @flattenable, flattenable, flatten, reconstruct, update!, modify,
       metaflatten, fieldnameflatten, parentnameflatten, fieldtypeflatten, parenttypeflatten

struct EmptyIgnore end

# Default behaviour when no flattentrait/use/ignore args are given
const USE = Real
const IGNORE = EmptyIgnore
const FLATTENTRAIT = flattenable

# Generalised nested struct walker
nested(T::Type, P::Type, expr_builder, expr_combiner) =
    expr_combiner(T, [Expr(:..., expr_builder(T, fn)) for fn in fieldnames(T)])
# Unknown parent type is replaced with Nothing
nested(T::Type, expr_builder, expr_combiner) =
    nested(T, Nothing, expr_builder, expr_combiner)

"""
    flatten(obj, [flattentrait::Function], [use::Type], [ignore::Type])

Flattening. Flattens a nested struct or tuple to a flat tuple.
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments

- `obj`: The target type to be reconstructed
- `flattentrait`: A function returning a Bool, such as a FielMetadata method. With the form:
  `f(::Type, ::Type{Val{:fieldname}}) = true`
- `use`: Type or `Union` of types to return in the tuple.
- `ignore`: Types or `Union` of types  to ignore completly. These are not reurned or recursed over.

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> foo = Foo(1, 2, 3)
Foo{Int64,Int64,Int64}(1, 2, 3)

julia> flatten(foo)
(1, 2, 3)

julia> nested = Foo(Foo(1, 2, 3), 4.0, 5.0)
Foo{Foo{Int64,Int64,Int64},Float64,Float64}(Foo{Int64,Int64,Int64}(1, 2, 3), 4.0, 5.0)

julia> flatten(nested)
(1, 2, 3, 4.0, 5.0)

To convert the tuple to a vector, simply use [flatten(x)...], or
using static arrays to avoid allocations: `SVector(flatten(x))`.
```

Fields can be excluded from flattening with the `flattenable(struct, field)`
method. These are easily defined using `@flattenable` from
[FieldMetadata.jl](https://github.com/rafaqz/FieldMetadata.jl), or defining your own
custom function with FieldMetadata, or manually with the form:

```julia
julia> import Flatten: flattenable

julia> @flattenable struct Partial{A,B,C}
           a::A | true
           b::B | true
           c::C | false
       end

julia> nestedpartial = Partial(Partial(1.0, 2.0, 3.0), 4, 5)
Partial{Partial{Float64,Float64,Float64},Int64,Int64}(Partial{Float64,Float64,Float64}(1.0, 2.0, 3.0), 4, 5)

julia> nestedpartial = Partial(Partial(1.0, 2.0, 3.0), 4, 5)
Partial{Partial{Float64,Float64,Float64},Int64,Int64}(Partial{Float64,Float64,Float64}(1.0, 2.0, 3.0), 4, 5)

julia> flatten(nestedpartial)
(1.0, 2.0, 4)
```

"""
function flatten end

flatten_builder(T, fname) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        flatten(getfield(obj, $(QuoteNode(fname))), flattentrait, use, ignore)
    else
        ()
    end
end

flatten_combiner(T, expressions) = Expr(:tuple, expressions...)

flatten_inner(T) = nested(T, flatten_builder, flatten_combiner)

flatten(obj) = flatten(obj, flattenable)
flatten(obj, args...) = flatten(obj, flattenable, args...)
flatten(obj, ft::Function) = flatten(obj, ft, USE)
flatten(obj, ft::Function, use) = flatten(obj, ft, use, IGNORE)
flatten(x::X, ft::Function, use::Type{U}, ignore::Type{I}) where {X,U,I} =
    if X <: I
        ()
    elseif X <: U
        (x,)
    else
        _flatten(x, ft, use, ignore)
    end

@generated _flatten(obj, flattentrait, use, ignore) = flatten_inner(obj)

"""
    reconstruct(obj, data, [flattentrait::Function], [use::Type], [ignore::Type])

Reconstruct an object from Tuple or Vector data and an existing object.
Data should be at least as long as the queried fields in the obj.
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments
- `obj`: The target type to be reconstructed
- `data`: Replacement data - an `AbstractArray`, `Tuple` or type that defines `getfield`.
- `flattentrait`: A function returning a Bool, such as a FielMetadata method. With the form:
  `f(::Type, ::Type{Val{:fieldname}}) = true`
- `use`: Type or `Union` of types to return in the tuple.
- `ignore`: Types or `Union` of types  to ignore completly. These are not reurned or recursed over.

# Examples

```julia
julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> reconstruct(Foo(1, 2, 3), (1, :two, 3.0))
Foo{Int64,Symbol,Float64}(1, :two, 3.0)
```
"""
function reconstruct end

reconstruct_builder(T, fname) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        val, n = reconstruct(x, data, flattentrait, use, ignore, n)
        val
    else
        (getfield(obj, $(QuoteNode(fname))),)
    end
end

reconstruct_combiner(T, expressions) = :($(Expr(:tuple, expressions...)), n)

reconstruct_inner(::Type{T}) where T =
    nested(T, reconstruct_builder, reconstruct_combiner)

# Run from first data index and extract the final return value from the nested tuple
reconstruct(obj, data) = reconstruct(obj, data, flattenable)
reconstruct(obj, data, args...) = reconstruct(obj, data, flattenable, args...)
reconstruct(obj, data, ft::Function) = reconstruct(obj, data, ft, USE)
reconstruct(obj, data, ft::Function, use) = reconstruct(obj, data, ft, use, IGNORE)
# Need to extract the final return value from the nested tuple
reconstruct(obj, data, ft::Function, use, ignore) =
    reconstruct(obj, data, ft, use, ignore, firstindex(data))[1][1]
# Return value unmodified
reconstruct(x::I, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} = (x, ), n
# Return value from data. Increment position counter - the returned n + 1 becomes n
reconstruct(x::U, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} =
    (data[n],), n + 1
@generated reconstruct(obj, data, flattentrait::Function, use, ignore, n) =
    quote
        args, n = $(reconstruct_inner(obj))
        if length(args) > 0
            (constructorof(typeof(obj))(args...),), n
        else
            (obj,), n
        end
    end

"""
    modify(func, obj, args...)

Modify field in a type with a function
"""
modify(f, obj, args...) = reconstruct(obj, map(f, flatten(obj, args...)), args...)

"""
    update!(obj, data, [flattentrait::Function], [use::Type], [ignore::Type])

Update a mutable object with a `Tuple` or `Vector` of data.
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments

- `obj`: The target type to be reconstructed
- `data`: Replacement data - an `AbstractArray`, `Tuple` or type that defines `getfield`.
- `flattentrait`: A function returning a Bool, such as a FielMetadat method. With the form:
  `f(::Type, ::Type{Val{:fieldname}}) = true`
- `use`: Types to return in the tuple.
- `ignore`: Types to ignore completly. These are not reurned or recursed over.

# Examples

```jldoctest
julia> using Flatten

julia> mutable struct MutableFoo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> mufoo = MutableFoo(1, 2, 3)
MutableFoo{Int64,Int64,Int64}(1, 2, 3)

julia> update!(mufoo, (2, 4, 6))
MutableFoo{Int64,Int64,Int64}(2, 4, 6)

julia> mufoo = MutableFoo(1, 2, :three)
MutableFoo{Int64,Int64,Symbol}(1, 2, :three)

julia> update!(mufoo, (:foo,), Symbol)
MutableFoo{Int64,Int64,Symbol}(1, 2, :foo)
```
"""
function update! end

update_builder(T, fname) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        val, n = update!(x, data, flattentrait, use, ignore, n)
        setfield!(obj, $(QuoteNode(fname)), val[1])
    end
    ()
end

# Identical to the reconstruct builder for tuples, they're immutable
update_builder(T::Type{<:Tuple}, fname) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        val, n = update!(x, data, flattentrait, use, ignore, n)
        val
    else
        (getfield(obj, $(QuoteNode(fname))),)
    end
end

update_combiner(T, expressions) = :($(Expr(:tuple, expressions...)); ((obj,), n))
update_combiner(T::Type{<:Tuple}, expressions) = reconstruct_combiner(T, expressions)

update_inner(::Type{T}) where T =
    nested(T, update_builder, update_combiner)

update!(obj, data) = update!(obj, data, flattenable)
update!(obj, data, args...) = update!(obj, data, flattenable, args...)
update!(obj, data, ft::Function) = update!(obj, data, ft, USE)
update!(obj, data, ft::Function, use) = update!(obj, data, ft, use, IGNORE)
update!(obj, data, ft::Function, use, ignore) = begin
    update!(obj, data, ft, use, ignore, firstindex(data))[1][1]
    obj
end
update!(x::I, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} = (x,), n
update!(x::U, data, ft::Function, use::Type{U}, ignore::Type{I}, n) where {U,I} = (data[n],), n + 1
@generated update!(obj, data, flattentrait::Function, use, ignore, n) =
    quote
        x, n = $(update_inner(obj))
        (obj,), n
    end

"""
    metaflatten(obj, func, [flattentrait::Function], [use::Type], [ignore::Type])

Metadata flattening. Flattens data attached to a field using a passed in function
Query types and flatten trait arguments are optional, but you must pass `use` to pass `ignore`.

# Arguments

- `obj`: The target type to be reconstructed
- `func`: A function with the form: `f(::Type, ::Type{Val{:fieldname}}) = metadata`
- `flattentrait`: A function returning a Bool, such as a FielMetadata method. With the form:
  `f(::Type, ::Type{Val{:fieldname}}) = true`
- `use`: Type or `Union` of types to return in the tuple.
- `ignore`: Types or `Union` of types  to ignore completly. These are not reurned or recursed over.

We can flatten this `@foobar` metadata:

```jldoctest
julia> using Flatten, FieldMetadata

julia> import Flatten: flattenable

julia> @metadata foobar :foo;

julia> @foobar struct Foo{A,B,C}
           a::A | :bar
           b::B | :foobar
           c::C | :foo
       end;

julia> @foobar struct Bar{X,Y}
           x::X | :foobar
           y::Y | :bar
       end;

julia> metaflatten(Foo(1, 2, Bar(3, 4)), foobar)
(:bar, :foobar, :foobar, :bar)
```
"""
function metaflatten end

metaflatten_builder(T, fname) = quote
    if flattentrait($T, Val{$(QuoteNode(fname))})
        x = getfield(obj, $(QuoteNode(fname)))
        metaflatten(x, func, flattentrait, use, ignore, $T, Val{$(QuoteNode(fname))})
    else
        ()
    end
end

metaflatten_inner(T::Type) =
    nested(T, metaflatten_builder, flatten_combiner)

metaflatten(obj, func::Function) = metaflatten(obj, func, flattenable)
metaflatten(obj, func::Function, args...) = metaflatten(obj, func, flattenable, args...)
metaflatten(obj, func::Function, ft::Function) = metaflatten(obj, func, ft, USE)
metaflatten(obj, func::Function, ft::Function, use) = metaflatten(obj, func, ft, use, IGNORE)
metaflatten(obj, func::Function, ft::Function, use, ignore) =
    metaflatten(obj, func::Function, ft::Function, use, ignore, Nothing, Val{:none})

metaflatten(x::I, func::Function, ft::Function, use::Type{U}, ignore::Type{I}, P, fname) where {U,I} = ()
metaflatten(x::U, func::Function, ft::Function, use::Type{U}, ignore::Type{I}, P, fname) where {U,I} =
    (func(P, fname),)
@generated metaflatten(obj, func::Function, flattentrait::Function, use, ignore, P, fname) =
    metaflatten_inner(obj)

# Helper functions to get field data with metaflatten

fieldname_meta(T, ::Type{Val{N}}) where N = N

"""jldoctest
    fieldnameflatten(obj, args...)

Flatten the field names of an object. Args are passed to [`metaflatten`](@ref).

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> fieldnameflatten(Foo(1, 2, 3))
(:a, :b, :c)
```
"""
fieldnameflatten(obj, args...) = metaflatten(obj, fieldname_meta, args...)

fieldtype_meta(T, ::Type{Val{N}}) where N = fieldtype(T, N)

"""
    fieldtypeflatten(obj, args...)

Flatten the field types of an object. Args are passed to [`metaflatten`](@ref).

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> fieldtypeflatten(Foo(1.0, :two, "Three"), Union{Real,String})
(Float64, String)
```
"""
fieldtypeflatten(obj, args...) = metaflatten(obj, fieldtype_meta, args...)

fieldparent_meta(T, ::Type{Val{N}}) where N = T.name.name

"""
    parentnameflatten(obj, args...)

Flatten the name of the parent type of an object. Args are passed to [`metaflatten`](@ref).

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> struct Bar{X,Y}
           x::X
           y::Y
       end

julia> parentnameflatten(Foo(1, 2, Bar(3, 4)))
(:Foo, :Foo, :Bar, :Bar)
```
"""
parentnameflatten(obj, args...) = metaflatten(obj, fieldparent_meta, args...)

fieldparenttype_meta(T, ::Type{Val{N}}) where N = T

"""
    parenttypeflatten(obj, args...)

Flatten the parent type of an object. Args are passed to [`metaflatten`](@ref).

# Examples

```jldoctest
julia> using Flatten

julia> struct Foo{A,B,C}
           a::A
           b::B
           c::C
       end

julia> struct Bar{X,Y}
           x::X
           y::Y
       end

julia> parenttypeflatten(Foo(1, 2, Bar(3, 4)))
(Foo{Int64,Int64,Bar{Int64,Int64}}, Foo{Int64,Int64,Bar{Int64,Int64}}, Bar{Int64,Int64}, Bar{Int64,Int64})
```
"""
parenttypeflatten(obj, args...) = metaflatten(obj, fieldparenttype_meta, args...)

end # module


text:
# # Entries
#
# The key ingredient for data representations are `AxisEntries`.
#
# ## The `AxisEntries` type
#
# An `AxisEntries` object is
# made of four components:
# - axis,
# - entries.

using SplitApplyPlot, CairoMakie
using SplitApplyPlot: CategoricalScale, ContinuousScale
resolution = (600, 600)
fig = Figure(; resolution)
N = 11
rg = range(1, 2, length=N)
ae = AxisEntries(
    Axis(fig[1, 1]),
    [
        Entry(
            Scatter,
            arguments(rg, cosh.(rg), color=1:N, marker=fill("b", N));
            markersize = 15
        ),
        Entry(
            Scatter,
            arguments(rg, sinh.(rg), color=1:N, marker=fill("c", N));
            markersize = 15
        ),
    ],
    arguments(
        ContinuousScale(identity, (0, 4)),
        ContinuousScale(identity, (0, 4)),
        color=ContinuousScale(identity, (1, N)),
        marker=CategoricalScale(["a", "b", "c"], [:circle, :utriangle, :dtriangle]), #scales
    ),
    arguments("x", "y", color="identity", marker="function"), #labels
)
plot!(ae)
fig


text:
using Dates, Tables
using DataStructures: OrderedDict

include("testfiles/testfiles.jl")

@testset "CSV.File" begin

for test in testfiles
    testfile(test...)
end

#test on non-existent file
@test_throws ArgumentError CSV.File("");

#test where datarow > headerrow
@test_throws ArgumentError CSV.File(joinpath(dir, "test_no_header.csv"); datarow=1, header=2);

#test bad types
f = CSV.File(joinpath(dir, "test_float_in_int_column.csv"); types=[Int, Int, Int], strict=true)
@test_throws CSV.Error (f |> columntable)

# Integer overflow; #100
@test_throws CSV.Error (CSV.File(joinpath(dir, "int8_overflow.csv"); types=[Int8], strict=true) |> columntable)

# #137
tbl = (a=[11,22], dt=[Date(2017,12,7), Date(2017,12,14)], dttm=[DateTime(2017,12,7), DateTime(2017,12,14)])
io = IOBuffer()
tbl |> CSV.write(io; delim='\t')
seekstart(io)
f = CSV.File(io; delim='\t', allowmissing=:auto)
@test (f |> columntable) == tbl

# #172
@test_throws ArgumentError CSV.File(joinpath(dir, "test_newline_line_endings.csv"), types=Dict(1=>Integer))

# #289
tmp = CSV.File(IOBuffer(" \"a, b\", \"c\" "), datarow=1) |> DataFrame
@test size(tmp) == (1, 2)
@test tmp.Column1[1] == "a, b"
@test tmp.Column2[1] == "c"

tmp = CSV.File(IOBuffer(" \"2018-01-01\", \"1\" ,1,2,3"), datarow=1) |> DataFrame
@test size(tmp) == (1, 5)
@test tmp.Column1[1] == Date(2018, 1, 1)
@test tmp.Column2[1] == 1
@test tmp.Column3[1] == 1
@test tmp.Column4[1] == 2
@test tmp.Column5[1] == 3

# #329
df = CSV.read(joinpath(dir, "test_types.csv"), types=Dict(:string=>Union{Missing,DateTime}), silencewarnings=true)
@test df.string[1] === missing

# #352
@test_throws ArgumentError first(CSV.File(joinpath(dir, "test_types.csv"))).a

# @time f = CSV.File(joinpath(dir, "pandas_zeros.csv"), allowmissing=:none);
# @time t = f |> columntable;
# @time t = Tables.buildcolumns(nothing, Tables.rows(f));
# using Profile
# Profile.clear()
# @profile f |> columntable;
# Profile.print(C=true)
# Profile.print()

# Test that using an OrderedDict for types doesn't error #373
df = CSV.read(joinpath(dir, "test_types.csv"), types=OrderedDict(:string=>Union{Missing,DateTime}))
@test df.string[1] === missing

end # testset


text:
## higher-order functions ##

import Base: +, -, div, mod, rem, &, |, xor
import SparseArrays: nnz

Base.map(f, d0::DArray, ds::AbstractArray...) = broadcast(f, d0, ds...)

function Base.map!(f::F, dest::DArray, src::DArray{<:Any,<:Any,A}) where {F,A}
    asyncmap(procs(dest)) do p
        remotecall_fetch(p) do
            map!(f, localpart(dest), makelocal(src, localindices(dest)...))
            return nothing
        end
    end
    return dest
end

function Base.reduce(f, d::DArray)
    results = asyncmap(procs(d)) do p
        remotecall_fetch(p) do
            return reduce(f, localpart(d))
        end
    end
    reduce(f, results)
end

function Base._mapreduce(f, op, ::IndexCartesian, d::DArray)
    results = asyncmap(procs(d)) do p
        remotecall_fetch((_f,_op,_d)->mapreduce(_f, _op, localpart(_d)), p, f, op, d)
    end

    reduce(op, results)
end
Base._mapreduce(f, op, ::IndexCartesian, d::SubDArray) = Base._mapreduce(f, op, IndexCartesian(), DArray(d))
# Base.mapreduce(f, opt::Union{typeof(|), typeof(&)}, d::DArray) = _mapreduce(f, opt, d)
# Base.mapreduce(f, opt::Function, d::DArray) = _mapreduce(f, opt, d)
# Base.mapreduce(f, opt, d::DArray) = _mapreduce(f, opt, d)

# mapreducedim
function Base.reducedim_initarray(A::DArray, region, v0, ::Type{R}) where {R}
    # Store reduction on lowest pids
    pids = A.pids[ntuple(i -> i in region ? (1:1) : (:), ndims(A))...]
    chunks = similar(pids, Future)
    @sync for i in eachindex(pids)
        @async chunks[i...] = remotecall_wait(() -> Base.reducedim_initarray(localpart(A), region, v0, R), pids[i...])
    end
    return DArray(chunks)
end
Base.reducedim_initarray(A::DArray, region, v0::T) where {T} = Base.reducedim_initarray(A, region, v0, T)

# Compute mapreducedim of each localpart and store the result in a new DArray
function mapreducedim_within(f, op, A::DArray, region)
    arraysize = [size(A)...]
    gridsize = [size(A.indices)...]
    arraysize[[region...]] = gridsize[[region...]]
    indx = similar(A.indices)

    for i in CartesianIndices(indx)
        indx[i] = ntuple(j -> j in region ? (i.I[j]:i.I[j]) : A.indices[i][j], ndims(A))
    end
    cuts = [i in region ? collect(1:arraysize[i] + 1) : A.cuts[i] for i in 1:ndims(A)]
    return DArray(next_did(), I -> mapreduce(f, op, localpart(A), dims=region),
        tuple(arraysize...), procs(A), indx, cuts)
end

# Compute mapreducedim accros the processes. This should be done after mapreducedim
# has been run on each localpart with mapreducedim_within. Eventually, we might
# want to write mapreducedim_between! as a binary reduction.
function mapreducedim_between!(f, op, R::DArray, A::DArray, region)
    asyncmap(procs(R)) do p
        remotecall_fetch(p, f, op, R, A, region) do f, op, R, A, region
            localind = [r for r = localindices(A)]
            localind[[region...]] = [1:n for n = size(A)[[region...]]]
            B = convert(Array, A[localind...])
            Base.mapreducedim!(f, op, localpart(R), B)
            nothing
        end
    end
    return R
end

function Base.mapreducedim!(f, op, R::DArray, A::DArray)
    lsize = Base.check_reducedims(R,A)
    if isempty(A)
        return copy(R)
    end
    region = tuple(collect(1:ndims(A))[[size(R)...] .!= [size(A)...]]...)
    if isempty(region)
        return copyto!(R, A)
    end
    B = mapreducedim_within(f, op, A, region)
    return mapreducedim_between!(identity, op, R, B, region)
end

## Some special cases
function Base._all(f, A::DArray, ::Colon)
    B = asyncmap(procs(A)) do p
        remotecall_fetch(p) do
            all(f, localpart(A))
        end
    end
    return all(B)
end

function Base._any(f, A::DArray, ::Colon)
    B = asyncmap(procs(A)) do p
        remotecall_fetch(p) do
            any(f, localpart(A))
        end
    end
    return any(B)
end

function Base.count(f, A::DArray)
    B = asyncmap(procs(A)) do p
        remotecall_fetch(p) do
            count(f, localpart(A))
        end
    end
    return sum(B)
end

function nnz(A::DArray)
    B = asyncmap(A.pids) do p
        remotecall_fetch(nnz∘localpart, p, A)
    end
    return reduce(+, B)
end

function Base.extrema(d::DArray)
    r = asyncmap(procs(d)) do p
        remotecall_fetch(p) do
            extrema(localpart(d))
        end
    end
    return reduce((t,s) -> (min(t[1], s[1]), max(t[2], s[2])), r)
end

Statistics._mean(A::DArray, region) = sum(A, dims = region) ./ prod((size(A, i) for i in region))

# Unary vector functions
(-)(D::DArray) = map(-, D)

map_localparts(f::Callable, d::DArray) = DArray(i->f(localpart(d)), d)
map_localparts(f::Callable, d1::DArray, d2::DArray) = DArray(d1) do I
    f(localpart(d1), localpart(d2))
end

function map_localparts(f::Callable, DA::DArray, A::Array)
    s = verified_destination_serializer(procs(DA), size(DA.indices)) do pididx
        A[DA.indices[pididx]...]
    end
    DArray(DA) do I
        f(localpart(DA), localpart(s))
    end
end

function map_localparts(f::Callable, A::Array, DA::DArray)
    s = verified_destination_serializer(procs(DA), size(DA.indices)) do pididx
        A[DA.indices[pididx]...]
    end
    DArray(DA) do I
        f(localpart(s), localpart(DA))
    end
end

function map_localparts!(f::Callable, d::DArray)
    asyncmap(procs(d)) do p
        remotecall_fetch((f,d)->(f(localpart(d)); nothing), p, f, d)
    end
    return d
end

# Here we assume all the DArrays have
# the same size and distribution
map_localparts(f::Callable, As::DArray...) = DArray(I->f(map(localpart, As)...), As[1])

function samedist(A::DArray, B::DArray)
    (size(A) == size(B)) || throw(DimensionMismatch())
    if (procs(A) != procs(B)) || (A.cuts != B.cuts)
        B = DArray(x->B[x...], A)
    end
    B
end

for f in (:+, :-, :div, :mod, :rem, :&, :|, :xor)
    @eval begin
        function ($f)(A::DArray{T}, B::DArray{T}) where T
            B = samedist(A, B)
            map_localparts($f, A, B)
        end
        ($f)(A::DArray{T}, B::Array{T}) where {T} = map_localparts($f, A, B)
        ($f)(A::Array{T}, B::DArray{T}) where {T} = map_localparts($f, A, B)
    end
end

function Base.mapslices(f, D::DArray{T,N,A}; dims) where {T,N,A}
    if !(dims isa AbstractVector)
        dims = [dims...]
    end
    if !all(t -> t == 1, size(D.indices)[dims])
        p = ones(Int, ndims(D))
        nondims = filter(t -> !(t in dims), 1:ndims(D))
        p[nondims] = defaultdist([size(D)...][[nondims...]], procs(D))
        DD = DArray(size(D), procs(D), p) do I
            return convert(A, D[I...])
        end
        return mapslices(f, DD, dims=dims)
    end

    refs = Future[remotecall((x,y,z)->mapslices(x,localpart(y),dims=z), p, f, D, dims) for p in procs(D)]

    DArray(reshape(refs, size(procs(D))))
end

function _ppeval(f, A...; dim = map(ndims, A))
    if length(dim) != length(A)
        throw(ArgumentError("dim argument has wrong length. length(dim) = $(length(dim)) but should be $(length(A))"))
    end
    narg = length(A)
    dimlength = size(A[1], dim[1])
    for i = 2:narg
        if dim[i] > 0 && dimlength != size(A[i], dim[i])
            throw(ArgumentError("lengths of broadcast dimensions must be the same. size(A[1], $(dim[1])) = $dimlength but size(A[$i], $(dim[i])) = $(size(A[i], dim[i]))"))
        end
    end
    dims = []
    idx  = []
    args = []
    for i = 1:narg
        push!(dims, ndims(A[i]))
        push!(idx, Any[Colon() for d in 1:dims[i]])
        if dim[i] > 0
            idx[i][dim[i]] = 1
            push!(args, view(A[i], idx[i]...))
        else
            push!(args, A[i])
        end
    end
    R1 = f(args...)
    ridx = Any[1:size(R1, d) for d in 1:ndims(R1)]
    push!(ridx, 1)
    Rsize = map(last, ridx)
    Rsize[end] = dimlength
    R = Array{eltype(R1)}(undef, Rsize...)

    for i = 1:dimlength
        for j = 1:narg
            if dim[j] > 0
                idx[j][dim[j]] = i
                args[j] = view(A[j], idx[j]...)
            else
                args[j] = A[j]
            end
        end
        ridx[end] = i
        R[ridx...] = f(args...)
    end

    return R
end

"""
     ppeval(f, D...; dim::NTuple)

Evaluates the callable argument `f` on slices of the elements of the `D` tuple.

#### Arguments
`f` can be any callable object that accepts sliced or broadcasted elements of `D`.
The result returned from `f` must be either an array or a scalar.

`D` has any number of elements and the elements can have any type. If an element
of `D` is a distributed array along the dimension specified by `dim`. If an
element of `D` is not distributed, the element is by default broadcasted and
applied on all evaluations of `f`.

`dim` is a tuple of integers specifying the dimension over which the elements
of `D` is slices. The length of the tuple must therefore be the same as the
number of arguments `D`. By default distributed arrays are slides along the
last dimension. If the value is less than or equal to zero the element are
broadcasted to all evaluations of `f`.

#### Result
`ppeval` returns a distributed array of dimension `p+1` where the first `p`
sizes correspond to the sizes of return values of `f`. The last dimension of
the return array from `ppeval` has the same length as the dimension over which
the input arrays are sliced.

#### Examples
```jl
addprocs(Sys.CPU_THREADS)

using DistributedArrays

A = drandn((10, 10, Sys.CPU_THREADS), workers(), [1, 1, Sys.CPU_THREADS])

ppeval(eigvals, A)

ppeval(eigvals, A, randn(10,10)) # broadcasting second argument

B = drandn((10, Sys.CPU_THREADS), workers(), [1, Sys.CPU_THREADS])

ppeval(*, A, B)
```
"""
function ppeval(f, D...; dim::NTuple = map(t -> isa(t, DArray) ? ndims(t) : 0, D))
    #Ensure that the complete DArray is available on the specified dims on all processors
    for i = 1:length(D)
        if isa(D[i], DArray)
            for idxs in D[i].indices
                for d in setdiff(1:ndims(D[i]), dim[i])
                    if length(idxs[d]) != size(D[i], d)
                        throw(DimensionMismatch(string("dimension $d is distributed. ",
                            "ppeval requires dimension $d to be completely available on all processors.")))
                    end
                end
            end
        end
    end

    refs = Future[remotecall((x, y, z) -> _ppeval(x, map(localpart, y)...; dim = z), p, f, D, dim) for p in procs(D[1])]

    # The array of Futures has to be reshaped for the DArray constructor to work correctly.
    # This requires a fetch and the DArray is also fetching so it might be better to modify
    # the DArray constructor.
    sd = [size(D[1].pids)...]
    nd = remotecall_fetch((r)->ndims(fetch(r)), refs[1].where, refs[1])
    DArray(reshape(refs, tuple([sd[1:nd - 1]; sd[end]]...)))
end


text:
# This file illustrates the calculation of power spectral density and
# total average power of a sinusoidal signal.

using FFTW
using Plots

# Construct the signal
f = 20
fs = 100
l = 2^10
ts = 1 / fs
t = (0 : l -1) * ts
x = cos.(2π * f * t)

# Compute fft
xf = fft(x) * ts
ff = 0 : fs / (l - 1) : fs

# Compute power
psd = abs.(xf).^2 / (l - 1) / ts
Δf = fs / (l - 1)
Pf = sum(psd) * Δf              # Average power calculated in time
Pt = sum(abs.(x).^2) / (l - 1)  # Average power calculated in frequency
@show Pf, Pt

# Plots
plt = plot(layout=1)
plot!(ff, psd, subplot=1)
xlabel!("Frequency [Hz]")
ylabel!("Ψ(f)")
title!("Power Spectal Density")
display(plt)


text:
## header to provide surface and context
using Cairo
using Colors

c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);

z = Array{RGB24}(2,2)
c1 = convert(RGB24,colorant"grey20")
c2 = convert(RGB24,colorant"grey80")
z[1,1] = c1
z[1,2] = c2
z[2,1] = c2
z[2,2] = c1

img = CairoImageSurface(z)
pattern = CairoPattern(img);
pattern_set_extend(pattern, Cairo.EXTEND_REPEAT);
pattern_set_filter(pattern, Cairo.FILTER_BILINEAR);

m = CairoMatrix(1/8.0,0,0,1/8.0,0,0);

set_matrix(pattern, m);
set_source(cr, pattern);
paint(cr);
restore(cr);

save(cr);

# 5 uses of set_source

# color
save(cr)
rectangle(cr,16,32,224,32)
set_source(cr, colorant"red4")
fill(cr)
restore(cr)

# color with alpha
save(cr)
rectangle(cr,16,72,224,32)
set_source(cr, alphacolor(colorant"blue",0.5))
fill(cr)
restore(cr)

# image from surface
save(cr)
rectangle(cr,16,112,224,32)
clip(cr)
new_path(cr)
s = read_from_png("data/mulberry.png")
set_source(cr,s)
paint_with_alpha(cr,0.6)
restore(cr)

# image from context (from surface)
save(cr)
rectangle(cr,16,152,224,32)
clip(cr)
new_path(cr)
set_source(cr,creategc(s))
paint(cr)
restore(cr)

# image for surface with offset
save(cr)
rectangle(cr,16,192,224,32)
clip(cr)
set_source(cr,s,0.0,40.0)
paint(cr)
restore(cr)

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_alpha_paint.png");


text:
# Pkg.clone("https://github.com/lejon/TSne.jl.git")
# Pkg.add("Colors")
# Pkg.add("MNIST")
# Pkg.add("HDF5")

using Forge.DSL
using TSne
using Colors
using MNIST
# using HDF5

const N = 2500

function rescale(A, dim::Integer=1)
    res = A .- mean(A, dim)
    res ./= map(x -> x > 0.0 ? x : 1.0, std(A, dim))
    res
end
data, labels = traindata()
data = convert(Matrix{Float64}, data[:, 1:N])'
X = rescale(data, 1)
Y = tsne(X, 2, 50, 3000, 20.0)

# f = h5open("tsne.h5", "w")
# write(f, "Y", Y)
# write(f, "labels", labels)
# close(f)
#
# f = h5open("tsne.h5")
# Y = read(f, "Y")
# labels = read(f, "labels")
# close(f)

window = create_window(1000, 800, "t-SNE Demo")
show_window(window)

chart = create_chart_2d()
set_axes_limits(chart, -70, 70, -80, 80)

plots = Dict()
[setindex!(plots, [], n) for n in 0:9]
for idx in 1:N
    (x, y) = Y[idx, :]
    label = Int(labels[idx])
    append!(plots[label], [x y;])
end

colors = distinguishable_colors(10)

for label in 0:9
    pos = plots[label]
    data = reshape(pos, Int(length(pos)/2), 2)
    plot = create_plot(chart, length(data)/2, Forge.FG_PLOT_SCATTER, Forge.FG_MARKER_CIRCLE)
    set_color(plot, colors[label+1])
    set_legend(plot, string(label))
    update_vertex_buffer(plot, data)
end

draw_chart(!isclosed, window, chart)
save_framebuffer(window, "tsne.png")


text:
function hpipm_data()
    py_qp_data = hp.hpipm_data()
    ju_qp_data = py2ju_obj(py_qp_data)
    return ju_qp_data
end


text:
#Demo 4: Multi-dimensional datasets
#-------------------------------------------------------------------------------

import Statistics
using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))

#==Input data
===============================================================================#
tbit = 1e-9 #Bit period
osr = 20 #samples per bit
nbit_Π = 5 #Π-pulse length, in number of bits
nsamples = 20

#==Computations
===============================================================================#
seq = 1.0*prbs(reglen=5, seed=1, nsamples=nsamples)
tΠ = DataF1(0:(tbit/osr):(nbit_Π*tbit)) #Time of a single pulse

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("amp", [.8, 1, 1.2])
	PSweep("tau", (tbit/5).*[.5, 1, 2])
	PSweep("offset", [-1, 0, 1])
]

#Generate data:
#-------------------------------------------------------------------------------
#NOTE: verbose... trying to test DataHR & "parameter" operations.

#Generate single pulse for all parametric corners:
Π = fill(DataHR, sweeplist) do amp, tau, offset
	return pulse(tΠ, Pole(1/tau,:rad), tpw=tbit)
end
amp = parameter(Π, "amp")
offset = parameter(Π, "offset")
pat = (pattern(seq, Π, tbit=tbit)-0.5)*amp + offset #Center data pattern
t = xval(pat)

#Do some DataHR-level transformations:
result = pat
tmax = tbit*(nsamples+nbit_Π-1)
skew = DataF1([0, tmax],[0, 0.5])
result = result + result + 4 + skew
result += Statistics.mean(result)

#==Generate plot
===============================================================================#
plot = push!(cons(:plot, vvst, title="PRBS Pattern", legend=false),
	cons(:wfrm, result, label="pat"),
)

pcoll = push!(cons(:plotcoll, title="Mulit-Dataset Tests"), plot)
	pcoll.ncolumns = 1

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
# DomainSets.jl
module DomainSets

# We use static vectors internally

using StaticArrays
using Base

using LinearAlgebra, Statistics
import LinearAlgebra: cross, ×, pinv

using IntervalSets

################################
## Exhaustive list of imports
################################

# Generated functions
import Base: @ncall

# Operator symbols
import Base: *, +, -, /, \, ^
import Base: |, &
import Base: ∪, ∩
import Base: ==, isapprox
import Base: ∘

# Set operations
import Base: intersect, union, setdiff, in, isempty, minimum, maximum,
    issubset

# Arrays
import Base: size, length, ndims, getindex, eltype, ndims, hash
import Base: inv
import Base: isreal
import Base: zero

# Types, promotions and conversions
import Base: convert, widen, promote_rule

# Display
import Base: show

# Various
import Base: Bool

# IntervalSets
import IntervalSets: (..), endpoints, Domain, AbstractInterval, TypedEndpointsInterval,
                        leftendpoint, rightendpoint, isleftopen, isleftclosed,
                        isrightopen, isrightclosed, isopen, isclosed,
                        infimum, supremum
export ..

################################
## Exhaustive list of exports
################################

## Utils

# from util/tensorproducts.jl
export cartesianproduct

## Spaces

# from spaces/space.jl
export spaceof, superspaceof, spacetype, origin, superspace, issubspace, subeltype
export isomorphic, ≅, embedded, ↪
# from spaces/basic_spaces.jl
export ℤ, ℚ, ℝ, ℂ, ℝ1, ℝ2, ℝ3, ℝ4

## Maps

# from maps/maps.jl
export AbstractMap, applymap, jacobian, linearize
export domaintype, codomaintype
export left_inverse, right_inverse
export apply_inverse, apply_left_inverse, apply_right_inverse
export image
# from maps/composite_map.jl
export CompositeMap
export ∘
# from maps/productmap.jl
export ProductMap
export tensorproduct
# from maps/basic_maps.jl
export IdentityMap, ConstantMap
# from maps/embedding_map.jl
export embedding_map, restriction_map, isomorphism_map
# from maps/affine_maps.jl
export AffineMap, Translation, LinearMap
export islinear
export linear_map, interval_map, scaling_map, rotation_map, translation_map
# from maps/coordinates.jl
export CartToPolarMap, PolarToCartMap

## Generic domains

# from generic/domain.jl
export Domain, EuclideanDomain
export dimension
export approx_in
export isclosed, iscompact
export boundary, ∂
export point_in_domain

# from generic/derived_domain.jl
export DerivedDomain
export superdomain

# from generic/productdomain.jl
export ProductDomain

# from generic/setoperations.jl
export UnionDomain, IntersectionDomain, DifferenceDomain
export TranslatedDomain

# from generic/mapped_domain.jl
export source, target

# from generic/arithmetics.jl
export rotate

export infimum, supremum

## Specific domains

# from domains/simple.jl
export EmptySpace, FullSpace
# from domains/interval.jl
export AbstractInterval, Interval, UnitInterval, ChebyshevInterval, OpenInterval, ClosedInterval
export leftendpoint, rightendpoint, isleftopen, isrightopen
export cardinality
export HalfLine, NegativeHalfLine
# from domains/simple.jl
export UnitHyperBall, UnitDisk, UnitBall, Cube, Simplex, UnitSimplex, UnitHyperSphere, Point
export cylinder
# from domains/circle.jl
export UnitCircle, UnitSphere
export ellipse, ellipse_shape
export parameterization, gradient

include("util/common.jl")
include("util/products.jl")

include("spaces/space.jl")
include("spaces/space_promotions.jl")
include("spaces/basic_spaces.jl")
include("spaces/productspace.jl")

include("maps/maps.jl")
include("maps/composite_map.jl")
include("maps/productmap.jl")
include("maps/basic_maps.jl")
include("maps/embedding_map.jl")
include("maps/affine_map.jl")
include("maps/coordinates.jl")

include("generic/domain.jl")
include("generic/lazy.jl")
include("generic/derived_domain.jl")
include("generic/productdomain.jl")
include("generic/setoperations.jl")
include("generic/mapped_domain.jl")
include("generic/promotion.jl")
include("generic/arithmetics.jl")

include("domains/trivial.jl")
include("domains/interval.jl")
include("domains/simple.jl")
include("domains/circle.jl")

end # module


text:
using Meshing
using GeometryTypes
using LinearAlgebra: dot, norm
using FileIO

sphere(v) = sqrt(sum(v.^2)) - 1

sdf = SignedDistanceField(sphere, HyperRectangle(Vec(-1,-1,-1.), Vec(2,2,2.)))

mesh_type = HomogenousMesh{Point{3,Float64}, Face{3,Int}}

mc = mesh_type(sdf, MarchingCubes())
mt = mesh_type(sdf, MarchingTetrahedra())
ns = mesh_type(sdf, NaiveSurfaceNets())

mcf = mesh_type(sphere, sdf.bounds, MarchingCubes(), samples=size(sdf))
mtf = mesh_type(sphere, sdf.bounds, MarchingTetrahedra(), samples=size(sdf))
nsf = mesh_type(sphere, sdf.bounds, NaiveSurfaceNets(), samples=size(sdf))

# shift
@assert length(mc.vertices) == length(mcf.vertices)
@assert length(mt.vertices) == length(mtf.vertices)
@assert length(ns.vertices) == length(nsf.vertices)

for i in eachindex(mt.vertices)
    mt.vertices[i] = mt.vertices[i] .+ Point(3,0,0)
    mtf.vertices[i] = mtf.vertices[i] .+ Point(3,3,0)
end
for i in eachindex(mc.vertices)
    mc.vertices[i] = mc.vertices[i] .+ Point(-3,0,0)
    mcf.vertices[i] = mcf.vertices[i] .+ Point(-3,3,0)
end
for i in eachindex(ns.vertices)
    #ns.vertices[i] = ns.vertices[i] .+ Point(0,0,0)
    nsf.vertices[i] = nsf.vertices[i] .+ Point(0,3,0)
end

mesh = merge(mc,mt,ns,mcf,mtf,nsf)

# save all in one PLY file
save("output_comparison.ply",mesh)


text:
using ApproxFun, LinearAlgebra, DualNumbers
import ApproxFun: bandwidth

#
# Consider the following Sturm--Liouville problem with polynomial p, q, and w.
#
# [(-𝒟)(p𝒟) + q] u = λ w u,  (u±bu')(±1) = 0.
#
# For which value of b, the Neumann part of the Robin boundary conditions,
# is λ₀ equal to 1?
#
# To answer this question, we use a dual number to observe the sensitivity in λ₀
# with respect to perturbations in b. The `realpart` of the objective function
#
#    f(λ(b)) = λ(b) - 1,
#
# extracts the difference, whereas the `dualpart` evaluates to ∂f/∂b (λ(b)).
#

n = 100
λ = dual(0.0, 0.0)
b = dual(0.0, 1.0)
v = zeros(Dual{Float64}, n)
v[1] = 1
v .+= 0.003.*randn(Float64, n)
for _ in 1:8
    d = Segment(-1..1)
    S = Ultraspherical(0.5, d)
    NS = NormalizedPolynomialSpace(S)
    D = Derivative(S)
    p = Fun(x->1+x^2, S)
    q = Fun(x->x^3, S)
    L = -D*(p*D)+q
    w = Fun(x->1+x/4, S)
    M = Multiplication(w, S)
    C = Conversion(domainspace(L), rangespace(L))
    B = [Evaluation(S, -1, 0) - b*Evaluation(S, -1, 1); Evaluation(S, 1, 0) + b*Evaluation(S, 1, 1)]
    QS = QuotientSpace(B)
    Q = Conversion(QS, S)
    D1 = Conversion(S, NS)
    D2 = Conversion(NS, S)
    R = D1*Q
    P = cache(PartialInverseOperator(C, (0, bandwidth(L, 1) + bandwidth(R, 1) + bandwidth(C, 2))))
    A = R'D1*P*L*D2*R
    B = R'D1*M*D2*R
    SA = Symmetric(A[1:n,1:n], :L)
    SB = Symmetric(B[1:n,1:n], :L)
    for _ = 1:7
        global λ = dot(v, SA*v)/dot(v, SB*v)
        ldiv!(ldlt!(SA-λ*SB), v)
        global v = v/norm(v)
        println("\tThis is λ: ", λ)
    end
    global b -= (realpart(λ)-1)/dualpart(λ)
    println("This is b: ", b)
end

max_stars_repo_path: Extras/Inverse eigenvalue problem.jl

max_stars_repo_name: Gyoshi/ApproxFunExamples

max_issues_repo_path: Extras/Inverse eigenvalue problem.jl

max_issues_repo_name: Gyoshi/ApproxFunExamples

max_forks_repo_path: Extras/Inverse eigenvalue problem.jl

max_forks_repo_name: Gyoshi/ApproxFunExamples


text:
function GammaMnsFunction(arg0::jint, arg1::jdouble, arg2::jint)
    return GammaMnsFunction((jint, jdouble, jint), arg0, arg1, arg2)
end

function get_derivative(obj::GammaMnsFunction, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getDerivative", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function get_value(obj::GammaMnsFunction, arg0::jint, arg1::jint, arg2::jint)
    return jcall(obj, "getValue", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end


text:
module Aerospace

## Structures
export
    ## Constants
    go,
    d2r,
    r2d,
    REARTH,
    WEII3,
    GO,
    AGRAV,
    GCONST,
    EARTH_MASS,
    MU,
    C20,
    FLATTENING,
    SMAJOR_AXIS,
    GW_CLONG,
    RGAS,
    KBOLTZ,
    R2D,
    D2R,
    PI,
    EPS,
    SMALL,
    ILARGE,
    LARGE,
    M2F,
    M2NMI,
    J2,
    wgs_aa,
    wgs_finv,
    wgs_ff,
    wgs_bb,
    wgs_e2,

    ## Structures
    Aero6DOF,
    Aero5DOF,
    QuadControl,
    WGS_Pos,
    Velocity_States,
    Euler_Angles,
    MassProps,
    Angle_Rates,
    Quaternions,
    Motors,
    NavStates,
    Actuator,

    ## Transformations
    TR1,
    TR2,
    TR3,
    TR_BN,
    TR_NE,
    TR_EI,
    TR_WB,
    TR_VN,
    TR_BV,

    # Atmosphere
    Atmosphere,

    ## WGS-84
    ECEF2LLH,
    LLH2ECEF,
    WGS84_GRAVITY,

    ## Quaternion
    QuatInit,
    QuatInit!,
    QuatDerivative,
    QuatToDCM,

    ## Functions
    Pause,
    CloseAll,
    Bound,

    ## Euler Angle
    ComputeEuler,

    ## Actuator Functions
    Initialize_Act,
    Update_Act,

    ## Airframe Functions
    InitPos,
    InitVel,
    InitEuler,
    InitAngleRates,
    Mass_Props,
    Airframe_Initialize,
    Airframe5_Initialize,
    Airframe_Update,
    Airframe5_Update,

    ## Interpolation Functions
    interp1,
    interp2,
    interp3

mutable struct WGS_Pos
    time::Float64
    latitude::Float64
    longitude::Float64
    altitude::Float64
    rm_ecef::Array{Float64,1}
    rm_eci::Array{Float64,1}
    Tr_ei::Array{Float64,2}
    Tr_ne::Array{Float64,2}
end

mutable struct Velocity_States
    time_sim::Float64
    vmag::Float64
    gamma::Float64
    azimuth::Float64
    vm_ned::Array{Float64,1}
    vm_ecef::Array{Float64,1}
    vm_eci::Array{Float64,1}
end

mutable struct Quaternions
    q0::Float64
    q1::Float64
    q2::Float64
    q3::Float64
end

mutable struct Euler_Angles
    roll::Float64
    pitch::Float64
    heading::Float64
    Tr_bn::Array{Float64,2}
end

mutable struct Angle_Rates
    wb_bi_body::Array{Float64,1}
    wb_bn_body::Array{Float64,1}
    wb_ei_eci::Array{Float64,1}
    wb_ei_ecef::Array{Float64,1}
    wb_ne_eci::Array{Float64,1}
end

mutable struct MassProps
    mass::Float64
    II::Array{Float64,2}
    cg::Array{Float64,1}
end

mutable struct Aero6DOF
    time::Float64
    dt::Float64
    position::WGS_Pos
    velocity::Velocity_States
    euler::Euler_Angles
    angle_rates::Angle_Rates
    accel_meas::Array{Float64,1}
    massprops::MassProps
    quat_bi::Quaternions
    xo::Array{Float64,1}
end

mutable struct Aero5DOF
    time::Float64
    dt::Float64
    position::WGS_Pos
    velocity::Velocity_States
    euler::Euler_Angles
    angle_rates::Angle_Rates
    accel_meas::Array{Float64,1}
    mass::Float64
    quat_bi::Quaternions
    xo::Array{Float64,1}
end

mutable struct NavStates
    time::Float64
    wgs_pos::WGS_Pos
    velocity::Velocity_States
    rate_estimate::Angle_Rates
    rate_meas::Array{Float64,1}
    accel_meas::Array{Float64,1}
    euler::Euler_Angles
    quat::Quaternions
end

mutable struct Actuator
    time::Float64
    dt::Float64
    wn::Float64
    zeta::Float64
    position::Float64
    rate::Float64
    pos_lim_pos::Float64
    pos_lim_neg::Float64
    rate_lim_pos::Float64
    rate_lim_neg::Float64
end

include("Utils.jl");
include("Actuator.jl");
include("Airframe6DOF.jl");
include("Airframe5DOF.jl");
include("Interpolation.jl");
using LinearAlgebra

end # module


text:
# Tabix
# =====
#
# Generic index for tab-delimited files.
#
# Li, Heng. "Tabix: fast retrieval of sequence features from generic
# TAB-delimited files." Bioinformatics 27.5 (2011): 718-719.
# Specification: http://samtools.github.io/hts-specs/tabix.pdf

include("index.jl")
include("reader.jl")

"""
    overlapchunks(index::Tabix, seqname::String, interval::UnitRange)
    overlapchunks(index::Tabix, seqid::Integer, interval::UnitRange)

Return chunks possibly overlapping the range specified by `seqname` (or `seqid`)
and `interval`; `seqid` and `interval` must be 1-based index and inclusive.

NOTE: Records within the returned chunks are not guaranteed to overlap the query
region. The user need to check the condition as:
    line = readline(stream)
    values = split(chomp(line), '\t')
    seqname = value[index.columns[1]]
    start = parse(Int, index.columns[2])
    stop = parse(Int, index.columns[3]) - 1  # needs -1 for the BED format
    if seqname == target_seqname && intersect(start:stop, target_interval)
        # the line overlaps the query region
    end
"""
function overlapchunks(index::Union{BAI,Tabix}, seqid::Integer, interval::UnitRange)
    if !(1 ≤ seqid ≤ endof(index.indexes))
        throw(ArgumentError("sequence id $(seqid) is out of range"))
    end

    if isempty(interval)
        return Chunk[]
    end

    binindex, linindex, pbin = index.indexes[seqid]
    bins = reg2bins(first(interval), last(interval))
    offset = linindex[cld(first(interval), LinearWindowSize)]
    ret = Chunk[]
    for bin in bins
        if haskey(binindex, bin)
            for chunk in binindex[bin]
                if chunk.stop > offset
                    push!(ret, chunk)
                end
            end
        end
    end

    sort!(ret)
    reduce!(ret)

    return ret
end

function overlapchunks(index::Tabix, seqname::AbstractString, interval::UnitRange)
    seqid = findfirst(index.names, seqname)
    if seqid == 0
        throw(ArgumentError("sequence name $(seqname) is not included in the index"))
    end
    return overlapchunks(index, seqid, interval)
end

# Merge chunks to minimize the number of seek operations.
function reduce!(chunks)
    @assert issorted(chunks)
    i = 1
    while i < endof(chunks)
        chunk = chunks[i]
        next = chunks[i+1]
        if chunk.stop < next.start
            # neither overlapping nor adjacent
            i += 1
            continue
        end
        chunks[i] = Chunk(chunk.start, max(chunk.stop, next.stop))
        deleteat!(chunks, i + 1)
    end
    return chunks
end


text:
using Kalman
using GaussianDistributions
using DynamicIterators

Mo(θ) = let Φ = θ,
    b = 0.0,
    Q = 1.0,

    H = 1.0,
    R = 1.0,

    E = LinearEvolution(Φ, Gaussian(b, Q))

    LinearObservation(E, H, R)
end

θ0 = 1.0
M0 = Mo(θ0)
x = 0.2
Y = [0 => 0.2]

Y = collecttrace(M0.P, Trace(Sample(0 => x), Y, endtime(1000)));

x0 = 1.0
P0 = 1.0

target(x) = DynamicIterators._lastiterate(bind(Y, Mo(x)), Start(Kalman.Filter(0 => Gaussian(x0, P0), 0.0)))
@show target(1.0)
@time target(1.0)
true


text:
abstract type AbstractResource end

struct Resources{T<:AbstractResource}
    render_ptr::Ptr{SDL.Renderer}
    d::Dict{String, T}
end
Resources{T}(window::Window) where {T} = Resources(window.render_ptr, Dict{String, T}())
Resources(window::Window) = Resources{AbstractResource}(window)

Base.getindex(r::Resources, key) = getindex(r.d, key)
Base.setindex!(r::Resources, key, val) = setindex!(r.d, key, val)
Base.keys(r::Resources) = keys(r.d)
Base.values(r::Resources) = values(r.d)
Base.delete!(r::Resources{T}, key::T) where {T<:AbstractResource} = (delete!(r.d, key); r)


text:
# Stolen from examples/polymake.jl
# Needs polymake (obviously)

# Returns all vectors v such that Av == b and Cv >= 0.

export solve_non_negative, solve_mixed, solve_ineq

nrows(A::Polymake.MatrixAllocated) = Int(size(A)[1])
ncols(A::Polymake.MatrixAllocated) = Int(size(A)[2])

function _polytope(; A::fmpz_mat=zero_matrix(FlintZZ, 1, 1), b::fmpz_mat=zero_matrix(FlintZZ, ncols(A), 1), C::fmpz_mat=zero_matrix(FlintZZ, 1, 1))
  if !iszero(A)
    bA = Matrix{BigInt}(hcat(-b, A))
    z = findall(i->!is_zero_row(bA, i), 1:nrows(bA))
    zbA = Matrix{BigInt}(bA[z, :])
  else
    zbA = Matrix{BigInt}(undef, 0, 0)
  end
  if !iszero(C)
    z = findall(i->!is_zero_row(C, i), 1:nrows(C))
    zI = Matrix{BigInt}(hcat(zero_matrix(FlintZZ, nrows(C), 1), C))[z, :]
  else
    zI = Matrix{BigInt}(undef, 0, 0)
  end
  if length(zbA) == 0
    p =  Polymake.polytope.Polytope(INEQUALITIES = zI)
  else
    if nrows(zI) == 0
      p = Polymake.polytope.Polytope(EQUATIONS = zbA)
    else
      p = Polymake.polytope.Polytope(EQUATIONS = zbA, INEQUALITIES = zI)
    end
  end
  return p
end

@doc Markdown.doc"""
    solve_ineq(A::fmpz_mat, b::fmpz_mat)

Solves $Ax<=b$, assumes finite set of solutions.
"""
function solve_ineq(A::fmpz_mat, b::fmpz_mat)
  p = _polytope(C = hcat(b, -A))
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

@doc Markdown.doc"""
    solve_non_negative(A::fmpz_mat, b::fmpz_mat)

Finds all solutions to $Ax = b$, $x>=0$. Assumes a finite set of solutions.
"""
function solve_non_negative(A::fmpz_mat, b::fmpz_mat)
  p = _polytope(A = A, b = b, C = identity_matrix(FlintZZ, ncols(A)))
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, nrows(inner) + nrows(out), ncols(A))
  for i=1:nrows(out)
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

@doc Markdown.doc"""
    solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)

Solves $Ax = b$ under $Cx >= 0$, assumes a finite solution set.
"""
function solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat)  # Ax == b && Cx >= 0
  p = _polytope(A = A, b = b, C = C)
  inner = p.INTERIOR_LATTICE_POINTS
  out = p.BOUNDARY_LATTICE_POINTS

  res = zero_matrix(FlintZZ, (nrows(inner) + nrows(out))::Int, ncols(A)::Int)
  for i=1:nrows(out)
    if out[i,1] != 1
      println("unbounded polytope!!")
      global last_in = (A, b, C)
    end
    @assert out[i,1] == 1
    for j=1:ncols(A)
      res[i,j] = out[i, j+1]
    end
  end
  for i=1:nrows(inner)
    if inner[i,1] != 1
      println("unbounded polytope!!")
      global last_in = (A, b, C)
    end
    @assert inner[i,1] == 1
    for j=1:ncols(A)
      res[i+nrows(out), j] = inner[i, j+1]
    end
  end
  return res
end

@doc Markdown.doc"""
    solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat, d::fmpz_mat)

Solves $Ax = b$ under $Cx >= d$, assumes a finite solution set.
"""
function solve_mixed(A::fmpz_mat, b::fmpz_mat, C::fmpz_mat, d::fmpz_mat)
  n = ncols(A)
  A = cat(A, identity_matrix(FlintZZ, ncols(d)), dims=(1,2))
  b = vcat(b, identity_matrix(FlintZZ, ncols(d)))
  C = [C -d; zero_matrix(FlintZZ, ncols(d), ncols(C)) identity_matrix(FlintZZ, ncols(d))]
  s = solve_mixed(A, b, C)
  return s[:, 1:n]
end


text:
function normal_model(t :: Trace, data :: Vector{Float64})
    loc = sample(t, :loc, Normal(0.0, 10.0))
    scale = sample(t, :scale, LogNormal())
    plate(t, observe, :obs, Normal(loc, scale), data)
end

function cat_model(t, dim, data)
    probs = sample(t, "probs", Dirichlet(ones(dim)))
    plate(t, observe, "data", Categorical(probs), data)
end

@testset "aic 1" begin
    data = [1, 2, 1, 3, 4, 3, 3, 3, 4, 2]
    dim = 4
    results = mh(cat_model; params = (dim, data), burn=1000, thin=100, num_iterations=21000)
    @time cat_aic = aic(results)
    @info "Computed AIC of cat model: $cat_aic"
    @info "Cat model has $(numparams(results.traces[1])) parameters"
end

@testset "hpds 1" begin
    data = randn(10) .* 2.0 .+ 4.0
    results = mh(normal_model; params = (data,), burn=1000, thin=100, num_iterations=21000)
    @time intervals = hpdi(results, 0.95, [:loc, :scale])
    @info "HPDI: $intervals"
    normal_aic = aic(results)
    @info "Computed AIC of normal model: $normal_aic"
    @info "Normal model has $(numparams(results.traces[1])) parameters"
end


text:
# This file was generated, do not modify it. # hide
levels(y) |> pprint


text:
using ResettableStacks
using Test
using Random

S = ResettableStack{}(Tuple{Float64,Float64,Float64})

push!(S,(0.5,0.4,0.3))
push!(S,(0.5,0.4,0.4))
reset!(S)
push!(S,(0.5,0.4,0.3))
@test S.data[1] == (0.5,0.4,0.3)

S = ResettableStack{}(Float64)
for i=1:10
  push!(S,i)
end
@test pop!(S) == 10

### Iterator tests
s = ResettableStacks.ResettableStack{}(Float64)
Random.seed!(100)
for i in 1:6
  push!(s,rand())
end
for (i,c) in enumerate(s)
  @test i==1 ? c == 0.6456910432314067 : true
  @test i==2 ? c == 0.9675998379215747 : true
  @test i==3 ? c == 0.06719317094984745 : true
  @test i==4 ? c == 0.6609109399808133 : true
  @test i==5 ? c == 0.19031281518127185 : true
  @test i==6 ? c == 0.2601250914736861 : true
end
reset!(s)
push!(s,rand())
for c in s
  @test c==0.5459681993995775
end


text:
# file to solve the kkt system
mutable struct NewtonSystem{T}
    J::AbstractMatrix{T}
    F::AbstractVector{T}
    d_result::AbstractVector{T}
    dx_aux::AbstractVector{T}
    dp_aux::AbstractVector{T}

    function NewtonSystem{T}(lpip_pb::LPIPLinearProblem{T}, params::Params) where T
        J = zeros(T, 2lpip_pb.n + lpip_pb.m, 2lpip_pb.n + lpip_pb.m)
        build_first_J!(J, lpip_pb.A, lpip_pb.variables.x, lpip_pb.variables.s, lpip_pb.m, lpip_pb.n)
        return new(J,
                   zeros(T, 2lpip_pb.n + lpip_pb.m) ,
                   zeros(T, 2lpip_pb.n + lpip_pb.m) ,
                   zeros(T, lpip_pb.m),
                   zeros(T, lpip_pb.n))
    end
end

function build_J_F!(newton_system::NewtonSystem{T}, lpip_pb::LPIPLinearProblem{T}, params::Params) where T
    mi = params.rho * dot(lpip_pb.variables.x, lpip_pb.variables.s) / lpip_pb.n

    # Build J
    build_J!(newton_system.J, lpip_pb.A, lpip_pb.variables.x, lpip_pb.variables.s, lpip_pb.m, lpip_pb.n)

    # Build F
    build_F!(newton_system.F, newton_system.dx_aux, newton_system.dp_aux,
             lpip_pb.A, lpip_pb.b, lpip_pb.c,
             lpip_pb.variables.x, lpip_pb.variables.s,
             lpip_pb.variables.p, mi, lpip_pb.m, lpip_pb.n)

    return
end

function build_first_J!(J::AbstractMatrix{T}, A::AbstractMatrix{T}, x::AbstractVector{T},
                        s::AbstractVector{T}, m::Int, n::Int) where T

    @inbounds for j in 1:n
        J[m + n + j, j] = s[j]
        J[m + j, m + n + j] = one(T)
        J[m + n + j, m + n + j] = x[j]
        for i in 1:m
            J[i, j] = A[i, j]
            J[m + j, n + i] = A[i, j]
        end
    end
    return
end

function build_J!(J::AbstractMatrix{T}, A::AbstractMatrix{T}, x::AbstractVector{T},
                  s::AbstractVector{T}, m::Int, n::Int) where T

    @inbounds for j in 1:n
        J[m + n + j, j] = s[j]
        J[m + n + j, m + n + j] = x[j]
    end
    return
end

function build_F!(F::AbstractVector{T}, dx_aux::AbstractVector{T}, dp_aux::AbstractVector{T},
                  A::AbstractMatrix{T}, b::AbstractVector{T},
                  c::AbstractVector{T}, x::AbstractVector{T}, s::AbstractVector{T},
                  p::AbstractVector{T}, mi::T, m::Int, n::Int) where T

    LinearAlgebra.BLAS.gemv!('N', 1.0, A, x, 0.0, dx_aux) # newton_system.d.dx_aux = Ax
    LinearAlgebra.BLAS.gemv!('T', 1.0, A, p, 0.0, dp_aux) # newton_system.d.dp_aux = A'p
    @. F[(n + m + 1):(n + m + n)] = x*s - mi
    @. F[(m + 1):(m + n)] = dp_aux + s - c
    @. F[1:m] = dx_aux - b
    return
end

function solve_kkt(newton_system::NewtonSystem{T}, lpip_pb::LPIPLinearProblem{T}, params::Params) where T
    build_J_F!(newton_system, lpip_pb, params)
    newton_system.d_result = newton_system.J \ - newton_system.F # solve the system
    return
end


text:
export
    LogReg,
    logistic,
    loglogistic,
    gradloglogistic

"""
    LogReg(X, y, b)

Data Model for a Logistic Regression with `X` the design matrix, `y` the
response in `(-1,1)` and `b` is the proxy for the Lipschitz constant in order
to use thinning in the BPS algorithm.
"""
struct LogReg
    n::Int # number of observations
    p::Int # number of dimensions
    X::Matrix{Real} # design matrix
    y::Vector{Real} # response
    b::Real  # proxy for N*L, L=Lipschitz constant
    # constructor
    LogReg(X,y,b) = new(size(X,1), size(X,2), X, y, b)
end

# robust log(1+exp(x)) (from JuliaStats/StatsFuns.jl/src/basicfuns.jl)
log1pexp(x::Real) = x < 18.0 ? log1p(exp(x)) : x < 33.3 ? x + exp(-x) : x

# logistic and related functions
logistic(x::Real) = 1.0 / (1.0 + exp(-x))
loglogistic(x::Real) = -log1pexp(-x)
gradloglogistic(x::Real) = logistic(-x)

"""
    loglik(lr, w)

Loglikelihood of a Logistic Regression model `lr` at `w`.
"""
loglik(lr::LogReg, w::Vector{<:Real}) = sum(loglogistic.(lr.y .* (lr.X * w)))

# helper function, to be used only here; the big term on the right hand side
# is used in all the gradient computations
glli(w, lr, i) = gradloglogistic(lr.y[i] * (dot(lr.X[i, :], w)))

"""
    gradloglik(lr, w)

Gradient of the loglikelihood of a Logistic Regression model `lr` at `w`
"""
gradloglik(lr::LogReg, w::Vector{<:Real}) =
    sum(lr.y[i] * glli(w,lr,i) * lr.X[i,:] for i ∈ 1:lr.n)

"""
    gradloglik_cv(lr, wstar)

Return a function that can computes a gradient of the loglikelihood for the
model `lr` around a point `wstar`.
"""
function gradloglik_cv(lr::LogReg, wstar::Vector{<:Real})
    # Unbiased estimate of ∇U_j(x) using control variates
    gll_star = gradloglik(lr, wstar)

    function gll_cv(w::Vector{<:Real}, i=rand(1:lr.n))
        tw, tws = glli(w, lr, i), glli(wstar, lr, i)
        # unbiased estimate of the gradient
        gll_star + lr.n * lr.y[i] * (tw * lr.X[i, :] - tws * lr.X[i,:])
    end
    return gll_cv
end


text:
# Kernel Functions =========================================================================

abstract type Kernel{T<:AbstractFloat} end

function string(κ::Kernel)
    args = [string(getvalue(getfield(κ,θ))) for θ in fieldnames(typeof(κ))]
    kernelname = typeof(κ).name.name
    string(kernelname, "(", join(args, ","), ")")
end

function show(io::IO, κ::Kernel)
    print(io, string(κ))
end

function basefunction(::Kernel)
    error("No base function specified for kernel")
end

@inline eltype(::Type{<:Kernel{E}}) where {E} = E
@inline eltype(κ::Kernel) = eltype(typeof(κ))

"""
    ismercer(κ::Kernel)

Returns `true` if kernel `κ` is a Mercer kernel; `false` otherwise.
"""
ismercer(::Kernel) = false

"""
    isnegdef(κ::Kernel)

Returns `true` if the kernel `κ` is a negative definite kernel; `false` otherwise.
"""
isnegdef(::Kernel) = false

"""
    isstationary(κ::Kernel)

Returns `true` if the kernel `κ` is a stationary kernel; `false` otherwise.
"""
isstationary(κ::Kernel) = isstationary(basefunction(κ))

"""
    isisotropic(κ::Kernel)

Returns `true` if the kernel `κ` is an isotropic kernel; `false` otherwise.
"""
isisotropic(κ::Kernel)  = isisotropic(basefunction(κ))

thetafieldnames(κ::Kernel) = fieldnames(typeof(κ))

gettheta(κ::Kernel{T}) where {T} = T[gettheta(getfield(κ,θ)) for θ in thetafieldnames(κ)]

function settheta!(κ::Kernel{T},v::Vector{T}) where {T}
    fields = thetafieldnames(κ)
    if length(fields) != length(v)
        throw(DimensionMismatch("Update vector has invalid length"))
    end
    for i in eachindex(fields)
        settheta!(getfield(κ, fields[i]), v[i])
    end
    return κ
end

function checktheta(κ::Kernel{T},v::Vector{T}) where {T}
    fields = thetafieldnames(κ)
    if length(fields) != length(v)
        throw(DimensionMismatch("Update vector has invalid length"))
    end
    for i in eachindex(fields)
        if !checktheta(getfield(κ, fields[i]), v[i])
            return false
        end
    end
    return true
end

function floattype(T_i::DataType...)
    T_max = promote_type(T_i...)
    T_max <: AbstractFloat ? T_max : Float64
end

# Mercer Kernels ===========================================================================

abstract type MercerKernel{T<:AbstractFloat} <: Kernel{T} end
@inline ismercer(::MercerKernel) = true

const mercer_kernels = [
    "exponential",
    "squaredexponential",
    "gammaexponential",
    "rationalquadratic",
    "gammarational",
    "matern",
    "linear",
    "polynomial",
    "exponentiated",
    "periodic"
    ]

for kname in mercer_kernels
    include(joinpath("kernelfunctions", "mercer", "$(kname).jl"))
end

# Negative Definite Kernels ================================================================

abstract type NegativeDefiniteKernel{T<:AbstractFloat} <: Kernel{T} end
@inline isnegdef(::NegativeDefiniteKernel) = true

const negdef_kernels = [
    "power",
    "log"
]

for kname in negdef_kernels
    include(joinpath("kernelfunctions", "negativedefinite", "$(kname).jl"))
end

# Other Kernels ============================================================================

const other_kernels = [
    "sigmoid"
]

for kname in other_kernels
    include(joinpath("kernelfunctions", "$(kname).jl"))
end

for κ in [
        ExponentialKernel,
        SquaredExponentialKernel,
        GammaExponentialKernel,
        RationalQuadraticKernel,
        GammaRationalKernel,
        MaternKernel,
        LinearKernel,
        PolynomialKernel,
        ExponentiatedKernel,
        PeriodicKernel,
        PowerKernel,
        LogKernel,
        SigmoidKernel
    ]
    κ_sym = nameof(κ)
    κ_args = [:(getvalue(κ.$(θ))) for θ in fieldnames(κ)]

    @eval begin
        function ==(κ1::$(κ_sym), κ2::$(κ_sym))
            mapreduce(θ -> getfield(κ1,θ) == getfield(κ2,θ), &, fieldnames(typeof(κ1)), init = true)
        end
    end

    @eval begin
        function convert(::Type{$(κ_sym){T}}, κ::$(κ_sym)) where {T}
            $(Expr(:call, :($(κ_sym){T}), κ_args...))
        end
    end

    κs = supertype(κ)
    while κs != Any
        @eval begin
            function convert(::Type{$(nameof(κs)){T}}, κ::$(κ_sym)) where {T}
                convert($(κ_sym){T}, κ)
            end
        end
        κs = supertype(κs)
    end
end


text:
import Pkg

TRAVIS_BUILD_DIR = ENV["TRAVIS_BUILD_DIR"]

pushfirst!(Base.LOAD_PATH,joinpath(TRAVIS_BUILD_DIR, "ci", "RegistryTestingTools",),);

import RegistryTestingTools

RegistryTestingTools.test_registry(TRAVIS_BUILD_DIR)


text:
@testset "Bitcoin address functions" begin
        @testset "address" begin
                secret = 888^3
                mainnet_address = "148dY81A9BmdpMhvYEVznrM45kWN32vSCN"
                testnet_address = "mieaqB68xDCtbUBYFoUNcmZNwk74xcBfTP"
                point = secret * Secp256k1.G
                @test point2address(point, true, false) == mainnet_address
                @test point2address(point, true, true) == testnet_address
                secret = 321
                mainnet_address = "1S6g2xBJSED7Qr9CYZib5f4PYVhHZiVfj"
                testnet_address = "mfx3y63A7TfTtXKkv7Y6QzsPFY6QCBCXiP"
                point = secret * Secp256k1.G
                @test point2address(point, false, false) == mainnet_address
                @test point2address(point, false, true) == testnet_address
                secret = 4242424242
                mainnet_address = "1226JSptcStqn4Yq9aAmNXdwdc2ixuH9nb"
                testnet_address = "mgY3bVusRUL6ZB2Ss999CSrGVbdRwVpM8s"
                point = secret * Secp256k1.G
                @test point2address(point, false, false) == mainnet_address
                @test point2address(point, false, true) == testnet_address
                @testset "P2PKH" begin
                        h160 = hex2bytes("74d691da1574e6b3c192ecfb52cc8984ee7b6c56")
                        want = "1BenRpVUFK65JFWcQSuHnJKzc4M8ZP8Eqa"
                        @test h160_2_address(h160, false, "P2PKH") == want
                        want = "mrAjisaT4LXL5MzE81sfcDYKU3wqWSvf9q"
                        @test h160_2_address(h160, true, "P2PKH") == want
                end
                @testset "P2SH" begin
                        h160 = hex2bytes("74d691da1574e6b3c192ecfb52cc8984ee7b6c56")
                        want = "3CLoMMyuoDQTPRD3XYZtCvgvkadrAdvdXh"
                        @test h160_2_address(h160, false, "P2SH") == want
                        want = "2N3u1R6uwQfuobCqbCgBkpsgBxvr1tZpe7B"
                        @test h160_2_address(h160, true, "P2SH") == want
                end
        end
        @testset "WIF" begin
                kp = KeyPair{:ECDSA}(big(2)^256-big(2)^199)
                expected = "L5oLkpV3aqBJ4BgssVAsax1iRa77G5CVYnv9adQ6Z87te7TyUdSC"
                @test wif(kp, true, false) == expected
                kp = KeyPair{:ECDSA}(big(2)^256-big(2)^201)
                expected = "93XfLeifX7Jx7n7ELGMAf1SUR6f9kgQs8Xke8WStMwUtrDucMzn"
                @test wif(kp, false, true) == expected
                kp = KeyPair{:ECDSA}(big"0x0dba685b4511dbd3d368e5c4358a1277de9486447af7b3604a69b8d9d8b7889d")
                expected = "5HvLFPDVgFZRK9cd4C5jcWki5Skz6fmKqi1GQJf5ZoMofid2Dty"
                @test wif(kp, false, false) == expected
                kp = KeyPair{:ECDSA}(big"0x1cca23de92fd1862fb5b76e5f4f50eb082165e5191e116c18ed1a6b24be6a53f")
                expected = "cNYfWuhDpbNM1JWc3c6JTrtrFVxU4AGhUKgw5f93NP2QaBqmxKkg"
                @test wif(kp, true, true) == expected
        end
end


text:
vis = Visualizer()
open(vis)

# Create Roadway
roadway_opts = HighwayRoadwayOptions()
roadway = build_roadway(roadway_opts)

# Create players
T = Float64
p = 2
model = UnicycleGame(p=p)

players = Vector{Player{T}}(undef, p)
players[1] = Player(model, roadway.lane[1],
	Q=Diagonal(SVector{model.ni[2],T}([0,0.1,1,0])),
    x0=VehicleState(0.0,  0.1, 0.0, 0.05),
    xf=VehicleState(4.0,  0.1, 0.0, 0.05),
	μ=0.0)
players[2] = Player(model, roadway.lane[3],
	Q=Diagonal(SVector{model.ni[2],T}([0,5,0.5,10])),
    x0=VehicleState(0.05,-0.10, 0.0, 0.05),
    xf=VehicleState(4.0,  0.10, 0.0, 0.05),
	μ=3.0)
# players[3] = Player(model, roadway.lane[3],
#     x0=VehicleState(0.2,-0.1, 0.0, 0.3),
#     xf=VehicleState(4.0, 0.1, 0.0, 0.3))

# Create Scenario
sce = Scenario(model, roadway, players)

# Create AutoProblem
N = 40
dt = 0.05
solver_opts = Options()
solver_opts.inner_print = false
solver_opts.outer_print = false
solver_opts.ϵ_dyn = 1e-2
solver_opts.ϵ_sta = 1e-3
solver_opts.ϵ_con = 1e-2
solver_opts.ϵ_opt = 1e-2
prob = GameProblem(N, dt, sce, solver_opts)
@time newton_solve!(prob)

plot!(model, prob.pdtraj.pr)
plot!(prob.stats)

# Visualize trajectory
set_scenario!(vis, sce)
set_camera!(vis)
@time set_state!(vis, model, sce, Algames.state(prob.pdtraj.pr[10]))
# @time set_traj!(vis, model, sce, prob.pdtraj.pr)

function simulate_game!(vis::Visualizer, prob::Algames.GameProblem, x0::SVx, sce::Scenario{T},
    player_id::Int, mpc_stats::MPCStatistics, mpc_opts::MPCOptions) where {T,SVx}
	N = prob.probsize.N
	m = prob.probsize.m
	pu = prob.probsize.pu
    # Initialization
    reset!(mpc_stats)
    prob.x0 = x0
    newton_solve!(prob)
    prob.opts.shift = 0
	no_input_ct = 0
	char = 'N'

    for j = 1:mpc_opts.M
		@show j
        # Solve
		set_env!(vis, model, player_id, prob.x0)
        set_state!(vis, prob.model, sce, prob.x0)

		c, empty_stream  = get_char()
		if empty_stream
			no_input_ct += 1
		else
			char = c
			no_input_ct = 0
		end
		if no_input_ct > 5
			char = 'N'
		end
		xi = standardize(prob.model, prob.x0)[player_id]
		fric = [-0.1*xi.θ, -0.4*xi.v]
		if char == 'N'
			ui = [ 0., 0.] + fric
		elseif char == 'A'
			ui = [ 0.,  0.3] + fric
		elseif char == 'B'
			ui = [ 0., -0.3] + fric
		elseif char == 'D'
			ui = [ 0.3, 0.] + fric
		elseif char == 'C'
			ui = [-0.3, 0.] + fric
		end

		dt_min = 0.01
		dt = @elapsed newton_solve!(prob)
		if dt < dt_min
			sleep(dt_min-dt)
		end
		dt = max(dt_min, dt)

		@show dt
        z0 = prob.pdtraj.pr[1]
        z = Algames.KnotPoint(Algames.state(z0), Algames.control(z0), dt)
        record!(mpc_stats, z)
        if mpc_stats.t > mpc_opts.h
            break
        end
        # Update solver
		# Insert player's control
		traj_rollout = deepcopy(prob.pdtraj.pr)
		for k = 1:N-1
			zk = traj_rollout[k]
			xk = Algames.state(zk)
			uk = Algames.control(zk)
			dtk = zk.dt
			uk = Vector(uk)
			uk[pu[player_id]] = ui
			uk = SVector{m,T}(uk)
			traj_rollout[k] = Algames.KnotPoint(xk,uk,dtk)
		end
		Algames.rollout!(Algames.RK3, prob.model, traj_rollout)

        x,s = first_order_rollout(traj_rollout, dt)
        prob.x0 = x
        prob.opts.shift = s
    end
    return nothing
end

mpc_stats = MPCStatistics()
mpc_opts = MPCOptions(h=50.0, M=100000)
# simulate_MPC!(vis, prob, get_state(sce), sce, mpc_stats, mpc_opts)

simulate_game!(vis, prob, get_state(sce), sce, 1, mpc_stats, mpc_opts)

histogram(mpc_stats.dt)
mpc_stats.dt
mpc_stats.t
@time set_traj!(vis, model, sce, mpc_stats.traj[1:30:end])
#

# MeshCat.convert_frames_to_video(
# 	"home/simon/Downloads/game.tar",
# 	"home/simon/Documents/game.mp4")

# MeshCat.convert_frames_to_video(
# 	"/home/simon/Downloads/box_10_steps.tar",
#  	"/home/simon/Documents/box_10_steps.mp4")


text:
"""
$(TYPEDEF)

Epigraph of real or complex matrix spectral norm (i.e. maximum singular value)
for a matrix (stacked column-wise) of `nrows` rows and `ncols` columns with
`nrows ≤ ncols`.

    $(FUNCTIONNAME){T, R}(nrows::Int, ncols::Int, use_dual::Bool = false)
"""
mutable struct EpiNormSpectral{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d1::Int
    d2::Int
    is_complex::Bool

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_aux_updated::Bool
    inv_hess_aux_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    W_svd
    s::Vector{T}
    U::Matrix{R}
    Vt::Matrix{R}
    mu::Vector{T}
    zeta::Vector{T}
    Urzi::Matrix{R}
    mrziVt::Matrix{R}
    cu::T
    Zu::T
    Uz::Matrix{R}
    sVt::Matrix{R}
    umzdd::Matrix{T}
    simdot::Matrix{T}

    w1::Matrix{R}
    w2::Matrix{R}
    s1::Vector{T}
    s2::Vector{T}
    U1::Matrix{R}
    U2::Matrix{R}
    U3::Matrix{R}
    VVt::Matrix{R}

    function EpiNormSpectral{T, R}(
        d1::Int,
        d2::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert 1 <= d1 <= d2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.is_complex = (R <: Complex)
        cone.dim = 1 + vec_length(R, d1 * d2)
        cone.d1 = d1
        cone.d2 = d2
        return cone
    end
end

reset_data(cone::EpiNormSpectral) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated =
    cone.inv_hess_aux_updated = cone.hess_fact_updated = false)

function setup_extra_data!(
    cone::EpiNormSpectral{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    (d1, d2) = (cone.d1, cone.d2)
    cone.mu = zeros(T, d1)
    cone.zeta = zeros(T, d1)
    cone.Urzi = zeros(R, d1, d1)
    cone.mrziVt = zeros(R, d1, d2)
    cone.Uz = zeros(R, d1, d1)
    cone.sVt = zeros(R, d1, d2)
    cone.umzdd = zeros(T, d1, d1)
    cone.simdot = zeros(T, d1, d1)
    cone.w1 = zeros(R, d1, d2)
    cone.w2 = zeros(R, d1, d2)
    cone.s1 = zeros(T, d1)
    cone.s2 = zeros(T, d1)
    cone.U1 = zeros(R, d1, d1)
    cone.U2 = zeros(R, d1, d1)
    cone.U3 = zeros(R, d1, d1)
    return cone
end

get_nu(cone::EpiNormSpectral) = 1 + cone.d1

function set_initial_point!(
    arr::AbstractVector{T},
    cone::EpiNormSpectral{T},
    ) where {T <: Real}
    arr .= 0
    arr[1] = sqrt(T(get_nu(cone)))
    return arr
end

function update_feas(cone::EpiNormSpectral{T}) where {T <: Real}
    @assert !cone.feas_updated
    u = cone.point[1]
    cone.feas_updated = true
    cone.is_feas = false
    (u > eps(T)) || return false
    @views W = vec_copyto!(cone.w1, cone.point[2:end])

    # fast bounds:
    # spec <= frob <= spec * rtd1
    # opinf / rtd2 <= spec
    # op1 / rtd1 <= spec
    # spec <= sqrt(op1 * opinf)
    frob = norm(W, 2)
    op1 = opnorm(W, 1)
    opinf = opnorm(W, Inf)
    rtd1 = sqrt(T(cone.d1))
    rtd2 = sqrt(T(cone.d2))

    # lower bounds
    lb = max(max(frob, op1) / rtd1, opinf / rtd2)
    (u - lb > eps(T)) || return false

    # upper bounds
    ub = min(frob, sqrt(op1 * opinf))
    if u - ub < eps(T)
        # use fast Cholesky-based feasibility check, rescale W*W' by inv(u)
        Z = mul!(cone.U1, W, W', -inv(u), false)
        @inbounds for i in 1:cone.d1
            Z[i, i] += u
        end
        isposdef(cholesky!(Hermitian(Z, :U), check = false)) || return false
    end

    # compute SVD and final feasibility check
    cone.W_svd = svd(W, full = false) # TODO in place
    cone.is_feas = (u - maximum(cone.W_svd.S) > eps(T))

    return cone.is_feas
end

function is_dual_feas(cone::EpiNormSpectral{T}) where {T <: Real}
    u = cone.dual_point[1]
    (u > eps(T)) || return false
    @views W = vec_copyto!(cone.w1, cone.dual_point[2:end])

    # fast bounds: frob <= nuc <= frob * rtd1
    frob = norm(W, 2)
    (u - sqrt(T(cone.d1)) * frob > eps(T)) && return true
    (u - frob > eps(T)) || return false

    # final feasibility check: nuc = tr(sqrt(W*W')), rescale W*W' by inv(u)
    mul!(cone.U1, W, W', inv(u), false)
    λ = eigvals!(Hermitian(cone.U1, :U))
    return (sqrt(u) - sum(sqrt(abs(λ_i)) for λ_i in λ) > eps(T))
end

function update_grad(cone::EpiNormSpectral{T}) where T
    @assert cone.is_feas
    u = cone.point[1]
    U = cone.U = cone.W_svd.U
    Vt = cone.Vt = cone.W_svd.Vt
    s = cone.s = cone.W_svd.S
    mu = cone.mu
    zeta = cone.zeta
    s1 = cone.s1
    w1 = cone.w1
    U1 = cone.U1
    g = cone.grad

    @. mu = s / u
    @. zeta = T(0.5) * (u - mu * s)
    cone.cu = (cone.d1 - 1) / u

    g[1] = cone.cu - sum(inv, zeta)

    @. s1 = mu / zeta
    mul!(U1, U, Diagonal(s1))
    mul!(w1, U1, Vt)
    @views vec_copyto!(g[2:end], w1)

    cone.grad_updated = true
    return cone.grad
end

function update_hess_aux(cone::EpiNormSpectral)
    @assert !cone.hess_aux_updated
    @assert cone.grad_updated
    s1 = cone.s1
    s2 = cone.s2

    @. s1 = sqrt(cone.zeta)
    @. s2 = cone.mu / s1
    mul!(cone.mrziVt, Diagonal(s2), cone.Vt)
    @. s2 = inv(s1)
    mul!(cone.Urzi, cone.U, Diagonal(s2))

    cone.hess_aux_updated = true
end

function update_hess(cone::EpiNormSpectral{T, R}) where {T, R}
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    d1 = cone.d1
    d2 = cone.d2
    isdefined(cone, :VVt) || (cone.VVt = zeros(R, d2, d2))
    u = cone.point[1]
    zeta = cone.zeta
    Urzi = cone.Urzi
    mrziVt = cone.mrziVt
    s1 = cone.s1
    w1 = cone.w1
    U1 = cone.U1
    H = cone.hess.data

    # u, u
    ui = inv(u)
    H[1, 1] = sum((inv(z_i) - ui) / z_i for z_i in zeta) - cone.cu / u

    # u, w
    @. s1 = -inv(zeta)
    mul!(U1, Urzi, Diagonal(s1))
    mul!(w1, U1, mrziVt)
    @views vec_copyto!(H[1, 2:end], w1)

    # w, w
    Zi = mul!(U1, Urzi, Urzi')
    ZiW = mul!(w1, Urzi, mrziVt)
    WZiWI = mul!(cone.VVt, mrziVt', mrziVt, T(0.5), false)
    @inbounds for i in 1:d2
        WZiWI[i, i] += ui
    end

    idx_incr = (cone.is_complex ? 2 : 1)
    for i in 1:d2, j in 1:d1
        c_idx = r_idx = 2 + idx_incr * ((i - 1) * d1 + j - 1)
        @inbounds for k in i:d2
            ZiWjk = T(0.5) * ZiW[j, k]
            WZiWIik = WZiWI[i, k]
            lstart = (i == k ? j : 1)
            @inbounds for l in lstart:d1
                term1 = Zi[l, j] * WZiWIik
                term2 = ZiW[l, i] * ZiWjk
                spectral_kron_element!(H, r_idx, c_idx, term1, term2)
                c_idx += idx_incr
            end
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiNormSpectral{T},
    ) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    d1 = cone.d1
    u = cone.point[1]
    zeta = cone.zeta
    Urzi = cone.Urzi
    mrziVt = cone.mrziVt
    r = w1 = cone.w1
    simU = w2 = cone.w2
    sim = cone.U1
    S1 = cone.U2
    @views S1diag = S1[diagind(S1)]

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views vec_copyto!(r, arr[2:end, j])

        pui = p / u
        mul!(simU, Urzi', r)
        mul!(sim, simU, mrziVt')
        @. S1 = T(0.5) * (sim + sim')
        @. S1diag -= p / zeta

        prod[1, j] = -sum((pui + real(S1[i, i])) / zeta[i] for i in 1:d1) -
            cone.cu * pui

        mul!(w2, Hermitian(S1, :U), mrziVt, true, inv(u))
        mul!(w1, Urzi, w2)
        @views vec_copyto!(prod[2:end, j], w1)
    end

    return prod
end

function update_inv_hess_aux(cone::EpiNormSpectral{T}) where T
    @assert !cone.inv_hess_aux_updated
    @assert cone.grad_updated
    u = cone.point[1]
    s = cone.s
    zeta = cone.zeta
    umzdd = cone.umzdd
    simdot = cone.simdot

    mul!(cone.sVt, Diagonal(s), cone.Vt)
    mul!(cone.Uz, cone.U, Diagonal(zeta))

    cone.Zu = -cone.cu + sum(inv, u - z_i for z_i in zeta)

    # umzdd = 0.5 * (u .+ mu * s')
    # simdot = zeta ./ (u .- mu * s')
    @inbounds for j in 1:cone.d1
        mu_j = cone.mu[j]
        z_j = zeta[j]
        for i in 1:(j - 1)
            mus_ij = mu_j * s[i]
            umzdd[i, j] = umzdd[j, i] = T(0.5) * (u + mus_ij)
            umus_ij = u - mus_ij
            simdot[i, j] = zeta[i] / umus_ij
            simdot[j, i] = z_j / umus_ij
        end
        umzdd[j, j] = u - z_j
        simdot[j, j] = T(0.5)
    end

    cone.inv_hess_aux_updated = true
end

function update_inv_hess(cone::EpiNormSpectral)
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    d1 = cone.d1
    d2 = cone.d2
    u = cone.point[1]
    U = cone.U
    Uz = cone.Uz
    sVt = cone.sVt
    umzdd = cone.umzdd
    s1 = cone.s1
    w1 = cone.w1
    w2 = cone.w2
    U1 = cone.U1
    U2 = cone.U2
    Hi = cone.inv_hess.data

    # u, u
    hiuu = Hi[1, 1] = u / cone.Zu

    # w, w
    Ut = copyto!(cone.U3, U')
    c_idx = 2
    @inbounds for j in 1:d2, i in 1:d1
        @views U_i = Ut[:, i]
        @views @. U1 = U_i * sVt[:, j]' / umzdd * cone.simdot

        @. U2 = U1 + U1'
        mul!(w2, U2, sVt, -1, false)
        @. @views w2[:, j] += u * U_i
        mul!(w1, Uz, w2)
        @views vec_copyto!(Hi[2:end, c_idx], w1)
        c_idx += 1

        if cone.is_complex
            U1 .*= im
            @. U2 = U1 + U1'
            mul!(w2, U2, sVt, -1, false)
            @. @views w2[:, j] += u * im * U_i
            mul!(w1, Uz, w2)
            @views vec_copyto!(Hi[2:end, c_idx], w1)
            c_idx += 1
        end
    end

    # u, w and w, w
    rthiuu = sqrt(hiuu)
    @inbounds for i in 1:d1
        s1[i] = rthiuu / umzdd[i, i]
    end
    mul!(U1, U, Diagonal(s1))
    mul!(w1, U1, sVt)
    @views Hiuwvec = Hi[1, 2:end]
    vec_copyto!(Hiuwvec, w1)
    @views mul!(Hi[2:end, 2:end], Hiuwvec, Hiuwvec', true, true)
    Hiuwvec .*= rthiuu

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormSpectral,
    )
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    d1 = cone.d1
    u = cone.point[1]
    zeta = cone.zeta
    umzdd = cone.umzdd
    sVt = cone.sVt
    r = w1 = cone.w1
    simU = w2 = cone.w2
    sim = cone.U1
    S1 = cone.U2
    @views S1diag = S1[diagind(S1)]

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views vec_copyto!(r, arr[2:end, j])

        mul!(simU, cone.U', r)
        mul!(sim, simU, sVt')

        c1 = u * (p + sum(real(sim[i, i]) / umzdd[i, i] for i in 1:d1)) / cone.Zu
        prod[1, j] = c1

        sim .*= cone.simdot
        @. S1 = sim + sim'
        @. S1diag -= c1 / zeta
        S1 ./= umzdd

        mul!(w2, Hermitian(S1, :U), sVt, -1, u)
        mul!(w1, cone.Uz, w2)
        @views vec_copyto!(prod[2:end, j], w1)
    end

    return prod
end

function dder3(cone::EpiNormSpectral{T}, dir::AbstractVector{T}) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    u = cone.point[1]
    zeta = cone.zeta
    Urzi = cone.Urzi
    mrziVt = cone.mrziVt
    r = w1 = cone.w1
    simU = w2 = cone.w2
    sim = cone.U1
    S1 = cone.U2
    S2 = cone.U3
    @views S1diag = S1[diagind(S1)]
    @views S2diag = S2[diagind(S2)]
    dder3 = cone.dder3

    p = dir[1]
    @views vec_copyto!(r, dir[2:end])

    pui = p / u
    mul!(simU, Urzi', r)
    mul!(sim, simU, mrziVt')
    @. S1 = T(-0.5) * (sim + sim')
    @. S1diag += p / zeta

    mul!(S2, simU, simU', T(-0.5) / u, false)
    @. S2diag += T(0.5) * p / zeta * pui
    mul!(S2, Hermitian(S1, :U), S1, -1, true)

    @inbounds dder3[1] = -sum((real(S1[i, i]) * pui + real(S2[i, i])) / zeta[i]
        for i in 1:cone.d1) - cone.cu * abs2(pui)

    mul!(w1, Hermitian(S2, :U), mrziVt)
    mul!(w1, Hermitian(S1, :U), simU, inv(u), true)
    mul!(w2, Urzi, w1)
    @views vec_copyto!(dder3[2:end], w2)

    return dder3
end


text:
using DifferentialEquations
using ParameterizedFunctions

ball! = @ode_def BallBounce begin
  dy =  v
  dv = -g
end g

function condition(u,t,integrator)
  u[1]
end

function affect!(integrator)
    integrator.u[2] = -integrator.p[2] * integrator.u[2]
end

bounce_cb = ContinuousCallback(condition,affect!)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
p = (9.8,0.9)
prob = ODEProblem(ball!,u0,tspan,p,callback=bounce_cb)

sol = solve(prob,Tsit5())
using Plots; gr()
plot(sol)

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

function condition_kick(u,t,integrator)
    t == 2
end

function affect_kick!(integrator)
    integrator.u[2] += 50
end

kick_cb = DiscreteCallback(condition_kick,affect_kick!)
u0 = [50.0,0.0]
tspan = (0.0,10.0)
p = (9.8,0.9)
prob = ODEProblem(ball!,u0,tspan,p,callback=kick_cb)

sol = solve(prob,Tsit5(),tstops=[2.0])
plot(sol)


text:
struct Resolution
    hours_in_period::Int

    function Resolution(hours_in_period::Integer)
        return new(hours_in_period)
    end
end

function num_periods(resolution::Resolution, horizon::Horizon)
    hours = num_hours(horizon)
    mod(hours, resolution.hours_in_period) == 0 || error("Time resolution with $(resolution.hours_in_period) not compatible with horizon of $hours hours.")
    return div(hours, resolution.hours_in_period)
end

function water_volume(resolution::Resolution, M::AbstractFloat)
    return M / resolution.hours_in_period
end

function marginal_production(resolution::Resolution, μ::AbstractFloat)
    return resolution.hours_in_period * μ
end

function water_flow_time(resolution::Resolution, R::Integer)
    minutes_in_periods = 60 * resolution.hours_in_period
    return floor(Int, R / minutes_in_periods)
end
function historic_flow(resolution, R::Integer)
    minutes_in_periods = 60 * resolution.hours_in_period
    return mod(R, minutes_in_periods)/minutes_in_periods
end
function overflow(resolution::Resolution, R::Integer)
    minutes_in_periods = 60 * resolution.hours_in_period
    return (1 - mod(R, minutes_in_periods)/minutes_in_periods)
end