AlgebraFixtures.hs

-- | The operand meshes and mesh observations both algebra slices are stated
-- over, built against the facade alone.
module Moonlight.Triangulation.AlgebraFixtures
  ( Mesh
  , PointMesh
  , meshOf
  , pointMeshOf
  , operands
  , separatedOperands
  , cocircularRing
  , collinearSites
  , latticeSites
  , pointsOf
  , dedupeAscending
  , siteList
  , siteSet
  , edgeKeys
  , assertMesh
  , assertMeshEquivalent
  , pairs
  , triples
  , advance
  , scramble
  , randomSites
  ) where

import Data.List (sortBy)
import Data.Ord (comparing)
import qualified Data.Set as Set
import qualified Data.Vector as V
import Data.Word (Word64)
import Moonlight.Triangulation
  ( DelaunayTriangulation
  , Point (Point)
  , VertexId
  , buildTriangulation
  , canonicalize
  , delaunay
  , mapVertices
  , numFaces
  , numUndirectedEdges
  , numVertices
  , undirectedEdges
  , undirectedEndpoints
  , unitElementDefaults
  , vertexPoint
  , vertices
  )
import Support (requireRight)

-- | The carrier. Geometry and nothing else: no vertex payload to need a
-- commutative combining rule, no element payloads to survive a rewrite that
-- destroys the elements they labelled.
type Mesh = DelaunayTriangulation ()

-- | The same carrier with its exact coordinate retained as the vertex
-- annotation. Annotation-preserving set laws use it to construct their
-- expected values without reaching below the public facade.
type PointMesh = DelaunayTriangulation (Point)

-- ── operands ─────────────────────────────────────────────────────────────────

-- | The site families the laws are exercised over.
--
-- Sized so that the cubic sweep in the associativity law stays cheap. Coverage
-- here is by /kind/ of degeneracy, not by count: what breaks a join is a
-- cocircular quad whose tie-break went the other way, or a mesh with no faces
-- at all, and neither of those becomes more likely at ten thousand sites.
families :: [(String, [Point])]
families =
  [ ("void", [])
  , ("single", [Point 3 (-7)])
  , ("pair", [Point 0 0, Point 4 1])
  , ("collinear", collinearSites)
  , ("cocircular", cocircularRing)
  , ("lattice", latticeSites)
  , ("repeated", concatMap (replicate 3) [Point 0 0, Point 5 0, Point 0 5, Point 5 5, Point 2 3])
  , ("scattered", randomSites 0xC0FFEEBABE 32)
  , ("extreme", [Point 1.0e-8 1.0e-8, Point 1.0e8 (-1.0e8), Point (-1.0e8) 1.0e8, Point 0 0, Point 1 1])
  ]

-- | Exactly cocircular lattice points on @x² + y² = 625@.
--
-- Trigonometric points would be cocircular only to within rounding, and the
-- rule this is here to exercise — the diagonal tie-break that fires when the
-- lifted quadrilateral is exactly flat — would then never fire at all.
cocircularRing :: [Point]
cocircularRing =
  [ Point (fromIntegral (signX * x)) (fromIntegral (signY * y))
  | (x, y) <- [(25, 0), (0, 25), (7, 24), (24, 7), (15, 20), (20, 15)] :: [(Int, Int)]
  , signX <- [1, -1]
  , signY <- [1, -1]
  ]

collinearSites :: [Point]
collinearSites = [Point (fromIntegral k) (2 * fromIntegral k - 1) | k <- [0 .. 9 :: Int]]

latticeSites :: [Point]
latticeSites = [Point (fromIntegral i) (fromIntegral j) | i <- [0 .. 4 :: Int], j <- [0 .. 4 :: Int]]

-- | The meshes every pairwise and triple law runs over.
--
-- Chosen for overlap structure rather than for size: disjoint operands, nested
-- ones, partially overlapping ones, the empty one, and — the case that matters
-- most — two meshes standing on the /same/ sites built in different orders, so
-- they are geometrically identical and structurally distinct.
operands :: IO [(String, Mesh)]
operands = do
  let sites name = maybe [] id (lookup name families)
      lattice = sites "lattice"
  void' <- meshOf "void" []
  single <- meshOf "single" (sites "single")
  collinear <- meshOf "collinear" (sites "collinear")
  ring <- meshOf "cocircular" (sites "cocircular")
  lower <- meshOf "lattice-lower" (take 15 lattice)
  upper <- meshOf "lattice-upper" (drop 10 lattice)
  repeated <- meshOf "repeated" (sites "repeated")
  scatterA <- meshOf "scattered" (sites "scattered")
  scatterB <- meshOf "scattered-scrambled" (scramble 0x5EED (sites "scattered"))
  pure
    [ ("void", void')
    , ("single", single)
    , ("collinear", collinear)
    , ("cocircular", ring)
    , ("lattice-lower", lower)
    , ("lattice-upper", upper)
    , ("repeated", repeated)
    , ("scattered", scatterA)
    , ("scattered-scrambled", scatterB)
    ]

-- | Operand pairs whose sites are separated by a vertical line, which is the
-- stratum a seam schedule is defined on and the one the old-edge identity is
-- sharpest over: with no shared site, every edge is unambiguously left, right
-- or cross.
separatedOperands :: IO [((String, Mesh), (String, Mesh))]
separatedOperands =
  traverse
    ( \(name, seed, count, shift) -> do
        let sites = randomSites seed count
        left <- meshOf (name <> "-left") sites
        right <- meshOf (name <> "-right") [Point (x + shift) y | Point x y <- sites]
        pure ((name <> "-left", left), (name <> "-right", right))
    )
    [ ("split-distant", 0x51DE1, 60, 1000)
    , ("split-near", 0x51DE2, 60, 3)
    , ("split-abutting", 0x51DE3, 60, 2.05)
    ]

-- ── construction ─────────────────────────────────────────────────────────────

meshOf :: String -> [Point] -> IO Mesh
meshOf label points = mapVertices (const ()) <$> pointMeshOf label points

pointMeshOf :: String -> [Point] -> IO PointMesh
pointMeshOf label points =
  buildTriangulation
    <$> requireRight ("build " <> label) (delaunay unitElementDefaults (V.fromList points))

pointsOf :: [(Double, Double)] -> [Point]
pointsOf keys = [Point x y | (x, y) <- keys]

dedupeAscending :: [(Double, Double)] -> [Point]
dedupeAscending sorted = [Point x y | (x, y) <- dropAdjacentDuplicates sorted]

dropAdjacentDuplicates :: Eq a => [a] -> [a]
dropAdjacentDuplicates (first : second : rest)
  | first == second = dropAdjacentDuplicates (second : rest)
  | otherwise = first : dropAdjacentDuplicates (second : rest)
dropAdjacentDuplicates rest = rest

-- ── observation ──────────────────────────────────────────────────────────────

siteKey :: Mesh -> VertexId -> (Double, Double)
siteKey mesh vertex = let Point x y = vertexPoint mesh vertex in (x, y)

siteList :: Mesh -> [(Double, Double)]
siteList mesh = [siteKey mesh vertex | vertex <- vertices mesh]

siteSet :: Mesh -> Set.Set (Double, Double)
siteSet = Set.fromList . siteList

edgeKeys :: Mesh -> Set.Set ((Double, Double), (Double, Double))
edgeKeys mesh =
  Set.fromList
    [ if left <= right then (left, right) else (right, left)
    | edge <- undirectedEdges mesh
    , let (from, to) = undirectedEndpoints mesh edge
    , let left = siteKey mesh from
    , let right = siteKey mesh to
    ]

-- | What a mesh looks like when two of them were supposed to be equal.
--
-- The interesting failure is the one where the coordinate-keyed edge sets
-- agree and the values do not: that is the join having become correct only up
-- to DCEL isomorphism, which is precisely what these laws exist to forbid, and
-- a report that only printed counts would hide it.
assertMesh :: String -> Mesh -> Mesh -> IO ()
assertMesh label expected actual
  | expected == actual = pure ()
  | otherwise = fail (label <> ": " <> report)
 where
  report
    | expectedEdges == actualEdges =
        "same geometry, different representation — "
          <> meshCounts expected
          <> " vs "
          <> meshCounts actual
    | otherwise =
        meshCounts expected
          <> " vs "
          <> meshCounts actual
          <> "; edges only in expected: "
          <> show (take 4 (Set.toList (Set.difference expectedEdges actualEdges)))
          <> "; only in actual: "
          <> show (take 4 (Set.toList (Set.difference actualEdges expectedEdges)))
  expectedEdges = edgeKeys expected
  actualEdges = edgeKeys actual

assertMeshEquivalent :: String -> Mesh -> Mesh -> IO ()
assertMeshEquivalent label expected actual = do
  canonicalExpected <-
    requireRight (label <> ": canonical expected") (canonicalize expected)
  canonicalActual <-
    requireRight (label <> ": canonical actual") (canonicalize actual)
  assertMesh label canonicalExpected canonicalActual

meshCounts :: Mesh -> String
meshCounts mesh =
  show (numVertices mesh)
    <> "v/"
    <> show (numUndirectedEdges mesh)
    <> "e/"
    <> show (numFaces mesh)
    <> "f"

-- ── combinatorics and pseudo-randomness ──────────────────────────────────────

pairs :: [a] -> [(a, a)]
pairs values = [(left, right) | left <- values, right <- values]

triples :: [a] -> [(a, a, a)]
triples values = [(a, b, c) | a <- values, b <- values, c <- values]

advance :: Word64 -> Word64
advance state = state * 6364136223846793005 + 1442695040888963407

randomWords :: Word64 -> Int -> [Word64]
randomWords seed count = take count (drop 1 (iterate advance seed))

-- | A deterministic permutation: decorate, sort by the key, discard it.
scramble :: Word64 -> [a] -> [a]
scramble seed values =
  map snd (sortBy (comparing fst) (zip (randomWords seed (length values)) values))

randomSites :: Word64 -> Int -> [Point]
randomSites seed count =
  [ Point (unitCoordinate first) (unitCoordinate second)
  | (first, second) <- take count (chunkPairs (randomWords seed (2 * count)))
  ]

chunkPairs :: [a] -> [(a, a)]
chunkPairs (first : second : rest) = (first, second) : chunkPairs rest
chunkPairs _ = []

unitCoordinate :: Word64 -> Double
unitCoordinate value = 2 * fromIntegral (value `div` 2048) / 9007199254740992 - 1