Merge pull request #69 from sandreza/sparse_operations

Sparse operations

Project.toml

@@ -1,7 +1,7 @@
 name = "MarkovChainHammer"
 uuid = "38c40fd0-bccb-4723-b30d-b2caea0ad8d9"
 authors = ["Andre Souza <[email protected]>"]
-version = "0.0.11"
+version = "0.0.12"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -13,15 +13,16 @@ ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-julia = "1.8, 1.9, 1.10"
 Distributions = "0.25"
 Documenter = "0.27"
 DocumenterTools = "0.1"
+PrecompileTools = "1.2"
 ProgressBars = "1.4"
-Revise = "3.4, 3.5"
-PrecompileTools = "1.2" 
 Reexport = "1.2"
+Revise = "3.4, 3.5"
+julia = "1.8, 1.9, 1.10"

src/TransitionMatrix/TransitionMatrix.jl

@@ -1,5 +1,6 @@
 @reexport module TransitionMatrix
 
 include("construct_from_data.jl")
+include("sparse_operations.jl")
 
 end # module TransitionMatrix

src/TransitionMatrix/sparse_operations.jl

@@ -0,0 +1,95 @@
+using SparseArrays, ProgressBars
+export sparse_perron_frobenius, sparse_generator
+
+function sparse_count_operator(markov_chain::Vector{S}, number_of_states::S, step::Int) where S
+    count_matrix = spzeros(typeof(markov_chain[1]), number_of_states, number_of_states)
+    for i in ProgressBar(0:step-1)
+        reduced_markov_chain = markov_chain[1+i:step:end]
+        for j in 1:length(reduced_markov_chain)-1
+            count_matrix[reduced_markov_chain[j+1], reduced_markov_chain[j]] += 1
+        end
+    end
+    return count_matrix
+end
+
+function sparse_count_operator(markov_chain::Vector{S}, number_of_states::S) where S
+    @info "computing sparse count matrix"
+    count_matrix = spzeros(S, number_of_states, number_of_states)
+    for i in ProgressBar(1:length(markov_chain)-1)
+        count_matrix[markov_chain[i+1], markov_chain[i]] += 1
+    end
+    return count_matrix
+end
+
+sparse_count_operator(markov_chain) = sparse_count_operator(markov_chain, maximum(markov_chain))
+
+"""
+sparse_perron_frobenius(markov_chain; step = 1)
+
+# Description 
+    Calculate the perron-frobenius matrix from a markov chain in sparse format
+
+# Arguments
+- `markov_chain::AbstractVector`: A vector of integers representing the state of a markov chain at each time step.
+
+# Keyword Arguments
+- `step::Integer=1`: The step size of the constructed operator.
+
+# Returns
+- `perron_frobenius_matrix::Matrix`: The perron-frobenius matrix of the markov chain in sparse format
+"""
+function sparse_perron_frobenius(partitions::Vector{S}; step = 1) where S
+    number_of_states = maximum(partitions)
+    count_matrix = sparse_count_operator(partitions, number_of_states, step)
+    number_of_states = maximum(partitions)
+    perron_frobenius_matrix = Float64.(count_matrix) 
+    normalization = sum(count_matrix, dims=1)
+    for i in ProgressBar(eachindex(normalization))
+        for j in perron_frobenius_matrix[:, i].nzind
+            perron_frobenius_matrix[j, i] /= normalization[i]
+        end
+    end
+    return perron_frobenius_matrix
+end
+
+"""
+sparse_generator(markov_chain; dt = 1)
+
+# Description 
+    Calculate the generator matrix from a markov chain in sparse format
+
+# Arguments
+- `markov_chain::AbstractVector`: A vector of integers representing the state of a markov chain at each time step.
+
+# Keyword Arguments
+- `dt::Real=1`: The time step of the constructed operator.
+
+# Returns
+- `generator_matrix::Matrix`: The generator matrix of the markov chain in sparse format
+"""
+function sparse_generator(partitions::Vector{S}; dt = 1) where S
+    number_of_states = maximum(partitions)
+    count_matrix = sparse_count_operator(partitions, number_of_states)
+    generator_matrix = Float64.(count_matrix) 
+    for i in 1:number_of_states
+        count_matrix[i, i] = 0.0
+    end
+    normalization = sum(count_matrix, dims=1)
+    holding_scale = 1 ./ mean.(holding_times(partitions, number_of_states; dt=dt))
+    # calculate generator and handle edge case where no transitions occur
+    for i in ProgressBar(eachindex(normalization))
+        for j in generator_matrix[:, i].nzind
+            generator_matrix[j, i] /= normalization[i]
+        end
+        generator_matrix[i, i] = -1.0
+        for j in generator_matrix[:, i].nzind
+            generator_matrix[j, i] *= holding_scale[i]
+        end
+        if normalization[i] == 0.0
+            for j in generator_matrix[:, i].nzind
+                generator_matrix[j, i] *= false
+            end
+        end
+    end
+    return generator_matrix
+end

src/Utils/histogram.jl

@@ -1,13 +1,15 @@
 export histogram
 
 """
-    histogram(array, bins=minimum([100, length(array)]), normalization=:uniform, custom_range=false)
+    histogram(array; bins=minimum([100, length(array)]), normalization=:uniform, custom_range=false)
 
 # Description
     Compute the histogram of an array. Useful for barplot in GLMakie.
 
 # Arguments
 - `array::AbstractArray`: Array to compute the histogram of.
+
+# Keyword Arguments
 - `bins::Integer`: Number of bins to use.
 - `normalization::AbstractArray`: Normalization to use. If `:uniform`, then the normalization is uniform.
 - `custom_range::Tuple`: Custom range to use. If `false`, then the range is computed from the data.

src/Utils/statistics.jl

@@ -53,12 +53,11 @@ end
 - `Q::Matrix`: The generator matrix or transition matrix.
 
 # Returns
-- `V⁻¹::Matrix`: The koopman modes of the generator matrix. Each row is a koopman mode
+- `W::Matrix`: The koopman modes of the generator matrix. Each column is a koopman mode
 """
 function koopman_modes(Q)
-    Λ, V = eigen(Q)
-    V⁻¹ = inv(V)
-    return V⁻¹
+    Λ, W = eigen(Q')
+    return W
 end
 
 """

test/test_TransitionMatrix.jl

@@ -1,5 +1,5 @@
 using MarkovChainHammer, Test, Revise, LinearAlgebra
-import MarkovChainHammer.TransitionMatrix: count_operator, generator, holding_times, perron_frobenius
+import MarkovChainHammer.TransitionMatrix: count_operator, generator, holding_times, perron_frobenius, sparse_count_operator
 import MarkovChainHammer.TransitionMatrix: symmetric_generator, symmetric_perron_frobenius
 
 @testset "Column Sum Consistency" begin
@@ -152,3 +152,22 @@ end
     @test all(generator_computed2 .== generator_computed3)
     @test all(generator_computed1 .== (0.5 * generator_computed4))
 end
+
+@testset "Hand Constructed: Default Case Sparse Check" begin
+    # default case. See all states in timeseries
+    timeseries = [1, 1, 1, 2, 2, 3, 3, 3, 2, 1]
+    count_operator_exact = [2.0 1.0 0.0; 1.0 1.0 1.0; 0.0 1.0 2.0]
+    generator_exact = [-1/2 1/3 0.0; 1/2 -2/3 1/3; 0.0 1/3 -1/3]
+    perron_frobenius_exact = [2/3 1/3 0.0; 1/3 1/3 1/3; 0.0 1/3 2/3]
+    perron_frobenius_exact_2step = [1/3 0 1/3; 2/3 0 1/3; 0 1 1/3]
+
+    count_operator_computed = sparse_count_operator(timeseries)
+    generator_computed = sparse_generator(timeseries)
+    perron_frobenius_computed = sparse_perron_frobenius(timeseries)
+    perron_frobenius_computed_2step = sparse_perron_frobenius(timeseries; step = 2)
+
+    @test all(count_operator_exact .== count_operator_computed)
+    @test all(generator_exact .== generator_computed)
+    @test all(perron_frobenius_exact .== perron_frobenius_computed)
+    @test all(perron_frobenius_exact_2step .== perron_frobenius_computed_2step)
+end

test/test_Utils.jl

@@ -26,7 +26,7 @@ end
     Q = [-1.0 2.0; 1.0 -2.0]
     p = steady_state(Q)
     km = koopman_modes(Q)
-    @test abs(km[1, :]' * p) < 100 * eps(1.0)
+    @test abs(km[:, 1]' * p) < 100 * eps(1.0)
 end
 
 @testset "utilities: decorrelation times" begin