Optimize transfer operator computation

[rusandris]

Initial illustrative PR (draft):

• replace visitors transition counting with a simple loop over the symbolic trajectory in transferoperator
• remove visitors from TransferOperatorApproximationRectangular
• remove overwriting L, use N to store number of symbols (alphabet size)
• lacks boundary conditions, stochasticity checking etc.

[rusandris]

Latest changes:

• add boundary conditions again, default is :none
• remove sort_idxs from TransferOperatorApproximationRectangular
• separate helper functions to utils.jl
• add logistic map invariant measure test
• rename some variables to hint towards more general usage, add comments

[Datseris]

I am on holiday now, can review seriously at end of September. In the meantime would be nice if you can update this with main branch so that tests can run.

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 72.85714% with 19 lines in your changes missing coverage. Please review.

Project coverage is 89.65%. Comparing base (9fcf9d7) to head (13a856f). Report is 6 commits behind head on main.

Files with missing lines	Patch %	Lines
src/outcome_spaces/transfer_operator/utils.jl	69.64%	17 Missing :warning:
...come_spaces/transfer_operator/transfer_operator.jl	85.71%	2 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #423      +/-   ##
==========================================
- Coverage   91.97%   89.65%   -2.33%     
==========================================
  Files          86       87       +1     
  Lines        2468     2475       +7     
==========================================
- Hits         2270     2219      -51     
- Misses        198      256      +58     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[rusandris]

Can you please post here the speed up versus the previous version?

Using the same benchmarks as before:

using DynamicalSystems
using BenchmarkTools
using ComplexityMeasures

#Henon system orbit
henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
henon = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
orbit, t = trajectory(henon, 10^7; Ttr = 10^4)
grid_size = 10 #10x10 grid binning

@btime begin
binning = RectangularBinning(grid_size)
to = ComplexityMeasures.transferoperator(orbit, binning);
P_cm = to.transfermatrix
end

Results:

#--------previous version---------------
1.749 s (1540 allocations: 828.16 MiB)
#--------new version---------------
777.273 ms (238 allocations: 152.62 MiB)

[Datseris]

@kahaaga this goes in, right?

[Datseris]

you mentioned you used some internal code somewhere, which i guess you can obtain by checking out an explicit commit?

[rusandris]

@kahaaga this goes in, right?

Any ideas on this PR? @kahaaga

[kahaaga]

I don't expect it to be breaking. My idea here was to avoid introducing an extra transition that's not present in the symbolic time series (hence changed the default from :circular to :none). But this definitely need to be discussed, and tests need to be added. Plus, if the symbolic time series is long enough, the differences in results are expected to be negligible

The circular boundary condition is a remnant from v1 of the package, where we had a simplex-based volume outcome space specifically dedicated to estimation from very short time series. In this case, not handling boundary conditions would be potentially detrimental to the estimation. For this outcome space, we'd approximate the transfer operator by volume overlaps between simplices, also using piecewise one-time-step linear projection of the simplex vertices. There would be an edge case where the projection of a simplex could lead to the projected simplex having a vertex outside the convex hull of the other simplices in the outcome space's volume. That would lead to a transition matrix where row sums didn't sum to 1, which would not be meaningful. The circular boundary condition remedied this, at the expense of some (assumed) minor bias.

However, the outcome space hasn't yet made it back in the package since we redesigned it. Therefore, a boundary condition of :none would not be breaking, I think. We can leave this as is, and deal with it properly in the next PR that deals with generalizing the transfer operator computation to all the outcome spaces.

So I've added a test for the logistic map (r=4) where the invariant measure is known analytically ρ(x) = 1/(π√x(1-x)). The test compares the probability of each bin using transferoperator and invariantmeasure with the probabilities obtained by integrating the analytical invariant measure over the same bins (10 bins were used). Intregrals is used for the numerical integration. It doesn't appear anywhere in the source code, only in the unit test (notice that the dependency is only added to test/Project.toml ).

Excellent! Thanks for adding this analytical example. The dependency is necessary for the tests, so keeping it in the test dependencies is fine.

Same deal, I think those comments were left by @kahaaga . My guess is it's related to how we remap the original symbols to unique indices (so a symbolic time series like [12,34,15] would become [1,2,3] internally) and use those to index into the matrix. Maybe we can discuss those at some later stage

We can deal with removing all the comments from the source code in the next PR which will probably introduce a more or less complete redesign of the interface for estimating the transfer operator. This PR is about improving speed of computation, not cleaning up the source code.

Using the same benchmarks as before:

The benchmarks are looking great. Thanks!

Otherwise, looks good to me. Before this can be merged, I'm just echoing @Datseris

• please export the transferoperator function in the docs.

[rusandris]

We can leave this as is, and deal with it properly in the next PR that deals with generalizing the transfer operator computation to all the outcome spaces.

Great, exactly my thoughts. Thanks!

please export the transferoperator function in the docs.

I've put it below the TransferOperator stuff, let me know if that is what you meant.

[rusandris]

The type it returns is not part of the public API.

You're right: only TransferOperator is exported, TransferOperatorApproximationRectangular isn't, neither was transferoperator, that's why I had to use it as ComplexityMeasures.transferoperator before.

I can export TransferOperatorApproximationRectangular as well for now, but I guess these types will have to be merged/redesigned when we move on to generalizing the transfer operator computation to all outcome spaces.

And I agree, it is a bit confusing to have all these around: TransferOperator, transferoperator, transfermatrix,TransferOperatorApproximationRectangular

[Datseris]

I am merging; it is best to not document the return type now since it will change in the near future.

[Datseris]

Merging and tagging this PR broke Associations.jl: it is no longer installable because it was using GroupSlices.

[rusandris]

Merging and tagging this PR broke Associations.jl: it is no longer installable because it was using GroupSlices.

I see. This is because src/outcome_spaces/transfer_operator/GroupSlices.jl is no longer included anywhere, and it had some exported stuff. It is tricky because the whole thing is a direct copy if an unregistered julia package: https://github.com/mcabbott/GroupSlices.jl

Maybe we could do one of the following things:

1. quickfix: include src/outcome_spaces/transfer_operator/GroupSlices.jl in ComplexityMeasures, even though it isn't used there
2. remove GroupSlices.jl from ComplexityMeasures, insert in Associations.jl
3. make registered package for GroupSlices, then it can be added to any package as a dependency

[Datseris]

Option 2 is the best; directly copy the code it needs there and use it there.

[Datseris]

but kristian probably knows best what is used and where an why