text:
abstract type AbstractTens{order,dim,T<:Number} <: AbstractArray{T,order} end

@pure getorder(::AbstractTens{getorder,dim,T}) where {getorder,dim,T} = getorder
@pure getdim(::AbstractTens{order,dim,T}) where {order,dim,T} = dim
@pure Base.eltype(::Type{AbstractTens{order,dim,T}}) where {order,dim,T} = T

"""
    Tens{order,dim,T,A<:AbstractArray}

Tensor type of any order defined by
- a multidata of components (of any type heriting from `AbstractArray`, e.g. `Tensor` or `SymmetricTensor`)
- a basis of `AbstractBasis` type
- a tuple of variances (covariant `:cov` or contravariant `:cont`) of length equal to the `order` of the tensor

# Examples
```julia
julia> ℬ = Basis(Sym[1 0 0; 0 1 0; 0 1 1]) ;

julia> T = Tens(metric(ℬ,:cov),ℬ,(:cov,:cov))
Tens{2, 3, Sym, SymmetricTensor{2, 3, Sym, 6}}
# data: 3×3 SymmetricTensor{2, 3, Sym, 6}:
 1  0  0
 0  2  1
 0  1  1
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 1  0  0
 0  1  0
 0  1  1
# var: (:cov, :cov)

julia> components(T,(:cont,:cov),b)
3×3 Matrix{Sym}:
 1  0  0
 0  1  0
 0  0  1
```
"""
struct Tens{order,dim,T,A<:AbstractArray} <: AbstractTens{order,dim,T}
    data::A
    basis::Basis
    var::NTuple{order,Symbol}
    function Tens(
        data::AbstractArray{T,order},
        basis::Basis{dim},
        var::NTuple{order,Symbol} = ntuple(_ -> :cont, Val(order)),
    ) where {order,dim,T}
        newdata = tensor_or_array(data)
        new{order,dim,T,typeof(newdata)}(newdata, basis, var)
    end
    Tens(data::AbstractArray, basis::RotatedBasis, args...) = TensRotated(data, basis)
    Tens(data::AbstractArray, basis::OrthogonalBasis, args...) =
        TensOrthogonal(data, basis, args...)
    Tens(
        data::AbstractArray{T},
        basis::OrthonormalBasis = CanonicalBasis{size(data)[1],T}(),
        args...,
    ) where {T} = TensOrthonormal(data, basis)
    Tens(data::AbstractArray, var::NTuple, basis::AbstractBasis = CanonicalBasis{size(data,1),eltype(data)}()) = Tens(data, basis, var)
    Tens(data::AbstractArray{T,0}, args...) where {T} = T(data[1])
    Tens(data::T, args...) where {T} = data
end

proj_tens(sym::Symbol, A::AbstractArray) = proj_tens(Val(sym), A)

for order ∈ (2, 4)
    @eval function best_sym_tens(
        t::AbstractTens{$order,dim,T},
        args...;
        proj = (:ISO,),
        ε = 1.e-6,
    ) where {dim,T}
        basis = relevant_OrthonormalBasis(getbasis(t))
        newt = change_tens(t, basis)
        for sym ∈ proj
            (projt, d, drel) = proj_tens(sym, newt)
            if d == zero(T) || drel < ε
                return projt, d, drel, sym
            end
        end
        return newt, zero(T), zero(T), :ANISO
    end
end

struct TensRotated{order,dim,T,A<:AbstractArray} <: AbstractTens{order,dim,T}
    data::A
    basis::RotatedBasis
    function TensRotated(
        data::AbstractArray{T,order},
        basis::RotatedBasis{dim},
    ) where {order,dim,T}
        newdata = tensor_or_array(data)
        new{order,dim,T,typeof(newdata)}(newdata, basis)
    end
end

struct TensCanonical{order,dim,T,A<:AbstractArray} <: AbstractTens{order,dim,T}
    data::A
    function TensCanonical(data::AbstractArray{T,order}) where {order,T}
        newdata = tensor_or_array(data)
        new{order,size(data, 1),T,typeof(newdata)}(newdata)
    end
end

struct TensOrthogonal{order,dim,T,A<:AbstractArray} <: AbstractTens{order,dim,T}
    data::A
    basis::OrthogonalBasis
    var::NTuple{order,Symbol}
    function TensOrthogonal(
        data::AbstractArray{T,order},
        basis::OrthogonalBasis{dim},
        var::NTuple{order,Symbol} = ntuple(_ -> :cont, Val(order)),
    ) where {order,dim,T}
        newdata = tensor_or_array(data)
        new{order,dim,T,typeof(newdata)}(newdata, basis, var)
    end
end

const TensOrthonormal{order,dim,T,A} =
    Union{TensRotated{order,dim,T,A},TensCanonical{order,dim,T,A}}
const AllTensOrthogonal{order,dim,T,A} = Union{TensOrthonormal{order,dim,T,A}, TensOrthogonal{order,dim,T,A}}
const TensVar{order,dim,T,A} = Union{Tens{order,dim,T,A},TensOrthogonal{order,dim,T,A}}
const TensBasis{order,dim,T,A} =
    Union{Tens{order,dim,T,A},TensRotated{order,dim,T,A},TensOrthogonal{order,dim,T,A}}
const TensArray{order,dim,T,A} = Union{
    Tens{order,dim,T,A},
    TensRotated{order,dim,T,A},
    TensCanonical{order,dim,T,A},
    TensOrthogonal{order,dim,T,A},
}

TensOrthonormal(data::AbstractArray, basis::RotatedBasis) = TensRotated(data, basis)
TensOrthonormal(data::AbstractArray, ::CanonicalBasis) = TensCanonical(data)

TensOrthonormal(t::TensOrthonormal) = t
function TensOrthonormal(t::TensOrthogonal{order,dim,T}) where {order,dim,T}
    m = Array(getarray(t))
    ℬ = getbasis(t)
    onℬ = relevant_OrthonormalBasis(ℬ)
    Λ = Dict(:cov => inv.(ℬ.λ), :cont => ℬ.λ)
    for ind ∈ CartesianIndices(m)
        m[ind] *= prod([Λ[getvar(t, i)][ind[i]] for i ∈ 1:order])
    end
    return Tens(m, onℬ)
end

@inline Base.size(t::TensArray) = size(getarray(t))
@inline Base.getindex(t::TensArray, ind...) = getindex(getarray(t), ind...)
@pure datatype(::TensArray{order,dim,T,A}) where {order,dim,T,A} = A

@inline Base.zero(t::AbstractTens) = Tens(zero.(getarray(t)), getbasis(t), getvar(t))

# This function aims at storing the table of components in the `Tensor` type whenever possible
tensor_or_array(tab::AbstractArray{T,1}) where {T} = Vec{size(tab, 1)}(tab)
for order ∈ (2, 4)
    @eval function tensor_or_array(tab::AbstractArray{T,$order}) where {T}
        dim = size(tab, 1)
        newtab = Tensor{$order,dim}(tab)
        if Tensors.issymmetric(newtab)
            newtab = convert(SymmetricTensor{$order,dim}, newtab)
        end
        return newtab
    end
    @eval tensor_or_array(tab::Tensor{$order,dim}) where {dim} = Tensors.issymmetric(tab) ? convert(SymmetricTensor{$order,dim}, tab) : tab
end
tensor_or_array(tab::Tensors.AllTensors) = tab
tensor_or_array(tab::AbstractArray) = tab

##############################
# Utility/Accessor Functions #
##############################

getarray(t::TensArray) = t.data
getbasis(t::TensBasis) = t.basis
getbasis(::TensCanonical{order,dim,T}) where {order,dim,T} = CanonicalBasis{dim,T}()
getvar(::TensOrthonormal{order}) where {order} = ntuple(_ -> :cov, Val(order))
getvar(::TensOrthonormal, i::Integer) = :cov
getvar(t::TensVar) = t.var
getvar(t::TensVar, i::Integer) = t.var[i]

#####################
# Display Functions #
#####################
for OP in (:show, :print, :display)
    @eval begin
        Base.$OP(U::FourthOrderTensor) = $OP(tomandel(U))

        # function Base.$OP(t::AbstractTens)
        #     $OP(typeof(t))
        #     print("→ array: ")
        #     $OP(getarray(t))
        #     print("→ basis: ")
        #     $OP(vecbasis(getbasis(t)))
        #     print("→ var: ")
        #     $OP(getvar(t))
        # end
        # function Base.$OP(t::TensOrthonormal)
        #     $OP(typeof(t))
        #     print("→ array: ")
        #     $OP(getarray(t))
        #     print("→ basis: ")
        #     $OP(vecbasis(getbasis(t)))
        # end

        Base.$OP(t::AbstractTens{order,dim,T}; vec = '𝐞', coords = ntuple(i -> i, dim)) where {order,dim,T} = printvec(t; vec= vec, coords = coords)
    end
end

function printvec(t::AbstractTens{order,dim,T}; vec = '𝐞', coords = ntuple(i -> i, dim)) where {order,dim,T}
    ind = CartesianIndices(t)
    ℬ = getbasis(t)
    firstprint = true
    s = ""
    for i ∈ ind
        x = t[i]
        if x ≠ zero(T)
            if !firstprint
                s *= " + "
            end
            if x ≠ one(T)
                s *= "(" * string(x) * ")"
            end
            j = Tuple(i)
            for k ∈ 1:order
                s *= strvecbasis(ℬ, coords[j[k]], invvar(getvar(t,k)); vec = vec)
                if k < order
                    s *= "⊗"
                end
            end
            firstprint = false
        end
    end
    if length(s)>0
        println(s)
    else
        println(0)
    end
end

########################
# Component extraction #
########################

"""
    components(t::AbstractTens{order,dim,T},ℬ::AbstractBasis{dim},var::NTuple{order,Symbol})
    components(t::AbstractTens{order,dim,T},ℬ::AbstractBasis{dim})
    components(t::AbstractTens{order,dim,T},var::NTuple{order,Symbol})

Extracts the components of a tensor for new variances and/or in a new basis

# Examples
```julia
julia> ℬ = Basis(Sym[0 1 1; 1 0 1; 1 1 0]) ;

julia> TV = Tens(Tensor{1,3}(i->symbols("v\$i",real=true)))
TensND.TensCanonical{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
 v₁
 v₂
 v₃
# basis: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# var: (:cont,)

julia> factor.(components(TV, ℬ, (:cont,)))
3-element Vector{Sym}:
 -(v1 - v2 - v3)/2
  (v1 - v2 + v3)/2
  (v1 + v2 - v3)/2

julia> components(TV, ℬ, (:cov,))
3-element Vector{Sym}:
 v₂ + v₃
 v₁ + v₃
 v₁ + v₂

julia> simplify.(components(TV, normalize(ℬ), (:cov,)))
3-element Vector{Sym}:
 sqrt(2)*(v2 + v3)/2
 sqrt(2)*(v1 + v3)/2
 sqrt(2)*(v1 + v2)/2

julia> TT = Tens(Tensor{2,3}((i,j)->symbols("t\$i\$j",real=true)))
TensND.TensCanonical{2, 3, Sym, Tensor{2, 3, Sym, 9}}
# data: 3×3 Tensor{2, 3, Sym, 9}:
 t₁₁  t₁₂  t₁₃
 t₂₁  t₂₂  t₂₃
 t₃₁  t₃₂  t₃₃
# basis: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# var: (:cont, :cont)

julia> components(TT, ℬ, (:cov,:cov))
3×3 Matrix{Sym}:
 t₂₂ + t₂₃ + t₃₂ + t₃₃  t₂₁ + t₂₃ + t₃₁ + t₃₃  t₂₁ + t₂₂ + t₃₁ + t₃₂
 t₁₂ + t₁₃ + t₃₂ + t₃₃  t₁₁ + t₁₃ + t₃₁ + t₃₃  t₁₁ + t₁₂ + t₃₁ + t₃₂
 t₁₂ + t₁₃ + t₂₂ + t₂₃  t₁₁ + t₁₃ + t₂₁ + t₂₃  t₁₁ + t₁₂ + t₂₁ + t₂₂

julia> factor.(components(TT, ℬ, (:cont,:cov)))
3×3 Matrix{Sym}:
 -(t12 + t13 - t22 - t23 - t32 - t33)/2  …  -(t11 + t12 - t21 - t22 - t31 - t32)/2
  (t12 + t13 - t22 - t23 + t32 + t33)/2      (t11 + t12 - t21 - t22 + t31 + t32)/2
  (t12 + t13 + t22 + t23 - t32 - t33)/2      (t11 + t12 + t21 + t22 - t31 - t32)/2
```
"""
components(t::AbstractTens) = getarray(t)

components(t::TensOrthonormal, ::NTuple) = getarray(t)

function components(
    t::TensOrthogonal{order,dim,T},
    var::NTuple{order,Symbol},
) where {order,dim,T}
    if var == getvar(t)
        return getarray(t)
    else
        m = Array(getarray(t))
        ℬ = getbasis(t)
        g_or_G = ntuple(i -> getvar(t, i) == var[i] ? I : metric(ℬ, var[i]), order)
        for ind ∈ CartesianIndices(m)
            m[ind] *= prod([g_or_G[i][ind[i], ind[i]] for i ∈ 1:order])
        end
        return m
    end
end

for B in (AbstractBasis, OrthogonalBasis, OrthonormalBasis)
    @eval function components(
        t::TensOrthogonal{order,dim,T},
        ℬ::$B{dim},
        var::NTuple{order,Symbol},
    ) where {order,dim,T}
        if ℬ == getbasis(t)
            return components(t, var)
        else
            return components(TensOrthonormal(t), ℬ, var)
        end
    end
end

function components(t::Tens{order,dim,T}, var::NTuple{order,Symbol}) where {order,dim,T}
    if var == getvar(t)
        return getarray(t)
    else
        m = Array(getarray(t))
        ec1 = ntuple(i -> i, Val(order))
        newcp = order + 1
        for i ∈ 1:order
            if getvar(t, i) ≠ var[i]
                g_or_G = metric(getbasis(t), var[i])
                ec2 = (i, newcp)
                ec3 = ntuple(j -> j ≠ i ? j : newcp, Val(order))
                m = einsum(EinCode((ec1, ec2), ec3), (m, g_or_G))
            end
        end
        return m
    end
end

function components(
    t::AbstractTens{order,dim,T},
    ℬ::AbstractBasis{dim},
    var::NTuple{order,Symbol},
) where {order,dim,T}
    if ℬ == getbasis(t)
        return components(t, var)
    else
        bb = Dict()
        for v1 ∈ (:cov, :cont), v2 ∈ (:cov, :cont)
            if v1 ∈ getvar(t) && v2 ∈ var
                bb[v1, v2] = vecbasis(getbasis(t), invvar(v1))' * vecbasis(ℬ, v2)
            end
        end
        m = Array(getarray(t))
        ec1 = ntuple(i -> i, Val(order))
        newcp = order + 1
        for i ∈ 1:order
            c = bb[getvar(t, i), var[i]]
            if c ≠ I
                ec2 = (i, newcp)
                ec3 = ntuple(j -> j ≠ i ? j : newcp, Val(order))
                m = einsum(EinCode((ec1, ec2), ec3), (m, c))
            end
        end
        return m
    end
end

function components(
    t::AbstractTens{order,dim,T},
    ℬ::OrthogonalBasis{dim},
    var::NTuple{order,Symbol},
) where {order,dim,T}
    if ℬ == getbasis(t)
        return components(t, var)
    else
        m = Array(components(t, relevant_OrthonormalBasis(ℬ)))
        Λ = Dict(:cont => inv.(ℬ.λ), :cov => ℬ.λ)
        for ind ∈ CartesianIndices(m)
            m[ind] *= prod([Λ[var[i]][ind[i]] for i ∈ 1:order])
        end
        return m
    end
end

components(t::AbstractTens{order,dim,T}, ℬ::AbstractBasis{dim}) where {order,dim,T} =
    components(t, ℬ, getvar(t))

components(t::AbstractTens{order,dim,T}, ℬ::OrthonormalBasis{dim}) where {order,dim,T} =
    components(t, ℬ, ntuple(_ -> :cont, Val(order)))

function components(
    t::TensOrthonormal{order,dim,T},
    ℬ::OrthonormalBasis{dim},
) where {order,dim,T}
    if ℬ == getbasis(t)
        return getarray(t)
    else
        bb = vecbasis(getbasis(t))' * vecbasis(ℬ)
        m = Array(getarray(t))
        ec1 = ntuple(i -> i, Val(order))
        newcp = order + 1
        for i ∈ 1:order
            if bb ≠ I
                ec2 = (i, newcp)
                ec3 = ntuple(j -> j ≠ i ? j : newcp, Val(order))
                m = einsum(EinCode((ec1, ec2), ec3), (m, bb))
            end
        end
        return m
    end
end

components(
    t::TensOrthonormal{order,dim,T},
    basis::OrthonormalBasis{dim},
    ::NTuple{order,Symbol},
) where {order,dim,T} = components(t, basis)

"""
    components_canon(t::AbstractTens)

Extracts the components of a tensor in the canonical basis
"""
components_canon(t::AbstractTens) =
    components(t, CanonicalBasis{getdim(t),eltype(t)}(), getvar(t))

components_canon(t::TensOrthonormal) = components(t, CanonicalBasis{getdim(t),eltype(t)}())

"""
    change_tens(t::AbstractTens{order,dim,T},ℬ::AbstractBasis{dim},var::NTuple{order,Symbol})
    change_tens(t::AbstractTens{order,dim,T},ℬ::AbstractBasis{dim})
    change_tens(t::AbstractTens{order,dim,T},var::NTuple{order,Symbol})

Rewrites the same tensor with components corresponding to new variances and/or to a new basis

```julia
julia> ℬ = Basis(Sym[0 1 1; 1 0 1; 1 1 0]) ;

julia> TV = Tens(Tensor{1,3}(i->symbols("v\$i",real=true)))
TensND.TensCanonical{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
 v₁
 v₂
 v₃
# basis: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# var: (:cont,)

julia> factor.(components(TV, ℬ, (:cont,)))
3-element Vector{Sym}:
 -(v1 - v2 - v3)/2
  (v1 - v2 + v3)/2
  (v1 + v2 - v3)/2

julia> ℬ₀ = Basis(Sym[0 1 1; 1 0 1; 1 1 1]) ;

julia> TV0 = change_tens(TV, ℬ₀)
Tens{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
     -v₁ + v₃
     -v₂ + v₃
 v₁ + v₂ - v₃
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 0  1  1
 1  0  1
 1  1  1
# var: (:cont,)
```
"""
function change_tens(t::AbstractTens, ℬ::AbstractBasis, newvar::NTuple)
    if ℬ == getbasis(t) && newvar == getvar(t)
        return t
    else
        return Tens(components(t, ℬ, newvar), ℬ, newvar)
    end
end

function change_tens(t::AbstractTens, newbasis::AbstractBasis)
    if newbasis == getbasis(t)
        return t
    else
        return Tens(components(t, newbasis, getvar(t)), newbasis, getvar(t))
    end
end

function change_tens(t::AbstractTens, newvar::NTuple)
    if newvar == getvar(t)
        return t
    else
        return Tens(components(t, getbasis(t), newvar), getbasis(t), newvar)
    end
end

"""
    change_tens_canon(t::AbstractTens{order,dim,T},var::NTuple{order,Symbol})

Rewrites the same tensor with components corresponding to the canonical basis

```julia
julia> ℬ = Basis(Sym[0 1 1; 1 0 1; 1 1 0]) ;

julia> TV = Tens(Tensor{1,3}(i->symbols("v\$i",real=true)), ℬ)
Tens{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
 v₁
 v₂
 v₃
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 0  1  1
 1  0  1
 1  1  1
# var: (:cont,)

julia> TV0 = change_tens_canon(TV)
TensND.TensCanonical{1, 3, Sym, Vec{3, Sym}}
# data: 3-element Vec{3, Sym}:
      v₂ + v₃
      v₁ + v₃
 v₁ + v₂ + v₃
# basis: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# var: (:cont,)
```
"""
change_tens_canon(t::AbstractTens) = change_tens(t, CanonicalBasis{getdim(t),eltype(t)}())

for OP in (:(simplify), :(factor), :(subs), :(diff))
    @eval SymPy.$OP(t::AbstractTens{order,dim,Sym}, args...; kwargs...) where {order,dim} =
        Tens($OP(getarray(t), args...; kwargs...), getbasis(t), getvar(t))
    @eval SymPy.$OP(t::AbstractArray{Sym}, args...; kwargs...) = $OP.(t, args...; kwargs...)
    @eval SymPy.$OP(t::Tensors.AllTensors{dim,Sym}, args...; kwargs...) where {dim} =
        Tensors.get_base(typeof(t))($OP.(Tensors.get_data(t), args...; kwargs...))
end

for OP in (:(trigsimp), :(expand_trig))
    @eval $OP(t::AbstractTens{order,dim}, args...; kwargs...) where {order,dim} =
        Tens($OP(getarray(t); kwargs...), getbasis(t), getvar(t))
    @eval $OP(t::AbstractArray{Sym}, args...; kwargs...) = sympy.$OP.(t, args...; kwargs...)
    @eval $OP(t::Tensors.AllTensors{dim,Sym}, args...; kwargs...) where {dim} =
        Tensors.get_base(typeof(t))(sympy.$OP.(Tensors.get_data(t), args...; kwargs...))
    @eval $OP(t::Sym, args...; kwargs...) = sympy.$OP.(t, args...; kwargs...)
end

# SymPy.simplify(t::AbstractTens{order,dim,Sym}; kwargs...) where {order,dim} = Tens(simplify(data(t); kwargs...), basis(t), getvar(t))
# SymPy.simplify(t::AbstractArray{Sym}; kwargs...) = simplify.(t; kwargs...)
# SymPy.simplify(t::Tensors.AllTensors{dim, Sym}; kwargs...) where {dim} = Tensors.get_base(typeof(t))(simplify.(Tensors.get_data(t); kwargs...))

# sympy.trigsimp(t::AbstractTens{order,dim,Sym}; kwargs...) where {order,dim} = Tens(trigsimp(data(t); kwargs...), basis(t), getvar(t))
# sympy.trigsimp(t::AbstractArray{Sym}; kwargs...) = trigsimp.(t; kwargs...)
# sympy.trigsimp(t::Tensors.AllTensors{dim, Sym}; kwargs...) where {dim} = Tensors.get_base(typeof(t))(trigsimp.(Tensors.get_data(t); kwargs...))

# SymPy.factor(t::AbstractTens{order,dim,Sym}; kwargs...) where {order,dim} = Tens(factor(data(t); kwargs...), basis(t), getvar(t))
# SymPy.factor(t::AbstractArray{Sym}; kwargs...) = factor.(t; kwargs...)
# SymPy.factor(t::Tensors.AllTensors{dim, Sym}; kwargs...) where {dim} = Tensors.get_base(typeof(t))(factor.(Tensors.get_data(t); kwargs...))

# SymPy.subs(t::AbstractTens{order,dim,Sym}, d...; kwargs...) where {order,dim} = Tens(subs(data(t), d...; kwargs...), basis(t), getvar(t))
# SymPy.subs(t::AbstractArray{Sym}, d...; kwargs...) = subs.(t, d...; kwargs...)
# SymPy.subs(t::Tensors.AllTensors{dim, Sym}, d...; kwargs...) where {dim} = Tensors.get_base(typeof(t))(subs.(Tensors.get_data(t), d...; kwargs...))

##############
# Operations #
##############

choose_best_basis(ℬ::AbstractBasis, ::AbstractBasis) = ℬ

choose_best_basis(ℬ::OrthonormalBasis, ::AbstractBasis) = ℬ
choose_best_basis(ℬ::OrthogonalBasis, ::AbstractBasis) = ℬ
choose_best_basis(::AbstractBasis, ℬ::OrthonormalBasis) = ℬ
choose_best_basis(::AbstractBasis, ℬ::OrthogonalBasis) = ℬ

choose_best_basis(::CanonicalBasis, ℬ::OrthonormalBasis) = ℬ
choose_best_basis(::CanonicalBasis, ℬ::OrthogonalBasis) = ℬ
choose_best_basis(ℬ::OrthonormalBasis, ::CanonicalBasis) = ℬ
choose_best_basis(ℬ::OrthogonalBasis, ::CanonicalBasis) = ℬ

choose_best_basis(::OrthonormalBasis, ℬ::OrthogonalBasis) = ℬ
choose_best_basis(ℬ::OrthogonalBasis, ::OrthonormalBasis) = ℬ
choose_best_basis(ℬ::CanonicalBasis, ::CanonicalBasis) = ℬ
choose_best_basis(ℬ::OrthonormalBasis, ::OrthonormalBasis) = ℬ
choose_best_basis(ℬ::OrthogonalBasis, ::OrthogonalBasis) = ℬ

function same_basis(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    ℬ = choose_best_basis(getbasis(t1), getbasis(t2))
    return change_tens(t1, ℬ), change_tens(t2, ℬ)
end

function same_basis_same_var(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    ℬ = choose_best_basis(getbasis(t1), getbasis(t2))
    return change_tens(t1, ℬ, getvar(t1)), change_tens(t2, ℬ, getvar(t1))
end

# same_basis(
#     t1::AbstractTens{order1,dim},
#     t2::AbstractTens{order2,dim},
# ) where {order1,order2,dim} = t1, change_tens(t2, basis(t1))

# same_basis_same_var(
#     t1::AbstractTens{order1,dim},
#     t2::AbstractTens{order2,dim},
# ) where {order1,order2,dim} = t1, change_tens(t2, basis(t1), getvar(t1))

for OP in (:(==), :(!=))
    @eval @inline function Base.$OP(
        t1::AbstractTens{order,dim},
        t2::AbstractTens{order,dim},
    ) where {order,dim}
        nt1, nt2 = same_basis_same_var(t1, t2)
        return $OP(getarray(nt1), getarray(nt2))
    end
end

for OP in (:+, :-)
    @eval @inline function Base.$OP(
        t1::AbstractTens{order,dim},
        t2::AbstractTens{order,dim},
    ) where {order,dim}
        nt1, nt2 = same_basis_same_var(t1, t2)
        return Tens($OP(getarray(nt1), getarray(nt2)), getbasis(nt1), getvar(nt1))
    end
    @eval @inline function Base.$OP(
        t1::AllTensOrthogonal{order,dim,T},
        t2::UniformScaling{T},
    ) where {order,dim,T<:SymPy.SymbolicObject}
        nt1 = TensOrthonormal(t1)
        return Tens($OP(getarray(nt1), t2), getbasis(nt1), getvar(nt1))
    end
    @eval @inline function Base.$OP(
        t1::UniformScaling{T},
        t2::AllTensOrthogonal{order,dim,T},
    ) where {order,dim,T<:SymPy.SymbolicObject}
        nt2 = TensOrthonormal(t2)
        return Tens($OP(t1, getarray(nt2)), getbasis(nt2), getvar(nt2))
    end
end

@inline Base.:-(t::AbstractTens) = Tens(.-(getarray(t)), getbasis(t), getvar(t))
@inline Base.:*(α::Number, t::AbstractTens) = Tens(α * getarray(t), getbasis(t), getvar(t))
@inline Base.:*(t::AbstractTens, α::Number) = Tens(α * getarray(t), getbasis(t), getvar(t))
@inline Base.:/(t::AbstractTens, α::Number) = Tens(getarray(t) / α, getbasis(t), getvar(t))

@inline Base.inv(t::AbstractTens{2}) =
    Tens(inv(getarray(t)), getbasis(t), (invvar(getvar(t, 2)), invvar(getvar(t, 1))))
@inline Base.inv(t::AbstractTens{4}) = Tens(
    inv(getarray(t)),
    getbasis(t),
    (
        invvar(getvar(t, 3)),
        invvar(getvar(t, 4)),
        invvar(getvar(t, 1)),
        invvar(getvar(t, 2)),
    ),
)

"""
    KM(t::AbstractTens{order,dim}; kwargs...)
    KM(t::AbstractTens{order,dim}, var::NTuple{order,Symbol}, b::AbstractBasis{dim}; kwargs...)

Writes the components of a second or fourth order tensor in Kelvin-Mandel notation

# Examples
```julia
julia> σ = Tens(SymmetricTensor{2,3}((i, j) -> symbols("σ\$i\$j", real = true))) ;

julia> KM(σ)
6-element Vector{Sym}:
         σ11
         σ22
         σ33
      √2⋅σ32
      √2⋅σ31
      √2⋅σ21

julia> C = Tens(SymmetricTensor{4,3}((i, j, k, l) -> symbols("C\$i\$j\$k\$l", real = true))) ;

julia> KM(C)
6×6 Matrix{Sym}:
         C₁₁₁₁     C₁₁₂₂     C₁₁₃₃  √2⋅C₁₁₃₂  √2⋅C₁₁₃₁  √2⋅C₁₁₂₁
         C₂₂₁₁     C₂₂₂₂     C₂₂₃₃  √2⋅C₂₂₃₂  √2⋅C₂₂₃₁  √2⋅C₂₂₂₁
         C₃₃₁₁     C₃₃₂₂     C₃₃₃₃  √2⋅C₃₃₃₂  √2⋅C₃₃₃₁  √2⋅C₃₃₂₁
      √2⋅C₃₂₁₁  √2⋅C₃₂₂₂  √2⋅C₃₂₃₃   2⋅C₃₂₃₂   2⋅C₃₂₃₁   2⋅C₃₂₂₁
      √2⋅C₃₁₁₁  √2⋅C₃₁₂₂  √2⋅C₃₁₃₃   2⋅C₃₁₃₂   2⋅C₃₁₃₁   2⋅C₃₁₂₁
      √2⋅C₂₁₁₁  √2⋅C₂₁₂₂  √2⋅C₂₁₃₃   2⋅C₂₁₃₂   2⋅C₂₁₃₁   2⋅C₂₁₂₁
```
"""
KM(t::Tensors.AllTensors; kwargs...) = tomandel(t; kwargs...)
KM(t::AbstractTens; kwargs...) = tomandel(getarray(t); kwargs...)

KM(
    t::AbstractTens{order,dim},
    b::AbstractBasis{dim},
    var::NTuple{order,Symbol},
    kwargs...,
) where {order,dim} = tomandel(tensor_or_array(components(t, b, var)); kwargs...)

KM(t::AbstractTens{order,dim}, b::AbstractBasis{dim}; kwargs...) where {order,dim} =
    tomandel(tensor_or_array(components(t, b)); kwargs...)

KM(t::AbstractArray; kwargs...) = KM(Tens(t); kwargs...)
KM(t::AbstractArray, b::AbstractBasis; kwargs...) = KM(Tens(t), b; kwargs...)

const select_type_KM = Dict(
    (6, 6) => SymmetricTensor{4,3},
    (9, 9) => Tensor{4,3},
    (3, 3) => SymmetricTensor{4,2},
    (4, 4) => Tensor{4,2},
    (6,) => SymmetricTensor{2,3},
    (9,) => Tensor{2,3},
    (3,) => SymmetricTensor{2,2},
    (4,) => Tensor{2,2},
)

"""
    invKM(v::AbstractVecOrMat; kwargs...)

Defines a tensor from a Kelvin-Mandel vector or matrix representation
"""
invKM(TT::Type{<:Tensors.AllTensors}, v::AbstractVecOrMat; kwargs...) =
    Tens(frommandel(TT, v; kwargs...))
invKM(v::AbstractVecOrMat; kwargs...) = invKM(select_type_KM[size(v)], v; kwargs...)

"""
    otimes(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a tensor product between two tensors

`(aⁱeᵢ) ⊗ (bʲeⱼ) = aⁱbʲ eᵢ⊗eⱼ`
"""
function Tensors.otimes(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    if t1 == t2
        return otimes(t1)
    else
        nt1, nt2 = same_basis(t1, t2)
        data = otimes(getarray(nt1), getarray(nt2))
        var = (getvar(nt1)..., getvar(nt2)...)
        return Tens(data, getbasis(nt1), var)
    end
end

function Tensors.otimes(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    if t1 == t2
        return otimes(t1)
    else
        nt1, nt2 = same_basis(t1, t2)
        data = otimes(getarray(nt1), getarray(nt2))
        return Tens(data, getbasis(nt1))
    end
end

Tensors.otimes(v::AbstractTens{1,dim}) where {dim} =
    Tens(otimes(getarray(v)), getbasis(v), (getvar(v)..., getvar(v)...))

Tensors.otimes(v::TensOrthonormal{1,dim}) where {dim} =
    Tens(otimes(getarray(v)), getbasis(v))

@inline function Tensors.otimes(S::SymmetricTensor{2,dim}) where {dim}
    return SymmetricTensor{4,dim}(@inline function (i, j, k, l)
        @inbounds S[i, j] * S[k, l]
    end)
end

Tensors.otimes(t::AbstractTens{2,dim}) where {dim} =
    Tens(otimes(getarray(t)), getbasis(t), (getvar(t)..., getvar(t)...))

Tensors.otimes(t::TensOrthonormal{2,dim}) where {dim} =
    Tens(otimes(getarray(t)), getbasis(t))

function scontract(
    t1::AbstractArray{T1,order1},
    t2::AbstractArray{T2,order2},
) where {T1,T2,order1,order2}
    ec1 = ntuple(i -> i, order1)
    ec2 = ntuple(i -> order1 - 1 + i, order2)
    ec3 = (ec1[begin:end-1]..., ec2[begin+1:end]...)
    return einsum(EinCode((ec1, ec2), ec3), (AbstractArray{T1}(t1), AbstractArray{T2}(t2)))
end

scontract(t1::AbstractArray{T1,1}, t2::AbstractArray{T2,1}) where {T1,T2} =
    dot(AbstractArray{T1}(t1), AbstractArray{T2}(t2))

for TT1 ∈ (Vec, SecondOrderTensor), TT2 ∈ (Vec, SecondOrderTensor)
    @eval scontract(S1::$TT1, S2::$TT2) = dot(S1, S2)
end

"""
    dot(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a contracted product between two tensors

`a ⋅ b = aⁱbⱼ`
"""
function LinearAlgebra.dot(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    var = (invvar(getvar(nt1)[end]), getvar(nt2)[begin+1:end]...)
    nt2 = change_tens(nt2, getbasis(nt2), var)
    data = scontract(getarray(nt1), getarray(nt2))
    var = (getvar(nt1)[begin:end-1]..., getvar(nt2)[begin+1:end]...)
    return Tens(data, getbasis(nt1), var)
end

function LinearAlgebra.dot(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = scontract(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

function LinearAlgebra.dot(t1::AbstractTens{1,dim}, t2::AbstractTens{1,dim}) where {dim}
    nt1, nt2 = same_basis(t1, t2)
    var = (invvar(getvar(nt1)[end]), getvar(nt2)[begin+1:end]...)
    nt2 = change_tens(nt2, getbasis(nt2), var)
    return scontract(getarray(nt1), getarray(nt2))
end

function LinearAlgebra.dot(
    t1::TensOrthonormal{1,dim},
    t2::TensOrthonormal{1,dim},
) where {dim}
    nt1, nt2 = same_basis(t1, t2)
    return scontract(getarray(nt1), getarray(nt2))
end

LinearAlgebra.dot(v::TensOrthonormal{1,dim}) where {dim} = dot(getarray(v))

LinearAlgebra.dot(t::TensOrthonormal{2,dim}) where {dim} = dot(getarray(t))

LinearAlgebra.norm(u::AbstractTens{1,dim}) where {dim} = √(dot(u, u))

LinearAlgebra.norm(t::AbstractTens{2,dim}) where {dim} = √(dot(t, t))

"""
    contract(t::AbstractTens{order,dim}, i::Integer, j::Integer)

Calculates the tensor obtained after contraction with respect to the indices `i` and `j`
"""
function contract(t::AbstractTens{order,dim}, i::Integer, j::Integer) where {order,dim}
    var = ntuple(k -> k == j ? invvar(getvar(t, i)) : getvar(t, k), order)
    nt = change_tens(t, getbasis(t), var)
    data = contract(getarray(nt), i, j)
    m = min(i, j)
    M = max(i, j)
    var = (getvar(nt)[1:m-1]..., getvar(nt)[m+1:M-1]..., getvar(nt)[M+1:order]...)
    return Tens(data, getbasis(nt), var)
end

function contract(t::AbstractTens{2,dim}, i::Integer, j::Integer) where {dim}
    var = ntuple(k -> k == j ? invvar(getvar(t, i)) : getvar(t, k), 2)
    nt = change_tens(t, getbasis(t), var)
    return contract(getarray(nt), i, j)
end

contract(t::TensOrthonormal{order,dim}, i::Integer, j::Integer) where {order,dim} =
    Tens(contract(getarray(t), i, j), getbasis(t))

contract(t::TensOrthonormal{2,dim}, i::Integer, j::Integer) where {dim} =
    contract(getarray(t), i, j)

LinearAlgebra.tr(t::AbstractTens{2}) = contract(t, 1, 2)

"""
    dcontract(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a double contracted product between two tensors

`𝛔 ⊡ 𝛆 = σⁱʲεᵢⱼ`
`𝛔 = ℂ ⊡ 𝛆`

# Examples
```julia
julia> 𝛆 = Tens(SymmetricTensor{2,3}((i, j) -> symbols("ε\$i\$j", real = true))) ;

julia> k, μ = symbols("k μ", real =true) ;

julia> ℂ = 3k * t𝕁() + 2μ * t𝕂() ;

julia> 𝛔 = ℂ ⊡ 𝛆
Tens{2, 3, Sym, Sym, SymmetricTensor{2, 3, Sym, 6}, CanonicalBasis{3, Sym}}
# data: 3×3 SymmetricTensor{2, 3, Sym, 6}:
 ε11*(k + 4*μ/3) + ε22*(k - 2*μ/3) + ε33*(k - 2*μ/3)                                              2⋅ε21⋅μ                                              2⋅ε31⋅μ
                                             2⋅ε21⋅μ  ε11*(k - 2*μ/3) + ε22*(k + 4*μ/3) + ε33*(k - 2*μ/3)                                              2⋅ε32⋅μ
                                             2⋅ε31⋅μ                                              2⋅ε32⋅μ  ε11*(k - 2*μ/3) + ε22*(k - 2*μ/3) + ε33*(k + 4*μ/3)
# var: (:cont, :cont)
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 1  0  0
 0  1  0
 0  0  1
```
"""
function Tensors.dcontract(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    var =
        (invvar(getvar(nt1)[end-1]), invvar(getvar(nt1)[end]), getvar(nt2)[begin+2:end]...)
    nt2 = change_tens(nt2, getbasis(nt2), var)
    data = Tensors.dcontract(getarray(nt1), getarray(nt2))
    var = (getvar(nt1)[begin:end-2]..., getvar(nt2)[begin+2:end]...)
    return Tens(data, getbasis(nt1), var)
end

function Tensors.dcontract(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = Tensors.dcontract(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

function Tensors.dcontract(t1::AbstractTens{2,dim}, t2::AbstractTens{2,dim}) where {dim}
    nt1, nt2 = same_basis(t1, t2)
    var =
        (invvar(getvar(nt1)[end-1]), invvar(getvar(nt1)[end]), getvar(nt2)[begin+2:end]...)
    nt2 = change_tens(nt2, getbasis(nt2), var)
    return Tensors.dcontract(getarray(nt1), getarray(nt2))
end

function Tensors.dcontract(
    t1::TensOrthonormal{2,dim},
    t2::TensOrthonormal{2,dim},
) where {dim}
    nt1, nt2 = same_basis(t1, t2)
    return Tensors.dcontract(getarray(nt1), getarray(nt2))
end

"""
    dotdot(v1::AbstractTens{order1,dim}, S::AbstractTens{orderS,dim}, v2::AbstractTens{order2,dim})

Defines a bilinear operator `𝐯₁⋅𝕊⋅𝐯₂`

# Examples
```julia
julia> n = Tens(Sym[0, 0, 1]) ;

julia> k, μ = symbols("k μ", real =true) ;

julia> ℂ = 3k * t𝕁() + 2μ * t𝕂() ;

julia> dotdot(n,ℂ,n) # Acoustic tensor
3×3 Tens{2, 3, Sym, Sym, Tensor{2, 3, Sym, 9}, CanonicalBasis{3, Sym}}:
 μ  0          0
 0  μ          0
 0  0  k + 4*μ/3
```
"""
function Tensors.dotdot(
    v1::AbstractTens{order1,dim},
    S::AbstractTens{orderS,dim},
    v2::AbstractTens{order2,dim},
) where {order1,orderS,order2,dim}
    nS, nv1 = same_basis(S, v1)
    nS, nv2 = same_basis(S, v2)
    var = (invvar(getvar(nS)[begin]),)
    nv1 = change_tens(nv1, getbasis(nv1), var)
    var = (invvar(getvar(nS)[end]),)
    nv2 = change_tens(nv2, getbasis(nv2), var)
    data = dotdot(nv1.data, nS.data, nv2.data)
    var = (getvar(nS)[begin+1], getvar(nS)[end-1])
    return Tens(data, getbasis(nS), var)
end

function Tensors.dotdot(
    v1::TensOrthonormal{order1,dim},
    S::TensOrthonormal{orderS,dim},
    v2::TensOrthonormal{order2,dim},
) where {order1,orderS,order2,dim}
    nS, nv1 = same_basis(S, v1)
    nS, nv2 = same_basis(S, v2)
    data = dotdot(nv1.data, nS.data, nv2.data)
    return Tens(data, getbasis(nS))
end

"""
    qcontract(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a quadruple contracted product between two tensors

`𝔸 ⊙ 𝔹 = AᵢⱼₖₗBⁱʲᵏˡ`

# Examples
```julia
julia> 𝕀 = t𝕀(Sym) ; 𝕁 = t𝕁(Sym) ; 𝕂 = t𝕂(Sym) ;

julia> 𝕀 ⊙ 𝕀
6

julia> 𝕁 ⊙ 𝕀
1

julia> 𝕂 ⊙ 𝕀
5

julia> 𝕂 ⊙ 𝕁
0
```
"""
function qcontract(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    var = (
        invvar(getvar(nt1)[end-3]),
        invvar(getvar(nt1)[end-2]),
        invvar(getvar(nt1)[end-1]),
        invvar(getvar(nt1)[end]),
        getvar(nt2)[begin+4:end]...,
    )
    nt2 = change_tens(nt2, getbasis(nt2), var)
    data = qcontract(getarray(nt1), getarray(nt2))
    var = (getvar(nt1)[begin:end-4]..., getvar(nt2)[begin+4:end]...)
    return Tens(data, getbasis(nt1), var)
end

function qcontract(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = qcontract(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

function qcontract(t1::AbstractTens{4,dim}, t2::AbstractTens{4,dim}) where {dim}
    nt1, nt2 = same_basis(t1, t2)
    var = (
        invvar(getvar(nt1)[end-3]),
        invvar(getvar(nt1)[end-2]),
        invvar(getvar(nt1)[end-1]),
        invvar(getvar(nt1)[end]),
        getvar(nt2)[begin+4:end]...,
    )
    nt2 = change_tens(nt2, getbasis(nt2), var)
    return qcontract(getarray(nt1), getarray(nt2))
end

function qcontract(t1::TensOrthonormal{4,dim}, t2::TensOrthonormal{4,dim}) where {dim}
    nt1, nt2 = same_basis(t1, t2)
    return qcontract(getarray(nt1), getarray(nt2))
end

"""
    otimesu(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a special tensor product between two tensors of at least second order

`(𝐚 ⊠ 𝐛) ⊡ 𝐩 = 𝐚⋅𝐩⋅𝐛 = aⁱᵏbʲˡpₖₗ eᵢ⊗eⱼ`
"""
function Tensors.otimesu(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = otimesu(getarray(nt1), getarray(nt2))
    var = (
        getvar(nt1)[begin:end-1]...,
        getvar(nt2)[begin],
        getvar(nt1)[end],
        getvar(nt2)[begin+1:end]...,
    )
    return Tens(data, getbasis(nt1), var)
end

function Tensors.otimesu(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = otimesu(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

function Tensors.otimesl(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = otimesl(getarray(nt1), getarray(nt2))
    var = (
        getvar(nt1)[begin:end-1]...,
        getvar(nt2)[begin+1],
        getvar(nt1)[end],
        getvar(nt2)[begin],
        getvar(nt2)[begin+2:end]...,
    )
    return Tens(data, getbasis(nt1), var)
end

function Tensors.otimesl(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = otimesl(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

otimesul(S1::SecondOrderTensor{dim}, S2::SecondOrderTensor{dim}) where {dim} =
    symmetric(otimesu(S1, S2))

"""
    otimesul(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a special tensor product between two tensors of at least second order

`(𝐚 ⊠ˢ 𝐛) ⊡ 𝐩 = (𝐚 ⊠ 𝐛) ⊡ (𝐩 + ᵗ𝐩)/2  = 1/2(aⁱᵏbʲˡ+aⁱˡbʲᵏ) pₖₗ eᵢ⊗eⱼ`
"""
function otimesul(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    var = (getvar(nt1)[end-1], getvar(nt1)[end], getvar(nt2)[begin+2:end]...)
    nt2 = change_tens(nt2, getbasis(nt2), var)
    data = otimesul(getarray(nt1), getarray(nt2))
    var = (getvar(nt1)..., getvar(nt2)...)
    return Tens(data, getbasis(nt1), var)
end

function otimesul(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = otimesul(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

"""
    sotimes(t1::AbstractTens{order1,dim}, t2::AbstractTens{order2,dim})

Defines a symmetric tensor product between two tensors

`(aⁱeᵢ) ⊗ˢ (bʲeⱼ) = 1/2(aⁱbʲ + aʲbⁱ) eᵢ⊗eⱼ`
"""
function sotimes(
    t1::AbstractTens{order1,dim},
    t2::AbstractTens{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    var = (getvar(nt1)[end], getvar(nt2)[begin+1:end]...)
    nt2 = change_tens(nt2, getbasis(nt2), var)
    data = sotimes(getarray(nt1), getarray(nt2))
    var = (getvar(nt1)..., getvar(nt2)...)
    return Tens(data, getbasis(nt1), var)
end

function sotimes(
    t1::AbstractTens{1,dim},
    t2::AbstractTens{1,dim},
) where {dim}
    if t1 == t2
        return otimes(t1)
    else
        nt1, nt2 = same_basis(t1, t2)
        var = (getvar(nt1)[end], getvar(nt2)[begin+1:end]...)
        nt2 = change_tens(nt2, getbasis(nt2), var)
        data = sotimes(getarray(nt1), getarray(nt2))
        var = (getvar(nt1)..., getvar(nt2)...)
        return Tens(data, getbasis(nt1), var)
    end
end

function sotimes(
    t1::TensOrthonormal{order1,dim},
    t2::TensOrthonormal{order2,dim},
) where {order1,order2,dim}
    nt1, nt2 = same_basis(t1, t2)
    data = sotimes(getarray(nt1), getarray(nt2))
    return Tens(data, getbasis(nt1))
end

function sotimes(
    t1::TensOrthonormal{1,dim},
    t2::TensOrthonormal{1,dim},
) where {dim}
    if t1 == t2
        return otimes(t1)
    else
        nt1, nt2 = same_basis(t1, t2)
        data = sotimes(getarray(nt1), getarray(nt2))
        return Tens(data, getbasis(nt1))
    end
end

Base.transpose(t::TensArray{order,dim,T,<:SecondOrderTensor}) where {order,dim,T} =
    Tens(transpose(getarray(t)), getbasis(t), (getvar(t)[2], getvar(t)[1]))

Base.transpose(t::TensOrthonormal{order,dim,T,<:SecondOrderTensor}) where {order,dim,T} =
    Tens(transpose(getarray(t)), getbasis(t))

Base.transpose(t::TensArray{order,dim,T,<:FourthOrderTensor}) where {order,dim,T} = Tens(
    Tensors.transpose(getarray(t)),
    getbasis(t),
    (getvar(t)[2], getvar(t)[1], getvar(t)[4], getvar(t)[3]),
)

Base.transpose(t::TensOrthonormal{order,dim,T,<:FourthOrderTensor}) where {order,dim,T} =
    Tens(Tensors.transpose(getarray(t)), getbasis(t))

Tensors.majortranspose(t::TensArray{order,dim,T,<:FourthOrderTensor}) where {order,dim,T} =
    Tens(
        majortranspose(getarray(t)),
        getbasis(t),
        (getvar(t)[3], getvar(t)[4], getvar(t)[1], getvar(t)[2]),
    )

Tensors.majortranspose(
    t::TensOrthonormal{order,dim,T,<:FourthOrderTensor},
) where {order,dim,T} = Tens(majortranspose(getarray(t)), getbasis(t))

Tensors.minortranspose(t::TensArray{order,dim,T,<:FourthOrderTensor}) where {order,dim,T} =
    Tens(
        minortranspose(getarray(t)),
        getbasis(t),
        (getvar(t)[2], getvar(t)[1], getvar(t)[4], getvar(t)[3]),
    )

Tensors.minortranspose(
    t::TensOrthonormal{order,dim,T,<:FourthOrderTensor},
) where {order,dim,T} = Tens(minortranspose(getarray(t)), getbasis(t))

function tensbasis(ℬ::AbstractBasis, i::Integer, ::Val{:cov})
    T = eltype(ℬ)
    t = [T(Int(j == i)) for j ∈ 1:getdim(ℬ)]
    return Tens(t, ℬ, (:cont,))
end
function tensbasis(ℬ::AbstractBasis, i::Integer, ::Val{:cont})
    T = eltype(ℬ)
    t = [T(Int(j == i)) for j ∈ 1:getdim(ℬ)]
    return Tens(t, ℬ, (:cov,))
end
tensbasis(ℬ::AbstractBasis, i::Integer, var = :cov) = tensbasis(ℬ, i, Val(var))
tensbasis(ℬ::AbstractBasis, var = :cov) = ntuple(i -> tensbasis(ℬ, i, Val(var)), getdim(ℬ))

export AbstractTens, Tens
export proj_tens, best_sym_tens
export getorder, arraytype, getdata, getarray, getbasis, getvar
export printvec
export components, components_canon, change_tens, change_tens_canon
export trigsimp, expand_trig
export KM, invKM
export getbasis, getvar
export tensbasis


text:
"""
    one_hot_encode(df::DataFrame;
                   cols::Vector{Symbol}=Vector{Symbol}(),
                   drop_original::Bool=false)

Utility function to perform one-hot-encoding in DataFrame. This will add new columns
with names `<original_col_name>_<value>`.

Following options can be passed to modify behavior.

* `cols` - Vector of Symbol to specify which columns to be encoded. Defaults to
  empty, which means all features will be encoded.
* `drop_original` - If true, this will drop the original feature set from resulting
  DataFrame. This defaults to false.

# Example

```jldoctest
julia> using RDatasets, FeatureSelectors

julia> titanic = dataset("datasets", "Titanic");

julia> first(one_hot_encode(titanic[:, [:Class, :Sex, :Age]]), 3)
3×11 DataFrame
 Row │ Class    Sex      Age      Class_1st  Class_2nd  Class_3rd  Class_Crew  ⋯
     │ String7  String7  String7  Bool       Bool       Bool       Bool        ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ 1st      Male     Child         true      false      false       false  ⋯
   2 │ 2nd      Male     Child        false       true      false       false
   3 │ 3rd      Male     Child        false      false       true       false
                                                               4 columns omitted

julia> first(one_hot_encode(titanic[:, [:Class, :Sex, :Age]], cols=[:Class], drop_original=true), 3)
3×6 DataFrame
 Row │ Sex      Age      Class_1st  Class_2nd  Class_3rd  Class_Crew
     │ String7  String7  Bool       Bool       Bool       Bool
─────┼───────────────────────────────────────────────────────────────
   1 │ Male     Child         true      false      false       false
   2 │ Male     Child        false       true      false       false
   3 │ Male     Child        false      false       true       false

```
"""
function one_hot_encode(
    df::DataFrame;
    cols::Vector{Symbol} = Vector{Symbol}(),
    drop_original::Bool = false,
)
    result_df = df
    for col_name in names(df)
        # Ignore if column not in cols
        if !(Symbol(col_name) in cols) && !isempty(cols)
            continue
        end
        tmp_df = DataFrame(Array(transpose(indicatormat(df[:, col_name]))), :auto)
        rename!(tmp_df, [Symbol("$(col_name)_$x") for x in sort(unique(df[:, col_name]))])
        result_df = drop_original ? hcat(result_df, tmp_df)[:, Not(col_name)] :
            hcat(result_df, tmp_df)
    end
    result_df
end

"""
    calculate_feature_importance(method::Function, X::DataFrame, y::Vector)

Calculate feature importance defined by `method`. Similar with `select_features`, this can
take `X` in `DataFrame` or splitted into `X_data` in `Matrix` and `X_features` in `Vector`.

This function will return `Dict` with feature names as key and scores as value.
"""
function calculate_feature_importance(
    method::Function,
    X_data::Matrix,
    X_features::Vector,
    y::Vector
)
    scores = method(X_data, y)
    Dict(feature_name=>score for (feature_name,score) in zip(X_features, scores))
end

calculate_feature_importance(
    method::Function,
    X::DataFrame,
    y::Vector
) = calculate_feature_importance(
    method,
    Matrix(X),
    names(X),
    y,
)

ext: jl


text:
######################
### Mutable states ###
######################

mutable struct ThreePhaseState
    i           :: Int
    window_size :: Int
    next_window :: Int
end

################
### Adaptors ###
################

# TODO: currently only StanNUTSAdaptor has the filed `n_adapts`. maybe we could unify all
# Acknowledgement: this adaption settings is mimicing Stan's 3-phase adaptation.
struct StanNUTSAdaptor <: AbstractCompositeAdaptor
    n_adapts    :: Int
    pc          :: AbstractPreConditioner
    ssa         :: StepSizeAdaptor
    init_buffer :: Int
    term_buffer :: Int
    state       :: ThreePhaseState
end

function StanNUTSAdaptor(n_adapts::Int, pc::AbstractPreConditioner, ssa::StepSizeAdaptor,
                         init_buffer::Int=75, term_buffer::Int=50, window_size::Int=25)
    next_window = init_buffer + window_size - 1
    return StanNUTSAdaptor(n_adapts, pc, ssa, init_buffer, term_buffer, ThreePhaseState(0, window_size, next_window))
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/windowed_adaptation.hpp
function is_in_window(tp::StanNUTSAdaptor)
    return (tp.state.i >= tp.init_buffer) &&
           (tp.state.i < tp.n_adapts - tp.term_buffer) &&
           (tp.state.i != tp.n_adapts)
end

function is_window_end(tp::StanNUTSAdaptor)
    return (tp.state.i == tp.state.next_window) &&
           (tp.state.i != tp.n_adapts)
end

is_final(tp::StanNUTSAdaptor) = tp.state.i == tp.n_adapts

function compute_next_window!(tp::StanNUTSAdaptor)
    if ~(tp.state.next_window == tp.n_adapts - tp.term_buffer - 1)
        tp.state.window_size *= 2
        tp.state.next_window = tp.state.i + tp.state.window_size
        if ~(tp.state.next_window == tp.n_adapts - tp.term_buffer - 1)
            next_window_boundary = tp.state.next_window + 2 * tp.state.window_size
            if (next_window_boundary >= tp.n_adapts - tp.term_buffer)
                tp.state.next_window = tp.n_adapts - tp.term_buffer - 1
            end
        end
    end
end

function adapt!(tp::StanNUTSAdaptor, θ::AbstractVector{<:Real}, α::AbstractFloat)
    tp.state.i += 1

    adapt!(tp.ssa, θ, α)

    # Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
    if is_in_window(tp)
        # We accumlate stats from θ online and only trigger the update of M⁻¹ in the end of window.
        is_update_M⁻¹ = is_window_end(tp)
        adapt!(tp.pc, θ, α, is_update_M⁻¹)
    end

    if is_window_end(tp)
        reset!(tp.ssa)
        reset!(tp.pc)
        # If window ends, compute next window
        compute_next_window!(tp)
    end

    if is_final(tp)
        finalize!(tp.ssa)
    end
end


text:
# Endowments

function init_test_endowments(lastWorkAge :: Int)
    n = 3;
    effGridV = collect(LinRange(1.0, 3.0, n));
    effProbV = collect(LinRange(1.0, 0.5, n));
    effProbV ./= sum(effProbV);

    trMatrix = zeros(n, n);
    for j = 1 : n
        prV = collect(LinRange(0.5 + j, 0.5, n));
        trMatrix[:, j] .= prV ./ sum(prV);
    end

    endow = Endowments(effGridV, effProbV, trMatrix);
    return endow
end

function validate_endowments(e :: Endowments)
    isValid = validate_efficiency_grid(e.effGridV)  &&
        validate_transition_matrix(e.trMatrix);
    return isValid
end

function validate_efficiency_grid(gridV :: Vector{Float64})
    isValid = all(gridV .> 0.0)  &&  all(diff(gridV) .> 0.0);
    return isValid
end

function validate_transition_matrix(trMatrix :: Matrix{Float64})
    isValid = all(trMatrix .>= 0.0)  &&  all(trMatrix .<= 1.0);
    n, m = size(trMatrix);
    isValid = isValid  &&  (n == m);
    for j = 1 : n
        isValid = isValid  &&  isapprox(sum(trMatrix[:, j]), 1.0);
    end
    return isValid
end

## -----------  Access

n_efficiencies(e :: Endowments) = length(e.effGridV);

# Probability of e' | e
prob_eprime(e :: Endowments, eIdx) = e.trMatrix[:, eIdx];

efficiency(e :: Endowments, eIdx) = e.effGridV[eIdx];

# -----------------


text:
# Automatically generated using Clang.jl

const LLVMBool = Cint
struct LLVMOpaqueMemoryBuffer end
const LLVMMemoryBufferRef = Ptr{LLVMOpaqueMemoryBuffer}
struct LLVMOpaqueContext end
const LLVMContextRef = Ptr{LLVMOpaqueContext}
struct LLVMOpaqueModule end
const LLVMModuleRef = Ptr{LLVMOpaqueModule}
struct LLVMOpaqueType end
const LLVMTypeRef = Ptr{LLVMOpaqueType}
struct LLVMOpaqueValue end
const LLVMValueRef = Ptr{LLVMOpaqueValue}
struct LLVMOpaqueBasicBlock end
const LLVMBasicBlockRef = Ptr{LLVMOpaqueBasicBlock}
struct LLVMOpaqueMetadata end
const LLVMMetadataRef = Ptr{LLVMOpaqueMetadata}
struct LLVMOpaqueBuilder end
const LLVMBuilderRef = Ptr{LLVMOpaqueBuilder}
struct LLVMOpaqueDIBuilder end
const LLVMDIBuilderRef = Ptr{LLVMOpaqueDIBuilder}
struct LLVMOpaqueModuleProvider end
const LLVMModuleProviderRef = Ptr{LLVMOpaqueModuleProvider}
struct LLVMOpaquePassManager end
const LLVMPassManagerRef = Ptr{LLVMOpaquePassManager}
struct LLVMOpaquePassRegistry end
const LLVMPassRegistryRef = Ptr{LLVMOpaquePassRegistry}
struct LLVMOpaqueUse end
const LLVMUseRef = Ptr{LLVMOpaqueUse}
struct LLVMOpaqueAttributeRef end
const LLVMAttributeRef = Ptr{LLVMOpaqueAttributeRef}
struct LLVMOpaqueDiagnosticInfo end
const LLVMDiagnosticInfoRef = Ptr{LLVMOpaqueDiagnosticInfo}

@cenum LLVMVerifierFailureAction::UInt32 begin
    LLVMAbortProcessAction = 0
    LLVMPrintMessageAction = 1
    LLVMReturnStatusAction = 2
end

const LLVMFatalErrorHandler = Ptr{Cvoid}

@cenum LLVMOpcode::UInt32 begin
    LLVMRet = 1
    LLVMBr = 2
    LLVMSwitch = 3
    LLVMIndirectBr = 4
    LLVMInvoke = 5
    LLVMUnreachable = 7
    LLVMAdd = 8
    LLVMFAdd = 9
    LLVMSub = 10
    LLVMFSub = 11
    LLVMMul = 12
    LLVMFMul = 13
    LLVMUDiv = 14
    LLVMSDiv = 15
    LLVMFDiv = 16
    LLVMURem = 17
    LLVMSRem = 18
    LLVMFRem = 19
    LLVMShl = 20
    LLVMLShr = 21
    LLVMAShr = 22
    LLVMAnd = 23
    LLVMOr = 24
    LLVMXor = 25
    LLVMAlloca = 26
    LLVMLoad = 27
    LLVMStore = 28
    LLVMGetElementPtr = 29
    LLVMTrunc = 30
    LLVMZExt = 31
    LLVMSExt = 32
    LLVMFPToUI = 33
    LLVMFPToSI = 34
    LLVMUIToFP = 35
    LLVMSIToFP = 36
    LLVMFPTrunc = 37
    LLVMFPExt = 38
    LLVMPtrToInt = 39
    LLVMIntToPtr = 40
    LLVMBitCast = 41
    LLVMAddrSpaceCast = 60
    LLVMICmp = 42
    LLVMFCmp = 43
    LLVMPHI = 44
    LLVMCall = 45
    LLVMSelect = 46
    LLVMUserOp1 = 47
    LLVMUserOp2 = 48
    LLVMVAArg = 49
    LLVMExtractElement = 50
    LLVMInsertElement = 51
    LLVMShuffleVector = 52
    LLVMExtractValue = 53
    LLVMInsertValue = 54
    LLVMFence = 55
    LLVMAtomicCmpXchg = 56
    LLVMAtomicRMW = 57
    LLVMResume = 58
    LLVMLandingPad = 59
    LLVMCleanupRet = 61
    LLVMCatchRet = 62
    LLVMCatchPad = 63
    LLVMCleanupPad = 64
    LLVMCatchSwitch = 65
end

@cenum LLVMTypeKind::UInt32 begin
    LLVMVoidTypeKind = 0
    LLVMHalfTypeKind = 1
    LLVMFloatTypeKind = 2
    LLVMDoubleTypeKind = 3
    LLVMX86_FP80TypeKind = 4
    LLVMFP128TypeKind = 5
    LLVMPPC_FP128TypeKind = 6
    LLVMLabelTypeKind = 7
    LLVMIntegerTypeKind = 8
    LLVMFunctionTypeKind = 9
    LLVMStructTypeKind = 10
    LLVMArrayTypeKind = 11
    LLVMPointerTypeKind = 12
    LLVMVectorTypeKind = 13
    LLVMMetadataTypeKind = 14
    LLVMX86_MMXTypeKind = 15
    LLVMTokenTypeKind = 16
end

@cenum LLVMLinkage::UInt32 begin
    LLVMExternalLinkage = 0
    LLVMAvailableExternallyLinkage = 1
    LLVMLinkOnceAnyLinkage = 2
    LLVMLinkOnceODRLinkage = 3
    LLVMLinkOnceODRAutoHideLinkage = 4
    LLVMWeakAnyLinkage = 5
    LLVMWeakODRLinkage = 6
    LLVMAppendingLinkage = 7
    LLVMInternalLinkage = 8
    LLVMPrivateLinkage = 9
    LLVMDLLImportLinkage = 10
    LLVMDLLExportLinkage = 11
    LLVMExternalWeakLinkage = 12
    LLVMGhostLinkage = 13
    LLVMCommonLinkage = 14
    LLVMLinkerPrivateLinkage = 15
    LLVMLinkerPrivateWeakLinkage = 16
end

@cenum LLVMVisibility::UInt32 begin
    LLVMDefaultVisibility = 0
    LLVMHiddenVisibility = 1
    LLVMProtectedVisibility = 2
end

@cenum LLVMDLLStorageClass::UInt32 begin
    LLVMDefaultStorageClass = 0
    LLVMDLLImportStorageClass = 1
    LLVMDLLExportStorageClass = 2
end

@cenum LLVMCallConv::UInt32 begin
    LLVMCCallConv = 0
    LLVMFastCallConv = 8
    LLVMColdCallConv = 9
    LLVMWebKitJSCallConv = 12
    LLVMAnyRegCallConv = 13
    LLVMX86StdcallCallConv = 64
    LLVMX86FastcallCallConv = 65
end

@cenum LLVMValueKind::UInt32 begin
    LLVMArgumentValueKind = 0
    LLVMBasicBlockValueKind = 1
    LLVMMemoryUseValueKind = 2
    LLVMMemoryDefValueKind = 3
    LLVMMemoryPhiValueKind = 4
    LLVMFunctionValueKind = 5
    LLVMGlobalAliasValueKind = 6
    LLVMGlobalIFuncValueKind = 7
    LLVMGlobalVariableValueKind = 8
    LLVMBlockAddressValueKind = 9
    LLVMConstantExprValueKind = 10
    LLVMConstantArrayValueKind = 11
    LLVMConstantStructValueKind = 12
    LLVMConstantVectorValueKind = 13
    LLVMUndefValueValueKind = 14
    LLVMConstantAggregateZeroValueKind = 15
    LLVMConstantDataArrayValueKind = 16
    LLVMConstantDataVectorValueKind = 17
    LLVMConstantIntValueKind = 18
    LLVMConstantFPValueKind = 19
    LLVMConstantPointerNullValueKind = 20
    LLVMConstantTokenNoneValueKind = 21
    LLVMMetadataAsValueValueKind = 22
    LLVMInlineAsmValueKind = 23
    LLVMInstructionValueKind = 24
end

@cenum LLVMIntPredicate::UInt32 begin
    LLVMIntEQ = 32
    LLVMIntNE = 33
    LLVMIntUGT = 34
    LLVMIntUGE = 35
    LLVMIntULT = 36
    LLVMIntULE = 37
    LLVMIntSGT = 38
    LLVMIntSGE = 39
    LLVMIntSLT = 40
    LLVMIntSLE = 41
end

@cenum LLVMRealPredicate::UInt32 begin
    LLVMRealPredicateFalse = 0
    LLVMRealOEQ = 1
    LLVMRealOGT = 2
    LLVMRealOGE = 3
    LLVMRealOLT = 4
    LLVMRealOLE = 5
    LLVMRealONE = 6
    LLVMRealORD = 7
    LLVMRealUNO = 8
    LLVMRealUEQ = 9
    LLVMRealUGT = 10
    LLVMRealUGE = 11
    LLVMRealULT = 12
    LLVMRealULE = 13
    LLVMRealUNE = 14
    LLVMRealPredicateTrue = 15
end

@cenum LLVMLandingPadClauseTy::UInt32 begin
    LLVMLandingPadCatch = 0
    LLVMLandingPadFilter = 1
end

@cenum LLVMThreadLocalMode::UInt32 begin
    LLVMNotThreadLocal = 0
    LLVMGeneralDynamicTLSModel = 1
    LLVMLocalDynamicTLSModel = 2
    LLVMInitialExecTLSModel = 3
    LLVMLocalExecTLSModel = 4
end

@cenum LLVMAtomicOrdering::UInt32 begin
    LLVMAtomicOrderingNotAtomic = 0
    LLVMAtomicOrderingUnordered = 1
    LLVMAtomicOrderingMonotonic = 2
    LLVMAtomicOrderingAcquire = 4
    LLVMAtomicOrderingRelease = 5
    LLVMAtomicOrderingAcquireRelease = 6
    LLVMAtomicOrderingSequentiallyConsistent = 7
end

@cenum LLVMAtomicRMWBinOp::UInt32 begin
    LLVMAtomicRMWBinOpXchg = 0
    LLVMAtomicRMWBinOpAdd = 1
    LLVMAtomicRMWBinOpSub = 2
    LLVMAtomicRMWBinOpAnd = 3
    LLVMAtomicRMWBinOpNand = 4
    LLVMAtomicRMWBinOpOr = 5
    LLVMAtomicRMWBinOpXor = 6
    LLVMAtomicRMWBinOpMax = 7
    LLVMAtomicRMWBinOpMin = 8
    LLVMAtomicRMWBinOpUMax = 9
    LLVMAtomicRMWBinOpUMin = 10
end

@cenum LLVMDiagnosticSeverity::UInt32 begin
    LLVMDSError = 0
    LLVMDSWarning = 1
    LLVMDSRemark = 2
    LLVMDSNote = 3
end

# the following isn't defined as a regular enum in llvm-c, so Clang.jl fails to pick it up
@cenum LLVMAttributeIndex::UInt32 begin
    LLVMAttributeReturnIndex = 0
    LLVMAttributeFunctionIndex = reinterpret(UInt32, Int32(-1))
end

const LLVMDiagnosticHandler = Ptr{Cvoid}
const LLVMYieldCallback = Ptr{Cvoid}

@cenum LLVMDIFlags::UInt32 begin
    LLVMDIFlagZero = 0
    LLVMDIFlagPrivate = 1
    LLVMDIFlagProtected = 2
    LLVMDIFlagPublic = 3
    LLVMDIFlagFwdDecl = 4
    LLVMDIFlagAppleBlock = 8
    LLVMDIFlagBlockByrefStruct = 16
    LLVMDIFlagVirtual = 32
    LLVMDIFlagArtificial = 64
    LLVMDIFlagExplicit = 128
    LLVMDIFlagPrototyped = 256
    LLVMDIFlagObjcClassComplete = 512
    LLVMDIFlagObjectPointer = 1024
    LLVMDIFlagVector = 2048
    LLVMDIFlagStaticMember = 4096
    LLVMDIFlagLValueReference = 8192
    LLVMDIFlagRValueReference = 16384
    LLVMDIFlagReserved = 32768
    LLVMDIFlagSingleInheritance = 65536
    LLVMDIFlagMultipleInheritance = 131072
    LLVMDIFlagVirtualInheritance = 196608
    LLVMDIFlagIntroducedVirtual = 262144
    LLVMDIFlagBitField = 524288
    LLVMDIFlagNoReturn = 1048576
    LLVMDIFlagMainSubprogram = 2097152
    LLVMDIFlagIndirectVirtualBase = 36
    LLVMDIFlagAccessibility = 3
    LLVMDIFlagPtrToMemberRep = 196608
end

@cenum LLVMDWARFSourceLanguage::UInt32 begin
    LLVMDWARFSourceLanguageC89 = 0
    LLVMDWARFSourceLanguageC = 1
    LLVMDWARFSourceLanguageAda83 = 2
    LLVMDWARFSourceLanguageC_plus_plus = 3
    LLVMDWARFSourceLanguageCobol74 = 4
    LLVMDWARFSourceLanguageCobol85 = 5
    LLVMDWARFSourceLanguageFortran77 = 6
    LLVMDWARFSourceLanguageFortran90 = 7
    LLVMDWARFSourceLanguagePascal83 = 8
    LLVMDWARFSourceLanguageModula2 = 9
    LLVMDWARFSourceLanguageJava = 10
    LLVMDWARFSourceLanguageC99 = 11
    LLVMDWARFSourceLanguageAda95 = 12
    LLVMDWARFSourceLanguageFortran95 = 13
    LLVMDWARFSourceLanguagePLI = 14
    LLVMDWARFSourceLanguageObjC = 15
    LLVMDWARFSourceLanguageObjC_plus_plus = 16
    LLVMDWARFSourceLanguageUPC = 17
    LLVMDWARFSourceLanguageD = 18
    LLVMDWARFSourceLanguagePython = 19
    LLVMDWARFSourceLanguageOpenCL = 20
    LLVMDWARFSourceLanguageGo = 21
    LLVMDWARFSourceLanguageModula3 = 22
    LLVMDWARFSourceLanguageHaskell = 23
    LLVMDWARFSourceLanguageC_plus_plus_03 = 24
    LLVMDWARFSourceLanguageC_plus_plus_11 = 25
    LLVMDWARFSourceLanguageOCaml = 26
    LLVMDWARFSourceLanguageRust = 27
    LLVMDWARFSourceLanguageC11 = 28
    LLVMDWARFSourceLanguageSwift = 29
    LLVMDWARFSourceLanguageJulia = 30
    LLVMDWARFSourceLanguageDylan = 31
    LLVMDWARFSourceLanguageC_plus_plus_14 = 32
    LLVMDWARFSourceLanguageFortran03 = 33
    LLVMDWARFSourceLanguageFortran08 = 34
    LLVMDWARFSourceLanguageRenderScript = 35
    LLVMDWARFSourceLanguageBLISS = 36
    LLVMDWARFSourceLanguageMips_Assembler = 37
    LLVMDWARFSourceLanguageGOOGLE_RenderScript = 38
    LLVMDWARFSourceLanguageBORLAND_Delphi = 39
end

@cenum LLVMDWARFEmissionKind::UInt32 begin
    LLVMDWARFEmissionNone = 0
    LLVMDWARFEmissionFull = 1
    LLVMDWARFEmissionLineTablesOnly = 2
end

const LLVMDisassembler_VariantKind_None = 0
const LLVMDisassembler_VariantKind_ARM_HI16 = 1
const LLVMDisassembler_VariantKind_ARM_LO16 = 2
const LLVMDisassembler_VariantKind_ARM64_PAGE = 1
const LLVMDisassembler_VariantKind_ARM64_PAGEOFF = 2
const LLVMDisassembler_VariantKind_ARM64_GOTPAGE = 3
const LLVMDisassembler_VariantKind_ARM64_GOTPAGEOFF = 4
const LLVMDisassembler_VariantKind_ARM64_TLVP = 5
const LLVMDisassembler_VariantKind_ARM64_TLVOFF = 6
const LLVMDisassembler_ReferenceType_InOut_None = 0
const LLVMDisassembler_ReferenceType_In_Branch = 1
const LLVMDisassembler_ReferenceType_In_PCrel_Load = 2
const LLVMDisassembler_ReferenceType_In_ARM64_ADRP = 0x0000000100000001
const LLVMDisassembler_ReferenceType_In_ARM64_ADDXri = 0x0000000100000002
const LLVMDisassembler_ReferenceType_In_ARM64_LDRXui = 0x0000000100000003
const LLVMDisassembler_ReferenceType_In_ARM64_LDRXl = 0x0000000100000004
const LLVMDisassembler_ReferenceType_In_ARM64_ADR = 0x0000000100000005
const LLVMDisassembler_ReferenceType_Out_SymbolStub = 1
const LLVMDisassembler_ReferenceType_Out_LitPool_SymAddr = 2
const LLVMDisassembler_ReferenceType_Out_LitPool_CstrAddr = 3
const LLVMDisassembler_ReferenceType_Out_Objc_CFString_Ref = 4
const LLVMDisassembler_ReferenceType_Out_Objc_Message = 5
const LLVMDisassembler_ReferenceType_Out_Objc_Message_Ref = 6
const LLVMDisassembler_ReferenceType_Out_Objc_Selector_Ref = 7
const LLVMDisassembler_ReferenceType_Out_Objc_Class_Ref = 8
const LLVMDisassembler_ReferenceType_DeMangled_Name = 9
const LLVMDisassembler_Option_UseMarkup = 1
const LLVMDisassembler_Option_PrintImmHex = 2
const LLVMDisassembler_Option_AsmPrinterVariant = 4
const LLVMDisassembler_Option_SetInstrComments = 8
const LLVMDisassembler_Option_PrintLatency = 16
const LLVMDisasmContextRef = Ptr{Cvoid}
const LLVMOpInfoCallback = Ptr{Cvoid}

struct LLVMOpInfoSymbol1
    Present::UInt64
    Name::Cstring
    Value::UInt64
end

struct LLVMOpInfo1
    AddSymbol::LLVMOpInfoSymbol1
    SubtractSymbol::LLVMOpInfoSymbol1
    Value::UInt64
    VariantKind::UInt64
end

const LLVMSymbolLookupCallback = Ptr{Cvoid}

@cenum LLVMByteOrdering::UInt32 begin
    LLVMBigEndian = 0
    LLVMLittleEndian = 1
end

struct LLVMOpaqueTargetData end
const LLVMTargetDataRef = Ptr{LLVMOpaqueTargetData}
struct LLVMOpaqueTargetLibraryInfotData end
const LLVMTargetLibraryInfoRef = Ptr{LLVMOpaqueTargetLibraryInfotData}
struct LLVMOpaqueTargetMachine end
const LLVMTargetMachineRef = Ptr{LLVMOpaqueTargetMachine}
const LLVMTarget = Cvoid
const LLVMTargetRef = Ptr{LLVMTarget}

@cenum LLVMCodeGenOptLevel::UInt32 begin
    LLVMCodeGenLevelNone = 0
    LLVMCodeGenLevelLess = 1
    LLVMCodeGenLevelDefault = 2
    LLVMCodeGenLevelAggressive = 3
end

@cenum LLVMRelocMode::UInt32 begin
    LLVMRelocDefault = 0
    LLVMRelocStatic = 1
    LLVMRelocPIC = 2
    LLVMRelocDynamicNoPic = 3
end

@cenum LLVMCodeModel::UInt32 begin
    LLVMCodeModelDefault = 0
    LLVMCodeModelJITDefault = 1
    LLVMCodeModelSmall = 2
    LLVMCodeModelKernel = 3
    LLVMCodeModelMedium = 4
    LLVMCodeModelLarge = 5
end

@cenum LLVMCodeGenFileType::UInt32 begin
    LLVMAssemblyFile = 0
    LLVMObjectFile = 1
end

struct LLVMOpaqueGenericValue end
const LLVMGenericValueRef = Ptr{LLVMOpaqueGenericValue}
struct LLVMOpaqueExecutionEngine end
const LLVMExecutionEngineRef = Ptr{LLVMOpaqueExecutionEngine}
struct LLVMOpaqueMCJITMemoryManager end
const LLVMMCJITMemoryManagerRef = Ptr{LLVMOpaqueMCJITMemoryManager}

struct LLVMMCJITCompilerOptions
    OptLevel::UInt32
    CodeModel::LLVMCodeModel
    NoFramePointerElim::LLVMBool
    EnableFastISel::LLVMBool
    MCJMM::LLVMMCJITMemoryManagerRef
end

const LLVMMemoryManagerAllocateCodeSectionCallback = Ptr{Cvoid}
const LLVMMemoryManagerAllocateDataSectionCallback = Ptr{Cvoid}
const LLVMMemoryManagerFinalizeMemoryCallback = Ptr{Cvoid}
const LLVMMemoryManagerDestroyCallback = Ptr{Cvoid}
const llvm_lto_t = Ptr{Cvoid}

@cenum llvm_lto_status::UInt32 begin
    LLVM_LTO_UNKNOWN = 0
    LLVM_LTO_OPT_SUCCESS = 1
    LLVM_LTO_READ_SUCCESS = 2
    LLVM_LTO_READ_FAILURE = 3
    LLVM_LTO_WRITE_FAILURE = 4
    LLVM_LTO_NO_TARGET = 5
    LLVM_LTO_NO_WORK = 6
    LLVM_LTO_MODULE_MERGE_FAILURE = 7
    LLVM_LTO_ASM_FAILURE = 8
    LLVM_LTO_NULL_OBJECT = 9
end

const llvm_lto_status_t = llvm_lto_status

@cenum LLVMLinkerMode::UInt32 begin
    LLVMLinkerDestroySource = 0
    LLVMLinkerPreserveSource_Removed = 1
end

struct LLVMOpaqueObjectFile end
const LLVMObjectFileRef = Ptr{LLVMOpaqueObjectFile}
struct LLVMOpaqueSectionIterator end
const LLVMSectionIteratorRef = Ptr{LLVMOpaqueSectionIterator}
struct LLVMOpaqueSymbolIterator end
const LLVMSymbolIteratorRef = Ptr{LLVMOpaqueSymbolIterator}
struct LLVMOpaqueRelocationIterator end
const LLVMRelocationIteratorRef = Ptr{LLVMOpaqueRelocationIterator}
struct LLVMOpaqueSharedModule end
const LLVMSharedModuleRef = Ptr{LLVMOpaqueSharedModule}
const LLVMOrcOpaqueJITStack = Cvoid
const LLVMOrcJITStackRef = Ptr{LLVMOrcOpaqueJITStack}
const LLVMOrcModuleHandle = UInt32
const LLVMOrcTargetAddress = UInt64
const LLVMOrcSymbolResolverFn = Ptr{Cvoid}
const LLVMOrcLazyCompileCallbackFn = Ptr{Cvoid}

@cenum LLVMOrcErrorCode::UInt32 begin
    LLVMOrcErrSuccess = 0
    LLVMOrcErrGeneric = 1
end

const LTO_API_VERSION = 21
const lto_bool_t = Bool

@cenum lto_symbol_attributes::UInt32 begin
    LTO_SYMBOL_ALIGNMENT_MASK = 31
    LTO_SYMBOL_PERMISSIONS_MASK = 224
    LTO_SYMBOL_PERMISSIONS_CODE = 160
    LTO_SYMBOL_PERMISSIONS_DATA = 192
    LTO_SYMBOL_PERMISSIONS_RODATA = 128
    LTO_SYMBOL_DEFINITION_MASK = 1792
    LTO_SYMBOL_DEFINITION_REGULAR = 256
    LTO_SYMBOL_DEFINITION_TENTATIVE = 512
    LTO_SYMBOL_DEFINITION_WEAK = 768
    LTO_SYMBOL_DEFINITION_UNDEFINED = 1024
    LTO_SYMBOL_DEFINITION_WEAKUNDEF = 1280
    LTO_SYMBOL_SCOPE_MASK = 14336
    LTO_SYMBOL_SCOPE_INTERNAL = 2048
    LTO_SYMBOL_SCOPE_HIDDEN = 4096
    LTO_SYMBOL_SCOPE_PROTECTED = 8192
    LTO_SYMBOL_SCOPE_DEFAULT = 6144
    LTO_SYMBOL_SCOPE_DEFAULT_CAN_BE_HIDDEN = 10240
    LTO_SYMBOL_COMDAT = 16384
    LTO_SYMBOL_ALIAS = 32768
end

@cenum lto_debug_model::UInt32 begin
    LTO_DEBUG_MODEL_NONE = 0
    LTO_DEBUG_MODEL_DWARF = 1
end

@cenum lto_codegen_model::UInt32 begin
    LTO_CODEGEN_PIC_MODEL_STATIC = 0
    LTO_CODEGEN_PIC_MODEL_DYNAMIC = 1
    LTO_CODEGEN_PIC_MODEL_DYNAMIC_NO_PIC = 2
    LTO_CODEGEN_PIC_MODEL_DEFAULT = 3
end

struct LLVMOpaqueLTOModule end
const lto_module_t = Ptr{LLVMOpaqueLTOModule}
struct LLVMOpaqueLTOCodeGenerator end
const lto_code_gen_t = Ptr{LLVMOpaqueLTOCodeGenerator}
struct LLVMOpaqueThinLTOCodeGenerator end
const thinlto_code_gen_t = Ptr{LLVMOpaqueThinLTOCodeGenerator}

@cenum lto_codegen_diagnostic_severity_t::UInt32 begin
    LTO_DS_ERROR = 0
    LTO_DS_WARNING = 1
    LTO_DS_REMARK = 3
    LTO_DS_NOTE = 2
end

const lto_diagnostic_handler_t = Ptr{Cvoid}

struct LTOObjectBuffer
    Buffer::Cstring
    Size::Csize_t
end

struct LLVMOpaquePassManagerBuilder end
const LLVMPassManagerBuilderRef = Ptr{LLVMOpaquePassManagerBuilder}


text:
#####
##### Shortcut for truncating uniform distributions.
#####

truncated(d::Uniform, l::T, u::T) where {T <: Real} = Uniform(max(l, d.a), min(u, d.b))


text:
using NearestNeighbors
using NearestNeighborDescent
using Distances
using MLDatasets
using Statistics

recall(true_ids, cand_ids) = mean(_recall(trues, cands) for (trues, cands) in zip(true_ids, cand_ids))
_recall(true_ids, cand_ids) = length(intersect(true_ids, cand_ids))/length(true_ids)

mnist_x = MNIST.convert2features(MNIST.traintensor(Float64))
fmnist_x = FashionMNIST.convert2features(FashionMNIST.traintensor(Float64))

mnist_y = MNIST.convert2features(MNIST.testtensor(Float64))
fmnist_y = FashionMNIST.convert2features(FashionMNIST.testtensor(Float64))

mnist_tree = KDTree(mnist_y; leafsize=50)
fmnist_tree = KDTree(fmnist_y; leafsize=50)

mnist_graph = nndescent(mnist_y, 20, Euclidean())
fmnist_graph = nndescent(fmnist_y, 20, Euclidean())

recalls = []

for (name, graph, tree, data) in [(:mnist, mnist_graph, mnist_tree, mnist_y), (:fmnist, fmnist_graph, fmnist_tree, fmnist_y)]

    nn_ids, _ = knn(tree, data, 21, true)
    true_ids = [v[2:end] for v in nn_ids]

    cand_nn_ids, _ = knn_matrices(graph)
    cand_ids = collect(eachcol(cand_nn_ids))

    push!(recalls, name=>recall(true_ids, cand_ids))
end

recalls


text:
# This is essentially the definition of sampling from a dirac delta.
Base.rand(RNG::AbstractRNG, P::Gaussian{U, Zero}) where U = P.μ

function kalman_sample(m::Model, inputs, rng::AbstractRNG, start_state,
                       N::Integer=length(inputs))
    # This code was optimized in 0.6, but sampling doesn't usually have to be
    # hyper-efficient, and it would look cleaner with a loop.
    result = accumulate(1:N; init=(start_state, nothing)) do v, t
        state, _ = v
        next_state = transition_mat(m, inputs, t) * state +
            rand(rng, transition_noise(m, inputs, t))
        return (next_state,
                observation_mat(m, inputs, t) * next_state +
                rand(rng, observation_noise(m, inputs, t)))
    end
    return (map(first, result), map(last, result))
end

################################################################################
# sample_and_recover

""" A data structure containing the results of `sample_and_recover`. """
struct RecoveryResults
    true_model       # the true model, passed by the user
    estimated_model  # the model estimated from maximum likelihood
    true_samples     # the samples generated from the true model
    estimated_state  # the P(state|observation, estimated_model) distributions
    obs              # the sampled observations
    optim            # the result from calling `optimize` (contains convergence details)
end
parameter_accuracy_ratios(rr::RecoveryResults) =
    [f=>getfield(rr.estimated_model, f) ./ getfield(rr.true_model, f)
     for f in parameters(rr.estimated_model)]

function Base.show(io::IO, rr::RecoveryResults)
    println(io, "Ratio of estimated/true parameters (1.0 is best): ")
    for (f, ratio) in parameter_accuracy_ratios(rr)
        println(io, "  ",
                f, " => ", round.(ratio*1000)/1000)
    end
    # This functionality was neat, but requiring Plots is not nice.
    # show(io, MIME"text/html"(),
    #      plot_hidden_state(1:length(rr.obs), rr.estimated_state;
    #                        true_samples=rr.true_samples))
end

""" See if we can recover the model parameters _and_ the true parameters using
data generated from the model.

Concretely, we sample observations and hidden state from `true_model` for the
given `inputs`, then call `optimize` on `true_model * fuzz_factor`.

If `start_model` isn't specified, we start from a model in the neighborhood of
`true_model` (with `fuzz_factor ~= 1.0` controlling how far we start).

We return a `RecoveryResults` object. See its definition for details. """
function sample_and_recover(true_model::Model, inputs, N=nothing; rng=GLOBAL_RNG,
                            parameters_to_optimize=parameters(true_model),
                            fuzz_factor=exp.(randn(rng, length(get_params(true_model, parameters_to_optimize)))),
                            initial_state::Gaussian=initial_state(true_model),
                            start_model=nothing)
    rng = rng isa AbstractRNG ? rng : MersenneTwister(rng::Integer)
    inputs2 = get_inputs(true_model, inputs)
    N = N === nothing ? length(inputs2) : N
    true_samples, obs =
        kalman_sample(true_model, inputs2, rng, rand(rng, initial_state), N)
    log_likelihood(true_model, inputs2, obs)
    if start_model === nothing
        true_params = get_params(true_model, parameters_to_optimize)
        start_model = set_params(true_model, (true_params .* fuzz_factor),
                                 parameters_to_optimize)
    end
    (best_model, o) = optimize(start_model, inputs, obs;
                               initial_state=initial_state,
                               parameters_to_optimize=parameters_to_optimize)
    estimated_state = kalman_smoother(best_model, get_inputs(best_model, inputs), obs,
                                      initial_state=initial_state)
    return RecoveryResults(true_model, best_model, true_samples, estimated_state, obs, o)
end


text:
using Libdl
cd(dirname(@__FILE__))

if (!ispath("usr"))
    run(`mkdir usr`)
end
if (!ispath("usr/lib"))
    run(`mkdir usr/lib`)
end

if (!ispath("TSSM"))
    run(`git clone git://github.com/HaraldHofstaetter/TSSM.git`)
else
    cd("TSSM")
    run(`git pull`)
    cd("..")
end

@info "Building TSSM standard version"
if occursin("juliabox", read(`uname -a`, String))
    run(`make -C TSSM/OPENMP JULIABOX_BUILD=1 all`)
    run(`mv TSSM/OPENMP/libtssm.$(Libdl.dlext) usr/lib`)
    run(`make -C TSSM/OPENMP clean`)
else
    run(`make -C TSSM/OPENMP all`)
    run(`mv TSSM/OPENMP/libtssm.$(Libdl.dlext) usr/lib`)
end

if "TSSM_DEBUG" in keys(ENV) && ENV["TSSM_DEBUG"]=="1"
    @info "Building TSSM debug version"
    if occursin("juliabox", read(`uname -a`))
        run(`make -C TSSM/DEBUG JULIABOX_BUILD=1 all`)
        run(`mv TSSM/DEBUG/libtssm.$(Libdl.dlext) usr/lib/libtssm_debug.$(Libdl.dlext)`)
        #run(`make -C TSSM/DEBUG clean`)
    else
        run(`make -C TSSM/DEBUG all`)
        run(`mv TSSM/DEBUG/libtssm.$(Libdl.dlext) usr/lib/libtssm_debug.$(Libdl.dlext)`)
    end
end

use_Float128 = false
#try
#    using Quadmath
#    use_Float128 = true
#end

if use_Float128
    @info "Building TSSM quadprecision version"
    if (   !ispath("usr/lib/libfftw3q.so.3.4.4")
        || !ispath("usr/lib/libfftw3q_omp.so.3.4.4"))
        if (!ispath("fftw-3.3.4"))
            download("http://www.fftw.org/fftw-3.3.4.tar.gz","fftw3.tgz")
            run(`tar xzf fftw3.tgz`)
            rm("fftw3.tgz")
        end
        install_dir = pwd()
        cd("fftw-3.3.4")
        run(`./configure --prefix=$install_dir/usr --enable-quad-precision --enable-openmp --enable-threads --enable-shared`)
        run(`make install`)
        cd("..")
        rm("fftw-3.3.4", recursive=true)
    end

    if searchindex(readstring(`uname -a`), "juliabox")>0
        run(`make -C TSSM/QUADPREC_OPENMP JULIABOX_BUILD=1 all`)
        run(`mv TSSM/QUADPREC_OPENMP/libtssmq.$(Libdl.dlext) usr/lib`)
        run(`make -C TSSM/QUADPREC_OPENMP clean`)
    else
        run(`make -C TSSM/QUADPREC_OPENMP all`)
        run(`mv TSSM/QUADPREC_OPENMP/libtssmq.$(Libdl.dlext) usr/lib`)
    end
end


text:
export nrmsd

"""
    nrmsd(I,Ireco)

computes the normalized root mean squared error of the image `Ireco`
with respect to the image `I`.
"""
function nrmsd(I,Ireco)
  N = length(I)

  # This is a little trick. We usually are not interested in simple scalings
  # and therefore "calibrate" them away
  alpha = (dot(vec(I),vec(Ireco))+dot(vec(Ireco),vec(I))) /
          (2*dot(vec(Ireco),vec(Ireco)))
  Ireco[:] .*= alpha

  RMS =  1.0/sqrt(N)*norm(vec(I)-vec(Ireco))
  NRMS = RMS/(maximum(abs.(I))-minimum(abs.(I)) )
  return NRMS
end


text:
# IupFill: Example
#  Uses the Fill element to horizontally centralize a button and to justify it to the left and right.

module fill_

using IUP

export
    fill

#------------------------example html/examples/C/fill.c ----------------------------------
function fill()
    IupOpen()    #Initializes IUP

    # Creates frame with left aligned button */
    frame_left = IupFrame(
        IupHbox(IupButton("Ok" , ""), IupFill())
    )

    # Sets frame's title */
    IupSetAttribute( frame_left, "TITLE", "Left aligned" );

    # Creates frame with centered button
    frame_center = IupFrame(
        IupHbox(IupFill(), IupButton("Ok", ""), IupFill())
    )

    # Sets frame's title
    IupSetAttribute(frame_center, "TITLE", "Centered")

    # Creates frame with right aligned button
    frame_right = IupFrame(IupHbox(IupFill(),IupButton("Ok", "")))

    # Sets frame's title
    IupSetAttribute(frame_right, "TITLE", "Right aligned")

    # Creates dialog with these three frames
    dialog = IupDialog(IupVbox(frame_left, frame_center, frame_right))

    # Sets dialog's size and title
    IupSetAttributes(dialog, "SIZE=120, TITLE=IupFill")

    IupShow(dialog)                    # Shows dialog in the center of the screen
    IupMainLoop()                      # Initializes IUP main loop
    IupClose()                         # And close it when ready
end

end        # module


text:
#####
##### _compare! - CPU Edition
#####

function Profiler._compare!(
        stats,
        f,
        opt::ILPOptimizer,
        backend::nGraph.Backend{nGraph.CPU};
        skip_run = false,
        cache = nothing,
        kw...,
    )

    fex, frame, _metadata = factory(backend, f, opt; cache = cache, kw...)
    GC.gc()

    # Reprofile to capture the move nodes.
    data = profile(fex; cache = cache)

    # Unpack some of the return values for logging
    creation_times = _metadata[:creation_times]
    optimization_times = _metadata[:optimization_times]

    # Get the predicted run time and then the actual run time
    nt = Dict(
        :creation_times => creation_times,
        :optimization_times => optimization_times,
        :predicted_runtime => predict(frame),
        :dram_limit => maxlimit(frame.modeltype),
        :tensor_size_map => Dict(nGraph.name(t) => sizeof(t) for t in tensors(data)),
        :config_map => Dict(nGraph.name(n) => getconfig(nGraph.Node(unx(n))) for n in nodes(data)),
        :ratio => getratio(opt),
    )

    nt_new = Dict(
        # Some statistics on nodes and tensors

        # Number of move nodes plus bytes moved around
        :num_move_nodes => count(_move_filter(), nodes(data)),
        :num_pmem_move_nodes => count(_move_filter(PMEM), nodes(data)),
        :num_dram_move_nodes => count(_move_filter(DRAM), nodes(data)),

        :bytes_moved => _count(inputs, sizeof, data; filt = _move_filter()),
        :bytes_moved_pmem => _count(inputs, sizeof, data; filt = _move_filter(PMEM)),
        :bytes_moved_dram => _count(inputs, sizeof, data; filt = _move_filter(DRAM)),

        :num_async_move_nodes => count(_async_filter(), nodes(data)),
        :num_pmem_async_move_nodes => count(_async_filter(PMEM), nodes(data)),
        :num_dram_async_move_nodes => count(_async_filter(DRAM), nodes(data)),

        :bytes_async_moved => _count(inputs, sizeof, data; filt = _async_filter()),
        :bytes_async_moved_pmem => _count(inputs, sizeof, data; filt = _async_filter(PMEM)),
        :bytes_async_moved_dram => _count(inputs, sizeof, data; filt = _async_filter(DRAM)),

        # Total number of kernels
        :num_kernels => count(hasprofile, nodes(data)),
        :num_input_tensors => _count(inputs, data; filt = hasprofile),
        :num_output_tensors => _count(outputs, data; filt = hasprofile),

        :num_dram_input_tensors => _count(
            x -> filter(!nGraph.is_persistent, inputs(x)),
            data; filt = hasprofile
        ),
        :num_dram_output_tensors => _count(
            x -> filter(!nGraph.is_persistent, outputs(x)),
            data; filt = hasprofile
        ),

        # Get the sizes of the input and output tensors
        :bytes_input_tensors => _count(inputs, sizeof, data; filt = hasprofile),
        :bytes_output_tensors => _count(outputs, sizeof, data; filt = hasprofile),

        :bytes_dram_input_tensors => _count(
            x -> filter(!nGraph.is_persistent, inputs(x)),
            sizeof,
            data;
            filt = hasprofile
        ),
        :bytes_dram_output_tensors => _count(
            x -> filter(!nGraph.is_persistent, outputs(x)),
            sizeof,
            data;
            filt = hasprofile
        ),

        # Info on global allocations
        :dram_alloc_size => nGraph.get_temporary_pool_size(fex.ex.ngraph_function),
        :pmem_alloc_size => nGraph.get_pmem_pool_size(fex.ex.ngraph_function),
        :move_time => estimate_move_time(fex, frame),
    )
    nt = merge(nt, nt_new)

    if !skip_run
        nt_new = Dict(
            :actual_runtime => gettime(fex),
            :kernel_times => nGraph.get_performance(fex.ex)
        )
        nt = merge(nt, nt_new)
    end

    push!(stats.runs, nt)

    return nothing
end

#####
##### _compare! - GPU Edition
#####
function Profiler._compare!(
        stats,
        f,
        opt::ILPOptimizer,
        backend::nGraph.Backend{nGraph.GPU};
        cache = nothing,
        kw...
    )

    # If we're just profiling, `factory` can return `nothing`.
    #
    # Hwere, we just perform a quick check for that actions so we don't get an indexing error.
    x = factory(backend, f, opt; cache = cache, kw...)
    isnothing(x) && return nothing

    fex, frame = x

    actual_time = gettime(fex)
    oracle_time = Profiler.fastest_time(fex, frame)

    @show oracle_time

    GC.gc()
    data = profile(fex; cache = cache)

    # Get the predicted run time and then the actual run time
    nt = Dict(
        :predicted_runtime => predict(frame),
        :dram_limit => maxlimit(frame.modeltype),
        :tensor_size_map => Dict(nGraph.name(t) => sizeof(t) for t in tensors(data)),

        # Number of move nodes plus bytes moved around
        :num_move_nodes => count(_move_filter(), nodes(data)),
        :num_pmem_move_nodes => count(_move_filter(PMEM), nodes(data)),
        :num_dram_move_nodes => count(_move_filter(DRAM), nodes(data)),

        :bytes_moved => _count(inputs, sizeof, data; filt = _move_filter()),
        :bytes_moved_pmem => _count(inputs, sizeof, data; filt = _move_filter(PMEM)),
        :bytes_moved_dram => _count(inputs, sizeof, data; filt = _move_filter(DRAM)),

        :num_async_move_nodes => count(_async_filter(), nodes(data)),
        :num_pmem_async_move_nodes => count(_async_filter(PMEM), nodes(data)),
        :num_dram_async_move_nodes => count(_async_filter(DRAM), nodes(data)),

        :bytes_async_moved => _count(inputs, sizeof, data; filt = _async_filter()),
        :bytes_async_moved_pmem => _count(inputs, sizeof, data; filt = _async_filter(PMEM)),
        :bytes_async_moved_dram => _count(inputs, sizeof, data; filt = _async_filter(DRAM)),

        # Total number of kernels
        :num_kernels => count(hasprofile, nodes(data)),
        :num_input_tensors => _count(inputs, data; filt = hasprofile),
        :num_output_tensors => _count(outputs, data; filt = hasprofile),

        :num_dram_input_tensors => _count(
            x -> filter(!nGraph.is_persistent, inputs(x)),
            data; filt = hasprofile
        ),
        :num_dram_output_tensors => _count(
            x -> filter(!nGraph.is_persistent, outputs(x)),
            data; filt = hasprofile
        ),

        # Get the sizes of the input and output tensors
        :bytes_input_tensors => _count(inputs, sizeof, data; filt = hasprofile),
        :bytes_output_tensors => _count(outputs, sizeof, data; filt = hasprofile),

        :bytes_dram_input_tensors => _count(
            x -> filter(!nGraph.is_persistent, inputs(x)),
            sizeof,
            data;
            filt = hasprofile
        ),
        :bytes_dram_output_tensors => _count(
            x -> filter(!nGraph.is_persistent, outputs(x)),
            sizeof,
            data;
            filt = hasprofile
        ),

        # Info on global allocations
        :dram_alloc_size => nGraph.get_temporary_pool_size(fex.ex.ngraph_function),
        :pmem_alloc_size => nGraph.get_pmem_pool_size(fex.ex.ngraph_function),
        :move_time => estimate_move_time(fex, frame),

        # Timing Breakdowns
        :actual_runtime => actual_time,
        :oracle_time => oracle_time,
        #:kernel_times => read_timing_data(fex.ex.ngraph_function)
    )

    push!(stats.runs, nt)
    return nothing
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
# ---
# title: 1449. Form Largest Integer With Digits That Add up to Target
# id: problem1449
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/form-largest-integer-with-digits-that-add-up-to-target/description/>
# hidden: true
# ---
#
# Given an array of integers `cost` and an integer `target`. Return the
# **maximum** integer you can paint under the following rules:
#
#   * The cost of painting a digit (i+1) is given by `cost[i]` (0 indexed).
#   * The total cost used must be equal to `target`.
#   * Integer does not have digits 0.
#
# Since the answer may be too large, return it as string.
#
# If there is no way to paint any integer given the condition, return "0".
#
#
#
# **Example 1:**
#
#
#
#     Input: cost = [4,3,2,5,6,7,2,5,5], target = 9
#     Output: "7772"
#     Explanation: The cost to paint the digit '7' is 2, and the digit '2' is 3. Then cost("7772") = 2*3+ 3*1 = 9. You could also paint "977", but "7772" is the largest number.
#     **Digit    cost**
#       1  ->   4
#       2  ->   3
#       3  ->   2
#       4  ->   5
#       5  ->   6
#       6  ->   7
#       7  ->   2
#       8  ->   5
#       9  ->   5
#
#
# **Example 2:**
#
#
#
#     Input: cost = [7,6,5,5,5,6,8,7,8], target = 12
#     Output: "85"
#     Explanation: The cost to paint the digit '8' is 7, and the digit '5' is 5. Then cost("85") = 7 + 5 = 12.
#
#
# **Example 3:**
#
#
#
#     Input: cost = [2,4,6,2,4,6,4,4,4], target = 5
#     Output: "0"
#     Explanation: It's not possible to paint any integer with total cost equal to target.
#
#
# **Example 4:**
#
#
#
#     Input: cost = [6,10,15,40,40,40,40,40,40], target = 47
#     Output: "32211"
#
#
#
#
# **Constraints:**
#
#   * `cost.length == 9`
#   * `1 <= cost[i] <= 5000`
#   * `1 <= target <= 5000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
mutable struct DocumentTermMatrix{T}
    dtm::SparseMatrixCSC{Int, Int}
    terms::Vector{T}
    column_indices::Dict{T, Int}
end

"""
    columnindices(terms::Vector{String})

Creates a column index lookup dictionary from a vector of terms.
"""
function columnindices(terms::Vector{T}) where T
    column_indices = Dict{T, Int}()
    for i in 1:length(terms)
        term = terms[i]
        column_indices[term] = i
    end
    column_indices
end

"""
    DocumentTermMatrix(crps::Corpus)
    DocumentTermMatrix(crps::Corpus, terms::Vector{String})
    DocumentTermMatrix(crps::Corpus, lex::AbstractDict)
    DocumentTermMatrix(dtm::SparseMatrixCSC{Int, Int},terms::Vector{String})

Represent documents as a matrix of word counts.

Allow us to apply linear algebra operations and statistical techniques.
Need to update lexicon before use.

# Examples
```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
                      StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> m = DocumentTermMatrix(crps)
A 2 X 6 DocumentTermMatrix

julia> m.dtm
2×6 SparseArrays.SparseMatrixCSC{Int64,Int64} with 10 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
  [1, 2]  =  2
  [2, 3]  =  2
  [1, 4]  =  1
  [2, 4]  =  1
  [1, 5]  =  1
  [2, 5]  =  1
  [1, 6]  =  1
  [2, 6]  =  1
```
"""
function DocumentTermMatrix(crps::Corpus, terms::Vector{T}) where T
    column_indices = columnindices(terms)

    m = length(crps)
    n = length(terms)

    rows = Array{Int}(undef, 0)
    columns = Array{Int}(undef, 0)
    values = Array{Int}(undef, 0)
    for i in 1:m
        doc = crps.documents[i]
        ngs = ngrams(doc)
        for ngram in keys(ngs)
            j = get(column_indices, ngram, 0)
            v = ngs[ngram]
            if j != 0
                push!(rows, i)
                push!(columns, j)
                push!(values, v)
            end
        end
    end
    if length(rows) > 0
        dtm = sparse(rows, columns, values, m, n)
    else
        dtm = spzeros(Int, m, n)
    end
    DocumentTermMatrix(dtm, terms, column_indices)
end
DocumentTermMatrix(crps::Corpus) = DocumentTermMatrix(crps, lexicon(crps))

DocumentTermMatrix(crps::Corpus, lex::AbstractDict) = DocumentTermMatrix(crps, sort(collect(keys(lex))))

DocumentTermMatrix(dtm::SparseMatrixCSC{Int, Int},terms::Vector{T}) where T = DocumentTermMatrix{T}(dtm, terms, columnindices(terms))

"""
    dtm(crps::Corpus)
    dtm(d::DocumentTermMatrix)
    dtm(d::DocumentTermMatrix, density::Symbol)

Creates a simple sparse matrix of DocumentTermMatrix object.

# Examples
```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
                      StringDocument("To become or not to become")])

julia> update_lexicon!(crps)

julia> dtm(DocumentTermMatrix(crps))
2×6 SparseArrays.SparseMatrixCSC{Int64,Int64} with 10 stored entries:
  [1, 1]  =  1
  [2, 1]  =  1
  [1, 2]  =  2
  [2, 3]  =  2
  [1, 4]  =  1
  [2, 4]  =  1
  [1, 5]  =  1
  [2, 5]  =  1
  [1, 6]  =  1
  [2, 6]  =  1

julia> dtm(DocumentTermMatrix(crps), :dense)
2×6 Array{Int64,2}:
 1  2  0  1  1  1
 1  0  2  1  1  1
```
"""
function dtm(d::DocumentTermMatrix, density::Symbol)
    if density == :sparse
        return d.dtm
    else
        return Matrix(d.dtm)
    end
end

function dtm(d::DocumentTermMatrix)
    return d.dtm
end

function dtm(crps::Corpus)
    dtm(DocumentTermMatrix(crps))
end

tdm(crps::DocumentTermMatrix, density::Symbol) = dtm(crps, density)' #'

tdm(crps::DocumentTermMatrix) = dtm(crps)' #'

tdm(crps::Corpus) = dtm(crps)' #'

##############################################################################
#
# Produce the signature of a DTM entry for a document
#
# TODO: Make this more efficieny by reusing column_indices
#
##############################################################################

function dtm_entries(d::AbstractDocument, lex::Dict{T, Int}) where T
    ngs = ngrams(d)
    indices = Array{Int}(undef, 0)
    values = Array{Int}(undef, 0)
    terms = sort(collect(keys(lex)))
    column_indices = columnindices(terms)

    for ngram in keys(ngs)
        if haskey(column_indices, ngram)
            push!(indices, column_indices[ngram])
            push!(values, ngs[ngram])
        end
    end
    return (indices, values)
end

"""
    dtv(d::AbstractDocument, lex::Dict{String, Int})

Produce a single row of a DocumentTermMatrix.

Individual documents do not have a lexicon associated with them,
we have to pass in a lexicon as an additional argument.

# Examples
```julia-repl
julia> dtv(crps[1], lexicon(crps))
1×6 Array{Int64,2}:
 1  2  0  1  1  1
```
"""
function dtv(d::AbstractDocument, lex::Dict{T, Int}) where T
    p = length(keys(lex))
    row = zeros(Int, 1, p)
    indices, values = dtm_entries(d, lex)
    for i in 1:length(indices)
        row[1, indices[i]] = values[i]
    end
    return row
end

function dtv(d::AbstractDocument)
    error("Cannot construct a DTV without a pre-existing lexicon")
end

##############################################################################
#
# The hash trick: use a hash function instead of a lexicon to determine the
# columns of a DocumentTermMatrix-like encoding of the data
#
##############################################################################
"""
    hash_dtv(d::AbstractDocument)
    hash_dtv(d::AbstractDocument, h::TextHashFunction)

Represents a document as a vector with N entries.

# Examples
```julia-repl
julia> crps = Corpus([StringDocument("To be or not to be"),
                      StringDocument("To become or not to become")])

julia> h = TextHashFunction(10)
TextHashFunction(hash, 10)

julia> hash_dtv(crps[1], h)
1×10 Array{Int64,2}:
 0  2  0  0  1  3  0  0  0  0

julia> hash_dtv(crps[1])
1×100 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
```
"""
function hash_dtv(d::AbstractDocument, h::TextHashFunction)
    p = cardinality(h)
    res = zeros(Int, 1, p)
    ngs = ngrams(d)
    for ng in keys(ngs)
        res[1, index_hash(ng, h)] += ngs[ng]
    end
    return res
end

hash_dtv(d::AbstractDocument) = hash_dtv(d, TextHashFunction())

"""
    hash_dtm(crps::Corpus)
    hash_dtm(crps::Corpus, h::TextHashFunction)

Represents a Corpus as a Matrix with N entries.
"""
function hash_dtm(crps::Corpus, h::TextHashFunction)
    n, p = length(crps), cardinality(h)
    res = zeros(Int, n, p)
    for i in 1:length(crps)
        doc = crps.documents[i]
        res[i, :] = hash_dtv(doc, h)
    end
    return res
end

hash_dtm(crps::Corpus) = hash_dtm(crps, hash_function(crps))

hash_tdm(crps::Corpus) = hash_dtm(crps)' #'

##############################################################################
#
# Produce entries for on-line analysis when DTM would not fit in memory
#
##############################################################################

mutable struct EachDTV
    crps::Corpus
end

start(edt::EachDTV) = 1

function next(edt::EachDTV, state::Int)
    return (dtv(edt.crps.documents[state], lexicon(edt.crps)), state + 1)
end

done(edt::EachDTV, state::Int) = state > length(edt.crps.documents)

mutable struct EachHashDTV
    crps::Corpus
end

start(edt::EachHashDTV) = 1

function next(edt::EachHashDTV, state::Int)
    (hash_dtv(edt.crps.documents[state]), state + 1)
end

done(edt::EachHashDTV, state::Int) = state > length(edt.crps.documents)

each_dtv(crps::Corpus) = EachDTV(crps)

each_hash_dtv(crps::Corpus) = EachHashDTV(crps)

##
## getindex() methods
##

Base.getindex(dtm::DocumentTermMatrix, k::AbstractString) = dtm.dtm[:, dtm.column_indices[k]]
Base.getindex(dtm::DocumentTermMatrix, i::Any) = dtm.dtm[i]
Base.getindex(dtm::DocumentTermMatrix, i::Any, j::Any) = dtm.dtm[i, j]


text:
# Julia program to find the sum of divisors of a number.

function sum_of_divisors(num)
    i = 1
    sum = 0
    while (i <= sqrt(num))
        if num % i == 0 && i == sqrt(num)
            sum = sum + i
        elseif num % i == 0
            sum = sum + i
            sum = sum + (num ÷ i)
        end
        i = i + 1
    end
    return sum
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
sum = sum_of_divisors(num)
print("The Sum of the Divisors of the given number is $sum.")

"""
Time Complexity - O(n^(0.5)), where 'n' given is the number
Space Complexity - O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE I

Enter the number: 256
The Sum of the Divisors of the given number is 511.

SAMPLE II

Enter the number: 5687
The Sum of the Divisors of the given number is 6384.
"""