From af42ceced64f5a1aede157b96033af79f8bec5ba Mon Sep 17 00:00:00 2001 From: Mikael Vejdemo-Johansson Date: Sun, 5 Nov 2023 13:17:21 -0500 Subject: [PATCH] Vietoris-Rips implemented and somewhat tested; symmetric Vietoris-Rips implemented; profiling. TODO: figure out whether/why the Vietoris-Rips streams are incorrectly sorted. --- .../tda4j/FiniteMetricSpace.kt | 40 +++-- .../org/appliedtopology/tda4j/Simplex.kt | 32 +++- .../appliedtopology/tda4j/SimplexStream.kt | 2 +- .../appliedtopology/tda4j/SymmetryGroup.kt | 138 ++++++++++++++++++ .../org/appliedtopology/tda4j/VietorisRips.kt | 36 +++-- .../appliedtopology/tda4j/ProfilingSpec.kt | 32 ++++ .../org/appliedtopology/tda4j/SimplexSpec.kt | 15 +- .../tda4j/SymmetryGroupSpec.kt | 84 +++++++++++ .../appliedtopology/tda4j/VietorisRipsSpec.kt | 18 ++- 9 files changed, 360 insertions(+), 37 deletions(-) create mode 100644 src/commonMain/kotlin/org/appliedtopology/tda4j/SymmetryGroup.kt create mode 100644 src/commonTest/kotlin/org/appliedtopology/tda4j/ProfilingSpec.kt create mode 100644 src/commonTest/kotlin/org/appliedtopology/tda4j/SymmetryGroupSpec.kt diff --git a/src/commonMain/kotlin/org/appliedtopology/tda4j/FiniteMetricSpace.kt b/src/commonMain/kotlin/org/appliedtopology/tda4j/FiniteMetricSpace.kt index d64d60f5..aebec1d3 100644 --- a/src/commonMain/kotlin/org/appliedtopology/tda4j/FiniteMetricSpace.kt +++ b/src/commonMain/kotlin/org/appliedtopology/tda4j/FiniteMetricSpace.kt @@ -12,21 +12,19 @@ interface FiniteMetricSpace { abstract fun contains(x: VertexT): Boolean companion object { - fun , Double : Comparable> - MaximumDistanceFiltrationValue(metricSpace: FiniteMetricSpace): - Filtered = Filtered { simplex: AbstractSimplex -> - ( - if (simplex.size <= 1) { - 0.0 - } else { - var distances = listOf(1.0, 2.0, 0.0) - var dist = simplex.vertices.flatMap({ v -> - simplex.vertices.filter({ it > v }).map({ w -> metricSpace.distance(v, w) }) - }) - distances.max() - } - ) as Double? - } + fun > MaximumDistanceFiltrationValue(metricSpace: FiniteMetricSpace): + Filtered = + Filtered { simplex: AbstractSimplex -> + ( + if (simplex.size <= 1) { + 0.0 + } else { + simplex.vertices.flatMap({ v -> + simplex.vertices.filter({ it > v }).map({ w -> metricSpace.distance(v, w) }) + }).maxOf { it ?: Double.POSITIVE_INFINITY } + } + ) + } } } @@ -54,3 +52,15 @@ class EuclideanMetricSpace(val points: Array>) : FiniteMetricSpace override fun contains(x: Int): Boolean = (0 <= x) and (x < size) } + +class HyperCube(val dimension: Int) : FiniteMetricSpace { + val top: Int = (1 shl dimension) - 1 + + override fun distance(x: Int, y: Int): Double? = (x xor y).countOneBits().toDouble() + + override val size: Int = (1 shl dimension) + + override val elements: Iterable = 0..top + + override fun contains(x: Int): Boolean = (x and top) == x +} diff --git a/src/commonMain/kotlin/org/appliedtopology/tda4j/Simplex.kt b/src/commonMain/kotlin/org/appliedtopology/tda4j/Simplex.kt index 5f1e807d..7e2447b2 100644 --- a/src/commonMain/kotlin/org/appliedtopology/tda4j/Simplex.kt +++ b/src/commonMain/kotlin/org/appliedtopology/tda4j/Simplex.kt @@ -18,6 +18,18 @@ open class AbstractSimplex> : Set { override fun containsAll(elements: Collection): Boolean = _simplex.containsAll(elements) + override fun equals(other: Any?): Boolean { + val that: AbstractSimplex? = other as? AbstractSimplex + if (that == null) return false + val mine = this.vertices + val thine = that.vertices + return mine.containsAll(thine) and thine.containsAll(mine) + } + + override fun hashCode(): Int { + return _simplex.hashCode() + } + override val size: Int get() = _simplex.size @@ -29,15 +41,22 @@ open class AbstractSimplex> : Set { return AbstractSimplex(vertices) } - val vertices: List = _simplex.sortedBy { it } + val vertices: List + get() = _simplex.sorted() + + override fun toString(): String = _simplex.joinToString( + ",", + "abstractSimplexOf(", + ")" + ) companion object { fun > compare(left: AbstractSimplex, right: AbstractSimplex): Int { - for (lr in left.vertices.sorted().padZip(right.vertices.sorted())) { + for (lr in left.vertices.padZip(right.vertices)) { val l = lr.first val r = lr.second - if (l == null) return +1 - if (r == null) return -1 + if (l == null) return -1 + if (r == null) return +1 val lv: VertexT = l!! val rv: VertexT = r!! if (lv.compareTo(rv) != 0) { @@ -54,4 +73,9 @@ typealias Simplex = AbstractSimplex fun >abstractSimplexOf(vararg elements: VertexT): AbstractSimplex = AbstractSimplex(elements.toList()) +fun >abstractSimplexOf(elements: Collection): AbstractSimplex = + AbstractSimplex(elements.toList()) + fun simplexOf(vararg elements: Int): Simplex = Simplex(elements.toList()) + +fun simplexOf(elements: Collection): Simplex = Simplex(elements.toList()) diff --git a/src/commonMain/kotlin/org/appliedtopology/tda4j/SimplexStream.kt b/src/commonMain/kotlin/org/appliedtopology/tda4j/SimplexStream.kt index af284ce2..2c899b84 100644 --- a/src/commonMain/kotlin/org/appliedtopology/tda4j/SimplexStream.kt +++ b/src/commonMain/kotlin/org/appliedtopology/tda4j/SimplexStream.kt @@ -6,7 +6,7 @@ fun interface Filtered, FiltrationT : Comparable, FiltrationT : Comparable> : Filtered, Iterable> { - val comparator: Comparator> = + open val comparator: Comparator> = compareBy> { filtrationValue(it) } .thenComparator({ a: AbstractSimplex, b: AbstractSimplex -> AbstractSimplex.compare(a, b) }) } diff --git a/src/commonMain/kotlin/org/appliedtopology/tda4j/SymmetryGroup.kt b/src/commonMain/kotlin/org/appliedtopology/tda4j/SymmetryGroup.kt new file mode 100644 index 00000000..dd9e0d27 --- /dev/null +++ b/src/commonMain/kotlin/org/appliedtopology/tda4j/SymmetryGroup.kt @@ -0,0 +1,138 @@ +package org.appliedtopology.tda4j + +interface SymmetryGroup> { + val elements: Collection + + fun action(g: GroupT): (VertexT) -> VertexT + + fun orbit(vertex: VertexT): Set = elements.map { action(it)(vertex) }.toSet() + + fun orbit(simplex: AbstractSimplex): Set> = + elements.map { g -> abstractSimplexOf(simplex.map { v -> action(g)(v) }) }.toSet() + + fun representative(vertex: VertexT): VertexT = orbit(vertex).min() + + fun representative(simplex: AbstractSimplex): AbstractSimplex = + orbit(simplex).minWith { left, right -> AbstractSimplex.compare(left, right) } + + fun isRepresentative(vertex: VertexT): Boolean = vertex == representative(vertex) + + fun isRepresentative(simplex: AbstractSimplex): Boolean = simplex == representative(simplex) +} + +class HyperCubeSymmetry(val elementCount: Int) : SymmetryGroup { + override val elements: Collection = (0..factorial(elementCount) - 1).toList() + + val permutations: List> = elements.map { + run { + val p = permutation(it) + (0..(1 shl elementCount) - 1).associateWith { bits -> + (0..elementCount - 1) + .map { ((bits and (1 shl it) shr it) shl p[it]) } + .fold(0, { x, y -> x or y }) + } + } + } + val knownRepresentatives: MutableSet> = HashSet() + + fun permutation(n: Int): List = buildList(elementCount) { + val source = (0..elementCount - 1).toMutableList() + var pos = n + + while (source.isNotEmpty()) { + val div = pos.floorDiv(factorial(source.size - 1)) + pos = pos.mod(factorial(source.size - 1)) + add(source.removeAt(div)) + } + } + + override fun action(g: Int): (Int) -> Int = { permutations[g][it] ?: it } + + override fun isRepresentative(simplex: AbstractSimplex): Boolean { + if (simplex in knownRepresentatives) return true + val isR = super.isRepresentative(simplex) + if (isR) knownRepresentatives.add(simplex) + return isR + } + + companion object { + fun factorial(n: Int, accum: Int = 1): Int { + return if (n <= 1) accum else factorial(n - 1, n * accum) + } + } +} + +class ExpandSequence> ( + val representatives: List>, + val symmetryGroup: SymmetryGroup<*, VertexT> +) : Sequence> { + val orbitSizes: List = representatives.map { symmetryGroup.orbit(it).size } + val orbitRanges: List = orbitSizes.scan(0) { x, y -> x + y } + + val size: Int = orbitSizes.sum() + + fun get(index: Int): AbstractSimplex { + val orbitIndex = orbitRanges.indexOfFirst { it >= index } - 1 + return symmetryGroup.orbit(representatives[orbitIndex]).toList()[index - orbitRanges[orbitIndex]] + } + + fun isEmpty(): Boolean = size == 0 + + override fun iterator(): Iterator> = sequence { + representatives.forEach({ t -> yieldAll(symmetryGroup.orbit(t)) }) + }.iterator() +} + +class SymmetricZomorodianIncremental>(val symmetryGroup: SymmetryGroup<*, VertexT>) : CliqueFinder { + override fun cliques( + metricSpace: FiniteMetricSpace, + maxFiltrationValue: Double, + maxDimension: Int + ): Sequence> { + val edges: List>> = + weightedEdges(metricSpace, maxFiltrationValue).sortedBy { it.first }.toList() + val lowerNeighbors = buildMap { + edges.forEach { + dvw -> + ( + getOrPut( + dvw.second.second, + defaultValue = { -> hashSetOf() } + ) as MutableSet + ).add(dvw.second.first) + } + } + + val V: MutableSet> = HashSet(metricSpace.size + edges.size) + + val tasks: ArrayDeque, Set>> = ArrayDeque(edges.size) + metricSpace.elements.forEach { vertex -> + tasks.addFirst( + Pair( + abstractSimplexOf(vertex), + lowerNeighbors.getOrElse(vertex, { -> emptySet() }) + ) + ) + } + + while (tasks.size > 0) { + val task: Pair, Set> = tasks.removeFirst() + val tau: AbstractSimplex = task.first + val N: Set = task.second + if (symmetryGroup.isRepresentative(tau)) V.add(tau) + if (tau.size < maxDimension) { + N.forEach { v: VertexT -> + run { + val sigma: AbstractSimplex = tau.plus(v) + val M: Set = + N.intersect(lowerNeighbors.get(v)?.asIterable() ?: emptySequence().asIterable()) + tasks.addFirst(Pair(sigma, M)) + } + } + } + } + + val filtered: Filtered = FiniteMetricSpace.MaximumDistanceFiltrationValue(metricSpace) + return ExpandSequence (V.sortedWith(VietorisRips.getComparator(filtered)), symmetryGroup) + } +} diff --git a/src/commonMain/kotlin/org/appliedtopology/tda4j/VietorisRips.kt b/src/commonMain/kotlin/org/appliedtopology/tda4j/VietorisRips.kt index a2524703..c2219d2a 100644 --- a/src/commonMain/kotlin/org/appliedtopology/tda4j/VietorisRips.kt +++ b/src/commonMain/kotlin/org/appliedtopology/tda4j/VietorisRips.kt @@ -1,6 +1,6 @@ package org.appliedtopology.tda4j -open fun interface CliqueFinder> { +fun interface CliqueFinder> { abstract fun cliques( metricSpace: FiniteMetricSpace, maxFiltrationValue: Double, @@ -10,12 +10,12 @@ open fun interface CliqueFinder> { fun weightedEdges( metricSpace: FiniteMetricSpace, maxFiltrationValue: Double - ): Sequence>> { + ): Sequence>> { return metricSpace.elements.flatMap({ v -> metricSpace.elements.filter({ v < it }) .map({ w -> Pair(metricSpace.distance(v, w), Pair(v, w)) }) .filter({ (d, _) -> (d != null) && (d <= maxFiltrationValue) }) - }) as Sequence>> + }).asSequence() } } @@ -31,6 +31,15 @@ open class VietorisRips>( cliqueFinder.cliques(metricSpace, maxFiltrationValue, maxDimension) override fun iterator(): Iterator> = simplices.iterator() + + override val comparator: Comparator> + get() = getComparator(this) + + companion object { + fun > getComparator(filtered: Filtered): Comparator> = + compareBy> { filtered.filtrationValue(it) } + .thenComparator { x: AbstractSimplex, y: AbstractSimplex -> AbstractSimplex.compare(x, y) } + } } class ZomorodianIncremental> : CliqueFinder { @@ -39,8 +48,8 @@ class ZomorodianIncremental> : CliqueFinder> { - val edges: List>> = - weightedEdges(metricSpace, maxFiltrationValue).sortedBy { it.first }.toList() + val edges: List>> = + weightedEdges(metricSpace, maxFiltrationValue).sortedBy { it.first ?: Double.POSITIVE_INFINITY }.toList() val lowerNeighbors = buildMap { edges.forEach { dvw -> @@ -56,7 +65,14 @@ class ZomorodianIncremental> : CliqueFinder> = HashSet(metricSpace.size + edges.size) val tasks: ArrayDeque, Set>> = ArrayDeque(edges.size) - lowerNeighbors.forEach { entry -> tasks.addFirst(Pair(abstractSimplexOf(entry.key), entry.value)) } + metricSpace.elements.forEach { vertex -> + tasks.addFirst( + Pair( + abstractSimplexOf(vertex), + lowerNeighbors.getOrElse(vertex, { -> emptySet() }) + ) + ) + } while (tasks.size > 0) { val task: Pair, Set> = tasks.removeFirst() @@ -65,15 +81,17 @@ class ZomorodianIncremental> : CliqueFinder - { + run { val sigma: AbstractSimplex = tau.plus(v) - val M: Set = N.intersect(lowerNeighbors.get(v)?.asIterable() ?: emptySequence().asIterable()) + val M: Set = + N.intersect(lowerNeighbors.get(v)?.asIterable() ?: emptySequence().asIterable()) tasks.addFirst(Pair(sigma, M)) } } } } - return V.sortedWith(Comparator> { a, b -> AbstractSimplex.compare(a, b) }).asSequence() + val filtered: Filtered = FiniteMetricSpace.MaximumDistanceFiltrationValue(metricSpace) + return V.sortedWith(VietorisRips.getComparator(filtered)).asSequence() } } diff --git a/src/commonTest/kotlin/org/appliedtopology/tda4j/ProfilingSpec.kt b/src/commonTest/kotlin/org/appliedtopology/tda4j/ProfilingSpec.kt new file mode 100644 index 00000000..41175d06 --- /dev/null +++ b/src/commonTest/kotlin/org/appliedtopology/tda4j/ProfilingSpec.kt @@ -0,0 +1,32 @@ +package org.appliedtopology.tda4j + +import io.kotest.core.spec.style.FunSpec + +val dimension: Int = 4 +val maxF = dimension.toDouble() +val maxD = 1 shl dimension +class ZomorodianProfilingSpec : FunSpec({ + val hc4 = HyperCube(dimension) + val hcs4 = HyperCubeSymmetry(dimension) + + test("Time VietorisRips construction") { + val vr = VietorisRips(hc4, maxF, maxD, ZomorodianIncremental()) + val sstream = vr.simplices + var size = 0 + sstream.forEach { size += 1 } + println(size) + } +}) + +class SymmetricZomorodianProfilingSpec : FunSpec({ + val hc4 = HyperCube(dimension) + val hcs4 = HyperCubeSymmetry(dimension) + + test("Time SymmetricVietorisRips construction") { + val vr = VietorisRips(hc4, maxF, maxD, SymmetricZomorodianIncremental(hcs4)) + val sstream = vr.simplices + var size = 0 + sstream.forEach { size += 1 } + println(size) + } +}) diff --git a/src/commonTest/kotlin/org/appliedtopology/tda4j/SimplexSpec.kt b/src/commonTest/kotlin/org/appliedtopology/tda4j/SimplexSpec.kt index c4e61eaf..1cf2bd56 100644 --- a/src/commonTest/kotlin/org/appliedtopology/tda4j/SimplexSpec.kt +++ b/src/commonTest/kotlin/org/appliedtopology/tda4j/SimplexSpec.kt @@ -1,13 +1,12 @@ -package org.appliedtopology.tda4j.org.appliedtopology.tda4j +package org.appliedtopology.tda4j import io.kotest.core.spec.style.FunSpec import io.kotest.matchers.booleans.shouldBeTrue import io.kotest.matchers.comparables.shouldBeGreaterThan +import io.kotest.matchers.equals.shouldBeEqual import io.kotest.matchers.equals.shouldNotBeEqual +import io.kotest.matchers.ints.shouldBeLessThan import io.kotest.matchers.types.shouldBeInstanceOf -import org.appliedtopology.tda4j.Chain -import org.appliedtopology.tda4j.Simplex -import org.appliedtopology.tda4j.simplexOf class SimplexSpec : FunSpec({ val simplex: Simplex = simplexOf(1, 2, 3) @@ -30,4 +29,12 @@ class SimplexSpec : FunSpec({ test("A simplex can remove a vertex") { (simplex - 2).equals(listOf(1, 3)).shouldBeTrue() } + + test("Two simplices can be equal") { + simplexOf(0, 1, 2).shouldBeEqual(simplexOf(0, 1, 2)) + } + + test("Lexicographic ordering works") { + Simplex.compare(simplexOf(0, 1), simplexOf(0, 2)).shouldBeLessThan(0) + } }) diff --git a/src/commonTest/kotlin/org/appliedtopology/tda4j/SymmetryGroupSpec.kt b/src/commonTest/kotlin/org/appliedtopology/tda4j/SymmetryGroupSpec.kt new file mode 100644 index 00000000..fdfd4221 --- /dev/null +++ b/src/commonTest/kotlin/org/appliedtopology/tda4j/SymmetryGroupSpec.kt @@ -0,0 +1,84 @@ +package org.appliedtopology.tda4j + +import io.kotest.core.spec.style.FunSpec +import io.kotest.matchers.collections.shouldContainAll +import io.kotest.matchers.collections.shouldContainExactly +import io.kotest.matchers.collections.shouldContainInOrder +import io.kotest.matchers.equals.shouldBeEqual + +class SymmetryGroupSpec : FunSpec({ + test("Factorial function computes correctly") { + (1..5).map({ HyperCubeSymmetry.factorial(it) }) + .shouldContainInOrder(listOf(1, 2, 6, 24, 120)) + } + + val hcs = HyperCubeSymmetry(3) + test("Permutations of 1,2,3 enumerate correctly") { + hcs.elements.map { hcs.permutation(it) }.shouldContainInOrder( + listOf( + listOf(0, 1, 2), + listOf(0, 2, 1), + listOf(1, 0, 2), + listOf(1, 2, 0), + listOf(2, 0, 1), + listOf(2, 1, 0) + ) + ) + } + + test("Permutations of 1,2,3 shift bits correctly") { + hcs.orbit(3).shouldContainExactly(listOf(3, 5, 6)) + } + + val hc2 = HyperCube(2) + val hcs2 = HyperCubeSymmetry(2) + val sstream = VietorisRips(hc2, 3.0, 3, SymmetricZomorodianIncremental(hcs2)) + val expandseq = sstream.simplices as ExpandSequence + + test("An orbit can have a representative simplex") { + hcs2.isRepresentative(abstractSimplexOf(0)) + } + + test("Orbits work as expected") { + hcs2.orbit(abstractSimplexOf(0, 1, 3)).shouldContainExactly( + abstractSimplexOf(0, 1, 3), + abstractSimplexOf(0, 2, 3) + ) + } + + test("Singleton orbits work as expected") { + hcs2.orbit(abstractSimplexOf(0, 1, 2)).shouldBeEqual(setOf(abstractSimplexOf(0, 1, 2))) + } + + test("SimplexStream should have the right representatives") { + expandseq.representatives.shouldContainAll( + abstractSimplexOf(0), + abstractSimplexOf(1), + abstractSimplexOf(3), + abstractSimplexOf(0, 1), + abstractSimplexOf(1, 3), + abstractSimplexOf(0, 3), + abstractSimplexOf(0, 1, 3), + abstractSimplexOf(0, 1, 2) + ) + } + + test("SimplexStream should have the right elements") { + sstream.simplices.toList().shouldContainAll( + abstractSimplexOf(0), + abstractSimplexOf(1), + abstractSimplexOf(2), + abstractSimplexOf(3), + abstractSimplexOf(0, 1), + abstractSimplexOf(0, 2), + abstractSimplexOf(1, 3), + abstractSimplexOf(2, 3), + abstractSimplexOf(0, 3), + abstractSimplexOf(1, 2), + abstractSimplexOf(0, 1, 2), + abstractSimplexOf(0, 1, 3), + abstractSimplexOf(0, 2, 3), + abstractSimplexOf(1, 2, 3) + ) + } +}) diff --git a/src/commonTest/kotlin/org/appliedtopology/tda4j/VietorisRipsSpec.kt b/src/commonTest/kotlin/org/appliedtopology/tda4j/VietorisRipsSpec.kt index d86ab2e0..bc5d7dd1 100644 --- a/src/commonTest/kotlin/org/appliedtopology/tda4j/VietorisRipsSpec.kt +++ b/src/commonTest/kotlin/org/appliedtopology/tda4j/VietorisRipsSpec.kt @@ -1,10 +1,20 @@ -package org.appliedtopology.tda4j.org.appliedtopology.tda4j +package org.appliedtopology.tda4j import io.kotest.core.spec.style.FunSpec -import io.kotest.matchers.booleans.shouldBeFalse +import io.kotest.matchers.collections.shouldContainAll class VietorisRipsSpec() : FunSpec({ - test("Need to do a test here") { - true.shouldBeFalse() + val hc2 = HyperCube(2) + val sstream = VietorisRips(hc2, 3.0, 3, ZomorodianIncremental()) + + test("HyperCube(2) generates the right simplices") { + sstream.simplices.toList().shouldContainAll( + simplexOf(0), simplexOf(1), simplexOf(2), simplexOf(3), + simplexOf(0, 1), simplexOf(0, 2), + simplexOf(1, 3), simplexOf(2, 3), + simplexOf(0, 3), simplexOf(1, 2), + simplexOf(0, 1, 3), simplexOf(0, 2, 3), + simplexOf(0, 1, 2), simplexOf(1, 2, 3) + ) } })