Support.hs

{-# LANGUAGE DataKinds #-}

-- | Input generation, canonical encoding, and the shared build helpers.
module Moonlight.Planar.Bench.SpadeCompare.Support where

import Data.List (sort)
import Data.Word (Word64)
import GHC.Float (castDoubleToWord64)
import Numeric (showHex)
import qualified Data.Vector as V
import Moonlight.Planar.BulkLoad (delaunay)
import Moonlight.Planar.Cdt (constraintEdges)
import Moonlight.Planar.Dcel (undirectedEndpoints, vertexPoint)
import Moonlight.Planar.Handles.Iterators.FixedIterators (undirectedEdges)
import Moonlight.Planar.Point (mkQueryPoint)
import Moonlight.Planar.Math (orient2d)
import Moonlight.Planar.Types (ConstrainedDelaunayTriangulation, DelaunayTriangulation, Triangulation, buildTriangulation, unitElementDefaults)
import Moonlight.Planar.Point (HasPosition (..), Point (..), QueryPoint)
import Moonlight.Planar.Telemetry (LocationHint (..), nearestNeighbor)
import Moonlight.Planar.Bench.SpadeCompare.Input (randomPoints)
import Moonlight.Planar.HintGenerator

-- | The same generated points flattened onto a sliver a picometre thick. No
-- coordinate becomes equal to another and no four points become cocircular, so
-- the input stays in general position and the load stays on the ordinary
-- face-building path; what collapses is the angular spread the bulk loaders
-- index their hulls by. That collapse is the whole of what this input
-- measures, and it does not reach the face-less location path — that one needs
-- exact collinearity, which 'exactlyCollinearPoints' supplies and gates.
nearCollinearPoints :: Word64 -> Int -> [Point]
nearCollinearPoints seed count =
  [Point x (y * 1.0e-12) | Point x y <- randomPoints seed count]

-- | Points on one exact line: @(i, 0)@, every coordinate an integer that a
-- 'Double' holds exactly. No three are in general position, so 'orient2d'
-- returns an exact zero on both sides and no face is ever built.
--
-- Uniqueness does not need general position, only one valid answer. A
-- collinear set admits exactly one triangulation — the chain of consecutive
-- segments, @n - 1@ of them and no triangle — so the canonical edge set is
-- fully determined and gates bit for bit against the strict side. Both
-- implementations state that same count as their own degenerate Euler
-- invariant, which is what makes this a shared answer rather than a
-- coincidence.
exactlyCollinearPoints :: Int -> [Point]
exactlyCollinearPoints count = [Point (fromIntegral index) 0 | index <- [0 .. count - 1]]

hex64 :: Double -> String
hex64 value =
  let raw = showHex (castDoubleToWord64 value) ""
   in replicate (16 - length raw) '0' <> raw

-- | Canonical edge set: each undirected edge as its two endpoint coordinates in
-- IEEE-754 bit patterns, endpoint-ordered then globally sorted. Formatting
-- cannot introduce a spurious disagreement and rounding cannot hide a real one.
canonicalEdges :: Triangulation mode vertex directed undirected face -> [String]
canonicalEdges triangulation =
  sort
    [ encodeCanonicalEdge first second
    | edge <- undirectedEdges triangulation
    , let (from, to) = undirectedEndpoints triangulation edge
          p = vertexPoint triangulation from
          q = vertexPoint triangulation to
          (first, second) = if key p <= key q then (p, q) else (q, p)
    ]
 where
  key (Point x y) = (x, y)

canonicalConstraintEdges
  :: ConstrainedDelaunayTriangulation vertex
  -> [String]
canonicalConstraintEdges triangulation =
  sort
    [ encodeCanonicalEdge first second
    | edge <- constraintEdges triangulation
    , let (from, to) = undirectedEndpoints triangulation edge
          p = vertexPoint triangulation from
          q = vertexPoint triangulation to
          (first, second) =
            if pointKey p <= pointKey q
              then (p, q)
              else (q, p)
    ]
 where
  pointKey (Point x y) = (x, y)

encodeCanonicalEdge :: Point -> Point -> String
encodeCanonicalEdge (Point ax ay) (Point bx by) =
  concatMap hex64 [ax, ay, bx, by]

coordHex :: Point -> String
coordHex (Point x y) = hex64 x <> hex64 y

-- | A cycle, rotated to begin at its smallest encoded element.
--
-- Both implementations keep one entry edge per vertex and per face — whichever
-- their builder wrote last — so the two walk the same cycle from different
-- places. Rotating pins the order, which the mesh determines, and declines to
-- pin the starting point, which it does not. The rotation compares the encoded
-- strings rather than the coordinates so that both sides order by the same
-- comparison: the bit pattern of a negative double does not sort in numeric
-- order, so a Haskell 'Ord' on 'Point' and a Rust @partial_cmp@ on a tuple
-- would each have to be argued to agree where two byte strings simply do.
rotateToSmallest :: [String] -> [String]
rotateToSmallest [] = []
rotateToSmallest ring =
  let !smallest = minimum ring
      (before, after) = break (== smallest) ring
   in after <> before

require :: Show failure => Either failure value -> IO value
require (Left failure) = fail (show failure)
require (Right value) = pure value

requireQueryPoints
  :: Traversable container
  => container Point
  -> IO (container QueryPoint)
requireQueryPoints = require . traverse mkQueryPoint

requireQueryChords
  :: Traversable container
  => container (Point, Point)
  -> IO (container (QueryPoint, QueryPoint))
requireQueryChords =
  require
    . traverse
      (\(from, to) -> (,) <$> mkQueryPoint from <*> mkQueryPoint to)

delaunayOf :: Int -> IO (DelaunayTriangulation Point)
delaunayOf count =
  buildTriangulation <$> require (delaunay unitElementDefaults (V.fromList (randomPoints 0x9e3779b97f4a7c15 count)))

-- | One line per query: the coordinates of the nearest vertex the hierarchy
-- routed the search to, or @none@.
nearestAnswer
  :: DelaunayTriangulation Point
  -> HierarchyHint
  -> QueryPoint
  -> String
nearestAnswer triangulation hierarchy query =
  let hint = case hierarchyHint hierarchy query of
        Just (VertexHint vertex) -> Just vertex
        _ -> Nothing
   in case nearestNeighbor triangulation hint query of
        Nothing -> "none"
        Just (vertex, _) ->
          let Point x y = vertexPoint triangulation vertex
           in hex64 x <> hex64 y

pointHex :: Point -> String
pointHex (Point x y) = hex64 x <> "," <> hex64 y

constraintPairs :: Int -> Int -> [(Int, Int)]
constraintPairs pointCount constraintCount =
  take constraintCount
    [ (a, b)
    | index <- [0 ..]
    , let a = index `mod` pointCount
          b = (index * 6151 + pointCount `quot` 2) `mod` pointCount
    , a /= b
    ]

-- | A closed square with no dangling constraint. The interior segment used by
-- the package's own benchmark is outside spade's documented "closed shape"
-- contract and makes it refuse the whole domain — not a workload the two can be
-- timed on.
refinementInput :: (V.Vector Point, V.Vector (Int, Int))
refinementInput =
  ( V.fromList [Point 0 0, Point 64 0, Point 64 64, Point 0 64]
  , V.fromList [(0, 1), (1, 2), (2, 3), (3, 0)]
  )