Skip to content

Commit

Permalink
Vietoris-Rips implemented and somewhat tested; symmetric Vietoris-Rip…
Browse files Browse the repository at this point in the history
…s implemented; profiling.

TODO: figure out whether/why the Vietoris-Rips streams are incorrectly sorted.
  • Loading branch information
Mikael Vejdemo-Johansson committed Nov 5, 2023
1 parent 8c27c69 commit af42cec
Show file tree
Hide file tree
Showing 9 changed files with 360 additions and 37 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -12,21 +12,19 @@ interface FiniteMetricSpace<VertexT> {
abstract fun contains(x: VertexT): Boolean

companion object {
fun <VertexT : Comparable<VertexT>, Double : Comparable<Double>>
MaximumDistanceFiltrationValue(metricSpace: FiniteMetricSpace<VertexT>):
Filtered<VertexT, Double> = Filtered { simplex: AbstractSimplex<VertexT> ->
(
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 <VertexT : Comparable<VertexT>> MaximumDistanceFiltrationValue(metricSpace: FiniteMetricSpace<VertexT>):
Filtered<VertexT, Double> =
Filtered { simplex: AbstractSimplex<VertexT> ->
(
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 }
}
)
}
}
}

Expand Down Expand Up @@ -54,3 +52,15 @@ class EuclideanMetricSpace(val points: Array<Array<Double>>) : FiniteMetricSpace

override fun contains(x: Int): Boolean = (0 <= x) and (x < size)
}

class HyperCube(val dimension: Int) : FiniteMetricSpace<Int> {
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<Int> = 0..top

override fun contains(x: Int): Boolean = (x and top) == x
}
32 changes: 28 additions & 4 deletions src/commonMain/kotlin/org/appliedtopology/tda4j/Simplex.kt
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,18 @@ open class AbstractSimplex<VertexT : Comparable<VertexT>> : Set<VertexT> {

override fun containsAll(elements: Collection<VertexT>): Boolean = _simplex.containsAll(elements)

override fun equals(other: Any?): Boolean {
val that: AbstractSimplex<VertexT>? = other as? AbstractSimplex<VertexT>
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

Expand All @@ -29,15 +41,22 @@ open class AbstractSimplex<VertexT : Comparable<VertexT>> : Set<VertexT> {
return AbstractSimplex<VertexT>(vertices)
}

val vertices: List<VertexT> = _simplex.sortedBy { it }
val vertices: List<VertexT>
get() = _simplex.sorted()

override fun toString(): String = _simplex.joinToString(
",",
"abstractSimplexOf(",
")"
)

companion object {
fun <VertexT : Comparable<VertexT>> compare(left: AbstractSimplex<VertexT>, right: AbstractSimplex<VertexT>): 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) {
Expand All @@ -54,4 +73,9 @@ typealias Simplex = AbstractSimplex<Int>
fun <VertexT : Comparable<VertexT>>abstractSimplexOf(vararg elements: VertexT): AbstractSimplex<VertexT> =
AbstractSimplex(elements.toList())

fun <VertexT : Comparable<VertexT>>abstractSimplexOf(elements: Collection<VertexT>): AbstractSimplex<VertexT> =
AbstractSimplex(elements.toList())

fun simplexOf(vararg elements: Int): Simplex = Simplex(elements.toList())

fun simplexOf(elements: Collection<Int>): Simplex = Simplex(elements.toList())
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ fun interface Filtered<VertexT : Comparable<VertexT>, FiltrationT : Comparable<F

abstract class SimplexStream<VertexT : Comparable<VertexT>, FiltrationT : Comparable<FiltrationT>> :
Filtered<VertexT, FiltrationT>, Iterable<AbstractSimplex<VertexT>> {
val comparator: Comparator<AbstractSimplex<VertexT>> =
open val comparator: Comparator<AbstractSimplex<VertexT>> =
compareBy<AbstractSimplex<VertexT>> { filtrationValue(it) }
.thenComparator({ a: AbstractSimplex<VertexT>, b: AbstractSimplex<VertexT> -> AbstractSimplex.compare(a, b) })
}
Expand Down
138 changes: 138 additions & 0 deletions src/commonMain/kotlin/org/appliedtopology/tda4j/SymmetryGroup.kt
Original file line number Diff line number Diff line change
@@ -0,0 +1,138 @@
package org.appliedtopology.tda4j

interface SymmetryGroup<GroupT, VertexT : Comparable<VertexT>> {
val elements: Collection<GroupT>

fun action(g: GroupT): (VertexT) -> VertexT

fun orbit(vertex: VertexT): Set<VertexT> = elements.map { action(it)(vertex) }.toSet()

fun orbit(simplex: AbstractSimplex<VertexT>): Set<AbstractSimplex<VertexT>> =
elements.map { g -> abstractSimplexOf(simplex.map { v -> action(g)(v) }) }.toSet()

fun representative(vertex: VertexT): VertexT = orbit(vertex).min()

fun representative(simplex: AbstractSimplex<VertexT>): AbstractSimplex<VertexT> =
orbit(simplex).minWith { left, right -> AbstractSimplex.compare(left, right) }

fun isRepresentative(vertex: VertexT): Boolean = vertex == representative(vertex)

fun isRepresentative(simplex: AbstractSimplex<VertexT>): Boolean = simplex == representative(simplex)
}

class HyperCubeSymmetry(val elementCount: Int) : SymmetryGroup<Int, Int> {
override val elements: Collection<Int> = (0..factorial(elementCount) - 1).toList()

val permutations: List<Map<Int, Int>> = 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<AbstractSimplex<Int>> = HashSet()

fun permutation(n: Int): List<Int> = 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<Int>): 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<VertexT : Comparable<VertexT>> (
val representatives: List<AbstractSimplex<VertexT>>,
val symmetryGroup: SymmetryGroup<*, VertexT>
) : Sequence<AbstractSimplex<VertexT>> {
val orbitSizes: List<Int> = representatives.map { symmetryGroup.orbit(it).size }
val orbitRanges: List<Int> = orbitSizes.scan(0) { x, y -> x + y }

val size: Int = orbitSizes.sum()

fun get(index: Int): AbstractSimplex<VertexT> {
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<AbstractSimplex<VertexT>> = sequence {
representatives.forEach({ t -> yieldAll(symmetryGroup.orbit(t)) })
}.iterator()
}

class SymmetricZomorodianIncremental<VertexT : Comparable<VertexT>>(val symmetryGroup: SymmetryGroup<*, VertexT>) : CliqueFinder<VertexT> {
override fun cliques(
metricSpace: FiniteMetricSpace<VertexT>,
maxFiltrationValue: Double,
maxDimension: Int
): Sequence<AbstractSimplex<VertexT>> {
val edges: List<Pair<Double?, Pair<VertexT, VertexT>>> =
weightedEdges(metricSpace, maxFiltrationValue).sortedBy { it.first }.toList()
val lowerNeighbors = buildMap {
edges.forEach {
dvw ->
(
getOrPut(
dvw.second.second,
defaultValue = { -> hashSetOf<VertexT>() }
) as MutableSet<VertexT>
).add(dvw.second.first)
}
}

val V: MutableSet<AbstractSimplex<VertexT>> = HashSet(metricSpace.size + edges.size)

val tasks: ArrayDeque<Pair<AbstractSimplex<VertexT>, Set<VertexT>>> = ArrayDeque(edges.size)
metricSpace.elements.forEach { vertex ->
tasks.addFirst(
Pair(
abstractSimplexOf(vertex),
lowerNeighbors.getOrElse(vertex, { -> emptySet() })
)
)
}

while (tasks.size > 0) {
val task: Pair<AbstractSimplex<VertexT>, Set<VertexT>> = tasks.removeFirst()
val tau: AbstractSimplex<VertexT> = task.first
val N: Set<VertexT> = task.second
if (symmetryGroup.isRepresentative(tau)) V.add(tau)
if (tau.size < maxDimension) {
N.forEach { v: VertexT ->
run {
val sigma: AbstractSimplex<VertexT> = tau.plus(v)
val M: Set<VertexT> =
N.intersect(lowerNeighbors.get(v)?.asIterable() ?: emptySequence<VertexT>().asIterable())
tasks.addFirst(Pair(sigma, M))
}
}
}
}

val filtered: Filtered<VertexT, Double> = FiniteMetricSpace.MaximumDistanceFiltrationValue(metricSpace)
return ExpandSequence<VertexT> (V.sortedWith(VietorisRips.getComparator(filtered)), symmetryGroup)
}
}
36 changes: 27 additions & 9 deletions src/commonMain/kotlin/org/appliedtopology/tda4j/VietorisRips.kt
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
package org.appliedtopology.tda4j

open fun interface CliqueFinder<VertexT : Comparable<VertexT>> {
fun interface CliqueFinder<VertexT : Comparable<VertexT>> {
abstract fun cliques(
metricSpace: FiniteMetricSpace<VertexT>,
maxFiltrationValue: Double,
Expand All @@ -10,12 +10,12 @@ open fun interface CliqueFinder<VertexT : Comparable<VertexT>> {
fun weightedEdges(
metricSpace: FiniteMetricSpace<VertexT>,
maxFiltrationValue: Double
): Sequence<Pair<Double, Pair<VertexT, VertexT>>> {
): Sequence<Pair<Double?, Pair<VertexT, VertexT>>> {
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<Pair<Double, Pair<VertexT, VertexT>>>
}).asSequence()
}
}

Expand All @@ -31,6 +31,15 @@ open class VietorisRips<VertexT : Comparable<VertexT>>(
cliqueFinder.cliques(metricSpace, maxFiltrationValue, maxDimension)

override fun iterator(): Iterator<AbstractSimplex<VertexT>> = simplices.iterator()

override val comparator: Comparator<AbstractSimplex<VertexT>>
get() = getComparator(this)

companion object {
fun <VertexT : Comparable<VertexT>> getComparator(filtered: Filtered<VertexT, Double>): Comparator<AbstractSimplex<VertexT>> =
compareBy<AbstractSimplex<VertexT>> { filtered.filtrationValue(it) }
.thenComparator { x: AbstractSimplex<VertexT>, y: AbstractSimplex<VertexT> -> AbstractSimplex.compare(x, y) }
}
}

class ZomorodianIncremental<VertexT : Comparable<VertexT>> : CliqueFinder<VertexT> {
Expand All @@ -39,8 +48,8 @@ class ZomorodianIncremental<VertexT : Comparable<VertexT>> : CliqueFinder<Vertex
maxFiltrationValue: Double,
maxDimension: Int
): Sequence<AbstractSimplex<VertexT>> {
val edges: List<Pair<Double, Pair<VertexT, VertexT>>> =
weightedEdges(metricSpace, maxFiltrationValue).sortedBy { it.first }.toList()
val edges: List<Pair<Double?, Pair<VertexT, VertexT>>> =
weightedEdges(metricSpace, maxFiltrationValue).sortedBy { it.first ?: Double.POSITIVE_INFINITY }.toList()
val lowerNeighbors = buildMap {
edges.forEach {
dvw ->
Expand All @@ -56,7 +65,14 @@ class ZomorodianIncremental<VertexT : Comparable<VertexT>> : CliqueFinder<Vertex
val V: MutableSet<AbstractSimplex<VertexT>> = HashSet(metricSpace.size + edges.size)

val tasks: ArrayDeque<Pair<AbstractSimplex<VertexT>, Set<VertexT>>> = 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<AbstractSimplex<VertexT>, Set<VertexT>> = tasks.removeFirst()
Expand All @@ -65,15 +81,17 @@ class ZomorodianIncremental<VertexT : Comparable<VertexT>> : CliqueFinder<Vertex
V.add(tau)
if (tau.size < maxDimension) {
N.forEach { v: VertexT ->
{
run {
val sigma: AbstractSimplex<VertexT> = tau.plus(v)
val M: Set<VertexT> = N.intersect(lowerNeighbors.get(v)?.asIterable() ?: emptySequence<VertexT>().asIterable())
val M: Set<VertexT> =
N.intersect(lowerNeighbors.get(v)?.asIterable() ?: emptySequence<VertexT>().asIterable())
tasks.addFirst(Pair(sigma, M))
}
}
}
}

return V.sortedWith(Comparator<AbstractSimplex<VertexT>> { a, b -> AbstractSimplex.compare(a, b) }).asSequence()
val filtered: Filtered<VertexT, Double> = FiniteMetricSpace.MaximumDistanceFiltrationValue(metricSpace)
return V.sortedWith(VietorisRips.getComparator(filtered)).asSequence()
}
}
32 changes: 32 additions & 0 deletions src/commonTest/kotlin/org/appliedtopology/tda4j/ProfilingSpec.kt
Original file line number Diff line number Diff line change
@@ -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<Int>(hc4, maxF, maxD, ZomorodianIncremental<Int>())
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<Int>(hc4, maxF, maxD, SymmetricZomorodianIncremental<Int>(hcs4))
val sstream = vr.simplices
var size = 0
sstream.forEach { size += 1 }
println(size)
}
})
Loading

0 comments on commit af42cec

Please sign in to comment.