text:
mutable struct GeoArray{T<:Union{Real, Union{Missing, Real}}} <: AbstractArray{T, 3}
    A::AbstractArray{T, 3}
    f::AffineMap
    crs::AbstractString
end
GeoArray(A::AbstractArray{T, 3} where T<:Union{Real, Union{Missing, Real}}) = GeoArray(A, geotransform_to_affine(SVector(0.,1.,0.,0.,0.,1.)), "")
GeoArray(A::AbstractArray{T, 2} where T<:Union{Real, Union{Missing, Real}}) = GeoArray(reshape(A, size(A)..., 1), geotransform_to_affine(SVector(0.,1.,0.,0.,0.,1.)), "")

Base.size(ga::GeoArray) = size(ga.A)
Base.IndexStyle(::Type{T}) where {T<:GeoArray} = IndexLinear()
Base.getindex(ga::GeoArray, i::Int) = getindex(ga.A, i)
Base.getindex(ga::GeoArray, I::Vararg{Int, 2}) = getindex(ga.A, I..., :)
Base.getindex(ga::GeoArray, I::Vararg{Int, 3}) = getindex(ga.A, I...)
Base.iterate(ga::GeoArray) = iterate(ga.A)
Base.length(ga::GeoArray) = length(ga.A)
Base.parent(ga::GeoArray) = ga.A
# Base.map(f, ga::GeoArray) = GeoArray(map(f, ga.A), ga.f, ga.crs)
# Base.convert(::Type{Array{T, 3}}, A::GeoArray{T}) where {T} = convert(Array{T,3}, ga.A)
Base.eltype(::Type{GeoArray{T}}) where {T} = T

function Base.show(io::IO, ga::GeoArray)
    print(io, "$(join(size(ga), "x")) $(typeof(ga.A)) with $(ga.f) and WKT $(ga.crs)")
end
function Base.show(ga::GeoArray)
    print("$(join(size(ga), "x")) $(typeof(ga.A)) with $(ga.f) and WKT $(ga.crs)")
end

# Generate upper left coordinates for specic index
function coords(ga::GeoArray, p::SVector{2, Int})
    ga.f(p.-1)
end
coords(ga::GeoArray, p::Vector{Int}) = coords(ga, SVector{2}(p))

# Generate center coordinates for specific index
function centercoords(ga::GeoArray, p::SVector{2, Int})
    ga.f(p.-0.5)
end
centercoords(ga::GeoArray, p::Vector{Int}) = centercoords(ga, SVector{2}(p))

# Convert coordinates back to indices
function indices(ga::GeoArray, p::SVector{2, Float64})
    map(x->round(Int, x), inv(ga.f)(p)::SVector{2, Float64}).+1
end
indices(ga::GeoArray, p::Vector{Float64}) = indices(ga, SVector{2}(p))
indices(ga::GeoArray, p::Tuple{Float64, Float64}) = indices(ga, SVector{2}(p))

# Overload indexing directly into GeoArray
# TODO This could be used for interpolation instead
function Base.getindex(ga::GeoArray, I::SVector{2, Float64})
    (i, j) = indices(ga, I)
    return ga[i, j, :]
end
Base.getindex(ga::GeoArray, I::Vararg{Float64, 2}) = Base.getindex(ga, SVector{2}(I))

# Generate coordinates for complete GeoArray
function coords(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [coords(ga, SVector{2}(i,j)) for i in 1:ui+1, j in 1:uj+1]
end
function centercoords(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [centercoords(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj]
end
function centercoordsnotmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [centercoords(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj if ~ismissing(ga.A[i, j])]
end
function centercoordsmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [centercoords(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj if ismissing(ga.A[i, j])]
end
function indexmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [[i,j] for i in 1:ui, j in 1:uj if ismissing(ga.A[i, j])]
end


text:
function Spectrum(w::Wave)

    amps = rfft(w.amps)
    fs = rfftfreq(length(w.amps), w.framerate)
    power = abs2.(amps)
    d = length(w.amps)

    return Spectrum(amps, fs, power,w.signal, w.framerate, d)

end

function Wave(spec::Spectrum)

    amps = irfft(spec.amps, spec.d)
    dt = inv(spec.framerate)
    t = range(0,step=dt,length = length(amps))
    #t = 0:dt:spec.d*dt
    return Wave(amps, t, spec.framerate, spec.signal)

end

# check type interference!
function +(w::Wave, wn::WhiteNoise)

    r = zeros(size(w.amps)) + wn

    return Wave(w.amps + r, w.time, w.framerate, w.signal)

end

function +(w::Wave, pn::PinkNoise)

    r = zeros(size(w.amps)) + pn

    return Wave(w.amps + r, w.time, w.framerate, w.signal)

end

function +(w::Wave, bn::BrownianNoise)

    r = zeros(size(w.amps)) + bn

    return Wave(w.amps + r, w.time, w.framerate, w.signal)

end

function +(w::Wave, gn::GaussianNoise)

    r = zeros(size(w.amps)) + gn

    return Wave(w.amps + r, w.time, w.framerate, w.signal)

end

function slope(spec::Spectrum)

    x = spec.fs[2:end] .|> log10
    y = spec.power[2:end] .|> log10
    X = [ones(length(x)) x]
    return X\y

end

Wave(wn::WhiteNoise;t=0.1,framerate=10_000) = wn(t,framerate)

Wave(gn::GaussianNoise;t=0.1,framerate=10_000) = gn(t,framerate)

Wave(bn::BrownianNoise;t=0.1,framerate=10_000) = bn(t,framerate)

Wave(pn::PinkNoise;t=0.1,framerate=10_000) = pn(t,framerate)

ext: jl


text:
# Autogenerated wrapper script for Zlib_jll for x86_64-linux-gnu
export libz

JLLWrappers.@generate_wrapper_header("Zlib")
JLLWrappers.@declare_library_product(libz, "libz.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libz,
        "lib/libz.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module PolynomialsTests

using Test

@testset "MonomialBases" begin include("MonomialBasesTests.jl") end

@testset "QGradMonomialBases" begin include("QGradMonomialBasesTests.jl") end

@testset "QCurlGradMonomialBases" begin include("QCurlGradMonomialBasesTests.jl") end

@testset "ChangeBasis" begin include("ChangeBasisTests.jl") end

end # module


text:
using Plots
push!(LOAD_PATH, "./src")
using Basis2DTri
using Basis1D
using LaTeXStrings
theme(:wong)
gr(aspect_ratio=1,legend=nothing,axis=nothing,border=nothing,ticks=nothing,
   markerstrokewidth=0,markersize=1)

r1D, w1D = gauss_lobatto_quad(0,0,4)

plot([-1;-1],[-1;1],seriescolor=:black,lw=3,label = nothing)
plot!([-1;1],[1;1],seriescolor=:black,lw=3,label = nothing)
plot!([1;1],[1;-1],seriescolor=:black,lw=3,label = nothing)
plot!([1;-1],[-1;-1],seriescolor=:black,lw=3,label = nothing)

count = 0
for i = (1:3:8*3).+8
    j = 10
    plot!([-1;-1].+i,[-1;2].+j,seriescolor=:black,lw=3,label = nothing)
    plot!([-1;2].+i,[2;2].+j,seriescolor=:black,lw=3,label = nothing)
    plot!([2;2].+i,[2;-1].+j,seriescolor=:black,lw=3,label = nothing)
    plot!([2;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    if (count < 2)
        scatter!([0.5;0.5].+i,[0.5;0.5].+j,markersize=12.5,markercolor=:darkred,markershape=:rect,label=nothing)
    elseif (count < 4)
        scatter!([0.5;0.5].+i,[0.5;0.5].+j,markersize=12.5,markercolor=:darkorange1,markershape=:rect,label=nothing)
    elseif (count < 6)
        scatter!([0.5;0.5].+i,[0.5;0.5].+j,markersize=12.5,markercolor=:aquamarine4,markershape=:rect,label=nothing)
    else
        scatter!([0.5;0.5].+i,[0.5;0.5].+j,markersize=12.5,markercolor=:royalblue1,markershape=:rect,label=nothing)
    end
    global count = count+1
end

for j = 0:2:2
    for i = 0:2:6
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    end
end

for j = 0:2:2
    for i = (0:2:6) .+ 12
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    end
end

for j = 0:2:2
    for i = (0:2:6) .+ 24
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    end
end

for j = 0:2:2
    for i = (0:2:6) .+ 36
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    end
end

############
#  Step 1
############

count = 0
for j = (2:-2:0) .- 12
    for i = 0:2:6
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        if (count < 2)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:darkred,markershape=:rect,label=nothing)
        elseif (count < 4)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:darkorange1,markershape=:rect,label=nothing)
        end
        global count = count + 1
    end
end

for j = (0:2:2) .- 12
    for i = (0:2:6) .+ 12
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    end
end

count = 0
for j = (2:-2:0) .- 12
    for i = (0:2:6) .+ 24
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        if (count < 2)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:aquamarine4,markershape=:rect,label=nothing)
        elseif (count < 4)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:royalblue1,markershape=:rect,label=nothing)
        end
        global count = count + 1
    end
end

for j = (0:2:2) .- 12
    for i = (0:2:6) .+ 36
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
    end
end

############
#  Step 2
############

count = 0
for j = (2:-2:0) .- 12*2
    for i = 0:2:6
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        if (count < 2)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:darkred,markershape=:rect,label=nothing)
        end
        global count = count + 1
    end
end

count = 0
for j = (2:-2:0) .- 12*2
    for i = (0:2:6) .+ 12
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        if (count < 2)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:darkorange1,markershape=:rect,label=nothing)
        end
        global count = count + 1
    end
end

count = 0
for j = (2:-2:0) .- 12*2
    for i = (0:2:6) .+ 24
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        if (count < 2)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:aquamarine4,markershape=:rect,label=nothing)
        end
        global count = count + 1
    end
end

count = 0
for j = (2:-2:0) .- 12*2
    for i = (0:2:6) .+ 36
        plot!([-1;-1].+i,[-1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([-1;1].+i,[1;1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;1].+i,[1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        plot!([1;-1].+i,[-1;-1].+j,seriescolor=:black,lw=3,label = nothing)
        if (count < 2)
            scatter!([0;0].+i,[0;0].+j,markersize=8,markercolor=:royalblue1,markershape=:rect,label=nothing)
        end
        global count = count + 1
    end
end

plot!([0;0],[0;0],seriescolor=:black,lw=3,label = nothing)
savefig("IO.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
import Compat.ASCIIString

function logsystem()
    logsys = LogSystem()

    register_log!(logsys, "nec_cards", ["card"], [ASCIIString], "cards",
        x->begin
            code = x[1]
            txt = cardtext(Val{Symbol(code)}, x)
            [txt]
        end)
    logsys
end

function cardtext(::Type{Val{:CM}}, args)
    #"CM", comment = args
    cardtext(args)
end

function cardtext(::Type{Val{:CE}}, args)
    #"CE", comment = args
    cardtext(args)
end

function cardtext(::Type{Val{:GW}}, args)
    #"GW", tag_id, segment_count, xw1, yw1, zw1, xw2, yw2, zw2, rad, radel, rrad = args
    rad = args[10]
    if rad == 0.0
        error("Tapering not implemented")
    end
    join(map(string, args[1:10]), "\t")
end

function cardtext(::Type{Val{:GE}}, args)
    #"GE", gpflag = args
    cardtext(args)
end

function cardtext(::Type{Val{:GN}}, args)
    #"GN", iperf, nradl, 0, 0, epse, sig, tmp3, tmp4, tmp5, tmp6 = args
    cardtext(args)
end

function cardtext(::Type{Val{:FR}}, args)
    #"FR", in_ifrq, in_nfrq, 0, 0, in_freq_mhz, in_del_freq = args
    cardtext(args)
end

function cardtext(::Type{Val{:EX}}, args)
    #"EX", extype, i2, i3, i4, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6
    cardtext(args)
end

function cardtext(::Type{Val{:RP}}, args)
    #"RP", calc_mode, n_theta, n_phi, output_format, normalization,
        #D, A, theta0, phi0, delta_theta, delta_phi, radial_distance,
        #gain_norm] = args
    cardtext(args)
end

cardtext(::Any, args) = cardtext(args)
cardtext(args) = join(map(string, args), "\t")


text:
"""
    bmm(A, B ; transA=false, transB=false)
Perform a batch matrix-matrix product of matrices stored in `A` and `B`. size(A,2) ==
size(B,1) and size(A)[3:end] and size(B)[3:end] must match.
If A is a (m,n,b...) tensor, B is a (n,k,b...) tensor, and the output is a (m,k,b...)
tensor.
"""
function bmm(A::AbstractArray{T}, B::AbstractArray{T}; transA::Bool = false, transB::Bool = false) where T
    sa, sb = size(A), size(B)
    m, k   = transA ? (sa[2],sa[1]) : (sa[1],sa[2])
    kb, n  = transB ? (sb[2],sb[1]) : (sb[1],sb[2])
    @assert kb == k && sa[3:end]==sb[3:end]
    a3, b3 = reshape(A,sa[1],sa[2],:), reshape(B,sb[1],sb[2],:)
    C = similar(A, m, n, size(a3,3))
    bmm!((transA ? 'T' : 'N'), (transB ? 'T' : 'N'), one(T), a3, b3, zero(T), C)
    reshape(C,m,n,sb[3:end]...)
end

function bmm(A::KnetArray{T}, B::KnetArray{T}; transA::Bool = false, transB::Bool = false) where {T}
    sa, sb = size(A), size(B)
    m, k   = transA ? (sa[2],sa[1]) : (sa[1],sa[2])
    kb, n  = transB ? (sb[2],sb[1]) : (sb[1],sb[2])
    @assert kb == k && sa[3:end]==sb[3:end]
    a3, b3 = reshape(A,sa[1],sa[2],:), reshape(B,sb[1],sb[2],:)
    C = similar(A, m, n, size(a3,3))
    bmm!((transA ? 'T' : 'N'), (transB ? 'T' : 'N'), one(T), a3, b3, zero(T), C)
    reshape(C,m,n,sb[3:end]...)
end

function bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T}) where {T}
    cublasop(c::Char)=(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    if ndims(A) != 3 || ndims(B) != 3
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    ma,ka,bsa = size(A)
    kb,nb,bsb = size(B)
    if bsa != bsb
        throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    bs = bsa
    if transA == 'N'
        m=ma; k=ka;
    else
        m=ka; k=ma;
    end
    if transB == 'N'
        k == kb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=nb; k==kb;
    else
        k == nb || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
        n=kb; k==nb;
    end
    (m == size(C,1) && n == size(C,2) && bs == size(C,3)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    lda,ldb,ldc= ma,kb,m
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T[alpha]; beta = T[beta]
    strideA, strideB, strideC = m*k, k*n, m*n
    if T<:Float64
        @cublas(cublasDgemmStridedBatched, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Clonglong, Ptr{T}, Cint, Clonglong, Ptr{T}, Ptr{T}, Cint, Clonglong, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    elseif T<:Float32
        @cublas(cublasSgemmStridedBatched, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Clonglong, Ptr{T}, Cint, Clonglong, Ptr{T}, Ptr{T}, Cint, Clonglong, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, strideA, B, ldb, strideB, beta, C, ldc, strideC, bs)
    else
        error("CUBLAS does not support $T")
    end
    return C
end

using Base: has_offset_axes, unsafe_convert
using LinearAlgebra: chkstride1, BlasInt
using LinearAlgebra.BLAS: libblas, @blasfunc

for (gemm, elty) in
    ((:dgemm_,:Float64),
     (:sgemm_,:Float32),)
    @eval begin
        function bmm!(transA::AbstractChar,
                               transB::AbstractChar,
                               alpha::($elty),
                               A::AbstractArray{$elty, 3},
                               B::AbstractArray{$elty, 3},
                               beta::($elty),
                               C::AbstractArray{$elty, 3})
            @assert !has_offset_axes(A, B, C)
            @assert size(A, 3) == size(B, 3) == size(C, 3) "batch size mismatch"
            m = size(A, transA == 'N' ? 1 : 2)
            ka = size(A, transA == 'N' ? 2 : 1)
            kb = size(B, transB == 'N' ? 1 : 2)
            n = size(B, transB == 'N' ? 2 : 1)
            if ka != kb || m != size(C,1) || n != size(C,2)
                throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))"))
            end
            chkstride1(A)
            chkstride1(B)
            chkstride1(C)

            ptrA = unsafe_convert(Ptr{$elty}, A)
            ptrB = unsafe_convert(Ptr{$elty}, B)
            ptrC = unsafe_convert(Ptr{$elty}, C)

            for k in 1:size(A, 3)
                ccall((@blasfunc($gemm), libblas), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                     Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt},
                     Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty},
                     Ref{BlasInt}),
                     transA, transB, m, n,
                     ka, alpha, ptrA, max(1,stride(A,2)),
                     ptrB, max(1,stride(B,2)), beta, ptrC,
                     max(1,stride(C,2)))

                ptrA += size(A, 1) * size(A, 2) * sizeof($elty)
                ptrB += size(B, 1) * size(B, 2) * sizeof($elty)
                ptrC += size(C, 1) * size(C, 2) * sizeof($elty)
            end

            C
        end
    end
end

@primitive bmm(x1,x2; transA::Bool=false, transB::Bool=false),dy,y (transA ? bmm(x2, dy; transA=transB , transB=true) :  bmm(dy, x2;  transA=false, transB=!transB) )    (transB ? bmm(dy,x1; transA=true , transB=transA) :  bmm(x1, dy;  transA=!transA , transB=false))
@zerograd  bmm!(transA::AbstractChar, transB::AbstractChar, alpha::Number, A, B, beta::Number, C)


text:
using Kinetic, Plots

# initialization
set, ctr, xface, yface, t = initialize("config.toml")

# solution algorithm
t = solve!(set, ctr, xface, yface, t)

# visualization
plot(set, ctr)


text:
module TestUtils

const __ATOL = 1e-9
const __RTOL = 1e-9

using Distances
using LinearAlgebra
using KernelFunctions
using Random
using Test

"""
    test_interface(
        k::Kernel,
        x0::AbstractVector,
        x1::AbstractVector,
        x2::AbstractVector;
        atol=__ATOL,
    )

Run various consistency checks on `k` at the inputs `x0`, `x1`, and `x2`.
`x0` and `x1` should be of the same length with different values, while `x0` and `x2` should
be of different lengths.

    test_interface([rng::AbstractRNG], k::Kernel, T::Type{<:AbstractVector}; atol=__ATOL)

`test_interface` offers certain types of test data generation to make running these tests
require less code for common input types. For example, `Vector{<:Real}`, `ColVecs{<:Real}`,
and `RowVecs{<:Real}` are supported. For other input vector types, please provide the data
manually.
"""
function test_interface(
    k::Kernel,
    x0::AbstractVector,
    x1::AbstractVector,
    x2::AbstractVector;
    atol=__ATOL,
    rtol=__RTOL,
)
    # Ensure that we have the required inputs.
    @assert length(x0) == length(x1)
    @assert length(x0) ≠ length(x2)

    # Check that kernelmatrix_diag basically works.
    @test kernelmatrix_diag(k, x0, x1) isa AbstractVector
    @test length(kernelmatrix_diag(k, x0, x1)) == length(x0)

    # Check that pairwise basically works.
    @test kernelmatrix(k, x0, x2) isa AbstractMatrix
    @test size(kernelmatrix(k, x0, x2)) == (length(x0), length(x2))

    # Check that elementwise is consistent with pairwise.
    @test kernelmatrix_diag(k, x0, x1) ≈ diag(kernelmatrix(k, x0, x1)) atol = atol

    # Check additional binary elementwise properties for kernels.
    @test kernelmatrix_diag(k, x0, x1) ≈ kernelmatrix_diag(k, x1, x0)
    @test kernelmatrix(k, x0, x2) ≈ kernelmatrix(k, x2, x0)' atol = atol

    # Check that unary elementwise basically works.
    @test kernelmatrix_diag(k, x0) isa AbstractVector
    @test length(kernelmatrix_diag(k, x0)) == length(x0)

    # Check that unary pairwise basically works.
    @test kernelmatrix(k, x0) isa AbstractMatrix
    @test size(kernelmatrix(k, x0)) == (length(x0), length(x0))
    @test kernelmatrix(k, x0) ≈ kernelmatrix(k, x0)' atol = atol

    # Check that unary elementwise is consistent with unary pairwise.
    @test kernelmatrix_diag(k, x0) ≈ diag(kernelmatrix(k, x0)) atol = atol

    # Check that unary pairwise produces a positive definite matrix (approximately).
    @test eigmin(Matrix(kernelmatrix(k, x0))) > -atol

    # Check that unary elementwise / pairwise are consistent with the binary versions.
    @test kernelmatrix_diag(k, x0) ≈ kernelmatrix_diag(k, x0, x0) atol = atol rtol = rtol
    @test kernelmatrix(k, x0) ≈ kernelmatrix(k, x0, x0) atol = atol rtol = rtol

    # Check that basic kernel evaluation succeeds and is consistent with `kernelmatrix`.
    @test k(first(x0), first(x1)) isa Real
    @test kernelmatrix(k, x0, x2) ≈ [k(xl, xr) for xl in x0, xr in x2]

    tmp = Matrix{Float64}(undef, length(x0), length(x2))
    @test kernelmatrix!(tmp, k, x0, x2) ≈ kernelmatrix(k, x0, x2)

    tmp_square = Matrix{Float64}(undef, length(x0), length(x0))
    @test kernelmatrix!(tmp_square, k, x0) ≈ kernelmatrix(k, x0)

    tmp_diag = Vector{Float64}(undef, length(x0))
    @test kernelmatrix_diag!(tmp_diag, k, x0) ≈ kernelmatrix_diag(k, x0)
    @test kernelmatrix_diag!(tmp_diag, k, x0, x1) ≈ kernelmatrix_diag(k, x0, x1)
end

function test_interface(
    rng::AbstractRNG, k::Kernel, ::Type{Vector{T}}; kwargs...
) where {T<:Real}
    return test_interface(
        k, randn(rng, T, 1001), randn(rng, T, 1001), randn(rng, T, 1000); kwargs...
    )
end

function test_interface(
    rng::AbstractRNG, k::MOKernel, ::Type{Vector{Tuple{T,Int}}}; dim_out=1, kwargs...
) where {T<:Real}
    return test_interface(
        k,
        [(randn(rng, T), rand(rng, 1:dim_out)) for i in 1:51],
        [(randn(rng, T), rand(rng, 1:dim_out)) for i in 1:51],
        [(randn(rng, T), rand(rng, 1:dim_out)) for i in 1:50];
        kwargs...,
    )
end

function test_interface(
    rng::AbstractRNG, k::Kernel, ::Type{<:ColVecs{T}}; dim_in=2, kwargs...
) where {T<:Real}
    return test_interface(
        k,
        ColVecs(randn(rng, T, dim_in, 1001)),
        ColVecs(randn(rng, T, dim_in, 1001)),
        ColVecs(randn(rng, T, dim_in, 1000));
        kwargs...,
    )
end

function test_interface(
    rng::AbstractRNG, k::Kernel, ::Type{<:RowVecs{T}}; dim_in=2, kwargs...
) where {T<:Real}
    return test_interface(
        k,
        RowVecs(randn(rng, T, 1001, dim_in)),
        RowVecs(randn(rng, T, 1001, dim_in)),
        RowVecs(randn(rng, T, 1000, dim_in));
        kwargs...,
    )
end

function test_interface(
    rng::AbstractRNG, k::Kernel, ::Type{<:Vector{Vector{T}}}; dim_in=2, kwargs...
) where {T<:Real}
    return test_interface(
        k,
        [randn(rng, T, dim_in) for _ in 1:1001],
        [randn(rng, T, dim_in) for _ in 1:1001],
        [randn(rng, T, dim_in) for _ in 1:1000];
        kwargs...,
    )
end

function test_interface(k::Kernel, T::Type{<:AbstractVector}; kwargs...)
    return test_interface(Random.GLOBAL_RNG, k, T; kwargs...)
end

function test_interface(rng::AbstractRNG, k::Kernel, T::Type{<:Real}; kwargs...)
    @testset "Vector{$T}" begin
        test_interface(rng, k, Vector{T}; kwargs...)
    end
    @testset "ColVecs{$T}" begin
        test_interface(rng, k, ColVecs{T}; kwargs...)
    end
    @testset "RowVecs{$T}" begin
        test_interface(rng, k, RowVecs{T}; kwargs...)
    end
    @testset "Vector{Vector{T}}" begin
        test_interface(rng, k, Vector{Vector{T}}; kwargs...)
    end
end

function test_interface(k::Kernel, T::Type{<:Real}=Float64; kwargs...)
    return test_interface(Random.GLOBAL_RNG, k, T; kwargs...)
end

end # module


text:
# This file is a part of UpROOT.jl, licensed under the MIT License (MIT).

"""
    TDirectory

`UpROOT.TDirectory` is a wrapper around Python objects of type
`uproot.rootio.ROOTDirectory`.

`TDirectory` behaves similar to a Julia `Dict`, it supports the functions
`keys`, `length` and `getindex`.

Limitations: Write access is not implemented yet.

[`TFile`](@ref) is defined as an alias for TDirectory in UpROOT.jl.
"""
struct TDirectory
    pyobj::PyObject

    function TDirectory(x::PyObject)
        if !pybuiltin(:isinstance)(x, uproot.rootio.ROOTDirectory)
            throw(ArgumentError("Python object $x is not a ROOT TDirectory"))
        end
        new(x)
    end
end

export TDirectory

Base.keys(tdir::TDirectory) = pyobj(tdir).keys()

Base.length(tdir::TDirectory) = length(keys(tdir))

Base.getindex(tdir::TDirectory, objname::AbstractString) = py2jl(get(pyobj(tdir), PyObject, objname))

"""
    TFile = TDirectory

The uproot Python package doesn't use separate types for the ROOT classes
`TFile` and `TDirectoy` - so in UpROOT.jl, so in UpROOT.jl, `TFile` is just
defined as an alias for [`TDirectory`](@ref).

Constructors:

    * `TFile(filename::AbstractString)::TFile`

Use `TFile(filename)` to open files instead of `TDirectory(filename)` for
increased clarity in your code.
"""
const TFile = TDirectory

export TFile

TFile(filename::AbstractString) = TDirectory(uproot.open(filename))


text:
"""
Computes the maximum multiroute flow (for any real number of routes)
between the source and target vertexes in a flow graph using the
[Extended Multiroute Flow algorithm](http://dx.doi.org/10.1016/j.disopt.2016.05.002).
If a number of routes is given, returns the value of the multiroute flow as
well as the final flow matrix, along with a multiroute cut if
Boykov-Kolmogorov max-flow algorithm is used as a subroutine.
Otherwise, it returns the vector of breaking points of the parametric
multiroute flow function.
Use a default capacity of 1 when the capacity matrix isn\'t specified.
Requires arguments:
- flow_graph::ADiGraph                    # the input graph
- source::Int                            # the source vertex
- target::Int                            # the target vertex
- capacity_matrix::AbstractMatrix{T}   # edge flow capacities
- flow_algorithm::AbstractFlowAlgorithm  # keyword argument for algorithm
- routes::Int                            # keyword argument for routes
"""
function emrf(
  flow_graph::ADiGraph,                   # the input graph
  source::Int,                           # the source vertex
  target::Int,                           # the target vertex
  capacity_matrix::AbstractMatrix{T},  # edge flow capacities
  flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
  routes::R = 0
  ) where {T<:AbstractFloat, R<:Real}
  breakingpoints = breakingPoints(flow_graph, source, target, capacity_matrix)
  if routes > zero(R)
    x, f = intersection(breakingpoints, routes)
    return maximum_flow(flow_graph, source, target, capacity_matrix,
                        algorithm = flow_algorithm, restriction = x)
  end
  return breakingpoints
end

"""
Output a set of (point,slope) that compose the restricted max-flow function.
One point by possible slope is enough (hence O(λ×max_flow) complexity).
Requires arguments:
- flow_graph::ADiGraph                    # the input graph
- source::Int                            # the source vertex
- target::Int                            # the target vertex
- capacity_matrix::AbstractMatrix{T}   # edge flow capacities
"""
function auxiliaryPoints(
  flow_graph::ADiGraph,                   # the input graph
  source::Int,                           # the source vertex
  target::Int,                           # the target vertex
  capacity_matrix::AbstractMatrix{T}   # edge flow capacities
  ) where T<:AbstractFloat
  # Problem descriptors
  λ = maximum_flow(flow_graph, source, target)[1] # Connectivity
  n = nv(flow_graph) # number of nodes
  r1, r2 = minmaxCapacity(capacity_matrix) # restriction left (1) and right (2)
  auxpoints = fill((0., 0.), λ + 1)

  # Initialisation of left side (1)
  f1, F1, cut1 = maximum_flow(flow_graph, source, target, capacity_matrix,
                 algorithm = BoykovKolmogorovAlgorithm(), restriction = r1)
  s1 = slope(flow_graph, capacity_matrix, cut1, r1) # left slope
  auxpoints[λ + 1 - s1] = (r1, f1) # Add left initial auxiliary point

  # Initialisation of right side (2)
  f2, F2, cut2 = maximum_flow(flow_graph, source, target, capacity_matrix,
                 algorithm = BoykovKolmogorovAlgorithm(), restriction = r2)
  s2 = slope(flow_graph, capacity_matrix, cut2, r2) # right slope
  auxpoints[λ + 1 - s2] = (r2, f2) # Add right initial auxiliary point

  # Loop if the slopes are distinct by at least 2
  if s1 > s2 + 1
    queue = [((f1, s1, r1), (f2, s2, r2))]

    while !isempty(queue)
      # Computes an intersection (middle) with a new associated slope
      (f1, s1, r1), (f2, s2, r2) = pop!(queue)
      r, expectedflow = intersection(r1, f1, s1, r2, f2, s2)
      f, F, cut = maximum_flow(flow_graph, source, target, capacity_matrix,
                  algorithm = BoykovKolmogorovAlgorithm(), restriction = r)
      s = slope(flow_graph, capacity_matrix, max.(cut, 1), r) # current slope
      auxpoints[λ + 1 - s] = (r, f)
      # If the flow at the intersection (middle) is as expected
      if expectedflow ≉ f # approximatively not equal (enforced by floating precision)
        # if the slope difference between (middle) and left is at least 2
        # push (left),(middle)
        if s1 > s + 1 && (r2, f2) ≉ (r, f)
          q = (f1, s1, r1), (f, s, r)
          push!(queue, q)
        end
        # if the slope difference between (middle) and right is at least 2
        # push (middle),(right)
        if s > s2 + 1 && (r1, f1) ≉ (r, f)
          q = (f, s, r), (f2, s2, r2)
          push!(queue, q)
        end
      end
    end
  end
  return auxpoints
end

"""
Calculates the breaking of the restricted max-flow from a set of auxiliary points.
Requires arguments:
- flow_graph::ADiGraph                    # the input graph
- source::Int                            # the source vertex
- target::Int                            # the target vertex
- capacity_matrix::AbstractMatrix{T}   # edge flow capacities
"""
function breakingPoints(
  flow_graph::ADiGraph,                   # the input graph
  source::Int,                           # the source vertex
  target::Int,                           # the target vertex
  capacity_matrix::AbstractMatrix{T}   # edge flow capacities
  ) where T<:AbstractFloat
  auxpoints = auxiliaryPoints(flow_graph, source, target, capacity_matrix)
  λ = length(auxpoints) - 1
  left_index = 1
  breakingpoints = Vector{Tuple{T, T, Int}}()

  for (id, point) in enumerate(auxpoints)
    if id == 1
      push!(breakingpoints, (0., 0., λ))
    else
      pleft = breakingpoints[left_index]
      if point[1] != 0
        x, y = intersection(pleft[1], pleft[2], pleft[3],
                            point[1], point[2], λ + 1 - id)
        push!(breakingpoints,(x, y, λ + 1 - id))
        left_index += 1
      end
    end
  end
  return breakingpoints
end

"""
Function to get the nonzero min and max function of a Matrix.

Note: this is more efficient than maximum() / minimum() / extrema()
since we have to ignore zero values.since we have to ignore zero values.

Requires argument:
- capacity_matrix::AbstractMatrix{T}   # edge flow capacities
"""
function minmaxCapacity(
  capacity_matrix::AbstractMatrix{T}    # edge flow capacities
  ) where T<:AbstractFloat
# Function to get the nonzero min and max function of a Matrix
# note: this is more efficient than maximum() / minimum() / extrema()
#       since we have to ignore zero values.
  cmin, cmax = typemax(T), typemin(T)
  for c in capacity_matrix
    if c > zero(T)
      cmin = min(cmin, c)
    end
    cmax = max(cmax, c)
  end
  return cmin, cmax
end

"""
Function to get the slope of the restricted flow. The slope is initialized at 0
and is incremented for each non saturated edge in the restricted min-cut.
Requires argument:
  flow_graph::ADiGraph,                   # the input graph
  capacity_matrix::AbstractMatrix{T},  # edge flow capacities
  cut::Vector{Int},                      # cut information for vertices
  restriction::T                         # value of the restriction
"""
function slope(
  flow_graph::ADiGraph,                   # the input graph
  capacity_matrix::AbstractMatrix{T},  # edge flow capacities
  cut::Vector{Int},                      # cut information for vertices
  restriction::T                         # value of the restriction
  ) where T<:AbstractFloat
  # Function to get the slope of the restricted flow
  slope = 0
  for e in edges(flow_graph)
      if cut[dst(e)] == 2 > cut[src(e)] &&
            capacity_matrix[src(e), dst(e)] > restriction
        slope += 1
    end
  end
  return slope
end

"""
Computes the intersection between:
1) Two lines
2) A set of segments and a linear function of slope k passing by the origin.
Requires argument:
1) - x1, y1, a1, x2, y2, a2::T<:AbstractFloat # Coordinates/slopes
2) - points::Vector{Tuple{T, T, Int}}         # vector of points with T<:AbstractFloat
   - k::R<:Real                             # number of routes (slope of the line)
"""
function intersection(
  x1::T,          # x coordinate of point 1
  y1::T,          # y coordinate of point 1
  a1::Int,        # slope passing by point 1
  x2::T,          # x coordinate of point 2
  y2::T,          # y coordinate of point 2
  a2::R           # slope passing by point 2
  ) where {T<:AbstractFloat, R<:Real}

  (a1 == a2) && return -1., -1. # result will be ignored in other intersection method
  b1 = y1 - a1 * x1
  b2 = y2 - a2 * x2
  x = (b2 - b1) / (a1 - a2)
  y = a1 * x + b1
  return x, y
end

# Compute the intersection between a set of segment and a line of slope k passing by the origin
function intersection(
  points::Vector{Tuple{T, T, Int}},  # vector of breaking points
  k::R                               # number of routes (slope of the line)
  ) where {T<:AbstractFloat, R<:Real}
  λ = points[1][1] # Connectivity

  # Loop over the segments (pair of breaking points)
  for (id, p) in enumerate(points[1:(end - 1)])
    if id == 1
      (k ≈ λ) && return points[2]
    else
      x, y = intersection(p[1], p[2], p[3], 0., 0., k)
      (p[1] ≤ x ≤ points[id + 1][1]) && return x, y
    end
  end
  p = points[end]
  return intersection(p[1], p[2], p[3], 0., 0., k)
end

"""
Redefinition of ≈ (isapprox) for a pair of points
Requires argument:
  a::Tuple{T, T},         # Point A with floating coordinates
  b::Tuple{T, T}          # Point B with floating coordinates
"""
function ≈(
  a::Tuple{T, T},         # Point A with floating coordinates
  b::Tuple{T, T}          # Point B with floating coordinates
  ) where T<:AbstractFloat
  return a[1] ≈ b[1] && a[2] ≈ b[2]
end


text:
# This file was generated, do not modify it. # hide
@load KNNRegressor
knn_model = KNNRegressor(K=10)