text:
#
# sparse.jl --
#
# Implement various format of compressed sparse linear operators.  The purpose
# of this module is mainly to convert between different formats and to provide
# LazyAlgebra wrappers to apply the corresponding linear mappings.  Julia's
# SparseArrays standard package only provides "Compressed Sparse Column" (CSC)
# format.
#
# See https://en.wikipedia.org/wiki/Sparse_matrix.
#
#------------------------------------------------------------------------------
#
# 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.
#
module SparseOperators

export
    CompressedSparseOperator,
    SparseOperator,
    SparseOperatorCOO,
    SparseOperatorCSC,
    SparseOperatorCSR,
    nrows,
    ncols,
    row_size,
    col_size,
    nonzeros,
    nnz

using StructuredArrays
using ZippedArrays

import LinearAlgebra

using ..LazyAlgebra
using ..Foundations
using ..LazyAlgebra: @certify

import .LazyAlgebra:
    MorphismType,
    apply!,
    vcreate,
    identical,
    coefficients,
    row_size,
    col_size,
    nrows,
    ncols,
    input_ndims,
    input_size,
    output_ndims,
    output_size

import SparseArrays
using SparseArrays: SparseMatrixCSC, nonzeros, nnz
if isdefined(SparseArrays, :AbstractSparseMatrixCSC)
    const AbstractSparseMatrixCSC{Tv,Ti} =
        SparseArrays.AbstractSparseMatrixCSC{Tv,Ti}
else
    const AbstractSparseMatrixCSC{Tv,Ti} = SparseArrays.SparseMatrixCSC{Tv,Ti}
end

import Base: getindex, setindex!, iterate
using Base: @propagate_inbounds

#------------------------------------------------------------------------------
# Convert to integer type suitable for indexing.
to_int(i::Int) = i
to_int(i::Integer) = Int(i)

# Convert to vector of indices.
to_indices(inds::AbstractVector{<:Integer}) = to_values(Int, inds)

# Convert to vector of values with given element type and make sure it is a
# fast vector.
to_values(vals::AbstractVector{T}) where {T} = to_values(T, vals)
to_values(::Type{Any}, vals::AbstractVector{T}) where {T} = to_values(T, vals)
to_values(::Type{T}, vals::Vector{T}) where {T} = vals
to_values(::Type{T}, vals::AbstractVector) where {T} = convert(Vector{T}, vals)
@inline to_values(::Type{T}, vals::AbstractVector{T}) where {T} =
    _to_values(T, vals, eachindex(vals))

@inline _to_values(::Type{T}, vals::AbstractVector, inds) where {T} =
    convert(Vector{T}, vals) # Convert because not a fast vector.

@inline function _to_values(::Type{T}, vals::AbstractVector,
                            inds::AbstractUnitRange{Int}) where {T}
    (first(inds) == 1 ? vals : convert(Vector{T}, vals))
end

# Union of types acceptable to define array size and methods to convert to
# canonic form.
const ArraySize = Union{Integer,Tuple{Vararg{Integer}}}
to_size(siz::Tuple{Vararg{Int}}) = siz
to_size(siz::Tuple{Vararg{Integer}}) = map(to_int, siz)
to_size(siz::Integer) = (to_int(siz),)

as_matrix(A::AbstractMatrix, nrows::Int, ncols::Int) = begin
    @certify size(A) == (nrows, ncols)
    return A
end
as_matrix(A::AbstractArray, nrows::Int, ncols::Int) =
    reshape(A, (nrows, ncols))

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

"""

`SparseOperator{T,M,N}` is the abstract type inherited by the sparse operator
types. Parameter `T` is the type of the elements.  Parameters `M` and `N` are
the number of dimensions of the *rows* and of the *columns* respectively.
Sparse operators are a generalization of sparse matrices in the sense that they
implement linear mappings which can be applied to `N`-dimensonal arguments to
produce `M`-dimensional results (as explained below).  See
[`GeneralMatrix`](@ref) for a similar generalization but for *dense* matrices.

See [`CompressedSparseOperator`](@ref) for usage of sparse operators
implementing compressed storage formats.

"""
abstract type SparseOperator{T,M,N} <: LinearMapping end

"""

`CompressedSparseOperator{F,T,M,N}` is an abstract sub-type of
`SparseOperator{T,M,N}` and is inherited by the concrete types implementing
sparse operators with compressed storage in format `F`.

Format `F` is specificed as a symbol and can be:

- `:COO` for *Compressed Sparse Coordinate* storage format.  This format is not
  the most efficient, it is mostly used as an intermediate for building a
  sparse operator in one of the following formats.

- `:CSC` for *Compressed Sparse Column* storage format.  This format is very
  efficient for applying the adjoint of the sparse operator.

- `:CSR` for *Compressed Sparse Row* storage format.  This format is very
  efficient for directly applying the sparse operator.

To construct (or convert to) a sparse operator with compressed storage format
`F`, you can call:

    CompressedSparseOperator{F}(args...; kwds...)
    CompressedSparseOperator{F,T}(args...; kwds...)
    CompressedSparseOperator{F,T,M}(args...; kwds...)
    CompressedSparseOperator{F,T,M,N}(args...; kwds...)

where given parameters `T`, `M` and `N`, arguments `args...` and optional
keywords `kwds...` will be passed to the concrete constructor
[`SparseOperatorCOO`](@ref), [`SparseOperatorCSC`](@ref) or
[`SparseOperatorCSR`](@ref) corresponding to the format `F`.

It is possible to use a compressed sparse operator `A` as an iterator:

```julia
for (Aij,i,j) in A # simple but slow for CSR and CSC
    ...
end
```

to retrieve the values `Aij` and respective row `i` and column `j` indices for
all the entries stored in `A`.  It is however more efficient to access them
according to their storage order which depends on the compressed format.

- If `A` is in CSC format:

  ```julia
  using LazyAlgebra.SparseMethods
  for j in each_col(A)        # loop over column index
      for k in each_off(A, j) # loop over structural non-zeros in this column
          i   = get_row(A, k) # get row index of entry
          Aij = get_val(A, k) # get value of entry
       end
  end
  ```

- If `A` is in CSR format:

  ```julia
  using LazyAlgebra.SparseMethods
  for i in each_row(A)        # loop over row index
      for k in each_off(A, i) # loop over structural non-zeros in this row
          j   = get_col(A, k) # get column index of entry
          Aij = get_val(A, k) # get value of entry
       end
  end
  ```

- If `A` is in COO format:

  ```julia
  using LazyAlgebra.SparseMethods
  for k in each_off(A)
       i   = get_row(A, k) # get row index of entry
       j   = get_col(A, k) # get column index of entry
       Aij = get_val(A, k) # get value of entry
  end
  ```

The low-level methods `each_row`, `each_col`, `each_off`, `get_row`, `get_col`
and `get_val` are not automatically exported by `LazyAlgebra`, this is the
purpose of the statement `using LazyAlgebra.SparseMethods`.

"""
abstract type CompressedSparseOperator{F,T,M,N} <: SparseOperator{T,M,N} end

struct SparseOperatorCSR{T,M,N,
                         V<:AbstractVector{T},
                         J<:AbstractVector{Int},
                         K<:AbstractVector{Int}
                         } <: CompressedSparseOperator{:CSR,T,M,N}
    m::Int                # equivalent number of rows of the operator
    n::Int                # number of columns of the operator
    vals::V               # values of entries
    cols::J               # linear column indices of entries
    offs::K               # row offsets in arrays of entries and column indices
    rowsiz::NTuple{M,Int} # dimensions of rows
    colsiz::NTuple{N,Int} # dimensions of columns

    # An inner constructor is defined to prevent Julia from providing a simple
    # outer constructor, it is not meant to be called directly as it does not
    # check whether arguments are correct.
    function SparseOperatorCSR{T,M,N,V,J,K}(m::Int,
                                            n::Int,
                                            vals::V,
                                            cols::J,
                                            offs::K,
                                            rowsiz::NTuple{M,Int},
                                            colsiz::NTuple{N,Int}) where {
                                                T,M,N,
                                                V<:AbstractVector{T},
                                                J<:AbstractVector{Int},
                                                K<:AbstractVector{Int}}
        new{T,M,N,V,J,K}(m, n, vals, cols, offs, rowsiz, colsiz)
    end
end

struct SparseOperatorCSC{T,M,N,
                         V<:AbstractVector{T},
                         I<:AbstractVector{Int},
                         K<:AbstractVector{Int}
                         } <: CompressedSparseOperator{:CSC,T,M,N}
    m::Int                # equivalent number of rows of the operator
    n::Int                # number of columns of the operator
    vals::V               # values of entries
    rows::I               # linear row indices of entries
    offs::K               # columns offsets in arrays of entries and row indices
    rowsiz::NTuple{M,Int} # dimensions of rows
    colsiz::NTuple{N,Int} # dimensions of columns

    # An inner constructor is defined to prevent Julia from providing a simple
    # outer constructor, it is not meant to be called directly as it does not
    # check whether arguments are correct.
    function SparseOperatorCSC{T,M,N,V,I,K}(m::Int,
                                            n::Int,
                                            vals::V,
                                            rows::I,
                                            offs::K,
                                            rowsiz::NTuple{M,Int},
                                            colsiz::NTuple{N,Int}) where {
                                                T,M,N,
                                                V<:AbstractVector{T},
                                                I<:AbstractVector{Int},
                                                K<:AbstractVector{Int}}
        new{T,M,N,V,I,K}(m, n, vals, rows, offs, rowsiz, colsiz)
    end
end

struct SparseOperatorCOO{T,M,N,
                         V<:AbstractVector{T},
                         I<:AbstractVector{Int},
                         J<:AbstractVector{Int}
                         } <: CompressedSparseOperator{:COO,T,M,N}
    m::Int                # equivalent number of rows of the operator
    n::Int                # number of columns of the operator
    vals::V               # values of entries
    rows::I               # linear row indices of entries
    cols::J               # linear column indices of entries
    rowsiz::NTuple{M,Int} # dimensions of rows
    colsiz::NTuple{N,Int} # dimensions of columns

    # An inner constructor is defined to prevent Julia from providing a simple
    # outer constructor, it is not meant to be called directly as it does not
    # check whether arguments are correct.
    function SparseOperatorCOO{T,M,N,V,I,J}(m::Int,
                                            n::Int,
                                            vals::V,
                                            rows::I,
                                            cols::J,
                                            rowsiz::NTuple{M,Int},
                                            colsiz::NTuple{N,Int}) where {
                                                T,M,N,
                                                V<:AbstractVector{T},
                                                I<:AbstractVector{Int},
                                                J<:AbstractVector{Int}}
        new{T,M,N,V,I,J}(m, n, vals, rows, cols, rowsiz, colsiz)
    end
end

# Unions of compressed sparse operators that can be considered as being in a
# given storage format.

const AnyCSR{T,M,N} = Union{CompressedSparseOperator{:CSR,T,M,N},
                            Adjoint{<:CompressedSparseOperator{:CSC,T,M,N}}}

const AnyCSC{T,M,N} = Union{CompressedSparseOperator{:CSC,T,M,N},
                            Adjoint{<:CompressedSparseOperator{:CSR,T,M,N}}}

const AnyCOO{T,M,N} = Union{CompressedSparseOperator{:COO,T,M,N},
                            Adjoint{<:CompressedSparseOperator{:COO,T,M,N}}}

#------------------------------------------------------------------------------
# Accessors and basic methods.

nrows(A::SparseOperator) = getfield(A, :m)
ncols(A::SparseOperator) = getfield(A, :n)
row_size(A::SparseOperator) = getfield(A, :rowsiz)
col_size(A::SparseOperator) = getfield(A, :colsiz)
output_size(A::SparseOperator) = row_size(A)
input_size(A::SparseOperator) = col_size(A)
output_ndims(A::SparseOperator{T,M,N}) where {T,M,N} = M
input_ndims(A::SparseOperator{T,M,N}) where {T,M,N} = N

Base.eltype(A::SparseOperator{T,M,N}) where {T,M,N} = T
Base.ndims(A::SparseOperator{T,M,N}) where {T,M,N} = M+N
Base.length(A::SparseOperator) = nrows(A)*ncols(A)
Base.size(A::SparseOperator) = (row_size(A)..., col_size(A)...)

# Use constructors to perform conversion (the first method is to resolve
# ambiguities).
Base.convert(::Type{T}, A::T) where {T<:SparseOperator} = A
Base.convert(::Type{T}, A) where {T<:SparseOperator} = T(A)

# FIXME: This cannot be considered as a *pure* trait as it does not only
#        depend on the type of the object.
MorphismType(A::SparseOperator) =
    (row_size(A) == col_size(A) ? Endomorphism() : Morphism())

coefficients(A::SparseOperator) = get_vals(A)

identical(A::T, B::T) where {T<:CompressedSparseOperator{:CSR}} =
    (get_vals(A) === get_vals(B) && get_cols(A) === get_cols(B) &&
     get_offs(A) === get_offs(B) &&
     row_size(A) == row_size(B) && col_size(A) == col_size(B))

identical(A::T, B::T) where {T<:CompressedSparseOperator{:CSC}} =
    (get_vals(A) === get_vals(B) && get_rows(A) === get_rows(B) &&
     get_offs(A) === get_offs(B) &&
     row_size(A) == row_size(B) && col_size(A) == col_size(B))

identical(A::T, B::T) where {T<:CompressedSparseOperator{:COO}} =
    (get_vals(A) === get_vals(B) && get_rows(A) === get_rows(B) &&
     get_cols(A) === get_cols(B) &&
     row_size(A) == row_size(B) && col_size(A) == col_size(B))

# Assume that a copy of a compressed sparse operator is to keep the same
# structure for the structural non-zeros but possibly change the values.  So
# only duplicate the value part.

Base.copy(A::SparseOperatorCSR{T,M,N}) where {T,M,N} =
    unsafe_csr(nrows(A), ncols(A), copy_vals(A), get_cols(A), get_offs(A),
               row_size(A), col_size(A))

Base.copy(A::SparseOperatorCSC{T,M,N}) where {T,M,N} =
    unsafe_csc(nrows(A), ncols(A), copy_vals(A), get_rows(A), get_offs(A),
               row_size(A), col_size(A))

Base.copy(A::SparseOperatorCOO{T,M,N}) where {T,M,N} =
    unsafe_coo(nrows(A), ncols(A), copy_vals(A), get_rows(A), get_cols(A),
               row_size(A), col_size(A))

# `findnz(A) -> I,J,V` yields the row and column indices and the values of the
# stored values in A`.
SparseArrays.findnz(A::SparseOperator) =
    (get_rows(A), get_cols(A), get_vals(A))

# Extend some methods in SparseArrays.  The "structural" non-zeros are the
# entries stored by the sparse structure, they may or not be equal to zero
# un-stored entries are always considered as being equal to zero.
SparseArrays.nonzeros(A::SparseOperator) = get_vals(A)
SparseArrays.nnz(A::SparseOperator) = length(nonzeros(A))

"""
    get_vals(A)

yields the array storing the values of the sparse operator `A`.  The returned
array is shared with `A`, call `copy_vals(A)` instead if you want to modify the
contents of the returned array with no side effects on `A`.

As a convenience, argument may also be the adjoint of a sparse operator:

    get_vals(A') -> get_vals(A)

which yields the **unmodified** values of `A`, hence the caller has to take the
conjugate of these values.  The method `get_val(A',k)` however takes care of
conjugating the values.

"""
get_vals(A::SparseOperator) = getfield(A, :vals)
get_vals(A::Adjoint{<:SparseOperator}) = get_vals(unveil(A))

"""
    copy_vals([T = eltype(A),] A) -> vals

yields a copy of the values of the entries in sparse operator `A` converted to
type `T`.  The result is a vector that is not shared by `A`, the caller may
thus modify its contents with no side effects on `A`.

"""
copy_vals(A::SparseOperator{T}) where {T} = copy_vals(T, A)
function copy_vals(::Type{T}, A::SparseOperator) where {T}
    vals = get_vals(A)
    copyto!(Vector{T}(undef, size(vals)), vals)
end

"""
    get_rows(A)

yields the row indices of the entries of the sparse operator `A`.  The returned
array may be shared with `A`, call `copy_rows(A)` instead if you want to modify
the contents of the returned array with no side effects on `A`.

"""
get_rows(A::SparseOperatorCSC) = getfield(A, :rows)
get_rows(A::SparseOperatorCOO) = getfield(A, :rows)
get_rows(A::CompressedSparseOperator{:CSR}) =
    copy_rows(A) # FIXME: yield an iterator
get_rows(A::Adjoint{<:SparseOperator}) = get_cols(unveil(A))

"""
    copy_rows(A) -> rows

yields a copy of the linear row indices of entries in sparse operator `A`.  The
result is a vector that is not shared by `A`, the caller may thus modify its
contents with no side effects on `A`.

"""
function copy_rows(A::SparseOperator)
    rows = get_rows(A)
    copyto!(Vector{Int}(undef, size(rows)), rows)
end
function copy_rows(A::CompressedSparseOperator{:CSR})
    rows = Vector{Int}(undef, length(get_vals(A)))
    @inbounds for i in each_row(A)
        @simd for k in each_off(A, i)
            rows[k] = i
        end
    end
    return rows
end

"""
    get_cols(A)

yields the column indices of the entries of the sparse operator `A`.  The
returned array may be shared with `A`, call `copy_cols(A)` instead if you want
to modify the contents of the returned array with no side effects on `A`.

"""
get_cols(A::SparseOperatorCSR) = getfield(A, :cols)
get_cols(A::SparseOperatorCOO) = getfield(A, :cols)
get_cols(A::Union{CompressedSparseOperator{:CSC},SparseMatrixCSC}) =
    copy_cols(A) # FIXME: yield an iterator
get_cols(A::Adjoint{<:SparseOperator}) = get_rows(unveil(A))

"""
    copy_cols(A) -> cols

yields a copy of the linear column indices of entries in sparse operator `A`.
The result is a vector that is not shared by `A`, the caller may thus modify
its contents with no side effects on `A`.

"""
function copy_cols(A::SparseOperator)
    cols = get_cols(A)
    copyto!(Vector{Int}(undef, size(cols)), cols)
end
function copy_cols(A::Union{CompressedSparseOperator{:CSC},SparseMatrixCSC})
    cols = Vector{Int}(undef, length(get_vals(A)))
    @inbounds for j in each_col(A)
        @simd for k in each_off(A, j)
            cols[k] = j
        end
    end
    return cols
end

"""
    get_offs(A)

yields the table of offsets of the sparse operator `A`.  Not all operators
extend this method.

!!! warning
    The interpretation of offsets depend on the type of `A`.  For instance,
    assuming `offs = get_offs(A)`, then the index range of the `j`-th column of
    a `SparseMatrixCSC` is `offs[j]:(offs[j+1]-1)` while the index range is
    `(offs[j]+1):offs[j+1]` for a `SparseOperatorCSC`.  For this reason,
    it is recommended to call [`each_off`](@ref) instead or to call `get_offs`
    with 2 arguments as shown below.

For a transparent usage of the offsets, the method should be called with 2
arguments:

    get_offs(A, i) -> k1, k2

which yields the offsets of the first and last elements in the arrays of values
and linear column indices for the `i`-th row of the sparse operator `A` stored
in a *Compressed Sparse Row* (CSR) format.  If `k2 < k1`, it means that the
`i`-th row is empty.  Calling `each_off(A,i)` directly yields `k1:k2`.

    get_offs(A, j) -> k1, k2

yields the offsets of the first and last elements in the arrays of values and
linear row indices for the `j`-th column of the sparse operator `A` stored in a
*Compressed Sparse Column* (CSC) format.  If `k2 < k1`, it means that the
`j`-th column is empty.  Calling `each_off(A,j)` directly yields `k1:k2`.

"""
get_offs(A::SparseOperatorCSR) = getfield(A, :offs)
get_offs(A::SparseOperatorCSC) = getfield(A, :offs)
get_offs(A::Adjoint{<:CompressedSparseOperator{:CSR}}) = get_offs(unveil(A))
get_offs(A::Adjoint{<:CompressedSparseOperator{:CSC}}) = get_offs(unveil(A))

@inline function get_offs(A::AnyCSR, i::Int)
    offs = get_offs(A)
    @boundscheck ((i < 1)|(i ≥ length(offs))) && out_of_range_row_index(A, i)
    return ((@inbounds offs[i] + 1),
            (@inbounds offs[i+1]))
end

@inline function get_offs(A::AnyCSC, j::Int)
    offs = get_offs(A)
    @boundscheck ((j < 1)|(j ≥ length(offs))) && out_of_range_column_index(A, j)
    return ((@inbounds offs[j] + 1),
            (@inbounds offs[j+1]))
end

@noinline out_of_range_row_index(A, i::Integer) =
    throw(ErrorException(string("out of range row index ", i,
                                " for sparse operator with ", nrows(A),
                                " rows")))

@noinline out_of_range_column_index(A, j::Integer) =
    throw(ErrorException(string("out of range column index ", j,
                                " for sparse operator with ", ncols(A),
                                " columns")))

"""

For a sparse operator `A` stored in a *Compressed Sparse Coordinate* (COO)
format, the call:

    each_off(A)

yields an iterator over the indices in the arrays of values and of linear row
and column indices for the `k`-th entry of `A`.

---

For a sparse operator `A` stored in a *Compressed Sparse Column* (CSC) format,
the call:

    each_off(A, j)

yields an iterator over the indices in the arrays of values and linear row
indices for the `j`-th column of `A`.

---

For a sparse operator `A` stored in a *Compressed Sparse Row* (CSR) format, the
call:

    each_off(A, i)

yields an iterator over the indices in the arrays of values and linear column
indices for the `i`-th row of `A`.

"""
@inline each_off(A::CompressedSparseOperator{:COO}) = Base.OneTo(nnz(A))
@inline each_off(A::Adjoint{<:CompressedSparseOperator{:COO}}) =
    each_off(unveil(A))

@inline @propagate_inbounds function each_off(A::AnyCSR, i::Int)
    k1, k2 = get_offs(A, i)
    return k1:k2
end

@inline @propagate_inbounds function each_off(A::AnyCSC, j::Int)
    k1, k2 = get_offs(A, j)
    return k1:k2
end

"""
    each_row(A)

yields an iterator over the linear row indices of the sparse operator `A`
stored in a *Compressed Sparse Row* (CSR) format, this includes the adjoint of
a sparse operator in *Compressed Sparse Column* (CSC) format.

"""
each_row(A::CompressedSparseOperator{:CSR}) = Base.OneTo(nrows(A))
each_row(A::Adjoint{<:CompressedSparseOperator{:CSC}}) = each_col(unveil(A))

"""
    each_col(A)

yields an iterator over the linear column indices of the sparse operator `A`
stored in a *Compressed Sparse Column* (CSC) format, this includes the adjoint
of a sparse operator in *Compressed Sparse Row* (CSR) format.

"""
each_col(A::CompressedSparseOperator{:CSC}) = Base.OneTo(ncols(A))
each_col(A::Adjoint{<:CompressedSparseOperator{:CSR}}) = each_row(unveil(A))

"""
    get_row(A, k) -> i

yields the linear row index of the `k`-th entry of the sparse operator `A`
stored in a *Compressed Sparse Column* (CSC) or *Coordinate* (COO) formats
(this includes adjoint of sparse operators in CSR format).

"""
@inline @propagate_inbounds function get_row(
    A::Union{CompressedSparseOperator{:COO},
             CompressedSparseOperator{:CSC},
             Adjoint{<:CompressedSparseOperator{:CSR}}}, k::Integer)
    get_rows(A)[k]
end

"""
    get_col(A, k) -> j

yields the linear column index of the `k`-th entry of the sparse operator `A`
stored in a *Compressed Sparse Row* (CSR) or *Coordinate* (COO) formats (this
includes adjoint of sparse operators in CSC format).

"""
@inline @propagate_inbounds function get_col(
    A::Union{
    CompressedSparseOperator{:COO},
        CompressedSparseOperator{:CSR},
        Adjoint{<:CompressedSparseOperator{:CSC}}}, k::Integer)
    get_cols(A)[k]
end

"""
    get_val(A, k) -> v

yields the value of the `k`-th entry of the sparse operator `A` stored in a
*Compressed Sparse Row* (CSR), *Compressed Sparse Column* (CSC) or *Coordinate*
(COO) format.

Argument may also be the adjoint of a sparse operator:

    get_val(A', k) -> conj(get_val(A, k))

"""
@inline @propagate_inbounds get_val(A, k::Integer) = get_vals(A)[k]
@inline @propagate_inbounds get_val(A::Adjoint, k::Integer) =
    conj(get_vals(A)[k])

"""
    set_val!(A, k, v) -> v

assigns `v` to the value of the `k`-th entry of the sparse operator `A` stored
in a *Compressed Sparse Row* (CSR), *Compressed Sparse Column* (CSC) or
*Coordinate* (COO) format.

"""
@inline @propagate_inbounds set_val!(A, k::Integer, v) = get_vals(A)[k] = v
@inline @propagate_inbounds set_val!(A::Adjoint, k::Integer, v) = begin
    get_vals(A)[k] = conj(v)
    return v
end

# Iterators to deliver (v,i,j).

@inline function Base.iterate(A::AnyCSR, state::Tuple{Int,Int,Int} = (0,0,0))
    i, k, kmax = state
    @inbounds begin
        k += 1
        while k > kmax
            if i ≥ nrows(A)
                return nothing
            end
            i += 1
            k, kmax = get_offs(A, i)
        end
        v = get_val(A, k)
        j = get_col(A, k)
        return ((v, i, j), (i, k, kmax))
    end
end

@inline function Base.iterate(A::AnyCSC, state::Tuple{Int,Int,Int} = (0,0,0))
    j, k, kmax = state
    @inbounds begin
        k += 1
        while k > kmax
            if j ≥ ncols(A)
                return nothing
            end
            j += 1
            k, kmax = get_offs(A, j)
        end
        v = get_val(A, k)
        i = get_row(A, k)
        return ((v, i, j), (j, k, kmax))
    end
end

@inline function Base.iterate(A::AnyCOO, state::Tuple{Int,Int} = (0, nnz(A)))
    k, kmax = state
    @inbounds begin
        if k < kmax
            k += 1
            return ((get_val(A, k), get_row(A, k), get_col(A, k)), (k, kmax))
        else
            return nothing
        end
    end
end

#------------------------------------------------------------------------------
# Extend methods for SparseMatrixCSC defined in SparseArrays.

nrows(A::SparseMatrixCSC) = getfield(A, :m)
ncols(A::SparseMatrixCSC) = getfield(A, :n)
get_vals(A::SparseMatrixCSC) = getfield(A, :nzval)
get_offs(A::SparseMatrixCSC) = getfield(A, :colptr)
get_rows(A::SparseMatrixCSC) = getfield(A, :rowval)
# get_cols is already done elsewhere.
row_size(A::SparseMatrixCSC) = (nrows(A),)
col_size(A::SparseMatrixCSC) = (ncols(A),)
each_col(A::SparseMatrixCSC) = Base.OneTo(ncols(A))

@inline @propagate_inbounds each_off(A::SparseMatrixCSC, j::Integer) =
    ((k1, k2) = get_offs(A, j); k1:k2)

# Provide a specific version of `get_offs(A,j)` because offsets have a slightly
# different definition than our CSC format.
@inline function get_offs(A::SparseMatrixCSC, j::Integer)
    offs = get_offs(A)
    @boundscheck ((j < 1)|(j ≥ length(offs))) && out_of_range_column_index(A, j)
    return ((@inbounds offs[j]),
            (@inbounds offs[j+1]-1))
end

#------------------------------------------------------------------------------
# Constructors.

"""

Sparse operators in *Compressed Sparse Coordinate* (COO) format store the
significant entries in no particular order, as a vector of values, a vector of
linear row indices and a vector of linear column indices. It is even possible
to have repeated entries. This format is very useful to build a sparse linear
operator. It can be converted to a more efficient format like *Compressed
Sparse Column* (CSC) or *Compressed Sparse Row* (CSR) for fast application of
the sparse linear mapping or of its adjoint.

A sparse operator in COO storage format can be constructed by providing all
necessary information:

    SparseOperatorCOO(vals, rows, cols, rowsiz, colsiz)

where `vals` is the vector of values of the sparse entries, `rows` and `cols`
are integer valued vectors with the linear row and column indices of the sparse
entries, `rowsiz` and `colsiz` are the sizes of the row and column dimensions.
The entries values and respective linear row and column indices of the `k`-th
sparse entry are given by `vals[k]`, `rows[k]` and `cols[k]`.  For efficiency
reasons, sparse operators are currently limited to *fast* arrays because they
can be indexed linearly with no loss of performances.  If `vals`, `rows` and/or
`cols` are not fast arrays, they will be automatically converted to linearly
indexed arrays.

A sparse operator in COO storage format can be directly constructed from a
2-dimensional Julia array `A`:

    SparseOperatorCOO(A, sel = (v,i,j) -> (v != zero(v)))

where optional argument `sel` is a selector function which is called as
`sel(v,i,j)` with `v`, `i` and `j` the value, the row and the column linear
indices for each entries of `A` and which is assumed to yield `true` for the
entries of `A` to be selected in the sparse structure and `false` for the
entries of `A` to discard.  The default selector is such that all non-zeros of
`A` are selected.

The element type, say `T`, for the sparse coefficients can be imposed by
rewriting the above examples as:

    SparseOperatorCOO{T}(args...)

A sparse operator in COO storage format implementing generalized matrix-vector
multiplication can also be directly constructed from a `L`-dimensional Julia
array (with `L ≥ 2`) `A` by:

    SparseOperatorCOO{T,M}(A[, sel])

with `M` the number of leading dimensions of `A` corresponding to the *rows* of
the operator, the trailing `N = L - M` dimensions being assumed to correspond
to the *columns* of the operator.  These dimensions are the size of,
respectively, the output and the input arrays when applying the operator.  The
parameter `N` may be specified (although it can be automatically determined):

    SparseOperatorCOO{T,M,N}(A[, sel])

provided the equality `M + N = ndims(A)` holds.

A last parameter `V` can be specified for the type of the vector to store the
coefficients of the sparse operator:

    SparseOperatorCOO{T,M,N,V}(args...)

provided `V` implements standard linear indexing.  The default is to take `V =
Vector{T}`.  As a special case, you can choose a uniform boolean vector from
the `StructuredArrays` package to store the sparse coefficients:

    SparseOperatorCOO{T,M,N,UniformVector{Bool}}(args...)

to get a compressed sparse operator in COO format whose values are an immutable
uniform vector of true values requiring no storage.  This is useful to only
store the sparse structure of the operator, that is the indices in COO format
of the sparse coefficients not their values.

The `SparseOperatorCOO` constructor can also be used to convert a sparse
operator in another storage format into the COO format.  In that case,
parameter `T` may also be specified to convert the type of the sparse
coefficients.

""" SparseOperatorCOO

"""

Sparse operators in *Compressed Sparse Column* (CSC) format store the
significant entries in a column-wise order, as a vector of values, a vector of
corresponding linear row indices and a vector of offsets indicating, for each
column, the range of indices in the vectors of values and of row indices.  This
storage format is very suitable for fast application of the operator, notably
its adjoint.

A sparse operator in CSC storage format can be constructed by providing all
necessary information:

    SparseOperatorCSC(vals, rows, offs, rowsiz, colsiz)

where `vals` is the vector of values of the sparse entries, `rows` is an
integer valued vector of the linear row indices of the sparse entries, `offs`
is a column-wise table of offsets in these arrays, `rowsiz` and `colsiz` are
the sizes of the row and column dimensions.  The entries values and respective
linear row indices of the `j`-th column are given by `vals[k]` and `rows[k]`
with `k ∈ offs[j]+1:offs[j+1]`.  The linear column index `j` is in the range
`1:n` where `n = prod(colsiz)` is the equivalent number of columns.  For
efficiency reasons, sparse operators are currently limited to *fast* arrays
because they can be indexed linearly with no loss of performances.  If `vals`,
`rows` and/or `offs` are not fast arrays, they will be automatically converted
to linearly indexed arrays.

A sparse operator in CSC storage format can be directly constructed from a
2-dimensional Julia array `A`:

    SparseOperatorCSC(A, sel = (v,i,j) -> (v != zero(v)))

where optional argument `sel` is a selector function which is called as
`sel(v,i,j)` with `v`, `i` and `j` the value, the row and the column linear
indices for each entries of `A` and which is assumed to yield `true` for the
entries of `A` to be selected in the sparse structure and `false` for the
entries of `A` to discard.  The default selector is such that all non-zeros of
`A` are selected.

The element type, say `T`, for the sparse coefficients can be imposed by
rewriting the above examples as:

    SparseOperatorCSC{T}(args...)

A sparse operator in CSC storage format implementing generalized matrix-vector
multiplication can also be directly constructed from a `L`-dimensional Julia
array (with `L ≥ 2`) `A` by:

    SparseOperatorCSC{T,M}(A[, sel])

with `M` the number of leading dimensions of `A` corresponding to the *rows* of
the operator, the trailing `N = L - M` dimensions being assumed to correspond
to the *columns* of the operator.  These dimensions are the size of,
respectively, the output and the input arrays when applying the operator.  The
parameter `N` may be specified (although it can be automatically determined):

    SparseOperatorCSC{T,M,N}(A[, sel])

provided the equality `M + N = ndims(A)` holds.

A last parameter `V` can be specified for the type of the vector to store the
coefficients of the sparse operator:

    SparseOperatorCSC{T,M,N,V}(args...)

provided `V` implements standard linear indexing.  The default is to take `V =
Vector{T}`.  As a special case, you can choose a uniform boolean vector from
the `StructuredArrays` package to store the sparse coefficients:

    SparseOperatorCSC{T,M,N,UniformVector{Bool}}(args...)

to get a compressed sparse operator in CSC format whose values are an immutable
uniform vector of true values requiring no storage.  This is useful to only
store the sparse structure of the operator, that is the indices in CSC format
of the sparse coefficients not their values.

The `SparseOperatorCSC` constructor can also be used to convert a sparse
operator in another storage format into the CSC format.  In that case,
parameter `T` may also be specified to convert the type of the sparse
coefficients.

""" SparseOperatorCSC

"""

Sparse operators in *Compressed Sparse Row* (CSR) format store the significant
entries in a row-wise order, as a vector of values, a vector of corresponding
linear column indices and a vector of offsets indicating, for each row, the
range of indices in the vectors of values and of column indices.  This storage
format is very suitable for fast application of the operator.

A sparse operator in CSR storage format can be constructed by providing all
necessary information:

    SparseOperatorCSR(vals, cols, offs, rowsiz, colsiz)

where `vals` is the vector of values of the sparse entries, `cols` is an
integer valued vector of the linear column indices of the sparse entries,
`offs` is a column-wise table of offsets in these arrays, `rowsiz` and `colsiz`
are the sizes of the row and column dimensions.  The entries values and
respective linear column indices of the `i`-th row are given by `vals[k]` and
`cols[k]` with `k ∈ offs[i]+1:offs[i+1]`.  The linear row index `i` is in the
range `1:m` where `m = prod(rowsiz)` is the equivalent number of rows.  For
efficiency reasons, sparse operators are currently limited to *fast* arrays
because they can be indexed linearly with no loss of performances.  If `vals`,
`cols` and/or `offs` are not fast arrays, they will be automatically converted
to linearly indexed arrays.

A sparse operator in CSR storage format can be directly constructed from a
2-dimensional Julia array `A`:

    SparseOperatorCSR(A, sel = (v,i,j) -> (v != zero(v)))

where optional argument `sel` is a selector function which is called as
`sel(v,i,j)` with `v`, `i` and `j` the value, the row and the column linear
indices for each entries of `A` and which is assumed to yield `true` for the
entries of `A` to be selected in the sparse structure and `false` for the
entries of `A` to discard.  The default selector is such that all non-zeros of
`A` are selected.

The element type, say `T`, for the sparse coefficients can be imposed by
rewriting the above examples as:

    SparseOperatorCSR{T}(args...)

A sparse operator in CSR storage format implementing generalized matrix-vector
multiplication can also be directly constructed from a `L`-dimensional Julia
array (with `L ≥ 2`) `A` by:

    SparseOperatorCSR{T,M}(A[, sel])

with `M` the number of leading dimensions of `A` corresponding to the *rows* of
the operator, the trailing `N = L - M` dimensions being assumed to correspond
to the *columns* of the operator.  These dimensions are the size of,
respectively, the output and the input arrays when applying the operator.  The
parameter `N` may be specified (although it can be automatically determined):

    SparseOperatorCSR{T,M,N}(A[, sel])

provided the equality `M + N = ndims(A)` holds.

A last parameter `V` can be specified for the type of the vector to store the
coefficients of the sparse operator:

    SparseOperatorCSR{T,M,N,V}(args...)

provided `V` implements standard linear indexing.  The default is to take `V =
Vector{T}`.  As a special case, you can choose a uniform boolean vector from
the `StructuredArrays` package to store the sparse coefficients:

    SparseOperatorCSR{T,M,N,UniformVector{Bool}}(args...)

to get a compressed sparse operator in CSR format whose values are an immutable
uniform vector of true values requiring no storage.  This is useful to only
store the sparse structure of the operator, that is the indices in CSR format
of the sparse coefficients not their values.

The `SparseOperatorCSR` constructor can also be used to convert a sparse
operator in another storage format into the CSR format.  In that case,
parameter `T` may also be specified to convert the type of the sparse
coefficients.

""" SparseOperatorCSR

# Make sparse operators callable.
@callable SparseOperatorCSR
@callable SparseOperatorCSC
@callable SparseOperatorCOO

SparseOperator(A::SparseOperator) = A
SparseOperator{T}(A::SparseOperator{T}) where {T} = A
SparseOperator{T,M}(A::SparseOperator{T,M}) where {T,M} = A
SparseOperator{T,M,N}(A::SparseOperator{T,M,N}) where {T,M,N} = A

# Change element type.
SparseOperator{T}(A::SparseOperator{<:Any,M,N}) where {T,M,N} =
    SparseOperator{T,M,N}(A)
SparseOperator{T,M}(A::SparseOperator{<:Any,M,N}) where {T,M,N} =
    SparseOperator{T,M,N}(A)
for F in (:SparseOperatorCSC, :SparseOperatorCSR, :SparseOperatorCOO)
    @eval begin
        SparseOperator{T,M,N}(A::$F{<:Any,M,N}) where {T,M,N} =
            $F{T,M,N}(A)
    end
end

for (fmt,func) in ((:CSC, :SparseOperatorCSC),
                   (:CSR, :SparseOperatorCSR),
                   (:COO, :SparseOperatorCOO),)
    F = Expr(:quote, Symbol(fmt))
    @eval begin
        CompressedSparseOperator{$F}(args...; kwds...) =
            $func(args...; kwds...)
        CompressedSparseOperator{$F,T}(args...; kwds...) where {T} =
            $func{T}(args...; kwds...)
        CompressedSparseOperator{$F,T,M}(args...; kwds...) where {T,M} =
            $func{T,M}(args...; kwds...)
        CompressedSparseOperator{$F,T,M,N}(args...; kwds...) where {T,M,N} =
            $func{T,M,N}(args...; kwds...)
    end
end

# Conversion without changing format (mostly for changing element type).
CompressedSparseOperator{Any}(A::CompressedSparseOperator) = A
CompressedSparseOperator{Any,T}(A::CompressedSparseOperator{F}) where {F,T} =
    CompressedSparseOperator{F,T}(A)
CompressedSparseOperator{Any,T,M}(A::CompressedSparseOperator{F}) where {F,T,M} =
    CompressedSparseOperator{F,T,M}(A)
CompressedSparseOperator{Any,T,M,N}(A::CompressedSparseOperator{F}) where {F,T,M,N} =
    CompressedSparseOperator{F,T,M,N}(A)

# Basic outer constructors return a fully checked structure.

function SparseOperatorCSR(vals::AbstractVector,
                           cols::AbstractVector{<:Integer},
                           offs::AbstractVector{<:Integer},
                           rowsiz::ArraySize,
                           colsiz::ArraySize)
    check_structure(unsafe_csr(to_values(vals),
                               to_indices(cols),
                               to_indices(offs),
                               to_size(rowsiz),
                               to_size(colsiz)))
end

function SparseOperatorCSC(vals::AbstractVector,
                           rows::AbstractVector{<:Integer},
                           offs::AbstractVector{<:Integer},
                           rowsiz::ArraySize,
                           colsiz::ArraySize)
    check_structure(unsafe_csc(to_values(vals),
                               to_indices(rows),
                               to_indices(offs),
                               to_size(rowsiz),
                               to_size(colsiz)))
end

function SparseOperatorCOO(vals::AbstractVector,
                           rows::AbstractVector{<:Integer},
                           cols::AbstractVector{<:Integer},
                           rowsiz::ArraySize,
                           colsiz::ArraySize)
    check_structure(unsafe_coo(to_values(vals),
                               to_indices(rows),
                               to_indices(cols),
                               to_size(rowsiz),
                               to_size(colsiz)))
end

@inline isnonzero(v::T, i::Integer, j::Integer) where {T} = (v != zero(T))

for CS in (:SparseOperatorCSR,
           :SparseOperatorCSC,
           :SparseOperatorCOO)
    @eval begin
        # Get rid of the M,N parameters, but keep/set T for conversion of
        # values.
        $CS{T,M,N}(A::SparseOperator{<:Any,M,N}) where {T,M,N} = $CS{T}(A)
        $CS{T,M}(A::SparseOperator{<:Any,M}) where {T,M} = $CS{T}(A)
        $CS(A::SparseOperator{T}) where {T} = $CS{T}(A)

        # Cases which do nothing (it makes sense that a constructor of an
        # immutable type be able to just return its argument if it is already
        # of the correct type).
        $CS{T}(A::$CS{T}) where {T} = A

        # Manage to call constructors of compressed sparse operator given a
        # regular Julia array with correct parameters and selector.
        $CS(A::AbstractMatrix{T}, args...; kwds...) where {T} =
            $CS{T,1,1}(A, args...; kwds...)
        $CS{Any}(A::AbstractMatrix{T}, args...; kwds...) where {T} =
            $CS{T,1,1}(A, args...; kwds...)
        $CS{T}(A::AbstractMatrix, args...; kwds...) where {T} =
            $CS{T,1,1}(A, args...; kwds...)
        $CS{Any,M}(A::AbstractArray{T}, args...; kwds...) where {T,M} =
            $CS{T,M}(A, args...; kwds...)
        $CS{Any,M,N}(A::AbstractArray{T}, args...; kwds...) where {T,M,N} =
            $CS{T,M,N}(A, args...; kwds...)
        $CS{Any,M,N,V}(A::AbstractArray{T}, args...; kwds...) where {T,M,N,V} =
            $CS{T,M,N,V}(A, args...; kwds...)
        function $CS{T,M}(A::AbstractArray{S,L}, args...; kwds...) where {S,T,L,M}
            1 ≤ M < L ||
                error("parameters M=$M and L=$L are not such that 1 ≤ M < L")
            $CS{T,M,L-M}(A, args...; kwds...)
        end
        $CS{T,M,N}(A::AbstractArray, args...; kwds...) where {T,M,N} =
            $CS{T,M,N,Vector{T}}(A, args...; kwds...)
    end
end

# Constructors that convert array of values.  Other fields have already been
# checked so do not check structure again.

SparseOperatorCSR{T}(A::SparseOperatorCSR{S,M,N}) where {S,T,M,N} =
    unsafe_csr(nrows(A), ncols(A), to_values(T, get_vals(A)),
               get_cols(A), get_offs(A), row_size(A), col_size(A))

SparseOperatorCSC{T}(A::SparseOperatorCSC{S,M,N}) where {S,T,M,N} =
    unsafe_csc(nrows(A), ncols(A), to_values(T, get_vals(A)),
               get_rows(A), get_offs(A), row_size(A), col_size(A))

SparseOperatorCOO{T}(A::SparseOperatorCOO{S,M,N}) where {S,T,M,N} =
    unsafe_coo(nrows(A), ncols(A), to_values(T, get_vals(A)),
               get_rows(A), get_cols(A), row_size(A), col_size(A))

# Constructors for CSR format similar to the basic ones but have parameters
# that may imply converting arguments.

function SparseOperatorCSR{T,M,N}(vals::AbstractVector,
                                  cols::AbstractVector{<:Integer},
                                  offs::AbstractVector{<:Integer},
                                  rowsiz::ArraySize,
                                  colsiz::ArraySize) where {T,M,N}
    check_row_ndims(rowsiz, M)
    check_column_ndims(colsiz, N)
    SparseOperatorCSR{T}(vals, cols, offs, rowsiz, colsiz)
end

function SparseOperatorCSR{T,M}(vals::AbstractVector,
                                cols::AbstractVector{<:Integer},
                                offs::AbstractVector{<:Integer},
                                rowsiz::ArraySize,
                                colsiz::ArraySize) where {T,M}
    check_row_ndims(rowsiz, M)
    SparseOperatorCSR{T}(vals, cols, offs, rowsiz, colsiz)
end

function SparseOperatorCSR{T}(vals::AbstractVector,
                              cols::AbstractVector{<:Integer},
                              offs::AbstractVector{<:Integer},
                              rowsiz::ArraySize,
                              colsiz::ArraySize) where {T}
    SparseOperatorCSR(to_values(T, vals), cols, offs, rowsiz, colsiz)
end

# Idem for CSC format.

function SparseOperatorCSC{T,M,N}(vals::AbstractVector,
                                  rows::AbstractVector{<:Integer},
                                  offs::AbstractVector{<:Integer},
                                  rowsiz::ArraySize,
                                  colsiz::ArraySize) where {T,M,N}
    check_row_ndims(rowsiz, M)
    check_column_ndims(colsiz, N)
    SparseOperatorCSC{T}(vals, rows, offs, rowsiz, colsiz)
end

function SparseOperatorCSC{T,M}(vals::AbstractVector,
                                rows::AbstractVector{<:Integer},
                                offs::AbstractVector{<:Integer},
                                rowsiz::ArraySize,
                                colsiz::ArraySize) where {T,M}
    check_row_ndims(rowsiz, M)
    SparseOperatorCSC{T}(vals, rows, offs, rowsiz, colsiz)
end

function SparseOperatorCSC{T}(vals::AbstractVector,
                              rows::AbstractVector{<:Integer},
                              offs::AbstractVector{<:Integer},
                              rowsiz::ArraySize,
                              colsiz::ArraySize) where {T}
    SparseOperatorCSC(to_values(T, vals), rows, offs, rowsiz, colsiz)
end

# Idem for COO format.

function SparseOperatorCOO{T,M,N}(vals::AbstractVector,
                                  rows::AbstractVector{<:Integer},
                                  cols::AbstractVector{<:Integer},
                                  rowsiz::ArraySize,
                                  colsiz::ArraySize) where {T,M,N}
    check_row_ndims(rowsiz, M)
    check_column_ndims(colsiz, N)
    SparseOperatorCOO{T}(vals, rows, cols, rowsiz, colsiz)
end

function SparseOperatorCOO{T,M}(vals::AbstractVector,
                                rows::AbstractVector{<:Integer},
                                cols::AbstractVector{<:Integer},
                                rowsiz::ArraySize,
                                colsiz::ArraySize) where {T,M}
    check_row_ndims(rowsiz, M)
    SparseOperatorCOO{T}(vals, rows, cols, rowsiz, colsiz)
end

function SparseOperatorCOO{T}(vals::AbstractVector,
                              rows::AbstractVector{<:Integer},
                              cols::AbstractVector{<:Integer},
                              rowsiz::ArraySize,
                              colsiz::ArraySize) where {T}
    SparseOperatorCOO(to_values(T, vals), rows, cols, rowsiz, colsiz)
end

check_row_ndims(rowsiz::ArraySize, M::Integer) =
    length(rowsiz) == M ||
        throw_dimension_mismatch("number of row dimensions is not M=$M")

check_column_ndims(colsiz::ArraySize, N::Integer) =
    length(colsiz) == N ||
        throw_dimension_mismatch("number of columns dimensions is not N=$N")

# Constructors of a sparse operator in various format given a regular Julia
# array and a selector function.  Julia arrays are usually in column-major
# order but this is not always the case, to handle various storage orders when
# extracting selected entries, we convert the input array into a equivalent
# "matrix", that is a 2-dimensional array.

function SparseOperatorCSR{T,M,N,V}(arr::AbstractArray{S,L},
                                    sel::Function = isnonzero) where {
                                        S,T,L,M,N,V<:AbstractVector{T}}
    # Get equivalent matrix dimensions.
    nrows, ncols, rowsiz, colsiz = get_equivalent_size(arr, Val(M), Val(N))

    # Convert into equivalent matrix.
    A = as_matrix(arr, nrows, ncols)

    # Count the number of selected entries.
    nvals = count_selection(A, sel)

    # Extract the selected entries and their column indices and count the
    # numver of selected entries per row.  The pseudo-matrix is walked in
    # row-major order.
    cols = Vector{Int}(undef, nvals)
    offs = Vector{Int}(undef, nrows + 1)
    k = 0
    if V <: UniformVector{Bool}
        # Just extract the structure, not the values.
        vals = V(true, nvals)
        @inbounds for i in 1:nrows
            offs[i] = k
            for j in 1:ncols
                if sel(A[i,j], i, j)
                    (k += 1) ≤ nvals || bad_selector()
                    cols[k] = j
                end
            end
        end
    else
        # Extract the structure and the values.
        vals = V(undef, nvals)
        @inbounds for i in 1:nrows
            offs[i] = k
            for j in 1:ncols
                v = A[i,j]
                if sel(v, i, j)
                    (k += 1) ≤ nvals || bad_selector()
                    vals[k] = v
                    cols[k] = j
                end
            end
        end
    end
    k == nvals || bad_selector()
    offs[end] = nvals

    # By construction, the sparse structure should be correct so just call the
    # "unsafe" constructor.
    return unsafe_csr(nrows, ncols, vals, cols, offs, rowsiz, colsiz)
end

function SparseOperatorCSC{T,M,N,V}(arr::AbstractArray{S,L},
                                    sel::Function = isnonzero) where {
                                        S,T,L,M,N,V<:AbstractVector{T}}
    # Get equivalent matrix dimensions.
    nrows, ncols, rowsiz, colsiz = get_equivalent_size(arr, Val(M), Val(N))

    # Convert into equivalent matrix.
    A = as_matrix(arr, nrows, ncols)

    # Count the number of selected entries.
    nvals = count_selection(A, sel)

    # Extract the selected entries and their row indices and count the numver
    # of selected entries per column.  The pseudo-matrix is walked in
    # column-major order.
    rows = Vector{Int}(undef, nvals)
    offs = Vector{Int}(undef, ncols + 1)
    k = 0
    if V <: UniformVector{Bool}
        # Just extract the structure, not the values.
        vals = V(true, nvals)
        @inbounds for j in 1:ncols
            offs[j] = k
            for i in 1:nrows
                if sel(A[i,j], i, j)
                    (k += 1) ≤ nvals || bad_selector()
                    rows[k] = i
                end
            end
        end
    else
        # Extract the structure and the values.
        vals = V(undef, nvals)
        @inbounds for j in 1:ncols
            offs[j] = k
            for i in 1:nrows
                v = A[i,j]
                if sel(v, i, j)
                    (k += 1) ≤ nvals || bad_selector()
                    vals[k] = v
                    rows[k] = i
                end
            end
        end
    end
    k == nvals || bad_selector()
    offs[end] = nvals

    # By construction, the sparse structure should be correct so just call the
    # "unsafe" constructor.
    return unsafe_csc(nrows, ncols, vals, rows, offs, rowsiz, colsiz)
end

function SparseOperatorCOO{T,M,N,V}(arr::AbstractArray{S,L},
                                    sel::Function = isnonzero) where {
                                        S,T,L,M,N,V<:AbstractVector{T}}
    # Get equivalent matrix dimensions.
    nrows, ncols, rowsiz, colsiz = get_equivalent_size(arr, Val(M), Val(N))

    # Convert into equivalent matrix.
    A = as_matrix(arr, nrows, ncols)

    # Count the number of selected entries.
    nvals = count_selection(A, sel)

    # Extract the selected entries and their row and column indices.  The
    # pseudo-matrix is walked in column-major order since most Julia arrays are
    # stored in that order.
    rows = Vector{Int}(undef, nvals)
    cols = Vector{Int}(undef, nvals)
    k = 0
    if V <: UniformVector{Bool}
        # Just extract the structure, not the values.
        vals = V(true, nvals)
        @inbounds for j in 1:ncols
            for i in 1:nrows
                if sel(A[i,j], i, j)
                    (k += 1) ≤ nvals || bad_selector()
                    vals[k] = v
                    rows[k] = i
                    cols[k] = j
                end
            end
        end
    else
        # Extract the structure and the values.
        vals = V(undef, nvals)
        @inbounds for j in 1:ncols
            for i in 1:nrows
                v = A[i,j]
                if sel(v, i, j)
                    (k += 1) ≤ nvals || bad_selector()
                    vals[k] = v
                    rows[k] = i
                    cols[k] = j
                end
            end
        end
    end
    k == nvals || bad_selector()

    # By construction, the sparse structure should be correct so just call the
    # "unsafe" constructor.
    return unsafe_coo(nrows, ncols, vals, rows, cols, rowsiz, colsiz)
end

"""
    unpack!(A, S; flatten=false) -> A

unpacks the non-zero coefficients of the sparse operator `S` into the array `A`
and returns `A`.  Keyword `flatten` specifies whether to only consider the
length of `A` instead of its dimensions.  In any cases, `A` must have as many
elements as `length(S)` and standard linear indexing.

Just call `Array(S)` to unpack the coefficients of a sparse operator `S`
without providing the destination array.

""" unpack!

# Convert to standard Julia arrays which are stored in column-major order,
# hence the stride is the equivalent number of rows.  For COO format, as
# duplicates are allowed, values must be combined by an operator.

Base.Array(A::SparseOperator{T,M,N}) where {T,M,N} = Array{T,M+N}(A)
Base.Array{T}(A::SparseOperator{<:Any,M,N}) where {T,M,N} = Array{T,M+N}(A)
function Base.Array{T,L}(A::SparseOperator{<:Any,M,N}) where {T,L,M,N}
    L == M + N || throw_incompatible_number_of_dimensions()
    return unpack!(Array{T}(undef, (row_size(A)..., col_size(A)...,)), A)
end

function prepare_unpack!(dst::AbstractArray,
                         src::SparseOperator,
                         flatten::Bool)
    is_fast_array(dst) || throw_non_standard_indexing("destination array")
    if flatten
        length(dst) == length(src) ||
            throw_incompatible_number_of_elements()
    else
        size(dst) == (row_size(src)..., col_size(src)...,) ||
            throw_incompatible_dimensions()
    end
    fill!(dst, zero(eltype(dst)))
end

function unpack!(B::AbstractArray{T,L},
                 A::SparseOperatorCSR{<:Any,M,N};
                 flatten::Bool = false) where {T,L,M,N}
    prepare_unpack!(B, A, flatten)
    m = nrows(A) # used as the "stride" in B
    @inbounds for i in each_row(A)
        for k in each_off(A, i)
            j = get_col(A, k)
            v = get_val(A, k)
            B[i + m*(j - 1)] = v
        end
    end
    return B
end

function unpack!(B::AbstractArray{T,L},
                 A::SparseOperatorCSC{<:Any,M,N};
                 flatten::Bool = false) where {T,L,M,N}
    prepare_unpack!(B, A, flatten)
    m = nrows(A) # used as the "stride" in B
    @inbounds for j in each_col(A)
        for k in each_off(A, j)
            i = get_row(A, k)
            v = get_val(A, k)
            B[i + m*(j - 1)] = v
        end
    end
    return B
end

unpack!(B::AbstractArray, A::SparseOperatorCOO; kwds...) =
    unpack!(B, A, +; kwds...)

unpack!(B::AbstractArray, A::SparseOperatorCOO{Bool}; kwds...) =
    unpack!(B, A, |; kwds...)

function unpack!(B::AbstractArray{T,L},
                 A::SparseOperatorCOO{<:Any,M,N},
                 op::Function; flatten::Bool = false) where {T,L,M,N}
    prepare_unpack!(B, A, flatten)
    m = nrows(A) # used as the "stride" in B
    @inbounds for k in each_off(A)
        i = get_row(A, k)
        j = get_col(A, k)
        v = get_val(A, k)
        l = i + m*(j - 1)
        B[l] = op(B[l], v)
    end
    return B
end

function check_new_shape(A::SparseOperator,
                         rowsiz::Tuple{Vararg{Int}},
                         colsiz::Tuple{Vararg{Int}})
    prod(rowsiz) == nrows(A) ||
        bad_size("products of row dimensions must be equal")
    prod(colsiz) == ncols(A) ||
        bad_size("products of column dimensions must be equal")
end

Base.reshape(A::SparseOperator, rowsiz::ArraySize, colsiz::ArraySize) =
    reshape(A, to_size(rowsiz), to_size(colsiz))

function Base.reshape(A::SparseOperatorCSR,
                      rowsiz::Tuple{Vararg{Int}},
                      colsiz::Tuple{Vararg{Int}})
    check_new_shape(A, rowsiz, colsiz)
    unsafe_csr(nrows(A), ncols(A), get_vals(A), get_cols(A), get_offs(A),
               rowsiz, colsiz)
end

function Base.reshape(A::SparseOperatorCSC,
                      rowsiz::Tuple{Vararg{Int}},
                      colsiz::Tuple{Vararg{Int}})
    check_new_shape(A, rowsiz, colsiz)
    unsafe_csc(nrows(A), ncols(A), get_vals(A), get_rows(A), get_offs(A),
               rowsiz, colsiz)
end

function Base.reshape(A::SparseOperatorCOO,
                      rowsiz::Tuple{Vararg{Int}},
                      colsiz::Tuple{Vararg{Int}})
    check_new_shape(A, rowsiz, colsiz)
    unsafe_coo(nrows(A), ncols(A), get_vals(A), get_rows(A), get_cols(A),
               rowsiz, colsiz)
end

# Convert from other compressed sparse formats.  For compressed sparse row and
# column (CSR and CSC) formats, the compressed sparse coordinate (COO) format
# is used as an intermediate representation and entries are sorted in
# row/column major order.  To avoid side-effects, they must be copied first.
# Unless values are converted, there is no needs to copy when converting to a
# compressed sparse coordinate (COO) format.

SparseOperatorCSR{T}(A::SparseOperator) where {T} =
    coo_to_csr!(copy_vals(T, A),
                copy_rows(A),
                copy_cols(A),
                row_size(A),
                col_size(A))

SparseOperatorCSC{T}(A::SparseOperator) where {T} =
    coo_to_csc!(copy_vals(T, A),
                copy_rows(A),
                copy_cols(A),
                row_size(A),
                col_size(A))

SparseOperatorCOO{T}(A::SparseOperator) where {T} =
    SparseOperatorCOO(copy_vals(T, A),
                      get_rows(A),
                      get_cols(A),
                      row_size(A),
                      col_size(A))

SparseOperatorCOO{T}(A::SparseOperator{T}) where {T} =
    SparseOperatorCOO(get_vals(A),
                      get_rows(A),
                      get_cols(A),
                      row_size(A),
                      col_size(A))

"""
    coo_to_csr!(vals, rows, cols, rowsiz, colsiz [, mrg]) -> A

yields the a compressed sparse operator in a CSR format given the components
`vals`, `rows` and `cols` in the COO format and the sizes `rowsiz` and `colsiz`
of the row and column dimensions.  Input arrays are modified in-place.
Optional argument `mrg` is a function called to merge values of entries with
the same row and column indices.

Input arrays must be regular Julia vectors to ensure type stability in case of
duplicates.

"""
function coo_to_csr!(vals::Vector{T},
                     rows::Vector{Int},
                     cols::Vector{Int},
                     rowsiz::NTuple{M,Int},
                     colsiz::NTuple{N,Int},
                     mrg::Function = (T <: Bool ? (|) : (+))) where {T,M,N}
    # Check row and column sizes.
    nrows = check_size(rowsiz, "row")
    ncols = check_size(colsiz, "column")

    # Check row and column indices.
    check_rows(rows, nrows)
    check_cols(cols, ncols)

    # Sort and merge entries in row-major order, resize arrays if needed and
    # compute offsets.
    nvals = sort_and_merge!(vals, rows, cols, mrg)
    if nvals < length(vals)
        vals = vals[1:nvals]
        cols = cols[1:nvals]
    end
    offs = compute_offsets(nrows, rows, nvals)

    # Since everything will have
    # been checked, we can call the unsafe constructor.
    return unsafe_csr(nrows, ncols, vals, cols, offs, rowsiz, colsiz)
end

"""
    coo_to_csc!(vals, rows, cols, rowsiz, colsiz [, mrg]) -> A

yields the a compressed sparse operator in a CSC format given the components
`vals`, `rows` and `cols` in the COO format and the sizes `rowsiz` and `colsiz`
of the row and column dimensions.  Input arrays are modified in-place.
Optional argument `mrg` is a function called to merge values of entries with
the same row and column indices.

Input arrays must be regular Julia vectors to ensure type stability in case of
duplicates.

"""
function coo_to_csc!(vals::Vector{T},
                     rows::Vector{Int},
                     cols::Vector{Int},
                     rowsiz::NTuple{M,Int},
                     colsiz::NTuple{N,Int},
                     mrg::Function = (T <: Bool ? (|) : (+))) where {T,M,N}
    # Check row and column sizes.
    nrows = check_size(rowsiz, "row")
    ncols = check_size(colsiz, "column")

    # Check row and column indices.
    check_rows(rows, nrows)
    check_cols(cols, ncols)

    # Sort and merge entries in column-major order, resize arrays if needed and
    # compute offsets.
    nvals = sort_and_merge!(vals, cols, rows, mrg)
    if nvals < length(vals)
        vals = vals[1:nvals]
        rows = rows[1:nvals]
    end
    offs = compute_offsets(ncols, cols, nvals)

    # Since everything will have
    # been checked, we can call the unsafe constructor.
    return unsafe_csc(nrows, ncols, vals, rows, offs, rowsiz, colsiz)
end

# "less-than" method for sorting entries in order, arguments are 3-tuples
# `(v,major,minor)` with `v` the entry value, `major` the major index and
# `minir` the minor index.
@inline major_order(a::Tuple{T,Int,Int},b::Tuple{T,Int,Int}) where {T} =
    ifelse(a[2] == b[2], a[3] < b[3], a[2] < b[2])

"""
    sort_and_merge!(vals, major, minor, mrg) -> nvals

sorts entries and merges duplicates in input arrays `vals`, `major` and
`minor`.  Entries consist in the 3-tuples `(vals[i],major[i],minor[i])`.  The
sorting order of the `i`-th entry is based the value of `major[i]` and, if
equal, on the value of `minor[i]`.  After sorting, duplicates entries, that is
those which have the same `(major[i],minor[i])`, are removed merging the
associated values in `vals` with the `mrg` function.  All operations are done
in-place, the number of unique entries is returned but inputs arrays are not
resized, only the `nvals` first entries are valid.

"""
function sort_and_merge!(vals::AbstractVector,
                         major::AbstractVector{Int},
                         minor::AbstractVector{Int},
                         mrg::Function)
    # Sort entries in order (this also ensures that all arrays have the same
    # dimensions).
    sort!(ZippedArray(vals, major, minor); lt=major_order)

    # Merge duplicates.
    j = 1
    @inbounds for k in 2:length(vals)
        if major[k] == major[j] && minor[k] == minor[j]
            vals[j] = mrg(vals[j], vals[k])
        elseif (j += 1) < k
            vals[j], major[j], minor[j] = vals[k], major[k], minor[k]
        end
    end
    return j
end

"""
    compute_offsets(n, inds, len=length(inds)) -> offs

yields the `n+1` offsets computed from the list of indices `inds[1:len]`.
Indices must be in non-increasing order an all in the range `1:n`.

"""
function compute_offsets(n::Int,
                         inds::AbstractVector{Int},
                         len::Int = length(inds))
    @certify len ≤ length(inds)
    @inbounds begin
        offs = Vector{Int}(undef, n + 1)
        i = 0
        for k in 1:len
            j = inds[k]
            j == i && continue
            ((j < i)|(j > n)) &&
                error(1 ≤ j ≤ n ?
                      "indices must be in non-increasing order" :
                      "out of bound indices")
            off = k - 1
            while i < j
                i += 1
                offs[i] = off
            end
        end
        while i ≤ n
            i += 1
            offs[i] = len
        end
    end
    return offs
end

# This error is due to the non-zeros selector not returning the same results in
# the two selection passes.
bad_selector() = throw_argument_error("inconsistent selector function")

@inline select_non_zeros(v::Bool, i::Int, j::Int) = v
@inline select_non_zeros(v::T, i::Int, j::Int) where {T} = (v != zero(T))

@inline select_non_zeros_in_diagonal(v::T, i::Int, j::Int) where {T} =
    ((i == j)|(v != zero(T)))

@inline select_non_zeros_in_lower_part(v::T, i::Int, j::Int) where {T} =
    ((i ≥ j)&(v != zero(T)))

@inline select_non_zeros_in_upper_part(v::T, i::Int, j::Int) where {T} =
    ((i ≤ j)&(v != zero(T)))

"""
    get_equivalent_size(A, Val(M), Val(N)) -> nrows, ncols, rowsiz, colsiz

yields equivalent matrix dimensions of array `A` assuming the *rows* account
for the `M` leading dimensions while the *columns* account for the other `N`
dimensions.

"""
function get_equivalent_size(A::AbstractArray{T,L},
                             ::Val{M}, ::Val{N}) where {T,L,M,N}
    @certify L == M + N
    @certify M ≥ 1
    @certify N ≥ 1
    eachindex(A) == 1:length(A) ||
        throw_argument_error("array must have standard linear indexing")
    siz = size(A)
    rowsiz = siz[1:M]
    colsiz = siz[M+1:end]
    nrows = prod(rowsiz)
    ncols = prod(colsiz)
    return nrows, ncols, rowsiz, colsiz
end

"""
    count_selection(A, sel, rowmajor=false) -> nvals

yields the number of selected entries in matrix `A` such that `sel(A[i,j],i,j)`
is `true` and with `i` and `j` the row and column indices.  if optinal argument
`rowmajor` is true, the array is walked in row-major order; otherwise (the
default), the array is walked in column-major order.

"""
function count_selection(A::AbstractMatrix, sel::Function,
                         rowmajor::Bool = false)
    nrows, ncols = size(A)
    nvals = 0
    if rowmajor
        # Walk the coefficients in row-major order.
        @inbounds for i in 1:nrows, j in 1:ncols
            if sel(A[i,j], i, j)
                nvals += 1
            end
        end
    else
        # Walk the coefficients in column-major order.
        @inbounds for j in 1:ncols, i in 1:nrows
            if sel(A[i,j], i, j)
                nvals += 1
            end
        end
    end
    return nvals
end

"""
    check_structure(A) -> A

checks the structure of the compressed sparse operator `A` throwing an
exception if there are any inconsistencies.

"""
function check_structure(A::CompressedSparseOperator{:CSR})
    check_size(A)
    check_vals(A)
    check_cols(A)
    check_offs(A)
    return A
end

function check_structure(A::CompressedSparseOperator{:CSC})
    check_size(A)
    check_vals(A)
    check_rows(A)
    check_offs(A)
    return A
end

function check_structure(A::CompressedSparseOperator{:COO})
    check_size(A)
    check_vals(A)
    check_rows(A)
    check_cols(A)
    return A
end

"""
    check_size(siz, id="array") -> len

checks the array size `siz` and returns the corresponding number of elements.
An `ArgumentError` is thrown if a dimension is invalid, using `id` to identify
the argument.

    check_size(A)

checks the validity of the row and column sizes in compressed sparse operator
`A` throwing an exception if there are any inconsistencies.

"""
function check_size(siz::NTuple{N,Int}, id::String="array") where {N}
    len = 1
    @inbounds for i in 1:N
        (dim = siz[i]) ≥ 0 || bad_dimension(dim, i, id)
        len *= dim
    end
    return len
end

function check_size(A::SparseOperator)
    check_size(row_size(A), "row") == nrows(A) ||
        throw_dimension_mismatch("incompatible equivalent number of rows and row size")
    check_size(col_size(A), "column") == ncols(A) ||
        throw_dimension_mismatch("incompatible equivalent number of columns and column size")
    nothing
end

@noinline bad_dimension(dim::Integer, i::Integer, id) =
    throw_argument_error("invalid ", i, ordinal_suffix(i), " ", id,
                         " dimension: ", dim)

"""
    check_vals(A)

checks the array of values in compressed sparse operator `A` throwing an
exception if there are any inconsistencies.

"""
function check_vals(A::SparseOperator)
    vals = get_vals(A)
    is_fast_array(vals) || throw_not_fast_array("array of values")
    length(vals) == nnz(A) || throw_argument_error("bad number of values")
    nothing
end

"""
    check_rows(A)

checks the array of linear row indices in the compressed sparse operator `A`
stored in a *Compressed Sparse Column* (CSC) or *Compressed Sparse Coordinate*
(COO) format.  Throws an exception in case of inconsistency.

    check_rows(rows, m)

check the array of linear row indices `rows` for being a fast vetor of values
in the range `1:m`.

"""
function check_rows(A::Union{<:CompressedSparseOperator{:CSC},
                             <:CompressedSparseOperator{:COO}})
    rows = get_rows(A)
    length(rows) == nnz(A) ||
        throw_argument_error("bad number of row indices")
    check_rows(rows, nrows(A))
    # FIXME: also check sorting for CompressedSparseOperator{:CSC}?
end

function check_rows(rows::AbstractVector{Int}, m::Int)
    is_fast_array(rows) || throw_not_fast_array("array of row indices")
    anyerror = false
    @inbounds @simd for k in eachindex(rows)
        i = rows[k]
        anyerror |= ((i < 1)|(i > m))
    end
    anyerror && error("out of range row indices")
    nothing
end

"""
    check_cols(A)

checks the array of linear column indices in the compressed sparse operator `A`
stored in a *Compressed Sparse Row* (CSR) or *Compressed Sparse Coordinate*
(COO) format.  Throws an exception in case of inconsistency.

    check_cols(cols, n)

check the array of linear column indices `cols` for being a fast vetor of
values in the range `1:n`.

"""
function check_cols(A::Union{<:CompressedSparseOperator{:CSR},
                             <:CompressedSparseOperator{:COO}})
    cols = get_cols(A)
    length(cols) == nnz(A) ||
        throw_argument_error("bad number of column indices")
    check_cols(cols, ncols(A))
    # FIXME: also check sorting for CompressedSparseOperator{:CSR}?
end

function check_cols(cols::AbstractVector{Int}, n::Int)
    is_fast_array(cols) || throw_not_fast_array("array of column indices")
    anyerror = false
    @inbounds @simd for k in eachindex(cols)
        j = cols[k]
        anyerror |= ((j < 1)|(j > n))
    end
    anyerror && error("out of range column indices")
    nothing
end

"""
    check_offs(A)

checks the array of offsets in the compressed sparse operator `A` stored in a
*Compressed Sparse Row* (CSR) or *Compressed Sparse Column* (CSC) format.
Throws an exception in case of inconsistency.

"""
function check_offs(A::T) where {T<:Union{CompressedSparseOperator{:CSR},
                                          CompressedSparseOperator{:CSC}}}
    offs = get_offs(A)
    is_fast_array(offs) || throw_not_fast_array("array of offsets")
    n = (T <: CompressedSparseOperator{:CSR} ? nrows(A) : ncols(A))
    length(offs) == n + 1 || error("bad number of offsets")
    offs[1] == 0 || error("bad initial offset")
    len = 0
    anyerrors = false
    @inbounds for i in 1:n
        k1, k2 = offs[i], offs[i+1]
        anyerrors |= (k2 < k1)
        len += ifelse(k2 > k1, k2 - k1, 0)
    end
    anyerrors && error("offsets must be non-decreasing")
    len == nnz(A) ||
        error("offsets incompatible with number of structural non-zeros")
    nothing
end

"""
    unsafe_csr([m, n,] vals, cols, offs, rowsiz, colsiz)

yields a compressed sparse operator in *Compressed Sparse Row* (CSR) format as
an instance of `SparseOperatorCSR`.  This method assumes that arguments are
correct, it just calls the inner constructor with suitable parameters.  This
method is mostly used by converters and outer constructors.

"""
function unsafe_csr(m::Integer, n::Integer,
                    vals::V, cols::J, offs::K,
                    rowsiz::NTuple{M,Int},
                    colsiz::NTuple{N,Int}) where {T,M,N,
                                                  V<:AbstractVector{T},
                                                  J<:AbstractVector{Int},
                                                  K<:AbstractVector{Int}}
    SparseOperatorCSR{T,M,N,V,J,K}(to_int(m), to_int(n),
                                   vals, cols, offs,
                                   rowsiz, colsiz)
end

function unsafe_csr(vals::AbstractVector,
                    cols::AbstractVector{Int},
                    offs::AbstractVector{Int},
                    rowsiz::NTuple{M,Int},
                    colsiz::NTuple{N,Int}) where {M,N}
    unsafe_csr(prod(rowsiz), prod(colsiz), vals, cols, offs, rowsiz, colsiz)
end

"""
    unsafe_csc([m, n,] vals, rows, offs, rowsiz, colsiz)

yields a compressed sparse operator in *Compressed Sparse Column* (CSC) format
as an instance of `SparseOperatorCSC`.  This method assumes that arguments are
correct, it just calls the inner constructor with suitable parameters.  This
method is mostly used by converters and outer constructors.

"""
function unsafe_csc(m::Integer, n::Integer,
                    vals::V, rows::I, offs::K,
                    rowsiz::NTuple{M,Int},
                    colsiz::NTuple{N,Int}) where {T,M,N,
                                                  V<:AbstractVector{T},
                                                  I<:AbstractVector{Int},
                                                  K<:AbstractVector{Int}}
    SparseOperatorCSC{T,M,N,V,I,K}(to_int(m), to_int(n),
                                   vals, rows, offs,
                                   rowsiz, colsiz)
end

function unsafe_csc(vals::AbstractVector,
                    rows::AbstractVector{Int},
                    offs::AbstractVector{Int},
                    rowsiz::NTuple{M,Int},
                    colsiz::NTuple{N,Int}) where {M,N}
    unsafe_csc(prod(rowsiz), prod(colsiz), vals, rows, offs, rowsiz, colsiz)
end

"""
    unsafe_coo([m, n,] vals, rows, cols, rowsiz, colsiz)

yields a compressed sparse operator in *Compressed Sparse Coordinate* (COO)
format as an instance of `SparseOperatorCOO`.  This method assumes that
arguments are correct, it just calls the inner constructor with suitable
parameters.  This method is mostly used by converters and outer constructors.

"""
function unsafe_coo(m::Integer, n::Integer,
                    vals::V, rows::I, cols::J,
                    rowsiz::NTuple{M,Int},
                    colsiz::NTuple{N,Int}) where {T,M,N,
                                                  V<:AbstractVector{T},
                                                  I<:AbstractVector{Int},
                                                  J<:AbstractVector{Int}}
    SparseOperatorCOO{T,M,N,V,I,J}(to_int(m), to_int(n),
                                   vals, rows, cols,
                                   rowsiz, colsiz)
end

function unsafe_coo(vals::AbstractVector,
                    rows::AbstractVector{Int},
                    cols::AbstractVector{Int},
                    rowsiz::NTuple{M,Int},
                    colsiz::NTuple{N,Int}) where {M,N}
    unsafe_coo(prod(rowsiz), prod(colsiz), vals, rows, cols, rowsiz, colsiz)
end

"""
    is_fast_array(A) -> bool

yields whether `A` is a *fast array* that is an array with standard linear
indexing.

"""
is_fast_array(A::AbstractArray) = is_fast_indices(eachindex(A))

"""
    is_fast_indices(inds) -> bool

yields whether `inds` is an iterator for *fast indices* that is linear indices
starting at `1`.

"""
is_fast_indices(inds::AbstractUnitRange{Int}) = (first(inds) == 1)
is_fast_indices(inds) = false

"""
    check_argument(A, siz, id="array")

checks whether array `A` has size `siz` and implements standard linear
indexing.  An exception is thrown if any of these do not hold.

"""
function check_argument(A::AbstractArray{<:Any,N},
                        siz::NTuple{N,Int},
                        id="array") where {N}
    IndexStyle(A) === IndexLinear() || throw_non_linear_indexing(id)
    inds = axes(A)
    @inbounds for i in 1:N
        first(inds[i]) == 1 || throw_non_standard_indexing(id)
        length(inds[i]) == siz[i] || throw_incompatible_dimensions(id)
    end
    nothing
end

"""
    throw_argument_error(args...)

throws an `ArgumentError` exception with a textual message made of `args...`.

"""
throw_argument_error(mesg::AbstractString) = throw(ArgumentError(mesg))
@noinline throw_argument_error(args...) = throw_argument_error(string(args...))

@noinline throw_not_fast_array(id) =
    throw_argument_error(id, " does not implement fast indexing")

@noinline throw_non_linear_indexing(id) =
    throw_argument_error(id, " does not implement linear indexing")

@noinline throw_non_standard_indexing(id) =
    throw_argument_error(id, " has non-standard indexing")

"""
    throw_assertion_error(args...)

throws an `AssertionError` exception with a textual message made of `args...`.

"""
throw_assertion_error(mesg::AbstractString) = throw(AssertionError(mesg))
@noinline throw_assertion_error(args...) = throw_assertion_error(string(args...))

"""
    throw_dimension_mismatch(args...)

throws a `DimensionMismatch` exception with a textual message made of
`args...`.

"""
throw_dimension_mismatch(mesg::AbstractString) =
    throw(DimensionMismatch(mesg))

@noinline throw_dimension_mismatch(args...) =
    throw_dimension_mismatch(string(args...))

@noinline throw_incompatible_dimensions(id) =
    throw_dimension_mismatch(id, " has incompatible dimensions")

@noinline throw_incompatible_dimensions() =
    throw_dimension_mismatch("incompatible dimensions")

@noinline throw_incompatible_number_of_dimensions() =
    throw_dimension_mismatch("incompatible number of dimensions")

@noinline throw_incompatible_number_of_elements() =
    throw_dimension_mismatch("incompatible number of elements")

"""
    ordinal_suffix(n) -> "st" or "nd" or "rd" or "th"

yields the ordinal suffix for integer `n`.

"""
function ordinal_suffix(n::Integer)
    if n > 0
        d = mod(n, 10)
        if d == 1
            return "st"
        elseif d == 2
            return "nd"
        elseif d == 3
            return "rd"
        end
    end
    return "th"
end

#------------------------------------------------------------------------------
# Apply operators.

"""
    dispatch_multipliers!(α, f, A, x, β, y) -> y

dispatch calls to function `f` as `f(α,A,x,β,axpy)` with `α`, `A`, `x`, `β` and
`y` the other arguments and where `axpy` is a function called with 4 scalar
arguments as `axpy(α,x,β,y)` to yield `α*x + β*y` but which is optimized
depending on the values of the multipliers `α` and `β`.  For instance, if `α=1`
and `β=0`, then `axpy(α,x,β,y)` just evaluates as `x`.

The `dispatch_multipliers!` method is a helper to apply a mapping `A` to an
argument `x` and store the result in `y`.  In pseudo-code, this amounts to
performing `y <- α*op(A)(x) + β*y` where `op(A)` denotes a variant of `A` which
usually depends on `f`.

"""
@inline function dispatch_multipliers!(α::Number, f::Function, A, x,
                                       β::Number, y)
    if α == 0
        vscale!(y, β)
    elseif α == 1
        if β == 0
            f(1, A, x, 0, y, axpby_yields_x)
        elseif β == 1
            f(1, A, x, 1, y, axpby_yields_xpy)
        else
            b = promote_multiplier(β, eltype(y))
            f(1, A, x, b, y, axpby_yields_xpby)
        end
    else
        a = promote_multiplier(α, eltype(A), eltype(x))
        if β == 0
            f(a, A, x, 0, y, axpby_yields_ax)
        elseif β == 1
            f(a, A, x, 1, y, axpby_yields_axpy)
        else
            b = promote_multiplier(β, eltype(y))
            f(a, A, x, b, y, axpby_yields_axpby)
        end
    end
    return y
end

# Generic version of `vcreate` for most compressed sparse operators.
#
# We assume that in-place operation is not possible and thus simply ignore the
# `scratch` flag.  Operators which can be applied in-place shall specialize
# this method.  We do not check the dimensions and indexing of `x` as this will
# be done when `apply!` is called.

function vcreate(::Type{P},
                 A::SparseOperator{Ta,M,N},
                 x::AbstractArray{Tx,N},
                 scratch::Bool) where {Ta,Tx,M,N,P<:Union{Direct,InverseAdjoint}}
    Ty = promote_type(Ta,Tx)
    return Array{Ty}(undef, row_size(A))
end

function vcreate(::Type{P},
                 A::SparseOperator{Ta,M,N},
                 x::AbstractArray{Tx,M},
                 scratch::Bool) where {Ta,Tx,M,N,P<:Union{Adjoint,Inverse}}
    Ty = promote_type(Ta,Tx)
    return Array{Ty}(undef, col_size(A))
end

# Apply a sparse linear mapping, and its adjoint, stored in Compressed Sparse
# Row (CSR) format.

function apply!(α::Number,
                ::Type{Direct},
                A::CompressedSparseOperator{:CSR,Ta,M,N},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{Ty,M}) where {Ta,Tx,Ty,M,N}
    check_argument(x, col_size(A))
    check_argument(y, row_size(A))
    dispatch_multipliers!(α, unsafe_apply_direct!, A, x, β, y)
end

function unsafe_apply_direct!(α::Number,
                              A::CompressedSparseOperator{:CSR,Ta,M,N},
                              x::AbstractArray{Tx,N},
                              β::Number,
                              y::AbstractArray{Ty,M},
                              axpby::Function) where {Ta,Tx,Ty,M,N}
    @inbounds for i in each_row(A)
        s = zero(promote_type(Ta, Tx))
        for k in each_off(A, i)
            j = get_col(A, k)
            v = get_val(A, k)
            s += v*x[j]
        end
        y[i] = axpby(α, s, β, y[i])
    end
end

function apply!(α::Number,
                ::Type{Adjoint},
                A::CompressedSparseOperator{:CSR,Ta,M,N},
                x::AbstractArray{Tx,M},
                scratch::Bool,
                β::Number,
                y::AbstractArray{Ty,N}) where {Ta,Tx,Ty,M,N}
    check_argument(x, row_size(A))
    check_argument(y, col_size(A))
    Tm = promote_type(Ta, Tx) # to promote multipliers
    β == 1 || vscale!(y, β)
    if α == 1
        @inbounds for i in each_row(A)
            q = promote_multiplier(x[i], Tm)
            if q != 0
                for k in each_off(A, i)
                    j = get_col(A, k)
                    v = get_val(A, k)
                    y[j] += q*conj(v)
                end
            end
        end
    elseif α != 0
        a = promote_multiplier(α, Tm)
        @inbounds for i in each_row(A)
            q = a*promote_multiplier(x[i], Tm)
            if q != 0
                for k in each_off(A, i)
                    j = get_col(A, k)
                    v = get_val(A, k)
                    y[j] += q*conj(v)
                end
            end
        end
    end
    return y
end

# Apply a sparse operator, and its adjoint, stored in Compressed Sparse Column
# (CSC) format.

function apply!(α::Number,
                ::Type{Direct},
                A::CompressedSparseOperator{:CSC,Ta,M,N},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Real,
                y::AbstractArray{Ty,M}) where {Ta,Tx,Ty,M,N}
    check_argument(x, col_size(A))
    check_argument(y, row_size(A))
    Tm = promote_type(Ta, Tx) # to promote multipliers
    β == 1 || vscale!(y, β)
    if α == 1
        @inbounds for j in each_col(A)
            q = promote_multiplier(x[j], Tm)
            if q != 0
                for k in each_off(A, j)
                    i = get_row(A, k)
                    v = get_val(A, k)
                    y[i] += q*v
                end
            end
        end
    elseif α != 0
        a = promote_multiplier(α, Tm)
        @inbounds for j in each_col(A)
            q = a*promote_multiplier(x[j], Tm)
            if q != 0
                for k in each_off(A, j)
                    i = get_row(A, k)
                    v = get_val(A, k)
                    y[i] += q*v
                end
            end
        end
    end
    return y
end

function apply!(α::Real,
                ::Type{Adjoint},
                A::CompressedSparseOperator{:CSC,Ta,M,N},
                x::AbstractArray{Tx,M},
                scratch::Bool,
                β::Real,
                y::AbstractArray{Ty,N}) where {Ta,Tx,Ty,M,N}
    check_argument(x, row_size(A))
    check_argument(y, col_size(A))
    dispatch_multipliers!(α, unsafe_apply_adjoint!, A, x, β, y)
end

function unsafe_apply_adjoint!(α::Number,
                               A::CompressedSparseOperator{:CSC,Ta,M,N},
                               x::AbstractArray{Tx,M},
                               β::Number,
                               y::AbstractArray{Ty,N},
                               axpby::Function) where {Ta,Tx,Ty,M,N}
    @inbounds for j in each_col(A)
        s = zero(promote_type(Ta, Tx))
        for k in each_off(A, j)
            i = get_row(A, k)
            v = get_val(A, k)
            s += conj(v)*x[i]
        end
        y[j] = axpby(α, s, β, y[j])
    end
    return y
end

# Apply a sparse operator, and its adjoint, stored in Compressed Sparse
# Coordinate (COO) format.

function apply!(α::Number,
                ::Type{Direct},
                A::CompressedSparseOperator{:COO,Ta,M,N},
                x::AbstractArray{Tx,N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{Ty,M}) where {Ta,Tx,Ty,M,N}
    check_argument(x, col_size(A))
    check_argument(y, row_size(A))
    β == 1 || vscale!(y, β)
    if α != 0
        V, I, J = get_vals(A), get_rows(A), get_cols(A)
        if α == 1
            @inbounds for k in eachindex(V, I, J)
                v, i, j = V[k], I[k], J[k]
                y[i] += x[j]*v
            end
        elseif α == -1
            @inbounds for k in eachindex(V, I, J)
                v, i, j = V[k], I[k], J[k]
                y[i] -= x[j]*v
            end
        else
            # The ordering of operations is to minimize the number of
            # operations in case `v` is complex while `α` and `x` are reals.
            alpha = promote_multiplier(α, Ta, Tx)
            @inbounds for k in eachindex(V, I, J)
                v, i, j = V[k], I[k], J[k]
                y[i] += (alpha*x[j])*v
            end
        end
    end
    return y
end

function apply!(α::Number,
                ::Type{Adjoint},
                A::CompressedSparseOperator{:COO,Ta,M,N},
                x::AbstractArray{Tx,M},
                scratch::Bool,
                β::Number,
                y::AbstractArray{Ty,N}) where {Ta,Tx,Ty,M,N}
    check_argument(x, row_size(A))
    check_argument(y, col_size(A))
    β == 1 || vscale!(y, β)
    if α != 0
        V, I, J = get_vals(A), get_rows(A), get_cols(A)
        if α == 1
            @inbounds for k in eachindex(V, I, J)
                v, i, j = V[k], I[k], J[k]
                y[j] += x[i]*conj(v)
            end
        elseif α == -1
            @inbounds for k in eachindex(V, I, J)
                v, i, j = V[k], I[k], J[k]
                y[j] -= x[i]*conj(v)
            end
        else
            # The ordering of operations is to minimize the number of
            # operations in case `v` is complex while `α` and `x` are reals.
            alpha = promote_multiplier(α, Ta, Tx)
            @inbounds for k in eachindex(V, I, J)
                v, i, j = V[k], I[k], J[k]
                y[j] += (alpha*x[i])*conj(v)
            end
        end
    end
    return y
end

end # module SparseOperators

# The following module is to facilitate using compressed sparse operators at a
# lower level than the exported API.
module SparseMethods

export
    CompressedSparseOperator,
    SparseOperator,
    SparseOperatorCOO,
    SparseOperatorCSC,
    SparseOperatorCSR,
    col_size,
    copy_cols,
    copy_rows,
    copy_vals,
    each_col,
    each_off,
    each_row,
    get_col,
    get_cols,
    get_offs,
    get_row,
    get_rows,
    get_val,
    get_vals,
    ncols,
    nnz,
    nonzeros,
    nrows,
    row_size,
    set_val!

import ..SparseOperators:
    CompressedSparseOperator,
    SparseOperator,
    SparseOperatorCOO,
    SparseOperatorCSC,
    SparseOperatorCSR,
    col_size,
    copy_cols,
    copy_rows,
    copy_vals,
    each_col,
    each_off,
    each_row,
    get_col,
    get_cols,
    get_offs,
    get_row,
    get_rows,
    get_val,
    get_vals,
    ncols,
    nnz,
    nonzeros,
    nrows,
    row_size,
    set_val!

end # module SparseMethods


text:
# adaptive step
# with noise
using Distributions
@everywhere using DifferentialEquations
using PyPlot
using StatsBase
include("../src/SMCABC_RKHS.jl")
include("Models.jl")
cd(@__DIR__)

# data generation for test---------
@everywhere function generate_dataset()
    tspan = (0.0, 1.0)
    t = linspace(0,1,1001)
    t_mark = linspace(0,1,51)
    y0 = [1.0, 0.0, 0.0]
    p = [4.0, 2.0, 40.0, 20.0]

    alg = CVODE_BDF(max_hnil_warns = 0, max_error_test_failures = 0)
    prob = ODEProblem(testmodel, y0, tspan, p)
    res = solve(prob, alg, saveat=t, maxiters=1e7,
                abstol = 1e-8, reltol = 1e-8)
    res_mark = solve(prob, alg, saveat=t_mark, maxiters=1e7,
                     abstol = 1e-8, reltol = 1e-8)
    dataset = [vcat(res_mark.u'...)]
    return dataset
end
@everywhere dataset = generate_dataset()
nrow, ncol = size(dataset[1])
# dataset = [dataset[1] .+ rand(Normal(0,0.05), nrow, ncol)]
#--------------------------------
srand(20180503)

N = 5000
ϵ₀ = Inf
Δϵ_tol = 1e-8
α = 0.98
target = 0.03
Δϵ_window = 100
save_interval = 50
parallel = true
@everywhere γ = 1.
@everywhere nu = 1.
output = joinpath(@__DIR__, "result_parallel1")
@everywhere lbs = [1e-5, 1e-5, 1e-5, 1e-5]
@everywhere ubs = [1e+5, 1e+5, 1e+5, 1e+5]
@everywhere f(param) = obj_func(dataset, testmodel, param, lbs, ubs)

run_smc(N, lbs, ubs, f, ϵ₀, Δϵ_tol, α, Δϵ_window, output;
        save_interval=save_interval, target=target, parallel=parallel, γ=γ, nu=nu)

# plot -----------------------------------------------
# nrow=4
# ncol=1
# figsize=(5,10)
# fig, axes = subplots(nrow, ncol, figsize=figsize);
#
# for i in 1:4
#     axes[i][:hist](est_params[:,i], bins=100, normed=true, color="b");
#     axes[i][:axvline](p[i], color="r");
# end
# fig[:savefig](joinpath(output, "param_dist.png"), dpi=300);
# PyPlot.plt[:close]()
#
#
# using PyCall
# @pyimport numpy as np
# tspan = (0.0, 1.0)
# t = linspace(0,1,1001)
# y0 = [1.0, 0.0, 0.0]
# p = [4., 2., 40., 20.]
#
# trajectories = zeros(size(est_params)[1], length(t), 3)
# for i in 1:size(est_params)[1]
#     param = est_params[i,:]
#     prob = ODEProblem(testmodel, y0, tspan, param)
#     res_temp = solve(prob, alg, saveat=t, maxiters=1e7,
#                 abstol = 1e-8, reltol = 1e-8)
#     trajectories[i,:,:] = np.array(res_temp.u)
# end
#
# ave = mean(trajectories, Dims(1))
# sd = std(trajectories, Dims(1))
# y1=ave.-2sd
# y2=ave.+2sd
#
# prob = ODEProblem(testmodel, y0, tspan, p)
# res = solve(prob, alg, saveat=t, maxiters=1e7,
#             abstol = 1e-8, reltol = 1e-8)
#
# for i in 1:3
#     plot(t, vec(ave[:,:,i]), color="red", ls="--", lw=1);
#     fill_between(t, vec(y1[:,:,i]), vec(y2[:,:,i]));
# end
#
# plot(t, res[1, :], lw=1, color="k", ls="--");
# plot(t, res[2, :], lw=1, color="k", ls="--");
# plot(t, res[3, :], lw=1, color="k", ls="--");
#
# plot(t_mark, dataset[1][:, 1], lw=1, color="b", marker="o");
# plot(t_mark, dataset[1][:, 2], lw=1, color="orange", marker="o");
# plot(t_mark, dataset[1][:, 3], lw=1, color="g", marker="o");
#
# savefig(joinpath(output, "predict.png"), dpi=300);
# PyPlot.plt[:close]()

ext: jl


text:
Base.size(layer::Union{Binary, Spin, Potts, Gaussian, ReLU, pReLU, xReLU}) = size(layer.θ)
Base.length(layer::Union{Binary, Spin, Potts, Gaussian, ReLU, pReLU, xReLU}) = length(layer.θ)

const _FieldLayers = Union{Binary, Spin, Potts}

"""
    energies(layer, x)

Energies of units in layer (not reduced over layer dimensions).
"""
function energies(layer::_FieldLayers, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    return -layer.θ .* x
end

function energy(layer::_FieldLayers, x::AbstractArray)
    @assert size(layer) == size(x)[1:ndims(layer)]
    xconv = activations_convert_maybe(layer.θ, x)
    if ndims(layer) == ndims(x)
        return -dot(layer.θ, x)
    else
        Eflat = -vec(layer.θ)' * flatten(layer, xconv)
        return reshape(Eflat, batch_size(layer, x))
    end
end

function ∂free_energies(layer::_FieldLayers, inputs::Union{Real,AbstractArray} = 0)
    return (; θ = -mean_from_inputs(layer, inputs))
end

struct FieldStats{A<:AbstractArray}
    μ::A
    function FieldStats(layer::_FieldLayers, x::AbstractArray; wts=nothing)
        μ = batchmean(layer, x; wts)
        return new{typeof(μ)}(μ)
    end
end
Base.size(stats::FieldStats) = size(stats.μ)
Base.ndims(stats::FieldStats) = size(stats.μ)

suffstats(layer::_FieldLayers, data::AbstractArray; wts = nothing) = FieldStats(layer, data; wts)

function ∂energy(layer::_FieldLayers, stats::FieldStats)
    @assert size(layer) == size(stats.μ)
    return (; θ = -stats.μ)
end

"""
    colors(layer)

Number of possible states of units in discrete layers.
"""
colors(layer::Union{Spin,Binary}) = 2
colors(layer::Potts) = size(layer, 1)

"""
    sitedims(layer)

Number of dimensions of layer, with special handling of Potts layer,
for which the first dimension doesn't count as a site dimension.
"""
sitedims(layer::AbstractLayer) = ndims(layer)
sitedims(layer::Potts) = ndims(layer) - 1

"""
    sitesize(layer)

Size of layer, with special handling of Potts layer,
for which the first dimension doesn't count as a site dimension.
"""
sitesize(layer::AbstractLayer) = size(layer)
sitesize(layer::Potts) = size(layer)[2:end]

function pReLU(layer::dReLU)
    abs_γp = abs.(layer.γp)
    abs_γn = abs.(layer.γn)
    γ = @. 2abs_γp * abs_γn / (abs_γp + abs_γn)
    η = @. (abs_γn - abs_γp) / (abs_γp + abs_γn)
    θ = @. (layer.θp * abs_γn + layer.θn * abs_γp) / (abs_γp + abs_γn)
    Δ = @. γ * (layer.θp - layer.θn) / (abs_γp + abs_γn)
    return pReLU(θ, γ, Δ, η)
end

function dReLU(layer::pReLU)
    γp = @. layer.γ / (1 + layer.η)
    γn = @. layer.γ / (1 - layer.η)
    θp = @. layer.θ + layer.Δ / (1 + layer.η)
    θn = @. layer.θ - layer.Δ / (1 - layer.η)
    return dReLU(θp, θn, γp, γn)
end

function xReLU(layer::dReLU)
    abs_γp = abs.(layer.γp)
    abs_γn = abs.(layer.γn)
    γ = @. 2abs_γp * abs_γn / (abs_γp + abs_γn)
    ξ = @. (abs_γn - abs_γp) / (abs_γp + abs_γn - abs(abs_γn - abs_γp))
    θ = @. (layer.θp * abs_γn + layer.θn * abs_γp) / (abs_γp + abs_γn)
    Δ = @. γ * (layer.θp - layer.θn) / (abs_γp + abs_γn)
    return xReLU(θ, γ, Δ, ξ)
end

function dReLU(layer::xReLU)
    ξp = @. (1 + abs(layer.ξ)) / (1 + max(2layer.ξ, 0))
    ξn = @. (1 + abs(layer.ξ)) / (1 - min(2layer.ξ, 0))
    γp = @. layer.γ * ξp
    γn = @. layer.γ * ξn
    θp = @. layer.θ + layer.Δ * ξp
    θn = @. layer.θ - layer.Δ * ξn
    return dReLU(θp, θn, γp, γn)
end

function xReLU(layer::pReLU)
    ξ = @. layer.η / (1 - abs(layer.η))
    return xReLU(layer.θ, layer.γ, layer.Δ, ξ)
end

function pReLU(layer::xReLU)
    η = @. layer.ξ / (1 + abs(layer.ξ))
    return pReLU(layer.θ, layer.γ, layer.Δ, η)
end

dReLU(layer::Gaussian) = dReLU(layer.θ, layer.θ, layer.γ, layer.γ)
pReLU(layer::Gaussian) = pReLU(dReLU(layer))
xReLU(layer::Gaussian) = xReLU(dReLU(layer))

#dReLU(layer::ReLU) = dReLU(layer.θ, zero(layer.θ), layer.γ, inf.(layer.γ))

# function pReLU(layer::ReLU)
#     θ = layer.θ
#     γ = 2layer.γ
#     η = one.(layer.γ)
#     Δ = zero.(layer.θ)
#     return pReLU(θ, γ, Δ, η)
# end

# function xReLU(layer::ReLU)
#     θ = layer.θ
#     γ = 2layer.γ
#     ξ = inf.(layer.γ)
#     Δ = zero.(layer.θ)
#     return xReLU(θ, γ, Δ, ξ)
# end


text:
__precompile__()
module SeisNoise

using Dates, DataFrames, DSP, FFTW, Glob, JLD2, LinearAlgebra, SeisIO
using Statistics, StatsBase, Interpolations, GLM, Plots, Distributed, LightXML

# import types first
include("Types/RawData.jl")
include("Types/FFTData.jl")
include("Types/CorrData.jl")
include("Types/show.jl")
include("Types/InputParams.jl")

# import pre and post processing tools
include("distance.jl")
include("io.jl")
include("ArrayFuncs.jl")
include("tools.jl")
include("slicing.jl")
include("filter.jl")
include("downsample.jl")
include("availability.jl")
include("phase_shift.jl")
include("stacking.jl")

# import  routines for doin' stuff
include("compute_fft.jl")
include("correlation.jl")
include("rotation.jl")
include("VelocityChange/MWCS.jl")
include("VelocityChange/Stretching.jl")
include("VelocityChange/Wavelets.jl")
include("Plotting/plotting.jl")

end # module


text:
# set random seed to promote repeatability in CI unit tests
using Random
Random.seed!(101)

for filename in (
        "jop_convolve.jl",
        "jop_envelope.jl",
        "jop_filter.jl")
    include(filename)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct BotsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_botsInfo(_api::BotsApi; bot=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/bots.info", ["slackAuth"])
    Swagger.set_param(_ctx.query, "bot", bot)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets information about a bot user.
Param: bot::String
Return: Dict{String, Any}
"""
function botsInfo(_api::BotsApi; bot=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_botsInfo(_api; bot=bot, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function botsInfo(_api::BotsApi, response_stream::Channel; bot=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_botsInfo(_api; bot=bot, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export botsInfo