text:
# ==================================== bernoulli_numbers(nmax) =================

@doc raw"""
    bernoulli_numbers(nmax [, T=Int])

Bernoulli numbers ``B_0,\ \cdots,\ B_{nmax}`` calculated by repetative use of the recurrence relation
```math
    B_n = - \frac{1}{n+1}\sum_{k=0}^{n-1}\frac{(n+1)!}{k!(n+1-k)}B_k.
```
Special numbers: ``B_0=1,\ B_1=-1/2,\ B_{2n+1}=0\ (\rm{for}\ n>1)``.
### Examples:
```
bernoulli_numbers(10)
11-element Vector{Rational{Int64}}:
  1//1
 -1//2
  1//6
  0//1
 -1//30
  0//1
  1//42
  0//1
 -1//30
  0//1
  5//66
```
"""
function bernoulli_numbers(nmax::Int; T=Int)

    B = Base.ones(Rational{T},nmax+1)

    for m = 2:nmax+1
        B[m] = m > 2 ? 0//1 : -1
        if Base.isodd(m)
            b = 1
            for j = 1:m-1
                B[m] -= B[j] * b
                b *= m+1-j
                b = b÷j      # binomial coefficients are integers
            end
        end
        B[m] = B[m] // m
    end

    return B

end
function bernoulli_numbers(nmax::Int)       # short argument: better performance

    B = Base.ones(Rational{Int},nmax+1)

    for m = 2:nmax+1
        B[m] = m > 2 ? 0//1 : -1
        if Base.isodd(m)
            b = 1
            for j = 1:m-1
                B[m] -= B[j] * b
                b *= m+1-j
                b = b÷j      # binomial coefficients are integers
            end
        end
        B[m] = B[m] // m
    end

    return B

end

# ==================================== faulhaber_polynom(p) ====================

@doc raw"""
    faulhaber_polynom(p [, T=Int])

Vector representation of the Faulhaber polynomial of degree ``p``,
```math
    F(n,p)=\frac{1}{p}\sum_{j=1}^{p}{\binom {p}{p-j}}B_{p-j}n^{j}.
```
``F(n,p)=`` `polynom(c,n)`, where ``c=[c_0,\cdots,\ c_p]`` is the coefficient vector, with
```math
    c_0=0,\ c_j=\frac{1}{p}{\binom {p}{p-j}}B_{p-j},
```
with ``j∈\{ 1,\cdots,\ p\}``. The ``B_0,\cdots,\ B_{p-1}`` are Bernoulli numbers
(but with ``B_1=+\frac{1}{2}`` rather than ``-\frac{1}{2}``).
### Example:
```
faulhaber_polynom(6)
7-element Vector{Rational{Int64}}:
  0//1
  0//1
 -1//12
  0//1
  5//12
  1//2
  1//6
```
"""
function faulhaber_polynom(k::Int; T=Int)

    k < 1 && return 0
    k > 1 || return 1//1

    P = CamiXon.pascal_triangle(k; T)[end][1:end-1]
    B = CamiXon.bernoulli_numbers(k-1; T); B[2]=-B[2]

    F = (B .* P)  // k

    F = Base.append!(F,0//1)   # add polynomial constant (zero in this case)

    return Base.reverse(F)     # reverse to standard order

end
function faulhaber_polynom(k::Int)       # short argument: better performance

    k < 1 && return 0
    k > 1 || return 1//1

    P = CamiXon.pascal_triangle(k)[end][1:end-1]
    B = CamiXon.bernoulli_numbers(k-1); B[2]=-B[2]

    F = (B .* P)  // k

    F = Base.append!(F,0//1)   # add polynomial constant (zero in this case)

    return Base.reverse(F)     # reverse to standard order

end

# =================================== faulhaber_summation(n,p;T) ===============

@doc raw"""
    faulhaber_summation(n, p [, T=Int])

Sum of powers of natural numbers ``1,\ \cdots,\ n``,
```math
    FS(n,p)=\sum_{k=1}^{n}k^{p}=F(n,p+1).
```
where ``F(n,p)`` is the Faulhamer polynomial of degree ``p``.
### Examples:
```
faulhaber_summation(5,1)
 15

faulhaber_summation(3,60; T=BigInt)
  42391158276369125018901280178
```
"""
function faulhaber_summation(n::Int, p::Int; T=Int)

    n ≠ 0 || return nothing

    F = CamiXon.faulhaber_polynom(p+1; T)
    o = 0
    for k=1:p+1
        for i=1:k
            F[k+1] *= n # avoid n^k in o = Base.sum([F[k+1]*n^k for k=1:p+1])
        end
        o += F[k+1]
    end

    Base.denominator(o) == 1 || error("jwError: Faulhaber sum failed")

    return Base.numerator(o)

end
function faulhaber_summation(n::Int, p::Int)   # short argument: better performance

    n ≠ 0 || return 0

    F = CamiXon.faulhaber_polynom(p+1)
    o = 0
    for k=1:p+1
        for i=1:k
            F[k+1] *= n # avoid n^k in o = Base.sum([F[k+1]*n^k for k=1:p+1])
        end
        o += F[k+1]
    end

    Base.denominator(o) == 1 || error("jwError: Faulhaber sum failed")

    return Base.numerator(o)

end

# =================================== harmonic number(n;T) ===============

@doc raw"""
    harmonic_number(n [, T=Int])

Sum of the reciprocals of the first ``n`` natural numbers
```math
    H_n=\sum_{k=1}^{n}\frac{1}{k}.
```
### Examples:
```
o = [harmonic_number(i) for i=1:10]; println(o)
 [1//1, 3//2, 11//6, 25//12, 137//60, 49//20, 363//140, 761//280, 7129//2520, 7381//2520]

harmonic_number(60; T=BigInt)
 15117092380124150817026911//3230237388259077233637600

harmonic_number(12) == harmonic_number(12, 1)
 true
```
"""
function harmonic_number(n::Int; T=Int)

    n ≠ 0 || return 0

    o::Base.Rational{T} = 0//1
    for j=1:n
        o += 1//j
    end

    return o

end
function harmonic_number(n::Int)            # short argument: better performance

    n ≠ 0 || return 0

    o = 0//1
    for j=1:n
        o += 1//j
    end

    return o

end

# =================================== harmonic number(n, p;T) ===============

@doc raw"""
    harmonic_number(n, p [, T=Int])

Sum of the ``p_{th}`` power of reciprocals of the first ``n`` numbers
```math
    H_{n,p}=\sum_{k=1}^{n}\frac{1}{k^p}.
```
### Examples:
```
harmonic_number(12, 3)
 25535765062457//21300003648000

harmonic_number(12, 5; T=BigInt)
 16971114472329088045481//16366888723117363200000

harmonic_number(12, -3) == faulhaber_summation(12, 3)
  true
```
"""
function harmonic_number(n::Int, p::Int; T=Int)

    n ≠ 0 || return 0

    if p > 0
        o::Base.Rational{T} = 0//1
        for j=1:n
            a = 1
            for i=1:p
                a *= j
            end
        o += 1//a
        end
    else
        p = -p
        F = CamiXon.faulhaber_polynom(p+1; T)
        o = 0
        for k=1:p+1
            for i=1:k
                F[k+1] *= n
            end
            o += F[k+1]
        end
        Base.denominator(o) == 1 || error("jwError: Faulhaber sum failed")
        o = Base.numerator(o)
    end

    return o

end
function harmonic_number(n::Int, p::Int)    # short argument: better performance

    n ≠ 0 || return 0

    if p > 0
        o = 0//1
        for j=1:n
            a = 1
            for i=1:p
                a *= j
            end
        o += 1//a
        end
    else
        p = -p
        F = CamiXon.faulhaber_polynom(p+1)
        o = 0
        for k=1:p+1
            for i=1:k
                F[k+1] *= n
            end
            o += F[k+1]
        end
        Base.denominator(o) == 1 || error("jwError: Faulhaber sum failed")
        o = Base.numerator(o)
    end

    return o

end

# ==================================== _canonical_partition(n, m) =======================

function _canonical_partition(n::Int, m::Int)

    o = Base.fill(m,Base.cld(n,m))                              # init partition
    o[Base.cld(n,m)]=((n%m)≠0 ? n%m : m)                        # adjust last element of partition

    return o

end

"""
    canonical_partitions(n; header=false, reverse=true)

The canonical partition in integers of the integer n

header=true : unit patition included in output
#### Examples:
```
canonical_partitions(6; header=true, reverse=false)
6-element Array{Array{Int64,1},1}:
 [6]
 [5, 1]
 [4, 2]
 [3, 3]
 [2, 2, 2]
 [1, 1, 1, 1, 1, 1]

canonical_partitions(6; header=true)
6-element Array{Array{Int64,1},1}:
 [1, 1, 1, 1, 1, 1]
 [2, 2, 2]
 [3, 3]
 [4, 2]
 [5, 1]
 [6]

canonical_partitions(6)
5-element Array{Array{Int64,1},1}:
 [1, 1, 1, 1, 1, 1]
 [2, 2, 2]
 [3, 3]
 [4, 2]
 [5, 1]
```
"""
function canonical_partitions(n::Int, m=0; header=true, reverse=true)

    h = header ? n : n-1

    if m == 0
        if reverse
            o = [_canonical_partition(n,m) for m=1:h]
        else
            o = [_canonical_partition(n,m) for m=h:-1:1]
        end
    elseif 0 < m <= n
        o = _canonical_partition(n,m)
    else
        o = nothing
    end

    return o

end

function _partition_count(n::Int,k::Int)

    (n<0)|(k<0)|(k>n) ? 0 : (k==n)|(k==1) ? 1 : _partition_count(n-k,k) + _partition_count(n-1,k-1)

end

function _partition(a::Array{Int,1}, n::Int, i::Int, cp::Array{Array{Array{Int,1},1},1})

    o = a[1:i-1]
    m = a[i]-1                                           # m: partition value
    ni = n - Base.sum(o)                                 # ni: sub-partition index at partition index i

    Base.append!(o,cp[ni][m])                            # complete partition by appending it to a

    return o

end

function _restricted_partitions(o::Array{Int,1}, n::Int, np::Int, cp::Array{Array{Array{Int,1},1},1})

    oo = [o]

    for p=1:np-1
        i = Base.findlast(x -> x > 1, oo[p])
        Base.append!(oo,[_partition(oo[p],n,i,cp)])
    end

    return oo

end

"""
    integer_partitions(n [,m]; transpose=false, count=false)

default              : The integer partitions of n

count=true           : The number of integer partitions of n

transpose=false/true : for m>0 restricted to partitions with maximum part/length m

definitions:

The integer partition of the positive integer n is a nonincreasing sequence of positive integers p1, p2,... pk whose sum is n.

The elements of the sequence are called the parts of the partition.
#### Examples:
```
integer_partitions(7)
15-element Array{Array{Int64,1},1}:
 [1, 1, 1, 1, 1, 1, 1]
 [2, 2, 2, 1]
 [3, 3, 1]
 [4, 3]
 [5, 2]
 [6, 1]
 [7]
 [2, 2, 1, 1, 1]
 [3, 2, 2]
 [4, 2, 1]
 [5, 1, 1]
 [2, 1, 1, 1, 1, 1]
 [3, 2, 1, 1]
 [4, 1, 1, 1]
 [3, 1, 1, 1, 1]

integer_partitions(7; count=true)
15

integer_partitions(7,4; count=true)
3

integer_partitions(7,4)
3-element Array{Array{Int64,1},1}:
 [4, 3]
 [4, 2, 1]
 [4, 1, 1, 1]

integer_partitions(7,4; transpose=true)
3-element Array{Array{Int64,1},1}:
 [2, 2, 2, 1]
 [3, 2, 1, 1]
 [4, 1, 1, 1]
```
"""
function integer_partitions(n::Int, m=0; transpose=false, count=false)

    cp = [canonical_partitions(m) for m=1:n]
    pc = [_partition_count(n,m)  for m=1:n]
    oo = [ones(Int,n)]

    np = m > 0 ? pc[m] : sum(pc)

    if !count

        if m == 0
            o = [_restricted_partitions(cp[n][p],n,pc[p],cp) for p=2:n]
            for p=1:n-1 append!(oo,o[p]) end
        else
            oo = _restricted_partitions(cp[n][m],n,pc[m],cp)
        end

        if transpose
            for p=1:np
                l = length(oo[p])
                s=max(oo[p][1],l)
                mat = zeros(Int,s,s)
                for j=1:l for i=1:oo[p][j] mat[i,j]=1 end end
                oo[p] = [sum(mat[i,:]) for i=1:oo[p][1]]
            end

        end

    end

    return count ? np : oo

end

# ===================================== log10_characteristic_power(x) ==============================================
"""
    log10_characteristic_power(x)

characteristic power-of-10 of the number x
#### Examples:
```
log10_characteristic_power.([3,30,300])
3-element Vector{Int64}:
 0
 1
 2
```
"""
log10_characteristic_power(x) = Base.round(Int,Base.floor(log10(x)))

# ==================================== log10_mantissa(x) ============================================================
"""
    log10_mantissa(x)

log10 mantissa of the number x
#### Examples:
```
log10_mantissa.([3,30,300])
3-element Vector{Float64}:
 0.47712125471966244
 0.4771212547196624
 0.4771212547196626
```
"""
log10_mantissa(x) = Base.log10(x)-Base.floor(Base.log10(x))

# ==================================== pascal_triangle(nmax)  ============

@doc raw"""
    pascal_triangle(nmax [, T=Int])

Pascal triangle of binomial coefficients ``\binom{n}{k}`` for ``n=0,\ 1,\ \cdots,\ nmax``
### Example:
```
pascal_triangle(5)
6-element Vector{Vector{Int64}}:
 [1]
 [1, 1]
 [1, 2, 1]
 [1, 3, 3, 1]
 [1, 4, 6, 4, 1]
 [1, 5, 10, 10, 5, 1]
```
"""
function pascal_triangle(nmax::Int; T=Int)

    nmax < 0 && error("jwError: nmax must be a non-negative integer")

    o = [Base.ones(T,n+1) for n=0:nmax]

    for n=2:nmax
        for k=1:n÷2
            o[n+1][k+1] = o[n][k+1] + o[n][k]
            o[n+1][n+1-k] = o[n+1][k+1]
        end
    end

    return o

end
function pascal_triangle(nmax::Int)

    nmax < 0 && error("jwError: nmax must be a non-negative integer")

    o = [Base.ones(Int,n+1) for n=0:nmax]

    for n=2:nmax
        for k=1:n÷2
            o[n+1][k+1] = o[n][k+1] + o[n][k]
            o[n+1][n+1-k] = o[n+1][k+1]
        end
    end

    return o

end

# ==================================== pascal_next(nmax)  ============

@doc raw"""
    pascal_next(nmax)

Next row of Pascal triangle
### Example:
```
a = [1, 4, 6, 4, 1]
pascal_next(a)
 [1, 5, 10, 10, 5, 1]
```
"""
function pascal_next(a::Vector{Int})

    n = Base.length(a) + 1
    o = Base.ones(Int,n)

    for k=1:n÷2
        o[k+1] = a[k+1] + a[k]
        o[n-k] = o[k+1]
    end

    return o

end

# ==================================== permutations_unique_count(p, i) =======================

@doc raw"""
    permutations_unique_count(p::Array{Array{Int64,1},1}, i::Int)

Number of unique permutations of the subarray ``p[i]``.
#### Example:
```
p = [[1,2,3],[2,3,1,4,3]]
permutations_unique_count(p,2)
 60
```
"""
function permutations_unique_count(p::Array{Array{Int64,1},1}, i::Int)

    o = Base.factorial(Base.length(p[i]))
    d = Base.Dict([(n,Base.count(x->x==n,p[i])) for n ∈ Base.unique(p[i])])

    for j ∈ Base.eachindex(Base.unique(p[i]))
        o = o ÷ Base.factorial(d[Base.unique(p[i])[j]])
    end

    return o

end

# ==================================== polynom(c,x) ============================================================

@doc raw"""
    polynom(coords,x)

Method to evaluate the function ``f(x)=polynom(c,x)``, where
``c=[c_0,\ \ldots,\ c_d]`` is the vector representation of a polynomial of degree ``d``.
```math
    p(c,x)=c_0 + c_1 x + \cdots + c_d x^d.
```
### Examples:
```
coords = ones(6)             # for polynomial of degree 5 with unit coefficients
f(x) = polynom(coords,x)
println([f(1.0),f(2.0)])     # values of polynomial for x = 1.0 and x = 2.0
 [6.0, 63.0]
```
"""
function polynom(coords::Vector{T}, x::T) where T<:Number

    k = Base.length(coords)
    X = Base.ones(T,k)

    for i=2:k
        X[i] = X[i-1] * x
    end

    return LinearAlgebra.dot(coords, X)

end

# ==================================== polynom_derivative(coords) ==============

@doc raw"""
    polynom_derivative(coords)

Vector representation of the first derivative of the polynomial `coords`,
```math
    p'(c,x)=c_1 + 2 c_2 x + \cdots + d c_d x^{d-1},
```
Polynomials of degree ``d`` are represented by a vector in a vector space of dimension ``d+1``.
The polynomial `coords` is specified by the coordinates vector ``c=[c_0,\ \ldots,\ c_d]``
consisting of the polynomial coefficients.
### Examples:
```
coords=[1,1,1,1,1]                 # vector representation of polynomial of degree d=4
polynom_derivative(coords)         # (first) derivative of polynomial `coords`
4-element Vector{Int64}:
 1
 2
 3
 4
```
"""
function polynom_derivative(coords::Vector{<:Number})

    k = Base.length(coords)
    k > 1 || return [0]

    return coords[2:end] .* Base.OneTo(k-1)

end

# =============================== polynom_derivative(coords[,deriv=0]) =========

@doc raw"""
    polynom_derivatives(coords[,deriv=0])

Vector representation of derivatives of the polynomial `coords`.

Polynomials of degree ``d`` are represented by a vector in a vector space of dimension ``d+1``.
The polynomial `coords` is specified by the coordinates vector ``c=[c_0,\ \ldots,\ c_d]``
consisting of the polynomial coefficients.

`deriv`: derivative of choice; `default`: collection of all (nontrivial) derivatives.

### Examples:
```
coords=[1,1,1,1,1]               # vector representation of a polynomial of degree d=4
polynom_derivatives(coords)      # `all' (nontrivial) derivatives of polynomial `coords`
5-element Vector{Vector{Int64}}:
 [1, 2, 3, 4]
 [2, 6, 12]
 [6, 24]
 [24]
 [0]

polynom_derivatives(coords; deriv=2)          # second derivative of polynomial `coords`
3-element Vector{Int64}:
  2
  6
 12
```
"""
function polynom_derivatives(coords::Vector{<:Number}; deriv=0)

    deriv < 0 && error("jwError: negative derivative not defined")

    k = deriv > 0 ? deriv+1 : Base.length(coords)

    coords = CamiXon.polynom_derivative(coords)

    deriv ≠ 1 ? o = [coords] : return coords

    for i=2:k-1
        coords = polynom_derivative(coords)
        Base.push!(o,coords)
    end

    deriv == 0 ? Base.push!(o,[0]) : return o[deriv]

    return o

end

# ==================================== polynom_power(coords, p) ================

@doc raw"""
    polynom_power(coords, p)

Vector representation of the polynomial `coords` raised to the power `p` which
results in a polynomial in a vector space of dimension ``p d + 1``.

Polynomials of degree ``d`` are represented by a vector in a vector space of dimension ``d+1``.
The polynomial `coords` is specified by the coordinates vector ``c=[c_0,\ \ldots,\ c_d]``
consisting of the polynomial coefficients.
### Examples:
```
coords=[1,1,1]             # vector representation of polynomial of degree ``d=2``
polynom_power(coords,2)
5-element Vector{Int64}:
 1
 2
 3
 2
 1
```
"""
function polynom_power(coords::Vector{<:Number}, power::Int)

    power >= 0 || error("jwError: negative powers not allowed")
    power == 2 && return polynom_product(coords, coords)
    power == 1 && return coords
    power == 0 && return [1]

    o = CamiXon.polynom_product(coords, coords)

    for i=1:power-2
        o = CamiXon.polynom_product(o, coords)
    end

    return o

end

# ==================================== polynom_powers(coords, pmax) ============

@doc raw"""
    polynom_powers(coords, pmax)

The polynomial `coords` raised to the powers 1,...,pmax  which
results in a collection of polynomials in vector spaces of dimension ``d+1`` tot ``p d + 1``.

Polynomials of degree ``d`` are represented by a vector in a vector space of dimension ``d+1``.
The polynomial `coords` is specified by the coordinates vector ``c=[c_0,\ \ldots,\ c_d]``
consisting of the polynomial coefficients.
### Examples:
```
coords=[1,1,1]                   # vector representation of polynomial of degree d=2
polynom_powers(coords,3)
3-element Vector{Vector{Int64}}:
 [1, 1, 1]
 [1, 2, 3, 2, 1]
 [1, 3, 6, 7, 6, 3, 1]
```
"""
function polynom_powers(coords::Vector{<:Number}, pmax::Int)

    pmax > 0 || error("jwError: minimum power included is unity")

    o = [coords]

    for i=1:pmax-1
        Base.push!(o,CamiXon.polynom_product(o[end], coords))
    end

    return o

end

# ==================================== polynom_primitive(coords) ====================

@doc raw"""
    polynom_primitive(coords)

Vector representation of the primitive of the polynomial `coords` which is a
polynomial in a vector space of dimension ``p d + 1``.
```math
    P(c,x)=c_{int} +c_0 x + \frac{1}{2} c_1 x^2 + \frac{1}{3} c_2 x^3 + \cdots + \frac{1}{d+1} c_d x^{d+1},
```
The constant of integration is set to zero, ``c_{int} = 0``.

Polynomials of degree ``d`` are represented by a vector in a vector space of dimension ``d+1``.
The polynomial `coords` is specified by the coordinates vector ``c=[c_0,\ \ldots,\ c_d]``
consisting of the polynomial coefficients.
### Examples:
```
coords=[1,1,1,1,1]         # vector representation of polynomial of degree ``d=4``
polynom_primitive(coords)
6-element Vector{Rational{Int64}}:
 0//1
 1//1
 1//2
 1//3
 1//4
 1//5
```
"""
function polynom_primitive(coords::Vector{<:Number})

    d = [1//p for p ∈ Base.eachindex(coords)]

    coords = coords .* d

    return Base.pushfirst!(coords,0)      # constant of integration equal to zero

end

# ==================================== polynom_product(a, b) ============================================================

@doc raw"""
    polynom_product(a::Vector{<:Number}, b::Vector{<:Number})

Vector representation of the product of two polynomials, ``a`` and ``b`` which
is a polynomial in a vector space of dimension ``d=m+n``,
```math
    p(c,x)=a_0b_0 + (a_0b_1 + b_0a_1)x + \cdots + a_n b_m x^{n+m}.
```
Polynomials of degree ``d`` are represented by a vector in a vector space of dimension ``d+1``
The polynomial `coords` is specified by the coordinates vector ``c=[c_0,\ \ldots,\ c_d]``
consisting of the polynomial coefficients.
####
```
a = [1,1]
b = [1,- 1]
o = polynomial_product(a, b); println(o)
 [1, 0, -1]
```
"""
function polynom_product(a::Vector{<:Number}, b::Vector{<:Number})

    n = Base.length(a)
    m = Base.length(b)

    if m ≥ n
        o = [Base.sum(a[1+j-i]*b[1+i] for i=0:j) for j=0:n-1]
        if m≠n Base.append!(o,[Base.sum(a[n-i]*b[1+i+j] for i=0:n-1) for j=1:m-n]) end
        Base.append!(o,[Base.sum(a[n-i]*b[1+i+j+m-n] for i=0:n-1-j) for j=1:n-1])
    else
        o = [Base.sum(b[1+j-i]*a[1+i] for i=0:j) for j=0:m-1]
        if m≠n Base.append!(o,[Base.sum(b[m-i]*a[1+i+j] for i=0:m-1) for j=1:n-m]) end
        Base.append!(o,[Base.sum(b[m-i]*a[1+i+j+n-m] for i=0:m-1-j) for j=1:m-1])
    end

    return o

end

# ==================================== polynom_product_expansion(a, b, p) ============================================================

@doc raw"""
    polynom_product_expansion(a::Vector{<:Number}, b::Vector{<:Number}, p::Int)

Vector representation of the product of two polynomials, ``a`` (of degree ``n``) and ``b`` (of degree ``m``), with ``m≤n``
truncated at the order ``p`` is a polynomial in a vector space of dimension ``d=p+1``. If ``ab`` is the `polynom_product`,
the `polynom_product_expansion` is ``ab[1:p+1]``
####
```
a = [1,-1,1]
b = [1,1,-1,1,1,1]
o = polynom_product(a, b); println(o)
 [1, 0, -1, 3, -1, 1, 0, 1]

o = expand_product(a, b, 4); println(o)
 [1, 0, -1, 3, -1]

```
"""
function polynom_product_expansion(a::Vector{<:Number}, b::Vector{<:Number}, p::Int)

    n = Base.length(a)
    m = Base.length(b)

    if m ≥ n
        o = [Base.sum(a[1+j-i]*b[1+i] for i=0:j) for j=0:min(n-1,p)]
        p+1 == length(o) && return o
        if m≠n Base.append!(o,[Base.sum(a[n-i]*b[1+i+j] for i=0:n-1) for j=1:min(m-n,p-n+1)]) end
        p+1 == length(o) && return o
        Base.append!(o,[Base.sum(a[n-i]*b[1+i+j+m-n] for i=0:n-1-j) for j=1:min(n-1,p-m+1)])
        p+1 == length(o) && return o
    else
        o = [Base.sum(b[1+j-i]*a[1+i] for i=0:j) for j=0:min(m-1,p)]
        p+1 == length(o) && return o
        Base.append!(o,[Base.sum(b[m-i]*a[1+i+j] for i=0:m-1) for j=1:min(n-m,p-m+1)])
        p+1 == length(o) && return o
        Base.append!(o,[Base.sum(b[m-i]*a[1+i+j+n-m] for i=0:m-1-j) for j=1:min(m-1,p-n+1)])
        p+1 == length(o) && return o
    end

    return o

end

# ...................................................... VectorRational .........................................................

@doc raw"""
    VectorRational

Object to decompose a vector of rational numbers

The fields are:
* `.num::Vector{Int}``: vector of normalized numerators
* `.den::Int`: common denominator
* `.val::Vector{Rational}`: vector of rational numbers (simplified = not normalized)
"""
struct VectorRational

    num::Vector{Int}
    den::Int
    val::Vector{Rational{Int64}}

end

# ==================================== analyzeVectorRational(vec) =======================

@doc raw"""
    normalize_VectorRational(vec::Vector{Rational{Int}})

Decompose vector of rational numbers.
#### Example:
```
v = [2//3,4//5]
normalize_VectorRational(v)
 VectorRational([10, 12], 15, Rational{Int64}[2//3, 4//5])
```
"""
function normalize_VectorRational(vec::Vector{Rational{Int}})

    val = Base.gcd(vec)
    den = val.den
    num = Base.convert(Vector{Int},(vec .* den))

    return CamiXon.VectorRational(num, den, vec)

end


text:
@enum LocTag begin
    O_Tag = 0
    B_Tag = 1
    I_Tag = 2
    E_Tag = 3
    S_Tag = 4
end

function toSym(t::LocTag)
    if t == O_Tag
        return :O
    elseif t == B_Tag
        return :B
    elseif t == I_Tag
        return :I
    elseif t == E_Tag
        return :E
    elseif t == S_Tag
        return :S
    end
end

function LocTag(s::AbstractString)
    if s == "O"
        return O_Tag
    elseif s == "B"
        return B_Tag
    elseif s == "I"
        return I_Tag
    elseif s == "E"
        return E_Tag
    elseif s == "S"
        return S_Tag
    else
        error("unknown loc tag: $s")
    end
end

struct NEType
    sym::Symbol
end

NEType(s::AbstractString="") = NEType(Symbol(s))

Base.isempty(t::NEType) = t.sym === Symbol()

Base.show(io::IO, t::NEType) = print(io, t.sym)

struct NETag
    loc::LocTag
    type::NEType
end

const prefix_form = Ref(r"^([BIES])-([^-]+)$")
const suffix_form = Ref(r"^([^-]+)-([BIES])$")

function parse_form(label::AbstractString)
    global prefix_form, suffix_form
    m = match(prefix_form[], label)
    if !isnothing(m)
        tag, type = m.captures
        return tag, type
    end

    m = match(suffix_form[], label)
    if !isnothing(m)
        type, tag = m.captures
        return tag, type
    end

    return nothing
end

get_tags(x::AbstractVector{<:AbstractString}) = map(NETag, x)
get_tags(x::AbstractVector{<:AbstractVector{<:AbstractString}}) = map(get_tags, x)

NETag() = NETag(O_Tag, NEType())

const nulltag = NETag()

function NETag(label::AbstractString)
    if length(label) == 1
        return NETag(LocTag(label), NEType())
    else
        p = parse_form(label)
        isnothing(p) && error("unknown ne tag: $label")
        tag, type = p
        return NETag(LocTag(tag), NEType(type))
    end
end

function Base.show(io::IO, t::NETag)
    print(io, toSym(t.loc))
    if !isempty(t.type)
        print(io, '-', t.type.sym)
    end
end

ext: jl


text:
using SparseArrays, BlockArrays, Base64
import BlockArrays: _BlockArray

function test_error_message(f, needle, expected = Exception)
    err = nothing
    try
        f()
    catch err
    end
    @test err isa expected
    @test occursin(needle, sprint(showerror ,err))
    return err
end

@testset "block constructors" begin
    ret = BlockArray{Float64}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1,Vector{Vector{Float64}}}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64}(undef, BlockArrays.BlockSizes(1:3))
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1}(undef, BlockArrays.BlockSizes(1:3))
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1,Vector{Vector{Float64}}}(undef, BlockArrays.BlockSizes(1:3))
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArrays._BlockArray([[0.0],[0.0,0.0],[0.0,0.0,0.0]], 1:3)
    @test Array(ret)  == zeros(6)

    ret = BlockArrays._BlockArray([[0.0],[0.0,0.0],[0.0,0.0,0.0]], BlockArrays.BlockSizes(1:3))
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float32}(undef_blocks, 1:3)
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float32,1}(undef_blocks, 1:3)
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float32,1,Vector{Vector{Float32}}}(undef_blocks, 1:3)
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float32}(undef_blocks, 1:3, 1:3)
    @test eltype(ret.blocks) == Matrix{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray(undef_blocks, Vector{Float32}, 1:3)
    @test eltype(ret) == Float32
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float64}(undef, 1:3, 1:3)
    fill!(ret, 0)
    @test Matrix(ret) == zeros(6,6)

    ret = PseudoBlockArray{Float64}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = PseudoBlockArray{Float64,1}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = PseudoBlockArray{Float64,1,Vector{Float64}}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = PseudoBlockArray{Float64}(undef, 1:3, 1:3)
    fill!(ret, 0)
    Matrix(ret) == zeros(6,6)

    @test_throws DimensionMismatch BlockArray([1,2,3],[1,1])

    @testset for sizes in [(1:3,), (1:3, 1:3), (1:3, 1:3, 1:3)]
        dims = sum.(sizes)
        A = BlockArray(copy(reshape(1:prod(dims), dims)), sizes...)
        @test mortar(A.blocks) == A
    end

    ret = mortar([spzeros(2), spzeros(3)])
    @test eltype(ret.blocks) <: SparseVector
    @test blocksizes(ret) == BlockArrays.BlockSizes([2, 3])

    ret = mortar(
        (spzeros(1, 3), spzeros(1, 4)),
        (spzeros(2, 3), spzeros(2, 4)),
        (spzeros(5, 3), spzeros(5, 4)),
     )
    @test Array(ret) == zeros(8, 7)
    @test eltype(ret.blocks) <: SparseMatrixCSC
    @test blocksizes(ret) == BlockArrays.BlockSizes([1, 2, 5], [3, 4])

    test_error_message("must have ndims consistent with ndims = 1") do
        mortar([ones(2,2)])
    end
    test_error_message("must have ndims consistent with ndims = 2") do
        mortar(reshape([ones(2), ones(2, 2)], (1, 2)))
    end
    test_error_message("size(blocks[2, 2]) (= (111, 222)) is incompatible with expected size: (2, 4)") do
        mortar(
            (zeros(1, 3), zeros(1, 4)),
            (zeros(2, 3), zeros(111, 222)),
        )
    end

    a_data = [1,2,3]
    a = BlockVector(a_data,[1,2])
    a[1] = 2
    @test a == [2,2,3]
    @test a_data == [1,2,3]
    a = BlockVector(a_data,BlockArrays.BlockSizes([1,2]))
    a[1] = 2
    @test a == [2,2,3]
    @test a_data == [1,2,3]
    a = PseudoBlockVector(a_data,[1,2])
    a[1] = 2
    @test a == [2,2,3]
    @test a_data == [2,2,3]
    a = PseudoBlockVector(a_data,BlockArrays.BlockSizes([1,2]))
    a[1] = 3
    @test a == [3,2,3]
    @test a_data == [3,2,3]

    a_data = [1 2; 3 4]
    a = BlockMatrix(a_data,[1,1],[2])
    a[1] = 2
    @test a == [2 2; 3 4]
    @test a_data == [1 2; 3 4]
    a = BlockMatrix(a_data,BlockArrays.BlockSizes([1,1],[2]))
    a[1] = 2
    @test a == [2 2; 3 4]
    @test a_data == [1 2; 3 4]
    a = PseudoBlockMatrix(a_data,[1,1],[2])
    a[1] = 2
    @test a == [2 2; 3 4]
    @test a_data == [2 2; 3 4]
    a = PseudoBlockMatrix(a_data,BlockArrays.BlockSizes([1,1],[2]))
    a[1] = 3
    @test a == [3 2; 3 4]
    @test a_data == [3 2; 3 4]
end

@testset "block indexing" begin
    BA_1 = BlockArray(undef_blocks, Vector{Float64}, [1,2,3])
    @test Base.IndexStyle(typeof(BA_1)) == IndexCartesian()

    a_1 = rand(2)
    BA_1[Block(2)] = a_1
    @test BA_1[BlockIndex(2, 1)] == a_1[1]
    @test BA_1[BlockIndex(2, 2)] == a_1[2]
    @test BA_1[Block(2)] == a_1
    @test BA_1[2] == a_1[1]

    @test_throws DimensionMismatch (BA_1[Block(3)] = rand(4))
    @test_throws BlockBoundsError blockcheckbounds(BA_1, 4)
    @test_throws BlockBoundsError BA_1[Block(4)]

    BA_2 = BlockArray(undef_blocks, Matrix{Float64}, [1,2], [3,4])
    a_2 = rand(1,4)
    BA_2[Block(1,2)] = a_2
    @test BA_2[Block(1,2)] == a_2
    BA_2[Block(1,2)] = a_2

    @test BA_2[1,5] == a_2[2]
    @test_throws DimensionMismatch BA_2[Block(1,2)] = rand(1,5)
end

@testset "misc block tests" begin
    for BlockType in (BlockArray, PseudoBlockArray)
        a_1 = rand(6)
        BA_1 = BlockType(a_1, [1,2,3])
        @test Array(BA_1) == a_1
        @test nblocks(BA_1) == (3,)
        @test nblocks(BA_1,1) == 3
        @test eltype(similar(BA_1, Float32)) == Float32
        q = rand(1)
        BA_1[Block(1)] = q
        BA_1[BlockIndex(3, 2)] = a_1[5]
        @test BA_1[Block(1)] == q
        if BlockType == PseudoBlockArray
            q2 = zero(q)
            getblock!(q2, BA_1, 1)
            @test q2 == q
            @test_throws DimensionMismatch getblock!(zeros(2), BA_1, 1)
            fill!(q2, 0)
            getblock!(q2, BA_1, 1)
            @test q2 == q
        end
        fill!(BA_1, 1.0)
        @test BA_1 == ones(size(BA_1))
        ran = rand(size(BA_1)...)
        copyto!(BA_1, ran)
        @test BA_1 == ran

        a_1_sparse = sprand(6, 0.9)
        BA_1_sparse = BlockType(a_1_sparse, [1,2,3])
        @test Array(BA_1_sparse) == a_1_sparse
        BA_1_sparse[4] = 3.0
        @test BA_1_sparse[4] == 3.0

        a_2 = rand(3, 7)
        BA_2 = BlockType(a_2, [1,2], [3,4])
        @test Array(BA_2) == a_2
        @test nblocks(BA_2) == (2,2)
        @test nblocks(BA_2, 1) == 2
        @test nblocks(BA_2, 2, 1) == (2, 2)
        BA_2[BlockIndex((2,1), (2,2))] = a_2[3,2]
        @test eltype(similar(BA_2, Float32)) == Float32
        q = rand(1,4)
        BA_2[Block(1,2)] = q
        @test_throws DimensionMismatch BA_2[Block(1,2)] = rand(1,5)
        @test BA_2[Block(1,2)] == q
        if BlockType == PseudoBlockArray
            q2 = zero(q)
            getblock!(q2, BA_2, 1, 2)
            @test q2 == q
            @test_throws DimensionMismatch getblock!(zeros(1,5), BA_2, 1, 2)
        end
        fill!(BA_2, 1.0)
        @test BA_2 == ones(size(BA_2))
        ran = rand(size(BA_2)...)
        copyto!(BA_2, ran)
        @test BA_2 == ran

        a_2_sparse = sprand(3, 7, 0.9)
        BA_2_sparse = BlockType(a_2_sparse, [1,2], [3,4])
        @test Array(BA_2_sparse) == a_2_sparse
        BA_2_sparse[1,2] = 3.0
        @test BA_2_sparse[1,2] == 3.0

        a_3 = rand(3, 7,4)
        BA_3 = BlockType(a_3, [1,2], [3,4], [1,2,1])
        @test Array(BA_3) == a_3
        @test nblocks(BA_3) == (2,2,3)
        @test nblocks(BA_3, 1) == 2
        @test nblocks(BA_3, 3, 1) == (3, 2)
        @test nblocks(BA_3, 3) == 3
        BA_3[BlockIndex((1,1,1), (1,1,1))] = a_3[1,1,1]
        @test eltype(similar(BA_3, Float32)) == Float32
        q = rand(1,4,2)
        BA_3[Block(1,2,2)] = q
        @test BA_3[Block(1,2,2)] == q
        if BlockType == PseudoBlockArray
            q3 = zero(q)
            getblock!(q3, BA_3, 1, 2, 2)
            @test q3 == q
            @test_throws DimensionMismatch getblock!(zeros(1,3,2), BA_3, 1, 2,2)
        end
        fill!(BA_3, 1.0)
        @test BA_3 == ones(size(BA_3))
        ran = rand(size(BA_3)...)
        copyto!(BA_3, ran)
        @test BA_3 == ran
    end
end

@testset "convert" begin
    # Could probably be DRY'd.
    A = PseudoBlockArray(rand(2,3), [1,1], [2,1])
    C = convert(BlockArray, A)
    @test C == A == BlockArray(A)
    @test eltype(C) == eltype(A)

    C = convert(BlockArray{Float32}, A)
    @test C ≈ A ≈ BlockArray(A)
    @test eltype(C) == Float32

    C = convert(BlockArray{Float32, 2}, A)
    @test C ≈ A ≈ BlockArray(A)
    @test eltype(C) == Float32

    A = BlockArray(rand(2,3), [1,1], [2,1])
    C = convert(PseudoBlockArray, A)
    @test C == A == PseudoBlockArray(A)
    @test eltype(C) == eltype(A)

    C = convert(PseudoBlockArray{Float32}, A)
    @test C ≈ A ≈ PseudoBlockArray(A)
    @test eltype(C) == Float32

    C = convert(PseudoBlockArray{Float32, 2}, A)
    @test C ≈ A ≈ PseudoBlockArray(A)
    @test eltype(C) == Float32
end

@testset "string" begin
    A = BlockArray(rand(4, 5), [1,3], [2,3]);
    buf = IOBuffer()
    Base.showerror(buf, BlockBoundsError(A, (3,2)))
    @test String(take!(buf)) == "BlockBoundsError: attempt to access 2×2-blocked 4×5 BlockArray{Float64,2,Array{Array{Float64,2},2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}} at block index [3,2]"

    A = PseudoBlockArray(rand(4, 5), [1,3], [2,3]);
    Base.showerror(buf, BlockBoundsError(A, (3,2)))
    @test String(take!(buf)) == "BlockBoundsError: attempt to access 2×2-blocked 4×5 PseudoBlockArray{Float64,2,Array{Float64,2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}} at block index [3,2]"
end

@testset "replstring" begin
    @test stringmime("text/plain",BlockArray(collect(reshape(1:16, 4, 4)), [1,3], [2,2])) == "2×2-blocked 4×4 BlockArray{Int64,2}:\n 1  5  │   9  13\n ──────┼────────\n 2  6  │  10  14\n 3  7  │  11  15\n 4  8  │  12  16"
    @test stringmime("text/plain",PseudoBlockArray(collect(reshape(1:16, 4, 4)), [1,3], [2,2])) == "2×2-blocked 4×4 PseudoBlockArray{Int64,2}:\n 1  5  │   9  13\n ──────┼────────\n 2  6  │  10  14\n 3  7  │  11  15\n 4  8  │  12  16"
    design = zeros(Int16,6,9);
    A = BlockArray(design,[6],[4,5])
    @test stringmime("text/plain",A) == "1×2-blocked 6×9 BlockArray{Int16,2}:\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0"
    A = PseudoBlockArray(design,[6],[4,5])
    @test stringmime("text/plain",A) == "1×2-blocked 6×9 PseudoBlockArray{Int16,2}:\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0"
end

@testset "AbstractVector{Int} blocks" begin
    A = BlockArray(ones(6,6), 1:3, 1:3)
    @test A[1,1] == 1
    @test A[Block(2,3)] == ones(2,3)

    A = BlockArray(undef_blocks, Matrix{Float64}, 1:3, 1:3)
    A[Block(2,3)] = ones(2,3)
    @test A[Block(2,3)] == ones(2,3)
end

@testset "Block arithmetic" begin
    @test +(Block(1)) == Block(1)
    @test -(Block(1)) == Block(-1)
    @test Block(2) + Block(1) == Block(3)
    @test Block(2) + 1 == Block(3)
    @test 2 + Block(1) == Block(3)
    @test Block(2) - Block(1) == Block(1)
    @test Block(2) - 1 == Block(1)
    @test 2 - Block(1) == Block(1)
    @test 2*Block(1) == Block(2)
    @test Block(1)*2 == Block(2)

    @test isless(Block(1), Block(2))
    @test !isless(Block(1), Block(1))
    @test !isless(Block(2), Block(1))
    @test Block(1) < Block(2)
    @test Block(1) ≤ Block(1)
    @test Block(2) > Block(1)
    @test Block(1) ≥ Block(1)
    @test min(Block(1), Block(2)) == Block(1)
    @test max(Block(1), Block(2)) == Block(2)

    @test +(Block(1,2)) == Block(1,2)
    @test -(Block(1,2)) == Block(-1,-2)
    @test Block(1,2) + Block(2,3) == Block(3,5)
    @test Block(1,2) + 1 == Block(2,3)
    @test 1 + Block(1,2) == Block(2,3)
    @test Block(2,3) - Block(1,2) == Block(1,1)
    @test Block(1,2) - 1 == Block(0,1)
    @test 1 - Block(1,2) == Block(0,-1)
    @test 2*Block(1,2) == Block(2,4)
    @test Block(1,2)*2 == Block(2,4)

    @test isless(Block(1,1), Block(2,2))
    @test isless(Block(1,1), Block(2,1))
    @test !isless(Block(1,1), Block(1,1))
    @test !isless(Block(2,1), Block(1,1))
    @test Block(1,1) < Block(2,1)
    @test Block(1,1) ≤ Block(1,1)
    @test Block(2,1) > Block(1,1)
    @test Block(1,1) ≥ Block(1,1)
    @test min(Block(1,2), Block(2,2)) == Block(1,2)
    @test max(Block(1,2), Block(2,2)) == Block(2,2)

    @test convert(Int, Block(2)) == 2
    @test convert(Float64, Block(2)) == 2.0

    @test_throws MethodError convert(Int, Block(2,1))
    @test convert(Tuple{Int,Int}, Block(2,1)) == (2,1)
    @test convert(Tuple{Float64,Int}, Block(2,1)) == (2.0,1)
end

@testset "Strided array interface" begin
    A = PseudoBlockArray{Float64}(undef, 1:3, 1:3)
    fill!(A, 1)
    @test strides(A) == (1, size(A,1))
    x = randn(size(A,2))
    y = similar(x)
    @test BLAS.gemv!('N', 2.0, A, x, 0.0, y) ≈ 2A*x
end

@testset "lmul!/rmul!" begin
    A = PseudoBlockArray{Float64}(undef, 1:3)
    @test fill!(A, NaN) === A
    @test all(isnan, lmul!(0.0, copy(A))) == all(isnan, lmul!(0.0, Array(A)))
    @test lmul!(false, copy(A)) == lmul!(false, Array(A))
    @test lmul!(0.0, A) === A
    @test fill!(A, NaN) === A
    @test all(isnan, rmul!(copy(A), 0.0)) == all(isnan, rmul!(Array(A), 0.0))
    @test rmul!(copy(A), false) == rmul!(Array(A), false)
    @test rmul!(A, 0.0) === A

    A = BlockArray{Float64}(undef, 1:3)
    @test fill!(A, NaN) === A
    @test all(isnan, lmul!(0.0, copy(A))) == all(isnan, lmul!(0.0, Array(A)))
    @test lmul!(false, copy(A)) == lmul!(false, Array(A))
    @test lmul!(0.0, A) === A
    @test fill!(A, NaN) === A
    @test all(isnan, rmul!(copy(A), 0.0)) == all(isnan, rmul!(Array(A), 0.0))
    @test rmul!(copy(A), false) == rmul!(Array(A), false)
    @test rmul!(A, 0.0) === A
end

@testset "copy" begin
    A = PseudoBlockArray(randn(6), 1:3)
    B = copy(A)
    @test typeof(A) == typeof(B)
    @test blocksizes(A) == blocksizes(B)
    @test A == B
    B[1] = 2
    @test B[1] == 2
    @test A[1] ≠ 2

    A = BlockArray(randn(6), 1:3)
    B = copy(A)
    @test typeof(A) == typeof(B)
    @test blocksizes(A) == blocksizes(B)
    @test A == B
    B[1] = 2
    @test B[1] == 2
    @test A[1] ≠ 2
end

@testset "const block size" begin
    N = 10
    # In the future this will be automated via mortar(..., Fill(2,N))
    A = mortar(fill([1,2], N), BlockArrays.BlockSizes((1:2:2N+1,)))
    B = PseudoBlockArray(vcat(fill([1,2], N)...), BlockArrays.BlockSizes((1:2:2N+1,)))
    @test A == vcat(A.blocks...) == B
    @test A[Block(1)] == B[Block(1)] == [1,2]
end

@testset "mixed blocks #92" begin
    A = BlockArray{Any}(undef_blocks, [1,2], [3,2])
    A.blocks[1,1] = fill(1.0,1,3)
    A.blocks[1,2] = fill(2,1,2)
    A.blocks[2,1] = fill(3//3,2,3)
    A.blocks[2,2] = fill(4.0f0,2,2)
    @test A[1,1] == 1.0
    @test A[2,1] == 1//1

    B = Matrix{Matrix}(undef,2,2)
    B[1,1] = fill(1.0,1,3)
    B[1,2] = fill(2,1,2)
    B[2,1] = fill(3//3,2,3)
    B[2,2] = fill(4.0f0,2,2)
    @test mortar(B) == A
    @test mortar(B)[1,1] == 1.0
end


text:
using LightGraphs
using MatrixDepot
function symmetrize(A)
    println("Symmetrizing ")
    tic()
    if !issymmetric(A)
        println(STDERR, "the matrix is not symmetric using A+A'")
        A = A + A'
    end
    if isa(A, Base.LinAlg.Symmetric)
        A.data.nzval = abs(A.data.nzval)
    else
        A.nzval = abs(A.nzval)
    end
    #= spA = abs(sparse(A)) =#
    toc()
    return A
end

function loadmat(matname)
    println("Reading MTX for $matname")
    tic()
    try
        info = matrixdepot(matname, :get)
    end
    A = matrixdepot(matname, :read)
    A = symmetrize(A)
    g = Graph(A)
    toc()
    println(STDERR, "size(A) = $(size(A))")
    return g
end

function vec2tree(v::Vector{Int})
    nv = length(v)
    I = Vector{Int}()
    J = Vector{Int}()
    for i in 1:nv
        if v[i] > 0
            push!(I, i)
            push!(J, v[i])
        end
    end
    #= @assert length(I) == length(J) =#
    ncolumns = max(maximum(J), maximum(I))
    nentries = length(J)
    n = ncolumns
    return sparse(I,J, fill(1, nentries), n,n)
end

names = ["Newman/football" , "Newman/cond-mat-2003", "SNAP/amazon0302", "SNAP/roadNet-CA"]
for matname in names
    g = loadmat(matname)
    seed = 1
    println("bfs_tree original")
    # woah bfs_tree allocates a lot of memory
    @time tdg = LightGraphs.bfs_tree(g, seed)
    println(tdg)
    # preallocating the output tree reduces the number of allocations.
    # much faster
    println("bfs_tree!(vector)")
    visitor = LightGraphs.TreeBFSVisitorVector(zeros(Int, nv(g)))
    @time tvec = LightGraphs.bfs_tree!(visitor, g, seed)
    println("converting to Sparse")
    @time m = vec2tree(visitor.tree)
    println("converting to DiGraph")
    @time h = DiGraph(m)
end


text:
#!/usr/bin/env julia

using Plots
pyplot()

julia1 = readdlm("julia_times-1.dat")
julia4 = readdlm("julia_times-4.dat")
python = readdlm("python_times.dat")

plot(xaxis = (:log,), yaxis = (:log,),
     xlabel = "Datapoints", ylabel = "Time (seconds)",
     size=(900, 600))
plot!(julia1[:,1], julia1[:,2], linewidth = 2, marker = (:auto,), color = :blue,   lab = "LombScargle.jl - single thread")
plot!(julia4[:,1], julia4[:,2], linewidth = 2, marker = (:auto,), color = :orange, lab = "LombScargle.jl - 4 threads")
plot!(python[:,1], python[:,2], linewidth = 2, marker = (:auto,), color = :green,  lab = "Astropy")
savefig("benchmarks.svg")
savefig("benchmarks.png")


text:
# ---
# title: 560. Subarray Sum Equals K
# id: problem560
# author: Indigo
# date: 2021-02-16
# difficulty: Medium
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/subarray-sum-equals-k/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and an integer `k`, return _the total number
# of continuous subarrays whose sum equals to`k`_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,1], k = 2
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3], k = 3
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2 * 104`
#   * `-1000 <= nums[i] <= 1000`
#   * `-107 <= k <= 107`
#
#
## @lc code=start
using LeetCode

function subarray_sum(nums::Vector{Int}, k::Int)
    mp = DefaultDict(0, 0 => 1)
    res, pre_sum = 0, 0
    for num in nums
        pre_sum += num
        res += mp[pre_sum - k]
        mp[pre_sum] += 1
    end
    res
end
## @lc code=end


text:
const Match = Tuple{Rule, Union{Nothing, Pattern}, Sub, Int64}
const MatchesBuf = Vector{Match}

function cached_ids(g::EGraph, p::Pattern)#::Vector{Int64}
    keys(g.classes)
end

# FIXME
function cached_ids(g::EGraph, p::PatTerm)
    # println(p)
    # println(get(g.symcache, p.head, []))
    keys(g.classes)
    # get(g.symcache, p.head, [])
end

# function cached_ids(g::EGraph, p::PatLiteral)
#     get(g.symcache, p.val, [])
# end

function search_rule!(g::EGraph, r::SymbolicRule, id::Int64, matches::MatchesBuf)
    for sub in ematch(g, r.left, id)
        push!(matches, (r, r.right, sub, id))
    end
end

function search_rule!(g::EGraph, r::DynamicRule, id::Int64, matches::MatchesBuf)
    for sub in ematch(g, r.left, id)
        push!(matches, (r, nothing, sub, id))
    end
end

function search_rule!(g::EGraph, r::BidirRule, id::Int64, matches::MatchesBuf)
    for sub in ematch(g, r.left, id)
        push!(matches, (r, r.right, sub, id))
    end
    for sub in ematch(g, r.right, id)
        push!(matches, (r, r.left, sub, id))
    end
end

function Base.show(io::IO, s::Sub)
    print(io, "Sub[")
    kvs = collect(s)
    n = length(kvs)
    for i ∈ 1:n
        print(io, kvs[i][1], " => ", kvs[i][2][1].id)
        if i < n
            print(io, ",")
        end
    end
    print(io, "]")
end

function search_rule!(g::EGraph, r::MultiPatRewriteRule, id::Int64, matches::MatchesBuf)
    buf = ematch(g, r.left, id)
    if isempty(buf)
        return
    end
    pats_todo = reverse(copy(r.pats))
    while !isempty(pats_todo)
        pat = pop!(pats_todo)
        # println("====================")
        # @show pat
        ids = cached_ids(g, pat)
        newbuf = SubBuf()
        while !isempty(buf)
            sub = pop!(buf)
            # @show sub
            isempty(sub) && continue
            for i ∈ ids
                ematch(g, pat, i, sub, newbuf)
            end
        end
        buf = copy(newbuf)
    end
    for sub in buf
        # println("FINALLY ", sub, " $id")
        push!(matches, (r, r.right, sub, id))
    end
end

function eqsat_search!(egraph::EGraph, theory::Vector{<:Rule},
        scheduler::AbstractScheduler)::MatchesBuf
    matches=MatchesBuf()
    for rule ∈ theory
        # don't apply banned rules
        if shouldskip(scheduler, rule)
            # println("skipping banned rule $rule")
            continue
        end

        ids = cached_ids(egraph, rule.left)

        for id ∈ ids
            search_rule!(egraph, rule, id, matches)
        end
    end
    return matches
end


text:
using BasisFunctions, BasisFunctions.Test, DoubleFloats

using Test
⊕ = BasisFunctions.:⊕
⊗ = BasisFunctions.:⊗

types = (Float64,LargeFloat)

function test_operators(T)
    @testset "identity operator" begin
        test_identity_operator(T) end

    @testset "diagonal operators" begin
        test_diagonal_operators(T) end

    @testset "multidiagonal operators" begin
        test_multidiagonal_operators(T) end

    @testset "test invertible operators" begin
        test_invertible_operators(T) end

    @testset "noninvertible operators" begin
        test_noninvertible_operators(T) end

    @testset "tensor operators" begin
        test_tensor_operators(T) end

    @testset "circulant operator" begin
        test_circulant_operator(T) end

    @testset "banded operators" begin
        test_banded_operator(T) end

    @testset "sparse operator" begin
        test_sparse_operator(T) end
end

function test_generic_operators(T)
    b1 = Fourier{T}(3)
    b2 = ChebyshevT{T}(4)
    b3 = ChebyshevT{T}(3)
    b4 = Legendre{T}(3)

    cartlist = [CartesianIndex(1,1),CartesianIndex(2,3)]
    operators = [
        ["Scaling operator", ScalingOperator(b1, T(2))],
        ["Zero operator", BF.ZeroOperator(b1, b2)],
        ["Diagonal operator", DiagonalOperator(b2, b2, rand(T, length(b2)))],
        ["Block operator", differentiation(b1⊕b2)],
        ["Coefficient scaling operator", BF.CoefficientScalingOperator(b1, 1, T(2))],
        ["Index restriction operator", IndexRestriction(b2, b1, 1:3) ],
        ["Index restriction operator (2)", IndexRestriction(b1⊗b2, cartlist) ],
        ["Index extension operator", IndexExtension(b1, b2, 1:3) ],
        ["Index extension operator (2)", IndexExtension(b1⊗b2, cartlist) ],
        ["Linearization operator", BasisFunctions.LinearizationOperator(b1⊗b2) ],
        ["DelinearizationOperator operator", BasisFunctions.DelinearizationOperator(b1⊗b2) ],
    ]

    for ops in operators
        @testset "$(ops[1])" begin
            test_generic_operator_interface(ops[2], T)
        end
    end
end

function test_tensor_operators(T)
    m1 = 3
    n1 = 4
    m2 = 10
    n2 = 24
    A1 = ArrayOperator(rand(T,m1,n1))
    A2 = ArrayOperator(rand(T,m2,n2))

    A_tp = TensorProductOperator(A1, A2)
    b = rand(T, n1, n2)
    c_tp = zeros(T, m1, m2)
    apply!(A_tp, c_tp, b)

    # Apply the operators manually, row by row and then column by column
    intermediate = zeros(T, m1, n2)
    result1 = zeros(T, m1)
    for i = 1:n2
        apply!(A1, result1, b[:,i])
        intermediate[:,i] = result1
    end
    c_manual = similar(c_tp)
    result2 = zeros(T, m2)
    for i = 1:m1
        apply!(A2, result2, reshape(intermediate[i,:], n2))
        c_manual[i,:] = result2
    end

    @test sum(abs.(c_tp-c_manual)) + 1 ≈ 1
end

function test_diagonal_operators(T)
    for SRC in (Fourier{T}(10), ChebyshevT{T}(11))
        operators = (BF.CoefficientScalingOperator(SRC, 3, rand(coefficienttype(SRC))),
            BF.AlternatingSignOperator(SRC), IdentityOperator(SRC),
            ScalingOperator(SRC, rand(coefficienttype(SRC))), ScalingOperator(SRC,3),
            DiagonalOperator(SRC, rand(coefficienttype(SRC), size(SRC))))
           # PseudoDiagonalOperator(SRC, map(coefficienttype(SRC), rand(size(SRC)))))
        for Op in operators
            m = matrix(Op)
            # Test in-place
            if isinplace(Op)
                coef_src = rand(src(Op))
                coef_dest_m = m * coef_src
                coef_dest = apply!(Op, coef_src)
                @test sum(abs.(coef_dest-coef_dest_m)) + 1 ≈ 1
            end
            # Test out-of-place
            coef_src = rand(src(Op))
            coef_dest = apply(Op, coef_src)
            coef_dest_m = m * coef_src
            @test sum(abs.(coef_dest-coef_dest_m)) + 1 ≈ 1
            # Test inverse
            I = inv(Op)
            @test abs(sum(abs.(I*(Op*coef_src)-coef_src))) + 1 ≈ 1
            @test abs(sum(abs.((I*Op)*coef_src-coef_src))) + 1 ≈ 1
            # Test transposes
            @test abs(sum(abs.(m'*coef_src-Op'*coef_src))) + 1 ≈ 1
            # Test Sum
            @test abs(sum(abs.((Op+Op)*coef_src-2*(Op*coef_src)))) + 1 ≈ 1
            # Test Equivalence to diagonal operator
            @test abs(sum(abs.(Op*coef_src-diag(Op).*coef_src))) + 1 ≈ 1
            # Make sure diagonality is retained
            @test isdiag(2*Op)
            @test isdiag(Op')
            @test isdiag(inv(Op))
            @test isdiag(Op*Op)
            # Make sure in_place is retained
            if isinplace(Op)
                @test isinplace(2*Op)
                @test isinplace(Op')
                @test isinplace(inv(Op))
                @test isinplace(Op*Op)
            end
        end
    end
end

function test_multidiagonal_operators(T)
    MSet = Fourier{T}(10)⊕ChebyshevT{T}(11)
    operators = (BF.CoefficientScalingOperator(MSet, 3, rand(coefficienttype(MSet))),
        BF.AlternatingSignOperator(MSet), IdentityOperator(MSet),
        ScalingOperator(MSet,2.0+2.0im), ScalingOperator(MSet, 3),
        DiagonalOperator(MSet, rand(coefficienttype(MSet),length(MSet))))
    for Op in operators
        # Test in-place
        coef_src = zeros(MSet)
        for i in eachindex(coef_src)
            coef_src[i] = rand(eltype(coef_src))
        end
        m = matrix(Op)
        coef_dest_m = m * linearize_coefficients(MSet,coef_src)
        coef_dest = linearize_coefficients(MSet, apply!(Op, coef_src))
        @test sum(abs.(coef_dest-coef_dest_m)) + 1 ≈ 1
        # Test out-of-place
        coef_src = zeros(MSet)
        for i in eachindex(coef_src)
            coef_src[i] = rand() * one(eltype(coef_src))
        end
        coef_dest = linearize_coefficients(MSet, apply(Op, coef_src))
        coef_dest_m = m * linearize_coefficients(MSet, coef_src)
        @test sum(abs.(coef_dest-coef_dest_m)) + 1 ≈ 1

        # Test inverse
        I = inv(Op)
        @test I*(Op*coef_src)≈coef_src
        @test (I*Op)*coef_src≈coef_src
        # Test transposes
        @test m'*linearize_coefficients(MSet,coef_src)≈linearize_coefficients(MSet,Op'*coef_src)
        # Test Sum
        @test (Op+Op)*coef_src≈2*(Op*coef_src)
        # Test Equivalence to diagonal operator
        @test linearize_coefficients(MSet,Op*coef_src)≈ diag(Op).*linearize_coefficients(MSet,coef_src)
        # Make sure diagonality is retained
        @test isdiag(2*Op) && isinplace(2*Op)
        @test isdiag(Op') && isinplace(Op')
        @test isdiag(inv(Op)) && isinplace(inv(Op))
        @test isdiag(Op*Op) && isinplace(Op*Op)
    end
end

function test_sparse_operator(ELT)
    S = SparseMatrixOperator(ArrayOperator(rand(ELT,4,4)))
    test_generic_operator_interface(S, ELT)
end

function test_banded_operator(ELT)
    a = [ELT(1),ELT(2),ELT(3)]
    H = HorizontalBandedOperator(Fourier{ELT}(6), Fourier{ELT}(3),a,3,2)
    test_generic_operator_interface(H, ELT)
    h = matrix(H)
    e = zeros(ELT,6,1)
    for i in 1:3
        e .= 0
        I = (2+(i-1)*3) .+ (1:3)
        for (j,k) in enumerate(mod.(I .- 1,6) .+ 1)
            e[k] = a[j]
        end
        @test e ≈ h[i,:]
    end

    V = VerticalBandedOperator(Fourier{ELT}(3), Fourier{ELT}(6),a,3,2)
    test_generic_operator_interface(V, ELT)
    v = matrix(V)
    e = zeros(ELT,6,1)
    for i in 1:3
        e .= 0
        I = (2+(i-1)*3) .+ (1:3)
        for (j,k) in enumerate(mod.(I .- 1,6) .+ 1)
            e[k] = a[j]
        end
        @test e ≈ v[:,i]
    end
end

function test_circulant_operator(ELT)
    n = 20
    for T in (ELT, complex(ELT))
        e1 = zeros(T,n); e1[1] = 1
        c = rand(T, n)
        C = CirculantOperator(c)
        m = matrix(C)
        for i in 1:n
            @test m[:,i] ≈ circshift(c,i-1)
        end
        invC = inv(C)
        @test BasisFunctions.eigvals(invC).*BasisFunctions.eigvals(C) ≈ ones(T,n)
        @test m'*e1 ≈ C'*e1

        sumC = invC+C
        @test (typeof(sumC) <:BasisFunctions.CirculantOperator)
        @test sumC*e1 ≈ m[:,1] + matrix(invC)[:,1]

        subC = invC-C
        @test (typeof(subC) <:BasisFunctions.CirculantOperator)
        @test subC*e1 ≈ -m[:,1] + matrix(invC)[:,1]

        mulC = invC*C
        @test (typeof(sumC) <:BasisFunctions.CirculantOperator)
        @test mulC*e1 ≈ (m*matrix(invC))[:,1]

        mulC = 2*C
        @test (typeof(sumC) <:BasisFunctions.CirculantOperator)
        @test mulC*e1 ≈ 2*C*e1

        mulC = C*2
        @test (typeof(sumC) <:BasisFunctions.CirculantOperator)
        @test mulC*e1 ≈ 2*C*e1
    end
end

function test_identity_operator(T)
    b = Fourier{T}(10)
    I = IdentityOperator(b)
    @test I isa IdentityOperator
    @test eltype(I) == Complex{T}
    @test src(I) == b
    I = IdentityOperator(Fourier{T}(10), ChebyshevT{T}(10))
    @test src(I) == Fourier{T}(10)
    @test dest(I) == ChebyshevT{T}(10)
end

function test_invertible_operators(T)
    for SRC in (Fourier{T}(10), ChebyshevT{T}(10))
        operators = (MultiplicationOperator(SRC, SRC, rand(coefficienttype(SRC),length(SRC),length(SRC))),
            CirculantOperator(rand(T, 10)),
            CirculantOperator(rand(complex(T), 10)))
        for Op in operators
            m = matrix(Op)
            # Test out-of-place
            coef_src = rand(src(Op))
            coef_dest = apply(Op, coef_src)
            coef_dest_m = m * coef_src
            @test sum(abs.(coef_dest-coef_dest_m)) + 1 ≈ 1
            # Test inverse
            I = inv(Op)
            @test abs(sum(abs.(I*(Op*coef_src)-coef_src))) + 1 ≈ 1
            @test abs(sum(abs.((I*Op)*coef_src-coef_src))) + 1 ≈ 1
            # Test transposes
            @test abs(sum(abs.(m'*coef_src-Op'*coef_src))) + 1 ≈ 1
            # Test Sum
            @test abs(sum(abs.((Op+Op)*coef_src-2*(Op*coef_src)))) + 1 ≈ 1
        end
    end
end

# FunctionOperator is currently not tested, since we cannot assume it is a linear operator.
function test_noninvertible_operators(T)
    for SRC in (Fourier{T}(10), ChebyshevT{T}(11))
        operators = (MultiplicationOperator(SRC, resize(SRC,length(SRC)+2), rand(coefficienttype(SRC),length(SRC)+2,length(SRC))),
            BF.ZeroOperator(SRC))
        for Op in operators
            m = matrix(Op)
            # Test out-of-place
            coef_src = rand(src(Op))
            coef_dest = apply(Op, coef_src)
            coef_dest_m = m * coef_src
            @test sum(abs.(coef_dest-coef_dest_m)) + 1 ≈ 1
            # Test transposes
            @test (sum(abs.(m'-matrix(Op')))) + 1 ≈ 1
            # Test Sum
            @test (sum(abs.((Op+Op)*coef_src-2*(Op*coef_src)))) + 1 ≈ 1
        end
    end
end

for T in types
    delimit(string(T))
    test_operators(T)
    test_generic_operators(T)
end


text:
using CUDA

function cuda_kernel(
    xyz::CuDeviceVector{T},
    forces::CuDeviceVector{T},
    energy::CuDeviceVector{T},
    N::Int,
    cutsq::T) where {T<:AbstractFloat}

    blocksize = blockDim().x

    # thread ID
    tid = (blockIdx().x-1) * blocksize + threadIdx().x

    if tid <= N
        id = threadIdx().x

        shmem = @cuDynamicSharedMem(T, 3*blocksize)

        ii = tid<<2 - (2 + tid) # AoS

        e  = T(0)
        fx = T(0)
        fy = T(0)
        fz = T(0)

        xi = xyz[ii]
        yi = xyz[ii+1]
        zi = xyz[ii+2]

        for i=1:blocksize:N
            k = ((i+id)<<1) + (i+id-2) # AoS
            if (i+id) <= N
                shmem[(id-1)*3+1] = xyz[k]
                shmem[(id-1)*3+2] = xyz[k+1]
                shmem[(id-1)*3+3] = xyz[k+2]

                sync_threads()

                #if i!=tid
                for j=1:blocksize
                    dx = shmem[(j-1)*3+1] - xi
                    dy = shmem[(j-1)*3+2] - yi
                    dz = shmem[(j-1)*3+3] - zi
                    dij_sq = dx*dx + dy*dy + dz*dz

                    if dij_sq <= cutsq && dij_sq > 1e-12
                        lj2 = T(1) / dij_sq
                        lj6 = lj2*lj2*lj2
                        e += lj6*lj6 - lj6

                        fc = T(24.0) * (lj6 - T(2.0) * lj6 * lj6) / dij_sq
                        fx += fc*dx
                        fy += fc*dy
                        fz += fc*dz
                    end

                end
                #end
                sync_threads()
            end
        end

        energy[tid]  = e/2
        forces[ii  ] = fx
        forces[ii+1] = fy
        forces[ii+2] = fz
    end

    return nothing
end

function cuda(state::State{F, T}, cutoff::T) where {F <: AbstractMatrix, T <: AbstractFloat}

    N = state.n
    x_gpu = CuArray(view(state.coords.values, :))
    f_gpu = similar(x_gpu)
    e_gpu = CuArray{T}(undef, N)

    total_threads = min(N, 512)
    threads       = (total_threads, )
    blocks        = ceil.(Int, N ./ threads)

    # calculate size of dynamic shared memory
    shmem = 3 * sum(threads) * sizeof(T)
    cutsq = convert(T, cutoff*cutoff)

    @cuda blocks=blocks threads=threads shmem=shmem cuda_kernel(x_gpu, f_gpu, e_gpu, N, cutsq)

    e = sum(e_gpu)
    return e
end


text:
# This file contains functions for the Profit Efficiency DEA model
"""
    ProfitDEAModel
An data structure representing a profit DEA model.
"""
struct ProfitDEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
end

"""
    deaprofit(X, Y, P)
Compute profit efficiency using data envelopment analysis model for
inputs `X`, outputs `Y`, price of inputs `W`, and price of outputs `P`.

# Examples
```jldoctest
julia> X = [1 1; 1 1; 0.75 1.5; 0.5 2; 0.5 2; 2 2; 2.75 3.5; 1.375 1.75];

julia> Y = [1 11; 5 3; 5 5; 2 9; 4 5; 4 2; 3 3; 4.5 3.5];

julia> P = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1];

julia> W = [2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1; 2 1];

julia> GxGydollar = 1 ./ (sum(P, dims = 2) + sum(W, dims = 2));

julia> Gx = repeat(GxGydollar, 1, 2);

julia> Gy = repeat(GxGydollar, 1, 2);

julia> deaprofit(X, Y, W, P, Gx, Gy)
Profit DEA Model
DMUs = 8; Inputs = 2; Outputs = 2
Returns to Scale = VRS
─────────────────────────────────────
   Profit     Technical    Allocative
─────────────────────────────────────
1     2.0   0.0           2.0
2     2.0  -5.41234e-16   2.0
3     0.0   0.0           0.0
4     2.0   0.0           2.0
5     2.0   0.0           2.0
6     8.0   6.0           2.0
7    12.0  12.0          -1.77636e-15
8     4.0   3.0           1.0
─────────────────────────────────────
```
"""
function deaprofit(X::Matrix, Y::Matrix, W::Matrix, P::Matrix, Gx::Matrix, Gy::Matrix)::ProfitDEAModel
    # Check parameters
    nx, m = size(X)
    ny, s = size(Y)

    nw, mw = size(W)
    np, sp = size(P)

    if nx != ny
        error("number of observations is different in inputs and outputs")
    end
    if nw != nx
        error("number of observations is different in input prices and inputs")
    end
    if np != ny
        error("number of observations is different in output prices and outputs")
    end
    if mw != m
        error("number of input prices and intputs is different")
    end
    if sp != s
        error("number of output prices and outputs is different")
    end
    if size(Gx) != size(X)
        error("size of inputs should be equal to size of inputs direction")
    end
    if size(Gy) != size(Y)
        error("size of outputs should be equal to size of outputs direction")
    end

    # Compute efficiency for each DMU
    n = nx

    Xefficient = zeros(n,m)
    Yefficient = zeros(n,s)
    pefficiency = zeros(n)
    plambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        w0 = W[i,:]
        p0 = P[i,:]

        # Create the optimization model
        deamodel = Model(with_optimizer(GLPK.Optimizer))
        @variable(deamodel, Xeff[1:m])
        @variable(deamodel, Yeff[1:s])
        @variable(deamodel, lambda[1:n] >= 0)

        @objective(deamodel, Max, (sum(p0[j] .* Yeff[j] for j in 1:s)) - (sum(w0[j] .* Xeff[j] for j in 1:m)))

        @constraint(deamodel, [j in 1:m], sum(X[t,j] * lambda[t] for t in 1:n) <= Xeff[j])
        @constraint(deamodel, [j in 1:s], sum(Y[t,j] * lambda[t] for t in 1:n) >= Yeff[j])

        @constraint(deamodel, sum(lambda) == 1)

        # Optimize and return results
        JuMP.optimize!(deamodel)

        Xefficient[i,:]  = JuMP.value.(Xeff)
        Yefficient[i,:]  = JuMP.value.(Yeff)
        plambdaeff[i,:] = JuMP.value.(lambda)

    end

    # Profit, technical and allocative efficiency
    maxprofit = sum(P .* Yefficient, dims = 2) .- sum(W .* Xefficient, dims = 2)

    pefficiency_num  = maxprofit .- ( sum(P .* Y, dims = 2) .- sum(W .* X, dims = 2))
    pefficiency_den = sum(P .* Gy, dims = 2) .+ sum(W .* Gx, dims = 2)
    pefficiency = vec( pefficiency_num ./ pefficiency_den )

    techefficiency = efficiency(deaddf(X, Y, Gx, Gy, rts = :VRS, slack = false))
    allocefficiency = pefficiency - techefficiency

    return ProfitDEAModel(n, m, s, pefficiency, plambdaeff, techefficiency, allocefficiency)

end

function deaprofit(X::Vector, Y::Matrix, W::Vector, P::Matrix, Gx::Vector, Gy::Matrix)::ProfitDEAModel
    X = X[:,:]
    W = W[:,:]
    Gx = Gx[:,:]
    return deaprofit(X, Y, W, P, Gx, Gy)
end

function deaprofit(X::Matrix, Y::Vector, W::Matrix, P::Vector, Gx::Matrix, Gy::Vector)::ProfitDEAModel
    Y = Y[:,:]
    P = P[:,:]
    Gy = Gy[:,:]
    return deaprofit(X, Y, W, P, Gx, Gy)
end

function deaprofit(X::Vector, Y::Vector, W::Vector, P::Vector, Gx::Vector, Gy::Vector)::ProfitDEAModel
    X = X[:,:]
    Y = Y[:,:]
    W = W[:,:]
    P = P[:,:]
    Gx = Gx[:,:]
    Gy = Gy[:,:]
    return deaprofit(X, Y, W, P, Gx, Gy)
end

function Base.show(io::IO, x::ProfitDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Profit DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = VRS")
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Profit", "Technical", "Allocative"], ["$i" for i in 1:n]))

    else

    end
end


text:
function prgp(cfg, X_train, y_train, X_test, y_test; rec=Record(), seed=1, verbose=true)
    Random.seed!(seed)
    init_id()
    init_ref_counter(cfg)
    init_entry_pool()

    data_train = Dataset(X_train, y_train)
    data_test  = Dataset(X_test, y_test)

    init_entries = cfg.create_population(cfg, data_train, data_test)
    add_entry!.(init_entries)

    pop = getfield.(init_entries, :index)
    new_pop = similar(pop)
    old_pop = similar(pop)

    avoid_overwrite.(pop)

    fitness = getfield.(init_entries, :train_fitness)

    register!(rec, "train_fitness", train_fitness, Float64)
    register!(rec, "test_fitness", test_fitness, Float64)

    train_rec = rec
    add_record!(rec, 1, pop, cfg)

    overwrite = false
    verbose && print_latest_fitness(rec, 1)

    for g in 2:cfg.num_generations
        overwrite = g > 2

        elite_idx = argmin(fitness)
        new_pop[1] = pop[elite_idx]
        avoid_overwrite(pop[elite_idx])
        for i in 2:cfg.population_size
            rand_num = rand()

            if rand_num < cfg.p_crossover
                p1 = tournament_selection(pop, fitness, cfg)
                p2 = tournament_selection(pop, fitness, cfg)
                offspring = projection_crossover(p1, p2, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            else
                p = tournament_selection(pop, fitness, cfg)
                offspring = projection_mutation(p, data_train, data_test, cfg, overwrite)
                new_pop[i] = offspring.index
            end
            # TODO: reproduction
        end

        update_fitness!(fitness, new_pop)
        add_record!(rec, g, pop, cfg)

        g > 2 && update_entry_pool!(old_pop)

        old_pop = copy(pop)
        pop = copy(new_pop)

        verbose && print_latest_fitness(rec, g)
    end

    return pop[argmin(fitness)], rec
end

optimal_coeff(v1, v2, t) = pinv([v1 v2])*t

function projection_crossover(t1, t2, data_train, data_test, cfg, overwrite)
    increment_ref(entry_pool[t1])
    increment_ref(entry_pool[t2])

    local r, offspring
    if overwrite
        offspring = projection_crossover_entry(entry_pool[t1], entry_pool[t2], data_train, data_test, cfg)
    else
        offspring = projection_crossover_entry(entry_pool[t1], entry_pool[t2], data_train, data_test, cfg, index=get_counter_size()+1)
    end

    add_entry!(offspring)
    set_ref(offspring)

    # to avoid overwritten
    # need to decrement at the end of current generation
    increment_ref(offspring)

    return offspring
end

function projection_crossover_entry(t1, t2, data_train, data_test, cfg; index=get_usable_index())

    r = optimal_coeff(t1.train_semantics, t2.train_semantics, data_train.target)

    train_semantics = Vector{Float64}(undef, length(data_train.target))
    @inbounds @simd for i in 1:length(train_semantics)
        train_semantics[i] = r[1] * t1.train_semantics[i] + r[2] * t2.train_semantics[i]
    end

    test_semantics = Vector{Float64}(undef, length(data_test.target))
    @inbounds @simd for i in 1:length(test_semantics)
        test_semantics[i] = r[1] * t1.test_semantics[i] + r[2] * t2.test_semantics[i]
    end

    train_fitness = cfg.fitness(train_semantics, data_train)
    test_fitness = cfg.fitness(test_semantics, data_test)

    id = get_new_id()

    return Entry(id, index, crossover, dummy_node(), (t1.id, t2.id, -1), (t1.index, t2.index, -1), train_semantics, test_semantics, train_fitness, test_fitness)
end

function projection_mutation(t, data_train, data_test, cfg, overwrite)
    increment_ref(entry_pool[t])

    local offspring, r1, r2
    if overwrite
        r1 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg)
        add_entry!(r1)
        set_ref(r1)
        increment_ref(r1)

        r2 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg)
        add_entry!(r2)
        set_ref(r2)
        increment_ref(r2)

        offspring = projection_mutation_entry(entry_pool[t], r1, r2,data_train, data_test, cfg)
    else
        r1 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg, index=get_counter_size()+1)
        add_entry!(r1)
        set_ref(r1)
        increment_ref(r1)

        r2 = Entry(cfg.create_rand_tree(cfg.pset, cfg.min_depth, cfg.max_depth), data_train, data_test, cfg, index=get_counter_size()+1)
        add_entry!(r2)
        set_ref(r2)
        increment_ref(r2)

        offspring = projection_mutation_entry(entry_pool[t], r1, r2, data_train, data_test, cfg, index=get_counter_size()+1)
    end

    add_entry!(offspring)
    set_ref(offspring)

    # to avoid overwritten
    # need to decrement at the end of current generation
    increment_ref(offspring)

    return offspring
end

function projection_mutation_entry(t, r1, r2, data_train, data_test, cfg;index=get_usable_index())

    coeff = optimal_coeff(t.train_semantics, r1.train_semantics-r2.train_semantics, data_train.target)

    train_semantics = Vector{Float64}(undef, length(data_train.target))
    @inbounds @simd for i in 1:length(train_semantics)
        train_semantics[i] = coeff[1]*t.train_semantics[i] + coeff[2]*(r1.train_semantics[i] - r2.train_semantics[i])
    end

    test_semantics = Vector{Float64}(undef, length(data_test.target))
    @inbounds @simd for i in 1:length(test_semantics)
        test_semantics[i] = coeff[1]*t.test_semantics[i] + coeff[2]*(r1.test_semantics[i] - r2.test_semantics[i])
    end

    train_fitness = cfg.fitness(train_semantics, data_train)
    test_fitness = cfg.fitness(test_semantics, data_test)

    id = get_new_id()

    return Entry(id, index, mutation, dummy_node(), (t.id, r1.id, r2.id), (t.index, r1.index, r2.index), train_semantics, test_semantics, train_fitness, test_fitness)
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
using IncrementalInference
using Test

@testset "Test ccolamd for constrained variable ordering" begin

fg = generateCanonicalFG_Kaess(graphinit=false)

vo = getEliminationOrder(fg, constraints=[:x3], ordering=:ccolamd)

@test vo[end] == :x3
@test length(vo) == length(ls(fg))

vo = getEliminationOrder(fg, constraints=[:l2], ordering=:ccolamd)

@test vo[end] == :l2

vo = getEliminationOrder(fg, constraints=[:x3;:l2], ordering=:ccolamd)

@test intersect(vo[end-1:end], [:x3;:l2]) |> length == 2

end


text:
cpyop(f::Function, x) = begin
    xo = C.PyObject_From(x)
    isnull(xo) && pythrow()
    r = f(xo)
    C.Py_DecRef(xo)
    r
end

cpyop(f::Function, x, y) = begin
    xo = C.PyObject_From(x)
    isnull(xo) && pythrow()
    yo = C.PyObject_From(y)
    isnull(yo) && (C.Py_DecRef(xo); pythrow())
    r = f(xo, yo)
    C.Py_DecRef(xo)
    C.Py_DecRef(yo)
    r
end

cpyop(f::Function, x, y, z) = begin
    xo = C.PyObject_From(x)
    isnull(xo) && pythrow()
    yo = C.PyObject_From(y)
    isnull(yo) && (C.Py_DecRef(xo); pythrow())
    zo = C.PyObject_From(z)
    isnull(yo) && (C.Py_DecRef(xo); C.Py_DecRef(yo); pythrow())
    r = f(xo, yo, zo)
    C.Py_DecRef(xo)
    C.Py_DecRef(yo)
    C.Py_DecRef(zo)
    r
end

cpyop(::Type{T}, f::Function, args...) where {T} = checknullconvert(T, cpyop(f, args...))

"""
    pyconvert(T, x) :: T

Convert Python object `x` to a `T`.
"""
pyconvert(::Type{T}, x) where {T} = checknullconvert(T, C.PyObject_From(x))
export pyconvert

"""
    pyis(x, y) :: Bool

Equivalent to `x is y` in Python.
"""
pyis(x::X, y::Y) where {X,Y} = begin
    if ispyreftype(X)
        xo = pyptr(x)
        isnull(xo) && pythrow()
    else
        xo = C.PyObject_From(x)
        isnull(xo) && pythrow()
        C.Py_DecRef(xo)
    end
    if ispyreftype(Y)
        yo = pyptr(y)
        isnull(yo) && pythrow()
    else
        yo = C.PyObject_From(y)
        isnull(yo) && pythrow()
        C.Py_DecRef(yo)
    end
    xo == yo
end
export pyis

"""
    pyhasattr(x, k) :: Bool

Equivalent to `hasattr(x, k)` in Python, returned as a `Bool`.
"""
pyhasattr(x, k) = checkm1(cpyop(C.PyObject_HasAttr, x, k)) != 0
pyhasattr(x, k::String) = checkm1(cpyop(xo -> C.PyObject_HasAttrString(xo, k), x)) != 0
pyhasattr(x, k::Symbol) = pyhasattr(x, string(k))
export pyhasattr

"""
    pygetattr([T=PyObject,] x, k) :: T

Equivalent to `x.k` or `getattr(x, k)` in Python.
"""
pygetattr(::Type{T}, x, k) where {T} = cpyop(T, C.PyObject_GetAttr, x, k)
pygetattr(::Type{T}, x, k::String) where {T} =
    cpyop(T, xo -> C.PyObject_GetAttrString(xo, k), x)
pygetattr(::Type{T}, x, k::Symbol) where {T} = pygetattr(T, x, string(k))
pygetattr(x, k) = pygetattr(PyObject, x, k)
export pygetattr

"""
    pysetattr(x, k, v)

Equivalent to `x.k = v` or `setattr(x, k, v)` in Python, but returns `x`.
"""
pysetattr(x, k, v) = (checkm1(cpyop(C.PyObject_SetAttr, x, k, v)); x)
pysetattr(x, k::String, v) =
    (checkm1(cpyop((xo, vo) -> C.PyObject_SetAttrString(xo, k, vo), x, v)); x)
pysetattr(x, k::Symbol, v) = pysetattr(x, string(k), v)
export pysetattr

"""
    pydir([T=PyObject,] x) :: T

Equivalent to `dir(x)` in Python.
"""
pydir(::Type{T}, x) where {T} = cpyop(T, C.PyObject_Dir, x)
pydir(x) = pydir(PyObject, x)
export pydir

"""
    pycall([T=PyObject,] f, args...; kwargs...) :: T

Equivalent to `f(*args, **kwargs)` in Python.
"""
pycall(::Type{T}, f, args...; opts...) where {T} =
    cpyop(T, fo -> C.PyObject_CallArgs(fo, args, opts), f)
pycall(f, args...; opts...) = pycall(PyObject, f, args...; opts...)
export pycall

"""
    pyrepr([T=PyObject,] x) :: T

Equivalent to `repr(x)` in Python.
"""
pyrepr(::Type{T}, x) where {T} = cpyop(T, C.PyObject_Repr, x)
pyrepr(x) = pyrepr(PyObject, x)
export pyrepr

"""
    pystr([T=PyObject,] x) :: T

Equivalent to `str(x)` in Python.
"""
pystr(::Type{T}, x) where {T} = cpyop(T, C.PyObject_Str, x)
pystr(::Type{T}, x::Union{String,SubString{String},Vector{Int8},Vector{UInt8}}) where {T} =
    checknullconvert(T, C.PyUnicode_From(x))
pystr(x) = pystr(PyObject, x)
export pystr

"""
    pybytes([T=PyObject,] x) :: T

Equivalent to `str(x)` in Python.
"""
pybytes(::Type{T}, x) where {T} = cpyop(T, C.PyObject_Bytes, x)
pybytes(
    ::Type{T},
    x::Union{Vector{Int8},Vector{UInt8},String,SubString{String}},
) where {T} = checknullconvert(T, C.PyBytes_From(x))
pybytes(x) = pybytes(PyObject, x)
export pybytes

"""
    pylen(x) :: Integer

Equivalent to `len(x)` in Python.
"""
pylen(x) = checkm1(cpyop(C.PyObject_Length, x))
export pylen

"""
    pycontains(x, v) :: Bool

Equivalent to `v in x` in Python.
"""
pycontains(x, v) = checkm1(cpyop(C.PySequence_Contains, x, v)) != 0
export pycontains

"""
    pyin(v, x) :: Bool

Equivalent to `v in x` in Python.
"""
pyin(v, x) = pycontains(x, v)
export pyin

"""
    pygetitem([T=PyObject,] x, k) :: T

Equivalent to `x[k]` or `getitem(x, k)` in Python.
"""
pygetitem(::Type{T}, x, k) where {T} = cpyop(T, C.PyObject_GetItem, x, k)
pygetitem(x, k) = pygetitem(PyObject, x, k)
export pygetitem

"""
    pysetitem(x, k, v)

Equivalent to `x[k] = v` or `setitem(x, k, v)` in Python, but returns `x`.
"""
pysetitem(x, k, v) = (checkm1(cpyop(C.PyObject_SetItem, x, k, v)); x)
export pysetitem

"""
    pydelitem(x, k)

Equivalent to `del x[k]` or `delitem(x, k)` in Python, but returns x.
"""
pydelitem(x, k) = (checkm1(cpyop(C.PyObject_DelItem, x, k)); x)
export pydelitem

"""
    pynone([T=PyObject]) :: T

Equivalent to `None` in Python.
"""
pynone(::Type{T}) where {T} = (checkm1(C.PyObject_Convert(C.Py_None(), T)); takeresult(T))
pynone() = pynone(PyObject)
export pynone

"""
    pybool([T=PyObject,] ...) :: T

Equivalent to `bool(...)` in Python.
"""
pybool(::Type{T}, x::Bool) where {T} = checknullconvert(T, C.PyBool_From(x))
pybool(::Type{T}, args...; kwargs...) where {T} =
    checknullconvert(T, C.PyObject_CallArgs(C.PyBool_Type(), args, kwargs))
pybool(args...; kwargs...) = pybool(PyObject, args...; kwargs...)
export pybool

"""
    pyint([T=PyObject,] ...) :: T

Equivalent to `int(...)` in Python.
"""
pyint(
    ::Type{T},
    x::Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128,BigInt},
) where {T} = checknullconvert(T, C.PyLong_From(x))
pyint(::Type{T}, x) where {T} = cpyop(T, C.PyNumber_Long, x)
pyint(::Type{T}, args...; kwargs...) where {T} =
    checknullconvert(T, C.PyObject_CallArgs(C.PyLong_Type(), args, kwargs))
pyint(args...; kwargs...) = pyint(PyObject, args...; kwargs...)
export pyint

"""
    pyfloat([T=PyObject,] ...) :: T

Equivalent to `float(...)` in Python.
"""
pyfloat(::Type{T}, x::Union{Float16,Float32,Float64}) where {T} =
    checknullconvert(T, C.PyFloat_From(x))
pyfloat(::Type{T}, x) where {T} = cpyop(T, C.PyNumber_Float, x)
pyfloat(::Type{T}, args...; kwargs...) where {T} =
    checknullconvert(T, C.PyObject_CallArgs(C.PyFloat_Type(), args, kwargs))
pyfloat(args...; kwargs...) = pyfloat(PyObject, args...; kwargs...)
export pyfloat

"""
    pyimport([T=PyObject,] name) :: T
    pyimport([T=PyObject,] name=>attr) :: T
    pyimport([T=PyObject,] name=>(attr,...)) :: Tuple{T,...}

Imports and returns the Python module `name`.

If additionally `attr` is given, the given attribute of the module is returned instead.
It may also be a tuple of attributes.

If several arguments are given, each one is imported and a tuple is returned.
"""
pyimport(::Type{T}, m) where {T} = cpyop(T, C.PyImport_Import, m)
pyimport(::Type{T}, m::String) where {T} = checknullconvert(T, C.PyImport_ImportModule(m))
pyimport(::Type{T}, m::Pair) where {T} = begin
    r = pyimport(PyObject, m[1])
    m[2] isa Tuple ? map(k -> pygetattr(T, r, k), m[2]) : pygetattr(T, r, m[2])
end
pyimport(::Type{T}, m1, m2, ms...) where {T} = map(m -> pyimport(T, m), (m1, m2, ms...))
pyimport(m1, ms...) = pyimport(PyObject, m1, ms...)
export pyimport

"""
    pytruth(x) :: Bool

The truthyness of `x`, equivalent to `bool(x)` or `not not x` in Python, or to `pybool(Bool, x)`.
"""
pytruth(x) = checkm1(cpyop(C.PyObject_IsTrue, x)) != 0
export pytruth

"""
    pyissubclass(x, y) :: Bool

Equivalent to `issubclass(x, y)` in Python.
"""
pyissubclass(x, y) = checkm1(cpyop(C.PyObject_IsSubclass, x, y)) != 0
export pyissubclass

"""
    pyisinstance(x, y) :: Bool

Equivalent to `isinstance(x, y)` in Python.
"""
pyisinstance(x, y) = checkm1(cpyop(C.PyObject_IsInstance, x, y)) != 0
export pyisinstance

"""
    pyhash(x) :: Integer

Equivalent to `hash(x)` in Python.
"""
pyhash(x) = checkm1(cpyop(C.PyObject_Hash, x))
export pyhash

"""
    pycompare([T=PyObject,] x, op, y) :: T

Equivalent to `x op y` in Python where `op` is one of: `=`, `≠`, `<`, `≤`, `>`, `≥`.
"""
pycompare(::Type{T}, x, ::typeof(==), y) where {T} = _pycompare(T, x, C.Py_EQ, y)
pycompare(::Type{T}, x, ::typeof(!=), y) where {T} = _pycompare(T, x, C.Py_NE, y)
pycompare(::Type{T}, x, ::typeof(<), y) where {T} = _pycompare(T, x, C.Py_LT, y)
pycompare(::Type{T}, x, ::typeof(<=), y) where {T} = _pycompare(T, x, C.Py_LE, y)
pycompare(::Type{T}, x, ::typeof(>), y) where {T} = _pycompare(T, x, C.Py_GT, y)
pycompare(::Type{T}, x, ::typeof(>=), y) where {T} = _pycompare(T, x, C.Py_GE, y)
pycompare(x, op, y) = pycompare(PyObject, x, op, y)
_pycompare(::Type{T}, x, op::Cint, y) where {T} =
    cpyop(T, (xo, yo) -> C.PyObject_RichCompare(xo, yo, op), x, y)
_pycompare(::Type{Bool}, x, op::Cint, y) =
    checkm1(cpyop((xo, yo) -> C.PyObject_RichCompareBool(xo, yo, op), x, y)) != 0
export pycompare

"""
    pyeq([T=PyObject,] x, y) :: T

Equivalent to `x == y` in Python.
"""
pyeq(::Type{T}, x, y) where {T} = pycompare(T, x, ==, y)
pyeq(x, y) = pycompare(x, ==, y)
export pyeq

"""
    pyne([T=PyObject,] x, y) :: T

Equivalent to `x != y` in Python.
"""
pyne(::Type{T}, x, y) where {T} = pycompare(T, x, !=, y)
pyne(x, y) = pycompare(x, !=, y)
export pyne

"""
    pyge([T=PyObject,] x, y) :: T

Equivalent to `x >= y` in Python.
"""
pyge(::Type{T}, x, y) where {T} = pycompare(T, x, >=, y)
pyge(x, y) = pycompare(x, >=, y)
export pyge

"""
    pygt([T=PyObject,] x, y) :: T

Equivalent to `x > y` in Python.
"""
pygt(::Type{T}, x, y) where {T} = pycompare(T, x, >, y)
pygt(x, y) = pycompare(x, >, y)
export pygt

"""
    pyle([T=PyObject,] x, y) :: T

Equivalent to `x <= y` in Python.
"""
pyle(::Type{T}, x, y) where {T} = pycompare(T, x, <=, y)
pyle(x, y) = pycompare(x, <=, y)
export pyle

"""
    pylt([T=PyObject,] x, y) :: T

Equivalent to `x < y` in Python.
"""
pylt(::Type{T}, x, y) where {T} = pycompare(T, x, <, y)
pylt(x, y) = pycompare(x, <, y)
export pylt

"""
    pyadd([T=PyObject,] x, y) :: T

Equivalent to `x + y` in Python.
"""
pyadd(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Add, x, y)
pyadd(x, y) = pyadd(PyObject, x, y)
export pyadd

"""
    pyiadd([T=typeof(x),] x, y) :: T

`x = pyiadd(x, y)` is equivalent to `x += y` in Python.
"""
pyiadd(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceAdd, x, y)
pyiadd(x, y) = pyiadd(typeof(x), x, y)
export pyiadd

"""
    pysub([T=PyObject,] x, y) :: T

Equivalent to `x - y` in Python.
"""
pysub(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Subtract, x, y)
pysub(x, y) = pysub(PyObject, x, y)
export pysub

"""
    pyisub([T=typeof(x),] x, y) :: T

`x = pyisub(x, y)` is equivalent to `x -= y` in Python.
"""
pyisub(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceSubtract, x, y)
pyisub(x, y) = pyisub(typeof(x), x, y)
export pyisub

"""
    pymul([T=PyObject,] x, y) :: T

Equivalent to `x * y` in Python.
"""
pymul(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Multiply, x, y)
pymul(x, y) = pymul(PyObject, x, y)
export pymul

"""
    pyimul([T=typeof(x),] x, y) :: T

`x = pyimul(x, y)` is equivalent to `x *= y` in Python.
"""
pyimul(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceMultiply, x, y)
pyimul(x, y) = pyimul(typeof(x), x, y)
export pyimul

"""
    pymatmul([T=PyObject,] x, y) :: T

Equivalent to `x @ y` in Python.
"""
pymatmul(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_MatrixMultiply, x, y)
pymatmul(x, y) = pymatmul(PyObject, x, y)
export pymatmul

"""
    pyimatmul([T=typeof(x),] x, y) :: T

`x = pyimatmul(x, y)` is equivalent to `x @= y` in Python.
"""
pyimatmul(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceMatrixMultiply, x, y)
pyimatmul(x, y) = pyimatmul(typeof(x), x, y)
export pyimatmul

"""
    pyfloordiv([T=PyObject,] x, y) :: T

Equivalent to `x // y` in Python.
"""
pyfloordiv(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_FloorDivide, x, y)
pyfloordiv(x, y) = pyfloordiv(PyObject, x, y)
export pyfloordiv

"""
    pyifloordiv([T=typeof(x),] x, y) :: T

`x = pyifloordiv(x, y)` is equivalent to `x //= y` in Python.
"""
pyifloordiv(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceFloorDivide, x, y)
pyifloordiv(x, y) = pyifloordiv(typeof(x), x, y)
export pyifloordiv

"""
    pytruediv([T=PyObject,] x, y) :: T

Equivalent to `x / y` in Python.
"""
pytruediv(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_TrueDivide, x, y)
pytruediv(x, y) = pytruediv(PyObject, x, y)
export pytruediv

"""
    pyitruediv([T=typeof(x),] x, y) :: T

`x = pyitruediv(x, y)` is equivalent to `x /= y` in Python.
"""
pyitruediv(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceTrueDivide, x, y)
pyitruediv(x, y) = pyitruediv(typeof(x), x, y)
export pyitruediv

"""
    pymod([T=PyObject,] x, y) :: T

Equivalent to `x % y` in Python.
"""
pymod(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Remainder, x, y)
pymod(x, y) = pymod(PyObject, x, y)
export pymod

"""
    pyimod([T=typeof(x),] x, y) :: T

`x = pyimod(x, y)` is equivalent to `x %= y` in Python.
"""
pyimod(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceRemainder, x, y)
pyimod(x, y) = pyimod(typeof(x), x, y)
export pyimod

"""
    pydivmod([T=PyObject,] x, y) :: T

Equivalent to `divmod(x, y)` in Python.
"""
pydivmod(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Divmod, x, y)
pydivmod(x, y) = pydivmod(PyObject, x, y)
export pydivmod

"""
    pylshift([T=PyObject,] x, y) :: T

Equivalent to `x << y` in Python.
"""
pylshift(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Lshift, x, y)
pylshift(x, y) = pylshift(PyObject, x, y)
export pylshift

"""
    pyilshift([T=typeof(x),] x, y) :: T

`x = pyilshift(x, y)` is equivalent to `x <<= y` in Python.
"""
pyilshift(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceLshift, x, y)
pyilshift(x, y) = pyilshift(typeof(x), x, y)
export pyilshift

"""
    pyrshift([T=PyObject,] x, y) :: T

Equivalent to `x >> y` in Python.
"""
pyrshift(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Rshift, x, y)
pyrshift(x, y) = pyrshift(PyObject, x, y)
export pyrshift

"""
    pyirshift([T=typeof(x),] x, y) :: T

`x = pyirshift(x, y)` is equivalent to `x >>= y` in Python.
"""
pyirshift(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceRshift, x, y)
pyirshift(x, y) = pyirshift(typeof(x), x, y)
export pyirshift

"""
    pyand([T=PyObject,] x, y) :: T

Equivalent to `x & y` in Python.
"""
pyand(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_And, x, y)
pyand(x, y) = pyand(PyObject, x, y)
export pyand

"""
    pyiand([T=typeof(x),] x, y) :: T

`x = pyiand(x, y)` is equivalent to `x &= y` in Python.
"""
pyiand(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceAnd, x, y)
pyiand(x, y) = pyiand(typeof(x), x, y)
export pyiand

"""
    pyxor([T=PyObject,] x, y) :: T

Equivalent to `x ^ y` in Python.
"""
pyxor(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Xor, x, y)
pyxor(x, y) = pyxor(PyObject, x, y)
export pyxor

"""
    pyixor([T=typeof(x),] x, y) :: T

`x = pyixor(x, y)` is equivalent to `x ^= y` in Python.
"""
pyixor(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceXor, x, y)
pyixor(x, y) = pyixor(typeof(x), x, y)
export pyixor

"""
    pyor([T=PyObject,] x, y) :: T

Equivalent to `x | y` in Python.
"""
pyor(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_Or, x, y)
pyor(x, y) = pyor(PyObject, x, y)
export pyor

"""
    pyior([T=typeof(x),] x, y) :: T

`x = pyior(x, y)` is equivalent to `x |= y` in Python.
"""
pyior(::Type{T}, x, y) where {T} = cpyop(T, C.PyNumber_InPlaceOr, x, y)
pyior(x, y) = pyior(typeof(x), x, y)
export pyior

"""
    pypow([T=PyObject,] x, y, [z]) :: T

Equivalent to `x**y` or `pow(x, y, z)` in Python.
"""
pypow(::Type{T}, x, y, z = C.PyObjectRef(C.Py_None())) where {T} =
    cpyop(T, C.PyNumber_Power, x, y, z)
pypow(x, y, z) = pypow(PyObject, x, y, z)
pypow(x, y) = pypow(PyObject, x, y)
export pypow

"""
    pyipow([T=typeof(x),] x, y, [z]) :: T

`x = pyipow(x, y)` is equivalent to `x **= y` in Python.
"""
pyipow(::Type{T}, x, y, z = C.PyObjectRef(C.Py_None())) where {T} =
    cpyop(T, C.PyNumber_InPlacePower, x, y, z)
pyipow(x, y, z) = pyipow(typeof(x), x, y, z)
pyipow(x, y) = pyipow(typeof(x), x, y)
export pyipow

"""
    pyneg([T=typeof(x),] x) :: T

Equivalent to `-x` in Python.
"""
pyneg(::Type{T}, x) where {T} = cpyop(T, C.PyNumber_Negative, x)
pyneg(x) = pyneg(typeof(x), x)
export pyneg

"""
    pypos([T=typeof(x),] x) :: T

Equivalent to `+x` in Python.
"""
pypos(::Type{T}, x) where {T} = cpyop(T, C.PyNumber_Positive, x)
pypos(x) = pypos(typeof(x), x)
export pypos

"""
    pyabs([T=typeof(x),] x) :: T

Equivalent to `abs(x)` in Python.
"""
pyabs(::Type{T}, x) where {T} = cpyop(T, C.PyNumber_Absolute, x)
pyabs(x) = pyabs(typeof(x), x)
export pyabs

"""
    pyinv([T=typeof(x),] x) :: T

Equivalent to `-x` in Python.
"""
pyinv(::Type{T}, x) where {T} = cpyop(T, C.PyNumber_Invert, x)
pyinv(x) = pyinv(typeof(x), x)
export pyinv

"""
    pyiter([T=PyObject] x) :: T

Equivalent to `iter(x)` in Python.
"""
pyiter(::Type{T}, x) where {T} = cpyop(T, C.PyObject_GetIter, x)
pyiter(x) = pyiter(PyObject, x)
export pyiter

"""
    pywith(f, o, d=nothing)

Equivalent to `with o as x: f(x)` in Python, where `x` is a `PyObject`.

On success, the value of `f(x)` is returned.
If an exception occurs but is suppressed then `d` is returned.
"""
function pywith(f, _o, d = nothing)
    o = PyObject(_o)
    t = pytype(o)
    exit = t.__exit__
    value = t.__enter__(o)
    exited = false
    try
        return f(value)
    catch err
        if err isa PyException
            exited = true
            if pytruth(exit(o, err.tref, err.vref, err.bref))
                rethrow()
            else
                return d
            end
        else
            rethrow()
        end
    finally
        exited || exit(o, pynone(), pynone(), pynone())
    end
end
export pywith

"""
    pytuple([T=PyObject,] [x]) :: T

Create a Python `tuple` from the elements of iterable `x`.

If `x` is a Python object, this is equivalent to `tuple(x)` in Python.
"""
pytuple(::Type{T}, x) where {T} = checknullconvert(
    T,
    ispyref(x) ? C.PyObject_CallNice(C.PyTuple_Type(), x) : C.PyTuple_FromIter(x),
)
pytuple(::Type{T}) where {T} = checknullconvert(T, C.PyTuple_New(0))
pytuple(x) = pytuple(PyObject, x)
pytuple() = pytuple(PyObject)
export pytuple

"""
    pylist([T=PyObject,] [x]) :: T

Create a Python `list` from the elements of iterable `x`.

If `x` is a Python object, this is equivalent to `list(x)` in Python.
"""
pylist(::Type{T}, x) where {T} = checknullconvert(
    T,
    ispyref(x) ? C.PyObject_CallNice(C.PyList_Type(), x) : C.PyList_FromIter(x),
)
pylist(::Type{T}) where {T} = checknullconvert(T, C.PyList_New(0))
pylist(x) = pylist(PyObject, x)
pylist() = pylist(PyObject)
export pylist

"""
    pycollist([T=PyObject,] x::AbstractArray) :: T

Create a nested Python `list`-of-`list`s from the elements of `x`. For matrices, this is a list of columns.
"""
pycollist(::Type{T}, x::AbstractArray) where {T} =
    ndims(x) == 0 ? pyconvert(T, x[]) :
    pylist(T, pycollist(PyRef, y) for y in eachslice(x; dims = ndims(x)))
pycollist(x::AbstractArray) = pycollist(PyObject, x)
export pycollist

"""
    pyrowlist([T=PyObject,] x::AbstractArray) :: T

Create a nested Python `list`-of-`list`s from the elements of `x`. For matrices, this is a list of rows.
"""
pyrowlist(::Type{T}, x::AbstractArray) where {T} =
    ndims(x) == 0 ? pyconvert(T, x[]) :
    pylist(T, pyrowlist(PyRef, y) for y in eachslice(x; dims = 1))
pyrowlist(x::AbstractArray) = pyrowlist(PyObject, x)
export pyrowlist

"""
    pyset([T=PyObject,] [x]) :: T

Create a Python `set` from the elements of iterable `x`.

If `x` is a Python object, this is equivalent to `set(x)` in Python.
"""
pyset(::Type{T}, x) where {T} = checknullconvert(
    T,
    ispyref(x) ? C.PyObject_CallNice(C.PySet_Type(), x) : C.PySet_FromIter(x),
)
pyset(::Type{T}) where {T} = checknullconvert(T, C.PySet_New(C_NULL))
pyset(x) = pyset(PyObject, x)
pyset() = pyset(PyObject)
export pyset

"""
    pyfrozenset([T=PyObject,] [x]) :: T

Create a Python `frozenset` from the elements of iterable `x`.

If `x` is a Python object, this is equivalent to `frozenset(x)` in Python.
"""
pyfrozenset(::Type{T}, x) where {T} = checknullconvert(
    T,
    ispyref(x) ? C.PyObject_CallNice(C.PyFrozenSet_Type(), x) : C.PyFrozenSet_FromIter(x),
)
pyfrozenset(::Type{T}) where {T} = checknullconvert(T, C.PyFrozenSet_New(C_NULL))
pyfrozenset(x) = pyfrozenset(PyObject, x)
pyfrozenset() = pyfrozenset(PyObject)
export pyfrozenset

"""
    pydict([T=PyObject,] [x]) :: T
    pydict([T=PyObject;] key=value, ...)

Create a Python `dict` from the given key-value pairs in `x` or keyword arguments.

If `x` is a Python object, this is equivalent to `dict(x)` in Python.
"""
pydict(::Type{T}, x) where {T} = checknullconvert(
    T,
    ispyref(x) ? C.PyObject_CallNice(C.PyDict_Type(), x) : C.PyDict_FromPairs(x),
)
pydict(::Type{T}; opts...) where {T} =
    checknullconvert(T, isempty(opts) ? C.PyDict_New() : C.PyDict_FromStringPairs(opts))
pydict(x) = pydict(PyObject, x)
pydict(; opts...) = pydict(PyObject; opts...)
export pydict

"""
    pyslice([T=PyObject,] [start,] stop, [step]) :: T

Equivalent to `slice(start, stop, step)` in Python (or `start:stop:step` while indexing).
"""
pyslice(::Type{T}, x) where {T} = cpyop(T, x -> C.PySlice_New(C_NULL, x, C_NULL), x)
pyslice(::Type{T}, x, y) where {T} = cpyop(T, (x, y) -> C.PySlice_New(x, y, C_NULL), x, y)
pyslice(::Type{T}, x, y, z) where {T} = cpyop(T, C.PySlice_New, x, y, z)
pyslice(x) = pyslice(PyObject, x)
pyslice(x, y) = pyslice(PyObject, x, y)
pyslice(x, y, z) = pyslice(PyObject, x, y, z)
export pyslice

"""
    pyellipsis([T=PyObject]) :: T

Equivalent to `Ellipsis` in Python (or `...` while indexing).
"""
pyellipsis(::Type{T}) where {T} = checknullconvert(T, C.PyEllipsis_New())
pyellipsis() = pyellipsis(PyObject)
export pyellipsis

"""
    pynotimplemented([T=PyObject]) :: T

Equivalent to `NotImplemented` in Python.
"""
pynotimplemented(::Type{T}) where {T} = checknullconvert(T, C.PyNotImplemented_New())
pynotimplemented() = pynotimplemented(PyObject)
export pynotimplemented

"""
    pymethod([T=PyObject,] x) :: T

Convert `x` to a Python instance method.
"""
pymethod(::Type{T}, x) where {T} = cpyop(T, C.PyInstanceMethod_New, x)
pymethod(x) = pymethod(PyObject, x)
export pymethod

"""
    pytype([T=PyObject,] x) :: T

Equivalent to `type(x)` in Python.
"""
pytype(::Type{T}, x) where {T} = cpyop(T, o -> (t = C.Py_Type(o); C.Py_IncRef(t); t), x)
pytype(x) = pytype(PyObject, x)
export pytype

"""
    pytype([T=PyObject,] name, bases, dict) :: T

Equivalent to `type(name, bases, dict)` in Python.
"""
pytype(::Type{T}, name, bases, dict) where {T} =
    @pyv `type($name, $(pytuple(bases)), $(dict isa NamedTuple ? pydict(;dict...) : pydict(dict)))`::T
pytype(name, bases, dict) = pytype(PyObject, name, bases, dict)
export pytype

### MULTIMEDIA DISPLAY

_py_mime_reprmethod(::MIME) = nothing
_py_mime_reprmethod(::MIME"text/plain") = "__repr__"
_py_mime_reprmethod(::MIME"text/html") = "_repr_html_"
_py_mime_reprmethod(::MIME"text/markdown") = "_repr_markdown_"
_py_mime_reprmethod(::MIME"text/json") = "_repr_json_"
_py_mime_reprmethod(::MIME"text/latex") = "_repr_latex_"
_py_mime_reprmethod(::MIME"application/javascript") = "_repr_javascript_"
_py_mime_reprmethod(::MIME"application/pdf") = "_repr_pdf_"
_py_mime_reprmethod(::MIME"image/jpeg") = "_repr_jpeg_"
_py_mime_reprmethod(::MIME"image/png") = "_repr_png_"
_py_mime_reprmethod(::MIME"image/svg+xml") = "_repr_svg_"

_py_mime_data(m::MIME, o) = begin
    @py ```
    o = $o
    m = $(string(m))
    r = $(_py_mime_reprmethod(m))
    data = None
    meta = None
    try:
        x = type(o)._repr_mimebundle_(o, include=[m])
        if isinstance(x, tuple):
            data = x[0][m]
            meta = x[1].get(m)
        else:
            data = x[m]
    except:
        pass
    if data is None and r is not None:
        try:
            x = getattr(type(o), r)(o)
            if isinstance(x, tuple):
                data = x[0]
                meta = x[1]
            else:
                data = x
        except:
            pass
    $(data::Union{PyRef, Nothing}) = data
    $(meta::Union{PyRef, Nothing}) = meta
    ```
    data, meta
end

_py_mime_showable(m::MIME, o) = begin
    data, meta = _py_mime_data(m, o)
    data !== nothing
end

_py_mime_show(io::IO, m::MIME, o) = begin
    data, meta = _py_mime_data(m, o)
    write(io, istextmime(m) ? pyconvert(String, data) : pyconvert(Vector{UInt8}, data))
    nothing
end

### IO

"""
    pytextio([T=PyObject], io::IO) :: T

Convert `io` to a Python text IO stream, specifically a `juliaaa.TextIOValue`.
"""
pytextio(::Type{T}, io::IO) where {T} = checknullconvert(T, C.PyJuliaTextIOValue_New(io))
pytextio(io::IO) = pytextio(PyObject, io)
export pytextio

"""
    pybufferedio([T=PyObject], io::IO) :: T

Convert `io` to a Python buffered byte IO stream, specifically a `juliaaa.BufferedIOValue`.
"""
pybufferedio(::Type{T}, io::IO) where {T} =
    checknullconvert(T, C.PyJuliaBufferedIOValue_New(io))
pybufferedio(io::IO) = pybufferedio(PyObject, io)
export pybufferedio

"""
    pyrawio([T=PyObject], io::IO) :: T

Convert `io` to a Python raw (unbuffered byte) IO stream, specifically a `juliaaa.RawIOValue`.
"""
pyrawio(::Type{T}, io::IO) where {T} = checknullconvert(T, C.PyJuliaRawIOValue_New(io))
pyrawio(io::IO) = pyrawio(PyObject, io)
export pyrawio