text:
if VERSION >= v"0.7.0-beta.0"
    using Test
    using Dates
    using Printf
else
    using Base.Test
end
using NCDatasets

sz = (4,5)
filename = tempname()
#filename = "/tmp/test-6.nc"

# The mode "c" stands for creating a new file (clobber)
NCDatasets.Dataset(filename,"c") do ds

    # define the dimension "lon" and "lat"
    ds.dim["lon"] = sz[1]
    ds.dim["lat"] = sz[2]

    v = NCDatasets.defVar(ds,"small",Float64,("lon","lat"))
    @test_throws NCDatasets.NetCDFError v[:] = zeros(sz[1]+1,sz[2])
    #@test_throws NCDatasets.NetCDFError v[1:sz[1],1:sz[2]] = zeros(sz[1]+1,sz[2])
    @test_throws NCDatasets.NetCDFError v[sz[1]+1,1] = 1
    @test_throws NCDatasets.NetCDFError v[-1,1] = 1

    # variables
    for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64,
              Char,String]
    #for T in []
        local data
        data, scalar_data =
            if T == String
                [Char(i+60) * Char(j+60) for i = 1:sz[1], j = 1:sz[2]], "abcde"
            else
                [T(i+2*j) for i = 1:sz[1], j = 1:sz[2]], T(100)
            end

        v = NCDatasets.defVar(ds,"var-$T",T,("lon","lat"))
        v[:,:] = data
        @test v[:,:] == data[:,:]

        @test v[2,:] == data[2,:]

        @test v[:,3] == data[:,3]

        @test v[2,3] == data[2,3]

        # ignore extra index
        @test v[2,3,1,1] == data[2,3,1,1]

        # ignore extra index
        @test v[2:3,3,1,1] == data[2:3,3,1,1]

        @test v[[1,2,3],:] == data[[1,2,3],:]

        # write scalar,
        v.var[:,:] = scalar_data
        @test all(v.var[:,:][:] .== scalar_data)

    end
end

# quick interface
Dataset(filename,"c") do ds
    data = Int32[i+3*j for i = 1:sz[1], j = 1:sz[2]]
    defVar(ds,"temp",data,("lon","lat"), attrib = [
        "units" => "degree_Celsius",
        "long_name" => "Temperature"
    ])
    @test ds["temp"][:] == data
    @test eltype(ds["temp"].var) == Int32
    @test ds.dim["lon"] == sz[1]
    @test ds.dim["lat"] == sz[2]

    # load in-place
    data2 = similar(data)
    NCDatasets.load!(ds["temp"].var,data2,:,:)
    @test data2 == data

    data2 = zeros(eltype(data),sz[1],2)
    NCDatasets.load!(ds["temp"].var,data2,:,1:2)
    @test data2 == data[:,1:2]
end


text:
const madsen = let res_init=zeros(3), jac_init=zeros(3,2), x_init=[3.0, 1.0]

    function res(x, r)
        r[1] = x[1] ^ 2 + x[2] ^ 2 + x[1] * x[2]
        r[2] = sin(x[1])
        r[3] = cos(x[2])
        return r
    end

    function jac(x, j)
        j[1,1] = 2.e0 * x[1] + x[2]
        j[1,2] = 2.e0*x[2] + x[1]
        j[2,1] = cos(x[1])
        j[2,2] = 0.e0
        j[3,1] = 0.e0
        j[3,2] = -sin(x[2])
        return j
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("madsen", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end

ext: jl


text:
using Test, Flux, SparseArrays, Mill
using Mill: BagConv, convsum, bagconv, legacy_bagconv, _convshift, ∇convsum, ArrayNode, BagNode, ∇wbagconv, ∇xbagconv

@testset "testing convolution shift" begin
    @test _convshift(2) == 0:1
    @test _convshift(3) == -1:1
    @test _convshift(4) == -1:2
    @test _convshift(5) == -2:2
end

@testset "testing matvec and vecvec products " begin
    W = randn(3, 4)
    Wt = Matrix(transpose(W))
    xs = sprand(4, 10, 0.5)
    x = Matrix(xs)

    o = zeros(3, 10)
    foreach(i -> Mill._addmatvec!(o, i, W, x, i), 1:10)
    @test o ≈ W*x
    fill!(o, 0)
    foreach(i -> Mill._addmatvec!(o, i, W, xs, i), 1:10)
    @test o ≈ W*x

    o = zeros(3, 1)
    foreach(i -> Mill._addmatvec!(o, 1, W, x, i), 1:10)
    @test o ≈ sum(W*x, dims = 2)
    fill!(o, 0)
    foreach(i -> Mill._addmatvec!(o, 1, W, xs, i), 1:10)
    @test o ≈ sum(W*x, dims = 2)

    o = zeros(3, 10)
    foreach(i -> Mill._addmattvec!(o, i, Wt, x, i), 1:10)
    @test o ≈ W*x

    o = zeros(3, 1)
    foreach(i -> Mill._addmattvec!(o, 1, Wt, x, i), 1:10)
    @test o ≈ sum(W*x, dims = 2)

    xs = sprand(10, 1, 0.5)
    r, s = randn(10,1), Matrix(xs)
    o = zeros(10, 10)
    Mill._addvecvect!(o, r, 1, s, 1)
    @test o ≈ r * transpose(s)
    fill!(o, 0)
    Mill._addvecvect!(o, r, 1, xs, 1)
    @test o ≈ r * transpose(s)
end

@testset "testing forward convolution & gradient" begin
    x = Float64.([1 10  100  1000  10000]);
    y = 2 .* x;
    z = 4 .* x;
    Δ = ones(1, 5)
    for bags in [AlignedBags([1:2,3:5]), ScatteredBags([[1,2],[3,4,5]])];

        @test convsum(bags, x) == x
        @test convsum(bags, x, y) == [21  10  2100  21000  10000]
        @test convsum(bags, x, y, z) == [42  21  4200  42100  21000]

        @test all(ngradient(x -> sum(convsum(bags,x, y, z)), x)[1] .== ∇convsum(Δ, bags, 3)[1])
        @test all(ngradient(y -> sum(convsum(bags,x, y, z)), y)[1] .== ∇convsum(Δ, bags, 3)[2])
        @test all(ngradient(z -> sum(convsum(bags,x, y, z)), z)[1] .== ∇convsum(Δ, bags, 3)[3])

        @test gradtest((a, b, c) -> convsum(bags, a, b, c), x, y, z)
    end
end

@testset "testing the convolution" begin
    xs = sprand(3, 15, 0.5)
    x = Matrix(xs)
    filters = randn(4, 3, 3)
    fs = [filters[:,:,i] for i in 1:3]
    for bags in [AlignedBags([1:1, 2:3, 4:6, 7:15]), ScatteredBags(collect.([1:1, 2:3, 4:6, 7:15]))]
        @test bagconv(x, bags, fs...) ≈ legacy_bagconv(x, bags, filters)
        @test bagconv(x, bags, fs...) ≈ bagconv(xs, bags, fs...)
        @test isapprox(ngradient(f -> sum(bagconv(x, bags, f, fs[2], fs[3])), fs[1])[1],  ∇wbagconv(ones(4, 15), x, bags, fs...)[1], rtol = 1e-6)
        @test isapprox(ngradient(f -> sum(bagconv(xs, bags, f, fs[2], fs[3])), fs[1])[1],  ∇wbagconv(ones(4, 15), xs, bags, fs...)[1], rtol = 1e-6)
        @test isapprox(ngradient(x -> sum(bagconv(x, bags, fs...)), x)[1],  ∇xbagconv(ones(4, 15), x, bags, fs...), rtol = 1e-5)
        @test gradtest((a, b, c) -> bagconv(x, bags, a, b, c), fs...)
        @test gradtest((a, b, c) -> bagconv(xs, bags, a, b, c), fs...)
        @test gradtest(x -> bagconv(x, bags, fs...), x)
        @test gradtest((x, a, b , c) -> bagconv(x, bags, a, b, c), x, fs...)
    end
end

@testset "testing convolution with ScatteredBags" begin
    xs = sprand(3, 7, 0.5)
    x = Matrix(xs)
    filters = randn(4, 3, 3)
    fs = [filters[:,:,i] for i in 1:3]
    baga = AlignedBags([1:3, 4:7])
    bags = ScatteredBags([[1,2,3],[4,5,6,7]])
    bagsp = ScatteredBags([[1,4,7],[2,3,5,6]])
    xp = x[:, [1,4,5,2,6,7,3]]

    @testset "Test that all versions of convolutions are equal" begin
        @test bagconv(x, baga, fs...) == bagconv(x, bags, fs...)
        @test x[:,bags[1]] == xp[:,bagsp[1]]
        @test x[:,bags[2]] == xp[:,bagsp[2]]
        @test bagconv(x, bags, fs...)[:,bags[1]] == bagconv(xp, bagsp, fs...)[:,bagsp[1]]
        @test bagconv(x, bags, fs...)[:,bags[2]] == bagconv(xp, bagsp, fs...)[:,bagsp[2]]
    end

    @testset "Test that gradient of scattered convolution is correct" begin
        @test isapprox(ngradient(f -> sum(bagconv(xp, bagsp, f, fs[2], fs[3])), fs[1])[1],  ∇wbagconv(ones(4, 15), xp, bagsp, fs...)[1], rtol = 1e-6)
        @test isapprox(ngradient(f -> sum(bagconv(x, bags, f, fs[2], fs[3])), fs[1])[1],  ∇wbagconv(ones(4, 15), xp, bagsp, fs...)[1], rtol = 1e-6)
        @test isapprox(ngradient(x -> sum(bagconv(x, bagsp, fs...)), xp)[1],  ∇xbagconv(ones(4, 15), xp, bagsp, fs...), rtol = 1e-5)
        @test gradtest((a, b, c) -> bagconv(x, bags, a, b, c), fs...)
        @test gradtest(x -> bagconv(x, bagsp, fs...), xp)
        @test gradtest((x, a, b , c) -> bagconv(x, bagsp, a, b, c), xp, fs...)
    end
end

@testset "testing convolution layer" begin
    xs = sprand(3, 15, 0.5)
    x = Matrix(xs)
    for bags in [AlignedBags([1:1, 2:3, 4:6, 7:15]), ScatteredBags(collect.([1:1, 2:3, 4:6, 7:15]))]
        ds = BagNode(ArrayNode(Float32.(x)), bags)

        m = BagConv(3, 4, 3, relu)
        @test length(params(m)) == 3
        @test size(m(x, bags)) == (4, 15)
        @test size(m(xs, bags)) == (4, 15)

        m = BagConv(3, 4, 1)
        @test size(m(x, bags)) == (4, 15)
        @test size(m(xs, bags)) == (4, 15)

        m = BagChain(BagConv(3, 4, 3, relu), BagConv(3, 4, 2))
        @test length(params(m)) == 5
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

module FMU2

include(joinpath(dirname(@__FILE__), "FMU2_content.jl"))

end # module


text:
"""
    function polarization(para, extK = DiagTree.getK(para.totalLoopNum, 1), subdiagram = false; name = :Π, resetuid = false)

    Generate polarization diagrams using Parquet Algorithm.

# Arguments
- `para`            : parameters. It should provide internalLoopNum, interactionTauNum, firstTauIdx
- `extK`            : basis of external loop.
- `subdiagram`      : a sub-vertex or not
- `name`            : name of the vertex
- `resetuid`        : restart uid count from 1

# Output
- A DataFrame with fields `:response`, `:diagram`, `:hash`.
- All polarization share the same external Tau index. With imaginary-time variables, they are extT = (para.firstTauIdx, para.firstTauIdx+1)
"""
function polarization(para, extK=DiagTree.getK(para.totalLoopNum, 1), subdiagram=false; name=:Π, resetuid=false)
    resetuid && uidreset()
    @assert para.diagType == PolarDiag
    @assert para.innerLoopNum >= 1
    # @assert length(extK) == para.totalLoopNum
    @assert length(extK) >= para.totalLoopNum "expect dim of extK>=$(para.totalLoopNum), got $(length(extK))"
    extK = extK[1:para.totalLoopNum]

    #polarization diagram is always proper
    para = reconstruct(para, filter=union(Proper, para.filter), transferLoop=extK)

    if (para.extra isa ParquetBlocks) == false
        para = reconstruct(para, extra=ParquetBlocks())
    end

    K = zero(extK)
    LoopIdx = para.firstLoopIdx
    K[LoopIdx] = 1.0
    @assert (K ≈ extK) == false
    t0 = para.firstTauIdx
    extT = para.hasTau ? (t0, t0 + 1) : (t0, t0)
    legK = [extK, K, K .- extK]

    polar = DataFrame()

    for (oVer3, oGin, oGout) in orderedPartition(para.innerLoopNum - 1, 3, 0)
        # ! Vertex3 must be in the first place, because we want to make sure that the bosonic extT of the vertex3 start with t0+1

        idx, maxLoop = findFirstLoopIdx([oVer3, oGin, oGout], LoopIdx + 1) # GGΓ3 consumes one internal loop
        @assert maxLoop <= para.totalLoopNum "maxLoop = $maxLoop > $(para.totalLoopNum)"
        Ver3Kidx, GinKidx, GoutKidx = idx

        if isValidG(para.filter, oGin) && isValidG(para.filter, oGout)
            if oVer3 == 0
                ######################## Π0 = GG #########################################
                gt0 = para.hasTau ? extT[2] + 1 : extT[1]
                idx, maxTau = findFirstTauIdx([oGin, oGout], [GreenDiag, GreenDiag], gt0, para.interactionTauNum)
                @assert maxTau <= para.totalTauNum "maxTau = $maxTau > $(para.totalTauNum)"
                GinTidx, GoutTidx = idx

                paraGin = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGin,
                    firstLoopIdx=GinKidx, firstTauIdx=GinTidx)
                paraGout = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGout,
                    firstLoopIdx=GoutKidx, firstTauIdx=GoutTidx)

                response = UpUp
                polarid = PolarId(para, response, k=extK, t=extT)
                gin = green(paraGin, K, (extT[1], extT[2]), true, name=:Gin)
                gout = green(paraGout, K .- extK, (extT[2], extT[1]), true, name=:Gout)
                @assert gin isa Diagram && gout isa Diagram "$gin or $gout is not a single diagram"

                sign = para.isFermi ? -1.0 : 1.0
                polardiag = Diagram(polarid, Prod(), [gin, gout], name=name, factor=sign)
                push!(polar, (response=response, extT=extT, diagram=polardiag))
            else
                ##################### composite polarization #####################################
                idx, maxTau = findFirstTauIdx([oVer3, oGin, oGout], [Ver3Diag, GreenDiag, GreenDiag], extT[2], para.interactionTauNum)
                @assert maxTau <= para.totalTauNum "maxTau = $maxTau > $(para.totalTauNum)"
                Ver3Tidx, GinTidx, GoutTidx = idx

                paraGin = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGin,
                    firstLoopIdx=GinKidx, firstTauIdx=GinTidx)
                paraGout = reconstruct(para, diagType=GreenDiag, innerLoopNum=oGout,
                    firstLoopIdx=GoutKidx, firstTauIdx=GoutTidx)

                paraVer3 = reconstruct(para, diagType=Ver3Diag, innerLoopNum=oVer3,
                    firstLoopIdx=Ver3Kidx, firstTauIdx=Ver3Tidx)
                ver3 = vertex3(paraVer3, legK, true)
                if isnothing(ver3) || isempty(ver3)
                    continue
                end
                if para.hasTau
                    @assert all(x -> x[1] == extT[2], ver3[:, :extT]) "The bosonic T must be firstTauIdx+1 if hasTau\n$ver3"
                    @assert all(x -> x[2] == ver3[1, :extT][2], ver3[:, :extT]) "The TinL must be firstTauIdx+2 if hasTau\n$ver3"
                end

                #transform extT coloum into extT for Vertex4 and the extT for Gin and Gout
                df = transform(ver3, :extT => ByRow(x -> [extT, (extT[1], x[2]), (x[3], extT[1])]) => [:extT, :GinT, :GoutT])

                groups = mergeby(df, [:response, :GinT, :GoutT, :extT], operator=Sum())

                for v3 in eachrow(groups)
                    response = v3[:response]
                    @assert response == UpUp || response == UpDown
                    #type: Instant or Dynamic
                    polarid = PolarId(para, response, k=extK, t=v3[:extT])
                    gin = green(paraGin, K, v3[:GinT], true, name=:Gin)
                    gout = green(paraGout, K .- extK, v3[:GoutT], true, name=:Gout)
                    @assert gin isa Diagram && gout isa Diagram

                    polardiag = Diagram(polarid, Prod(), [gin, gout, v3[:diagram]], name=name)
                    push!(polar, (response=response, extT=v3[:extT], diagram=polardiag))
                end
            end
        end
    end

    # for (oGinL, oGoutL, oGinR, oGoutR, ver4) in orderedPartition(para.innerLoopNum - 1, 5, 0)
    # end
    if isempty(polar)
        return DataFrame(response=[], extT=[], diagram=[])
    end

    # legK = [extK, K, K, extK]
    Factor = 1 / (2π)^para.loopDim
    polar = mergeby(polar, [:response, :extT]; name=name, factor=Factor,
        getid=g -> PolarId(para, g[1, :response], k=extK, t=extT)
    )
    return polar
end


text:
#All time related types and functions

import Base: convert, isless, sleep, +, -, *, ==
export Time, Duration, Rate, to_sec, to_nsec, get_rostime, rossleep

#Time type definitions
abstract type AbstractTime end

"""
    Time(secs, nsecs), Time(), Time(t::Real)

Object representing an absolute time from a fixed past reference point at nanosecond precision.

Basic arithmetic can be performed on combinations of `Time` and `Duration` objects that make sense.
For example, if `t::Time` and `d::Duration`, `t+d` will be a `Time`, `d+d` a `Duration`, `t-d` a
`Time`, `d-d` a `Duration`, and `t-t` a `Duration`.
"""
struct Time <: AbstractTime
    secs::Int32
    nsecs::Int32
    function Time(s::Real,n::Real)
        cs, cns = _canonical_time(s,n)
        new(cs, cns)
    end
end
Time() = Time(0,0)
Time(t::Real) = Time(t,0)

"""
    Duration(secs, nsecs), Duration(), Duration(t::Real)

Object representing a relative period of time at nanosecond precision.

Basic arithmetic can be performed on combinations of `Time` and `Duration` objects that make sense.
For example, if `t::Time` and `d::Duration`, `t+d` will be a `Time`, `d+d` a `Duration`, `t-d` a
`Time`, `d-d` a `Duration`, and `t-t` a `Duration`.
"""
struct Duration <: AbstractTime
    secs::Int32
    nsecs::Int32
    function Duration(s::Real,n::Real)
        cs, cns = _canonical_time(s,n)
        new(cs, cns)
    end
end
Duration() = Duration(0,0)
Duration(t::Real) = Duration(t,0)

#Enforce 0 <= nsecs < 1e9
function _canonical_time(secs, nsecs)
    nsec_conv = convert(Int32, 1_000_000_000)
    secs32  = floor(Int32, secs)
    nsecs32 = floor(Int32, mod(secs,1)*1e9 + nsecs)

    addsecs = div(nsecs32, nsec_conv)
    crnsecs = rem(nsecs32, nsec_conv)
    if crnsecs < 0
        addsecs -= one(Int32)
        crnsecs += nsec_conv
    end
    (secs32 + addsecs, crnsecs)
end

+(t1::Time,     t2::Duration) = Time(    t1.secs+t2.secs, t1.nsecs+t2.nsecs)
+(t1::Duration, t2::Time)     = Time(    t1.secs+t2.secs, t1.nsecs+t2.nsecs)
+(t1::Duration, t2::Duration) = Duration(t1.secs+t2.secs, t1.nsecs+t2.nsecs)
-(t1::Time,     t2::Duration) = Time(    t1.secs-t2.secs, t1.nsecs-t2.nsecs)
-(t1::Duration, t2::Duration) = Duration(t1.secs-t2.secs, t1.nsecs-t2.nsecs)
-(t1::Time,     t2::Time)     = Duration(t1.secs-t2.secs, t1.nsecs-t2.nsecs)
*(td::Duration, tf::Real)     = Duration(tf*td.secs     , tf*td.nsecs)
*(tf::Real, td::Duration)     = Duration(tf*td.secs     , tf*td.nsecs)

#PyObject conversions
convert(::Type{Time},     o::PyObject) = Time(    o[:secs],o[:nsecs])
convert(::Type{Duration}, o::PyObject) = Duration(o[:secs],o[:nsecs])
convert(::Type{PyObject}, t::Time)     = __rospy__[:Time](    t.secs,t.nsecs)
convert(::Type{PyObject}, t::Duration) = __rospy__[:Duration](t.secs,t.nsecs)

#Real number conversions
"""
    to_sec(t)

Return the value of a ROS time object in absolute seconds (with nanosecond precision)
"""
to_sec(t::T) where {T <: AbstractTime} = t.secs + 1e-9*t.nsecs

"""
    to_nsec(t)

Return the value of a ROS time object in nanoseconds as an integer.
"""
to_nsec(t::T) where {T <: AbstractTime} = 1_000_000_000*t.secs + t.nsecs
convert(::Type{Float64}, t::T) where {T <: AbstractTime} = to_sec(t)

#Comparisons
==(t1::T, t2::T) where {T <: AbstractTime} = (t1.secs == t2.secs) && (t1.nsecs == t2.nsecs)
isless(t1::T, t2::T) where {T <: AbstractTime} = to_nsec(t1) < to_nsec(t2)

"""
    Rate(hz::Real), Rate(d::Duration)

Used to allow a loop to run at a fixed rate. Construct with a frequency or `Duration` and use with
`rossleep` or `sleep`. The rate object will record execution time of other work in the loop and
modify the sleep time to compensate, keeping the loop rate as consistent as possible.
"""
mutable struct Rate
    duration::Duration
    last_time::Time
end
Rate(d::Duration) = Rate(d, get_rostime())
Rate(hz::Real) = Rate(Duration(1.0/hz), get_rostime())

"""
    get_rostime()

Return the current ROS time as a `Time` object.
"""
function get_rostime()
    t = try
        __rospy__[:get_rostime]()
    catch ex
        error(pycall(pybuiltin("str"), PyAny, ex.val))
    end
    convert(Time, t)
end

"""
    RobotOS.now()

Return the current ROS time as a `Time` object.
"""
now() = get_rostime()

"""
    rossleep(t)

Sleep and process callbacks for a number of seconds implied by the type and value of `t`, which may
be a real-value, a `Duration` object, or a `Rate` object.
"""
function rossleep(td::Duration)
    #Busy sleep loop needed to allow both julia and python async activity
    tnsecs = to_nsec(td)
    t0 = time_ns()
    while time_ns()-t0 < tnsecs
        yield()                   #Allow julia callback loops to run
        __rospy__[:sleep](0.001)  #Allow rospy comm threads to run
    end
end
rossleep(t::Real) = rossleep(Duration(t))

function rossleep(r::Rate)
    ctime = get_rostime()
    if r.last_time > ctime
        r.last_time = ctime
    end
    elapsed = ctime - r.last_time
    rossleep(r.duration - elapsed)
    r.last_time += r.duration

    if ctime - r.last_time > r.duration*2
        r.last_time = ctime
    end
end

"""
    sleep(t::Duration), sleep(t::Rate)

Call `rossleep` with a `Duration` or `Rate` object. Use `rossleep` to specify sleep time directly.
"""
sleep(t::Duration) = rossleep(t)
sleep(t::Rate) = rossleep(t)


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module RiordanSquares

using Nemo, SeriesExpansion, Triangles

export ModuleRiordanSquares
export RiordanProduct, RiordanSquare
export T039599, T116392, T172094, T321620, T321621
export T321623, T321624, T322942

"""

The Riordan product is a map a, b ↦ [a, b] associating two formal power series a, b with a lower triangular matrix [a, b]. The Riordan square is the case a = b of the Riordan product. Formally we can describe the Riordan square as a transform RS: Z[[x]] ↦ Mat[Z] which maps power series over the integers to (lower triangular) integer matrices.

* RiordanProduct, RiordanSquare
* T039599, T116392, T172094, T321620, T321621, T321623, T321624, T322942

[Introduction to the Riordan Square](http://luschny.de/math/seq/RiordanSquare.html)
"""
const ModuleRiordanSquares = ""

"""

Return the Riordan array associated with the generating functions a and b.
"""
function RiordanProduct(a, b, dim, expo = false)
    A = Coefficients(a, dim)
    B = b === nothing ? A : Coefficients(b, dim)
    M = identity_matrix(QQ, dim)
    for k ∈ 1:dim
        M[k, 1] = A[k]
    end

    for k ∈ 2:dim, m ∈ k+1:dim
        M[m, k] = sum(M[j+1, k-1] * B[m-j] for j ∈ k-2:m-2)
    end
    #expo ? ExponentialWeights(M) : M
    toΔ(M)
end

"""

Return the Riordan array (Riordan product) ``a \\times a``.
"""
RiordanSquare(a, n, expo = false) = RiordanProduct(a, nothing, n, expo)

"""

The Riordan square of the Catalan numbers.
"""
T039599(n) = RiordanSquare(G000108, n)

"""

The Riordan square of the central trinomial.
"""
T116392(n) = RiordanSquare(G002426, n)

"""

The Riordan square of the little Schröder numbers.
"""
T172094(n) = RiordanSquare(G001003, n)

"""

The Riordan square of the Riordan numbers (with 1 prepended).
"""
T321620(n) = RiordanSquare(G005043, n)

"""

The Riordan square of the Motzkin numbers.
"""
T321621(n) = RiordanSquare(G001006, n)

"""

The Riordan square of the large Schröder numbers.
"""
T321623(n) = RiordanSquare(G006318, n)

"""

The Riordan square of the Lucas numbers.
"""
T321624(n) = RiordanSquare(G000032, n)

"""

The Riordan square of the Jacobsthal numbers.
"""
T322942(n) = RiordanSquare(G001045, n)

#"""
#The Riordan square of the number of rooted bicubic maps.
#(Does not yet exist in the OEIS.)
#"""
# Tx(n) = RiordanSquare(G000257, n)

#START-TEST-########################################################

using Test

function test()
    @testset "RiordanSqr" begin
        a = T039599(5)
        b = [1, 1, 1, 2, 3, 1, 5, 9, 5, 1, 14, 28, 20, 7, 1]
        @test all(a .== b)

        a = T116392(5)
        b = [1, 1, 1, 3, 4, 1, 7, 13, 7, 1, 19, 42, 32, 10, 1]
        @test all(a .== b)

        a = T172094(5)
        b = [1, 1, 1, 3, 4, 1, 11, 17, 7, 1, 45, 76, 40, 10, 1]
        @test all(a .== b)

        a = T321620(5)
        b = [1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1]
        @test all(a .== b)

        a = T321621(5)
        b = [1, 1, 1, 2, 3, 1, 4, 8, 5, 1, 9, 21, 18, 7, 1]
        @test all(a .== b)

        ###############################################
        a = T321623(5)
        b = T321623(5)
        @test all(a .== b)

        a = T321624(5)
        b = [1, 1, 1, 3, 4, 1, 4, 10, 7, 1, 7, 24, 26, 10, 1]
        @test all(a .== b)

        a = T322942(5)
        b = [1, 1, 1, 1, 2, 1, 3, 5, 3, 1, 5, 12, 10, 4, 1]
        @test all(a .== b)
    end
end

function demo()
    println("\nT039599")
    T039599(5) |> ShowAsΔ

    println("\nT116392")
    T116392(5) |> ShowAsΔ

    println("\nT172094")
    T172094(5) |> ShowAsΔ

    println("\nT321620")
    T321620(5) |> ShowAsΔ

    println("\nT321621")
    T321621(5) |> ShowAsΔ

    println("\nT321623")
    T321623(5) |> ShowAsΔ

    println("\nT321624")
    T321624(5) |> ShowAsΔ

    println("\nT322942")
    T322942(5) |> ShowAsΔ
end

function perf() end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
macro diff!(ex)
    f_name = ex.args[1]
    df_name = Symbol(String(f_name) * "′")
    quote
        function $(esc(df_name))(x)
            x_var = Variable(x)
            y = $(esc(ex.args[1]))(x_var)
            backward!(y)
            return x_var.grad.values
        end
    end
end


text:
####################################################################
function gpmodelVI(x, flux, σ2, idx; logs=logs, maxiter = 10, seed = 1, show_trace = true)
####################################################################

    #--------------------------------------------------------
    # Auxiliary
    #--------------------------------------------------------

    rg = MersenneTwister(seed)

    JITTER = 1e-8

    Σdiag = Diagonal(σ2)

    numF  = length(unique(idx))

    @assert(all(sort(unique(idx)) .== collect(1:numF)))

    N = length(x)

    @assert(N == length(flux) == length(idx) == length(σ2))

    @printf("Running gpmodelVIexact_grad with %d and %d filters\n", N, numF)

    #--------------------------------------------------------
    # Calculate kernel matrix and keep fixed!
    #--------------------------------------------------------

    dist2 = pairwise(SqEuclidean(), reshape(x, 1, N), dims=2)

    K1 = zeros(N, N)

    calculatekernelmatrix!(K1, dist2, rbf, [0.0; logs])

    K1 = K1 + JITTER * I

    K = PDMat((K1 + K1') / 2) # This helps us speed calculations

    #--------------------------------------------------------------------
    # Define functions for calling lower bound.
    # Actual calculation takes place in function "complete_lower_bound"
    #--------------------------------------------------------------------

    v²ₐ = v²ᵦ = 1000.0

    # Convenient calls to lower bound

    function lb(μα, Σα, μβ, Σβ, μf, Σf)

        complete_lower_bound(; y=flux, idx=idx, S=Σdiag, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    end

    function lb_grad(μα, Σα, μβ, Σβ, μf, Σf)

        complete_lower_bound_grad(; y=flux, idx=idx, S=Σdiag, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    end

    # Unpacks vectorised parameters

    function unpack(param)

        @assert(length(param) == numF*4)

        local μα =               param[1+0*numF:1*numF]

        # the parametrisation used here has implications on the calculation
        # of the gradient, see complete_lower_bound.jl and exact_expectations_grad.jl
        local Σα = Diagonal(exp.(param[1+1*numF:2*numF]))

        local μβ =               param[1+2*numF:3*numF]

        # same comment as above applies here
        local Σβ = Diagonal(exp.(param[1+3*numF:4*numF]))

        return μα, Σα, μβ, Σβ

    end

    # Constructs necessary matrices

    function construct_matrices(μα, Σα, μβ, Σβ)

        local Ā  = Diagonal(μα[idx])

        local Σa = Diagonal(Σα.diag[idx])

        local b̄  = μβ[idx]

        return Ā, Σa, b̄

    end

    # Calculate optimal posterior q(f)

    function optimal_qf(Ā, Σa, b̄)

        # local Σf = ((Ā*inv(Σdiag)*Ā + inv(Σdiag)*Σa) + inv(K))\I

        local invDiagPart = (Ā*(Σdiag\Ā) + Σdiag\Σa)\I #
                                                       # This is faster/more stable
        local Σf = invDiagPart*((invDiagPart + K)\K)   #

        local μf = Σf * Ā * (Σdiag \ (flux - b̄))

        return μf, Σf

    end

    function fg!(F,G, param)

      # do common computations here

      local μα, Σα, μβ, Σβ = unpack(param)

      local Ā, Σa, b̄ = construct_matrices(μα, Σα, μβ, Σβ)

      local μf, Σf = optimal_qf(Ā, Σa, b̄)

      if G != nothing
        # code to compute gradient here
        # writing the result to the vector G
        local gradα, gradΣα, gradβ, gradΣβ = lb_grad(μα, Σα, μβ, Σβ, μf, Σf)

        copyto!(G, -1.0 * [gradα; gradΣα; gradβ; gradΣβ])

      end

      if F != nothing
        # value = ... code to compute objective function
        local value = -1.0 * lb(μα, Σα, μβ, Σβ, μf, Σf)
        return value
      end

    end

    #--------------------------------------------------------
    # Initialise GP hyperparameters with random search
    # and few iterations
    #--------------------------------------------------------

    a0 = [std(flux[findall(i.==idx)])  for i in unique(idx)]

    b0 = [mean(flux[findall(i.==idx)]) for i in unique(idx)]

    initsolutions   = [[a0.+3*randn(numF); log.(3*rand()*ones(numF)); b0.+3*randn(numF); log.(3*rand()*ones(numF))] for i=1:100]

    initialfitness  = @showprogress "\t init random search " map(x->fg!(0.0, nothing, x), initsolutions)

    solution = initsolutions[argmin(initialfitness)]

    # Use the LBFGS optimiser

    opt = Optim.Options(show_trace = show_trace, iterations = maxiter, show_every = 50, f_tol=0.0, g_tol=1e-8)

    result = Optim.optimize(Optim.only_fg!(fg!), solution, LBFGS(), opt)

    # retrieve result and return posterior of parameters describing
    # fitted line as object of type Distributions.MvNormal

    μα, Σα, μβ, Σβ = unpack(result.minimizer)

    return MvNormal([μα; μβ], Diagonal([Σα.diag; Σβ.diag]))

end


text:
# ------------------------------------------------------------------------------
#
#  Gmsh Python tutorial 16
#
#  Constructive Solid Geometry, OpenCASCADE geometry kernel
#
# ------------------------------------------------------------------------------
# Instead of constructing a model in a bottom-up fashion with Gmsh's built-in
# geometry kernel, starting with version 3 Gmsh allows you to directly use
# alternative geometry kernels. Here we will use the OpenCASCADE kernel.
import gmsh

gmsh.initialize()

gmsh.model.add("t16")

# Let's build the same model as in `t5.py', but using constructive solid
# geometry.
# We first create two cubes:
gmsh.model.occ.addBox(0,0,0, 1,1,1, 1)
gmsh.model.occ.addBox(0,0,0, 0.5,0.5,0.5, 2)
# We apply a boolean difference to create the "cube minus one eigth" shape:
gmsh.model.occ.cut([(3,1)], [(3,2)], 3)

# Boolean operations with OpenCASCADE always create new entities. By default the
# extra arguments `removeObject' and `removeTool' in `cut()' are set to `True',
# which will delete the original entities.

# We then create the five spheres:
x = 0; y = 0.75; z = 0; r = 0.09

holes = []
for t in 1:5
    global x, z
    x += 0.166
    z += 0.166
    gmsh.model.occ.addSphere(x,y,z,r, 3 + t)
    t = (3, 3 + t)
    push!(holes, t)
end

# If we had wanted five empty holes we would have used `cut()' again. Here we
# want five spherical inclusions, whose mesh should be conformal with the mesh
# of the cube: we thus use `fragment()', which intersects all volumes in a
# conformal manner (without creating duplicate interfaces):
ov = gmsh.model.occ.fragment([(3,3)], holes)
gmsh.model.occ.synchronize()

# Creating entities using constructive solid geometry is very powerful, but can
# lead to practical issues for e.g. setting mesh sizes at points, or identifying
# boundaries.

# To identify points or other bounding entities you can take advantage of the
# `getEntities()', `getBoundary()' and `getEntitiesInBoundingBox()' functions:
lcar1 = .1
lcar2 = .0005
lcar3 = .055

# Assign a mesh size to all the points:
ov = gmsh.model.getEntities(0);
gmsh.model.mesh.setSize(ov, lcar1);

# Override this constraint on the points of the five spheres:
ov = gmsh.model.getBoundary(holes, false, false, true);
gmsh.model.mesh.setSize(ov, lcar3);

# Select the corner point by searching for it geometrically:
eps = 1e-3
ov = gmsh.model.getEntitiesInBoundingBox(0.5-eps, 0.5-eps, 0.5-eps,
                                         0.5+eps, 0.5+eps, 0.5+eps, 0)
gmsh.model.mesh.setSize(ov, lcar2)

gmsh.model.mesh.generate(3)

gmsh.write("t16.msh")

gmsh.finalize()


text:
# import GPUArrays
import Knet
import ProgressMeter
import ValueHistories

function KnetModel(
        ;
        name::AbstractString = "",
        predict_function_source::AbstractString = "",
        loss_function_source::AbstractString = "",
        predict_function::Function = identity,
        loss_function::Function = identity,
        losshyperparameters::AbstractDict = Dict(),
        optimizationalgorithm::Symbol = :nothing,
        optimizerhyperparameters::AbstractDict = Dict(),
        minibatchsize::Integer = 0,
        modelweights::AbstractArray = [],
        isclassificationmodel::Bool = false,
        isregressionmodel::Bool = false,
        maxepochs::Integer = 0,
        printlosseverynepochs::Integer = 0,
        )
    optimizersymbol2type=Dict()
    optimizersymbol2type[:Sgd] = Knet.Sgd
    optimizersymbol2type[:Momentum] = Knet.Momentum
    optimizersymbol2type[:Nesterov] = Knet.Nesterov
    optimizersymbol2type[:Rmsprop] = Knet.Rmsprop
    optimizersymbol2type[:Adagrad] = Knet.Adagrad
    optimizersymbol2type[:Adadelta] = Knet.Adadelta
    optimizersymbol2type[:Adam] = Knet.Adam
    optimizersymbol2type=fix_type(optimizersymbol2type)
    modelweightoptimizers = Knet.optimizers(
        modelweights,
        optimizersymbol2type[optimizationalgorithm];
        optimizerhyperparameters...
        )
    # lastepoch = 0
    # lastiteration = 0
    history = ValueHistories.MVHistory()
    ValueHistories.push!(
        history,
        :epoch_at_iteration,
        0,
        0,
        )
    losshyperparameters = fix_type(losshyperparameters)
    optimizerhyperparameters = fix_type(optimizerhyperparameters)
    result = KnetModel(
        name,
        isclassificationmodel,
        isregressionmodel,
        predict_function_source,
        loss_function_source,
        predict_function,
        loss_function,
        losshyperparameters,
        optimizationalgorithm,
        optimizerhyperparameters,
        minibatchsize,
        maxepochs,
        printlosseverynepochs,
        modelweights,
        modelweightoptimizers,
        history,
        )
    return result
end

"""
"""
function set_max_epochs!(
        x::KnetModel,
        new_max_epochs::Integer,
        )
    x.maxepochs = new_max_epochs
    return nothing
end

"""
"""
function get_underlying(
        x::KnetModel;
        saving::Bool = false,
        loading::Bool = false,
        )
    result = (x.modelweights, x.modelweightoptimizers,)
    return result
end

"""
"""
function get_history(
        x::KnetModel;
        saving::Bool = false,
        loading::Bool = false,
        )
    result = x.history
    return result
end

function parse_functions!(estimator::KnetModel)
    estimator.predict_function = eval(
        Meta.parse(
            strip(
                estimator.predict_function_source
                )
            )
        )
    estimator.loss_function = eval(
        Meta.parse(
            strip(
                estimator.loss_function_source
                )
            )
        )
    return nothing
end

"""
"""
function fit!(
        estimator::KnetModel,
        training_features_array::AbstractArray,
        training_labels_array::AbstractArray,
        tuning_features_array::Union{Nothing, AbstractArray} = nothing,
        tuning_labels_array::Union{Nothing, AbstractArray} = nothing,
        )
    has_tuning_data::Bool =
        !is_nothing(tuning_labels_array) && !is_nothing(tuning_features_array)
    training_features_array = Float64.(training_features_array)
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        training_labels_array = Int.(training_labels_array)
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        training_labels_array = Float64.(training_labels_array)
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
    training_data = Knet.minibatch(
        training_features_array,
        training_labels_array,
        estimator.minibatchsize,
        )
    loss_function_gradient = Knet.grad(
        estimator.loss_function,
        2,
        )
    all_iterations_so_far, all_epochs_so_far = ValueHistories.get(
        estimator.history,
        :epoch_at_iteration,
        )
    last_iteration = all_iterations_so_far[end]
    last_epoch = all_epochs_so_far[end]
    @info(
        string(
            "Starting to train Knet model. Max epochs: ",
            estimator.maxepochs,
            ".",
            )
        )
    training_lossbeforetrainingstarts = estimator.loss_function(
        estimator.predict_function,
        estimator.modelweights,
        training_features_array,
        training_labels_array;
        estimator.losshyperparameters...
        )
    if has_tuning_data
        tuning_lossbeforetrainingstarts = estimator.loss_function(
           estimator.predict_function,
           estimator.modelweights,
           tuning_features_array,
           tuning_labels_array;
           estimator.losshyperparameters...
           )
    end
    if (estimator.printlosseverynepochs) > 0
        if has_tuning_data
            @info(
                string(
                    "Epoch: ",
                    last_epoch,
                    ". Loss (training set): ",
                    training_lossbeforetrainingstarts,
                    ". Loss (tuning set): ",
                    tuning_lossbeforetrainingstarts,
                    ".",
                    )
                )
        else
            @info(
                string(
                    "Epoch: ",
                    lastepoch,
                    ". Loss: ",
                    lossbeforetrainingstarts,
                    "."
                    )
                )
        end
    end
    while last_epoch < estimator.maxepochs
        for (x_training, y_training) in training_data
            grads = loss_function_gradient(
                estimator.predict_function,
                estimator.modelweights,
                x_training,
                y_training;
                estimator.losshyperparameters...
                )
            Knet.update!(
                estimator.modelweights,
                grads,
                estimator.modelweightoptimizers,
                )
            last_iteration += 1
            training_currentiterationloss = estimator.loss_function(
                estimator.predict_function,
                estimator.modelweights,
                x_training,
                y_training;
                estimator.losshyperparameters...
                )
            ValueHistories.push!(
                estimator.history,
                :training_loss_at_iteration,
                last_iteration,
                training_currentiterationloss,
                )
        end # end for
        last_epoch += 1
        ValueHistories.push!(
            estimator.history,
            :epoch_at_iteration,
            last_iteration,
            last_epoch,
            )
        training_currentepochloss = estimator.loss_function(
            estimator.predict_function,
            estimator.modelweights,
            training_features_array,
            training_labels_array;
            estimator.losshyperparameters...
            )
        ValueHistories.push!(
            estimator.history,
            :training_loss_at_epoch,
            last_epoch,
            training_currentepochloss,
            )
        if has_tuning_data
            tuning_currentepochloss = estimator.loss_function(
                estimator.predict_function,
                estimator.modelweights,
                tuning_features_array,
                tuning_labels_array;
                estimator.losshyperparameters...
                )
            ValueHistories.push!(
                estimator.history,
                :tuning_loss_at_epoch,
                last_epoch,
                tuning_currentepochloss,
                )
        end
        printlossthisepoch = (estimator.printlosseverynepochs > 0) &&
            ( (last_epoch == estimator.maxepochs) ||
                ( (last_epoch %
                    estimator.printlosseverynepochs) == 0 ) )
        if printlossthisepoch
            if has_tuning_data
                @info(
                   string(
                       "Epoch: ",
                       last_epoch,
                       ". Loss (training set): ",
                       training_currentepochloss,
                       ". Loss (tuning set): ",
                       tuning_currentepochloss,
                       ".",
                       ),
                   )
            else
                @info(
                   string(
                       "Epoch: ",
                       last_epoch,
                       ". Loss: ",
                       training_currentepochloss,
                       ".",
                       ),
                   )
            end
        end
    end # end while
    @info(string("Finished training Knet model."))
    return estimator
end

"""
"""
function predict(
        estimator::KnetModel,
        featuresarray::AbstractArray,
        )
    if estimator.isclassificationmodel
        probabilitiesassoc = predict_proba(
            estimator,
            featuresarray,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc
            )
        return predictionsvector
    elseif estimator.isregressionmodel
        output = estimator.predict_function(
            estimator.modelweights,
            featuresarray;
            )
        outputtransposed = transpose(output)
        result = convert(Array, outputtransposed)
        return result
    else
        error("unable to predict")
    end
end

"""
"""
function predict(
        estimator::KnetModel,
        featuresarray::AbstractArray,
        positive_class::Integer,
        threshold::AbstractFloat,
        )
    if estimator.isclassificationmodel
        probabilitiesassoc = predict_proba(
            estimator,
            featuresarray,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc,
            positive_class,
            threshold,
            )
        return predictionsvector
    else
        error("can only use the `threshold` argument with classification models")
    end
end

"""
"""
function predict_proba(
        estimator::KnetModel,
        featuresarray::AbstractArray,
        )
    if estimator.isclassificationmodel
        output = estimator.predict_function(
            estimator.modelweights,
            featuresarray;
            probabilities = true,
            )
        outputtransposed = transpose(output)
        numclasses = size(outputtransposed, 2)
        result = Dict()
        for i = 1:numclasses
            result[i] = outputtransposed[:, i]
        end
        return result
    elseif estimator.isregressionmodel
        error("predict_proba is not defined for regression models")
    else
        error("unable to predict")
    end
end

"""
"""
function single_labelmulticlassdataframeknetclassifier_Knet(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        name::AbstractString = "",
        predict_function_source::AbstractString = "",
        loss_function_source::AbstractString = "",
        losshyperparameters::AbstractDict = Dict(),
        optimizationalgorithm::Symbol = :nothing,
        optimizerhyperparameters::AbstractDict = Dict(),
        minibatchsize::Integer = 0,
        modelweights::AbstractArray = [],
        maxepochs::Integer = 0,
        printlosseverynepochs::Integer = 0,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} =
            nothing,
        )
    label_names = [single_label_name]
    label_levels = Dict()
    label_levels[single_label_name] = single_label_levels
    label_levels = fix_type(label_levels)
    dftransformer_index = 1
    dftransformer_transposefeatures = true
    dftransformer_transposelabels = true
    dftransformer = MutableDataFrame2ClassificationKnetTransformer(
        feature_names,
        label_names,
        label_levels,
        dftransformer_index;
        transposefeatures = dftransformer_transposefeatures,
        transposelabels = dftransformer_transposelabels,
        )
    knetestimator = KnetModel(
        ;
        name = name,
        predict_function_source = predict_function_source,
        loss_function_source = loss_function_source,
        losshyperparameters = losshyperparameters,
        optimizationalgorithm = optimizationalgorithm,
        optimizerhyperparameters = optimizerhyperparameters,
        minibatchsize = minibatchsize,
        modelweights = modelweights,
        isclassificationmodel = true,
        isregressionmodel = false,
        maxepochs = maxepochs,
        printlosseverynepochs = printlosseverynepochs,
        )
    predprobalabelfixer =
            ImmutablePredictProbaSingleLabelInt2StringTransformer(
                        1,
                        single_label_levels
                        )
    predictlabelfixer =
            ImmutablePredictionsSingleLabelInt2StringTransformer(
                        1,
                        single_label_levels
                        )
    probapackager = ImmutablePackageSingleLabelPredictProbaTransformer(
        single_label_name,
        )
    predpackager = ImmutablePackageSingleLabelPredictionTransformer(
        single_label_name,
        )
    finalpipeline = dftransformer |>
                    knetestimator |>
                    predprobalabelfixer |>
                    predictlabelfixer |>
                    probapackager |>
                    predpackager
    finalpipeline.name = name
    if !is_nothing(feature_contrasts)
        set_feature_contrasts!(finalpipeline, feature_contrasts)
    end
    return finalpipeline
end

"""
"""
function single_labelmulticlassdataframeknetclassifier(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        package::Symbol = :none,
        name::AbstractString = "",
        predict_function_source::AbstractString = "",
        loss_function_source::AbstractString = "",
        losshyperparameters::AbstractDict = Dict(),
        optimizationalgorithm::Symbol = :nothing,
        optimizerhyperparameters::AbstractDict = Dict(),
        minibatchsize::Integer = 0,
        modelweights::AbstractArray = [],
        maxepochs::Integer = 0,
        printlosseverynepochs::Integer = 0,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} =
            nothing,
        )
    if package == :Knet
        result = single_labelmulticlassdataframeknetclassifier_Knet(
            feature_names,
            single_label_name,
            single_label_levels;
            name = name,
            predict_function_source = predict_function_source,
            loss_function_source = loss_function_source,
            losshyperparameters = losshyperparameters,
            optimizationalgorithm = optimizationalgorithm,
            optimizerhyperparameters = optimizerhyperparameters,
            minibatchsize = minibatchsize,
            modelweights = modelweights,
            maxepochs = maxepochs,
            printlosseverynepochs = printlosseverynepochs,
            feature_contrasts = feature_contrasts
            )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end

"""
"""
function single_labeldataframeknetregression_Knet(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        name::AbstractString = "",
        predict_function_source::AbstractString = "",
        loss_function_source::AbstractString = "",
        losshyperparameters::AbstractDict = Dict(),
        optimizationalgorithm::Symbol = :nothing,
        optimizerhyperparameters::AbstractDict = Dict(),
        minibatchsize::Integer = 0,
        modelweights::AbstractArray = [],
        maxepochs::Integer = 0,
        printlosseverynepochs::Integer = 0,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} =
            nothing,
        )
    label_names = [single_label_name]
    dftransformer_index = 1
    dftransformer_transposefeatures = true
    dftransformer_transposelabels = true
    dftransformer = MutableDataFrame2RegressionKnetTransformer(
        feature_names,
        label_names;
        transposefeatures = true,
        transposelabels = true,
        )
    knetestimator = KnetModel(
        ;
        name = name,
        predict_function_source = predict_function_source,
        loss_function_source = loss_function_source,
        losshyperparameters = losshyperparameters,
        optimizationalgorithm = optimizationalgorithm,
        optimizerhyperparameters = optimizerhyperparameters,
        minibatchsize = minibatchsize,
        modelweights = modelweights,
        isclassificationmodel = false,
        isregressionmodel = true,
        maxepochs = maxepochs,
        printlosseverynepochs = printlosseverynepochs,
        )
    predpackager = ImmutablePackageMultiLabelPredictionTransformer(
        [single_label_name,],
        )
    finalpipeline = dftransformer |>
                    knetestimator |>
                    predpackager
    finalpipeline.name = name
    if !is_nothing(feature_contrasts)
        set_feature_contrasts!(finalpipeline, feature_contrasts)
    end
    return finalpipeline
end

"""
"""
function single_labeldataframeknetregression(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        package::Symbol = :none,
        name::AbstractString = "",
        predict_function_source::AbstractString = "",
        loss_function_source::AbstractString = "",
        losshyperparameters::AbstractDict = Dict(),
        optimizationalgorithm::Symbol = :nothing,
        optimizerhyperparameters::AbstractDict = Dict(),
        minibatchsize::Integer = 0,
        modelweights::AbstractArray = [],
        maxepochs::Integer = 0,
        printlosseverynepochs::Integer = 0,
        feature_contrasts::Union{Nothing, AbstractFeatureContrasts} =
            nothing,
        )
    if package == :Knet
        result = single_labeldataframeknetregression_Knet(
            feature_names,
            single_label_name;
            name = name,
            predict_function_source = predict_function_source,
            loss_function_source = loss_function_source,
            losshyperparameters = losshyperparameters,
            optimizationalgorithm = optimizationalgorithm,
            optimizerhyperparameters = optimizerhyperparameters,
            minibatchsize = minibatchsize,
            modelweights = modelweights,
            maxepochs = maxepochs,
            printlosseverynepochs = printlosseverynepochs,
            feature_contrasts = feature_contrasts
            )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end


text:
using Test
import GDAL
import ArchGDAL as AG

@testset "test_rasterattrtable.jl" begin
    @testset "Testing Raster Attribute Tables" begin
        AG.createRAT() do rat
            @test AG.ncolumn(rat) == 0
            @test AG.nrow(rat) == 0
            @test AG.changesarewrittentofile(rat) == false

            AG.createcolumn!(rat, "col1", AG.GFT_Integer, AG.GFU_Generic)
            @test AG.ncolumn(rat) == 1
            @test AG.columnname(rat, 0) == "col1"
            @test AG.columnusage(rat, 0) == AG.GFU_Generic
            @test AG.columntype(rat, 0) == AG.GFT_Integer
            @test AG.findcolumnindex(rat, AG.GFU_Generic) == 0
            @test AG.findcolumnindex(rat, AG.GFU_Red) == -1

            AG.createcolumn!(rat, "col2", AG.GFT_Real, AG.GFU_MinMax)
            AG.createcolumn!(rat, "col3", AG.GFT_String, AG.GFU_PixelCount)
            @test AG.ncolumn(rat) == 3
            @test AG.nrow(rat) == 0
            AG.setrowcount!(rat, 5)
            @test AG.nrow(rat) == 5

            @test AG.attributeio!(
                rat,
                AG.GF_Read,
                0,
                0,
                5,
                Array{Cint}(undef, 5),
            ) == fill(0, 5)
            @test AG.attributeio!(
                rat,
                AG.GF_Read,
                0,
                0,
                5,
                Array{Float64}(undef, 5),
            ) == fill(0, 5)
            @test AG.attributeio!(
                rat,
                AG.GF_Read,
                1,
                0,
                5,
                Array{Float64}(undef, 5),
            ) == fill(0, 5)

            @test AG.asstring(rat, 2, 0) == "0"
            @test AG.asint(rat, 2, 0) == 0
            @test AG.asdouble(rat, 2, 0) == 0

            AG.setvalue!(rat, 2, 0, "2")
            @test AG.asstring(rat, 2, 0) == "2"
            @test AG.asint(rat, 2, 0) == 2
            @test AG.asdouble(rat, 2, 0) == 2

            AG.setvalue!(rat, 2, 0, 3)
            @test AG.asstring(rat, 2, 0) == "3"
            @test AG.asint(rat, 2, 0) == 3
            @test AG.asdouble(rat, 2, 0) == 3

            AG.setvalue!(rat, 2, 0, 4.5)
            @test AG.asstring(rat, 2, 0) == "4"
            @test AG.asint(rat, 2, 0) == 4
            @test AG.asdouble(rat, 2, 0) == 4

            @test AG.asstring(rat, 2, 1) == "0"
            @test AG.asstring(rat, 2, 2) == ""

            @test AG.asstring(rat, 0, 2) == ""
            @test AG.asstring(rat, 2, 2) == ""
            @test AG.asstring(rat, 4, 2) == ""
            AG.attributeio!(rat, AG.GF_Write, 2, 0, 5, fill("abc", 5))
            @test AG.asstring(rat, 0, 2) == "abc"
            @test AG.asstring(rat, 2, 2) == "abc"
            @test AG.asstring(rat, 4, 2) == "abc"

            AG.clone(rat) do ratclone
                @test AG.asstring(ratclone, 0, 2) == "abc"
                @test AG.asstring(ratclone, 2, 2) == "abc"
                @test AG.asstring(ratclone, 4, 2) == "abc"
                @test AG.ncolumn(ratclone) == 3
                @test AG.nrow(ratclone) == 5
                @test AG.findcolumnindex(ratclone, AG.GFU_Generic) == 0
                @test AG.findcolumnindex(ratclone, AG.GFU_Red) == -1
            end

            AG.setlinearbinning!(rat, 0, 10)
            @test (@test_logs (
                :warn,
                "There is no linear binning information.",
            ) AG.getlinearbinning(rat)) == (0, 10)
            AG.setlinearbinning!(rat, -1.5, 12.0)
            @test (@test_logs (
                :warn,
                "There is no linear binning information.",
            ) AG.getlinearbinning(rat)) == (-1.5, 12.0)

            @test AG.findrowindex(rat, 0) == 0
            @test AG.findrowindex(rat, -1) == 0
            @test AG.findrowindex(rat, -1.5) == 0
            @test AG.findrowindex(rat, 7.5) == 0
            @test AG.findrowindex(rat, 12) == 1
            @test AG.findrowindex(rat, 13) == 1
        end
    end

    @testset ("Testing Color Tables") begin
        AG.createcolortable(AG.GPI_RGB) do ct
            @test sprint(print, ct) == "ColorTable[GPI_RGB]"
            @test AG.paletteinterp(ct) == AG.GPI_RGB
            AG.clone(ct) do ct1
                @test sprint(print, ct1) == "ColorTable[GPI_RGB]"
            end
            AG.clone(AG.ColorTable(C_NULL)) do ct2
                @test sprint(print, ct2) == "NULL ColorTable"
            end
            AG.createRAT(ct) do rat
                @test sprint(print, AG.toColorTable(rat, 0)) ==
                      "ColorTable[GPI_RGB]"
            end
        end
    end
end


text:
function handle_drop(files::Vector{String})
    for f in files
        try
            glplot(load(f))
        catch e
            warn(e)
        end
    end
end

toolbar_area(pa, toolbar_width) = SimpleRectangle(0, 0, toolbar_width, pa.h)
function viewing_area(area_l, area_r)
    SimpleRectangle(area_l.x+area_l.w, 0, area_r.x-area_l.w, area_r.h)
end
function edit_rectangle(visible, area, tarea)
    w = visible ? 70mm : 1.5mm
    x = area.w-w
    SimpleRectangle(x, 0, w, area.h)
end

layout_pos_ho(i) = map(icon_size) do ip
    SimpleRectangle{Float32}(0, i*ip + i*2, ip, ip)
end
layout_pos_ver(i, border) = map(icon_size) do ip
    SimpleRectangle{Float32}(i*ip + i*border, 0, ip, ip)
end

function edit_item_area(la, item_height, left_gap)
    y = la.y-item_height-2
    return SimpleRectangle(left_gap, y, la.w, item_height)
end

function save_record(frames)
    path = joinpath(homedir(), "Desktop")
    GLVisualize.create_video(frames, "test.webm", path, 0)
end

const _compute_callbacks = []
register_compute(f) = push!(_compute_callbacks, f)
poll_reactive() = (Base.n_avail(Reactive._messages) > 1) && Reactive.run_till_now()
function glplot_renderloop(window, compute_s, record_s)
    was_recording = false
    frames = []
    i = 0
    Reactive.stop()
    yield()
    while isopen(window)
        tic()
        GLWindow.pollevents()
        record = !value(record_s)
        if isready(Reactive._messages)
            Reactive.run_till_now()
            GLWindow.render_frame(window)
            GLWindow.swapbuffers(window)
        end
        if record
            push!(frames, screenbuffer(window))
        elseif was_recording && !record
            save_record(frames)
            frames = []
            gc()
        end
        yield()
        diff = (1/60) - toq()
        while diff >= 0.001
            tic()
            sleep(0.001) # sleep for the minimal amount of time
            diff -= toq()
        end
        was_recording = record
    end
    GLWindow.destroy!(window)
    GLVisualize.cleanup_old_screens()
end

function get_dpi(window)
    monitor = GLFW.GetPrimaryMonitor()
    props = GLWindow.MonitorProperties(monitor)
    props.dpi[1]# we do not start fiddling with differently scaled xy dpi's
end

immutable Millimeter
end
const mm = Millimeter()
function Base.:(*)(x::Millimeter, y::Number)
    round(Int, y * pixel_per_mm)
end
function Base.:(*)(x::Number, y::Millimeter)
    round(Int, x * pixel_per_mm)
end

function init()

    w = glscreen("GLPlot")

    preserve(map(handle_drop, w.inputs[:dropped_files]))

    global const pixel_per_mm = get_dpi(w)/25.4

    global const icon_size = Signal(10mm)
    w.inputs[:key_pressed] = const_lift(GLAbstraction.singlepressed,
        w.inputs[:mouse_buttons_pressed],
        GLFW.MOUSE_BUTTON_LEFT
    )
    button_pos = map(w.area) do a
        Point2f0[(0, a.h/2)]
    end
    button_width = 1.5mm
    edit_screen_show_button = visualize(
        (SimpleRectangle{Float32}(0, 0, button_width, button_width*2), button_pos),
        offset=Vec2f0(0, -button_width),
        color=RGBA{Float32}(0.6,0.6,0.6,1)
    )
    tarea = map(toolbar_area, w.area, icon_size)

    show_edit_screen = toggle(edit_screen_show_button, w, false)
    edit_screen_area = map(edit_rectangle,
        show_edit_screen, w.area, tarea
    )

    global const viewing_screen = Screen(w,
        name=Symbol("Viewing Screen"),
        area=map(viewing_area, tarea, edit_screen_area),
        color=RGBA{Float32}(1,1,1,1)
    )
    global const toolbar_screen = Screen(w, area=tarea)
    global const edit_screen = Screen(
        w, area=edit_screen_area,
        color=RGBA{Float32}(0.9,0.9,0.9,1)
    )
    _view(edit_screen_show_button, edit_screen, camera=:fixed_pixel)
    play_record, record_sig = toggle_button(imload("record.png"), imload("break.png"), w)
    compute_record, compute_sig = toggle_button(imload("play.png"), imload("break.png"), w)
    persp_ortho, persp_ortho_toggle_sig = toggle_button(imload("perspective.png"), imload("ortho.png"), w)
    persp_ortho_sig = map(persp_ortho_toggle_sig) do isp
        isp && return GLAbstraction.PERSPECTIVE
        GLAbstraction.ORTHOGRAPHIC
    end
    cube = cubecamera(viewing_screen, persp_ortho_sig)

    image_names = ["center", "screenshot"]
    tools = Matrix{BGRA{U8}}[imload("$name.png") for name in image_names]
    center_b, center_s = button(tools[1], w)
    screenshot_b, screenshot_s = button(tools[2], w)

    tools = [center_b, screenshot_b, play_record, persp_ortho, compute_record]
    tools_robjs = Any[]

    i = 0
    for tool in tools
        robj = layout!(layout_pos_ho(i), visualize(tool))
        _view(robj, toolbar_screen, camera=:fixed_pixel)
        push!(tools_robjs, robj.children[])
        i += 1
    end
    preserve(map(center_s) do pressed
        if pressed
            center!(viewing_screen)
        end
        nothing
    end)
    preserve(map(screenshot_s) do pressed
        if pressed
            screenshot(viewing_screen, path=joinpath(homedir(), "Desktop", "glplot.png"))
        end
        nothing
    end)
    rot = cube.children[][:model]

    cube.children[][:model] = map(rot, icon_size) do r, ip
        half = ip/2
        translationmatrix(Vec3f0(half,i*ip + i*2 + half,0))*r*scalematrix(Vec3f0(half))
    end
    _view(cube, toolbar_screen, camera=:fixed_pixel)

    GLVisualize.add_screen(viewing_screen)
    @materialize scroll, mouseposition = edit_screen.inputs
    should_scroll = map(mouseposition) do mb
        isinside(value(w.area), mb...)
    end
    scroll = filterwhen(should_scroll, value(scroll), scroll)
    edit_screen.inputs[:menu_scroll] = foldp(0, scroll) do v0, s
        v0+(ceil(Int, s[2])*15)
    end
    global _renderloop_task = @async glplot_renderloop(w, compute_sig, record_sig)
    viewing_screen
end

function block()
    global _renderloop_task
    wait(_renderloop_task)
end


text:
module CGAL

using CxxWrap
using CxxWrap.CxxWrapCore:
    CxxBaseRef,
    IsCxxType,

    cpp_trait_type,
    reference_type_union

using Requires
using libcgal_julia_jll

function __init__()
   @initcxx
   @require(Khepri="c487d830-b105-11e8-2c9d-150ba7e503f8",
            include("glue/khepri.jl"))
end

libcgal_julia() = haskey(ENV, "JLCGAL_EXACT_CONSTRUCTIONS") ?
                      libcgal_julia_exact :
                      libcgal_julia_inexact

macro singleton(name::Symbol, cname::Union{Symbol,Nothing}=nothing)
    if isnothing(cname)
        cname = Symbol(uppercase(replace(string(name), r"(?<!)(?=[A-Z])" => "_")))
    end

    quote
        struct $name end
        const $cname = $name()
        export $cname
    end |> esc
end

# Origin constants
@singleton Origin
@singleton NullVector

# Transformation tags
@singleton IdentityTransformation IDENTITY
@singleton Reflection
@singleton Rotation
@singleton Scaling
@singleton Translation

@wrapmodule get(ENV, "JLCGAL_LIBPATH", libcgal_julia())

iscxxtype(T::Type) = _iscxxtype(cpp_trait_type(T))
_iscxxtype(x) = false
_iscxxtype(::Type{IsCxxType}) = true

include("kernel.jl")
include("algebra.jl")
include("global_kernel_functions.jl")
include("convex_hull_2.jl")
include("principal_component_analysis.jl")
include("polygon_2.jl")
include("straight_skeleton_2.jl")
include("triangulation_2.jl")
include("voronoi_diagram_2.jl")

include("visual.jl")

end # module CGAL


text:
__precompile__(true)
module SugarBLAS

export  @blas!
export  @scale!, @axpy!, @copy!, @ger!, @syr!, @syrk!,
        @her!, @herk!, @gbmv!, @sbmv!, @gemm!, @gemv!,
        @symm!, @symv!

include("Match/Match.jl")
using .Match

import Base: copy, abs, -

copy(s::Symbol) = s
-(expr) = Expr(:call, :-, expr)

function abs(ast)
    if @match(ast, -ast) | (ast == 0)
        ast
    else
        Expr(:call, :(-), ast)
    end
end

substracts(expr) = false
substracts(expr::Expr) = (expr.head == :call) & (expr.args[1] == :-)

char(s::Symbol) = string(s)[1]
function char(expr::Expr)
    (expr.head == :quote) || error("char doesn't support $(expr.head)")
    char(expr.args[1])
end

isempty(nl::Nullable) = nl.isnull

wrap(expr::Symbol) = QuoteNode(expr)
function wrap(expr::Expr)
    head = QuoteNode(expr.head)
    func = string(expr.args[1])
    :(Expr($head, parse($func), $(expr.args[2:end]...)))
end

function expand(expr::Expr)
    @match(expr, A += B) && return :($A = $A + $B)
    @match(expr, A -= B) && return :($A = $A - $B)
    expr
end

macro call(expr::Expr)
    esc(:(esc($(wrap(expr)))))
end

macro switch(expr::Expr)
    (expr.head == :block) || error("@switch statement must be followed by `begin ... end`")
    lines = filter(expr::Expr -> expr.head != :line, expr.args)
    failproof(s) = s
    failproof(s::Char) = string("'",s,"'")
    line = lines[1]
    exec = "if $(line.args[1])\n$(failproof(line.args[2]))\n"
    for line in lines[2:end-1]
        (line.head == :line) && continue
        line.head == :(=>) || error("Each condition must be followed by `=>`")
        exec *= "elseif $(line.args[1])\n$(failproof(line.args[2]))\n"
    end
    line = lines[end]
    exec *= (line.args[1] == :otherwise) && ("else\n$(failproof(line.args[2]))\n")
    exec *= "end"
    esc(parse(exec))
end

###########
# Mutable #
###########

#Must be ordered from most to least especific formulas
macro blas!(expr::Expr)
    expr = expand(expr)
    @switch begin
        @match(expr, X *= a)        => @call scale!(a,X)
        @match(expr, X = a*X)       => @call scale!(a,X)
        @match(expr, Y = Y - a*X)   => @call Base.LinAlg.axpy!(-a,X,Y)
        @match(expr, Y = Y - X)     => @call Base.LinAlg.axpy!(-1.0,X,Y)
        @match(expr, Y = a*X + Y)   => @call Base.LinAlg.axpy!(a,X,Y)
        @match(expr, Y = X + Y)     => @call Base.LinAlg.axpy!(1.0,X,Y)
        @match(expr, X = Y)         => @call copy!(X, Y)
        otherwise                   => error("No match found")
    end
end

macro copy!(expr::Expr)
    @switch begin
        @match(expr, X = Y) => @call copy!(X,Y)
        otherwise           => error("No match found")
    end
end

macro scale!(expr::Expr)
    @switch begin
        @match(expr, X *= a)    => @call scale!(a,X)
        @match(expr, X = a*X)   => @call scale!(a,X)
        otherwise               => error("No match found")
    end
end

macro axpy!(expr::Expr)
    expr = expand(expr)
    @switch begin
        @match(expr, Y = Y - a*X)   => @call(Base.LinAlg.axpy!(-a,X,Y))
        @match(expr, Y = Y - X)     => @call Base.LinAlg.axpy!(-1.0,X,Y)
        @match(expr, Y = a*X + Y)   => @call Base.LinAlg.axpy!(a,X,Y)
        @match(expr, Y = X + Y)     => @call Base.LinAlg.axpy!(1.0,X,Y)
        otherwise                   => error("No match found")
    end
end

macro ger!(expr::Expr)
    expr = expand(expr)
    f = @switch begin
        @match(expr, A = alpha*x*y' + A)    => identity
        @match(expr, A = A - alpha*x*y')    => (-)
        otherwise                           => error("No match found")
    end
    @call Base.LinAlg.BLAS.ger!(f(alpha),x,y,A)
end

macro syr!(expr::Expr)
    expr = expand(expr)
    @match(expr, A[uplo] = right) || error("No match found")
    c = char(uplo)
    f = @switch begin
        @match(right, alpha*x*x.' + Y)  => identity
        @match(right, Y - alpha*x*x.')  => (-)
        otherwise                       => error("No match found")
    end
    (@match(Y, Y[uplo]) && (Y == A)) || (Y == A) || error("No match found")
    @call Base.LinAlg.BLAS.syr!(c,f(alpha),x,A)
end

macro syrk!(expr::Expr)
    expr = expand(expr)
    @match(expr, C[uplo] = right) || error("No match found")
    c = char(uplo)
    f = @switch begin
        @match(right, alpha*X*Y + D)    => identity
        @match(right, D - alpha*X*Y)    => (-)
        otherwise                       => error("No match found")
    end
    trans = @switch begin
        @match(X, A.') && (Y == A)  => 'T'
        @match(Y, A.') && (X == A)  => 'N'
        otherwise                   => error("No match found")
    end
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[uplo]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.syrk!(c,trans,f(alpha),A,beta,C)
end

macro her!(expr::Expr)
    expr = expand(expr)
    @match(expr, A[uplo] = right) || error("No match found")
    c = char(uplo)
    f = @switch begin
        @match(right, alpha*x*x' + Y)   => identity
        @match(right, Y - alpha*x*x')   => (-)
        otherwise                       => error("No match found")
    end
    (@match(Y, Y[uplo]) && (Y == A)) || (Y == A) || error("No match found")
    @call Base.LinAlg.BLAS.her!(c,f(alpha),x,A)
end

macro herk!(expr::Expr)
    expr = expand(expr)
    @match(expr, C[uplo] = right) || error("No match found")
    c = char(uplo)
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*X*Y + D)    => identity
        @match(right, D - alpha*X*Y)    => (-)
        otherwise                       => error("No match found")
    end
    trans = @switch begin
        @match(X, A') && (Y == A)   =>  'T'
        @match(Y, A') && (X == A)   =>  'N'
        otherwise                   =>  error("No match found")
    end
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[uplo]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.herk!(c,trans,f(alpha),A,beta,C)
end

macro gbmv!(expr::Expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*Y*x + w)    => identity
        @match(right, w - alpha*Y*x)    => (-)
        otherwise                       => error("No match found")
    end
    trans = @match(Y, Y') ? 'T' : 'N'
    @match(Y, A[kl:ku,h=m])
    @match(w, beta*w) || (beta = 1.0)
    (@match(w, w[uplo]) && (y == w)) || (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.gbmv!(trans,m,abs(kl),ku,f(alpha),A,x,beta,y)
end

macro sbmv!(expr::Expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*A[0:k,uplo]*x + w)  => identity
        @match(right, w - alpha*A[0:k,uplo]*x)  => (-)
        otherwise                               => error("No match found")
    end
    c = char(uplo)
    @match(w, beta*w) || (beta = 1.0)
    (@match(w, w[uplo]) && (y == w)) || (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.sbmv!(c,k,f(alpha),A,x,beta,y)
end

macro gemm!(expr::Expr)
    expr = expand(expr)
    @match(expr, C = right) || error("No match found")
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*A*B + D)    => identity
        @match(right, D - alpha*A*B)    => (-)
        otherwise                       => error("No match found")
    end
    tA = @match(A, A') ? 'T' : 'N'
    tB = @match(B, B') ? 'T' : 'N'
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[uplo]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.gemm!(tA,tB,f(alpha),A,B,beta,C)
end

macro gemv!(expr::Expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*A*x + w)  => identity
        @match(right, w - alpha*A*x)  => (-)
        otherwise                         => error("No match found")
    end
    tA = @match(A, A') ? 'T' : 'N'
    @match(w, beta*w) || (beta = 1.0)
    (@match(w, w[uplo]) && (y == w)) || (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.gemv!(tA,f(alpha),A,x,beta,y)
end

macro symm!(expr::Expr)
    expr = expand(expr)
    @match(expr, C[uplo] = right) || error("No match found")
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*A*B + D)    => identity
        @match(right, D - alpha*A*B)    => (-)
        otherwise                       => error("No match found")
    end
    c = char(uplo)
    side = @switch begin
        @match(A, A[symm]) && (symm.args[1] == :symm)   => 'L'
        @match(B, B[symm]) && (symm.args[1] == :symm)   => 'R'
        otherwise                                       => error("No match found")
    end
    @match(D, beta*D) || (beta = 1.0)
    (@match(D, D[uplo]) && (C == D)) || (C == D) || error("No match found")
    @call Base.LinAlg.BLAS.symm!(side,c,f(alpha),A,B,beta,C)
end

macro symv!(expr::Expr)
    expr = expand(expr)
    @match(expr, y = right) || error("No match found")
    f = @switch begin #Right hand side must match one of these
        @match(right, alpha*A[uplo]*x + w)  => identity
        @match(right, w - alpha*A[uplo]*x)  => (-)
        otherwise                         => error("No match found")
    end
    c = char(uplo)
    @match(w, beta*w) || (beta = 1.0)
    (@match(w, w[uplo]) && (y == w)) || (y == w) || error("No match found")
    @call Base.LinAlg.BLAS.symv!(c,f(alpha),A,x,beta,y)
end

end


text:
import MeshArrays
MeshArrays.GRID_LLC90_download()

import OceanStateEstimation
OceanStateEstimation.ECCOdiags_download()
OceanStateEstimation.ECCOdiags_add("release4")

#OceanStateEstimation.ECCOdiags_add("interp_coeffs")
OceanStateEstimation.Downloads.download(
  "https://zenodo.org/record/5784905/files/interp_coeffs_halfdeg.jld2",
  joinpath(OceanStateEstimation.ECCOdiags_path,"interp_coeffs_halfdeg.jld2");
  timeout=60000.0)

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
using GrayCoding
using Documenter

DocMeta.setdocmeta!(GrayCoding, :DocTestSetup, :(using GrayCoding); recursive=true)
push!(LOAD_PATH,"../src/")
makedocs(;
    modules=[GrayCoding],
    authors="Nivedita Rethnakar et al.",
    repo="https://github.com/nivupai/GrayCoding.jl/blob/{commit}{path}#{line}",
    sitename="GrayCoding.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nivupai.github.io/GrayCoding.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Algebra of Gray Codes" => "algebra.md",
        "Applications" => [
            "List of Applications" => "applications.md",
            "Quantum Algorithms and Circuits" => "quantum.md",
            "QAM modulation in WiFi communications" => "wireless.md",
            "DNA Coding" => "dna.md"
        ],
        "Tutorials" => "tutorials.md",
    ],
)

deploydocs(;
    repo="github.com/nivupai/GrayCoding.jl",
    devbranch="main",
)


text:
# Precursor of the RF excitation code with feature directions, most of the comments below are already implemented

using Koma
using Koma: γ, Q, RF_fun, get_grads, Un, cross
using Plots, LaTeXStrings
gr(size = (800,800))

# Example
G = 30e-3 # 30 mT/m
T = 3e-3; # 3 ms
B1 = 15e-6; # 15 μT

N = 300
RFs = RF_fun(t->B1*sinc.(4(t-T/2)/T),T,N) #sinc pulse of duration T [s]
B1e = RFs.A
α = π / 2
ΔT = RFs.T
c = α / (2π*γ*sum(B1e.*ΔT))
c = c.re
RFs.A .= RFs.A * c

plot(real.(RFs.A), linewidth=5, dpi=200)

## In general, Bz = Δω(x)/γ + G ⋅ x
B1e = RFs.A
ΔT = RFs.T
Bx, By = real.(B1e), imag.(B1e)

# Plot
df = .2 #.02
ff = -1:df:1
colors = range(HSV(0,1,1), stop=HSV(-360,1,1), length=length(ff)) # inverse rotation

anim = @animate for j ∈ reverse(2:N)
plot()
for (i,f) = enumerate(ff)
Bz = zeros(size(B1e)) .+ f*.01e-3 # Bz = 0 in the rotating frame implies a non-selective excitation

φ = -2π*γ * ΔT .* sqrt.(abs.(B1e).^2 .+ abs.(Bz).^2) # angle of rotation
φ[φ.==0] .= 1e-17; # removes problems when dividing by φ
n =  2π*γ * ΔT .* [Bx By Bz]./abs.(φ) # axis of rotation
Qs = [Q(φ[i],n[i,:]) for i=1:N] # hard-pulse approximation for every RF element in RFs
Qt = *(Qs...) # Total rotation matrix

Mxy, Mz = 0. + 0. *im, 1. # [0,0,1]
M0 = Mag(Mxy, Mz) # Initial magnetization
Mt = Qt*M0 # Final magnetization
Mz = [(*(Qs[1:i]...)*M0)[2] for i=2:N-j]
Mxy = [(*(Qs[1:i]...)*M0)[1] for i=2:N-j]

# Plots
plot!(real.(Mxy),imag.(Mxy),real.(Mz),
    linewidth=5,xlabel="x",ylabel="y",zlabel="z",
    xlims=(-1,1),ylims=(-1,1),zlims=(-1,1),dpi=200,label=:none,
    linecolor=colors[i],subplot=1)
end

plot!(real.(RFs.A)[1:N-j], linewidth=5, dpi=200,
    xlims=(0,N),ylims=(-.5c*B1,1.1c*B1), ticks = nothing,
    inset_subplots=bbox(0.05, 0.05, 0.5, 0.25, :bottom, :right),
    subplot=2,label=:none)

end
#gif(anim, "anim_fps60.gif", fps = 60)

## TODO
# - A function that translates from Phantom and Sequence to a rotation matrix for each spin, something like:
# function get_rf_rotations(obj::Phantom, seq::Sequence, t)
#     Δt = t[2]-t[1]
#     Gx, Gy = get_grads(seq,t)
#     B1e = get_rfs(seq,t) #TODO: uniformly sample RF pulses with respect to the simulation sampling period Δt
#     for spin ∈ obj # this an example, it will not work
#         φ = -2π*γ * Δt .* sqrt.(abs.(B1e).^2 .+ abs.(Bz).^2) # angle of rotation
#         φ[φ.==0] .= 1e-17; # removes problems when dividing by φ
#         n =  2π*γ * Δt .* [Bx By Bz]./abs.(φ) # axis of rotation
#         Qtotal = prod([Q(φ[i],n[i,:]) for i=1:N]) # hard-pulse approximation for every RF element in RFs
#     end
#     Qtotal
# end
# - Obtaining rotation matrices efficiently: as we would NEVER acquire during an RF pulse we do not need
#   to store intermidiate steps. This is due to the fact that we need to change between Rx and Tx with a switch (hardware limitation).
# - I suggest to simulate RF blocks separatly, and then to use run_sim2D_spin() in between.
#   The only thing we need to add is a way to follow the T1 decay during the non-RF parts of the simulation.
#   The magnetization [Mx + iMx, Mz] should decay as follows (the first part is already in run_sim2D_spin()):
#                    [exp(-Δt/T2)(Mx + iMx), exp(-Δt/T1) Mz + ρ (1 - exp(-Δt/T1))].
# - Movement of the phantom during excitation? (using obj.x + obj.ux(obj.x,obj.y))
# - Effects of T1 and T2 during an RF block: For each Δt we do RF and then T1 & T2 decay.
# - Plotting function Koma.plot_grads(seq) in Display.jl is not currently plotting RF pulses
# - Pre-defined RF waveforms.
# - Shinnar-Le Roux pulse desing?


text:
using PyBaMM
using Test
using SafeTestsets
using SparseArrays, LinearAlgebra

@safetestset "Compare models with PyBaMM" begin include("test_full_models.jl") end
@safetestset "Test loss functions" begin include("test_loss.jl") end


text:
using PATHSolver
using Base.Test

M = [0  0 -1 -1 ;
     0  0  1 -2 ;
     1 -1  2 -2 ;
     1  2 -2  4 ]

q = [2; 2; -2; -6]

myfunc(x) = M*x + q

n = 4
lb = zeros(n)
ub = 100*ones(n)

path_options(   "convergence_tolerance 1e-2",
                "output yes",
                "time_limit 3600"      )

status, z, f = solveMCP(myfunc, lb, ub)

@show status
@show z
@show f

@test isapprox(z, [2.8, 0.0, 0.8, 1.2])
@test status == :Solved

println("-------------------------------------------------------")

M = [0  0 -1 -1 ;
     0  0  1 -2 ;
     1 -1  2 -2 ;
     1  2 -2  4 ]

q = [2; 2; -2; -6]

myfunc(x) = M*x + q

n = 4
lb = zeros(n)
ub = 100*ones(n)

path_options(   "convergence_tolerance 1e-2",
                "output no",
                "time_limit 3600"      )

status, z, f = solveMCP(myfunc, lb, ub)

@show status
@show z
@show f

@test isapprox(z, [2.8, 0.0, 0.8, 1.2])
@test status == :Solved

println("-------------------------------------------------------")

function elemfunc(x)
    val = similar(x)
    val[1] = -x[3]-x[4] + q[1]
    val[2] = x[3] -2x[4] + q[2]
    val[3] = x[1]-x[2]+2x[3]-2x[4] + q[3]
    val[4] = x[1]+2x[2]-2x[3]+4x[4] + q[4]
    return val
end

n = 4
lb = zeros(n)
ub = 100*ones(n)

status, z, f = solveMCP(elemfunc, lb, ub)

@show status
@show z
@show f

@test isapprox(z, [2.8, 0.0, 0.8, 1.2])
@test status == :Solved


text:
# test compute fft

T = Float32
fs = 100.  # sampling rate
N = 2^14   # number of points
Nwin = 10
freqmin = 10. # minimum frequency
freqmax = 20. # maximum frequency
corners = 4   # number of corners
zerophase = true # do zerophase filtering
starttime = d2u(DateTime(Date(now()))) # starttime for SeisChanel / SeisData

# test process raw
@testset "process raw" begin

    S = SeisData(2)
    S.x[1] = rand(T,N)
    S.id[1] = "EO.SYT.02.LLE"
    S.fs[1] = fs
    S.t[1] = [1 starttime * 1e6;N 0]
    S.x[2] = rand(T,N)
    S.id[2] = "MZ.KTFD.08.ENZ"
    S.fs[2] = fs
    S.t[2] = [1 (starttime + 0.004) * 1e6;N 0]

    # test without downsampling
    Snew = process_raw(S,fs)
    @test size(Snew.x[1]) == size(S.x[1])
    @test eltype(Snew.x[1]) == eltype(S.x[1])
    @test size(Snew.x[2]) == size(S.x[2])
    @test eltype(Snew.x[2]) == eltype(S.x[2])
    @test Snew.fs == S.fs
    @test Snew.t[1] == S.t[1]
    @test Snew.t[2][1,2] == S.t[2][1,2] - 0.004 * 1e6 # test phase shift

    # test with integer sampling rate
    SInt = process_raw(S,Int(fs))
    @test size(SInt.x[1]) == size(S.x[1])
    @test eltype(SInt.x[1]) == eltype(S.x[1])
    @test size(SInt.x[2]) == size(S.x[2])
    @test eltype(SInt.x[2]) == eltype(S.x[2])
    @test SInt.fs == S.fs
    @test SInt.t[1] == S.t[1]
    @test SInt.t[2][1,2] == S.t[2][1,2] - 0.004 * 1e6 # test phase shift

    # test with downsampling
    Snew = process_raw(S,20.)
    @test size(Snew.x[1]) != size(S.x[1])
    @test eltype(Snew.x[1]) == eltype(S.x[1])
    @test size(Snew.x[2]) != size(S.x[2])
    @test eltype(Snew.x[2]) == eltype(S.x[2])
    @test Snew.fs == [20.,20.]
    @test Snew.t[1][2,1] == length(Snew.x[1])
    @test Snew.t[2][2,1] == length(Snew.x[2])
    @test Snew.t[1] != S.t[1]
    @test Snew.t[2][1,2] == S.t[2][1,2] - 0.004 * 1e6 # test phase shift

    # test without ϕshift
    Sphi = process_raw(S,20.,ϕshift=false)
    @test Sphi.x[1] == Snew.x[1]
    @test Sphi.x[2] != Snew.x[2]
    @test Sphi.t[2] == Snew.t[2]

    # test in-place
    process_raw!(S,20.)
    @test Snew == S

    # test with SeisChannel
    C = SeisChannel()
    C.x = rand(T,N)
    C.fs = fs
    C.t = [1 (starttime + 0.004) * 1e6;N 0]

    # test no downsample
    Cnew = process_raw(C,fs)
    @test size(Cnew.x) == size(C.x)
    @test eltype(Cnew.x) == eltype(C.x)
    @test Cnew.fs == C.fs
    @test Cnew.t[1,2] == C.t[1,2] - 0.004 * 1e6 # test phase shift

    # test with int sampling rate
    CInt = process_raw(C,Int(fs))
    @test size(CInt.x) == size(C.x)
    @test eltype(CInt.x) == eltype(C.x)
    @test CInt.fs == C.fs
    @test CInt.t[1,2] == C.t[1,2] - 0.004 * 1e6 # test phase shift

    # test downsample
    Cnew = process_raw(C,20.)
    @test size(Cnew.x) != size(C.x)
    @test eltype(Cnew.x) == eltype(C.x)
    @test Cnew.fs == 20.
    @test Cnew.t[1,2] == C.t[1,2] - 0.004 * 1e6 # test phase shift

    # test in-place
    process_raw!(C,20.)
    @test Cnew == C
end

# test FFTData creation with rfft
@testset "rfft" begin

    cc_len = 10.24
    cc_step = 5.12
    C = SeisChannel()
    C.x = rand(T,N)
    C.fs = fs
    C.id = "MZ.KTFD.08.ENZ"
    C.t = [1 (starttime + 0.004) * 1e6;N 0]
    R = RawData(C,cc_len,cc_step)
    R.gain = 1.234234e6
    F = rfft(R)
    @test F.name == C.id == R.name
    @test F.id == R.id
    @test F.loc == R.loc
    @test F.fs == R.fs
    @test F.gain == R.gain
    @test F.freqmin == R.freqmin
    @test F.freqmax == F.freqmax
    @test F.cc_len == R.cc_len == cc_len
    @test F.cc_step == R.cc_step == cc_step
    @test F.whitened == R.whitened
    @test F.time_norm == R.time_norm
    @test F.resp == R.resp
    @test F.t == R.t
    @test length(F.t) == size(F.fft,2)
end

@testset "phase" begin

    # test with all zero row
    A = rand(T,N,Nwin)
    A[:,1] .= T(0)
    p = phase(A)
    @test all(.!SeisNoise.isweird.(p)) # test no NaNs

    # test phase
    A = rand(T,N,Nwin)
    p = phase(A)
    @test size(A) == size(p)
    @test all(p[N ÷ 2 + 2 + isodd(N),1] .== T(0))

    # test phase RawData
    R = RawData()
    R.x = A
    R.fs = fs
    F = phase(R)
    @test isa(F,FFTData)
    @test all(.!SeisNoise.isweird.(F.fft))
    @test F.fft == p
end

@testset "mute" begin

    # generate random signal with zero mean
    A = rand(T,N,Nwin)
    A .-= mean(A,dims=1)
    A[1000:1200,:] .*= 1000
    M = mute(A)
    @test size(M) == size(A)
    @test eltype(M) == eltype(A)

    # test begin , middle, end changes
    @test all(M[1:500,:] .== A[1:500,:])
    @test median(M[1000:1200,:]) == T(0)
    @test all(M[end-500:end,:] .== A[end-500:end,:])

    # test in-place
    mute!(A)
    @test A == M

    # test RawData
    R = RawData()
    A = rand(T,N,Nwin)
    A .-= mean(A,dims=1)
    A[1000:1200,:] .*= 1000
    R.x = A
    M = mute(R)
    @test isa(M,RawData)
    @test size(M.x) == size(A)
    @test eltype(M.x) == eltype(A)

    # test in-place
    mute!(R)
    @test R.x == M.x
end

@testset "clip/clamp" begin
    Nspikes = 100
    A = rand(T,N,Nwin) .- 0.5
    ind = rand(1:N,Nspikes)
    A[ind,:] .= randn(Nspikes,Nwin) .* 100000
    A[:,1] .= T(0)

    # test clipping
    B = clip(A,3)
    @test size(B) == size(A)
    @test eltype(B) == eltype(A)
    @test B != A
    @test all(maximum(abs.(B),dims=1) .<= maximum(abs.(A),dims=1))

    # test in-place
    clip!(A,3)
    @test A == B

    # test 1D array
    A = rand(T,N) .- 0.5
    ind = rand(1:N,Nspikes)
    A[ind] .= randn(Nspikes) .* 100000
    B = clip(A,3)
    @test size(B) == size(A)
    @test eltype(B) == eltype(A)
    @test all(maximum(abs.(B),dims=1) .<= maximum(abs.(A),dims=1))

    # test RawData
    Nspikes = 100
    A = rand(T,N,Nwin) .- 0.5
    ind = rand(1:N,Nspikes)
    A[ind,:] .= randn(Nspikes,Nwin) .* 100000
    A[:,1] .= T(0)
    R = RawData()
    R.x = A
    Rclip = clip(R,3)
    @test isa(Rclip,RawData)
    @test size(Rclip.x) == size(R.x)
    @test eltype(Rclip.x) == eltype(R.x)
    @test Rclip.x != R.x
    @test all(maximum(abs.(Rclip.x),dims=1) .<= maximum(abs.(R.x),dims=1))

    # test in-place
    clip!(R,3)
    @test R.x == Rclip.x

    # test clamping w RawData
    A = rand(T,N,Nwin) .- 0.5
    ind = rand(1:N,Nspikes)
    A[ind,:] .= randn(Nspikes,Nwin) .* 100000
    A[:,1] .= T(0)
    R = RawData()
    R.x = deepcopy(A)
    Rclamp = clamp(R,π)
    @test isa(Rclamp,RawData)
    @test size(Rclamp.x) == size(R.x)
    @test eltype(Rclamp.x) == eltype(R.x)
    @test all(abs.(Rclamp.x) .<= π)
    @test Rclamp.x != R.x

    # test in place
    clamp!(R,π)
    @test R.x == Rclamp.x

    # test clamping with two values
    R.x = deepcopy(A)
    Rclamp = clamp(R,-0.75,0.5)
    @test isa(Rclamp,RawData)
    @test size(Rclamp.x) == size(R.x)
    @test eltype(Rclamp.x) == eltype(R.x)
    @test all(Rclamp.x .>= -0.75)
    @test all(Rclamp.x .<= 0.5)

    # test switching hi and lo
    @test_throws ArgumentError clamp(R,0.5,-0.75)

    # test in-place
    clamp!(R,-0.75,0.5)
    @test R.x == Rclamp.x
end

@testset "remove_amp" begin
    # create array with std == 1
    A = rand(T,N,Nwin) .- 0.5
    A ./= std(A,dims=1)

    # add amplitude spikes
    Nspikes = 100
    cols = unique(rand(1:Nwin,Nwin÷2))
    ind = rand(1:N,Nspikes)
    R = RawData()
    R.x = deepcopy(A)
    R.x[ind,cols] .= randn(Nspikes,length(cols)) .* 100000
    R.t = collect(1.:10.)
    Ramp = remove_amp(R)
    @test isa(Ramp,RawData)
    @test size(Ramp.x,2) == Nwin - length(cols)
    @test length(Ramp.t) == Nwin - length(cols)

    # test in-place
    remove_amp!(R)
    @test R.x == Ramp.x

    # add amplitude spikes to all columns
    R.x = deepcopy(A)
    R.x[ind,:] .= randn(Nspikes) .* 100000
    @test_throws ErrorException remove_amp(R)

    # test with zeros
    cols = unique(rand(1:Nwin,Nwin÷2))
    R.x = deepcopy(A)
    R.x[:,cols] .= T(0)
    R.t = collect(1.:10.)
    Ramp = remove_amp(R)
    @test isa(Ramp,RawData)
    @test size(Ramp.x,2) == Nwin - length(cols)
    @test length(Ramp.t) == Nwin - length(cols)

    # test with all zeros
    R.x .= zero(T)
    @test_throws ErrorException remove_amp(R)
end

@testset "one-bit" begin
    A = rand(T,N,Nwin) .- T(0.5)
    R = RawData()
    R.x = A

    # one-bit RawData
    R1 = onebit(R)
    @test isa(R1,RawData)
    @test maximum(R1.x) == T(1)
    @test size(R1.x) == size(A)
    @test eltype(R1.x) == eltype(A)
    @test minimum(R1.x) == T(-1)
    @test maximum(R1.x) == T(1)

    # test in-place
    onebit!(R)
    @test R.x == R1.x

end

@testset "nonzero" begin
    A = rand(T,N,Nwin)
    cols = unique(rand(1:Nwin,Nwin÷2))
    A[:,cols] .= T(0)
    ind = SeisNoise.nonzero(A)
    # test no intersect between cols and ind
    @test setdiff(cols,ind) == cols
    # test ind are Int
    @test all(isa.(ind,Int))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

content =
    """
    Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor
    incididunt ut labore et dolore magna aliqua. Arcu non sodales neque sodales.
    Placerat orci nulla pellentesque dignissim enim sit amet venenatis. Mauris
    augue neque gravida in fermentum et sollicitudin. Amet venenatis urna cursus
    eget. Enim praesent elementum facilisis leo vel fringilla est. Vitae sapien
    pellentesque habitant morbi tristique. Ornare lectus sit amet est placerat in.
    Leo urna molestie at elementum eu facilisis. Aliquam vestibulum morbi blandit
    cursus risus at ultrices. Id aliquet lectus proin nibh. Facilisi etiam
    dignissim diam quis enim lobortis scelerisque fermentum. Pretium lectus quam id
    leo in vitae turpis massa sed. Elementum facilisis leo vel fringilla est.
    Vulputate ut pharetra sit amet aliquam. Quis enim lobortis scelerisque
    fermentum dui faucibus in ornare. Cursus turpis massa tincidunt dui ut.

    A arcu cursus vitae congue mauris rhoncus. Tellus rutrum tellus pellentesque
    eu. Fringilla phasellus faucibus scelerisque eleifend donec pretium. Aliquam
    etiam erat velit scelerisque. Volutpat lacus laoreet non curabitur gravida.
    Felis imperdiet proin fermentum leo vel orci. Viverra tellus in hac habitasse
    platea dictumst vestibulum rhoncus est. Ullamcorper dignissim cras tincidunt
    lobortis feugiat vivamus. Sit amet luctus venenatis lectus. Odio facilisis
    mauris sit amet massa vitae tortor condimentum. Purus sit amet volutpat
    consequat mauris nunc congue. Enim nunc faucibus a pellentesque sit amet. Purus
    non enim praesent elementum facilisis leo vel fringilla est.
    """ |> strip

let p = Pager(content)
    @test p.pagesize == 10
    @test length(p.lines) == 22
    @test startswith(content, p.lines[1])
    @test endswith(content, p.lines[end])
    buffer = IOBuffer()
    TerminalMenus.printmenu(buffer, p, 1)
    str = String(take!(buffer))
    @test contains(str, "(10 / 22)  45%")
    @test endswith(str, "leo in vitae turpis massa sed. Elementum facilisis leo vel fringilla est.")
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
abstract type AbstractClick{T,N} end

struct ClickPointer{T,N} <:  AbstractClick{T,N}
    series::AbstractArray{T,N}
    interval::ClosedInterval
    threshold::Real
end

series(c::ClickPointer) = c.series
interval(c::ClickPointer) = c.interval
left(c::ClickPointer) = interval(c).left
right(c::ClickPointer) = interval(c).right
SampledSignals.samplerate(c::ClickPointer) = samplerate(series(c))
samples(c::ClickPointer) = Array(series(c)[left(c)*s..right(c)*s])
times(c::ClickPointer) = left(c):(1/samplerate(c)):right(c)

struct Click{T,N} <: AbstractClick{T,N}
    samples::AbstractArray{T,N}
    samplerate::Real
    left::Real
    right::Real
end
Click(c::ClickPointer) = Click(samples(c), samplerate(c), left(c),  right(c))

left(c::Click) = c.left
right(c::Click) = c.right
SampledSignals.samplerate(c::Click) = c.samplerate
samples(c::Click) = c.samples
times(c::Click) = left(c):(1/samplerate(c)):right(c)

"""
    teager(x::AbstractArray)

Calculate the Teager energy operator of a sampled signal `x`, in the form of a
`Vector` or n x 1 `Array`.  The Teager energy at sample \$x_i\$ is defined as

```math
x_i^2 - x_{i-1} x_{i+1},
```

so `x` must be at least 3 samples long.
"""
function teager(x::AbstractArray)
    @assert length(x) >= 3 "Not enough measurements to compute Teager energy"
    output = zeros(size(x))
    for i in 2:length(x)-1
        output[i] = x[i]^2 - x[i-1] * x[i+1]
    end
    output[1] = output[2]
    output[end] = output[end-1]
    return output
end

"""
    merge_intervals(ints::Vector{I}, sort_first=true) where {I <: ClosedInterval}

Combine intervals which overlap.  Takes a vector `ints` of `ClosedIntervals`. If
`sort_first=true` (the default), `ints` is sorted by the left boundary of each
interval.  The only time `sort_first` should be `false` is in the case where `ints`
is already in ascending order.
"""
function merge_intervals(ints::Vector{I}, sort_first=true) where {I <: ClosedInterval}
    if length(ints) < 1
        return I[]
    end

    if sort_first
        sorted_ints = sort(ints, by = x -> x.left)
    end
    current_interval = sorted_ints[1]
    merged_intervals = I[]
    i = 2
    while i <= length(sorted_ints)
        if isempty(intersect(current_interval, sorted_ints[i]))
            push!(merged_intervals, current_interval)
            current_interval = sorted_ints[i]
        else
            current_interval = union(current_interval, sorted_ints[i])
        end
        i += 1
    end
    push!(merged_intervals, current_interval)
    return merged_intervals
end

"""
    detect_clicks(x::AbstractArray{T,N}, samplerate, thresh, span)
    detect_clicks(x::AbstractSampleBuf{T,N}, thresh, span)

Detect clicks in a sampled signal using the Teager energy operator. `x` can be an
`AbstractArray` or `AbstractSampleBuf` (from SampledSignals.jl).  If it is the
former, a samplerate must be supplied as well.

`thresh` is the detection threshold for the Teager energy.  `span` is the time
before and after each sample above the Teager threshold to include in the click.
The duration of the shortest click will therefore be `2 * thresh`.

Returns a `Vector{ClickPointer}`.
"""
function detect_clicks(x::AbstractArray{T,N}, samplerate, thresh, span)  where {T, N}
    y = SampleBuf(Array(x), samplerate)
    return detect_clicks(y, thresh, span)
end

function detect_clicks(x::AbstractSampleBuf{T,N}, thresh, span) where {T, N}
    nchannels(x) == 1 && ArgumentError("The SampleBuff cannot have more than one channel.")
    x_teager = teager(x)
    high_teager_idx = findall(x -> x > thresh, x_teager)
    times = domain(x)
    t0, t1 = times[1], times[end]
    intervals = [max(t0, times[i]-span)..min(times[i]+span, t1) for i in high_teager_idx]
    intervals = merge_intervals(intervals)
    return [ClickPointer(x, i, thresh) for i in intervals]
end
#
# function peak_lead_ratio(click::ClickPointer, nsamples_lead)
#     t, x = get_data(click)
#     x² = abs2.(x)
#     return mean(x²[nsamples_lead:end]) / mean(x²[1:nsamples_lead])
# end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using LinearAlgebra
dot2(v) = dot(v, v)
using Plots

n = 9
z = exp(2π*im/n)
t = range(0, 2π; length=101)
v0 = normalize([z^k for k in 0:n-1])
#v0 = normalize!(rand(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))
plot(t, f.(t); label="f(t)", ylim=(0, 1))
plot!(t, g.(t); label="g(t)", ls=:auto)

# %%
using LinearAlgebra
dot2(v) = dot(v, v)
using Plots

n = 9
z = exp(2π*im/n)
t = range(0, 2π; length=101)
#v0 = normalize([z^k for k in 0:n-1])
v0 = normalize!(rand(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))
plot(t, f.(t); ylim=(0, 1))
plot!(t, g.(t); ls=:auto)

# %%
using LinearAlgebra
using Plots

n = 10
v0 = normalize(randn(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))

s0 = sum(z^2 for z in v0)
n0 = norm(v0)^2
F(t) = (n0 + real(exp(im*2t) * s0))/2
G(t) = (n0 - real(exp(im*2t) * s0))/2

t = range(0, 2π; length=201)
P = plot(t, f.(t); label="f(t)", lw=2)
plot!(t, F.(t); label="F(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(t, g.(t); label="g(t)", lw=2)
plot!(t, G.(t); label="G(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(; legend=:outertopright)

# %%
using LinearAlgebra
using Plots

n = 10
v0 = normalize(randn(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))

s0 = sum(z^2 for z in v0)
n0 = norm(v0)^2
F(t) = (n0 + real(exp(im*2t) * s0))/2
G(t) = (n0 - real(exp(im*2t) * s0))/2

t = range(0, 2π; length=201)
P = plot(t, f.(t); label="f(t)", lw=2)
plot!(t, F.(t); label="F(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(t, g.(t); label="g(t)", lw=2)
plot!(t, G.(t); label="G(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(; legend=:outertopright)

# %%
using LinearAlgebra
using Plots

n = 10
v0 = normalize(randn(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))

s0 = sum(z^2 for z in v0)
n0 = norm(v0)^2
F(t) = (n0 + real(exp(im*2t) * s0))/2
G(t) = (n0 - real(exp(im*2t) * s0))/2

t = range(0, 2π; length=201)
P = plot(t, f.(t); label="f(t)", lw=2)
plot!(t, F.(t); label="F(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(t, g.(t); label="g(t)", lw=2)
plot!(t, G.(t); label="G(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(; legend=:outertopright)

# %%
using LinearAlgebra
using Plots

n = 10
v0 = normalize(randn(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))

s0 = sum(z^2 for z in v0)
n0 = norm(v0)^2
F(t) = (n0 + real(exp(im*2t) * s0))/2
G(t) = (n0 - real(exp(im*2t) * s0))/2

t = range(0, 2π; length=201)
P = plot(t, f.(t); label="f(t)", lw=2)
plot!(t, F.(t); label="F(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(t, g.(t); label="g(t)", lw=2)
plot!(t, G.(t); label="G(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(; legend=:outertopright)

# %%
using LinearAlgebra
using Plots

n = 10
v0 = normalize(randn(ComplexF64, n))
v(t) = exp(im*t)*v0
f(t) = real(dot2(real(v(t))))
g(t) = real(dot2(imag(v(t))))

s0 = sum(z^2 for z in v0)
n0 = norm(v0)^2
F(t) = (n0 + real(exp(im*2t) * s0))/2
G(t) = (n0 - real(exp(im*2t) * s0))/2

t = range(0, 2π; length=201)
P = plot(t, f.(t); label="f(t)", lw=2)
plot!(t, F.(t); label="F(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(t, g.(t); label="g(t)", lw=2)
plot!(t, G.(t); label="G(t)", ls=:auto, lw=2, ylim=(0, 1))
plot!(; legend=:outertopright)

# %%

max_stars_repo_path: 0009/squared norm of the real part of a complex vector.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using MatrixChainMultiply
using Test, LinearAlgebra

# case on readme
function readme_example(use_arrayfire=false)
  if use_arrayfire
    T = AFArray{Float64}
  else
    T = Float64
  end
  a = rand(T, 1000, 1000)
  b = rand(T, 1000, 100)
  c = rand(T, 100, 500)
  d = rand(T, 500)

  # get results
  result1 = *(a,b,c,d)
  result2 = matrixchainmultiply(a,b,c,d)
  difference = result1 - result2

  # do again to get timings
  t1 = time()
  *(a,b,c,d);
  stdtime = time() - t1
  t2 = time()
  matrixchainmultiply(a,b,c,d);
  chaintime = time() - t2
  speedup = stdtime / chaintime
  print("Speedup = $speedup")

  difference, speedup
end

difference, speedup = readme_example()

@test norm(difference) < 1e-5
@test speedup > 1.0