text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

using YAActL, Test

struct MyMsg <: Message
    x::Int
end

struct MySource <: Message
    x::Int
    from::Link
end

function writeLk(lk::Link, src=true, slp=false)
    for i in 1:4
        put!(lk.chn, MyMsg(i))
        yield()
        slp && sleep(rand()*0.1)
    end
    if src
        put!(lk.chn, MySource(10, USR))
        yield()
        slp && sleep(rand()*0.1)
    end
    for i in 5:8
        put!(lk.chn, MyMsg(i))
        yield()
        slp && sleep(rand()*0.1)
    end
end

function readLk(lk::Link)
    buf = Int[]
    while isready(lk.chn)
        push!(buf, take!(lk.chn).x)
    end
    buf
end

myLink = Link(10)

@test YAActL._match(MyMsg(1), nothing, nothing)
@test YAActL._match(MySource(1,USR), nothing, nothing)

@test YAActL._match(MyMsg(1), MyMsg, nothing)
@test !YAActL._match(MyMsg(1), Response, nothing)

@test !YAActL._match(MyMsg(1), nothing, USR)
@test YAActL._match(MySource(1,USR), nothing, USR)
@test !YAActL._match(MySource(1,USR), nothing, myLink)

@test !YAActL._match(MyMsg(1), MyMsg, USR)
@test !YAActL._match(MySource(1,USR), MyMsg, USR)
@test YAActL._match(MySource(1,USR), MySource, USR)
@test !YAActL._match(MySource(1,USR), MySource, myLink)

writeLk(myLink)
@test length(myLink.chn.data) == 9
msg = receive!(myLink, MySource, USR)
@test msg == MySource(10, USR)
@test length(myLink.chn.data) == 8
@test readLk(myLink) == collect(1:8)
@test length(myLink.chn.data) == 0

msg = receive!(myLink, MySource, USR, timeout=1)
@test msg == Timeout()

writeLk(myLink, false)
@test length(myLink.chn.data) == 8
msg = receive!(myLink, MySource, USR, timeout=0)
@test msg == Timeout()
@test length(myLink.chn.data) == 8
@test readLk(myLink) == collect(1:8)
@test length(myLink.chn.data) == 0

@async writeLk(myLink, true, true)
msg = receive!(myLink, MySource, USR)
@test msg == MySource(10, USR)
sleep(1)
@test length(myLink.chn.data) == 8
@test readLk(myLink) == collect(1:8)
@test length(myLink.chn.data) == 0

comtest(msg::MySource) = nothing
function comtest(msg::Request)
    res = (1, (2), (3,4,5), [6,7,8])
    if msg.x in 1:4
        send!(msg.from, Response(res[msg.x], self()))
    else
        send!(msg.from, Response("test", self()))
    end
    return nothing
end

A = Actor(comtest)
send!(A, MySource(1, USR))
msg = receive!(USR, timeout=1)
@test msg == Timeout()
send!(A, Request(10, USR))
msg = receive!(USR, timeout=1)
@test msg == Response("test", A)

res = request!(A, MySource(1, USR), timeout=1)
@test res == Timeout()
res = request!(A, Request(1, USR), full=true, timeout=1)
@test res == Response(1, A)
@test request!(A, Request(1, USR), timeout=1) == 1
@test request!(A, Request(2, USR), timeout=1) == 2
@test request!(A, Request(3, USR), timeout=1) == (3,4,5)
@test request!(A, Request(4, USR), timeout=1) == [6,7,8]
@test request!(A, Request(99, USR), timeout=1) == "test"


text:
import Base.findall

#=
function findall(f, a::Array{T, N}) where {T, N}
    j = 1
    b = Vector{Int}(undef, length(a))
    @inbounds for i in eachindex(a)
        @inbounds if f(a[i])
            b[j] = i
            j += 1
        end
    end
    resize!(b, j-1)
    sizehint!(b, length(b))
    return b
end
=#

function findall(t::Union{AbstractString,Regex}, s::AbstractString; overlap::Bool=false)
    found = UnitRange{Int}[]
    i, e = firstindex(s), lastindex(s)
    while true
        r = findnext(t, s, i)
        isnothing(r) && return found
        push!(found, r)
        j = overlap || isempty(r) ? first(r) : last(r)
        j > e && return found
        @inbounds i = nextind(s, j)
    end
end

ext: jl


text:
# Undergraduate Student: Arturo Burgos
# Professor: João Rodrigo Andrade
# Federal University of Uberlândia - UFU, Fluid Mechanics Laboratory - MFLab, Block 5P, Uberlândia, MG, Brazil

using Printf # Just the same as import a module
using Statistics
using CPUTime

println("Choose a N value: ")
N = chomp(readline()) # Here I get the char from the user
n = parse(Int32,N)  # Here I convert the char into an integer(32bits)

A = rand(1:20000,n,n,3);

println("\n\nThe A[1,1,1] is: ",A[1,1,1])
println("\nThe A[1,1,2] is: ",A[1,1,2])
println("\nThe A[1,1,3] is: ",A[1,1,3])
println("\n\nThe A[n,n,1] is: ",A[n,n,1])
println("The A[n,n,2] is: ",A[n,n,2])
println("The A[n,n,3] is: ",A[n,n,3])

function loop(mtx)

    for i = 1:n
        for j=1:n
            mtx[i, j, 1] = mtx[i, j, 2]
            mtx[i, j, 3] = mtx[i, j, 1]
        end
    end
    sleep(1)
    return mtx # retornar algo à variavel alocada, no caso "A"
end

A = loop(A);

println("\n\nThe new A[1,1,1] is: ",A[1,1,1])
println("\nThe new A[1,1,2] is: ",A[1,1,2])
println("\nThe new A[1,1,3] is: ",A[1,1,3])
println("\n\nThe A[n,n,1] is: ",A[n,n,1])
println("The A[n,n,2] is: ",A[n,n,2])
println("The A[n,n,3] is: ",A[n,n,3])

@time @CPUtime loop(A)

max_stars_repo_name: bangyen/Comparing-Languages

max_issues_repo_name: bangyen/Comparing-Languages

max_forks_repo_name: bangyen/Comparing-Languages


text:
module resnet
    using LoopVectorization

    mutable struct ResNet
        layers::Array{Any}
        add_layer::Any
        activate::Any
        initialize::Any
        update::Any
        num_layer::Int64

        loss::Float64

        default_input_shape::Tuple

        function ResNet()
            new(Any[Hidden_Input_Layer(Float32[])], ResNet_add_layer, activate_ResNet, init_ResNet, update_ResNet, 0, 0.0, ())
        end
    end

    mutable struct Hidden_Input_Layer
        output::Array{Float32}
    end
    mutable struct Hidden_Output_Layer
        δ::Array{Float32}
    end

    function ResNet_add_layer(model::ResNet, layer::Any; args...)
        push!(model.layers, layer(;input_shape=model.default_input_shape, args...))
        model.num_layer += 1
        model.default_input_shape = model.layers[end].output_shape
    end

    function init_ResNet(model::ResNet, mini_batch::Int64)
        if length(model.layers)==model.num_layer+1
            push!(model.layers, Hidden_Output_Layer(zeros(Float32, model.layers[end].output_shape..., mini_batch)))
        end
        for i in 2:length(model.layers)-1
            model.layers[i].initialize(model.layers[i], mini_batch)
        end
    end

    function activate_ResNet(model::ResNet, data::Array{Float32})
        model.layers[1].output = data
        for i in 2:length(model.layers)-1
            model.layers[i].activate(model.layers[i], model.layers[i-1].output)
        end
    end

    function update_ResNet(model::ResNet, optimizer::String, α::Float64, parameters...)
        for i in length(model.layers)-1:-1:2
            model.layers[i].update(model.layers[i], optimizer, model.layers[i-1].output, model.layers[i+1].δ, α, parameters)
        end
    end
end


text:
A = [-1 42 23; 44 -5 3; 2 -5 3] #do not change this line
B = [2 3 1; 4 45 -31; 22 12 10] #do not change this line

function get_module_3_excersise_1_author_name()
    return "YOUR NAME HERE"
end

function problem_1_1()
    C = [] # put result in C
    return C
end

function problem_1_2()
    C = [] # put result in C
    return C
end

function problem_1_3()
    C = [] # put result in C
    return C
end


text:
"""
This file contains a set of miscellaneous functions that would be useful not
only in Climate packages but also for other miscellaneous functions such as
trying to lighten colours, or setting the labels of a plot, or making the axis
symmetric in nature, etc.

"""


text:
#! /usr/bin/julia
using TimeSeries
using HTTP
root="/home/jls/data/2020-Corona/raw/RKI/"

function download_rki(source)
    r = HTTP.request("GET", source)
    String(r.body)
end

function parse_rki(pattern,file)
    m=match(pattern,file).match
    all=collect(eachmatch(r">([\d\.]+)<",m))
    confirmed=parse(Float64,replace(all[1][1],"." => ""))
    deaths=parse(Float64,replace(all[3][1],"." => ""))
    return confirmed,deaths
end

function date_rki(file)
    datum=replace(split(match(r"(Stand:.*),",file).match)[2],"," =>"")
    dd=parse(Int,match(r"^(\d+)",datum)[1])
    mm=parse(Int,match(r"\.(\d+)\.",datum)[1])
    yyyy=parse(Int,match(r"(\d+)$",datum)[1])
    d=Date(yyyy,mm,dd)
end

RKI="https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Fallzahlen.html"
rki_file=download_rki(RKI);
d=date_rki(rki_file)

laender=Dict([
    ("Baden-Würtemberg"    ,r"Baden.*Bayern")
    ("Bavaria"             ,r"Bayern.*Berlin")
    ("Berlin"              ,r"Berlin.*Brandenburg")
    ("Rheinland-Pfalz"     ,r"Rhein.*Saar")
    ("Nordrhein-Westphalen",r"Nordrhein-West.*Rhein")
])

for land in laender
    datafile=land[1]*".csv"
#    println(datafile)
    if isfile(root*datafile)
        dataarray=readtimearray(root*datafile,delim=',')
        LastDate=timestamp(dataarray)[end]
    else
        dataarray=TimeArray()
        LastDate=Date(2000)
    end
#    println(LastDate)
    confirmed,deaths=parse_rki(land[2],rki_file)
    #    println(datafile," ",LastDate," $confirmed, $deaths ")
    if d > LastDate
        dataarray=update(dataarray,d ,[confirmed deaths])
        writetimearray(dataarray,root*datafile)
    end
end


text:
using SMLMFrameConnection
using Test

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


text:
module Krotov

include("parametrization.jl")
export SquareParametrization,
    TanhParametrization,
    TanhSqParametrization,
    LogisticParametrization,
    LogisticSqParametrization

include("result.jl")
include("workspace.jl")
include("optimize.jl")

import QuantumControlBase

"""
```julia
opt_result = optimize(problem; method=:krotov, kwargs...)
```

optimizes [`problem`](@ref QuantumControlBase.ControlProblem) using Krotov's
method, see [`Krotov.optimize_krotov`](@ref).
"""
QuantumControlBase.optimize(problem, method::Val{:krotov}) = optimize_krotov(problem)
QuantumControlBase.optimize(problem, method::Val{:Krotov}) = optimize_krotov(problem)

end


text:
function e006(N)
    return difference(N)
end

function difference(N)
    sum_squares = begin
        s = 0
        for i = 1:N
            s += i ^ 2
        end
        s
    end

    square_sum = begin
        s = 0
        for i = 1:N
            s += i
        end
        s ^ 2
    end

    return square_sum - sum_squares
end

function solution006()
    N = 100
    answer = e006(N)
    println("The answer to problem 006 is $answer")
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
# FIELDS:
#==============================================================================#
"""
    abstract type Field end

Abstract type at the top of the type hierarchy for denoting fields over which vector spaces
(or more generally, linear categories) can be defined. Two common fields are `ℝ` and `ℂ`,
representing the field of real or complex numbers respectively.
"""
abstract type Field end

struct RealNumbers <: Field end
struct ComplexNumbers <: Field end

const ℝ = RealNumbers()
const ℂ = ComplexNumbers()

Base.show(io::IO, ::RealNumbers) = print(io, "ℝ")
Base.show(io::IO, ::ComplexNumbers) = print(io, "ℂ")

Base.in(::Any, ::Field) = false
Base.in(::Real, ::RealNumbers) = true
Base.in(::Number, ::ComplexNumbers) = true

Base.issubset(::Type, ::Field) = false
Base.issubset(::Type{<:Real}, ::RealNumbers) = true
Base.issubset(::Type{<:Number}, ::ComplexNumbers) = true
Base.issubset(::RealNumbers, ::RealNumbers) = true
Base.issubset(::RealNumbers, ::ComplexNumbers) = true
Base.issubset(::ComplexNumbers, ::RealNumbers) = false
Base.issubset(::ComplexNumbers, ::ComplexNumbers) = true

# VECTOR SPACES:
#==============================================================================#
"""
    abstract type VectorSpace end

Abstract type at the top of the type hierarchy for denoting vector spaces, or, more
accurately, 𝕜-linear categories. All instances of subtypes of VectorSpace will
represent objects in 𝕜-linear monoidal categories.
"""
abstract type VectorSpace end

"""
    field(V::VectorSpace) -> Field

Return the field type over which a vector space is defined.
"""
function field end
field(V::VectorSpace) = field(typeof(V))

# Basic vector space methods
#----------------------------
"""
    space(a) -> VectorSpace

Return the vector space associated to object `a`.
"""
function space end

"""
    dim(V::VectorSpace) -> Int

Return the total dimension of the vector space `V` as an Int.
"""
function dim end

"""
    dual(V::VectorSpace) -> VectorSpace

Return the dual space of `V`; also obtained via `V'`. This should satisfy
`dual(dual(V)) == V`. It is assumed that `typeof(V) == typeof(V')`.
"""
function dual end

# convenience definitions:
Base.adjoint(V::VectorSpace) = dual(V)

"""
    isdual(V::ElementarySpace) -> Bool

Return wether an ElementarySpace `V` is normal or rather a dual space. Always returns
`false` for spaces where `V == dual(V)`.
"""
function isdual end

# Hierarchy of elementary vector spaces
#---------------------------------------
"""
    abstract type ElementarySpace{𝕜} <: VectorSpace end

Elementary finite-dimensional vector space over a field `𝕜` that can be used as the index
space corresponding to the indices of a tensor. ElementarySpace is a super type for all
vector spaces (objects) that can be associated with the individual indices of a tensor,
as hinted to by its alias IndexSpace.

Every elementary vector space should respond to the methods [`conj`](@ref) and
[`dual`](@ref), returning the complex conjugate space and the dual space respectively. The
complex conjugate of the dual space is obtained as `dual(conj(V)) === conj(dual(V))`. These
different spaces should be of the same type, so that a tensor can be defined as an element
of a homogeneous tensor product of these spaces.
"""
abstract type ElementarySpace{𝕜} <: VectorSpace end
const IndexSpace = ElementarySpace

field(::Type{<:ElementarySpace{𝕜}}) where {𝕜} = 𝕜

"""
    oneunit(V::S) where {S<:ElementarySpace} -> S

Return the corresponding vector space of type `S` that represents the trivial
one-dimensional space, i.e. the space that is isomorphic to the corresponding field. Note
that this is different from `one(V::S)`, which returns the empty product space
`ProductSpace{S,0}(())`.
"""
Base.oneunit(V::ElementarySpace) = oneunit(typeof(V))

"""
    ⊕(V1::S, V2::S, V3::S...) where {S<:ElementarySpace} -> S

Return the corresponding vector space of type `S` that represents the direct sum sum of the
spaces `V1`, `V2`, ... Note that all the individual spaces should have the same value for
[`isdual`](@ref), as otherwise the direct sum is not defined.
"""
function ⊕ end
⊕(V1, V2, V3, V4...) = ⊕(⊕(V1, V2), V3, V4...)

"""
    ⊗(V1::S, V2::S, V3::S...) where {S<:ElementarySpace} -> S

Create a [`ProductSpace{S}(V1, V2, V3...)`](@ref) representing the tensor product of several
elementary vector spaces. For convience, Julia's regular multiplication operator `*` applied
to vector spaces has the same effect.

The tensor product structure is preserved, see [`fuse`](@ref) for returning a single
elementary space of type `S` that is isomorphic to this tensor product.
"""
function ⊗ end
⊗(V1, V2, V3, V4...) = ⊗(⊗(V1, V2), V3, V4...)

# convenience definitions:
Base.:*(V1::VectorSpace, V2::VectorSpace) = ⊗(V1, V2)

"""
    fuse(V1::S, V2::S, V3::S...) where {S<:ElementarySpace} -> S
    fuse(P::ProductSpace{S}) where {S<:ElementarySpace} -> S

Return a single vector space of type `S` that is isomorphic to the fusion product of the
individual spaces `V1`, `V2`, ..., or the spaces contained in `P`.
"""
function fuse end
fuse(V::ElementarySpace) = V
fuse(V1::VectorSpace, V2::VectorSpace, V3::VectorSpace...) =
    fuse(fuse(fuse(V1), fuse(V2)), V3...)
    # calling fuse on V1 and V2 will allow these to be `ProductSpace`

"""
    flip(V::S) where {S<:ElementarySpace} -> S

Return a single vector space of type `S` that has the same value of [`isdual`](@ref) as
`dual(V)`, but yet is isomorphic to `V` rather than to `dual(V)`. The spaces `flip(V)` and
`dual(V)` only differ in the case of [`GradedSpace{I}`](@ref).
"""
function flip end

"""
    conj(V::S) where {S<:ElementarySpace} -> S

Return the conjugate space of `V`. This should satisfy `conj(conj(V)) == V`.

For `field(V)==ℝ`, `conj(V) == V`. It is assumed that `typeof(V) == typeof(conj(V))`.
"""
Base.conj(V::ElementarySpace{ℝ}) = V

"""
    abstract type InnerProductSpace{𝕜} <: ElementarySpace{𝕜} end

Abstract type for denoting vector spaces with an inner product, thus a canonical mapping from
`dual(V)` to `V` (for `𝕜 ⊆ ℝ`) or from `dual(V)` to `conj(V)` (otherwise). This mapping
is provided by the metric, but no further support for working with metrics is currently
implemented.
"""
abstract type InnerProductSpace{𝕜} <: ElementarySpace{𝕜} end

"""
    abstract type EuclideanSpace{𝕜} <: InnerProductSpace{𝕜} end

Abstract type for denoting real or complex spaces with a standard Euclidean inner product
(i.e. orthonormal basis, and the metric is identity), such that the dual space is naturally
isomorphic to the conjugate space `dual(V) == conj(V)` (in the complex case) or even to
the space itself `dual(V) == V` (in the real case), also known as the category of
finite-dimensional Hilbert spaces ``FdHilb``. In the language of categories, this subtype
represents dagger or unitary categories, and support an adjoint operation.
"""
abstract type EuclideanSpace{𝕜} <: InnerProductSpace{𝕜} end # 𝕜 should be ℝ or ℂ

dual(V::EuclideanSpace) = conj(V)
isdual(V::EuclideanSpace{ℝ}) = false

"""
    sectortype(a) -> Type{<:Sector}

Return the type of sector over which object `a` (e.g. a representation space or a tensor) is
defined. Also works in type domain.
"""
sectortype(V::VectorSpace) = sectortype(typeof(V))
sectortype(::Type{<:ElementarySpace}) = Trivial

"""
    hassector(V::VectorSpace, a::Sector) -> Bool

Return whether a vector space `V` has a subspace corresponding to sector `a` with non-zero
dimension, i.e. `dim(V, a) > 0`.
"""
hassector(V::ElementarySpace, ::Trivial) = dim(V) != 0
Base.axes(V::ElementarySpace, ::Trivial) = axes(V)

struct TrivialOrEmptyIterator
    isempty::Bool
end
Base.IteratorSize(::TrivialOrEmptyIterator) = Base.HasLength()
Base.IteratorEltype(::TrivialOrEmptyIterator) = Base.HasEltype()
Base.isempty(V::TrivialOrEmptyIterator) = V.isempty
Base.length(V::TrivialOrEmptyIterator) = isempty(V) ? 0 : 1
Base.eltype(::TrivialOrEmptyIterator) = Trivial
function Base.iterate(V::TrivialOrEmptyIterator, state = true)
    return isempty(V) == state ? nothing : (Trivial(), false)
end

"""
    sectors(V::ElementarySpace)

Return an iterator over the different sectors of `V`.
"""
sectors(V::ElementarySpace) = TrivialOrEmptyIterator(dim(V) == 0)
dim(V::ElementarySpace, ::Trivial) =
    sectortype(V) == Trivial ? dim(V) : throw(SectorMismatch())

# Composite vector spaces
#-------------------------
"""
    abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end

Abstract type for composite spaces that are defined in terms of a number of elementary
vector spaces of a homogeneous type `S<:ElementarySpace{𝕜}`.
"""
abstract type CompositeSpace{S<:ElementarySpace} <: VectorSpace end

spacetype(S::Type{<:ElementarySpace}) = S
spacetype(V::ElementarySpace) = typeof(V) # = spacetype(typeof(V))
spacetype(::Type{<:CompositeSpace{S}}) where S = S
spacetype(V::CompositeSpace) = spacetype(typeof(V)) # = spacetype(typeof(V))

field(P::Type{<:CompositeSpace}) = field(spacetype(P))
sectortype(P::Type{<:CompositeSpace}) = sectortype(spacetype(P))

# make ElementarySpace instances behave similar to ProductSpace instances
blocksectors(V::ElementarySpace) = sectors(V)
blockdim(V::ElementarySpace, c::Sector) = dim(V, c)

# Specific realizations of ElementarySpace types
#------------------------------------------------
# spaces without internal structure
include("cartesianspace.jl")
include("complexspace.jl")
include("generalspace.jl")

# space with internal structure corresponding to the irreducible representations of
# a group, or more generally, the simple objects of a fusion category.
include("gradedspace.jl")

# Specific realizations of CompositeSpace types
#-----------------------------------------------
# a tensor product of N elementary spaces of the same type S
include("productspace.jl")
# deligne tensor product
include("deligne.jl")

# Other examples might include:
# symmetric and antisymmetric subspace of a tensor product of identical vector spaces
# ...

# HomSpace: space of morphisms
#------------------------------
include("homspace.jl")

# Partial order for vector spaces
#---------------------------------
"""
    isisomorphic(V1::VectorSpace, V2::VectorSpace)
    V1 ≅ V2

Return if `V1` and `V2` are isomorphic, meaning that there exists isomorphisms from `V1` to
`V2`, i.e. morphisms with left and right inverses.
"""
function isisomorphic(V1::VectorSpace, V2::VectorSpace)
    spacetype(V1) == spacetype(V2) || return false
    for c in union(blocksectors(V1), blocksectors(V2))
        if blockdim(V1, c) != blockdim(V2, c)
            return false
        end
    end
    return true
end

"""
    ismonomorphic(V1::VectorSpace, V2::VectorSpace)
    V1 ≾ V2

Return whether there exist monomorphisms from `V1` to `V2`, i.e. 'injective' morphisms with
left inverses.
"""
function ismonomorphic(V1::VectorSpace, V2::VectorSpace)
    spacetype(V1) == spacetype(V2) || return false
    for c in blocksectors(V1)
        if blockdim(V1, c) > blockdim(V2, c)
            return false
        end
    end
    return true
end

"""
    isepimorphic(V1::VectorSpace, V2::VectorSpace)
    V1 ≿ V2

Return whether there exist epimorphisms from `V1` to `V2`, i.e. 'surjective' morphisms with
right inverses.
"""
function isepimorphic(V1::VectorSpace, V2::VectorSpace)
    spacetype(V1) == spacetype(V2) || return false
    for c in blocksectors(V2)
        if blockdim(V1, c) < blockdim(V2, c)
            return false
        end
    end
    return true
end

# unicode alternatives
const ≅ = isisomorphic
const ≾ = ismonomorphic
const ≿ = isepimorphic

≺(V1::VectorSpace, V2::VectorSpace) = V1 ≾ V2 && !(V1 ≿ V2)
≻(V1::VectorSpace, V2::VectorSpace) = V1 ≿ V2 && !(V1 ≾ V2)

"""
    infimum(V1::ElementarySpace, V2::ElementarySpace, V3::ElementarySpace...)

Return the infimum of a number of elementary spaces, i.e. an instance `V::ElementarySpace`
such that `V ≾ V1`, `V ≾ V2`, ... and no other `W ≻ V` has this property. This requires
that all arguments have the same value of `isdual( )`, and also the return value `V` will
have the same value.
"""
infimum(V1::ElementarySpace, V2::ElementarySpace, V3::ElementarySpace...) =
    infimum(infimum(V1, V2), V3...)

"""
    supremum(V1::ElementarySpace, V2::ElementarySpace, V3::ElementarySpace...)

Return the supremum of a number of elementary spaces, i.e. an instance `V::ElementarySpace`
such that `V ≿ V1`, `V ≿ V2`, ... and no other `W ≺ V` has this property. This requires
that all arguments have the same value of `isdual( )`, and also the return value `V` will
have the same value.
"""
supremum(V1::ElementarySpace, V2::ElementarySpace, V3::ElementarySpace...) =
    supremum(supremum(V1, V2), V3...)

import Base: min, max


text:
using NFFT

N = 1024
M = N*N
m = 4
σ = 2.0
pre=NFFT.LINEAR
T=Float64
storeDeconvolutionIdx=false
fftflags=NFFT.FFTW.ESTIMATE

x = T.(rand(2,M) .- 0.5)

@time p = NFFTPlan(x, (N,N); m, σ, window=:kaiser_bessel,
             precompute=pre, fftflags, storeDeconvolutionIdx)


text:
module FoodCachingPlotting

using DataFrames, DataFramesMeta, Unitful, Random, StatsBase, Statistics
import PGFPlotsX, PlotlyJS, FoodCachingExperiments, FoodCachingModels
import PGFPlotsX: @pgf, Coordinates, Plot, Axis, Options
import FoodCachingExperiments: Experiment, EXPERIMENTS, CLAYTON0103_EXPERIMENTS,
       asdataframe

export plot_experiment, plot_compare, boxplot

function __init__()
    push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepackage{longtable}\\pgfplotsset{width=5cm}")
    PGFPlotsX.CLASS_OPTIONS[1] = "crop = false"
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"""
    \usepackage{xcolor}
    \definecolor{mblue}{HTML}{0992F2}
    \definecolor{morange}{HTML}{ff7f0e}
    \definecolor{mgreen}{HTML}{2ca02c}
    \definecolor{mred}{HTML}{d62728}
    \pagestyle{empty}
    """)
end

function plotdf(::Val{:pgf}, df, x, y; yerr, group, kwargs...)
    df = copy(df)
    gs = groupby(df, group)
    legend_entries = (x -> join(x, " ")).(collect(keys(gs.keymap)))
    xlabs = union(getproperty(df, x))
    try
        xlabs = Int.(xlabs)
    catch
    end
    xcoords = haskey(kwargs, :xtick) ? Dict(x => i for (i, x) in zip(kwargs[:xtick], xlabs)) : Dict(x => i - 1 for (i, x) in enumerate(xlabs))
    df.xcoords = map(x -> xcoords[x], getproperty(df, x))
    ymax = maximum(getproperty(df, y) .+ (yerr === nothing ? 0 : getproperty(df, yerr)))
    ymin = minimum(getproperty(df, y) .- (yerr === nothing ? 0 : getproperty(df, yerr)))
    ymax += .05 * (ymax - ymin)
    ymin -= .05 * (ymax - ymin)
    options = merge(@pgf({xtick = 0:length(xlabs)-1,
              xticklabels = xlabs, ymax = ymax, ymin = ymin,
               "no marks", xlabel = string(x),
               "error bars/error bar style = very thick", very_thick,
               axis_line_style = {draw = "none"},
               "error bars/error mark = |",
               legend_to_name = randstring(),
               legend_entries = legend_entries, legend_columns = 2,
               legend_style = {draw = "none", very_thick},
               "error bars/error mark options = {very thick, mark size = 3}",
               grid = "major", tick_style = "gray!50",
               "error bars/y explicit",
               "error bars/y dir=both"}), Options(kwargs...))
    Axis(options,
         [Plot(Coordinates(getproperty(d, :xcoords),
                           getproperty(d, y),
                           yerror = yerr === nothing ? nothing : getproperty(d, yerr)))
          for d in gs]...)
end
ticks(x::AbstractVector{Int}) = sort(union(x))
ticks(x) = typeof(x[1]) <: Int ? sort(union(x)) : 0:length(union(x))-1
function plotdf(::Val{:plotly}, df, x, y; yerr, group,
                  color = [:red, :blue, :green, :orange, :cyan, :brown],
                  kwargs...)
    layout = merge(Dict(:xaxis => Dict(:title => string(x),
                                       :zeroline => false,
                                       :tickvals => ticks(getproperty(df, x))),
                        :yaxis => Dict(:zeroline => false),
                        :legend => Dict(:xanchor => "center", :yanchor => "top",
                                        :x => .5, :y => -.3)), kwargs)
    p = PlotlyJS.Plot(df, PlotlyJS.Layout(; layout...); x, y, group, yerr)
    for (i, t) in enumerate(p.data)
        t[:error_y] = Dict(:array => t[:yerr])
        t[:line] = Dict(:color => color[i])
    end
    PlotlyJS.plot(p)
end

function plotdf(df, x, y = :μ;
                yerr = hasproperty(df, :sem) ? :sem : nothing,
                group = setdiff(Symbol.(names(df)), [x, y, yerr, :n]),
                backend = :plotly,
                kwargs...)
    plotdf(Val(backend), df, x, y; yerr, group, renamekw(Val(backend), kwargs)...)
end
function plotdfgroupby(df, group, x, y = :μ; backend = :plotly, kwargs...)
    [plotdf(d, x, y; backend, kwargs...) for d in groupby(df, group)]
end

plot_experiment(e::Symbol; kwargs...) = plot_experiment(EXPERIMENTS[e]; kwargs...)
plot_experiment(e::Symbol, data; kwargs...) = plot_experiment(EXPERIMENTS[e], data; kwargs...)
function plot_experiment(e::Symbol, data::AbstractVector; kwargs...)
    if e == :Clayton0103
        dfs = asdataframe(e, data)
        [plot_experiment(ei, di[1]; kwargs...)
         for (ei, di) in zip(CLAYTON0103_EXPERIMENTS, dfs)]
    else
        plot_experiment(e, asdataframe(e, data)[1]; kwargs...)
    end
end
plot_experiment(e::Experiment; kwargs...) = plot_experiment(e, e.data; kwargs...)

function renamekw(::Val{:plotly}, kwargs)
    rename = Dict(:ylabel => :yaxis_title,
                  :xtick => :xaxis_tickvals,
                  :xticklabels => :xaxis_ticktext)
    [(x -> haskey(rename, x) ? rename[x] : x)(k) => v for (k, v) in kwargs]
end
renamekw(::Any, kwargs) = kwargs
pushorappend!(x, v::AbstractVector) = append!(x, v)
pushorappend!(x, v) = push!(x, v)
function share_yaxis(::Val{:plotly}, plots)
    ps = hcat(plots...)
    y = ps.plot.data[1][:yaxis]
    for p in ps.plot.data
        PlotlyJS.restyle!(p, yaxis = y,
                          legendgroup = p[:name],
                          showlegend = p[:xaxis] == "x1")
    end
    PlotlyJS.plot(ps.plot)
end
function share_yaxis(::Val{:pgf}, plots)
    ymax = maximum([p.options.dict["ymax"] for p in plots])
    ymin = minimum([p.options.dict["ymin"] for p in plots])
    for p in plots
        p.options.dict["ymax"] = ymax
        p.options.dict["ymin"] = ymin
    end
    plots
end
arrange(::Any, ps) = ps
function arrange(::Val{:pgf}, ps)
    m = length(ps[1])
    d = PGFPlotsX.TikzDocument(raw"\begin{longtable}{" * "c"^m * "}")
    for row in ps
        for fig in row
            push!(d, PGFPlotsX.TikzPicture(fig))
            push!(d, " & ")
        end
        pop!(d.elements); push!(d, "\\\\")
        push!(d, "\\multicolumn{$m}{c}{\\ref{$(row[1].options.dict["legend_to_name"])}}\\\\")
    end
    push!(d, raw"\end{longtable}")
end
function plot_compare(e::Symbol, data1, data2 = nothing;
                      backend = :plotly, title1 = "Data",
                      title2 = data2 === nothing ? "Model" : "Model mean",
                      title3 = "Model closest",
                      color = [:red, :blue, :green, :orange, :cyan, :brown])
    if e == :Clayton0103
        dfs1 = asdataframe(e, data1)
        dfs2 = data2 === nothing ? fill(nothing, length(dfs1)) : asdataframe(e, data2)
        return [plot_compare(ei, di1[1], di2[1]; backend, title1, title2, title3, color)
                for (ei, di1, di2) in zip(CLAYTON0103_EXPERIMENTS, dfs1, dfs2)]
    end
    plots = []
    if data2 !== nothing
        pushorappend!(plots, plot_experiment(e, data2; title = title3, backend))
    end
    pushorappend!(plots, plot_experiment(e; title = title1, backend))
    pushorappend!(plots, plot_experiment(e, data1; title = title2, backend))
    m = 2 + (data2 !== nothing)
    n = length(plots) ÷ m
    ps = [share_yaxis(Val(backend), plots[collect(i:n:end)]) for i in 1:n]
    arrange(Val(backend), ps)
end

plot_experiment(::Experiment{:deKort07_exp1}, data; kwargs...) =
    plotdf(data, :trial; kwargs...)

plot_experiment(::Experiment{:deKort07_exp2}, data; kwargs...) =
    plotdf(data, :trial; kwargs...)

plot_experiment(::Experiment{:deKort07_exp3}, data; kwargs...) =
    plotdf(data, :trial; kwargs...)

plot_experiment(::Experiment{:deKort07_exp4}, data; kwargs...) =
    plotdf(data, :tray; kwargs...)

plot_experiment(::Experiment{:Raby07_planningforbreakfast}, data; kwargs...) =
    plotdf(data, :tray; kwargs...)

plot_experiment(::Experiment{:Raby07_breakfastchoice}, data; kwargs...) =
    plotdf(data, :tray; kwargs...)

function plot_experiment(::Experiment{:deKort05}, data; kwargs...)
    [plotdfgroupby(@select(data, :group, :set, :RI, :μ, :sem), :group, :set;
                   kwargs...);
     plotdf(@select(@where(data, :set .== 8), :group, :RI, :set, :firstinspection),
            :RI, :firstinspection;
            ymin = 0, ymax = 1, xtick = [4, 28],
            ytick = [0, 1], ylabel = "f.i. waxworm", kwargs...)]
end

plot_experiment(::Experiment{:Correia07_exp1}, data; kwargs...) =
    plotdf(data, :prefed; ylabel = "average items eaten", kwargs...)

plot_experiment(::Experiment{:Correia07_exp2}, data; kwargs...) =
    plotdfgroupby(data, :action, :trial; xtick = 1:3, kwargs...)

function plot_experiment(e::Experiment{:Cheke11_specsat}, data; kwargs...)
    data = hasproperty(data, :μ) ? data : FoodCachingExperiments.summarize(e, data)
    plotdf(data, :prefed; xtick = 0:1, xticklabels = ["peanut", "suet"], kwargs...)
end

function plot_experiment(e::Experiment{:Cheke11_planning}, data; kwargs...)
    data = hasproperty(data, :μ) ? data : FoodCachingExperiments.summarize(e, data)
    plotdf(data, :trial; xtick = 1:3, kwargs...)
end

plot_experiment(::Experiment{:Clayton99A_exp1}, data; kwargs...) =
    plotdfgroupby(data, :action, :tray; kwargs...)

plot_experiment(::Experiment{:Clayton99A_exp2}, data; kwargs...) =
    plotdf(data, :tray; kwargs...)

function plot_experiment(::Experiment{:Clayton99B_exp1}, data; kwargs...)
    [plotdfgroupby(@select(data, :group, :trial, :RI, :action, :foodtype, :μ, :sem),
                   [:group, :RI], :trial; kwargs...);
     plotdf(@select(@where(data, :trial .== 5,
                           :action .== "inspect", :foodtype .== "peanut"),
                    :group, :RI, :firstinspection),
            :RI, :firstinspection; ylabel = "f.i. peanut", ymin = 0, ymax = 1,
            kwargs...)]
end

function plot_experiment(::Experiment{:Clayton99B_exp2}, data; kwargs...)
    [plotdfgroupby(@select(@where(data, :action .== "inspect"),
                           :group, :trial, :order, :foodtype, :μ, :sem),
                   [:trial, :group], :order; kwargs...);
     plotdfgroupby(@select(@where(data, :action .== "inspect", :foodtype .== "peanut"),
                           :group, :trial, :order, :firstinspection),
                   [:trial], :order, :firstinspection;
                   ylabel = "f.i. peanuts", ymin = 0, ymax = 1, kwargs...)]
end

plot_experiment(::Experiment{:Clayton99C_exp1}, data; kwargs...) =
    plotdfgroupby(data, :foodtype, :group; kwargs...)

plot_experiment(::Experiment{:Clayton99C_exp2}, data; kwargs...) =
    plotdfgroupby(data, :action, :group; kwargs...)

plot_experiment(::Experiment{:Clayton99C_exp3}, data; kwargs...) =
    plotdfgroupby(data, :action, :group; xticklabel_style = @pgf({rotate = 45}),
                  kwargs...)

function plot_experiment(::Experiment{:Clayton01_exp1}, data; kwargs...)
    [plotdfgroupby(@select(@where(data, :set .== 5, :action .== "inspect"),
                           :foodtype, :group, :RI, :trial, :μ, :sem),
                     :RI, :foodtype; kwargs...);
     plotdfgroupby(@select(@where(data, :foodtype .== "peanut",
                                  :action .== "inspect", 0 .< :set .<= 6),
                           :group, :RI, :trial, :set, :firstinspection),
                   :RI, :set, :firstinspection; ylabel = "f.i. peanut",
                   ymin = 0, ymax = 1, kwargs...);
     plotdf(@select(@where(data, :foodtype .== "other", :action .== "inspect",
                           :set .> 4, :trial .== "pc"),
                    :group, :RI, :firstinspection),
            :RI, :firstinspection; ylabel = "f.i. cricket",
            ymin = 0, ymax = 1, kwargs...);
    ]
end

function plot_experiment(::Experiment{:Clayton01_exp2}, data; kwargs...)
    [plotdfgroupby(@select(@where(data, :set .== 6, :action .== "inspect"),
                           :foodtype, :group, :RI, :trial, :μ, :sem),
                   :RI, :foodtype; kwargs...);
     plotdf(@select(@where(data, :foodtype .== "mealworm",
                           :action .== "inspect", :set .== 6),
                    :group, :RI, :trial, :set, :firstinspection),
            :RI, :firstinspection;
            ylabel = "f.i. mealworm", ymin = 0, ymax = 1, kwargs...)]
end

function plot_experiment(::Experiment{:Clayton01_exp3}, data; kwargs...)
    ps = Vector{Any}(undef, 3)
    ps[1] = plotdf(@select(@where(data, :action .== "inspect"),
                           :foodtype, :group, :trial, :μ, :sem),
                     :foodtype; ylabel = "inspect", kwargs...)
    ps[2] = plotdf(@select(@where(data, :action .== "cache"),
                    :foodtype, :group, :μ, :sem),
            :foodtype; ylabel = "cache", kwargs...)
    ps[3] = plotdf(@select(@where(data, :foodtype .== "mealworm",
                                  :action .== "inspect"),
                           :group, :trial, :firstinspection),
                   :trial, :firstinspection; ylabel = "f.i. mealworm",
                   ymin = 0, ymax = 1, kwargs...)
    ps
end

function plot_experiment(::Experiment{:Clayton01_exp4}, data; kwargs...)
    ps = Vector{Any}(undef, 2)
    ps[1] = plotdf(@select(@where(data, :action .== "inspect"),
                           :foodtype, :group, :RI, :μ, :sem),
                     :foodtype; kwargs...);
    ps[2] = plotdf(@select(@where(data, :foodtype .== "peanut",
                                  :action .== "inspect"),
                           :group, :RI, :firstinspection),
                   :RI, :firstinspection; ylabel = "f.i. peanut",
                   ymin = 0, ymax = 1, kwargs...)
    ps
end

function plot_experiment(::Experiment{:Clayton03_exp1}, data; kwargs...)
    ps = Vector{Any}(undef, 2)
    ps[1] = plotdf(@select(@where(data, :foodtype .== "other", :action .== "inspect",
                           :set .> 4, :trial .== "pc"),
                    :group, :RI, :firstinspection),
            :RI, :firstinspection; ylabel = "f.i. cricket",
            ymin = 0, ymax = 1, kwargs...)
    ps[2] = plotdf(@select(@where(dropmissing(data, :foodtype, disallowmissing = true),
                           :foodtype .== "cricketproportion"),
                    :group, :RI, :μ, :sem),
            :RI; ylabel = "f.i. cricket", ymin = 0, ymax = 1, kwargs...)
    ps
end

function plot_experiment(exp::Experiment{:Clayton03_exp2}, data; kwargs...)
    [plotdfgroupby(@select(@where(data, :action .== "inspect"),
                          :group, :foodtype, :condition, :RI, :μ, :sem),
                  [:group, :RI], :foodtype; kwargs...);
     plotdfgroupby(@select(@where(data, :action .== "inspect", :foodtype .== "peanut"),
                           :group, :condition, :RI, :firstinspection),
                   [:group], :RI, :firstinspection; ylabel = "f.i. crickets",
                   ymin = 0, ymax = 1, kwargs...);
     plotdf(@select(@where(data, :action .== "cache"),
                    :group, :foodtype, :μ, :sem),
            :foodtype; kwargs...)]
end

function plot_experiment(e::Experiment{:Clayton0103}; kwargs...)
    plot_experiment(e, [FoodCachingExperiments.EXPERIMENTS[x].data
                        for x in FoodCachingExperiments.CLAYTON0103_EXPERIMENTS];
                    kwargs...)
end
function plot_experiment(::Experiment{:Clayton0103}, data; kwargs...)
    vcat([plot_experiment(e, d; kwargs...)
          for (e, d) in zip(FoodCachingExperiments.CLAYTON0103_EXPERIMENTS, data)]...)
end

plot_experiment(::Experiment{:Clayton05_exp1}, data; kwargs...) =
    plotdfgroupby(data, :group, :trial; kwargs...)

plot_experiment(::Experiment{:Clayton05_exp2}, data; kwargs...) =
    plotdfgroupby(data, :group, :trial; kwargs...)

plot_experiment(exp::Experiment{:Clayton05_exp3}, data; kwargs...) =
    plotdfgroupby(@where(data, :n .< 80), :group, :trial; kwargs...)

plot_experiment(exp::Experiment{:Clayton05_exp4}, data; kwargs...) =
    plotdf(data, :trial; kwargs...)

function trayprocessing(trays, i, f = nothing)
    trays = filter(x -> x.position == i, trays)
    length(trays) > 1 && @warn "multiple trays at position $i"
    isempty(trays) && return 0.
    if f === nothing
        trays[1].closed ? .75 : 1.5
    else
        FoodCachingModels.countfooditems(trays[1], f)
    end
end
function preprocess(field, t, v, foodtypes)
    if field == :ismdpresent
        [(t, v, "maintenance diet")]
    elseif field ∈ (:actions_eat, :actions_cache, :actions_inspect, :actions_other)
        [(t, first.(v), replace(string(field), "actions_" => ""))]
    elseif field ∈ (:eatableitems, :cacheableitems)
        [(t, FoodCachingModels.countfooditems.(v, f), "$field $f")
         for f in foodtypes]
    elseif field ∈ (:hunger, :stomach)
        [(t, getindex.(v, i), "$field $(foodtypes[i])")
         for i in eachindex(foodtypes)]
    elseif field == :trays
        trays = union((x -> x.position).(vcat(v...)))
        [[(t, trayprocessing.(v, i), "tray $i presence")
          for i in trays]...;
         [(t, trayprocessing.(v, i, f), "tray $i $f")
          for i in trays, f in foodtypes]...]
    elseif field == :trayweights
        ks = union(keys.(union(v))...)
        [(t, (x -> haskey(x, k) ? x[k] : missing).(v), "$(k[1]) $(k[2])")
         for k in ks]
    end
end
function traces(model; timeunit = 1.0u"hr")
    data = model.tracker.data
    foodtypes = union((x -> x.id).(vcat(data[:eatableitems].v...)))
    plotsdict = Dict()
    for (field, x) in data
        plotsdict[field] = preprocess(field,
                                      uconvert.(NoUnits, x.t ./ timeunit),
                                      x.v,
                                      foodtypes)
    end
    actions = intersect(keys(plotsdict),
                        [:actions_eat, :actions_cache,
                         :actions_inspect, :actions_other])
    plotsdict[:actions] = vcat([plotsdict[a] for a in actions]...)
    for a in actions
        pop!(plotsdict, a)
    end
    plotsdict
end
function titles(k)
    k == :trays && return "cached items in trays"
    k == :ismdpresent && return "maintentance diet availability"
    k == :cacheableitems && return "cacheable items"
    k == :eatableitems && return "eatable items"
    k == :trayweights && return "tray weights"
    string(k)
end
function plotexperiment(model; order = nothing)
    plotsdict = traces(model)
    if order !== nothing
        plotsdict = sort(plotsdict, by = x -> findfirst(a -> x == a, order))
    end
    data = PlotlyJS.GenericTrace{Dict{Symbol,Any}}[]
    layout = Dict{Symbol, Any}(:xaxis_domain => [0., 1.])
    layout[:annotations] = Dict{Symbol, Any}[]
    n = length(plotsdict)
    h0 = .04
    h = (1 - (n-1)*h0)/n
    i = n + 1
    for (k, v) in plotsdict
        i -= 1
        for (t, d, li) in v
            push!(data, PlotlyJS.scatter(x = t, y = d,
                                         mode = k == :actions ? "markers" : "lines",
                                         name = li, xaxis = "x",
                                         yaxis = "y$i"))
        end
        layout[Symbol(:yaxis, i, :_domain)] = [(i - 1)*h, i*h] .+ (i - 1)*h0
        push!(layout[:annotations],
              Dict(:yanchor => "bottom",
                   :xanchor => "center",
                   :y => (i-.5)*h + (i - 1)*h0,
                   :bgcolor => "#ffffff",
                   :font => Dict(),
                   :showarrow => false,
                   :yref => "paper",
                   :text => titles(k),
                   :xref => "paper",
                   :x => 0.8))
    end
    l = PlotlyJS.Layout(; xaxis = PlotlyJS.attr(title="time [hr]"), layout...)
    PlotlyJS.plot(data, l)
end

function prepare_boxplot(x; options = @pgf({}))
    opt = @pgf {mark_size = ".5pt",
                "every mark/.append style = {fill = white, line width = .2pt}"}
    opt = merge(opt, options)
    uq = quantile(x, .75)
    lq = quantile(x, .25)
    IQR = uq - lq
    lw = minimum(filter(x -> x > lq - 1.5 * IQR, x))
    uw = maximum(filter(x -> x < uq + 1.5 * IQR, x))
    bopt = @pgf {median = median(x), lower_quartile = lq,
                 upper_quartile = uq, lower_whisker = lw,
                 upper_whisker = uw}
    opt, bopt, filter(x -> x > uq + 1.5 * IQR || x < lq - 1.5 * IQR, x)
end
function boxplot(data; labels, axisopt = @pgf({}), options = @pgf({}))
    plots = []
    for (i, group) in enumerate(data)
        N = length(group)
        box_extend = .6/N
        offsets = N == 1 ? [0.] : range(-.35 + box_extend/2,
                                        stop = .35 - box_extend/2,
                                        length = N)
        for (k, x) in enumerate(group)
            opt, bopt, outliers = prepare_boxplot(x; options = options[k])
            bopt = merge(@pgf({box_extend = box_extend, draw_position =
                               offsets[k] + i}), bopt)
            opt = merge(@pgf({boxplot_prepared = bopt, thick}), opt)
            p = @pgf Plot(opt, Coordinates(fill(offsets[k] + i,
                                                          length(outliers)),
                                                         outliers))
            push!(plots, p)
        end
    end
    Axis(merge(@pgf({"boxplot/draw direction = y",
                     xtick = 1:length(labels),
                     ymode = "log",
                     ylabel = raw"$\Delta$",
                     xticklabels = labels,
                     xticklabel_style = {rotate = 70, anchor = "east"}}),
                axisopt), plots...)
end

end


text:
# SPDX-License-Identifier: MIT
# Author: Pierre Lairez

@genpredicate function orient(u :: 2, v :: 2, w :: 2)
    u = u - w
    v = v - w

    group!(u...)
    group!(v...)

    ext(u, v)
end
@doc """

    orient(p :: 2, q :: 2, r :: 2) -> Int

* Return 1 if `r` is on the left of the oriented line defined by `p` and `q`.
* Return –1 if `r` is on the right.
* Return 0 if `r` is on the line or if `p == q`.

""" orient(::NTuple{2, Float64}, ::NTuple{2, Float64}, ::NTuple{2, Float64})

@genpredicate function incircle(p :: 2, q :: 2, r :: 2, a :: 2)

    qp = q - p
    rp = r - p
    ap = a - p
    aq = a - q
    rq = r - q

    group!(qp...)
    group!(ap...)
    group!(rp...)
    group!(rq..., aq...)

    ext(qp, ap)*inp(rp, rq) - ext(qp, rp)*inp(ap, aq)
end
@doc """
    incircle(a :: 2, b :: 2, c :: 2, p :: 2) -> Int

Assume that `a`, `b` and `c` define a counterclockwise triangle.

* Return 1 if `p` is strictly inside the circumcircle of this triangle.
* Return –1 if `p` is outside.
* Return 0 if `p` is on the circle.

If the triangle is oriented clockwise, the signs are reversed.
If `a`, `b` and `c` are collinear, this degenerate to an orientation test.

If two of the four arguments are equal, return 0.

""" incircle(::NTuple{2, Float64}, ::NTuple{2, Float64}, ::NTuple{2, Float64}, ::NTuple{2, Float64})

@genpredicate function incircle2(a :: 2, b :: 2, c :: 2, p :: 2)
    a = a - p
    b = b - p
    c = c - p

    group!(a..., b..., c...)
    abs2(a)*ext(b, c) + abs2(b)*ext(c, a) + abs2(c)*ext(a, b)
end

@genpredicate function closest(p :: 2, q :: 2, a :: 2)
    qp = q - p
    pa = p - a
    qa = q - a

    group!(qp...)
    group!(pa..., qa...)

    inp(qp, pa+qa)
end
@doc """
    closest(p :: 2, q :: 2, a :: 2) -> Int

* Return 1 if `a` is closer to `p` than to `q`.
* Return –1 if `a` is closer to `q` than to `p`.
* Return 0 is `a` is equaly close to both.
""" closest(::NTuple{2, Float64}, ::NTuple{2, Float64}, ::NTuple{2, Float64})

"""
    sameside(p :: 2, a :: 2, b :: 2)

Assume that the three arguments are collinear, on some line L

* Return 1 if `a` and `b` are on the same side of `p` on L
* Return -1 if `a` and `b` are on different sides
* Return 0 if `a == p` or `b == p`.
"""
function sameside(p :: Tuple{Float64, Float64}, a :: Tuple{Float64, Float64}, b :: Tuple{Float64, Float64})
    if a[1] > p[1] && b[1] > p[1] || a[1] < p[1] && b[1] < p[1]
        return 1
    elseif a[1] < p[1] && b[1] > p[1] || a[1] > p[1] && b[1] < p[1]
        return -1
    elseif a[2] > p[2] && b[2] > p[2] || a[2] < p[2] && b[2] < p[2]
        return 1
    elseif a[2] < p[2] && b[2] > p[2] || a[2] > p[2] && b[2] < p[2]
        return -1
    else
        return 0
    end
end

function sameside(p, a, b)
    sameside(coord(p), coord(a), coord(b))
end

"""
    meet(p :: 2, q :: 2, a :: 2, b :: 2)

* Return 1 if the open line segments `(p, q)` and `(a, b)` meet in a single point.
* Return 0 if the the closed line segments `[p, a]` and `[a, b]` meet in one or several points.
* Return –1 otherwise.
"""
function meet(p, q, a, b)
    pqa = orient(p, q, a)
    pqb = orient(p, q, b)
    abp = orient(a, b, p)
    abq = orient(a, b, q)

    if opposite_signs(pqa, pqb) && opposite_signs(abp, abq)
        return 1
    elseif pqa & pqb == 1 || abp & abq == 1
        return -1
    elseif pqa == 0 && pqb == 0
        # all four points are collinear
        if sameside(p, a, b) == 1 && sameside(q, a, b) == 1 && sameside(a, p, q) == 1 && sameside(b, p, q) == 1
            return -1
        else
            return 0
        end
    else
        return 0
    end
end

@genpredicate function parallelorder(a :: 2, b :: 2, p :: 2, q :: 2)
    δ = b - a
    qp = q - p

    group!(δ...)
    group!(qp...)
    ext(δ, qp)
end
@doc """
    parallelorder(a :: 2, b :: 2, p :: 2, q :: 2) -> Int

Consider the oriented line defined by `a` and `b`
and the parallel lines passing through `p` and `q` respectively, with the same orientation.

* return 1 if the line passing through `p` is left of the line passing through `q`.
* return -1 in the reverse situation.
* return 0 if `a` and `b` are equal or if the parallel lines passing through `p` and `q` are equal.

This is a robust version of to `orient(b-a, q-p, 0)`.
Note also that `orient(a, b, c) == parallelorder(a, b, a, c)`.
""" parallelorder(:: NTuple{2, Float64}, :: NTuple{2, Float64}, :: NTuple{2, Float64}, :: NTuple{2, Float64})

"""
    rotation(pts :: AbstractVector)

Gives the number of rotation of the polygonal path defined by the elements of `pts`.
"""
function rotation(pts :: AbstractVector{T}) where T
    u = rand(SVector{2, Float64})
    origin = (0.0, 0.0)

    n = length(pts)
    @assert n >= 3

    pts = copy(pts)
    push!(pts, pts[1], pts[2])

    r = 0
    o1 = parallelorder(origin, u, pts[1], pts[2])
    for i in 1:n
        o2 = parallelorder(origin, u, pts[i+1], pts[i+2])
        if opposite_signs(o1, o2)
            ro = parallelorder(pts[i], pts[i+1], pts[i+1], pts[i+2])
            @assert ro != 0
            if ro > 0 && o1 < 0
                r += 1
            elseif ro < 0 && o1 > 0
                r -= 1
            end
        end
        o1 = o2
    end

    return r
end

@genpredicate function intersectorder(a :: 2, b :: 2, pa :: 2, pb :: 2, qa :: 2, qb :: 2)
    pp = pb - pa
    qq = qb - qa
    δ = b - a
    p0 = pa - a
    q0 = qa - a

    group!(pp...)
    group!(qq...)
    group!(δ...)
    group!(p0..., q0...)

    det(ext(p0, pp), ext(q0, qq),
        ext(δ, pp), ext(δ, qq))
end
@doc """
    intersectorder(a :: 2, b :: 2, pa :: 2, pb :: 2, qa :: 2, qb :: 2) -> Int

Consider the oriented line *L* defined by `a` and `b`, the line *P* defined by
`pa` and `pb` and the line *Q* defined by `qa` and `qb`.

Assumes that `parallelorder(a, b, pa, pb)` and `parallelorder(a, b, qa, qb)` have the same sign.
    Otherwise, the result has the opposite sign.

* return -1 if the intersection of *P* with *L* comes before the intersection of *Q* with *L*, following the orientation of *L*.
* return 1 in the reverse situation
* return 0 in case of equality or degeneracy.
""" intersectorder(:: NTuple{2, Float64}, :: NTuple{2, Float64}, :: NTuple{2, Float64}, :: NTuple{2, Float64}, :: NTuple{2, Float64}, :: NTuple{2, Float64})


text:
using DReal

a = [Var(Float64, -5.0, 5.0) for i = 1:14]
p = ForallVar(Float64, -2.0, 2.0)
r = ForallVar(Float64, -2.0, 2.0)
pest = ForallVar(Float64, -2.0, 2.0)
i = ForallVar(Float64, -2.0, 2.0)

c1 = 0.41328
c2 = 200.0
c3 = -0.366
c4 = 0.08979
c5 = -0.0337
c6 = 0.0001
c7 = 2.821
c8 = -0.05231
c9 = -0.05231
c9 = 0.10299
c10 = -0.00063
c11 = 1.0
c12 = 14.7
c13 = 0.9
c14 = 0.4
c15 = 0.4
c16 = 1.0
u = 23.0829

# Derivatives
pdot = c1 * ( 2*u*sqrt(p/c11 - (p/c11)^2) - (c3 + c4*c2*p + c5 * c2 *
  p^2 + c6*c2^2 * p))

rdot = 4 * ((c3 + c4*c2*p + c5*c2*p^2 + c6*c2^2*p)/
                (c13*(c3+c4*c2*pest^2 + c5*c2*pest^2 + c6*c2^2*pest)*
                     (1+i+c14*(r-c16)))-r)

pestdot = c1 * ((2*u*sqrt(p/c11 -
  (p/c11)^2))-c13*(c3+c4*c2*pest+c5*c2*pest^2+c6*c2^2*pest))

idot = c15*(r-c16)

V = (a[1]*p^2 + a[2]*r^2 + a[3]*pest^2 + a[4]*i^2 + a[5]*p*r +
a[6]*r*pest + a[7] * r*i + a[8] * pest*i + a[9]*p*i + a[10]*p +
a[11]*r + a[12]*pest + a[13]*i + a[14])^2

common = (a[14] + a[13]*i + a[4]*i^2 + a[10]*p +
   a[9]*i*p + a[1]*p^2 + a[12]*pest + a[8]*i*pest + a[3]*pest^2 + a[11]*r +
   a[7]*i*r + a[5]*p*r + a[6]*pest*r + a[2]*r^2)

# Partial Derivatives of V
dvbydp = 2 * (a[10] + a[9]*i + 2*a[1]*p + a[5]*r)
dvbydpest = 2*(a[12] + a[8]*i + 2*a[3]*pest + a[6]*r)
dvbydi = 2*(a[13] + 2*a[4]*i + a[9]*p + a[8]*pest + a[7]*r)
dvbydr = 2*(a[11] + a[7]*i + a[5]*p + a[6]*pest + 2*a[2]*r)

add!(common*(dvbydp * pdot + dvbydr * rdot + dvbydpest * pestdot + dvbydi * idot) < 0)
add!(V >= 0.0)
add!(p^2>0.001)
add!(r^2>0.001)
add!(pest^2 > 0.001)
add!(i^2 >0.001)


text:
using FastGaussQuadrature
using Polynomials

using DecFP: Dec128

function getCoefficientsGLRK1(T=Float64)
    a = reshape(Array([1//2]), (1,1))
    b = Vector([1//1])
    c = Vector([1//2])
    o = 2

    CoefficientsRK(T, :GLRK1, o, a, b, c)
end

function getCoefficientsGLRK2(T=Float64)
    a = @dec128 [
          [1/4       1/4-√3/6]
          [1/4+√3/6  1/4     ]
        ]
    b = @dec128 [1/2,      1/2     ]
    c = @dec128 [1/2-√3/6, 1/2+√3/6]
    o = 4

    CoefficientsRK(T, :GLRK2, o, a, b, c)
end

function getCoefficientsGLRK3(T=Float64)
    a = @dec128 [
          [5/36         2/9-√15/15  5/36-√15/30]
          [5/36+√15/24  2/9         5/36-√15/24]
          [5/36+√15/30  2/9+√15/15  5/36       ]
        ]
    b = @dec128 [5/18,        4/9,        5/18       ]
    c = @dec128 [1/2-√15/10,  1/2,        1/2+√15/10 ]
    o = 6

    CoefficientsRK(T, :GLRK3, o, a, b, c)
end

function getCoefficientsGLRK4(T=Float64)
    a = reshape(Array{Dec128}(@dec128 [
         [(-√30/144 + 1/8)
          (-13*√42/336 - √10/16 - √3/12 + √(-60*√30 + 450)/144 + √35/28 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√35/28 - √10/16 - √3/12 + √(-60*√30 + 450)/144 + 13*√42/336 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√35/14 - √(60*√30 + 450)/144 - √42/168 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
         [(-√(60*√30 + 450)/144 - √3/12 + √10/16 + √35/28 + 13*√42/336 + √(2*√30 + 15)/8)/√(2*√30 + 15)
          (√30/144 + 1/8)
          (-√35/14 + √42/168 + √(-60*√30 + 450)/144 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-13*√42/336 - √35/28 - √(60*√30 + 450)/144 - √3/12 + √10/16 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
         [(-√10/16 - √(60*√30 + 450)/144 + √3/12 + √35/28 + 13*√42/336 + √(2*√30 + 15)/8)/√(2*√30 + 15)
          (-√42/168 + √(-60*√30 + 450)/144 + √(-2*√30 + 15)/8 + √35/14)/√(-2*√30 + 15)
          (√30/144 + 1/8)
          (-13*√42/336 - √35/28 - √10/16 - √(60*√30 + 450)/144 + √3/12 + √(2*√30 + 15)/8)/√(2*√30 + 15)];
         [(-√(60*√30 + 450)/144 + √42/168 + √35/14 + √(2*√30 + 15)/8)/√(2*√30 + 15)
          (-13*√42/336 + √(-60*√30 + 450)/144 + √3/12 + √10/16 + √35/28 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√35/28 + √(-60*√30 + 450)/144 + √3/12 + √10/16 + 13*√42/336 + √(-2*√30 + 15)/8)/√(-2*√30 + 15)
          (-√30/144 + 1/8)]
        ]), (4,4))'

    b = @dec128 [
          49/(216+12*√30),
          49/(216-12*√30),
          49/(216-12*√30),
          49/(216+12*√30)
        ]

    c = @dec128 [
          1/2-√(3/7+2*√30/35)/2,
          1/2-√(3/7-2*√30/35)/2,
          1/2+√(3/7-2*√30/35)/2,
          1/2+√(3/7+2*√30/35)/2
        ]

    o = 8

    CoefficientsRK(T, :GLRK4, o, a, b, c)
end

function getCoefficientsGLRK5(T=Float64)
    a = reshape(Array{Dec128}(@dec128 [
         [-13*√70/3600 + 161/1800
          -59*√(20*√70 + 350)/19440 - 113*√(14*√70 + 245)/34020 - 11*√(-20*√70 + 350)/6480 + 4*√(-14*√70 + 245)/2835 + 13*√70/3600 + 161/1800
          -92*√(14*√70 + 245)/8505 + 4*√(20*√70 + 350)/1215 + 32/225
          -59*√(20*√70 + 350)/19440 - 113*√(14*√70 + 245)/34020 - 4*√(-14*√70 + 245)/2835 + 11*√(-20*√70 + 350)/6480 + 13*√70/3600 + 161/1800
          -2*√(14*√70 + 245)/315 - 13*√70/3600 + √(20*√70 + 350)/360 + 161/1800];
         [-113*√(-14*√70 + 245)/34020 - 13*√70/3600 + 4*√(14*√70 + 245)/2835 + 11*√(20*√70 + 350)/6480 + 59*√(-20*√70 + 350)/19440 + 161/1800
          13*√70/3600 + 161/1800
          -92*√(-14*√70 + 245)/8505 - 4*√(-20*√70 + 350)/1215 + 32/225
          -2*√(-14*√70 + 245)/315 - √(-20*√70 + 350)/360 + 13*√70/3600 + 161/1800
          -11*√(20*√70 + 350)/6480 - 113*√(-14*√70 + 245)/34020 - 13*√70/3600 - 4*√(14*√70 + 245)/2835 + 59*√(-20*√70 + 350)/19440 + 161/1800];
         [-23*√(20*√70 + 350)/5760 - 13*√70/3600 + 161/1800 + 11*√(14*√70 + 245)/1440
          13*√70/3600 + 23*√(-20*√70 + 350)/5760 + 11*√(-14*√70 + 245)/1440 + 161/1800
          32/225
          -11*√(-14*√70 + 245)/1440 - 23*√(-20*√70 + 350)/5760 + 13*√70/3600 + 161/1800
          -11*√(14*√70 + 245)/1440 - 13*√70/3600 + 161/1800 + 23*√(20*√70 + 350)/5760];
         [-59*√(-20*√70 + 350)/19440 - 13*√70/3600 + 4*√(14*√70 + 245)/2835 + 113*√(-14*√70 + 245)/34020 + 11*√(20*√70 + 350)/6480 + 161/1800
          13*√70/3600 + √(-20*√70 + 350)/360 + 2*√(-14*√70 + 245)/315 + 161/1800
          4*√(-20*√70 + 350)/1215 + 92*√(-14*√70 + 245)/8505 + 32/225
          13*√70/3600 + 161/1800
          -59*√(-20*√70 + 350)/19440 - 11*√(20*√70 + 350)/6480 - 13*√70/3600 - 4*√(14*√70 + 245)/2835 + 113*√(-14*√70 + 245)/34020 + 161/1800];
         [-√(20*√70 + 350)/360 - 13*√70/3600 + 161/1800 + 2*√(14*√70 + 245)/315
          -11*√(-20*√70 + 350)/6480 + 4*√(-14*√70 + 245)/2835 + 13*√70/3600 + 113*√(14*√70 + 245)/34020 + 59*√(20*√70 + 350)/19440 + 161/1800
          -4*√(20*√70 + 350)/1215 + 32/225 + 92*√(14*√70 + 245)/8505
          -4*√(-14*√70 + 245)/2835 + 11*√(-20*√70 + 350)/6480 + 13*√70/3600 + 113*√(14*√70 + 245)/34020 + 59*√(20*√70 + 350)/19440 + 161/1800
          -13*√70/3600 + 161/1800]
        ]), (5,5))'

    b = @dec128 [
          1/(4600/729+1300*√70/5103),
          1/(4600/729-1300*√70/5103),
          64/225,
          1/(4600/729-1300*√70/5103),
          1/(4600/729+1300*√70/5103)
        ]

    c = @dec128 [
          1/2-√(5/9+2*√70/63)/2,
          1/2-√(5/9-2*√70/63)/2,
          1/2,
          1/2+√(5/9-2*√70/63)/2,
          1/2+√(5/9+2*√70/63)/2
        ]

    o = 10

    CoefficientsRK(T, :GLRK5, o, a, b, c)
end

function getCoefficientsGLRK6(T=Float64)
    a = reshape(Array{Dec128}(@dec128 [
            ["+0.0428311230947925862600740355431832233817056253710109955996588599497364402",
             "-0.0147637259971974124753725910605206514420001878003408840631515394063317448",
             "+0.00932505070647775119143888450800314858828816322281640330116548464376180630",
             "-0.00566885804948351190092125641621650656214351346713436514621339322133564885",
             "+0.00285443331509933513092928583011602153367099528247654439080898011903375405",
             "-0.000812780171264762112299135651562540066903951464973606518543200348171357290"];
            ["+0.0926734914303788631865122917633203161433485534502333679517368248345204102",
             "+0.0901903932620346518924583784594290279153804731866863705724348100592850095",
             "-0.0203001022932395859524940805242724601067259796418104502230488160090068057",
             "+0.0103631562402464237307199458065599778725040265335003155034035217850387283",
             "-0.00488719292803767146341420376578964407137592740187386736001668777722125914",
             "+0.00135556105548506177551787075080010874364457167567922814822699418129216854"];
            ["+0.0822479226128438738077716511411289215554396236655129640997987178532804788",
             "+0.196032162333245006055759781563801382788760315019545751706035725357898777",
             "+0.116978483643172761847467585997387748702913901442302633827906329990978550",
             "-0.0204825277456560976298590118654006438219872235078237396303198827234514190",
             "+0.00798999189966233579720442148030827079362758194981571971201211474213850758",
             "-0.00207562578486633419359528915758164772805951363936400480613058750798236207"];
            ["+0.0877378719744515067137433602439480944914707643813859960054483074074552416",
             "+0.172390794624406967987712335438549785037133364423557021432857505376431512",
             "+0.254439495032001621324794183860176141227815026392429007286132542705408521",
             "+0.116978483643172761847467585997387748702913901442302633827906329990978551",
             "-0.0156513758091757022708430246449433269579993686461730105611661052393287579",
             "+0.00341432357674129871237641994523752520797162707650902709951900204619240137"];
            ["+0.0843066851341001107446302003355663380197666790663427630510907257181807085",
             "+0.185267979452106975248330960684647699902136873775246608504886307895791278",
             "+0.223593811046099099964215226188215519533323776351104952152409138196918370",
             "+0.254257069579585109647429252519047957512553782526415717878861475990963909",
             "+0.0901903932620346518924583784594290279153804731866863705724348100592850092",
             "-0.00701124524079369066636422067695386937993730270821137675241910493504753033"];
            ["+0.0864750263608499346324472067379289868303152022069955977178609202476442335",
             "+0.177526353208969968653987471088742034297089951090896196754060639999536273",
             "+0.239625825335829035595856428410992003967971316351739632802026053203292742",
             "+0.224631916579867772503496287486772348817539639661788864354647175338195299",
             "+0.195144512521266716260289347979378707272761134173713625208021159524901768",
             "+0.0428311230947925862600740355431832233817056253710109955996588599497364397"]
        ]), (6,6))'

    b = @dec128 [
            "0.0856622461895851725201480710863664467634112507420219911993177198994728803",
            "0.180380786524069303784916756918858055830760946373372741144869620118570019",
            "0.233956967286345523694935171994775497405827802884605267655812659981957101",
            "0.233956967286345523694935171994775497405827802884605267655812659981957101",
            "0.180380786524069303784916756918858055830760946373372741144869620118570019",
            "0.0856622461895851725201480710863664467634112507420219911993177198994728804"
        ]

    c = @dec128 [
            "0.0337652428984239860938492227530026954326171311438550875637251917366932495891",
            "0.169395306766867743169300202490047326496775717802414964592736647073908251704",
            "0.380690406958401545684749139159644032290694684929989324909302417712862532843",
            "0.619309593041598454315250860840355967709305315070010675090697582287137467157",
            "0.830604693233132256830699797509952673503224282197585035407263352926091748261",
            "0.966234757101576013906150777246997304567382868856144912436274808263306750376"
        ]

    o = 12

    CoefficientsRK(T, :GLRK6, o, a, b, c)
end

function getCoefficientsGLRK(s::Int; T=Float64, high_precision=true)

    if high_precision == true
        if s == 1
            return getCoefficientsGLRK1(T)
        elseif s == 2
            return getCoefficientsGLRK2(T)
        elseif s == 3
            return getCoefficientsGLRK3(T)
        elseif s == 4
            return getCoefficientsGLRK4(T)
        elseif s == 5
            return getCoefficientsGLRK5(T)
        elseif s == 6
            return getCoefficientsGLRK6(T)
        end
    end

    # order
    o = 2s

    # obtain Gauss-Legendre nodes and weights
    q = GaussLegendreQuadrature(s, T)

    # create Lagrange polynomial
    vdm = vandermonde_matrix_inverse(nodes(q))

    # compute monomial basis functions and corresponding integrals
    poly_ints = []
    for i in 1:s
        y = zeros(s)
        y[i] = 1
        mon = *(vdm, y)
        push!(poly_ints, polyint(Poly(mon)))
    end

    # compute Runge-Kutta coefficients
    a = zeros(s,s)
    for i in 1:s
        for j in 1:s
            a[i,j] = poly_ints[j](nodes(q)[i])
        end
    end

    CoefficientsRK(Symbol("GLRK", s), o, a, weights(q), nodes(q))
end


text:
if basename(pwd()) == "aoc"
    cd("2017/19")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using MeshCat
using ConstrainedDynamicsVis

convertvideo("test_files/test.tar")
@test ispath("test_files/test.mp4")
rm("test_files/test.mp4")
@test true


text:
using Documenter, MakieLayout, Makie, Animations

# don't open windows while generating animations
Makie.AbstractPlotting.inline!(true)

makedocs(
    sitename="MakieLayout.jl",
    pages = [
        "index.md",
        "GridLayout" => "grids.md",
        "LAxis" => "laxis.md",
        "How layouting works" => "layouting.md"
    ],
    format = Documenter.HTML(
            prettyurls = get(ENV, "CI", nothing) == "true"
        )
    )

# deploydocs(
#     repo = "github.com/jkrumbiegel/Animations.jl.git",
# )


text:
@testset "BSplineSpace" begin
    Random.seed!(42)

    P1 = BSplineSpace(2, Knots([1, 3, 5, 6, 8, 9]))
    @test bsplinesupport(P1, 2) == 3..8
    @test dim(P1) == 3
    @test properdim(P1) == 3
    @test isproper(P1) == true

    P2 = BSplineSpace(1, Knots([1, 3, 3, 3, 8, 9]))
    @test isproper(P2) == false
    @test properdim(P2) == 3

    i = 2
    k = Knots([5, 12, 13, 13, 14])
    p = 2
    P = BSplineSpace(p, k)
    @test bsplinesupport(P) == [5..13, 12..14]
    @test bsplinesupport(P, i) == 12..14

    @test isproper(BSplineSpace(2, Knots([1, 3, 5, 6, 8, 9])))
    @test !isproper(BSplineSpace(1, Knots([1, 3, 3, 3, 8, 9])))

    @test dim(BSplineSpace(2, Knots([1, 3, 5, 6, 8, 9]))) == 3

    P1 = BSplineSpace(1, Knots([1, 3, 5, 8]))
    P2 = BSplineSpace(1, Knots([1, 3, 5, 6, 8, 9]))
    P3 = BSplineSpace(2, Knots([1, 1, 3, 3, 5, 5, 8, 8]))
    @test P1 ⊆ P2
    @test P1 ⊆ P3
    @test P2 ⊈ P3

    n1, n2, n3 = dim(P1), dim(P2), dim(P3)
    A12 = changebasis(P1, P2)
    A13 = changebasis(P1, P3)
    Δ12 = [bsplinebasis(P1, i, t) - sum(A12[i, j] * bsplinebasis(P2, j, t) for j in 1:n2) for i in 1:n1, t in 8 * rand(10)]
    Δ13 = [bsplinebasis(P1, i, t) - sum(A13[i, j] * bsplinebasis(P3, j, t) for j in 1:n3) for i in 1:n1, t in 8 * rand(10)]
    @test norm(Δ12) < 1e-14
    @test norm(Δ13) < 1e-14

    P4 = BSplineSpace(1, Knots([1, 2, 3, 4, 5]))
    P5 = BSplineSpace(2, Knots([-1, 0.3, 2, 3, 3, 4, 5.2, 6]))
    @test P4 ⊑ P4
    @test P4 ⊑ P5
    @test P5 ⊒ P4
    @test P5 ⋢ P4
    @test P4 ⋣ P5

    P4_ = BSplineSpace(degree(P4) - 1, knots(P4)[2:end-1])
    P5_ = BSplineSpace(degree(P5) - 1, knots(P5)[2:end-1])
    @test P4_ ⊑ P5_

    n4, n5 = dim(P4), dim(P5)
    A45 = changebasis(P4, P5)
    Δ45 = [bsplinebasis(P4, i, t) - sum(A45[i, j] * bsplinebasis(P5, j, t) for j in 1:n5) for i in 1:n1, t in 2 * rand(10) .+ 2]
    @test norm(Δ45) < 1e-14
end


text:
# Interpolation
import Interpolations: interpolate

"""
    interpolate(xval:: Real, c1:: Curve)

Returns the interpolated or extrapolated y-value of the curve for a given x-value.
"""
function interpolate(xval:: Real, c1:: Curve)
    if isnothing(c1.etp)
        return c1.y[begin]
    else
        c = c1.etp(c1.logx ? log(xval) : xval)
        return c1.logy ? exp(c) : c
    end
end

interpolate(xval:: Tenor, c1:: Curve) = interpolate(get_days(xval), c1)

interpolate(xval:: AbstractString, c1:: Curve) = interpolate(Tenor(xval), c1)

"""
    interpolate(xval:: AbstractArray{T} where T, c1:: Curve; kwargs...):: Curve

Interpolates the curve onto the given array and returns a new Curve instance on the interpolated points.

The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
interpolate(xval:: AbstractArray{T} where T, c1:: Curve; kwargs...):: Curve = Curve(xval, interpolate.(xval, c1); kwargs...)

"""
    interpolate(c0:: Curve, c1:: Curve; kwargs...):: Curve

Interpolates the Curve ˋc1ˋ on the x-axis points of the Curve ˋc0ˋ and returns a new Curve instance on the interpolated
points.

The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
interpolate(c0:: Curve, c1:: Curve; kwargs...):: Curve = interpolate(c0.x, c1; kwargs...)

# Basic operations

import Base: +, -, *, /, ^, exp, log, sin, cos, tan, length, ==, ≈

length(c1:: Curve):: Int = length(c1.y)

# comparisons
==(c1:: Curve, c2:: Curve) = c1.x == c2.x && c1.y == c2.y && c1.etp == c2.etp &&
    c1.logx == c2.logx && c1.logy == c2.logy

≈(c1:: Curve, c2:: Curve) = c1.x ≈ c2.x && c1.y ≈ c2.y && c1.etp ≈ c2.etp &&
    c1.logx == c2.logx && c1.logy == c2.logy

#=
For all calculations on curves, the curve settings (interpolator, extrapolator, log axis settings) are taken as default
unless specified otherwise (they are not copied from input curves). This is to avoid issues e.g. with log-scales
on potentially negative values.
=#

operations = (:+, :-, :*, :/, :^)
for op in operations
    # Define operations with scalars
    @eval $op(c1:: Curve, a:: Number; kwargs...) = Curve(c1.x, $op.(c1.y, a); sort=false, kwargs...)
    @eval $op(a:: Number, c1:: Curve; kwargs...) = Curve(c1.x, $op.(a, c1.y); sort=false, kwargs...)

    # Operations with Multiple Curves
    @eval function $op(c1:: Curve, c2:: Curve; kwargs...)
        if c1.x == c2.x # fast mode if no interpolation is required
            x = c1.x
            y = $op.(c1.y, c2.y)
        else
            y_c1_grid = $op.(c1.y, interpolate.(c1.x, c2))
            y_c2_grid = $op.(interpolate.(c2.x, c1), c2.y)
            x, y = mergexy(c1.x, y_c1_grid, c2.x, y_c2_grid)
            x, y = uniquexy(x, y)
        end
        Curve(x, y; sort=true, kwargs...)
    end

    # propagate missing
    @eval $op(c1:: Curve, c2:: Missing; kwargs...) = missing
    @eval $op(c1:: Missing, c2:: Curve; kwargs...) = missing
end

# Define operations where Curve is the only argument
operations = (:exp, :log, :sin, :cos, :tan)
for op in operations
    @eval $op(c1:: Curve; kwargs...) = Curve(c1.x, $op.(c1.y); sort=false, kwargs...)
end

# first and last point functions

"""
    firstpoint(c1:: Curve; dims=1):: Real

Returns the first point of the X-axis (if `dims=1`) or the Y-axis (if `dims=2`).
"""
function firstpoint(c1:: Curve; dims=1):: Real
    if dims == 1
        return c1.x[begin]
    elseif dims == 2
        return c1.y[begin]
    else
        error("invalid dimension for Curve objects, only 1 and 2 supported")
    end
end

"""
    lastpoint(c1:: Curve; dims=1):: Real

Returns the last point of the X-axis (if `dims=1`) or the Y-axis (if `dims=2`).
"""
function lastpoint(c1:: Curve; dims=1):: Real
    if dims == 1
        return c1.x[end]
    elseif dims == 2
        return c1.y[end]
    else
        error("invalid dimension for Curve objects, only 1 and 2 supported")
    end
end

import Base: first, last, filter

"""
    first(c1:: Curve, n:: Integer; kwargs...)

Returns a new curve containing the first n points of the previous curve.
`n` must be at least 1.

The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
function first(c1:: Curve, n:: Integer; kwargs...)
    n < 1 && error("`n` must be at least 1 for definition of a curve")
    n = min(n, length(c1))
    Curve(c1.x[begin:n], c1.y[begin:n]; sort=false, kwargs...)
end

"""
    last(c1:: Curve, n:: Integer; kwargs...)

Returns a new curve containing the last n points of the previous curve.
`n` must be at least 1.

The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
function last(c1:: Curve, n:: Integer; kwargs...)
    n < 1 && error("`n` must be at least 1 for definition of a curve")
    n = min(n, length(c1))
    Curve(c1.x[end-n+1:end], c1.y[end-n+1:end]; sort=false, kwargs...)
end

"""
    filter(f:: Function, c1::Curve; axis:: Symbol = :x, kwargs...)

Filters curve points according to the given (boolean return) function.

The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
function filter(f:: Function, c1::Curve; axis:: Symbol = :x, kwargs...)
    if axis==:x
        mask = f.(c1.x)
    elseif axis==:y
        mask = f.(c1.y)
    else
        error("axis must be :x or :y, the value $axis is not allowed")
    end
    count(mask) < 1 && error("less than 1 point remaining")
    Curve(c1.x[mask], c1.y[mask]; sort=false, kwargs...)
end

# Helper functions for concatination

"""
Removes duplicates in the array x and the entries in the array y which have the same index as the x-duplicates.

Note that it does not check if the y-values for x-value duplicates are the same!
"""
@inline function uniquexy(x:: AbstractArray, y:: AbstractArray)
    @inbounds begin
        idx = unique(z -> x[z], 1:length(x))
        x[idx], y[idx]
    end
end

"""
Merges two pairs of arrays x, y into a combined (along x-values) array pair.
"""
mergexy(x1:: AbstractArray, y1:: AbstractArray, x2:: AbstractArray, y2:: AbstractArray) = vcat(x1, x2), vcat(y1, y2)

"""
    concat(c1:: Curve, c2:: Curve; drop_dup=true, kwargs...)

Merges two curves.

Element type is inferred by promotion. Dulicate points are dropped by default, unless ˋdrop_dup=falseˋ is set.
The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
function concat(c1:: Curve, c2:: Curve; kwargs...)
    x_all, y_all = mergexy(c1.x, c1.y, c2.x, c2.y)
    x_all, y_all = uniquexy(x_all, y_all)
    return Curve(x_all, y_all; sort=true, kwargs...)
end

concat(c1:: Curve, cx:: Curve...) = concat(c1, concat(cx[begin], cx[begin+1:end]...))

# Applying Functions

"""
    apply(f:: Function, c1:: Curve; axis:: Symbol = :xy, kwargs...)

If ˋaxisˋ is ˋ:xyˋ (default): applies a 2-argument function ˋf(x,y)=zˋ to each entry of the Curve.

If ˋaxisˋ is ˋ:xˋ or ˋ:yˋ: applies a 1-argument function ˋf(x)=zˋ to a single axis of the Curve.

The output Curve is constructed using default settings, alternative settings can be passed to the Curve constructor
using the ˋkwargs...ˋ
"""
function apply(f:: Function, c1:: Curve; axis:: Symbol = :xy, kwargs...)
    if axis==:xy
        x = c1.x
        y = f.(c1.x, c1.y)
        return Curve(x, y; sort=false, kwargs...)
    elseif axis==:x
        x = f.(c1.x)
        y = c1.y
        return Curve(x, y; sort=true, kwargs...)
    elseif axis==:y
        x = c1.x
        y = f.(c1.y)
        return Curve(x, y; sort=false, kwargs...)
    else
        error("axis must be :xy, :x or :y, the value $axis is not allowed")
    end
end


text:
function do_unsupported()
    a = 1
    Dict(key + 1 => value + 1 for (key, value) in Dict())
    b = Bool(a)
    println(b ? ("True") : ("False"))
end

if abspath(PROGRAM_FILE) == @__FILE__
    do_unsupported()
end


text:
function sum_values(x,y)
    return x+y
end


text:
### A bunch of helper functions to apply stuff with environments to mps
function _applymps0site(Λ::Matrix{T},
                        envL::Array{T, 3},
                        envR::Array{T, 3}) where {T<:Number}
    @tensor A[d,d'] := envL[u,m,d] * Λ[u, u'] * envR[u',m,d']
    A
end

function _mpsupdateright(renv::Array{T, 3},
                         mat::Array{T, 3},
                         hmpo::Array{T,4}) where {T<:Number}
    @tensor R[u,m,d] :=
        ((renv[u',m',d'] * mat[u,o,u']) * hmpo[m,o',m',o]) * conj(mat)[d,o',d']
    R
end

function _mpsupdateright(renv::SymTensor{S,T,3},
                         mat::SymTensor{S,T,3},
                         hmpo::SymTensor{S,T,4}) where {S,T}

    R = contract(contract(contract(renv, (-1,3,4), mat, (1,2,-1)),
                          (1, -1, -2, 4), hmpo, (3,2,-2,-1)),
                 (1,-1,2,-2), dual(mat), (3,-1,-2))
    R
end

# function _mpsupdateright(renv::SymTensor{Tv, 4},
#                          mat::SymTensor{Tv, 3},
#                          hmpo::SymTensor{Tv, 4}) where {Tv<:Number}

#     # println(renv.legs)
#     # println(mat.legs)
#     # println(hmpo.legs)
#     # renv is W_l x M_d x W_r x M_u
#     R = contract(contract(contract(renv, (1,2,3,-1), mat, (4,5,-1)),
#                           (-1,3,4,5,-2), hmpo, (1,2,-1,-2)),
#                  (1, -1,-2,3,4), invlegs(conj(mat)), (2, -1,-2))
#     R
# end

function _mpsupdateleft(lenv::Array{T, 3},
                        mat::Array{T, 3},
                        hmpo::Array{T,4}) where {T<:Number}
    @tensor L[u,m,d] :=
        ((lenv[u',m',d'] * mat[u',o,u]) * hmpo[m',o',m,o]) * conj(mat)[d',o',d]
    L
end

function _mpsupdateleft(lenv::SymTensor{S,T, 3},
                        mat::SymTensor{S,T, 3},
                        hmpo::SymTensor{S,T, 4}) where {S,T}

    L = contract(contract(contract(lenv, (-1,3,4), mat, (-1,2,1)),
                          (1, -1, -2, 4), hmpo, (-2,2,3,-1)),
                 (1,-1,2,-2), dual(mat), (-2,-1,3))

    L
end

# function _mpsupdateleft(lenv::SymTensor{Tv, 4},
#                         mat::SymTensor{Tv, 3},
#                         hmpo::SymTensor{Tv, 4}) where {Tv<:Number}

#     L = contract(contract(contract(lenv, (1,2,3,-1), mat, (-1,4,5)),
#                           (1,2,-1,-2,5), hmpo, (-1,3,4,-2)),
#                  (1,-1,-2,3,4), invlegs(conj(mat)), (-1,-2,2))
#     L
# end

function _applymps1site(v,
                        envL::Array{T,3},
                        envR::Array{T,3},
                        hmpo::Array{T,4}) where {T}

    @tensor v[l,o,r] := ((envL[l',ml,l] * v[l',o',r']) *
                         hmpo[ml,o,mr,o']) * envR[r',mr,r]
    v
end

function _applymps2site(v,
                        envL::Array{T,3},
                        envR::Array{T,3},
                        hmpoL::Array{T,4},
                        hmpoR::Array{T,4}) where {T}

    @tensor v[l,o1,o2,r] := (((envL[l',ml,l] * v[l',o1',o2',r']) *
                              hmpoL[ml,o1,mm,o1']) *
                             hmpoR[mm,o2,mr,o2']) *
                             envR[r',mr,r]
    v
end

function _applymps2site(v::SymTensor{S,T,4},
                        envL::SymTensor{S,T,3},
                        envR::SymTensor{S,T,3},
                        hmpoL::SymTensor{S,T,4},
                        hmpoR::SymTensor{S,T,4}) where {S,T}

    v = contract(contract(contract(contract(
        envL, (-1,2,1), v, (-1,3,4,5)),
                                   (1,-1,-2,4,5), hmpoL, (-1,2,3,-2)),
                          (1,2,-1,-2,5), hmpoR, (-1,3,4,-2)),
                 (1,2,3,-1,-2), envR, (-2,-1,4))
    v
end

# function _applymps2site(v::SymTensor{Tv,4},
#                         envL::SymTensor{Tv,4},
#                         envR::SymTensor{Tv,4},
#                         hmpoL::SymTensor{Tv,4},
#                         hmpoR::SymTensor{Tv,4}) where {Tv<:Number}

#     v = contract(contract(contract(contract(
#         envL, (1,2,3,-1), v, (-1,4,5,6)),
#                                    (1,2,-1,-2,5,6), hmpoL, (-1,3,4,-2)),
#                           (1,2,3,-1,-2,6), hmpoR, (-1,4,5,-2)),
#                  (-1,1,2,3,-2,-3), envR, (-2,4,-1,-3))
#     v
# end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
macro linklibrary_f2p(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:f_to_p, n, :e, es))

  inputtype = n == 32 ? Float64 : Float32
  outputtype = Symbol(:UInt, n)

  esc(quote
    f2p = (v) -> ccall( ($fnname, $path), $outputtype, ($inputtype,), v)
  end)
end
#spot unit testing

@linklibrary_f2p(8,0)
@test f2p(-32.5f0) == 0x82
@test f2p(-48.0f0) == 0x82
@test reinterpret(UInt8, Posit{8,0}(-48.0f0)) == 0x82

@linklibrary_f2p(16, 0)
@test f2p(0.0f0)  == 0x0000
@test f2p(Inf32)  == 0x8000
@test f2p(1.0f0)  == 0x4000
@test f2p(-1.0f0) == 0xc000
@test f2p(2.0f0)  == 0b0110_0000_0000_0000
@test f2p(1.5f0)  == 0b0101_0000_0000_0000

@linklibrary_f2p(16, 1)
@test f2p(0.0f0)  == 0x0000
@test f2p(Inf32)  == 0x8000
@test f2p(-1.0f0) == 0xc000
@test f2p(1.0f0)  == 0x4000
@test f2p(4.0f0)  == 0x6000
@test f2p(2.0f0)  == 0x5000
@test f2p(3.0f0)  == 0x5800
@test f2p(0.25f0) == 0x2000
@test f2p(0.50f0) == 0x3000

@linklibrary_f2p(16, 2)
@test f2p(1.0f0)  == 0x4000
@test f2p(16.0f0) == 0x6000
@test f2p(2.0f0)  == 0x4800
@test f2p(4.0f0)  == 0x5000
@test f2p(8.0f0)  == 0x5800

@linklibrary_f2p(32, 2)
@test f2p(1.0)  == 0x4000_0000
@test f2p(1.5)  == 0x4400_0000
@test f2p(2.0)  == 0x4800_0000
@test f2p(4.0)  == 0x5000_0000
@test f2p(8.0)  == 0x5800_0000
@test f2p(16.0) == 0x6000_0000

#comprehensive testing on 8,0;8,1;8,2
@linklibrary_p2f(8,0)
@linklibrary_f2p(8,0)
for idx=0x00:0xFF
  @test f2p(p2f(idx)) == idx
end
@linklibrary_p2f(8,1)
@linklibrary_f2p(8,1)
for idx=0x00:0xFF
  @test f2p(p2f(idx)) == idx
end
@linklibrary_p2f(8,2)
@linklibrary_f2p(8,2)
for idx=0x00:0xFF
  @test f2p(p2f(idx)) == idx
end

#comprehensive testing on 16,0; 16,1; 16,2
@linklibrary_p2f(16,0)
@linklibrary_f2p(16,0)
for idx=0x0000:0xFFFF
  @test f2p(p2f(idx)) == idx
end
@linklibrary_p2f(16,1)
@linklibrary_f2p(16,1)
for idx=0x0000:0xFFFF
  @test f2p(p2f(idx)) == idx
end
@linklibrary_p2f(16,2)
@linklibrary_f2p(16,2)
for idx=0x0000:0xFFFF
  @test f2p(p2f(idx)) == idx
end

#randomized testing on 32,0; 32,1; 32,2; 32,3
@linklibrary_p2f(32,0)
@linklibrary_f2p(32,0)
for idx=1:5000
  v = rand(UInt32)
  @test f2p(p2f(v)) == v
end
@linklibrary_p2f(32,1)
@linklibrary_f2p(32,1)
for idx=1:5000
  v = rand(UInt32)
  @test f2p(p2f(v)) == v
end
@linklibrary_p2f(32,2)
@linklibrary_f2p(32,2)
for idx=1:5000
  v = rand(UInt32)
  @test f2p(p2f(v)) == v
end
@linklibrary_p2f(32,3)
@linklibrary_f2p(32,3)
for idx=1:5000
  v = rand(UInt32)
  @test f2p(p2f(v)) == v
end

#spot test conversions of extreme values

@linklibrary_f2p(8,0)
#overflow goes to maxreal
@test f2p(prevfloat(Inf32))  == 0x7F
@test f2p(nextfloat(-Inf32)) == 0x81
#underflow doesn't go to zero
@test f2p(eps(0.0f0))        == 0x01
@test f2p(-eps(0.0f0))       == 0xFF


text:
# Creates a table of forecast times based on 1000 posterior draws
# run for 2600 days created in ttv_posterior.jl

using JLD2; using DelimitedFiles; using Printf; using Statistics

@load "../../data/T1_timing_posterior.jld2" n nplanet t0 h rstar tt1 count1 ngrid elements_grid tt_grid

# Now, try to improve the fit to the data:
# First, find the indices that match for each planet:
data_out = zeros(sum(count1[2:end]),4)
# Create a timings table:
io = open("timings_forecast.tex","w")
pname = ["b","c","d","e","f","g","h"]
itrans = 1
for ip=1:nplanet
  ittv = 1
  for j=1:count1[ip+1]
    # First, match the observation with a transit time:
    # Now compute the mean and standard deviation of the posterior:
    data_out[itrans,1] = ip
    data_out[itrans,2] = j
    data_out[itrans,3] = mean(tt_grid[ip+1,j,:])
    data_out[itrans,4] = std(tt_grid[ip+1,j,:])
    dstring = string(pname[ip]," & ",@sprintf("%i",j)," & ",@sprintf("%12.6f",data_out[itrans,3])," & ",@sprintf("%9.6f",data_out[itrans,4]),"\\cr")
    println(io,dstring)
    # Now, create output for timings table:
    global itrans += 1
  end
end
close(io)
writedlm("times_forecast.txt",data_out,',')


text:
using .Cornerback
using GLMakie

# Receiver runs the chosen route.

# x, y = fade(x, y, speed, δt)
# x, y = post(x, y, speed, δt)

function make_figure(obs)
    # Create figure
    fig = Figure(); display(fig)
    ax = Axis(fig[1, 1])

    # Plot observables
    scatter!(ax, obs)

    # Plot static parameters
    ax.aspect = DataAspect()
    xlims!(ax, -1, 20)
    ylims!(ax, -1, 40)

    return fig, ax
end

function main()
    # Receiver parameters
    speed = 10.0
    δt = 0.01
    x = 0.0
    y = 0.0
    pos = Observable(Point2f(x, y))

    fig, ax = make_figure(pos)

    for i in 1:10000
        x, y = post(x, y, speed, δt; cut = 10.0, angle = 60.0)
        pos[] = Point2f(x, y)
        sleep(0.01)
    end

end

main()


text:
# -*- coding: utf-8 -*-

include("dx9adl_demo_0001.jl")

println("Testing DirectX")
d9 = DirectX.connect(64 * 16, 64 * 9)
@test d9 != nothing
@test 0 == Demo_0001.demo_0001(d9)
@test DirectX.close(d9)
# @test ture != false
@test false != true
println("ok")


text:
using MarketData, TradeModels, Quandl, Distributions
import Base.show, TimeSeries.findall

tinv(df::Int, pval::Float64) = quantile(TDist(df), pval) # t table
tinv(df::Int, pval::FloatRange{Float64}) = quantile(TDist(df), pval) # t table

include(Pkg.dir("FinancialSeries/src/dev/byvolume.jl"))

crude = quandl("CHRIS/CME_CL1", rows=1000)
myvol = 825399
c = byvolume(crude, myvol)

function findall(b::BitArray)
    rownums = zeros(Int, sum(b))
    j = 1
    for i in 1:length(b)
      if b[i]
        rownums[j] = i
        j+=1
      end
    end
    rownums
end

# Phi = phi_0 to phi_p or length of p + 1
# y = all the ys starting at the p index, because we need pre values
# [yp+1 yp+2 ... end]
# the number of rows in y == number of rows in X
# X is a big matrix with the first row of 1s and the number of columns after that equal to p for a total of p+1 columns
# the number of columns in X must match the number of rows in Phi, because that's how it lines up
# to estimate Phi  :
##     inv(X'X) * X' * y

ret = percentchange(cl, method="log")
ret = flipud(ret.values)

#myfit = ar(c["Settle"],6)


text:
module Watchtower

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Watchtower" => Ahorn.EntityPlacement(
        Maple.Towerviewer
    ),
)

Ahorn.nodeLimits(entity::Maple.Towerviewer) = 0, -1

sprite = "objects/lookout/lookout05.png"

function Ahorn.selection(entity::Maple.Towerviewer)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny, jx=0.5, jy=1.0))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Towerviewer)
    px, py = Ahorn.position(entity)

    for node in get(entity.data, "nodes", [])
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py - 8, nx, ny - 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny, jx=0.5, jy=1.0)

        px, py = nx, ny
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Towerviewer, room::Maple.Room)
    x, y = Ahorn.position(entity)

    Ahorn.drawSprite(ctx, sprite, x, y, jx=0.5, jy=1.0)
end

end


text:
# Some track and camera properties
    f = 40e-3
    # pixelspermeter = 181818.0

    trackprop = TrackProperties(track_gauge = 1.435)

    # horizontal curvature
    β = 4e-6

    # vertical curvature
    α = -3.2e-6

    left_track  = [
        [x, -trackprop.track_gauge/2.0 + β*x^2, α*x^2]
    for x = LinRange(4.0,40.0,40)];

    right_track = [
        [x, trackprop.track_gauge/2.0 + β*x^2, α*x^2]
    for x = LinRange(4.0,40.0,40)];

    cameraposition_reference = [0.0,-0.1,-2.4];
    ψθφ_ref = [-1.0,-5.0,1.0] .* (pi/180.0);
    camera_reference = VideoCamera(cameraposition_reference;
        ψθφ = ψθφ_ref,
        focal_length = f
        # , pixelspermeter = pixelspermeter
    )

@testset "mapping track to image" begin

    # test two different ways to calculate the us for a straight track
    left_track_ref  = [
        [x, -trackprop.track_gauge/2.0, 0.0]
    for x = LinRange(4.0,40.0,40)]

    right_track_ref = [
        [x, trackprop.track_gauge/2.0, 0.0]
    for x = LinRange(4.0,40.0,40)]

    left_ref_uvs = camera_image(camera_reference, left_track_ref);
    right_ref_uvs = camera_image(camera_reference, right_track_ref);

    uLs_ref = [uv[1] for uv in left_ref_uvs];
    vLs_ref = [uv[2] for uv in left_ref_uvs];
    uRs_ref = [uv[1] for uv in right_ref_uvs];
    vRs_ref = [uv[2] for uv in right_ref_uvs];

    # plot(uLs_ref, vLs_ref, yflip = true)
    # plot!(uRs_ref, vRs_ref, yflip = true)

    uLs0 = left_track_image_u(camera_reference, trackprop, vLs_ref);
    uRs0 = right_track_image_u(camera_reference, trackprop, vRs_ref);

    @test maximum(abs.(uLs0 - uLs_ref) ./ abs.(uLs_ref)) < 1e-12
    @test maximum(abs.(uRs0 - uRs_ref) ./ abs.(uRs_ref)) < 1e-12

    # Test that linear distortions is the min solution
    left_uvs = camera_image(camera_reference, left_track);
    right_uvs = camera_image(camera_reference, right_track);

    uLs = [uv[1] for uv in left_uvs];
    vLs = [uv[2] for uv in left_uvs];
    uRs = [uv[1] for uv in right_uvs];
    vRs = [uv[2] for uv in right_uvs];

    # use same vs for reference track
    uLs0 = left_track_image_u(camera_reference, trackprop, vLs);
    uRs0 = right_track_image_u(camera_reference, trackprop, vRs);

    duLdβ = v_to_dudβ(camera_reference, -trackprop.track_gauge / (2));
    duLdα = v_to_dudα(camera_reference, -trackprop.track_gauge / (2));
    duLdβs = duLdβ.(vLs);
    duLdαs = duLdα.(vLs);

    duRdβ = v_to_dudβ(camera_reference, trackprop.track_gauge / (2));
    duRdα = v_to_dudα(camera_reference, trackprop.track_gauge / (2));
    duRdβs = duRdβ.(vRs);
    duRdαs = duRdα.(vRs);

    # The local minimum choice for β and α should be the values used for the ref track.
    βs = LinRange(0.0, β, 20);
    αs = LinRange(0.0, α, 20);

    uLs_arr = [duLdβs .* β + duLdαs .* α + uLs0 for β in βs];
    uRs_arr = [duRdβs .* β + duRdαs .* α + uRs0 for α in αs];

    m, i = findmin([norm(us - uLs) for us in uLs_arr])
    @test m / norm(uLs) < 1e-4
    @test norm(βs[i] - β) / abs(β) < 1e-5

    m, i = findmin([norm(us - uRs) for us in uRs_arr])
    @test m / norm(uRs) < 1e-4
    @test norm(αs[i] - α) / abs(α) < 1e-5
end

@testset "camera distortions" begin

    N = 1000;
    dYs = (rand(N) .- 0.5) .* 0.1;
    dZs = (rand(N) .- 0.5) .* 0.1;
    dXYZs = [[0,dYs[i],dZs[i]] for i in 1:N];
    dψθφs = [[(rand() - 0.5), (rand() - 0.5), 0.0] .* (0.5 * pi/180.0) for i in 1:N];

    # curvatures
    βs = (rand(N) .- 0.5) .* 1e-5;
    αs = (rand(N) .- 0.5) .* 1e-5;

    δs = [ [-dXYZs[i][2:3]; dψθφs[i][1:2]; αs[i]; βs[i]] for i in 1:N];

    choose_distortions = [:Y,:Z,:ψ,:θ,:α,:β]

    data = map(1:N) do i

        left_track  = [
            [x, -trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        right_track = [
            [x, trackprop.track_gauge/2.0 + βs[i] * x^2, αs[i] *x^2]
        for x = LinRange(4.0,40.0,40)];

        camera = VideoCamera(cameraposition_reference + dXYZs[i];
            focal_length = f
            , ψθφ = ψθφ_ref + dψθφs[i]
        );

        left_uvs = camera_image(camera, left_track);
        right_uvs = camera_image(camera, right_track);

        vLs = [uv[2] for uv in left_uvs];
        uLs = [uv[1] for uv in left_uvs];
        vRs = [uv[2] for uv in right_uvs];
        uRs = [uv[1] for uv in right_uvs];

        uLs_ref = left_track_image_u(camera_reference, trackprop, vLs);
        uRs_ref = right_track_image_u(camera_reference, trackprop, vRs);

        ΔuL = left_track_Δu(camera_reference, trackprop, vLs; choose_distortions=choose_distortions)
        ΔuR = right_track_Δu(camera_reference, trackprop, vRs; choose_distortions=choose_distortions)

        nL = maximum(abs.(uLs - (uLs_ref + ΔuL * δs[i]))) / mean(abs.(uLs))
        sL = sum(uLs - (uLs_ref + ΔuL * δs[i])) / sum(abs.(uLs))

        nR = maximum(abs.(uRs - (uRs_ref + ΔuR * δs[i]))) / mean(abs.(uRs))
        sR = sum(uRs - (uRs_ref + ΔuR * δs[i])) / sum(abs.(uRs))
        [nL,nR,sL,sR]
    end;

    nLs = [d[1] for d in data];
    nRs = [d[2] for d in data];
    sLs = [d[3] for d in data];
    sRs = [d[4] for d in data];

    @test maximum(nLs) < 0.02
    @test mean(nLs) < 0.002
    @test abs(mean(sLs)) < 0.001

    @test maximum(nRs) < 0.02
    @test mean(nRs) < 0.005
    @test abs(mean(sRs)) < 0.001

    # plot(uLs,vLs)
    # scatter!(uLs_ref,vLs)
    # plot!(uLs_ref + ΔuL * δ,vLs,linestyle = :dash, linewidth = 2)
    # #
    # plot(uRs,vRs)
    # scatter!(uRs_ref,vRs)
    # plot!(uRs_ref + ΔuR * δ,vRs,linestyle = :dash, linewidth = 2)

end


text:
function base_finder(seq)
    maximum( s -> parse(Int,s[end]), seq ) + 1
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas