text:
#
# lgemm.jl -
#
# Lazily Generalized Matrix-Matrix mutiplication.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2017-2020 Éric Thiébaut.
#

"""
# Lazily Generalized Matrix-Matrix mutiplication

```julia
lgemm([α=1,] [transA='N',] A, [transB='N',] B, Nc=2) -> C
```

yields `C = α*op(A)*op(B)` where `op(A)` is `A` if `transA` is `'N'`,
`transpose(A)` if `transA` is `'T'` and `adjoint(A)` if `transA` is `'C'` and
similarly for `op(B)` and `transB`.  The expression `op(A)*op(B)` is a
matrix-matrix multiplication but interpreted by grouping consecutive dimensions
of `op(A)` and `op(B)` (see [`lgemv`](@ref) for explanations).  Argument `Nc`
specifies the number of dimensions of the result.

The in-place version is called by:

```julia
lgemm!([α=1,] [transA='N',] A,  [transB='N',] B, [β=0,] C) -> C
```

and overwrites the contents of `C` with `α*op(A)*op(B) + β*C`.  Note that `C`
must not be aliased with `A` nor `B`.

The multipliers `α` and `β` must be both specified or omitted, they can be any
scalar numbers but are respectively converted to `promote_eltype(A,B)` and
`eltype(C)` which may throw an `InexactError` exception.

See also: [`lgemv`](@ref), [`LinearAlgebra.BLAS.gemm`](@ref),
[`LinearAlgebra.BLAS.gemm!`](@ref).

"""
function lgemm(α::Number,
               transA::Char,
               A::AbstractArray{<:Floats},
               transB::Char,
               B::AbstractArray{<:Floats},
               Nc::Integer=2)
    return _lgemm(Implementation(Val(:lgemm), α, transA, A, transB, B, Int(Nc)),
                  α, transA, A, transB, B, Int(Nc))
end

function lgemm!(α::Number,
                transA::Char,
                A::AbstractArray{<:Floats},
                transB::Char,
                B::AbstractArray{<:Floats},
                β::Number,
                C::AbstractArray{<:Floats})
    return _lgemm!(Implementation(Val(:lgemm), α, transA, A, transB, B, β, C),
                   α, transA, A, transB, B, β, C)
end

@doc @doc(lgemm) lgemm!

# Best implementations for lgemm and lgemm!

for (atyp, eltyp) in ((Real,   BlasReal),
                      (Number, BlasComplex))
    @eval begin
        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::Matrix{T},
                                transB::Char,
                                B::Matrix{T},
                                Nc::Int) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::Matrix{T},
                                transB::Char,
                                B::Matrix{T},
                                β::$atyp,
                                C::Matrix{T}) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::AbstractMatrix{T},
                                transB::Char,
                                B::AbstractMatrix{T},
                                Nc::Int) where {T<:$eltyp}
            return (is_flat_array(A, B) ? Blas() : Basic())
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::AbstractMatrix{T},
                                transB::Char,
                                B::AbstractMatrix{T},
                                β::$atyp,
                                C::AbstractMatrix{T}) where {T<:$eltyp}
            return (is_flat_array(A, B, C) ? Blas() : Basic())
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::Array{T},
                                transB::Char,
                                B::Array{T},
                                Nc::Int) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::Array{T},
                                transB::Char,
                                B::Array{T},
                                β::$atyp,
                                C::Array{T}) where {T<:$eltyp}
            return Blas()
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::AbstractArray{T},
                                transB::Char,
                                B::AbstractArray{T},
                                Nc::Int) where {T<:$eltyp}
            return (is_flat_array(A, B) ? Blas() : Basic())
        end

        function Implementation(::Val{:lgemm},
                                α::$atyp,
                                transA::Char,
                                A::AbstractArray{T},
                                transB::Char,
                                B::AbstractArray{T},
                                β::$atyp,
                                C::AbstractArray{T}) where {T<:$eltyp}
            return (is_flat_array(A, B, C) ? Blas() : Basic())
        end
    end
end

function Implementation(::Val{:lgemm},
                        α::Number,
                        transA::Char,
                        A::AbstractMatrix,
                        transB::Char,
                        B::AbstractMatrix,
                        Nc::Int)
    return Basic()
end

function Implementation(::Val{:lgemm},
                        α::Number,
                        transA::Char,
                        A::AbstractMatrix,
                        transB::Char,
                        B::AbstractMatrix,
                        β::Number,
                        C::AbstractMatrix)
    return Basic()
end

function Implementation(::Val{:lgemm},
                        α::Number,
                        transA::Char,
                        A::AbstractArray,
                        transB::Char,
                        B::AbstractArray,
                        Nc::Int)
    return (is_flat_array(A, B) ? Linear() : Generic())
end

function Implementation(::Val{:lgemm},
                        α::Number,
                        transA::Char,
                        A::AbstractArray,
                        transB::Char,
                        B::AbstractArray,
                        β::Number,
                        C::AbstractArray)
    return (is_flat_array(A, B, C) ? Linear() : Generic())
end

"""

```julia
_lgemm_ndims(Na, Nb, Nc) -> Ni, Nj, Nk
```

yields the number of consecutive indices involved in meta-indices `i`, `j` and
`k` for `lgemm` and `lgemm!` methods which compute:

```
C[i,j] = α⋅sum_k op(A)[i,k]*op(B)[k,j] + β⋅C[i,j]
```

Here `Na`, `Nb` and `Nc` are the respective number of dimensions of `A`, `B`
and `C`.  Calling this method also check the compatibility of the number of
dimensions.

"""
function _lgemm_ndims(Na::Int, Nb::Int, Nc::Int)
    #
    # The relations are:
    #
    #     Na = Ni + Nk         Ni = (Na + Nc - Nb)/2
    #     Nb = Nk + Nj   <==>  Nj = (Nc + Nb - Na)/2
    #     Nc = Ni + Nj         Nk = (Nb + Na - Nc)/2
    #
    # Note: Ni ≥ 1 and Nj ≥ 1 and Nk ≥ 1 implies that Na ≥ 2 and Nb ≥ 2
    # and Nc ≥ 2.
    #
    Li = Na + Nc - Nb
    Lj = Nc + Nb - Na
    Lk = Nb + Na - Nc
    if Li > 0 && iseven(Li) && Lj > 0 && iseven(Lj) && Lk > 0 && iseven(Lk)
        return Li >> 1, Lj >> 1, Lk >> 1
    end
    incompatible_dimensions()
end

# BLAS implementations for (generalized) matrices.
# Linear indexing is assumed (this must have been checked before).

function _lgemm(::Blas,
                α::Number,
                transA::Char,
                A::AbstractArray{T},
                transB::Char,
                B::AbstractArray{T},
                Nc::Int) where {T<:BlasFloat}
    m, n, p, shape = _lgemm_dims(transA, A, transB, B, Nc)
    return _blas_lgemm!(m, n, p, convert(T, α), transA, A,
                        transB, B, zero(T), Array{T}(undef, shape))
end

function _lgemm!(::Blas,
                 α::Number,
                 transA::Char,
                 A::AbstractArray{T},
                 transB::Char,
                 B::AbstractArray{T},
                 β::Number,
                 C::AbstractArray{T}) where {T<:BlasFloat}
    m, n, p = _lgemm_dims(transA, A, transB, B, C)
    return _blas_lgemm!(m, n, p, convert(T, α), transA, A,
                        transB, B, convert(T, β), C)
end

# Julia implementations for (generalized) matrices.
# Linear indexing is assumed (this must have been checked before).

function _lgemm(::Linear,
                α::Number,
                transA::Char,
                A::AbstractArray{<:Floats},
                transB::Char,
                B::AbstractArray{<:Floats},
                Nc::Int)
    # FIXME: check and improve type stability
    m, n, p, shape = _lgemm_dims(transA, A, transB, B, Nc)
    T = _lgemm_type(α, A, B)
    return _linear_lgemm!(m, n, p,
                          promote_multiplier(α, A, B), transA, A,
                          transB, B,
                          promote_multiplier(0, T),
                          Array{T}(undef, shape))
end

function _lgemm!(::Linear,
                 α::Number,
                 transA::Char,
                 A::AbstractArray{<:Floats},
                 transB::Char,
                 B::AbstractArray{<:Floats},
                 β::Number,
                 C::AbstractArray{<:Floats})
    m, n, p = _lgemm_dims(transA, A, transB, B, C)
    return _linear_lgemm!(m, n, p,
                          promote_multiplier(α, A, B), transA, A,
                          transB, B,
                          promote_multiplier(β, C), C)
end

# Julia implementations for any kind of abstract matrices.

function _lgemm(::Basic,
                α::Number,
                transA::Char,
                A::AbstractMatrix{<:Floats},
                transB::Char,
                B::AbstractMatrix{<:Floats},
                Nc::Int)
    I, J, K = _lgemm_indices(transA, A, transB, B, Nc)
    T = _lgemm_type(α, A, B)
    return _generic_lgemm!(I, J, K,
                           promote_multiplier(α, A, B), transA, A,
                           transB, B,
                           promote_multiplier(0, T),
                           similar(Array{T}, (I, J)))
end

function _lgemm!(::Basic,
                 α::Number,
                 transA::Char,
                 A::AbstractMatrix{<:Floats},
                 transB::Char,
                 B::AbstractMatrix{<:Floats},
                 β::Number,
                 C::AbstractMatrix{<:Floats})
    I, J, K = _lgemm_indices(transA, A, transB, B, C)
    return _generic_lgemm!(I, J, K,
                           promote_multiplier(α, A, B), transA, A,
                           transB, B,
                           promote_multiplier(β, C), C)
end

# Generic Julia implementation.

function _lgemm(::Generic,
                α::Number,
                transA::Char,
                A::AbstractMatrix{<:Floats},
                transB::Char,
                B::AbstractMatrix{<:Floats},
                Nc::Int)
    I, J, K = _lgemm_indices(transA, A, transB, B, Nc)
    T = _lgemm_type(α, A, B)
    return _generic_lgemm!(I, J, K,
                           promote_multiplier(α, A, B), transA, A,
                           transB, B,
                           promote_multiplier(0, T),
                           similar(Array{T}, (I, J)))
end

function _lgemm!(::Generic,
                 α::Number,
                 transA::Char,
                 A::AbstractMatrix{<:Floats},
                 transB::Char,
                 B::AbstractMatrix{<:Floats},
                 β::Number,
                 C::AbstractMatrix{<:Floats})
    I, J, K = _lgemm_indices(transA, A, transB, B, C)
    return _generic_lgemm!(I, J, K,
                           promote_multiplier(α, A, B), transA, A,
                           transB, B,
                           promote_multiplier(β, C), C)
end

function _lgemm(::Generic,
                α::Number,
                transA::Char,
                A::AbstractArray{<:Floats},
                transB::Char,
                B::AbstractArray{<:Floats},
                Nc::Int)
    I, J, K = _lgemm_indices(transA, A, transB, B, Nc)
    T = _lgemm_type(α, A, B)
    return _generic_lgemm!(cartesian_indices(I),
                           cartesian_indices(J),
                           cartesian_indices(K),
                           promote_multiplier(α, A, B), transA, A,
                           transB, B,
                           promote_multiplier(0, T),
                           similar(Array{T}, (I..., J...)))
end

function _lgemm!(::Generic,
                 α::Number,
                 transA::Char,
                 A::AbstractArray{<:Floats},
                 transB::Char,
                 B::AbstractArray{<:Floats},
                 β::Number,
                 C::AbstractArray{<:Floats})
    I, J, K = _lgemm_indices(transA, A, transB, B, C)
    return _generic_lgemm!(cartesian_indices(I),
                           cartesian_indices(J),
                           cartesian_indices(K),
                           promote_multiplier(α, A, B), transA, A,
                           transB, B,
                           promote_multiplier(β, C), C)
end

#
# Call low-level BLAS version.  The differences with LinearAlgebra.BLAS.gemm!
# are that inputs are assumed to be flat arrays (see is_flat_array) and that
# multipliers are automatically converted.
#
for (f, T) in ((:dgemm_, Float64),
               (:sgemm_, Float32),
               (:zgemm_, ComplexF64),
               (:cgemm_, ComplexF32))
    @eval begin
        #
        # FORTRAN prototype:
        #
        #     SUBROUTINE ${pfx}GEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,
        #                           C,LDC)
        #         ${T} ALPHA,BETA
        #         INTEGER M,N,K,LDA,LDB,LDC
        #         CHARACTER TRANSA,TRANSB
        #         ${T} A(LDA,*),B(LDB,*),C(LDC,*)
        #
        # Scalar arguments, α and β, can just be `Number` and integer arguments
        # can just be `Integer` but we want to keep the signature strict
        # because it is a low-level private method.
        #
        function _blas_lgemm!(m::Int, n::Int, p::Int, α::($T),
                              transA::Char, A::AbstractArray{$T},
                              transB::Char, B::AbstractArray{$T},
                              β::($T), C::AbstractArray{$T})
            lda = (transA == 'N' ? m : p)
            ldb = (transB == 'N' ? p : n)
            ldc = m
            ccall((@blasfunc($f), libblas), Cvoid,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ref{BlasInt}, Ref{$T}, Ptr{$T}, Ref{BlasInt},
                   Ptr{$T}, Ref{BlasInt}, Ref{$T}, Ptr{$T},
                   Ref{BlasInt}),
                  transA, transB, m, n, p, α, A, lda, B, ldb, β, C, ldc)
            return C
        end
    end
end
#
# This method is based on reference BLAS level 2 routine ZGEMM of LAPACK
# (version 3.7.0).  It is assumed that arrays A, B and C are flat arrays which
# can be linearly indexed.  Arguments are assumed to be correct, no checking is
# done.
#
function _linear_lgemm!(m::Int,
                        n::Int,
                        p::Int,
                        α::Number,
                        transA::Char,
                        A::AbstractArray{Ta,Na},
                        transB::Char,
                        B::AbstractArray{Tb,Nb},
                        β::Number,
                        C::AbstractArray{Tc,Nc}) where {Ta<:Floats,Na,
                                                        Tb<:Floats,Nb,
                                                        Tc<:Floats,Nc}
    #
    # Quick return if possible.
    #
    if m > 0 && n > 0 && (β != 1 || (p > 0 && α != 0))
        T = promote_type(Ta, Tb)
        if α == 0
            #
            # Quick computations when α = 0.
            #
            if β == 0
                @inbounds @simd for k in eachindex(C)
                    C[k] = zero(Tc)
                end
            elseif β != 1
                @inbounds @simd for k in eachindex(C)
                    C[k] *= β
                end
            end
        elseif transB == 'N'
            if transA == 'N'
                #
                # Form  C := α*A*B + β*C.
                #
                @inbounds for j in 1:n
                    if β == 0
                        @simd for i in 1:m
                            C[m*(j - 1) + i] = zero(Tc)
                        end
                    elseif β != 1
                        @simd for i in 1:m
                            C[m*(j - 1) + i] *= β
                        end
                    end
                    for k in 1:p
                        temp = α*B[p*(j - 1) + k]
                        if temp != zero(temp)
                            @simd for i in 1:m
                                C[m*(j - 1) + i] += temp*A[m*(k - 1) + i]
                            end
                        end
                    end
                end
            elseif Ta <: Real || transA == 'T'
                #
                # Form  C := α*A^T*B + β*C
                #
                @inbounds for j in 1:n
                    for i in 1:m
                        temp = zero(T)
                        @simd for k in 1:p
                            temp += A[p*(i - 1) + k]*B[p*(j - 1) + k]
                        end
                        C[m*(j - 1) + i] = (β == 0 ? α*temp : α*temp + β*C[m*(j - 1) + i])
                    end
                end
            else
                #
                # Form  C := α*A^H*B + β*C.
                #
                @inbounds for j in 1:n
                    for i in 1:m
                        temp = zero(T)
                        @simd for k in 1:p
                            temp += conj(A[p*(i - 1) + k])*B[p*(j - 1) + k]
                        end
                        C[m*(j - 1) + i] = (β == 0 ? α*temp : α*temp + β*C[m*(j - 1) + i])
                    end
                end
            end
        elseif transA == 'N'
            if Tb <: Real || transB == 'T'
                #
                # Form  C := α*A*B^T + β*C
                #
                @inbounds for j in 1:n
                    if β == 0
                        @simd for i in 1:m
                            C[m*(j - 1) + i] = zero(Tc)
                        end
                    elseif β != one
                        @simd for i in 1:m
                            C[m*(j - 1) + i] *= β
                        end
                    end
                    for k in 1:p
                        temp = α*B[n*(k - 1) + j]
                        if temp != zero(temp)
                            @simd for i in 1:m
                                C[m*(j - 1) + i] += temp*A[m*(k - 1) + i]
                            end
                        end
                    end
                end
            else
                #
                # Form  C := α*A*B^H + β*C.
                #
                @inbounds for j in 1:n
                    if β == 0
                        @simd for i in 1:m
                            C[m*(j - 1) + i] = zero(Tc)
                        end
                    elseif β != one
                        @simd for i in 1:m
                            C[m*(j - 1) + i] *= β
                        end
                    end
                    for k in 1:p
                        temp = α*conj(B[n*(k - 1) + j])
                        if temp != zero(temp)
                            @simd for i in 1:m
                                C[m*(j - 1) + i] += temp*A[m*(k - 1) + i]
                            end
                        end
                    end
                end
            end
        elseif Ta <: Real || transA == 'T'
            if Tb <: Real || transB == 'T'
                #
                # Form  C := α*A^T*B^T + β*C
                #
                @inbounds for j in 1:n
                    for i in 1:m
                        temp = zero(T)
                        @simd for k in 1:p
                            temp += A[p*(i - 1) + k]*B[n*(k - 1) + j]
                        end
                        C[m*(j - 1) + i] = (β == 0 ? α*temp : α*temp + β*C[m*(j - 1) + i])
                    end
                end
            else
                #
                # Form  C := α*A^T*B^H + β*C
                #
                @inbounds for j in 1:n
                    for i in 1:m
                        temp = zero(T)
                        @simd for k in 1:p
                            temp += A[p*(i - 1) + k]*conj(B[n*(k - 1) + j])
                        end
                        C[m*(j - 1) + i] = (β == 0 ? α*temp : α*temp + β*C[m*(j - 1) + i])
                    end
                end
            end
        else
            if Tb <: Real || transB == 'T'
                #
                # Form  C := α*A^H*B^T + β*C
                #
                @inbounds for j in 1:n
                    for i in 1:m
                        temp = zero(T)
                        @simd for k in 1:p
                            temp += conj(A[p*(i - 1) + k])*B[n*(k - 1) + j]
                        end
                        C[m*(j - 1) + i] = (β == 0 ? α*temp : α*temp + β*C[m*(j - 1) + i])
                    end
                end
            else
                #
                # Form  C := α*A^H*B^H + β*C.
                #
                @inbounds for j in 1:n
                    for i in 1:m
                        temp = zero(T)
                        @simd for k in 1:p
                            temp += conj(A[p*(i - 1) + k])*conj(B[n*(k - 1) + j])
                        end
                        C[m*(j - 1) + i] = (β == 0 ? α*temp : α*temp + β*C[m*(j - 1) + i])
                    end
                end
            end
        end
    end
    return C
end
#
# This method is based on reference BLAS level 2 routine ZGEMM of LAPACK
# (version 3.7.0).  Arguments A, B and C are any kind of arrays indexed by the
# Cartesian indices in I, J, and K.  Arguments are assumed to be correct, no
# checking is done.
#
function _generic_lgemm!(I, J, K,
                         α::Number,
                         transA::Char,
                         A::AbstractArray{Ta,Na},
                         transB::Char,
                         B::AbstractArray{Tb,Nb},
                         β::Number,
                         C::AbstractArray{Tc,Nc}) where {Ta<:Floats,Na,
                                                         Tb<:Floats,Nb,
                                                         Tc<:Floats,Nc}
    #
    # Quick return if possible.
    #
    if length(I) > 0 && length(J) > 0 && (β != 1 || (length(K) > 0 && α != 0))
        T = promote_type(Ta, Tb)
        if α == 0
            #
            # Quick computations when  α = 0.
            #
            if β == 0
                @inbounds @simd for k in eachindex(C)
                    C[k] = zero(Tc)
                end
            elseif β != 1
                @inbounds @simd for k in eachindex(C)
                    C[k] *= β
                end
            end
        elseif transB == 'N'
            if transA == 'N'
                #
                # Form  C := α*A*B + β*C.
                #
                @inbounds for j in J
                    if β == 0
                        @simd for i in I
                            C[i,j] = zero(Tc)
                        end
                    elseif β != 1
                        @simd for i in I
                            C[i,j] *= β
                        end
                    end
                    for k in K
                        temp = α*B[k,j]
                        if temp != zero(temp)
                            @simd for i in I
                                C[i,j] += temp*A[i,k]
                            end
                        end
                    end
                end
            elseif Ta <: Real || transA == 'T'
                #
                # Form  C := α*A^T*B + β*C
                #
                @inbounds for j in J
                    for i in I
                        temp = zero(T)
                        @simd for k in K
                            temp += A[k,i]*B[k,j]
                        end
                        C[i,j] = (β == 0 ? α*temp : α*temp + β*C[i,j])
                    end
                end
            else
                #
                # Form  C := α*A^H*B + β*C.
                #
                @inbounds for j in J
                    for i in I
                        temp = zero(T)
                        @simd for k in K
                            temp += conj(A[k,i])*B[k,j]
                        end
                        C[i,j] = (β == 0 ? α*temp : α*temp + β*C[i,j])
                    end
                end
            end
        elseif transA == 'N'
            if Tb <: Real || transB == 'T'
                #
                # Form  C := α*A*B^T + β*C
                #
                @inbounds for j in J
                    if β == 0
                        @simd for i in I
                            C[i,j] = zero(Tc)
                        end
                    elseif β != one
                        @simd for i in I
                            C[i,j] *= β
                        end
                    end
                    for k in K
                        temp = α*B[j,k]
                        if temp != zero(temp)
                            @simd for i in I
                                C[i,j] += temp*A[i,k]
                            end
                        end
                    end
                end
            else
                #
                # Form  C := α*A*B^H + β*C.
                #
                @inbounds for j in J
                    if β == 0
                        @simd for i in I
                            C[i,j] = zero(Tc)
                        end
                    elseif β != one
                        @simd for i in I
                            C[i,j] *= β
                        end
                    end
                    for k in K
                        temp = α*conj(B[j,k])
                        if temp != zero(temp)
                            @simd for i in I
                                C[i,j] += temp*A[i,k]
                            end
                        end
                    end
                end
            end
        elseif Ta <: Real || transA == 'T'
            if Tb <: Real || transB == 'T'
                #
                # Form  C := α*A^T*B^T + β*C
                #
                @inbounds for j in J
                    for i in I
                        temp = zero(T)
                        @simd for k in K
                            temp += A[k,i]*B[j,k]
                        end
                        C[i,j] = (β == 0 ? α*temp : α*temp + β*C[i,j])
                    end
                end
            else
                #
                # Form  C := α*A^T*B^H + β*C
                #
                @inbounds for j in J
                    for i in I
                        temp = zero(T)
                        @simd for k in K
                            temp += A[k,i]*conj(B[j,k])
                        end
                        C[i,j] = (β == 0 ? α*temp : α*temp + β*C[i,j])
                    end
                end
            end
        else
            if Tb <: Real || transB == 'T'
                #
                # Form  C := α*A^H*B^T + β*C
                #
                @inbounds for j in J
                    for i in I
                        temp = zero(T)
                        @simd for k in K
                            temp += conj(A[k,i])*B[j,k]
                        end
                        C[i,j] = (β == 0 ? α*temp : α*temp + β*C[i,j])
                    end
                end
            else
                #
                # Form  C := α*A^H*B^H + β*C.
                #
                @inbounds for j in J
                    for i in I
                        temp = zero(T)
                        @simd for k in K
                            temp += conj(A[k,i])*conj(B[j,k])
                        end
                        C[i,j] = (β == 0 ? α*temp : α*temp + β*C[i,j])
                    end
                end
            end
        end
    end
    return C
end
#
# This method yields the type of the elements for the result of lgemm.
#
_lgemm_type(α::Real, A::AbstractArray, B::AbstractArray) =
    promote_eltype(A, B)
_lgemm_type(α::Complex, A::AbstractArray, B::AbstractArray) =
    complex(promote_eltype(A, B))
#
# This method yields the lengths of the meta-dimensions for lgemm assuming
# linear indexing and check arguments.
#
@inline function _lgemm_dims(transA::Char,
                             A::AbstractMatrix,
                             transB::Char,
                             B::AbstractMatrix,
                             Nc::Int)
    Nc == 2 || incompatible_dimensions()
    lda = size(A, 1)
    ldb = size(B, 1)
    if transA == 'N'
        m, p = lda, size(A, 2)
    elseif transA == 'T' || transA == 'C'
        p, m = lda, size(A, 2)
    else
        invalid_transpose_character()
    end
    if transB == 'N'
        ldb == p || incompatible_dimensions()
        n = size(B, 2)
    elseif transB == 'T' || transB == 'C'
        size(B, 2) == p || incompatible_dimensions()
        n = ldb
    else
        invalid_transpose_character()
    end
    return m, n, p, (m, n)
end
#
# Idem for general arrays.
#
function _lgemm_dims(transA::Char,
                     A::AbstractArray{<:Any,Na},
                     transB::Char,
                     B::AbstractArray{<:Any,Nb},
                     Nc::Int) where {Na,Nb}
    nota = (transA == 'N')
    notb = (transB == 'N')
    if ((!nota && transA != 'T' && transA != 'C') ||
        (!notb && transB != 'T' && transB != 'C'))
        invalid_transpose_character()
    end
    Ni, Nj, Nk = _lgemm_ndims(Na, Nb, Nc)
    @inbounds begin
        p = 1
        if nota
            if notb
                # C[i,j] = sum_k A[i,k]*B[k,j]
                for d in 1:Nk
                    dim = size(A, d + Ni)
                    size(B, d) == dim || incompatible_dimensions()
                    p *= dim
                end
                shape = ntuple(d -> d ≤ Ni ? size(A, d) : size(B, d + (Nk - Ni)), Nc)
            else
                # C[i,j] = sum_k A[i,k]*B[j,k]
                for d in 1:Nk
                    dim = size(A, d + Ni)
                    size(B, d + Nj) == dim || incompatible_dimensions()
                    p *= dim
                end
                shape = ntuple(d -> d ≤ Ni ? size(A, d) : size(B, d - Ni), Nc)
            end
        else
            if notb
                # C[i,j] = sum_k A[k,i]*B[k,j]
                for d in 1:Nk
                    dim = size(A, d)
                    size(B, d) == dim || incompatible_dimensions()
                    p *= dim
                end
                shape = ntuple(d -> d ≤ Ni ? size(A, d + Nk) : size(B, d + (Nk - Ni)), Nc)
            else
                # C[i,j] = sum_k A[k,i]*B[j,k]
                for d in 1:Nk
                    dim = size(A, d)
                    size(B, d + Nj) == dim || incompatible_dimensions()
                    p *= dim
                end
                shape = ntuple(d -> d ≤ Ni ? size(A, d + Nk) : size(B, d - Ni), Nc)
            end
        end
        m = 1
        for d in 1:Ni
            m *= shape[d]
        end
        n = 1
        for d in Ni+1:Ni+Nj
            n *= shape[d]
        end
       return m, n, p, shape
    end
end
#
# This method yields the lengths of the meta-dimensions for lgemm! assuming
# linear indexing and check arguments.
#
@inline function _lgemm_dims(transA::Char,
                             A::AbstractMatrix,
                             transB::Char,
                             B::AbstractMatrix,
                             C::AbstractMatrix)
    m, n = size(C, 1), size(C, 2)
    adim1, adim2 = size(A, 1), size(A, 2)
    bdim1, bdim2 = size(B, 1), size(B, 2)
    if transA == 'N'
        adim1 == m || incompatible_dimensions()
        p = adim2
    elseif transA == 'T' || transA == 'C'
        adim2 == m || incompatible_dimensions()
        p = adim1
    else
        invalid_transpose_character()
    end
    if transB == 'N'
        (bdim1 == p && bdim2 == n) || incompatible_dimensions()
    elseif transB == 'T' || transB == 'C'
        (bdim1 == n && bdim2 == p) || incompatible_dimensions()
    else
        invalid_transpose_character()
    end
    return m, n, p
end
#
# Idem for general arrays.
#
function _lgemm_dims(transA::Char,
                     A::AbstractArray{<:Any,Na},
                     transB::Char,
                     B::AbstractArray{<:Any,Nb},
                     C::AbstractArray{<:Any,Nc}) where {Na,Nb,Nc}
    nota = (transA == 'N')
    notb = (transB == 'N')
    if ((!nota && transA != 'T' && transA != 'C') ||
        (!notb && transB != 'T' && transB != 'C'))
        invalid_transpose_character()
    end
    Ni, Nj, Nk = _lgemm_ndims(Na, Nb, Nc)
    @inbounds begin
        m = p = 1
        if nota
            for d in 1:Ni
                dim = size(C, d)
                size(A, d) == dim || incompatible_dimensions()
                m *= dim
            end
            if notb
                # C[i,j] = sum_k A[i,k]*B[k,j]
                for d in 1:Nk
                    dim = size(A, d + Ni)
                    size(B, d) == dim || incompatible_dimensions()
                    p *= dim
                end
            else
                # C[i,j] = sum_k A[i,k]*B[j,k]
                for d in 1:Nk
                    dim = size(A, d + Ni)
                    size(B, d + Nj) == dim || incompatible_dimensions()
                    p *= dim
                end
            end
        else
            for d in 1:Ni
                dim = size(C, d)
                size(A, d + Nk) == dim || incompatible_dimensions()
                m *= dim
            end
            if notb
                # C[i,j] = sum_k A[k,i]*B[k,j]
                for d in 1:Nk
                    dim = size(A, d)
                    size(B, d) == dim || incompatible_dimensions()
                    p *= dim
                end
            else
                # C[i,j] = sum_k A[k,i]*B[j,k]
                for d in 1:Nk
                    dim = size(A, d)
                    size(B, d + Nj) == dim || incompatible_dimensions()
                    p *= dim
                end
            end
        end
        n = 1
        if notb
            for d in 1:Nj
                dim = size(C, d + Ni)
                size(B, d + Nk) == dim || incompatible_dimensions()
                n *= dim
            end
        else
            for d in 1:Nj
                dim = size(C, d + Ni)
                size(B, d) == dim || incompatible_dimensions()
                n *= dim
            end
        end
        return m, n, p
    end
end
#
# This method yields the indices of the meta-dimensions for lgemm and check
# arguments.
#
@inline function _lgemm_indices(transA::Char,
                                A::AbstractMatrix,
                                transB::Char,
                                B::AbstractMatrix,
                                Nc::Int)
    Nc == 2 || incompatible_dimensions()
    if transA == 'N'
        I, K = axes(A, 1), axes(A, 2)
    elseif transA == 'T' || transA == 'C'
        K, I = axes(A, 1), axes(A, 2)
    else
        invalid_transpose_character()
    end
    if transB == 'N'
        axes(B, 1) == K || incompatible_dimensions()
        J = axes(B, 2)
    elseif transB == 'T' || transB == 'C'
        axes(B, 2) == K || incompatible_dimensions()
        J = axes(B, 1)
    else
        invalid_transpose_character()
    end
    return I, J, K
end
#
# Idem for general arrays.
#
function _lgemm_indices(transA::Char,
                        A::AbstractArray{<:Any,Na},
                        transB::Char,
                        B::AbstractArray{<:Any,Nb},
                        Nc::Int) where {Na,Nb}
    nota = (transA == 'N')
    notb = (transB == 'N')
    if ((!nota && transA != 'T' && transA != 'C') ||
        (!notb && transB != 'T' && transB != 'C'))
        invalid_transpose_character()
    end
    Ni, Nj, Nk = _lgemm_ndims(Na, Nb, Nc)
    @inbounds begin
        if nota
            I = ntuple(d -> axes(A, d), Ni)
            K = ntuple(d -> axes(A, d + Ni), Nk)
        else
            I = ntuple(d -> axes(A, d + Nk), Ni)
            K = ntuple(d -> axes(A, d), Nk)
        end
        if notb
            for d in 1:Nk
                axes(B, d) == K[d] || incompatible_dimensions()
            end
            J = ntuple(d -> axes(B, d + Nk), Nj)
        else
            for d in 1:Nk
                axes(B, d + Nj) == K[d] || incompatible_dimensions()
            end
            J = ntuple(d -> axes(B, d), Nj)
        end
        return I, J, K
    end
end
#
# This method yields the indices of the meta-dimensions for lgemm! and check
# arguments.
#
@inline function _lgemm_indices(transA::Char,
                                A::AbstractMatrix,
                                transB::Char,
                                B::AbstractMatrix,
                                C::AbstractMatrix)
    I, J = axes(C, 1), axes(C, 2)
    if transA == 'N'
        axes(A, 1) == I || incompatible_dimensions()
        K = axes(A, 2)
    elseif transA == 'T' || transA == 'C'
        axes(A, 2) == I || incompatible_dimensions()
        K = axes(A, 1)
    else
        invalid_transpose_character()
    end
    if transB == 'N'
        (axes(B, 1) == K && axes(B, 2) == J) || incompatible_dimensions()
    elseif transB == 'T' || transB == 'C'
        (axes(B, 1) == J && axes(B, 2) == K) || incompatible_dimensions()
    else
        invalid_transpose_character()
    end
    return I, J, K
end
#
# Idem for general arrays.
#
function _lgemm_indices(transA::Char,
                        A::AbstractArray{<:Any,Na},
                        transB::Char,
                        B::AbstractArray{<:Any,Nb},
                        C::AbstractArray{<:Any,Nc}) where {Na,Nb,Nc}
    nota = (transA == 'N')
    notb = (transB == 'N')
    if ((!nota && transA != 'T' && transA != 'C') ||
        (!notb && transB != 'T' && transB != 'C'))
        invalid_transpose_character()
    end
    Ni, Nj, Nk = _lgemm_ndims(Na, Nb, Nc)
    @inbounds begin
        I = ntuple(d -> axes(C, d), Ni)
        J = ntuple(d -> axes(C, d + Ni), Nj)
        if nota
            for d in 1:Ni
                axes(A, d) == I[d] || incompatible_dimensions()
            end
            K = ntuple(d -> axes(A, d + Ni), Nk)
        else
            for d in 1:Ni
                axes(A, d + Nk) == I[d] || incompatible_dimensions()
            end
            K = ntuple(d -> axes(A, d), Nk)
        end
        if notb
            for d in 1:Nk
                axes(B, d) == K[d] || incompatible_dimensions()
            end
            for d in 1:Nj
                axes(B, d + Nk) == J[d] || incompatible_dimensions()
            end
        else
            for d in 1:Nk
                axes(B, d + Nj) == K[d] || incompatible_dimensions()
            end
            for d in 1:Nj
                axes(B, d) == J[d] || incompatible_dimensions()
            end
        end
        return I, J, K
    end
end


text:
import Escher: Elem, Tile, render
import Base: convert

type VegaLiteTile <: Tile
  json::AbstractString
  svg::Bool
  actions::Bool
end

function VegaLiteTile(vis::AbstractString)
  VegaLiteTile(vis, SVG, SAVE_BUTTONS)
end

convert(::Type{Tile}, v::VegaLiteVis) = VegaLiteTile(JSON.json(v.vis))

function render(plot::VegaLiteTile, state)
  attr = Dict(:json=>plot.json)
  plot.svg     && ( attr[:svg]     = "" )
  plot.actions && ( attr[:actions] = "" )

  Elem(:"vega-lite-plot", attributes = attr)
end

ext: jl


text:
"""
This module defines a type called `PFWI` that is to be constructed
before performing
and kind of inversion.
The following functionality is currently added in this module:
* a simple RTM
* full-waveform inversion
* source and receiver filter inversion
Once a `PFWI` object is constructed, the following routines update certain fields
of the `PFWI` after performing the inversion.
* `xfwi!` : updates the initial subsurface models in `PFWI` using its observed data
* `wfwi!` : updates the source and receiver filters in `PFWI` using its observed data
"""

include("X.jl")

# add data inverse covariance matrix here later
struct LS end # just cls inversion

struct LS_prior # cls inversion including prior
	# add data inverse covariance matrix here later, instead of just Float64
	pdgls::Float64
	# pdgls
	pmgls::Misfits.P_gls{Float64}
end
LS_prior(α1::Float64, Q::LinearMap)=LS_prior(α1, Misfits.P_gls(Q))
struct Migr end
struct Migr_FD end

# create a timer object, used for inversion, see TimerOutputs.jl
global const fwi_to = TimerOutput();

mutable struct PFWI{Tmod, Tdatamisfit, Tinv}
	attrib_inv::Tinv
	"model inversion variable"
	mx::X{Float64,1}
	mxm::X{Float64,1}
	"forward modelling parameters for"
	paf::PFdtd{Tmod}
	"base source wavelet"
	srcwav::SrcWav
	"adjoint source functions"
	adjsrc::SrcWav
	"acquisition ageometry"
	ageom::AGeom
	"acquisition ageometry for adjoint propagation"
	adjageom::AGeom
	"Seismic model on modeling grid"
	mediumm::Medium
	"background Seismic model"
	mediumm0::Medium
	"Seismic on inversion grid"
	mediumi::Medium
	"initialize on this model"
	mod_initial::Medium
	"prior model on the inversion grid"
	priori::Medium
	"gradient Seismic model on modeling grid"
	gmediumm::Medium
	"gradient Seismic on inversion grid"
	gmediumi::Medium
	"parameterization"
	parameterization::Vector{Symbol}
	"model preconditioning"
	mprecon::AbstractArray{Float64,2}
	"compute data misfit"
	paTD::Tdatamisfit
	paminterp::Interpolation.Kernel_2D{Float64}
	verbose::Bool
	"synthetic or field inversion"
	attrib::Symbol
end

Base.show(io::Base.IO, pa::PFWI)=nothing

"""
```julia
pa=SeisInvExpt(attrib_mod, attrib_inv; srcwav, ageom, tgrid, mediumm, mediumm_obs, dobs)
```

# Arguments

* `attrib_mod::Union{Fdtd,FdtdBorn}` : choose modeling attribute
* `attrib_inv::Union{LS,LS_prior,Migr,Migr_FD}` : choose inversion attribute
* `obj::Union{LS,LS_prior}` : which objective function?
  * `=LS()` : least-squares inversion
  * `=LS_prior([1.0, 0.5])` regularized least-squares inversion
  * `=Migr()` : migration image
  * `=Migr_FD()` : finite-difference gradient (used for testing only)

# Keyword Arguments
* `srcwav::SrcWav` : source wavelets
* `ageom::AGeom` : acquisition
* `tgrid::StepRangeLen` : modelling time grid
* `mediumm::Medium` : an instance of `Medium` (used to generate initial model on the inversion grid)
* `dobs::Data` : observed data (required only when `attrib==:field`)
* `mediumm_obs::Medium` : medium used for generating *observed data* (required only when `attrib==:synthetic`)

# Optional Keyword Arguments

* `igrid` : inversion grid (defaults to modeling grid)
* `igrid_interp_scheme` : interpolation scheme
  * `=:B1` linear
  * `=:B2` second order
* `mediumm0` : background seismic model for Born modelling and inversion
* `parameterization::Vector{Symbol}` : subsurface parameterization
  * `=[:χvp, :χrho, :null]` for vp and density contrast inversion
  * `=[:null, :χrho, :null]` for density-only contrast inversion
  * `=[:χvp, :null, :null]` for vp-only contrast  inversion
* `tgrid_obs::StepRangeLen` : time grid for observed data
* `attrib::Symbol=:synthetic` : an attribute to control class
  * `=:synthetic` synthetic data inversion
  * `=:field` field data inversion
"""
SeisInvExpt(m1::Union{Fdtd, FdtdBorn},m2::Union{LS,LS_prior,Migr,Migr_FD};args1...)=PFWI(m1,m2;args1...) # change this when you have more inversion experiments

"""
this method constructs prior term with Q=α*I
* `ninv` : number of inversion variables, use `xfwi_ninv`
* `α` : scalar
"""
function LS_prior(ninv::Int, α=[1.0, 0.5])
	Q=LinearMap(
	     (y,x)->LinearAlgebra.mul!(y,x,α[2]),
	     (y,x)->LinearAlgebra.mul!(y,x,α[2]),
		      ninv, ninv, ismutating=true, issymmetric=true)
	return LS_prior(α[1],Misfits.P_gls(Q))
end

include("getprop.jl")
include("core.jl")
include("prior.jl")
include("xfwi.jl")
include("xwfwi.jl")
include("gallery.jl")

"""
Convert the data `Data` to `SrcWav` after time reversal.
"""
function update_adjsrc!(adjsrc, δdat::Data, adjageom)
	for i in 1:length(adjageom)
		for field in names(δdat[i].d)[1]
			wav=adjsrc[i].d[field]
			dat=δdat[i].d[field]
			for ir in 1:δdat[i][:nr]
				nt=length(adjsrc[i].grid)
				for it in 1:nt
					@inbounds wav[it,ir]=value_adjsrc(dat[nt-it+1,ir],eval(field)())
				end
			end
		end
	end
	return nothing
end

value_adjsrc(s, ::P) = s
value_adjsrc(s, ::Vx) = -1.0*s # see adj tests
value_adjsrc(s, ::Vz) = -1.0*s

function generate_adjsrc(fields, tgrid, adjageom)
	adjsrc=SrcWav(tgrid,adjageom,fields)
	return adjsrc
end

"""
Constructor for `PFWI`

# Arguments

* `attrib_mod::Union{Fdtd, FdtdBorn}` : modeling attribute
* `attrib_inv::Union{LS,LS_prior,Migr,Migr_FD} : inversion attribute

# Keyword Arguments
* `srcwav::SrcWav` : source time functions
* `ageom::AGeom` : acquisition ageometry
* `tgrid::StepRangeLen` : modeling time grid
* `mediumm::Medium` : seismic model on modeling mesh
* `mediumm_obs::Medium=mediumm` : actual seismic model to generate *observed* data

# Optional Keyword Arguments
* `tgrid_obs::StepRangeLen` : time grid for observed data
* `igrid::Vector{StepRangeLen}=mediumm.mgrid` : inversion grid if different from the modelling grid, i.e., `mediumm.mgrid`
* `igrid_interp_scheme` : interpolation scheme
  * `=:B1` linear
  * `=:B2` second order
* `mprecon_factor::Float64=1` : factor to control model preconditioner, always greater than 1
  * `=1` means the preconditioning is switched off
  * `>1` means the preconditioning is switched on, larger the mprecon_factor, stronger the applied preconditioner becomes
* `dobs::Data` : observed data
* `dprecon::Data=ones(Data)` : data preconditioning, defaults to one
* `tlagssf_fracs=0.0` : maximum lagtime of unknown source filter
* `tlagrf_fracs=0.0` : maximum lagtime of unknown receiver filter
* `srcwav_obs::SrcWav=srcwav` : source wavelets to generate *observed data*; can be different from `srcwav`
* `mediumm0::Medium=fill!(mediumm,0.0)` : background seismic model for Born modelling and inversion
* `parameterization::Vector{Symbol}` : subsurface parameterization
  * `=[:χvp, :χrho, :null]` for velocity and density inversion
  * `=[:null, :χrho, :null]` for density-only inversion
  * `=[:χvp, :null, :null]` for velocity-only inversion
* `verbose::Bool=false` : print level on STDOUT during inversion
* nworker : number of workers (input nothing to use all available)
"""
function PFWI(
	      attrib_mod::Union{Fdtd, FdtdBorn},
	      attrib_inv::Union{LS,LS_prior,Migr,Migr_FD};
	      srcwav::SrcWav=nothing,
	      ageom::AGeom=nothing,
	      tgrid::StepRangeLen=nothing,
	      mediumm::Medium=nothing,
	      # other optional
	      tgrid_obs::StepRangeLen=deepcopy(tgrid),
	      rfields=[:P],
	      igrid=nothing,
	      igrid_interp_scheme::Symbol=:B2,
	      mprecon_factor::Float64=1.0,
	      dobs::Data=Data(tgrid_obs, ageom, rfields),
	      dprecon=nothing,
	      tlagssf_fracs=0.0,
	      tlagrf_fracs=0.0,
	      srcwav_obs::SrcWav=deepcopy(srcwav),
	      mediumm_obs::Medium=mediumm,
	      mediumm0=nothing,
	      parameterization::Vector{Symbol}=[:χvp, :χrho, :null],
	      verbose::Bool=false,
	      attrib::Symbol=:synthetic,
	      nworker=nothing,
	      )

	# make a copy of initial model
	mediumm=deepcopy(mediumm)

	if(igrid===nothing)
		igrid=deepcopy(mediumm.mgrid)
	end
	# igrid has to truncated because the gradient evaluation
	# is inaccurate on boundaries
	mg=mediumm.mgrid
	amin=max(igrid[1][1],mg[1][3])
	amax=min(igrid[1][end],mg[1][end-2])
	if(length(igrid[1])==1)
		grid1=range(amin,step=min(step(igrid[1]), mg[1][end-2]-amin),length=1)
	else
		grid1=range(amin,stop=amax,length=length(igrid[1]))
	end
	amin=max(mg[2][3],igrid[2][1])
	amax=min(igrid[2][end],mg[2][end-2])
	if(length(igrid[2])==1)
		grid2=range(amin,step=min(step(igrid[2]), mg[2][end-2]-amin),length=1)
	else
		grid2=range(amin,stop=amax,length=length(igrid[2]))
	end
	igrid=[grid1, grid2]

	# create mediumi according to igrid and interpolation of mediumm
	mediumi = Medium(igrid, mediumm);
	# update bounds
	update!(mediumi,mediumm)
	interp_spray!(mediumm, mediumi, :interp)

	# save the initial model, as mediumi will be changed during the iterations
	mod_initial=deepcopy(mediumi)

	# allocate prior inputs
	priori=similar(mediumi)

	# use default background model mediumm0
	if(mediumm0 === nothing)
		mediumm0=deepcopy(mediumm) # some dummy
		fill!(mediumm0)
	end
	isequal(mediumm0, mediumm_obs) && (@warn "mediumm0 == mediumm_obs")

	if(attrib == :synthetic)
		isequal(mediumm, mediumm_obs) && error("initial model same as actual model, zero misfit?")
	end

	# ageom ageometry for adjoint propagation
	adjageom = AdjAGeom(ageom)

	# generate adjoint sources
	adjsrc=generate_adjsrc(rfields, tgrid, adjageom)

	# generating forward and adjoint modelling engines
	# to generate modelled data, border values, etc.
	# most of the parameters given to this are dummy
	paf=SeisForwExpt(attrib_mod,npw=2, medium=mediumm,
		ageom=[ageom, adjageom], srcwav=[srcwav, adjsrc],
		rfields=rfields,
		sflags=[3, 2], rflags=[1, 1],
		backprop_flag=1,
		tgrid=tgrid, gmodel_flag=true, verbose=verbose, illum_flag=false, nworker=nworker)

	gmediumm=Medium(mediumm.mgrid,[:χvp,:χrho])
	gmediumi=Medium(mediumi.mgrid,[:χvp,:χrho])

	#println("added fake precon")
	#dprecon=deepcopy(dobs)
	#fill!(dprecon, 1.0)
	#Data.taper!(dprecon, 20.)

	# check dprecon
	if(!(dprecon===nothing))
		(!(isapprox(dprecon,dobs))) && error("invalid dprecon used")
	end

	# create Parameters for data misfit
	#coup=Coupling.TD_delta(dobs.tgrid, tlagssf_fracs, tlagrf_fracs, rfields, ageom)
	# choose the data misfit
	#	(iszero(tlagssf_fracs))
 	# tlagssf_fracs==[0.0] | tlagssf_fracs=[0.0])
	# paTD=Data.P_misfit(Data.TD_zeros(rfields,tgrid,ageom),dobs,w=dprecon,coup=coup, func_attrib=optims[1]);

	if((attrib == :synthetic))
		if(iszero(dobs))
			Random.randn!(dobs) # dummy dobs, update later
		end
	else
		if(iszero(dobs))
			error("input observed data is zero")
		end
	end

	paTD=VNamedD_misfit(Data(tgrid,ageom,rfields),dobs,w=dprecon);

	paminterp=Interpolation.Kernel([mediumi.mgrid[2], mediumi.mgrid[1]], [mediumm.mgrid[2], mediumm.mgrid[1]], igrid_interp_scheme)

	mx=X(prod(length.(mediumi.mgrid))*count(parameterization.≠:null),2)
	mxm=X(prod(length.(mediumm.mgrid))*count(parameterization.≠:null),2)

	# put mediumm as a vector, according to parameterization in mxm.x
	copyto!(mxm.x,mediumm,parameterization)

	pa = PFWI(
		  attrib_inv,
	     mx, mxm,
	     paf,
	     deepcopy(srcwav),
	     SrcWav(tgrid,adjageom, rfields),
	     deepcopy(ageom),
	     adjageom,
	     deepcopy(mediumm), deepcopy(mediumm0), mediumi, mod_initial,
	     priori,
	     gmediumm,gmediumi,
	     parameterization,
	     zeros(2,2), # dummy, update mprecon later
	     paTD,
	     paminterp,
	     verbose, attrib)

	# generate observed data if attrib is synthetic
	if((attrib == :synthetic))
		update_observed_data!(pa, mediumm_obs)
	end

	iszero(dobs) && ((attrib == :real) ? error("input observed data for real data inversion") : error("problem generating synthetic observed data"))

	# generate modelled data
	F!(pa, nothing)

#	build_mprecon!(pa, Array(pa.paf.c.illum_stack), mprecon_factor)
	pa.paf.c.illum_flag=false # switch off illum flag for speed

	# default weights are 1.0
	fill!(pa.mx.w, 1.0)

	# update priori in pa to defaults
	update_prior!(pa)

	return pa
end

"""
update the *synthetic* observed data in `PFWI`
* allocated memory, don't use in inner loops
"""
function update_observed_data!(pa::PFWI, mediumm_obs)
	# save mediumm of pa to put it back later
	mediumm_copy=deepcopy(pa.mediumm)

	# change mediumm in pa to actual model
	copyto!(pa.mediumm, mediumm_obs)

	# update models in the forward engine and do modelling
	# FdtdBorn: mediumm0 will be used as a background model for Born modelling and mediumm will be perturbed model
	# Fdtd: mediumm will be *the* model
	F!(pa, nothing)

	# get observed data
	copyto!(pa.paTD.y, pa.paTD.x)

	# put back model and sources
	copyto!(pa.mediumm, mediumm_copy)

	initialize!(pa.paf.c)  # clear everything
end

"""
Return the number of inversion variables for PFWI corresponding to `PFWI`.
This number of inversion variables depend on the size of inversion mesh.
"""
function xfwi_ninv(pa::PFWI)
	length(pa.mx.x)
end

"""
Return the number of inversion variables for source and receiver filter inversion
corresponding to `PFWI`.
This number depends on the maximum lagtimes of the filters.
"""
function wfwi_ninv(pa::PFWI)
	return  pa.paTD.coup.tgridssf.nx
end

"build a model precon"
function build_mprecon!(pa,illum::Array{Float64}, mprecon_factor=1.0)
	"illum cannot be zero when building mprecon"
	(any(illum .<= 0.0)) && error("illum cannot be negative or zero")
	(mprecon_factor < 1.0)  && error("invalid mprecon_factor")

	illumi = zeros(length(pa.mediumi.mgrid[1]), length(pa.mediumi.mgrid[2]))
	Interpolation.interp_spray!(illumi, illum, pa.paminterp, :spray)
	(any(illumi .<= 0.0)) && error("illumi cannot be negative or zero")

	# use log to fix scaling of illum
	# (note that illum is not the exact diagonal of the Hessian matrix, so I am doing tricks here)
	# larger the mprecon_factor, stronger the preconditioning
	# division by maximum is safe because illum is non zero
	minillumi = minimum(illumi)
	maxillumi = maximum(illumi)
	illumi -= minillumi; # remove mean
	(maxillumi ≠ 0.0) && (illumi ./= maxillumi) # normalize using maxillumi
	#illumi = 1.0 + log.(1. + (mprecon_factor-1) * illumi) # illumi now lies between 1.0 and some value depending on how large illumi is
	illumi = (1. + (mprecon_factor-1) * illumi) # illumi now lies between 1.0 and some value depending on how large illumi is

	illumi = vcat([vec(illumi) for i in 1:count(pa.parameterization.≠:null)]...)

	length(illumi) == xfwi_ninv(pa) ? nothing : error("something went wrong with creating mprecon")
	pa.mprecon = spdiagm(illumi,(0),xfwi_ninv(pa), xfwi_ninv(pa))
end

# convert x to mediumm by interpolation
# x is on the sparse grid, and is changed by parameterization
function x_to_mediumm!(pa, x)

	# put sparse x into mediumi (useless step, just for visualization)
	copyto!(pa.mediumi, x, pa.parameterization)

	# get x on dense grid
	Interpolation.interp_spray!(x, pa.mxm.x, pa.paminterp, :interp, count(pa.parameterization.≠:null))

	# put dense x into mediumm, according to parameterization
	copyto!(pa.mediumm,pa.mxm.x,pa.parameterization)
end

function δx_to_δmods!(pa, δx)
	# get x on dense grid
	Interpolation.interp_spray!(δx, pa.mxm.x, pa.paminterp, :interp, count(pa.parameterization.≠:null))

	# reparameterize accordingly to get [δKI, δrhoI]
	copyto!(pa.paf.c.δmod, pa.mxm.x, pa.mediumm0, pa.parameterization)

	# input [δKI, δrhoI] to the modeling engine
	update_δmods!(pa.paf.c, pa.paf.c.δmod)
end

"""
convert the gradient output from Fdtd to gx
"""
function spray_gradient!(gx,  pa::PFWI{Fdtd,T1,T2}) where {T1,T2}

	# apply chain rule on gmediumm to get gradient w.r.t. dense x
	chainrule!(pa.mxm.gx, pa.paf.c.gradient, pa.mediumm, pa.parameterization)

	# spray gradient w.r.t. dense x to the space of sparse x
	Interpolation.interp_spray!(gx, pa.mxm.gx, pa.paminterp, :spray, count(pa.parameterization.≠:null))

	# visualize gmediumi here?
	# visualize_gx!(pa.gmediumm, pa.mediumm, pa.gmediumi, pa.mediumi, gx, pa)
end
function spray_gradient!(gx, pa::PFWI{FdtdBorn,T1,T2}) where {T1,T2}

	# apply chain rule on gmediumm to get gradient w.r.t. dense x (note that background model is used for Born modeling here)
	chainrule!(pa.mxm.gx, pa.paf.c.gradient, pa.mediumm, pa.parameterization)

	# spray gradient w.r.t. dense x to the space of sparse x
	Interpolation.interp_spray!(gx, pa.mxm.gx, pa.paminterp, :spray, count(pa.parameterization.≠:null))

	# visualize gmediumi here?
	# visualize_gx!(pa.gmediumm, pa.mediumm, pa.gmediumi, pa.mediumi, gx, pa)
end

"""
Return bound vectors for the `Seismic` model,
depeding on paramaterization
"""
function Seismic_xbound!(lower_x, upper_x, pa)
	nznx = prod(length.(pa.mediumi.mgrid))

	bound1 = similar(lower_x)
	# create a Seismic model with minimum possible values
	modbound = deepcopy(pa.mediumi);
	copyto!(modbound[:vp], fill(modbound.bounds[:vp][1], size(modbound[:vp])))
	if (:vs ∈ names(modbound.m)[1])
	#	copyto!(modbound[:vs], fill(modbound.bounds[:vs][1], size(modbound[:vs])))
	end
	copyto!(modbound[:rho], fill(modbound.bounds[:rho][1], size(modbound[:rho])))

	copyto!(bound1,modbound,pa.parameterization)

	bound2=similar(upper_x)
	# create a Seismic model with maximum possible values
	copyto!(modbound[:vp], fill(modbound.bounds[:vp][2], size(modbound[:vp])))
	if (:vs ∈ names(modbound.m)[1])
	#	copyto!(modbound[:vs], fill(modbound.bounds[:vs][2], size(modbound[:vs])))
	end
	copyto!(modbound[:rho], fill(modbound.bounds[:rho][2], size(modbound[:rho])))

	copyto!(bound2,modbound,pa.parameterization)

	# this sorting operation is important because
	lower_x[:] = min.(bound1, bound2)
	upper_x[:] = max.(bound1, bound2)

	return modbound
end

# gmediumi, gmediumm <-- gx (just for visualizing gradient)
function visualize_gx!(gmediumm, mediumm, gmediumi, mediumi, gx, pa)
	# chain rule depending on re-parameterization
	chainrule!(gmediumi, mediumi, gx, pa.parameterization, 1)

	# get gradient after interpolation (just for visualization, not exact)
#	interp_spray!(gmediumi, gmediumm, :interp, :B2, [:χvp,:χrho], pa=pa.paminterp)
end

"""
Modify the input acquisition ageometry
such that the adjoint source time functions can
be propagated from the receiver positions.
The number of supersources will remain the same.
All the recievers will be fired as simultaneous sources.
"""
function AdjAGeom(ageomin::AGeom)
	ageomout = deepcopy(ageomin);
	for i in 1:length(ageomin)
		ageomout[i].sx = ageomin[i].rx; ageomout[i].sz = ageomin[i].rz;
		ageomout[i].ns = ageomin[i].nr;

	#	ageomout.rx = ageomin.sx; ageomout.rz = ageomin.sz;
	#	ageomout.nr = ageomin.ns;

		ageomout[i].rx = ageomin[i].rx; ageomout[i].rz = ageomin[i].rz;
		ageomout[i].nr = ageomin[i].nr;
	end

	return ageomout
end

"""
Return functional and gradient of the LS objective
"""
function func_grad_Coupling!(storage, x::Vector{Float64},pa::PFWI)

	pa.verbose ? println("computing gradient...") : nothing

	# x to w
	Coupling_x!(x, pa, -1)

	if(storage === nothing)
		# compute misfit
		f = Data.error!(pa.paTD)
		return f
	else
		f = Data.error!(pa.paTD, :dJssf)
		Coupling_gx!(storage, pa.paTD.dJssf)
		return storage
	end
end

"""
Convert coupling functions to x and vice versa
"""
function Coupling_x!(x::Vector{Float64},
		   pa::PFWI,
		   flag::Int64)
	if(flag ==1) # convert w to x
		iss=1 # take only first source
		for ifield=1:length(pa.paTD.y.fields)
			w=pa.paTD.coup.ssf[iss,ifield]
			for i in eachindex(x)
				x[i]=w[i] # take only first source
			end
		end
	elseif(flag == -1) # convert x to w
		for iss=1:pa.paTD.y.ageom.nss, ifield=1:length(pa.paTD.y.fields)
			w=pa.paTD.coup.ssf[iss,ifield]
			for i in eachindex(x)
				w[i]=x[i]
			end
		end
	else
		error("invalid flag")
	end
end

"""
Convert dJssf to gradient vector
"""
function Coupling_gx!(gx, dJssf)
	gx[:]=0.0
	for ifield=1:size(dJssf,2), iss=1:size(dJssf,1)
		dwav=dJssf[iss,ifield]
		for i in eachindex(gx)
			gx[i]+=dwav[i]
		end
	end
end


text:
using Random
using LinearAlgebra
using Base.Threads
using IterativeSolvers
using TensorOperations
import Base.isempty
import Base.eltype
import Base.copy

"""
TT constructor of a tensor
``C[\\mu_1,…,\\mu_d] = A_1[1][\\mu_1]*⋯*A_d[d][\\mu_d] \\in \\mathbb{K}^{n_1 \\times ... \\times n_d}``
The following properties are stored
	* ttv_vec: the TT cores A_k as a list of 3-order tensors ``(A_1,...,A_d)`` where ``A_k = A_k[\\mu_k,\\alpha_{k-1},\\alpha_k]``, ``1 \\leq \\alpha_{k-1} \\leq r_{k-1}``, ``1 \\leq \\alpha_{k} \\leq r_{k}``, ``1 \\leq \\mu_k \\leq n_k``
	* ttv_dims: the dimension of the tensor along each mode
	* ttv_rks: the TT ranks ``(r_0,...,r_d)`` where ``r_0=r_d=1``
	* ttv_ot: the orthogonality of the TT where
		* ttv_ot[i] = 1 iff ``A_i`` is right-orthogonal *i.e.* ``\\sum_{\\mu_i} A_i[\\mu_i]^T A_i[\\mu_i] = I_{r_i}``
		* ttv_ot[i] = -1 iff ``A_i`` is left-orthogonal *i.e.* ``\\sum_{\\mu_i} A_i[\\mu_i] A_i[\\mu_i]^T = I_{r_{i-1}}``
		* ttv_ot[i] = 0 if nothing is known
"""
struct ttvector{T<:Number}
	ttv_vec :: Array{Array{T,3},1}
	ttv_dims :: Array{Int64,1}
	ttv_rks :: Array{Int64,1}
	ttv_ot :: Array{Int64,1}
end

Base.eltype(::ttvector{T}) where T<:Number = T
"""
Vidal representation of TT vector
"""
struct tt_vidal{T<:Number}
	#Vidal representation of a higher-order tensor
	#C[μ_1,…,μ_L] = core[μ_1] * Diagonal(Σ_1) * … * Diagonal(Σ_{L-1}) * core[μ_L]
	#Cores are orthogonal and Σ are the higher-order singular values
	core :: Array{Array{T,3},1}
	Σ :: Array{Array{Float64,1},1}
	dims :: Array{Int64,1}
	rks :: Array{Int64,1}
end

"""
TT constructor of a matrix ``M[i_1,..,i_d;j_1,..,j_d] \\in \\mathbb{K}^{n_1 \\cdots n_d \\times n_1 \\cdots n_d}`` where
``MM[i_1,..,i_d;j_1,..,j_d] = A_1[i_1,j_1] ... A_L[i_d,j_d]``
The following properties are stored
	* tto_vec: the TT cores A_k as a list of 4-order tensors ``(A_1,...,A_d)``
	* ttv_dims: the dimension of the tensor along each mode
	* ttv_rks: the TT ranks ``(r_0,...,r_d)`` where ``r_0=r_d=1``
"""
struct ttoperator{T<:Number}
	tto_vec :: Array{Array{T,4},1}
	tto_dims :: Array{Int64,1}
	tto_rks :: Array{Int64,1}
	tto_ot :: Array{Int64,1}
end

Base.eltype(::ttoperator{T}) where T<:Number = T

"""
returns a zero ttvector with dimensions `dims` and ranks `rks`
"""
function zeros_tt(dims,rks;T=Float64,ot=zeros(Int,length(dims)))
	@assert length(dims)+1==length(rks) "Dimensions and ranks are not compatible"
	d = length(dims)
	tt_vec = Array{Array{T,3}}(undef,d)
	for i in 1:d
		tt_vec[i] = zeros(T,dims[i],rks[i],rks[i+1])
	end
	return ttvector{T}(tt_vec,dims,copy(rks),copy(ot))
end

#returns partial isometry Q ∈ R^{n x m}
function rand_orthogonal(n,m;T=Float64)
    N = max(n,m)
    q,r = qr(rand(T,N,N))
    return Matrix(q)[1:n,1:m]
end

"""
Returns a random TTvector with dimensions `dims` and ranks `rks`
"""
function rand_tt(dims,rks;T=Float64)
	@assert length(dims)+1==length(rks) "Dimensions and ranks are not compatible"
	d = length(dims)
	tt_vec = Array{Array{T,3}}(undef,d)
	for i in 1:d
		tt_vec[i] = randn(T,dims[i],rks[i],rks[i+1])
	end
	return ttvector{T}(tt_vec,dims,copy(rks),zeros(Int,d))
end

function Base.copy(x_tt::ttvector{T}) where T<:Number
	y_tt = zeros_tt(x_tt.ttv_dims,x_tt.ttv_rks;T=T,ot=x_tt.ttv_ot)
	@threads for i in 1:length(x_tt.ttv_dims)
		y_tt.ttv_vec[i] = copy(x_tt.ttv_vec[i])
	end
	return y_tt
end

"""
TT decomposition by the Hierarchical SVD algorithm
	* Oseledets, I. V. (2011). Tensor-train decomposition. *SIAM Journal on Scientific Computing*, 33(5), 2295-2317.
	* Schollwöck, U. (2011). The density-matrix renormalization group in the age of matrix product states. *Annals of physics*, 326(1), 96-192.
The *root* of the TT decomposition is at index *i.e.* ``A_i`` for ``i < index`` are left-orthogonal and ``A_i`` for ``i > index`` are right-orthogonal. Singular values lower than tol are discarded.
"""
function ttv_decomp(tensor::Array{T};index=1,tol=1e-12) where T<:Number
	# Decomposes a tensor into its tensor train with core matrices at i=index
	dims = collect(size(tensor)) #dims = [n_1,...,n_d]
	n_max = maximum(dims)
	d = length(dims)
	ttv_vec = Array{Array{T}}(undef,d)
	# ttv_ot[i]= -1 if i < index
	# ttv_ot[i] = 0 if i = index
	# ttv_ot[i] = 1 if i > index
	ttv_ot = -ones(Int64,d)
	ttv_ot[index] = 0
	if index < d
		ttv_ot[index+1:d] = ones(d-index)
	end
	rks = ones(Int, d+1)
	tensor_curr = tensor
	# Calculate ttv_vec[i] for i < index
	for i = 1 : (index - 1)
		# Reshape the currently left tensor
		tensor_curr = reshape(tensor_curr, Int(rks[i] * dims[i]), :)
		# Perform the singular value decomposition
		u, s, v = svd(tensor_curr)
		# Define the i-th rank
		rks[i+1] = length(s[s .>= tol])
		# Initialize ttv_vec[i]
		ttv_vec[i] = zeros(T,dims[i],rks[i],rks[i+1])
		# Fill in the ttv_vec[i]
		for x = 1 : dims[i]
			ttv_vec[i][x, :, :] = u[(rks[i]*(x-1) + 1):(rks[i]*x), :]
		end
		# Update the currently left tensor
		tensor_curr = Diagonal(s[1:rks[i+1]])*v'[1:rks[i+1],:]
	end

	# Calculate ttv_vec[i] for i > index
	if index < d
		for i = d : (-1) : (index + 1)
			# Reshape the currently left tensor
			tensor_curr = reshape(tensor_curr, :, dims[i] * rks[i+1])
			# Perform the singular value decomposition
			u, s, v = svd(tensor_curr)
			# Define the (i-1)-th rank
			rks[i]=length(s[s .>= tol])
			# Initialize ttv_vec[i]
			ttv_vec[i] = zeros(T, dims[i], rks[i], rks[i+1])
			# Fill in the ttv_vec[i]
			i_vec = zeros(Int,rks[i+1])
			for x = 1 : dims[i]
				i_vec = dims[i]*((1:rks[i+1])-ones(Int, rks[i+1])) + x*ones(Int,rks[i+1])
				ttv_vec[i][x, :, :] = v'[1:rks[i],i_vec] #(rks[i+1]*(x-1)+1):(rks[i+1]*x)
			end
			# Update the currently left tensor
			tensor_curr = u[:,1:rks[i]]*Diagonal(s[1:rks[i]])
		end
	end
	# Calculate ttv_vec[i] for i = index
	# Reshape the currently left tensor
	tensor_curr = reshape(tensor_curr, Int(dims[index]*rks[index]),:)
	# Initialize ttv_vec[i]
	ttv_vec[index] = zeros(T, dims[index], rks[index], rks[index+1])
	# Fill in the ttv_vec[i]
	for x = 1 : dims[index]
		ttv_vec[index][x, :, :] =
		tensor_curr[Int(rks[index]*(x-1) + 1):Int(rks[index]*x), 1:rks[index+1]]
	end

	# Define the return value as a ttvector
	return ttvector{T}(ttv_vec, dims, rks, ttv_ot)
end

"""
Returns the tensor corresponding to x_tt
"""
function ttv_to_tensor(x_tt :: ttvector{T}) where T<: Number
	d = length(x_tt.ttv_dims)
	r_max = maximum(x_tt.ttv_rks)
	# Initialize the to be returned tensor
	tensor = zeros(T, x_tt.ttv_dims...)
	# Fill in the tensor for every t=(x_1,...,x_d)
	for t in CartesianIndices(tensor)
		curr = ones(T,r_max)
		a = collect(Tuple(t))
		for i = d:-1:1
			curr[1:x_tt.ttv_rks[i]] = x_tt.ttv_vec[i][a[i],:,:]*curr[1:x_tt.ttv_rks[i+1]]
		end
		tensor[t] = curr[1]
	end
	return tensor
end

#returns the Vidal representation of a TT
function tt_to_vidal(x_tt::ttvector{T};tol=1e-14) where T<:Number
	d = length(x_tt.ttv_dims)
	core = Array{Array{T,3},1}(undef,d)
	Σ = Array{Array{Float64,1},1}(undef,d-1)
	y_tt = orthogonalize(x_tt)
	y_rks = copy(y_tt.ttv_rks)
	for j in 1:d-1
		A = zeros(y_tt.ttv_dims[j],y_rks[j],y_tt.ttv_dims[j+1],y_rks[j+2])
		@tensor A[a,b,c,d] = y_tt.ttv_vec[j][a,z,b]*y_tt.ttv_vec[j+1][z,d,c]
		u,s,v = svd(reshape(A,size(A,1)*size(A,2),:))
		Σ[j] = s[s.>tol]
		y_rks[j+1] = length(Σ[j])
		core[j] = permutedims(reshape(u[:,s.>tol],y_rks[j],y_tt.ttv_dims[j],:),[1 3 2])
		y_tt.ttv_vec[j+1] = reshape(v[:,s.>tol],:,y_rks[j+2],y_tt.ttv_dims[j+1])
	end
	core[d] = y_tt.ttv_vec[d]
	return tt_vidal(core,Σ,y_tt.ttv_dims,y_rks)
end

"""
Transforms a ttoperator into a ttvector
"""
function tto_to_ttv(A::ttoperator{T}) where T<:Number
	d = length(A.tto_dims)
	xtt_vec = Array{Array{T,3},1}(undef,d)
	A_rks = A.tto_rks
	for i in 1:d
		xtt_vec[i] = reshape(A.tto_vec[i],A.tto_dims[i]^2,A_rks[i],A_rks[i+1])
	end
	return ttvector{T}(xtt_vec,A.tto_dims.^2,A.tto_rks,A.tto_ot)
end

"""
Transforms a ttvector (coming from a ttoperator) into a ttoperator
"""
function ttv_to_tto(x::ttvector{T}) where T<:Number
	d = length(x.ttv_dims)
	@assert(isqrt.(x.ttv_dims).^2 == x.ttv_dims, DimensionMismatch)
	Att_vec = Array{Array{T,4},1}(undef,d)
	x_rks = x.ttv_rks
	A_dims = zeros(Int64,d)
	for i in 1:d
		A_dims[i] = isqrt(x.ttv_dims[i])
		Att_vec[i] = reshape(x.ttv_vec[i],A_dims[i],A_dims[i],x_rks[i],x_rks[i+1])
	end
	return ttoperator{T}(Att_vec,A_dims,x.ttv_rks,x.ttv_ot)
end

"""
Returns the TT decomposition of a matrix using the HSVD algorithm
"""
function tto_decomp(tensor::Array{T}; index=1) where T<:Number
	# Decomposes a tensor operator into its tensor train
	# with core matrices at i=index
	# The tensor is given as tensor[x_1,...,x_d,y_1,...,y_d]
	d = Int(ndims(tensor)/2)
	tto_dims = collect(size(tensor))[1:d]
	n_max = maximum(tto_dims)
	dims_sq = tto_dims.^2
	# The tensor is reorder  into tensor[x_1,y_1,...,x_d,y_d],
	# reshaped into tensor[(x_1,y_1),...,(x_d,y_d)]
	# and decomposed into its tensor train with core matrices at i= index
	index_sorted = reshape(Transpose(reshape(1:(2*d),:,2)),1,:)
	ttv = ttv_decomp(reshape(permutedims(tensor,index_sorted),(dims_sq[1:(end-1)]...), :); index=index)
	# Define the array of ranks [r_0=1,r_1,...,r_d]
	rks = ttv.ttv_rks
	r_max = maximum(rks)
	# Initialize tto_vec
	tto_vec = Array{Array{T}}(undef,d)
	# Fill in tto_vec
	for i = 1:d
		# Initialize tto_vec[i]
		tto_vec[i] = zeros(T, tto_dims[i], tto_dims[i], rks[i], rks[i+1])
		# Fill in tto_vec[i]
		tto_vec[i] = reshape(ttv.ttv_vec[i], tto_dims[i], tto_dims[i], :, rks[i+1])
	end
	return ttoperator{T}(tto_vec, tto_dims, rks, ttv.ttv_ot)
end

function tto_to_tensor(tto :: ttoperator{T}) where T<:Number
	d = length(tto.tto_dims)
	# Define the array of ranks [r_0=1,r_1,...,r_d]
	rks = tto.tto_rks
	r_max = maximum(rks)
	# The tensor has dimensions [n_1,...,n_d,n_1,...,n_d]
	dims = zeros(Int, 2*d)
	dims[1:d] = tto.tto_dims
	dims[(d+1):(2*d)] = tto.tto_dims
	tensor = zeros(T,dims...)
	# Fill in the tensor for every t=(x_1,...,x_d,y_1,...,y_d)
	for t in CartesianIndices(tensor)
		curr = ones(T,r_max)
		a = collect(Tuple(t))
		for i = d:-1:1
			curr[1:rks[i]] = tto.tto_vec[i][t[i], t[d + i], :, :]*curr[1:rks[i+1]]
		end
		tensor[t] = curr[1]
	end
	return tensor
end

#TTO representation of the identity matrix
function id_tto(d;n_dim=2,T=Float64)
	dims = n_dim*ones(Int64,d)
	A = Array{Array{T,4},1}(undef,d)
	for j in 1:d
		A[j] = zeros(2,2,1,1)
		A[j][:,:,1,1] = Matrix{T}(I,2,2)
	end
	return ttoperator{T}(A,dims,ones(Int64,d+1),zeros(d))
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
module RobotDynamics

using Rotations
using StaticArrays
using LinearAlgebra
using ForwardDiff
using FiniteDiff
using UnsafeArrays
using RecipesBase
using SparseDiffTools
using SparseArrays

using Rotations: skew
using StaticArrays: SUnitRange

export
    AbstractModel,
    DynamicsExpansion,
    dynamics,
    jacobian!,
    discrete_dynamics,
    discrete_jacobian!,
    linearize,
    linearize!,
    state_dim,
    control_dim,
    state_diff_size,
    rollout!

# rigid bodies
export
    LieGroupModel,
    RigidBody,
    RBState,
    orientation,
    linear_velocity,
    angular_velocity

# linear model
export
    LinearModel,
    linear_dynamics,
    LinearizedModel,
    linearize_and_discretize!,
    discretize,
    discretize!,
    update_trajectory!

# knotpoints
export
    AbstractKnotPoint,
    KnotPoint,
    StaticKnotPoint,
    Traj,
    state,
    control,
    states,
    controls,
    set_states!,
    set_controls!

# integration
export
    QuadratureRule,
    RK2,
    RK3,
    RK4,
    HermiteSimpson,
    PassThrough,
    Exponential

include("rbstate.jl")
include("jacobian.jl")
include("knotpoint.jl")
include("model.jl")
include("liestate.jl")
include("rigidbody.jl")
include("integration.jl")
include("trajectories.jl")
include("linearmodel.jl")
include("linearization.jl")
include("plot_recipes.jl")

end # module


text:
sparseaccess(A::SparseMatrixCSC) = A
sparseaccess(A::SubArray{<:Any,2,<:SparseMatrixCSC{<:Any,<:Any},<:Tuple{I,<:Any},false}) where
    I<:Union{Base.Slice{<:Base.OneTo},AbstractUnitRange} = A

sparseaccess(A::SubArray) = XSubArray(A)

# the actual nonzero entries belonging to the view
nzvalview(S::SparseMatrixCSCView) = view(S.parent.nzval, first(nzrange(S, first(axes(S, 2)))):last(nzrange(S, last(axes(S, 2)))))

#
nonzeros(S::SubArray{<:Any,2,<:SparseMatrixCSC}) = nonzeros(S.parent)
nonzeros(S::XSubArray{<:Any,<:Integer,2}) = nonzeros(S.sub)

nnz(S::SparseMatrixCSCView) = nzrange(S.parent, last(S.indices[2]))[end] - nzrange(S.parent, first(S.indices[2]))[1] + 1
nzrange(S::SparseMatrixCSCView, col::Integer) = nzrange(S.parent, S.indices[2][col])
function nzrange(S::SubArray{<:Any,2,<:SparseMatrixCSC,<:Tuple{I,J},false}, i::Integer) where {I<:AbstractUnitRange,J<:AbstractVector{<:Integer}}
    A = S.parent
    r = nzrange(A, S.indices[2][i])
    rvA = rowvals(A)
    m = axes(A, 1)
    r1 = first(r)
    r2 = last(r)
    ri = S.indices[1]
    i1 = first(ri)
    i2 = last(ri)
    if i1 > first(m)
        r1 = searchsortedfirst(rvA, i1, r1, r2, Forward)
    end
    if i2 < last(m)
        r2 = searchsortedlast(rvA, i2, r1, r2, Forward)
    end
    r1:r2
end

function nzrange(xS::XSubArray{<:Any,<:Integer,2,<:SparseMatrixCSC,<:Tuple{I,J},false}, i::Integer) where {I<:AbstractVector{<:Integer},J<:AbstractVector{<:Integer}}
    S = xS.sub
    A = S.parent
    r = nzrange(A, S.indices[2][i])
    rvA = rowvals(A)
    m = axes(A, 1)
    rangestrip(rvA, r, m, xS)
end

function rangestrip(rvA, r::AbstractVector{Int}, m::AbstractUnitRange, S::XSubArray{<:Any,<:Integer,2,<:Any,<:Tuple{I,<:Any}}) where I<:StepRange
    r1 = first(r)
    r2 = last(r)
    ri = S.sub.indices[1]
    i1 = first(ri)
    i2 = last(ri)
    minmax(i1, i2)
    if i1 > first(m)
        r1 = searchsortedfirst(rvA, i1, r1, r2, Forward)
    end
    if i2 < last(m)
        r2 = searchsortedlast(rvA, i2, r1, r2, Forward)
    end
    nzr = Vector{eltype(ri)}(undef, r2 - r1 + 1)
    rv = rowvals(S)
    k = 0
    for i in r1:r2
        j = rv[i]
        if j != 0
            nzr[k+=1] = i
        end
    end
    resize!(nzr, k)
    step(ri) < 0 ? reverse(nzr) : nzr
end

function rangestrip(rvA, r::AbstractVector{Int}, m::AbstractUnitRange, S::XSubArray{<:Any,<:Integer,2,<:Any,<:Tuple{I,<:Any}}) where I<:AbstractVector{<:Integer}
    r1 = first(r)
    r2 = last(r)
    ri = S.sub.indices[1]
    i1, i2 = extrema(ri)
    if i1 > first(m)
        r1 = searchsortedfirst(rvA, i1, r1, r2, Forward)
    end
    if i2 < last(m)
        r2 = searchsortedlast(rvA, i2, r1, r2, Forward)
    end
    vind = rowvals(S)[r1:r2]
    p = sortperm(vind)
    k = searchsortedlast(vind[p], 0)
    (r1:r2)[p][k+1:end]
end

rowvals(S::SparseMatrixCSCView) = S.parent.rowval
rowvals(xs::XSubArray) = xs.rowval

struct RowIndexVector{U,T,D} <: AbstractVector{Int}
    rvp::T
    ind::U
    dict::D
end

function rowvals(S::SubArray{<:Any,2,<:Any,<:Tuple{<:Base.Slice,<:AbstractVector},false})
    rowvals(S.parent)
end
function rowvals(S::SubArray{<:Any,2,<:Any,<:Tuple{<:UnitRange,<:AbstractVector},false})
    RowIndexVector(rowvals(S.parent), S.indices[1], nothing)
end
function rowvals(S::SubArray{<:Any,2,<:Any,<:Tuple{<:OrdinalRange,<:AbstractVector},false})
    RowIndexVector(rowvals(S.parent), S.indices[1], nothing)
end
function rowvals(S::SubArray{<:Any,2,<:Any,<:Tuple{<:AbstractVector,<:AbstractVector},false})
    ind = S.indices[1]
    RowIndexVector(rowvals(S.parent), ind, Dict(reverse.(enumerate(ind))))
end

import Base: getindex, size
size(rv::RowIndexVector) = size(rv.rvp)
size(s::XSubArray) = size(s.sub)
getindex(s::XSubArray, ind...) = getindex(s.sub, ind...)

function getindex(rv::RowIndexVector{<:UnitRange}, k::Integer)
    i = rv.rvp[k]
    ind = rv.ind
    a = first(ind)
    b = last(ind)
    ( a <= i <= b ) ? ( i - a + 1 ) : 0
end
function getindex(rv::RowIndexVector{<:StepRange}, k::Integer)
    i = rv.rvp[k]
    ind = rv.ind
    a = first(ind)
    b = last(ind)
    a <= i <= b || return 0
    s = step(ind)
    i -= a
    rem(i, s) == 0 ? 0 : div(i, s)
end
function getindex(rv::RowIndexVector{<:AbstractVector}, k::Integer)
    i = rv.rvp[k]
    ind = rv.ind
    a, b = extrema(ind)
    get(rv.dict, i, 0)
end

SparseArrays.nzrange(X::SparseVector, i) = 1:length(X.nzind)
SparseArrays.rowvals(X::SparseVector) = X.nzind


text:
"""
    crandnToep(t)

Generate `x ~ N_C(0,T)` with Toeplitz cov. `T` that has `t` as
its first row.
"""
function crandnToep(t)
    C = toeplitzHe(t)
    y = sqrtm(C)*crandn(length(t))
end

"""
    toeplitz(t::Vector)

Generate an `n` x `n` Toeplitz matrix `C` with generating vector `t` of length `2n-1`.
`C[1,1] = t[n]`
`C[1,n] = t[2n-1]`
`C[n,1] = t[1]`
"""
function toeplitz(t::Vector)
    n = Int(floor(length(t)/2)+1)

    C = [t[j-i+n] for i=1:n,j=1:n]
end

function applyToep(t,x)
    x2 = [x;zeros(x)]
    t2 = [conj(t);0;t[end:-1:2]]
    y = ifft(fft(t2,1).*fft(x2,1))
    y = y[1:length(x)]
end

"""
    toeplitzHe(t::Vector)

Generate an `n` x `n` Hermitian Toeplitz matrix `C` with generating vector `t` of length `n`.
The vector `t` defines the first row of `C` and `t[1]` must be real.
"""
function toeplitzHe(t::Vector)
    t_full = [ conj(t[end:-1:2]); t[1:end] ]
    C = toeplitz(t_full)
end

"""
    best_circulant_approximation(t::Vector)

If `t` is the length `n` generating vector of a Hermitian Toeplitz matrix `T`, then
`c` is the length `n` generating vector of a Hermitian circulant matrix `C`, which
is the best approximation of `T`.
"""
function best_circulant_approximation(t::Vector)
    n = length(t);
    lins = ((n:-2:-n+2)./n);
    t1 = t;
    t2 = cat(1,[0], conj(reverse(t[2:end])));
    c = 0.5*n*real(ifft(lins .*(t1-t2) + t1 + t2));
end

"""
    best_circulant_approximation(t::Matrix)

If each `t[:,m]` is the length `n` generating vector of a Hermitian Toeplitz matrix `T_m`, then
each `c[:,m]` is the length `n` generating vector of a Hermitian circulant matrix `C_m`, which
is the best approximation of `T_m`.
"""
function best_circulant_approximation(t::Matrix)
    n = size(t,1);
    m = size(t,2);
    lins = ((n:-2:-n+2)./n);
    t1 = t;
    t2 = cat(1,zeros(1,m), conj(flipdim(t[2:end,:],1)));
    c = 0.5*n*real(ifft(lins .*(t1-t2) + t1 + t2,1));
end

max_stars_repo_name: tum-msv/learning-mmse-est


text:
using ACEdocs
using Test

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

max_stars_repo_name: ACEsuit/ACEdocs

max_issues_repo_name: ACEsuit/ACEdocs

max_forks_repo_name: ACEsuit/ACEdocs


text:
eachindex(arg::MultiValue) = eachindex(1:prod(size(arg)))

CartesianIndices(arg::MultiValue) = CartesianIndices(size(arg))

LinearIndices(arg::MultiValue) = LinearIndices(size(arg))

getindex(arg::VectorValue, i::Integer) = arg.data[i]

function getindex(arg::TensorValue{D},i::Integer,j::Integer) where D
  index = _2d_tensor_linear_index(D,i,j)
  arg.data[index]
end

function getindex(arg::SymTensorValue{D},i::Integer,j::Integer) where D
  index = _2d_sym_tensor_linear_index(D,i,j)
  arg.data[index]
end

function getindex(arg::SymFourthOrderTensorValue{D},i::Integer,j::Integer,k::Integer,l::Integer) where D
  index = _4d_sym_tensor_linear_index(D,i,j,k,l)
  arg.data[index]
end

function getindex(arg::ThirdOrderTensorValue{D1,D2},i::Integer,j::Integer,k::Integer) where {D1,D2}
  index = _3d_tensor_linear_index(D1,D2,i,j,k)
  arg.data[index]
end

getindex(arg::VectorValue, ci::CartesianIndex{1})              = getindex(arg,ci[1])
getindex(arg::TensorValue,ci::CartesianIndex{2})               = getindex(arg,ci[1],ci[2])
getindex(arg::SymTensorValue,ci::CartesianIndex{2})            = getindex(arg,ci[1],ci[2])
getindex(arg::ThirdOrderTensorValue,ci::CartesianIndex{3})     = getindex(arg,ci[1],ci[2],ci[3])
getindex(arg::SymFourthOrderTensorValue,ci::CartesianIndex{4}) = getindex(arg,ci[1],ci[2],ci[3],ci[4])

getindex(arg::MultiValue, i::Integer) = getindex(arg, CartesianIndices(arg)[i])
getindex(arg::TensorValue, i::Integer) = arg.data[i]
getindex(arg::ThirdOrderTensorValue, i::Integer) = arg.data[i]

@inline iterate(arg::MultiValue)        = iterate(arg.data)
@inline iterate(arg::MultiValue, state) = iterate(arg.data, state)

data_index(::Type{<:VectorValue},i) = i
data_index(::Type{<:TensorValue{D}},i,j) where D = _2d_tensor_linear_index(D,i,j)
data_index(::Type{<:SymTensorValue{D}},i,j) where D = _2d_sym_tensor_linear_index(D,i,j)
data_index(::Type{<:ThirdOrderTensorValue{D1,D2}},i,j,k) where {D1,D2} = _3d_tensor_linear_index(D1,D2,i,j,k)
data_index(::Type{<:SymFourthOrderTensorValue{D}},i,j,k,l) where D = _4d_sym_tensor_linear_index(D,i,j,k,l)

_symmetric_index_gaps(i::Integer) = i*(i-1)÷2

_2d_tensor_linear_index(D,i,j) = ((j-1)*D)+i

_3d_tensor_linear_index(D1,D2,i,j,k) = (k-1)*D1*D2+(j-1)*D1+i

function _2d_sym_tensor_linear_index(D,i,j)
  _j=min(i,j)
  _i=max(i,j)
  index=_2d_tensor_linear_index(D,_i,_j)-_symmetric_index_gaps(_j)
  index
end

#function _4d_sym_tensor_linear_index(D,i,j,k,l)
#  _j=min(i,j)
#  _i=max(i,j)
#  _l=min(k,l)
#  _k=max(k,l)
#  block_length=_symmetric_index_gaps(D+1)
#  element_index=_2d_tensor_linear_index(D,_i,_j)-_symmetric_index_gaps(_j)
#  block_index=_2d_tensor_linear_index(D,_l,_k)-_symmetric_index_gaps(_l)
#  index=(block_index-1)*block_length+element_index
#  index
#end

function _4d_sym_tensor_linear_index(D,i,j,k,l)
  block_length = (D*(D+1))÷2
  block_index = _2d_sym_tensor_linear_index(D,i,j)
  element_index = _2d_sym_tensor_linear_index(D,k,l)
  index=(block_index-1)*block_length+element_index
  index
end


text:
module TestAbstractTypes

using Test

const niters = 1 #TODO 100

function generator end
export generator

################################################################################

function generator(::Type{Char})
    function gen(channel::Channel{Char})
        yield!(x) = put!(channel, x)
        while true
            case = rand(1:20)
            if case == 1
                yield!('a')
            elseif case == 2
                yield!("β")
            elseif case == 3
                yield!(' ')
            elseif case == 4
                yield!('\0')
            elseif case == 5
                yield!(Char(127))
            elseif case == 6
                yield!(Char(128))
            elseif case ≤ 10
                yield!(Char(rand(0:127)))
            else
                len = rand(0:20)
                yield!(rand(Char))
            end
        end
    end
    return Channel{Char}(gen)
end

function testAbstractChar(::Type{T}, ch::AbstractChannel{T}) where {T}
    @testset "testAbstractChar($T)" begin
        @test T <: AbstractChar

        next!() = take!(ch)::T

        for iter in 1:niters
            x = next!()
            y = next!()

            if x == y
                @test cmp(x, y) == 0
            elseif x < y
                @test cmp(x, y) == -1
            elseif x > y
                @test cmp(x, y) == 1
            else
                @test false
            end

            xp = codepoint(x)
            @test xp isa Integer
            xc = T(xp)
            @test xc isa T
            @test xc == x
        end
    end
end
export testAbstractChar

################################################################################

function generator(::Type{String})
    function gen(channel::Channel{String})
        yield!(x) = put!(channel, x)
        while true
            case = rand(1:20)
            if case == 1
                yield!("")
            elseif case == 2
                yield!("a")
            elseif case == 3
                yield!("abc")
            elseif case == 4
                yield!("β")
            elseif case == 5
                yield!(" ")
            elseif case == 6
                yield!("\0")
            elseif case == 7
                yield!("\0abc")
            elseif case ≤ 10
                yield!(string(rand(Char)))
            elseif case ≤ 19
                len = rand(0:20)
                yield!(String(rand(Char, len)))
            else
                len = rand(10000:20000)
                yield!(String(rand(Char, len)))
            end
        end
    end
    return Channel{String}(gen)
end

function generator(::Type{SubString})
    function gen(channel::Channel{SubString})
        source = generator(String)
        while true
            str = take!(source)::String
            len = ncodeunits(str)
            if len == 0
                put!(channel, SubString(str))
            else
                i = thisind(str, rand(1:len))
                j = thisind(str, rand(1:len))
                if i ≤ j
                    put!(channel, SubString(str, i, j))
                end
            end
        end
    end
    return Channel{SubString}(gen)
end

function testAbstractString(::Type{T}, ch::AbstractChannel{T}) where {T}
    @testset "testAbstractString($T)" begin
        @test T <: AbstractString

        next!() = take!(ch)::T

        unit = T("")

        @test unit isa T
        @test length(unit) == 0

        for iter in 1:niters
            x = next!()
            y = next!()
            z = next!()

            # Equality
            @test unit == unit
            @test x == x

            # Unit
            @test unit * unit == unit
            @test unit * x == x
            @test x * unit == x

            # Associativity
            @test (x * y) * z == x * (y * z)

            # Exponentiation
            @test x^0 == unit
            @test x^1 == x
            @test x^2 == x * x
            n = rand(3:10)
            @test x^n == *((x for i in 1:n)...)

            # Lengths
            @test length(x) isa Int
            @test length(x * y) == length(x) + length(y)
            @test isempty(x) == (length(x) == 0)
            @test ncodeunits(x) isa Int
            @test ncodeunits(x * y) == ncodeunits(x) + ncodeunits(y)
            @test isempty(x) == (ncodeunits(x) == 0)

            # Elements
            # Find all possible indices
            inds = Int[]
            ind = firstindex(x)
            while ind ≤ lastindex(x)
                push!(inds, ind)
                ind = nextind(x, ind)
            end
            @test length(inds) == length(x)
            if !isempty(x)
                pos = inds[rand(1:length(inds))]
                @test eltype(x) isa Type
                @test x[pos] isa eltype(x)
                @test (x * y)[pos] == x[pos]
                @test (y * x)[ncodeunits(y) + pos] == x[pos]
            end

            # codeunit, codeunits
            # firstindex, lastindex, thisind, nextind, prevind, isvalid
            # cmp
        end

        close(ch)
    end
end
export testAbstractString

end


text:
"""
# module Shells

# Examples

```jldoctest
julia>
```
"""
module Shells

export Shell

struct Shell{T <: AbstractString}
    path::T
end

end


text:
function Emma()
        println("Willkommen Emma zum Rechentrainer")
        println("============================")
        println("Welche Rechenart möchtest du trainiren?")
        println("a für Addition")
        println("s für Subtraktion")
        println("m für Multiplikation")
        println("d für Divison")
        println("f für Finde den Platzhalter")
        println("r für Division mit Rest")

        v = chomp(readline())
        n = 10
        zeitLimit = 5 * n
        result = zeros(Int, n)
        result2 = zeros(Int, n)
        resultComp = zeros(Int, n)
        #flush(stdin)

        if v=="a"

                print("Addition:\n")
                a = rand(1:500, n)
                b = rand(1:500, n)

                startTime = Dates.Time(now())

                for i in (1:n)
                        print(string(a[i], " + ", b[i], " = "))
                        result[i] = parse(Int, chomp(readline()))
                        resultComp[i] = a[i] + b[i]
                end

        elseif v=="s"
                println("Subtraktion\n")

                a = rand(1:50, n)
                b = rand(1:50, n)

                startTime = Dates.Time(now())

                for i in (1:n)
                        if a[i] > b[i]
                                resultComp[i] = a[i] - b[i]
                                print(string(a[i], " - ", b[i], " = "))

                        else
                                resultComp[i] = b[i] - a[i]
                                print(string(b[i], " - ", a[i], " = "))

                        end

                        result[i] = parse(Int, chomp(readline()))

                end

        elseif v=="m"
                println("Multiplikation\n")
                a = rand(1:10, n)
                b = rand(1:10, n)

                startTime = Dates.Time(now())

                for i in (1:n)
                        print(string(a[i], " * ", b[i], " = "))
                        result[i] = parse(Int, chomp(readline()))
                        resultComp[i] = a[i] * b[i]
                end

        elseif v=="d"
                println("Divison")
                a = rand(1:10, n)
                b = rand(1:10, n)

                startTime = Dates.Time(now())

                for i in (1:n)
                        resultComp[i] = a[i] * b[i]
                        print(string(resultComp[i], " : ", a[i], " = "))
                        result[i] = parse(Int, chomp(readline()))
                end

        elseif v=="f"
                println("Finde den Platzhalter:")
                a = rand(1:10, n)
                b = rand(1:10, n)

                startTime = Dates.Time(now())

                for i in (1:n)
                        resultComp[i] = a[i] * b[i]
                        print(string(a[i], " * x = ", resultComp[i], "\t x = "))
                        result[i] = parse(Int, chomp(readline()))
                end

        elseif v=="r"
                println("Bestimme das Ergebnis der Division mit Rest:")
                a = rand(10:100, n)
                b = rand(1:9, n)
                c = zeros(Int64, n)

                startTime = Dates.Time(now())

                for i in (1:n)
                        c[i] = a[i] % b[i]
                        resultComp[i] = (a[i] - c[i]) / b[i]
                        print(string(a[i], " : ", b[i], " =  ?\t"))
                        result[i] = parse(Int, chomp(readline()))
                        print("Rest?\t")
                        result2[i] = parse(Int, chomp(readline()))
                end

        else
                println("Abbruch!!")
        end

    endTime = Dates.Time(now())
    #timeElapsed = Int64(round(Float64(endTime-startTime)/1000000000))
    timeElapsed = convert(Int64, round(Dates.value(endTime - startTime)/1000000000))

    if timeElapsed < zeitLimit

        zeitAnsage = string("Sehr gut, du warst um $(zeitLimit - timeElapsed) Sekunden schneller als erlaubt!")

    else

        zeitAnsage = string("Schade!! Du warst um $(timeElapsed - zeitLimit) Sekunden zu langsam.")

    end

        if v=="a"||v=="s"||v=="m"||v=="t"

                text2 = string("Du hast $(sum(result .== resultComp)) von $n Rechnungaufgaben richtig gelöst!")
                println(text2)

        elseif v=="d"||v=="f"
                text2 = string("Du hast $(sum(result .== b)) von $n Rechnungaufgaben richtig gelöst!")
                #text2 = @sprintf "Du hast %i von %i Rechnungaufgaben richtig gelöst!" sum(result .== resultComp) n;
                println(text2)

        elseif v=="r"

                text2 = string("Du hast $(sum(result .== resultComp)) von $n Quotienten gefunden!")
                text3 = string("Du hast $(sum(result2 .== c)) von $n Resten der Division gefunden!")
                #text2 = @sprintf "Du hast %i von %i Quotienten gefunden!" sum(result .== resultComp) n;
                #text3 = @sprintf "Du hast %i von %i Resten der Division gefunden!" sum(result2 .== c) n;

                println(text2)
                println(text3)

        end

        println(zeitAnsage)

end


text:
module MyGame
    using GameDataManager

    export foo, bar

    function foo()
        println("Hello World, I will return a XLSXTable")
        return GameDataManager.loadtable("items")
    end

    function bar(what)

    end
end


text:
module Plankton

export plankton_advection!
export plankton_diffusion!
export plankton_update!
export generate_individuals, individuals
export find_inds!, find_NPT!, acc_counts!, calc_par!

using CUDA
using StructArrays
using Random
using KernelAbstractions: @index, @kernel, Event, MultiEvent, wait
using KernelAbstractions.Extras.LoopInfo: @unroll

using PlanktonIndividuals.Architectures: device, Architecture, GPU, CPU, rng_type
using PlanktonIndividuals.Grids
using PlanktonIndividuals.Diagnostics

using PlanktonIndividuals: AbstractMode, CarbonMode, QuotaMode, MacroMolecularMode

#####
##### generate individuals of multiple species
#####
function generate_individuals(params::Dict, arch::Architecture, Nsp, N, maxN, g::AbstractGrid, mode::AbstractMode)
    plank_names = Symbol[]
    plank_data=[]
    for i in 1:Nsp
        name = Symbol("sp"*string(i))
        plank = construct_plankton(arch, i, params, maxN, mode::AbstractMode)
        generate_plankton!(plank, N, g, arch, mode)
        push!(plank_names, name)
        push!(plank_data, plank)
    end
    planks = NamedTuple{Tuple(plank_names)}(plank_data)
    return individuals(planks)
end

include("Advection/Advection.jl")
include("QuotaMode/QuotaMode.jl")
include("CarbonMode/CarbonMode.jl")
include("MacroMolecularMode/MacroMolecularMode.jl")
include("utils.jl")

using .Advection
import .Quota
import .Carbon
import .MacroMolecular

#####
##### some workarounds for function names
#####
construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN, mode::MacroMolecularMode) =
    MacroMolecular.construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN)

construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN, mode::QuotaMode) =
    Quota.construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN)

construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN, mode::CarbonMode) =
    Carbon.construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN)

generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture, mode::MacroMolecularMode) =
    MacroMolecular.generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture)

generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture, mode::QuotaMode) =
    Quota.generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture)

generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture, mode::CarbonMode) =
    Carbon.generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture)

plankton_update!(plank, nuts, proc, p, plk, diags_spcs, ΔT, t, arch::Architecture, mode::MacroMolecularMode) =
    MacroMolecular.plankton_update!(plank, nuts, proc, p, plk, diags_spcs, ΔT, t, arch::Architecture, mode::AbstractMode)
plankton_update!(plank, nuts, proc, p, plk, nothing::Nothing, ΔT, t, arch::Architecture, mode::MacroMolecularMode) =
    MacroMolecular.plankton_update!(plank, nuts, proc, p, plk, nothing::Nothing, ΔT, t, arch::Architecture, mode::AbstractMode)

plankton_update!(plank, nuts, proc, p, plk, diags_spcs, ΔT, t, arch::Architecture, mode::QuotaMode) =
    Quota.plankton_update!(plank, nuts, proc, p, plk, diags_spcs, ΔT, t, arch::Architecture, mode::AbstractMode)
plankton_update!(plank, nuts, proc, p, plk, nothing::Nothing, ΔT, t, arch::Architecture, mode::QuotaMode) =
    Quota.plankton_update!(plank, nuts, proc, p, plk, nothing::Nothing, ΔT, t, arch::Architecture, mode::AbstractMode)

plankton_update!(plank, nuts, proc, p, plk, diags_spcs, ΔT, t, arch::Architecture, mode::CarbonMode) =
    Carbon.plankton_update!(plank, nuts, proc, p, plk, diags_spcs, ΔT, t, arch::Architecture, mode::AbstractMode)
plankton_update!(plank, nuts, proc, p, plk, nothing::Nothing, ΔT, t, arch::Architecture, mode::CarbonMode) =
    Carbon.plankton_update!(plank, nuts, proc, p, plk, nothing::Nothing, ΔT, t, arch::Architecture, mode::AbstractMode)

end


text:
@testset "Corrupted/incorrect files" begin
    @test_nowarn readc3d(joinpath(datadir, "sample18", "bad_parameter_section.c3d"))
end


text:
"""
    power_table(sim, alpha = 0.05)

Returns a row table with fields `coefname` and `power`  based on the proportion of
simulated p-values less than alpha, for `sim`, the output of `parametricbootstrap`.

"""
function power_table(sim::MixedModelBootstrap, alpha = 0.05)
# `se` represents the Bernoulli standard error of the number of trials below signifance using the (1 - alpha) for the probability.
# q = 1 - alpha
# se = sqrt(alpha * q / nsim)
    nsim = length(sim.objective)
    keyvec = keys(first(sim.bstr).β)
    # if you don't initialize everything, then coefficents with zero power will be missing
    dd = Dict(coef => 0 for coef in keyvec)
    for row in sim.coefpvalues
        coef = row.coefname
        if row.p < alpha
            dd[coef] += 1
        end
    end
    return [ (; :coefname => string(key), :power => dd[key] / nsim) for key in keyvec ]
end


text:
using CairoMakie, Distributions, Random
Random.seed!(124)
cmap = to_colormap(:linear_wcmr_100_45_c42_n256)
alpha = rand(length(cmap))
cmap_alpha = RGBAf.(cmap, alpha)

set_theme!(theme_dark())
fig = Figure(resolution = (600,400))
ax = Axis(fig[1,1]; )
hist!(rand(Normal(2.6,0.4), 1000), normalization = :pdf, offset = -1,
    color = :values, colormap = :plasma, direction = :x, fillto = -0.5)
hist!(rand(Normal(2.6,0.4), 1000), normalization = :pdf, offset = 1,
    color = :grey90, direction = :x, scale_to = -0.5)
hist!(rand(Normal(0.2,0.2), 1000), normalization = :pdf, offset = 4,
    strokewidth = 1, color = :transparent, strokecolor = :white,
    direction = :y, scale_to = -1)
hist!(rand(Normal(0.2,0.2), 1000), normalization = :pdf, color = :values,
    colormap = cmap_alpha, strokewidth = 1, strokecolor = :black,
    bar_labels = :values, label_color = :white, label_size = 12,
    label_formatter=x-> round(x, digits=1))
hidedecorations!(ax; grid = false)
save("hist.png", fig, px_per_unit = 2)
display(fig)