Generalize TransferOperator

src/outcome_spaces/transfer_operator/transfer_operator.jl

@@ -97,22 +97,8 @@ below some threshold.
 See also: [`RectangularBinning`](@ref), [`FixedRectangularBinning`](@ref),
 [`invariantmeasure`](@ref).
 """
-struct TransferOperator{R<:AbstractBinning, RNG} <: OutcomeSpace
-    binning::R
-    warn_precise::Bool
-    rng::RNG
-    function TransferOperator(b::R; 
-            rng::RNG = Random.default_rng(), 
-            warn_precise = true) where {R <: AbstractBinning, RNG}
-        return new{R, RNG}(b, warn_precise, rng)
-    end
-end
-function TransferOperator(ϵ::Union{Real,Vector}; 
-        rng = Random.default_rng(), warn_precise = true)
-    return TransferOperator(RectangularBinning(ϵ), warn_precise, rng)
-end
-
 
+struct TransferOperator <: ProbabilitiesEstimator end
 
 """
     TransferOperatorApproximationRectangular(to, binning::RectangularBinning, mini,
@@ -135,13 +121,11 @@ points that visits `bins[i]`.
 
 See also: [`RectangularBinning`](@ref).
 """
-struct TransferOperatorApproximationRectangular{
-        T<:Real,
-        BINS,
-        E}
+#should it include an encoding as well?
+struct TransferOperatorApproximation{T<:Real,OC<:OutcomeSpace} <: ProbabilitiesEstimator
     transfermatrix::AbstractArray{T, 2}
-    encoding::E
-    bins::BINS
+    outcome_space::OC
+    outcomes
 end
 
 """
@@ -151,34 +135,29 @@ Approximate the transfer operator given a set of sequentially ordered points sub
 rectangular partition given by the `binning`.
 The keywords `boundary_condition = :none, warn_precise = true` are as in [`TransferOperator`](@ref).
 """
-function transferoperator(pts::AbstractStateSpaceSet{D, T},
-        binning::Union{FixedRectangularBinning, RectangularBinning};
+function transferoperator(o::OutcomeSpace,x;
         boundary_condition = :none, 
-        warn_precise = true) where {D, T<:Real}
-
-    L = length(pts)
-    if warn_precise && !binning.precise
+        warn_precise = true)     
+    
+    #warning (only when used with some kind of binning)
+    if warn_precise && typeof(o) <: ValueBinning  && !o.binning.precise
         @warn "`binning.precise == false`. You may be getting points outside the binning."
     end
-    encoding = RectangularBinEncoding(binning, pts)
 
-    # The L points visits a total of N bins, which are the following bins (referenced
-    # here as cartesian coordinates, not absolute bins):
-    outcomes = map(pᵢ -> encode(encoding, pᵢ), pts)
-    #sort_idxs = sortperm(visited_bins)
-    #sort!(visited_bins) # see todo on github
+    outcomes = codify(o,x)
+    L = length(outcomes)
 
-    # There are N=length(unique(visited_bins)) unique bins.
-    # Which of the unqiue bins does each of the L points visit?
-    visits_whichbin,unique_outcomes = inds_in_terms_of_unique(outcomes, false) # set to true when sorting is fixed
+    # There are L number of outcomes
+    # turn the time series of outcomes into a sequence of unique indices of outcomes
+    unique_indices,unique_outcomes = inds_in_terms_of_unique(outcomes, false) # set to true when sorting is fixed
     N = length(unique_outcomes)
 
     #apply boundary conditions (default is :none)
     if boundary_condition == :circular
-        append!(visits_whichbin, [1])
+        append!(unique_indices, [1])
         L += 1
     elseif boundary_condition == :random
-        append!(visits_whichbin, [rand(rng, 1:length(visits_whichbin))])
+        append!(unique_indices, [rand(rng, 1:length(unique_indices))])
         L += 1
     elseif boundary_condition != :none
         error("Boundary condition $(boundary_condition) not implemented")
@@ -189,16 +168,15 @@ function transferoperator(pts::AbstractStateSpaceSet{D, T},
 
 	#count transitions in Q, assuming symbols from 1 to N
 	for i in 1:(L - 1)
-        Q[visits_whichbin[i],visits_whichbin[i+1]] += 1.0
+        Q[unique_indices[i],unique_indices[i+1]] += 1.0
 	end
 
     #normalize Q (not strictly necessary) and fill P by normalizing rows of Q
     Q .= Q./sum(Q)
-    P = calculate_transition_matrix(Q)
+    P = normalize_transition_matrix(Q)
 
-    unique!(outcomes)
-    return TransferOperatorApproximationRectangular(
-        P, encoding, outcomes)
+    return TransferOperatorApproximation(
+        P, o, unique_outcomes)
 end
 
 """
@@ -216,7 +194,7 @@ struct InvariantMeasure{T}
 end
 
 function invariantmeasure(iv::InvariantMeasure)
-    return iv.ρ, iv.to.bins
+    return iv.ρ, iv.to.outcomes
 end
 
 
@@ -273,7 +251,7 @@ The element `ρ[i]` is the probability of visitation to the box `bins[i]`.
 
 See also: [`InvariantMeasure`](@ref).
 """
-function invariantmeasure(to::TransferOperatorApproximationRectangular;
+function invariantmeasure(to::TransferOperatorApproximation;
         N::Int = 200, tolerance::Float64 = 1e-8, delta::Float64 = 1e-8,
         rng = Random.default_rng())
 
@@ -284,7 +262,7 @@ function invariantmeasure(to::TransferOperatorApproximationRectangular;
     =#
     Ρ = rand(rng, Float64, 1, size(to.transfermatrix, 1))
     Ρ = Ρ ./ sum(Ρ, dims = 2)
-
+    
     #=
     # Start estimating the invariant distribution. We could either do this by
     # finding the left-eigenvector of M, or by repeated application of M on Ρ
@@ -293,7 +271,6 @@ function invariantmeasure(to::TransferOperatorApproximationRectangular;
     # iterations.
     =#
     distribution = Ρ * to.transfermatrix
-
     distance = norm(distribution - Ρ) / norm(Ρ)
 
     check = floor(Int, 1 / delta)
@@ -319,7 +296,6 @@ function invariantmeasure(to::TransferOperatorApproximationRectangular;
                 distribution = distribution ./ colsum_distribution
             end
         end
-
         distance = norm(distribution - Ρ) / norm(Ρ)
     end
     distribution = dropdims(distribution, dims = 1)
@@ -357,25 +333,19 @@ end
 
 # Explicitly extend `probabilities` because we can skip the decoding step, which is 
 # expensive.
-function probabilities(est::TransferOperator, x::Array_or_SSSet)
-    to = transferoperator(StateSpaceSet(x), est.binning; 
-        warn_precise = est.warn_precise)
-    return Probabilities(invariantmeasure(to; rng = est.rng).ρ)
+function probabilities(pest::TransferOperator, o::OutcomeSpace, x::Array_or_SSSet;kwargs...)
+    to = transferoperator(o,x)
+    return Probabilities(invariantmeasure(to; kwargs...).ρ)
 end
 
-function probabilities_and_outcomes(est::TransferOperator, x::Array_or_SSSet)
-    to = transferoperator(StateSpaceSet(x), est.binning; 
-        warn_precise = est.warn_precise)
-    probs = invariantmeasure(to; rng = est.rng).ρ
-
-    # Note: bins are *not* sorted. They occur in the order of first appearance, according
-    # to the input time series. Taking the unique bins preserves the order of first
-    # appearance.
-    bins = to.bins
-    unique!(bins)
-    outs = decode.(Ref(to.encoding), bins) # coordinates of the visited bins
+function probabilities_and_outcomes(pest::TransferOperator, o::OutcomeSpace, x::Array_or_SSSet)
+    to = transferoperator(o,x)
+    probs = probabilities(pest, o, x)
+
+    #doesn't work for ValueBinning outcome space
+    outs = decode.(Ref(o.encoding), to.outcomes) # coordinates of the visited bins
     probs = Probabilities(probs, (outs, ))
-    return probs, outcomes(probs)
+    return probs, to.outcomes
 end
 
 outcome_space(est::TransferOperator, x) = outcome_space(est.binning, x)

src/outcome_spaces/transfer_operator/utils.jl

@@ -51,15 +51,15 @@ end
 inds_in_terms_of_unique(x::AbstractStateSpaceSet) = inds_in_terms_of_unique(x.data)
 
 
-function calculate_transition_matrix(S::SparseMatrixCSC;verbose=true)
+function normalize_transition_matrix(S::SparseMatrixCSC;verbose=true)
 	S_returned = deepcopy(S)
-	calculate_transition_matrix!(S_returned,verbose=verbose)
+	normalize_transition_matrix!(S_returned,verbose=verbose)
 	return S_returned
 end
 
 #normalize each row of S (sum is 1) to get p_ij trans. probabilities
 #by looping through CSC sparse matrix efficiently
-function calculate_transition_matrix!(S::SparseMatrixCSC;verbose=true)
+function normalize_transition_matrix!(S::SparseMatrixCSC;verbose=true)
 
     stochasticity = true
 

test/outcome_spaces/implementations/transfer_operator.jl

@@ -2,6 +2,39 @@ using ComplexityMeasures, Test
 using Random: MersenneTwister
 using DynamicalSystemsBase
 using Integrals
+
+#test for all count-based outcome spaces
+
+@testset begin "Count-based outcome spaces" 
+
+    #simple 1d rand 
+    x = rand(100)
+    x_ue = rand([0.0:0.1:1.0;],100) #for unique elements
+
+    #def count-based 1d outcome spaces 
+    outcome_spaces = [BubbleSortSwaps(),AmplitudeAwareOrdinalPatterns(),OrdinalPatterns(),
+        WeightedOrdinalPatterns(),CosineSimilarityBinning(),Dispersion(),
+        SequentialPairDistances(x),UniqueElements(),ValueBinning(RectangularBinning(5))] 
+
+    #build transferoperator from every outcomespace
+    for ocs in outcome_spaces
+        if ocs isa UniqueElements
+            to = transferoperator(ocs,x_ue) #unique elements 
+        else 
+            to = transferoperator(ocs,x)
+        end
+
+        #test if transition matrix is normalized
+        sum_rows = sum(to.transfermatrix;dims=2) 
+        @test all( isapprox.(1.0, sum_rows ;atol = 1e-3))
+    end
+
+    #leave out spatial methods for now:
+    #SpatialBubbleSortSwaps,SpatialDispersion,SpatialOrdinalPatterns
+    #not trivial how to implement transferoperator!
+
+end
+
 @test TransferOperator(RectangularBinning(3)) isa TransferOperator
 
 D = StateSpaceSet(rand(MersenneTwister(1234), 100, 2))
@@ -37,16 +70,17 @@ binnings = [
     orbit,t = trajectory(ds, 10^7; Ttr = 10^4)
 
     #--------------estimate invariant measure----------
-    b = RectangularBinning(10, true)
-    to = ComplexityMeasures.transferoperator(orbit, b)
+    b = ValueBinning(FixedRectangularBinning(range(0,1;length=11),1,true))
+    to = transferoperator(b,orbit)
     iv = invariantmeasure(to)
     p,outcomes = invariantmeasure(iv)
+    @show outcomes
 
     #order from leftmost bin to rightmost bin  
     p_bins = p[sortperm(outcomes)]
 
     #----------compute probability for each bin analytically----------
-    bin_ranges = to.encoding.ranges[1]
+    bin_ranges = to.outcome_space.binning.ranges[1]
     ρ_bins = zeros(10)
 
     for i in 1:length(bin_ranges)-1
@@ -64,8 +98,8 @@ end
 
 @testset "Binning test $i" for i in eachindex(binnings)
     b = binnings[i]
-    to = ComplexityMeasures.transferoperator(D, b)
-    @test to isa ComplexityMeasures.TransferOperatorApproximationRectangular
+    to = ComplexityMeasures.transferoperator(ValueBinning(b),D)
+    @test to isa ComplexityMeasures.TransferOperatorApproximation
 
     iv = invariantmeasure(to)
     @test iv isa InvariantMeasure
@@ -74,12 +108,11 @@ end
     @test p isa Probabilities
     @test bins isa Vector{Int}
 
-    o = TransferOperator(binnings[i])
-    @test probabilities(o, D) isa Probabilities
-    @test probabilities_and_outcomes(o, D) isa Tuple{Probabilities, Vector{SVector{2, Float64}}}
+    @test probabilities(TransferOperator(), ValueBinning(b) , D) isa Probabilities
+    @test probabilities_and_outcomes(TransferOperator(), ValueBinning(b), D) isa Tuple{Probabilities, Vector{SVector{2, Float64}}}
 
     # Test that gives approximately same entropy as ValueBinning:
-    abs(information(TransferOperator(b), D) - information(ValueBinning(b), D) ) < 0.1 # or something like that
+    abs(information(Shannon(), p) - information(ValueBinning(b), D) ) < 0.1 # or something like that
 end
 
 # Warn if we're not using precise binnings.