BulkLoad.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}

-- | Generation. @delaunay@ builds a mesh from a whole site set by circle sweep;
-- the insertion verbs extend an existing mesh one site at a time.
module Moonlight.Triangulation.BulkLoad
  ( empty
  , clear
  , delaunay
  , DuplicatePayloadPolicy (..)
  , delaunayFromCoordinates
  , insert
  , insertAt
  , insertMany
  ) where

import Control.Monad (forM_)
import Control.Monad.ST (ST, runST)
import qualified Data.IntSet as IntSet
import Data.Primitive.PrimArray
  ( MutablePrimArray
  , newPrimArray
  , unsafeFreezePrimArray
  , writePrimArray
  )
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as MUV
import Data.Word (Word32)
import Moonlight.Triangulation.Dcel (numInnerFaces, numVertices)
import Moonlight.Triangulation.Handles.HandleDefs (DirectedEdgeId (..), FaceId (..), VertexId (..))
import Moonlight.Triangulation.Internal.BoxedPaged (boxedFromVector, boxedUpdate, emptyBoxedPaged)
import Moonlight.Triangulation.Internal.Capacity (ensureCapacity)
import Moonlight.Triangulation.Insertion (insertExistingVertexAtLocation, insertVertexAtPoint)
import Moonlight.Triangulation.Internal.Location (MutableLocation (..))
import Moonlight.Triangulation.Internal.Mutable
import Moonlight.Triangulation.Internal.OperationState
  ( Counter (..)
  , OperationState
  , addCounter
  , freezeBuildStats
  , newOperationState
  , setCounter
  )
import Moonlight.Triangulation.Internal.CircleSweep (circleSweepInsert)
import Moonlight.Triangulation.Internal.PointIndex
  ( MutablePointIndex
  , emptyPointIndex
  , newMutablePointIndex
  , resolveMutablePoint
  , seedMutablePointIndex
  )
import Moonlight.Triangulation.Math (canonicalCoordinate, validatePoint)
import Moonlight.Triangulation.PointLocation (locatePointWithHint)
import Moonlight.Triangulation.Internal.Probe (Probe (..))
import Moonlight.Triangulation.Internal.Representation (Triangulation (..))
import Moonlight.Triangulation.Internal.PackedIndex (noIndex)
import Moonlight.Triangulation.Internal.Paged (TransactionShape (DenseTransaction, LocalTransaction), emptyPaged, fromVector)
import Moonlight.Triangulation.Internal.Transaction (runTransaction)
import Moonlight.Triangulation.Types

-- | The vertexless triangulation: the outer face and nothing else. This is the
-- canonical origin of the type — bulk loading, incremental insertion and
-- refinement all agree with growing this value.
empty
  :: ElementDefaults directed undirected face
  -> Triangulation mode vertex directed undirected face
empty defaults@ElementDefaults{defaultDirectedEdgeData, defaultUndirectedEdgeData, defaultFaceData} =
  Triangulation
    { triPointX = emptyPaged
    , triPointY = emptyPaged
    , triPointIndex = emptyPointIndex
    , triVertexOut = emptyPaged
    , triVertexData = emptyBoxedPaged Nothing
    , triHalfTopology = emptyPaged
    , triDirectedData = emptyBoxedPaged (Just defaultDirectedEdgeData)
    , triUndirectedData = emptyBoxedPaged (Just defaultUndirectedEdgeData)
    , triFaceEdge = fromVector noIndex (U.singleton noIndex)
    , triFaceData = boxedFromVector (Just defaultFaceData) (V.singleton defaultFaceData)
    , triConstraint = emptyPaged
    , triConstraintCount = 0
    , triConstraintEdges = IntSet.empty
    , triElementDefaults = defaults
    }

-- | Discard every vertex while retaining the element defaults the
-- triangulation was built with.
clear
  :: Triangulation mode vertex directed undirected face
  -> Triangulation mode vertex directed undirected face
clear = empty . triElementDefaults

-- | How a canonical bulk source combines payloads whose exact coordinates
-- coincide. Geometry identity is settled independently by the point index.
data DuplicatePayloadPolicy vertex
  = KeepFirstPayload
  | CombineDuplicatePayload !(vertex -> vertex -> vertex)

-- | Build a finite Delaunay DCEL while preserving the first input payload at
-- every duplicate position. The returned mapping relates every input slot to
-- the canonical stored vertex.
delaunay
  :: forall vertex directed undirected face
   . HasPosition vertex
  => ElementDefaults directed undirected face
  -> V.Vector vertex
  -> Either BuildError (BuildResult 'Unconstrained vertex directed undirected face)
delaunay defaults input = do
  validateVertices input
  buildDelaunayFromSource
    defaults
    (V.length input)
    (position . (input V.!))
    (input V.!)
    KeepFirstPayload

-- | Canonical construction from separate geometry and annotation sources.
-- The coordinate vector remains the only geometry in ingress; payloads never
-- acquire a fabricated 'HasPosition' instance merely to reach the loader.
delaunayFromCoordinates
  :: forall vertex directed undirected face
   . ElementDefaults directed undirected face
  -> V.Vector (Point)
  -> V.Vector vertex
  -> DuplicatePayloadPolicy vertex
  -> Either BuildError (BuildResult 'Unconstrained vertex directed undirected face)
delaunayFromCoordinates defaults coordinates payloads duplicatePolicy
  | coordinateCount /= payloadCount =
      Left (CoordinatePayloadCountMismatch coordinateCount payloadCount)
  | otherwise = do
      V.iforM_ coordinates (\index point -> validatePoint (Just index) point)
      buildDelaunayFromSource
        defaults
        coordinateCount
        (coordinates V.!)
        (payloads V.!)
        duplicatePolicy
 where
  !coordinateCount = V.length coordinates
  !payloadCount = V.length payloads

buildDelaunayFromSource
  :: forall vertex directed undirected face
   . ElementDefaults directed undirected face
  -> Int
  -> (Int -> Point)
  -> (Int -> vertex)
  -> DuplicatePayloadPolicy vertex
  -> Either BuildError (BuildResult 'Unconstrained vertex directed undirected face)
buildDelaunayFromSource defaults inputCount pointAtInput payloadAtInput duplicatePolicy = do
  ensureCapacity inputCount
  runST $ do
    mutable <- newMutableDcel defaults inputCount
    operation <- newOperationState (halfEdgeCapacity mutable)
    table <- newMutablePointIndex inputCount
    mapping <- newPrimArray inputCount
    ingressed <- ingress mutable operation table mapping 0 0 0
    case ingressed of
      Left failure -> pure (Left failure)
      Right (sumX, sumY) -> do
        unique <- pointCount mutable
        inserted <-
          if unique == 0
            then pure (Right 0)
            else do
              let !scale = recip (fromIntegral unique)
              arena <-
                fillRadialArena
                  mutable
                  (sumX * scale)
                  (sumY * scale)
                  (\index -> pure (fromIntegral index))
                  unique
              circleSweepInsert mutable operation arena
        case inserted of
          Left failure -> pure (Left failure)
          Right seedCount -> do
            setCounter operation CounterSpatialSeedPoints seedCount
            frozenOutcome <- freezeTriangulation mutable
            case frozenOutcome of
              Left failure -> pure (Left failure)
              Right frozen -> do
                mapped <- unsafeFreezePrimArray mapping
                stats <- freezeBuildStats operation
                pure
                  ( Right
                      BuildResult
                        { buildTriangulation = frozen
                        , buildInputVertices = mapped
                        , buildStats = stats
                        }
                  )
 where
  -- One indexed ingress loop: read the position once, claim it against the
  -- transient table, write the input mapping, and accumulate the sort centre
  -- over the vertices that are actually new. No list, no decorated vector.
  ingress
    :: forall s
     . MutableDcel s vertex directed undirected face
    -> OperationState s
    -> MutablePointIndex s
    -> MutablePrimArray s Word32
    -> Int
    -> Double
    -> Double
    -> ST s (Either BuildError (Double, Double))
  ingress mutable operation table mapping !index !sumX !sumY
    | index >= inputCount = pure (Right (sumX, sumY))
    | otherwise = do
        addCounter operation CounterInputPoints 1
        let !vertexData = payloadAtInput index
        claimed <- claimPosition mutable table (pointAtInput index) vertexData
        case claimed of
          Left failure -> pure (Left failure)
          Right (vertex, fresh) -> do
            writePrimArray mapping index (fromIntegral vertex)
            if fresh
              then do
                addCounter operation CounterUniquePoints 1
                x <- readPointX mutable vertex
                y <- readPointY mutable vertex
                ingress mutable operation table mapping (index + 1) (sumX + x) (sumY + y)
              else do
                case duplicatePolicy of
                  KeepFirstPayload -> pure ()
                  CombineDuplicatePayload combine -> do
                    resident <- vertexDataAt mutable vertex
                    writeVertexData mutable vertex (combine resident vertexData)
                addCounter operation CounterDuplicatePoints 1
                ingress mutable operation table mapping (index + 1) sumX sumY

-- | Insert or replace a vertex payload. A payload at an existing position is
-- overwritten without changing topology.
--
-- The published mesh is independent of the one passed in. A singleton below
-- ten thousand resident sites copies densely; larger bases publish through
-- copy-on-write pages. A caller inserting a sequence wants one
-- 'Moonlight.Triangulation.Session.withSession' over
-- 'Moonlight.Triangulation.Session.insertVertex' instead — see 'insertAt'.
insert
  :: HasPosition vertex
  => Triangulation mode vertex directed undirected face
  -> vertex
  -> Either BuildError (InsertionResult mode vertex directed undirected face)
insert triangulation vertexData = insertAt triangulation (position vertexData) vertexData

-- | Insert at a stated point. 'insert' is this with the point read out of the
-- payload, which is what a caller holding only a payload wants; a caller that
-- computed the point — a constraint split, a Steiner refinement — wants to say
-- so rather than build a payload and hope the round trip through 'HasPosition'
-- returns what it started with.
--
-- This is one shaped transaction over a single insertion. Replacing a fold of
-- it with one session is sound because the two agree on every mesh and differ
-- only in how many intermediate meshes they publish. A fold publishes @k@
-- meshes and pays a thaw for each, so it runs in Θ(n·k); the session pays one
-- thaw and runs in O(k·log n) expected.
insertAt
  :: Triangulation mode vertex directed undirected face
  -> Point
  -> vertex
  -> Either BuildError (InsertionResult mode vertex directed undirected face)
insertAt triangulation rawPoint vertexData = do
  queryPoint <- validatePoint Nothing rawPoint
  case locatePointWithHint triangulation Nothing queryPoint of
    (OnVertex resident, walked) ->
      Right (replaceResidentPayload triangulation resident walked vertexData)
    (located, walked) -> do
      let transactionShape =
            if numVertices triangulation < 10_000
              then DenseTransaction
              else LocalTransaction
      ((vertex, disposition), frozen, stats) <-
        runTransaction
          id
          transactionShape
          triangulation
          1
          (\mutable operation -> do
             addCounter operation CounterInputPoints 1
             inserted <-
               case located of
                 -- A frozen degenerate-line location points at an arbitrary visible
                 -- segment, while the line extension interpreter requires a terminal
                 -- edge. The frozen section carries no terminal witness, so retain
                 -- the existing mutable line locator for this one non-lawful case.
                 OutsideConvexHull (Just _)
                   | numInnerFaces triangulation == 0 ->
                       insertVertexAtPoint @'ProbeOff mutable operation Nothing (queryPointValue queryPoint) vertexData
                 _ -> do
                   capacityOutcome <- ensurePointCapacity mutable 1
                   case capacityOutcome of
                     Left failure -> pure (Left failure)
                     Right () -> do
                       vertex <- appendVertex mutable (queryPointValue queryPoint) vertexData
                       let thawedLocation =
                             case located of
                               EmptyTriangulation -> MutableEmpty
                               OnEdge (DirectedEdgeId raw) -> MutableOnEdge (fromIntegral raw)
                               InFace (FaceId raw) -> MutableInFace (fromIntegral raw)
                               -- The frozen locator emits no edge only for a
                               -- singleton mesh. Its mutable interpreter ignores
                               -- this sentinel while constructing the second vertex.
                               OutsideConvexHull Nothing -> MutableOutsideHull 0
                               OutsideConvexHull (Just (DirectedEdgeId raw)) -> MutableOutsideHull (fromIntegral raw)
                       ((vertex, Inserted) <$) <$> insertExistingVertexAtLocation @'ProbeOff mutable operation vertex thawedLocation
             case inserted of
               Left failure -> pure (Left failure)
               Right (vertex, disposition) -> do
                 case disposition of
                   Inserted -> addCounter operation CounterUniquePoints 1
                   AlreadyPresent -> do
                     writeVertexData mutable vertex vertexData
                     addCounter operation CounterExistingPoints 1
                     addCounter operation CounterDuplicatePoints 1
                 pure (Right (vertex, disposition))
          )
      pure
        InsertionResult
          { insertionTriangulation = frozen
          , insertionVertex = VertexId (fromIntegral vertex)
          , insertionDisposition = disposition
          , insertionStats = withFrozenLocationStats walked stats
          }

-- | Publish a payload replacement without opening a transaction.
--
-- A position already resident changes exactly one thing: the payload slot the
-- vertex already occupies. No coordinate, no half-edge, no constraint flag and
-- no face record differs, so the five unboxed planes are the ones the argument
-- already holds rather than copies taken out of it — which is what a thaw costs
-- and what this exists to refuse. They are immutable values; nothing reached
-- from here is a mutable buffer, and 'boxedUpdate' materializes a fresh page
-- for the one it rewrites, leaving the argument's own directory intact.
--
-- The location counters are the frozen walk's, not a thawed walk's. They
-- describe the walk that actually ran.
replaceResidentPayload
  :: Triangulation mode vertex directed undirected face
  -> VertexId
  -> LocationStats
  -> vertex
  -> InsertionResult mode vertex directed undirected face
replaceResidentPayload triangulation resident@(VertexId raw) walked vertexData =
  InsertionResult
    { insertionTriangulation =
        triangulation
          { triVertexData =
              boxedUpdate (fromIntegral raw) vertexData (triVertexData triangulation)
          }
    , insertionVertex = resident
    , insertionDisposition = AlreadyPresent
    , insertionStats =
        withFrozenLocationStats
          walked
          emptyBuildStats
            { statInputPoints = 1
            , statExistingPoints = 1
            , statDuplicatePoints = 1
            }
    }

-- | Add the frozen locator's observation to the local topology interpreter's
-- operation-owned counters. Direct frozen-site insertion contributes no mutable
-- walk; the degenerate fallback contributes its real mutable walk rather than
-- having it erased from the published result.
withFrozenLocationStats :: LocationStats -> BuildStats -> BuildStats
withFrozenLocationStats walked stats =
  stats
    { statLocationWalkSteps = locationWalkSteps walked + statLocationWalkSteps stats
    , statLocationMaxWalk = max (locationWalkSteps walked) (statLocationMaxWalk stats)
    , statLocationFallbacks = (if locationUsedFallback walked then 1 else 0) + statLocationFallbacks stats
    }

-- | Apply a batch in one page transaction. Duplicate positions are processed
-- in input order, so their last payload wins exactly as repeated 'insert'
-- calls would, while topology is inserted only once per new position.
insertMany
  :: forall mode vertex directed undirected face
   . HasPosition vertex
  => Triangulation mode vertex directed undirected face
  -> V.Vector vertex
  -> Either BuildError (BuildResult mode vertex directed undirected face)
insertMany triangulation input = do
  validateVertices input
  ensureCapacity (numVertices triangulation + V.length input)
  runST $ do
    mutable <-
      thawTriangulationDense
        (numVertices triangulation + V.length input)
        triangulation
    operation <- newOperationState (halfEdgeCapacity mutable)
    table <- newMutablePointIndex (numVertices triangulation + V.length input)
    seeded <- seedPointTable mutable table
    case seeded of
      Left failure -> pure (Left failure)
      Right () -> do
        mapping <- newPrimArray (V.length input)
        freshBuffer <- MUV.new (V.length input)
        filled <- fill mutable operation table mapping freshBuffer 0 0 0 0
        case filled of
          Left failure -> pure (Left failure)
          Right (sumX, sumY, freshCount) -> do
            inserted <-
              if freshCount == 0
                then pure (Right 0)
                else do
                  let !scale = recip (fromIntegral freshCount)
                  arena <-
                    fillRadialArena
                      mutable
                      (sumX * scale)
                      (sumY * scale)
                      (MUV.unsafeRead freshBuffer)
                      freshCount
                  circleSweepInsert mutable operation arena
            case inserted of
              Left failure -> pure (Left failure)
              Right seedCount -> do
                setCounter operation CounterSpatialSeedPoints seedCount
                frozenOutcome <- freezeTriangulation mutable
                case frozenOutcome of
                  Left failure -> pure (Left failure)
                  Right frozen -> do
                    mapped <- unsafeFreezePrimArray mapping
                    stats <- freezeBuildStats operation
                    pure
                      ( Right
                          BuildResult
                            { buildTriangulation = frozen
                            , buildInputVertices = mapped
                            , buildStats = stats
                            }
                      )
 where
  fill
    :: forall s
     . MutableDcel s vertex directed undirected face
    -> OperationState s
    -> MutablePointIndex s
    -> MutablePrimArray s Word32
    -> MUV.MVector s Word32
    -> Int
    -> Double
    -> Double
    -> Int
    -> ST s (Either BuildError (Double, Double, Int))
  fill mutable operation table mapping freshBuffer !index !sumX !sumY !freshCount
    | index >= V.length input = pure (Right (sumX, sumY, freshCount))
    | otherwise = do
        addCounter operation CounterInputPoints 1
        let !vertexData = input V.! index
        claimed <- claimPosition mutable table (position vertexData) vertexData
        case claimed of
          Left failure -> pure (Left failure)
          Right (vertex, fresh) -> do
            writePrimArray mapping index (fromIntegral vertex)
            if fresh
              then do
                addCounter operation CounterUniquePoints 1
                MUV.unsafeWrite freshBuffer freshCount (fromIntegral vertex)
                x <- readPointX mutable vertex
                y <- readPointY mutable vertex
                fill
                  mutable
                  operation
                  table
                  mapping
                  freshBuffer
                  (index + 1)
                  (sumX + x)
                  (sumY + y)
                  (freshCount + 1)
              else do
                writeVertexData mutable vertex vertexData
                addCounter operation CounterExistingPoints 1
                addCounter operation CounterDuplicatePoints 1
                fill mutable operation table mapping freshBuffer (index + 1) sumX sumY freshCount

-- | One packed radial record per swept vertex — the derived sort fields and
-- the vertex handle, nothing else — filled straight from the coordinate
-- arenas and consumed in place by the sweep. The squared distance is stated
-- against the ingress-accumulated centre, in the widened comparison format.
fillRadialArena
  :: MutableDcel s vertex directed undirected face
  -> Double
  -> Double
  -> (Int -> ST s Word32)
  -> Int
  -> ST s (MUV.MVector s (Double, Double, Double, Word32))
fillRadialArena mutable centerX centerY lookupId count = do
  arena <- MUV.new count
  forM_ [0 .. count - 1] $ \index -> do
    raw <- lookupId index
    x <- readPointX mutable (fromIntegral raw)
    y <- readPointY mutable (fromIntegral raw)
    let !wideX = x
        !wideY = y
        !deltaX = centerX - wideX
        !deltaY = centerY - wideY
    MUV.unsafeWrite arena index (deltaX * deltaX + deltaY * deltaY, wideX, wideY, raw)
  pure arena

-- | Claim a position for the vertex the arena would append next, or answer the
-- vertex already holding it. The claim is written before the append, so the two
-- must stay adjacent: nothing may consume a vertex slot in between.
claimPosition
  :: MutableDcel s vertex directed undirected face
  -> MutablePointIndex s
  -> Point
  -> vertex
  -> ST s (Either BuildError (Int, Bool))
claimPosition mutable table rawPoint vertexData =
  case rawPoint of
    Point x y -> do
      let !canonicalX = canonicalCoordinate x
          !canonicalY = canonicalCoordinate y
      candidate <- pointCount mutable
      owner <-
        resolveMutablePoint
          table
          (readPointX mutable)
          (readPointY mutable)
          canonicalX
          canonicalY
          candidate
      case owner of
        Left failure -> pure (Left failure)
        Right (Just existing) -> pure (Right (existing, False))
        Right Nothing -> do
          vertex <- appendVertexCoordinates mutable canonicalX canonicalY vertexData
          pure (Right (vertex, True))
-- | Enter the positions a batch inherits from the triangulation it extends, so
-- that an input repeating one of them maps to the vertex already there.
seedPointTable
  :: MutableDcel s vertex directed undirected face
  -> MutablePointIndex s
  -> ST s (Either BuildError ())
seedPointTable mutable table = do
  existing <- pointCount mutable
  seedMutablePointIndex
    table
    existing
    (readPointX mutable)
    (readPointY mutable)

validateVertices :: HasPosition vertex => V.Vector vertex -> Either BuildError ()
validateVertices vertices =
  V.iforM_ vertices (\index vertexData -> validatePoint (Just index) (position vertexData))