Refactoring TransferOperator

[April-Hannah-Lena]

This PR brings a new version of the TransferOperator, based on the SparseMatrixCSC format. It is as fast as the previous version when calculating eigenvalues, and includes new features:

• Does not require copying the underlying data of the BoxSet, as was required (multiple times) before
• Ability to act on BoxFuns via multiplication
• "On the fly" extension of the support

The docstring explains hopefully how the TransferOperator works internally:

----------

TransferOperator(map::BoxMap, support::BoxSet)

Discretization of the Perron-Frobenius operator, or transfer operator. Implemented as a sparse matrix with the same linear indices as support, e.g. if

julia> B = BoxSet(partition, [3,10,30])
  Boxset over [...] partition

julia> T = TransferOperator(boxmap, B)
  TransferOperator over [...] BoxSet

for some partition and boxmap, then

julia> axes(T)
  ([3, 10, 30], [3, 10, 30])

Fields:

• boxmap: SampledBoxMap map which relates to the transfers.
• support: BoxSet which contains keys for the already calculated transfers. Effectively, these are row/column pointers, i.e. the first column of T.mat contains transfer weights FROM box B_1, where B_1 is the first box of support.
• variant_set: BoxSet which contains keys for potential boxes lying outside of support, i.e. it could be that support is not an invariant set. In this case, variant_set contains the boxes which were not in support, but whose preimage lies in support.
• mat: SparseMatrixCSC containing transfer weights. The index T.mat[i,j] represents the transfer weight FROM the j'th box in support TO the i'th box in support. If support is not invariant, then the matrix will be tall. In this case the rows are counted in support and then in variant_set.

support -->
  |     .   .   .   .   .
  |     .   .   .   .   .
  v     .   .   .   .   .
        .   .  mat  .   .
        .   .   .   .   .
 var-   .   .   .   .   .
 ian-   .   .   .   .   .
 t_set  .   .   .   .   .     
  |     .   .   .   .   .
  |     .   .   .   .   .
  v     .   .   .   .   .

It is important to note that TranferOperator is only supported over the box set B, but if one lets a TranferOperator act on a BoxFun, then the support B is extended "on the fly" to include the support of the BoxFun.

Methods Implemented:

:(==), axes, size, eltype, getindex, setindex!, SparseArrays.sparse, Arpack.eigs, LinearAlgebra.mul! #, etc ...

Implementation detail:

The reader may have noticed that the matrix representation depends on the order of boxes in support. For this reason an OrderedSet is used. BoxSets using regular Sets will be copied and converted to OrderedSets.

----------------

The one consideration to be made is the code complexity for multiplication:

        function LinearAlgebra.mul!(y::BoxFun, g::$type, x::BoxFun)
            checkbounds(Bool, g, keys(x.vals)) || throw(BoundsError(g, x))
            zer = zero(eltype(y))
            map!(x -> zer, values(y.vals))
            rows = rowvals($gmap.mat)
            vals = nonzeros($gmap.mat)
            m, n = size($gmap)
            # less overall checks need to be made with the square version
            if m == n
                # iterate over columns with the keys that the columns represent
                for (col_j, j) in enumerate($gmap.support.set)
                    for k in nzrange($gmap.mat, col_j)
                        w = vals[k]
                        row_i = rows[k]
                        # grab the key that this row represents
                        i = $gmap.support.set.dict.keys[row_i]
                        y[$ind1] = @muladd y[$ind1] + $func(w) * x[$ind2]
                    end
                end
            else
                # iterate over columns with the keys that the columns represent
                for (col_j, j) in enumerate($gmap.support.set)
                    for k in nzrange($gmap.mat, col_j)
                        w = vals[k]
                        row_i = rows[k]
                        # grab the key that this row represents
                        if row_i > n
                            i = $gmap.variant_set.set.dict.keys[row_i - n]
                        else
                            i = $gmap.support.set.dict.keys[row_i]
                        end
                        y[$ind1] = @muladd y[$ind1] + $func(w) * x[$ind2]
                    end
                end
            end
            return y
        end

Multiplying the first eigenmeasure of the Lorenz attractor takes ~6ms on average.

using GAIO

# the Lorenz system
const σ, ρ, β = 10.0, 28.0, 0.4
v((x,y,z)) = (σ*(y-x), ρ*x-y-x*z, x*y-β*z)
f(x) = rk4_flow_map(v, x)

center, radius = (0,0,25), (30,30,30)
P = BoxPartition(Box(center, radius), (128,128,128))
F = BoxMap(f, P, no_of_points=200)

x = (sqrt(β*(ρ-1)), sqrt(β*(ρ-1)), ρ-1)         # equilibrium
W = unstable_set!(F, P[x])

T = TransferOperator(F, W)
(λ, ev) = eigs(T)

μ = ev[1]
@time T * μ

The code can be simplified, however:

        function LinearAlgebra.mul!(y::BoxFun, g::$type, x::BoxFun)
            checkbounds(Bool, g, keys(x.vals)) || throw(BoundsError(g, x))
            zer = zero(eltype(y))
            map!(x -> zer, values(y.vals))
            rows = rowvals($gmap.mat)
            vals = nonzeros($gmap.mat)
            m, n = size($gmap)
            # iterate over columns with the keys that the columns represent
            for (col_j, j) in enumerate($gmap.support.set)
                for k in nzrange($gmap.mat, col_j)
                    w = vals[k]
                    row_i = rows[k]
                    # grab the key that this row represents
                    if row_i > n
                        i = $gmap.variant_set.set.dict.keys[row_i - n]
                    else
                        i = $gmap.support.set.dict.keys[row_i]
                    end
                    y[$ind1] = @muladd y[$ind1] + $func(w) * x[$ind2]
                end
            end
            return y
        end

This version adds roughly one ms to the run time, i.e. multiplication takes ~7ms.

!!! A note on the chain recurrent set: currently it is broken because I am setting up boxgraph.jl. The code is already there, I first wrote it months ago. But I'm currently adapting it to the new TransferOperator version. This shouldn't take nearly as long as the TransferOperator update, since I don't really have to come up with anything new for this one :)

[April-Hannah-Lena]

boxgraph.jl is now ready. I can push it, or wait to merge this PR first so that there are less changes all at once.

[gaioguy]

Awesome, works as expected for me! I vote on favor of the shorter code ... Yes, let's merge this first.