make ArnoldiMethod and KrylovKit optional extensions, update Readme

Project.toml

@@ -1,22 +1,28 @@
 name = "PCCA"
 uuid = "f48fc343-7e38-490c-be15-e66d68689cd5"
 authors = ["Alexander Sikorski"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
-ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
-KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+
+[extensions]
+ArnoldiExt = "ArnoldiMethod"
+KrylovExt = "KrylovKit"
+
 [compat]
 ArnoldiMethod = "0.2"
 Arpack = "0.5"
-KrylovKit = "0.5"
+KrylovKit = "0.5, 0.6"
 Optim = "1"
-julia = "1.6"
+julia = "1.9"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -4,3 +4,34 @@
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://axsk.github.io/PCCA.jl/dev)
 [![Build Status](https://github.com/axsk/PCCA.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/axsk/PCCA.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/axsk/PCCA.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/axsk/PCCA.jl)
+
+A [KISS](https://en.wikipedia.org/wiki/KISS_principle) style implementation of PCCA+ [1,2] with support for non-reversible systems [3].
+For a similar python implementation see also the [cmdtools](https://github.com/zib-cmd/cmdtools/) package.
+
+## Basic usage
+
+```julia
+using PCCA
+
+P=rand(10,10)
+P = P ./ sum(P, dims=2) # row stochastic matrix
+
+# basic PCCA+ clustering with 2 clusters (using no weighting and the ISA initial guess only)
+chi = pcca(P, 2)        # 
+
+using KrylovKit
+using SparseArray
+P = sprand(100,100, 0.1)
+P = P ./ sum(P, dims=2) # sparse row stochastic matrix
+
+# solve the PCCA+ problem weighted with the stationary density 
+# and optimize for crispness, using the KrylovKit.jl eigensolver
+chi = pcca(P, 2; pi=:auto, optimize=true, solver=KrylovSolver())
+```
+
+For sparse matrix support, add either the `ArnoldiMethod.jl` or `KrylovKit.jl` and pass the corresponding `ArnoldiSolver()` or `KrylovSolver()` as a solver.
+
+## References
+[1] 2006 - M. Weber: Meshless Methods in Conformation Dynamics
+[2] 2013 - S. Röblitz, M. Weber: Fuzzy Spectral Clustering by PCCA+
+[3] 2018 - K. Fackeldey, A. Sikorski, M. Weber: Spectral Clustering for Non-Reversible Markov Chains

ext/ArnoldiExt.jl

@@ -0,0 +1,12 @@
+module ArnoldiExt
+
+import PCCA
+import ArnoldiMethod
+
+function PCCA.schurvecs(T, n, israte, ::PCCA.ArnoldiSolver)
+    which = israte ? ArnoldiMethod.LR() : ArnoldiMethod.LM()
+    Q = ArnoldiMethod.partialschur(T; nev=n, which)[1].Q
+    Q[:, 1:n]
+end
+
+end # module

ext/KrylovExt.jl

@@ -0,0 +1,15 @@
+module KrylovExt
+import PCCA
+import KrylovKit
+
+
+
+
+function PCCA.schurvecs(T, n, israte, ::PCCA.KrylovSolver)
+    which = israte ? :LR : :LM
+    R, Qs, = KrylovKit.schursolve(T, rand(size(T, 1)), n, which, KrylovKit.Arnoldi())
+    Q = reduce(hcat, Qs)
+    Q[:, 1:n]
+end
+
+end # module

src/PCCA.jl

@@ -3,6 +3,6 @@ module PCCA
 include("pccap.jl")
 include("schur.jl")
 
-export pcca
+export pcca, KrylovSolver, ArnoldiSolver, BaseSolver
 
 end

src/pccap.jl

@@ -1,6 +1,20 @@
 import Arpack: eigs
 using LinearAlgebra
 
+
+""" pcca(T, n; pi=nothing, optimize=false, solver=BaseSolver())
+
+performs the pcca clustering on the transition matrix `T` with `n` clusters.
+Here `T` can be either a stochastic propagator with row sum 1 or a rate matrix with row sum 0.
+
+- `pi` is a density for weighting the result. `pi=:auto` uses the stationary density
+- `optimize` uses the optimization from Roeblitz (2013) to improve crispness
+- `solver` is the solver to use for computing the schur decomposition. 
+   The BaseSolver() is built in. ArnoldiSolver() and KrylovSolver() require ArnoldiMethod.jl and KrylovKit.jl respectively and provide support for sparse matrices.
+
+Returns the membership matrix `chi`.
+"""
+
 function pcca(T::AbstractMatrix, n::Integer; pi=nothing, optimize=false, solver=BaseSolver())
     israte = isratematrix(T)
     if pi == :auto

src/schur.jl

@@ -1,9 +1,6 @@
 ### Schur Eigenspace solvers
 
-import KrylovKit
-import ArnoldiMethod
-
-# TODO: both ArnoldiSolver and KrylovSolver cut the subspace, this should warn or error
+# TODO: there should be at least warnings when cutting conjugate eigenvalues
 
 """ Uses Julia's built in LinearAlgebra.schur solver. This has no support for sparse matrices """
 struct BaseSolver end
@@ -28,19 +25,6 @@ function schurvectors(T, pi::Nothing, n, israte, solver)
     return X
 end
 
-function schurvecs(T, n, israte, ::ArnoldiSolver)
-    which = israte ? ArnoldiMethod.LR() : ArnoldiMethod.LM()
-    Q = ArnoldiMethod.partialschur(T; nev=n, which)[1].Q
-    Q[:, 1:n]
-end
-
-function schurvecs(T, n, israte, ::KrylovSolver)
-    which = israte ? :LR : :LM
-    R, Qs, = KrylovKit.schursolve(T, rand(size(T, 1)), n, which, KrylovKit.Arnoldi())
-    Q = reduce(hcat, Qs)
-    Q[:, 1:n]
-end
-
 using SparseArrays
 function schurvecs(T, n, israte, ::BaseSolver)
     issparse(T) && error("The `BaseSolver` does not support sparse matrices")
@@ -49,7 +33,6 @@ function schurvecs(T, n, israte, ::BaseSolver)
     return Q
 end
 
-
 # select the schurvectors corresponding to the n largest real part of eigenvalues
 function selclusters!(S, n, israte)
     sortby = israte ? real.(S.values) : abs.(S.values)

test/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -1,6 +1,9 @@
 using PCCA
 using Test
 
+using ArnoldiMethod
+using KrylovKit
+
 @testset "PCCA.jl" begin
     # Write your tests here.
     x = rand(10, 10)
@@ -15,8 +18,6 @@ using Test
         P ./= sum(P, dims=2)
     end
 
-
-
     @testset "Reversible = $rev" for rev in [true, false]
         P = [randomstochasticmatrix(3 + mod(i, 12), rev) for i in 1:10]
         @testset "Method $method" for method in [PCCA.BaseSolver, PCCA.ArnoldiSolver, PCCA.KrylovSolver]
@@ -30,5 +31,4 @@ using Test
             end
         end
     end
-
 end