From 97ce9d9a207502d8250f89c5dc9c350c8211f215 Mon Sep 17 00:00:00 2001 From: rusandris Date: Sat, 18 Jan 2025 20:06:12 +0200 Subject: [PATCH 1/4] first redesign mods --- .../transfer_operator/transfer_operator.jl | 44 +++++++++---------- 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/src/outcome_spaces/transfer_operator/transfer_operator.jl b/src/outcome_spaces/transfer_operator/transfer_operator.jl index bb131956..dd5f52dc 100644 --- a/src/outcome_spaces/transfer_operator/transfer_operator.jl +++ b/src/outcome_spaces/transfer_operator/transfer_operator.jl @@ -97,6 +97,7 @@ below some threshold. See also: [`RectangularBinning`](@ref), [`FixedRectangularBinning`](@ref), [`invariantmeasure`](@ref). """ +#= struct TransferOperator{R<:AbstractBinning, RNG} <: OutcomeSpace binning::R warn_precise::Bool @@ -112,7 +113,7 @@ function TransferOperator(ϵ::Union{Real,Vector}; return TransferOperator(RectangularBinning(ϵ), warn_precise, rng) end - +=# """ TransferOperatorApproximationRectangular(to, binning::RectangularBinning, mini, @@ -135,13 +136,13 @@ points that visits `bins[i]`. See also: [`RectangularBinning`](@ref). """ -struct TransferOperatorApproximationRectangular{ - T<:Real, - BINS, - E} + +struct TransferOperator <: ProbabilitiesEstimator end + +struct TransferOperatorApproximation{T<:Real,OC<:OutcomeSpace} <: ProbabilitiesEstimator transfermatrix::AbstractArray{T, 2} - encoding::E - bins::BINS + outcome_space::OC + outcomes end """ @@ -151,22 +152,18 @@ 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) + warn_precise = true) + #= if warn_precise && !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? @@ -197,8 +194,8 @@ function transferoperator(pts::AbstractStateSpaceSet{D, T}, P = calculate_transition_matrix(Q) unique!(outcomes) - return TransferOperatorApproximationRectangular( - P, encoding, outcomes) + return TransferOperatorApproximation( + P, o, outcomes) end """ @@ -273,7 +270,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()) @@ -357,10 +354,9 @@ 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) From 67dd672bd413bd56ec6da72738d20c00b18b9289 Mon Sep 17 00:00:00 2001 From: rusandris Date: Sun, 19 Jan 2025 20:24:40 +0200 Subject: [PATCH 2/4] some clear-ups, renamings --- .../transfer_operator/transfer_operator.jl | 29 ++++++++++--------- src/outcome_spaces/transfer_operator/utils.jl | 6 ++-- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/outcome_spaces/transfer_operator/transfer_operator.jl b/src/outcome_spaces/transfer_operator/transfer_operator.jl index dd5f52dc..1b857d9e 100644 --- a/src/outcome_spaces/transfer_operator/transfer_operator.jl +++ b/src/outcome_spaces/transfer_operator/transfer_operator.jl @@ -155,27 +155,26 @@ The keywords `boundary_condition = :none, warn_precise = true` are as in [`Trans function transferoperator(o::OutcomeSpace,x; boundary_condition = :none, warn_precise = true) - #= - if warn_precise && !binning.precise + + #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) - =# 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") @@ -186,16 +185,15 @@ function transferoperator(o::OutcomeSpace,x; #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 TransferOperatorApproximation( - P, o, outcomes) + P, o, unique_outcomes) end """ @@ -274,6 +272,10 @@ function invariantmeasure(to::TransferOperatorApproximation; N::Int = 200, tolerance::Float64 = 1e-8, delta::Float64 = 1e-8, rng = Random.default_rng()) + + @show N + @show tolerance + @show delta TO = to.transfermatrix #= # Start with a random distribution `Ρ` (big rho). Normalise it so that it @@ -320,6 +322,7 @@ function invariantmeasure(to::TransferOperatorApproximation; distance = norm(distribution - Ρ) / norm(Ρ) end distribution = dropdims(distribution, dims = 1) + @show counter # Do the last normalisation and check colsum_distribution = sum(distribution) diff --git a/src/outcome_spaces/transfer_operator/utils.jl b/src/outcome_spaces/transfer_operator/utils.jl index 0bc39707..887fe34f 100644 --- a/src/outcome_spaces/transfer_operator/utils.jl +++ b/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 From e0b0ae2cba01c55e54070883ec595de739ea9f7e Mon Sep 17 00:00:00 2001 From: rusandris Date: Mon, 20 Jan 2025 21:39:15 +0200 Subject: [PATCH 3/4] clear up some docs, comments --- .../transfer_operator/transfer_operator.jl | 24 +------------------ 1 file changed, 1 insertion(+), 23 deletions(-) diff --git a/src/outcome_spaces/transfer_operator/transfer_operator.jl b/src/outcome_spaces/transfer_operator/transfer_operator.jl index 1b857d9e..88ee72ef 100644 --- a/src/outcome_spaces/transfer_operator/transfer_operator.jl +++ b/src/outcome_spaces/transfer_operator/transfer_operator.jl @@ -97,23 +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, @@ -137,8 +122,6 @@ points that visits `bins[i]`. See also: [`RectangularBinning`](@ref). """ -struct TransferOperator <: ProbabilitiesEstimator end - struct TransferOperatorApproximation{T<:Real,OC<:OutcomeSpace} <: ProbabilitiesEstimator transfermatrix::AbstractArray{T, 2} outcome_space::OC @@ -272,10 +255,6 @@ function invariantmeasure(to::TransferOperatorApproximation; N::Int = 200, tolerance::Float64 = 1e-8, delta::Float64 = 1e-8, rng = Random.default_rng()) - - @show N - @show tolerance - @show delta TO = to.transfermatrix #= # Start with a random distribution `Ρ` (big rho). Normalise it so that it @@ -322,7 +301,6 @@ function invariantmeasure(to::TransferOperatorApproximation; distance = norm(distribution - Ρ) / norm(Ρ) end distribution = dropdims(distribution, dims = 1) - @show counter # Do the last normalisation and check colsum_distribution = sum(distribution) From e55411d82f97206a5a602987b18b4cb6557ee132 Mon Sep 17 00:00:00 2001 From: rusandris Date: Mon, 15 Sep 2025 16:03:26 +0300 Subject: [PATCH 4/4] revisit transferop and tests --- .../transfer_operator/transfer_operator.jl | 27 ++++------ .../implementations/transfer_operator.jl | 51 +++++++++++++++---- 2 files changed, 52 insertions(+), 26 deletions(-) diff --git a/src/outcome_spaces/transfer_operator/transfer_operator.jl b/src/outcome_spaces/transfer_operator/transfer_operator.jl index 88ee72ef..06fa6f1f 100644 --- a/src/outcome_spaces/transfer_operator/transfer_operator.jl +++ b/src/outcome_spaces/transfer_operator/transfer_operator.jl @@ -121,7 +121,7 @@ points that visits `bins[i]`. See also: [`RectangularBinning`](@ref). """ - +#should it include an encoding as well? struct TransferOperatorApproximation{T<:Real,OC<:OutcomeSpace} <: ProbabilitiesEstimator transfermatrix::AbstractArray{T, 2} outcome_space::OC @@ -194,7 +194,7 @@ struct InvariantMeasure{T} end function invariantmeasure(iv::InvariantMeasure) - return iv.ρ, iv.to.bins + return iv.ρ, iv.to.outcomes end @@ -262,7 +262,7 @@ function invariantmeasure(to::TransferOperatorApproximation; =# Ρ = 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 Ρ @@ -271,7 +271,6 @@ function invariantmeasure(to::TransferOperatorApproximation; # iterations. =# distribution = Ρ * to.transfermatrix - distance = norm(distribution - Ρ) / norm(Ρ) check = floor(Int, 1 / delta) @@ -297,7 +296,6 @@ function invariantmeasure(to::TransferOperatorApproximation; distribution = distribution ./ colsum_distribution end end - distance = norm(distribution - Ρ) / norm(Ρ) end distribution = dropdims(distribution, dims = 1) @@ -340,19 +338,14 @@ function probabilities(pest::TransferOperator, o::OutcomeSpace, x::Array_or_SSSe 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) diff --git a/test/outcome_spaces/implementations/transfer_operator.jl b/test/outcome_spaces/implementations/transfer_operator.jl index 78e78d7d..f6eb81bf 100644 --- a/test/outcome_spaces/implementations/transfer_operator.jl +++ b/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.