moonlight-triangulation-1.3.0.2: bench/alpha/Main.hs

-- | Measure the exact alpha shortcut against recomputing homology at every
-- critical radius. Delaunay construction is deliberately outside every timed
-- lane: both algorithms consume the same already-built geometry.
module Main (main) where

import BenchSupport
  ( randomPoints
  , requireRight
  , timedProjection
  , timedValue
  )
import Data.Map.Strict qualified as Map
import Data.Vector qualified as Vector
import Moonlight.Homology.Boundary
  ( degreeCardinality
  , maxHomologicalDegree
  , restrictComplex
  )
import Moonlight.Homology.Chain
  ( HomologicalDegree (..)
  , HomologyFailure
  )
import Moonlight.Homology.Persistence
  ( FilteredFiniteChainComplex
  , criticalBettiVectors
  , filteredBaseComplex
  , filteredCellBirths
  , filteredCriticalValues
  , mod2PersistentPairsWithCriticalBettiTable
  )
import Moonlight.Homology.Topology (freeBettiVector)
import Moonlight.Triangulation.Alpha
  ( AlphaBirth
  , alphaFiltration
  )
import Moonlight.Triangulation.BulkLoad (delaunayGeometry)
import Moonlight.Triangulation.CellComplex (filteredAlphaComplex)
import Moonlight.Triangulation.Types (DelaunayTriangulation)

type PreparedAlpha = FilteredFiniteChainComplex AlphaBirth Int

main :: IO ()
main = benchmarkAlphaFiltration 64

benchmarkAlphaFiltration :: Int -> IO ()
benchmarkAlphaFiltration siteCount = do
  triangulation <-
    requireRight
      (delaunayGeometry (Vector.fromList (randomPoints 0x9e3779b97f4a7c15 siteCount)))
  filtered <-
    timedProjection
      (benchmarkLabel siteCount "alpha-filtered-complex")
      preparationReceipt
      (prepareAlpha triangulation)
  let thresholds = filteredCriticalValues filtered
      maximumDegree = maxHomologicalDegreeValue filtered
  persistentProfile <-
    timedValue
      (benchmarkLabel siteCount "persistence-critical-profile")
      ( do
          (_, bettiTable) <-
            requireRight (mod2PersistentPairsWithCriticalBettiTable filtered)
          pure (criticalBettiVectors bettiTable)
      )
  repeatedProfile <-
    timedValue
      (benchmarkLabel siteCount "repeated-threshold-homology")
      (requireRight (traverse (bettiVectorAtThreshold maximumDegree filtered) thresholds))
  if persistentProfile == repeatedProfile
    then
      putStrLn
        ( benchmarkLabel siteCount "agreement"
            <> ": critical-radii="
            <> show (length thresholds)
            <> " profiles="
            <> show (length persistentProfile)
            <> " prepared-cells="
            <> show (preparationReceipt filtered)
        )
    else fail (benchmarkLabel siteCount "profile-mismatch")

benchmarkLabel :: Int -> String -> String
benchmarkLabel siteCount suffix =
  "alpha/n" <> show siteCount <> "/" <> suffix

prepareAlpha
  :: DelaunayTriangulation ()
  -> IO PreparedAlpha
prepareAlpha triangulation = do
  filtration <- requireRight (alphaFiltration triangulation)
  requireRight (filteredAlphaComplex filtration)

preparationReceipt :: PreparedAlpha -> (Int, Int, [Int])
preparationReceipt filtered =
  let finite = filteredBaseComplex filtered
      HomologicalDegree maximumDegree = maxHomologicalDegree finite
   in ( length (filteredCriticalValues filtered)
      , Map.size (filteredCellBirths filtered)
      , fmap
          (degreeCardinality finite . HomologicalDegree)
          [0 .. maximumDegree]
      )

maxHomologicalDegreeValue
  :: FilteredFiniteChainComplex filtration coefficient
  -> Int
maxHomologicalDegreeValue filtered =
  case maxHomologicalDegree (filteredBaseComplex filtered) of
    HomologicalDegree degreeValue -> degreeValue

bettiVectorAtThreshold
  :: Int
  -> FilteredFiniteChainComplex AlphaBirth Int
  -> AlphaBirth
  -> Either HomologyFailure [Int]
bettiVectorAtThreshold maximumDegree filtered threshold = do
  restricted <-
    restrictComplex
      (Map.keysSet (Map.filter (<= threshold) (filteredCellBirths filtered)))
      (filteredBaseComplex filtered)
  pure (padBettiVector maximumDegree (freeBettiVector restricted))

padBettiVector :: Int -> [Int] -> [Int]
padBettiVector maximumDegree bettiVector =
  bettiVector
    <> replicate
      (max 0 (maximumDegree + 1 - length bettiVector))
      0