PowerDiagram.hs

{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveAnyClass #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DerivingStrategies #-}

-- | Exact regular topology and power cells. Labelled cell dispositions are
-- authoritative; regular topology and planar layers are derived views.
module Moonlight.Triangulation.PowerDiagram
  ( PowerWeight
  , PowerWeightError (..)
  , powerWeight
  , powerWeightExact
  , PowerSite
  , powerSite
  , powerSiteLabel
  , powerSitePosition
  , powerSiteWeight
  , PowerCellDisposition (..)
  , RegularSiteDisposition (..)
  , RegularFace
  , regularFaceLabels
  , regularFaceDualPoint
  , PowerDualEdge (..)
  , RegularEdge
  , regularEdgeLabels
  , regularEdgeDual
  , RegularTriangulation
  , regularTriangulation
  , regularSiteDisposition
  , regularFaces
  , regularEdges
  , regularNeighbours
  , RegularTriangulationReceipt (..)
  , BoundedPowerDiagram
  , boundedPowerDiagram
  , powerCellDisposition
  , powerCellDispositions
  , powerDiagramPlanarLayer
  , PowerDiagramError (..)
  , RegularTopologyError (..)
  , PowerDiagramReceipt (..)
  , powerDiagramInputSites
  , powerDiagramPeakIntermediateBitGrowth
  , powerDiagramFinalCoordinateBitGrowth
  , AffineForm (..)
  , UpperEnvelopeError (..)
  , upperEnvelope
  ) where

import Control.DeepSeq (NFData)
import Data.Bifunctor (first)
import qualified Data.Foldable as Foldable
import qualified Data.List as List
import Data.List.NonEmpty (NonEmpty (..))
import qualified Data.List.NonEmpty as NonEmpty
import qualified Data.Map.Strict as Map
import Data.Map.Strict (Map)
import Data.Ord (comparing)
import Data.Set (Set)
import qualified Data.Set as Set
import qualified Data.Vector as Vector
import GHC.Generics (Generic)
import Moonlight.Triangulation.Exact
  ( ExactClipDisposition (..)
  , ExactClipError
  , ExactClipReceipt (..)
  , ExactAffineLine
  , ExactClosedHalfPlane
  , ExactHalfPlaneError
  , ExactPoint
  , ExactRay
  , ExactRetainedPolygon
  , ExactSegment
  , exactClipRetainedPolygon
  , exactClosedHalfPlane
  , exactPointCoordinates
  , exactPointBitWidth
  , exactPointFromQueryPoint
  , exactRetainedPolygon
  , oppositeExactAffineLine
  )
import Moonlight.Triangulation.Internal.ExactRational
  ( ExactRational
  , exactRationalBitWidth
  , exactRationalFromFiniteDouble
  )
import Moonlight.Triangulation.Internal.BoundaryCycle
  ( consecutivePairs )
import Moonlight.Triangulation.Internal.RegularTriangulation
  ( DistinctSlopeGenerators (..)
  , ExactPowerGenerator (..)
  , GeneratorDualGeometry (..)
  , GeneratorRegularEdge (..)
  , GeneratorRegularFace (..)
  , GeneratorRegularReceipt (..)
  , GeneratorRegularTopology (..)
  , RegularGeneratorDisposition (..)
  , RegularTopologyError (..)
  , exactGeneratorAxis
  , regularGeneratorTopology
  )
import Moonlight.Triangulation.Internal.Minkowski.Convex
  ( convexHullPolygon
  , convexPolygonComponent
  , convexPolygonFromRetained
  )
import Moonlight.Triangulation.Internal.Overlay.Types
  ( OverlayCell (..)
  , OverlayCellGeometry (..)
  , OverlayError
  , OverlayResult (..)
  )
import Moonlight.Triangulation.Internal.Region.Publication
  ( planarLayerFromAdmittedComponents
  )
import Moonlight.Triangulation.Math (mkQueryPoint)
import Moonlight.Triangulation.Minkowski
  ( ConvexPolygon
  , MinkowskiError
  , convexPolygonPoints
  )
import Moonlight.Triangulation.Overlay
  ( overlayLayers )
import Moonlight.Triangulation.Region
  ( PlanarLayer
  , PolygonComponent
  , exactLoopPoints
  , polygonOuterLoop
  )
import Moonlight.Triangulation.Types
  ( NonFiniteValue
  , Point
  , PointValidationError
  , QueryPoint
  , classifyNonFinite
  , queryPointValue
  )

-- | An admitted signed additive power offset.  Power distance is
-- @||x-p||^2-w@, so negative values are lawful and this is deliberately not a
-- squared-radius refinement.
newtype PowerWeight = PowerWeight ExactRational
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | The sole obstruction to admitting a signed binary64 power offset.
data PowerWeightError
  = PowerWeightNonFinite !NonFiniteValue
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | Admit a finite binary64 power offset exactly.
powerWeight :: Double -> Either PowerWeightError PowerWeight
powerWeight value =
  case classifyNonFinite value of
    Just obstruction -> Left (PowerWeightNonFinite obstruction)
    Nothing -> Right (PowerWeight (exactRationalFromFiniteDouble value))

-- | Exact rational value of an admitted power offset.
powerWeightExact :: PowerWeight -> ExactRational
powerWeightExact (PowerWeight value) = value

-- | One labelled, admitted weighted site.  Construction validates and
-- canonicalizes the binary64 position once.
data PowerSite label = PowerSite
  { powerSiteLabel :: !label
  , powerSiteQueryPoint :: !QueryPoint
  , powerSiteWeight :: !PowerWeight
  }
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

powerSite
  :: label
  -> Point
  -> PowerWeight
  -> Either (PowerDiagramError label) (PowerSite label)
powerSite label point weight =
  PowerSite label <$> first (PowerSitePositionInvalid label) (mkQueryPoint point) <*> pure weight

powerSitePosition :: PowerSite label -> Point
powerSitePosition = queryPointValue . powerSiteQueryPoint

-- | Exactly one authoritative result for each submitted label.
data PowerCellDisposition label
  = PublishedPowerCell !ConvexPolygon
  | LowerDimensionalPowerCell !(NonEmpty ExactPoint)
  | EmptyPowerCell
  | CoincidentEquivalentTo !label
  | CoincidentDominatedBy !label
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | Visibility of one submitted label in the exact regular subdivision.
data RegularSiteDisposition label
  = RegularSiteVisible
  | RegularSiteLowerDimensional
  | RegularSiteHidden
  | RegularSiteCoincidentEquivalentTo !label
  | RegularSiteCoincidentDominatedBy !label
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | One oriented regular face and its exact weighted-dual vertex.
data RegularFace label = RegularFace !label !label !label !ExactPoint
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

regularFaceLabels :: RegularFace label -> (label, label, label)
regularFaceLabels (RegularFace firstLabel secondLabel thirdLabel _) =
  (firstLabel, secondLabel, thirdLabel)

regularFaceDualPoint :: RegularFace label -> ExactPoint
regularFaceDualPoint (RegularFace _ _ _ dualPoint) = dualPoint

-- | Exact weighted Voronoi geometry dual to one regular edge.
data PowerDualEdge
  = BoundedPowerDual !ExactSegment
  | UnboundedPowerDual !ExactRay
  | FullLinePowerDual !ExactAffineLine
  | CollapsedPowerDual !ExactPoint
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | One unordered regular edge and its exact dual geometry.
data RegularEdge label = RegularEdge !label !label !PowerDualEdge
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

regularEdgeLabels :: RegularEdge label -> (label, label)
regularEdgeLabels (RegularEdge firstLabel secondLabel _) =
  (firstLabel, secondLabel)

regularEdgeDual :: RegularEdge label -> PowerDualEdge
regularEdgeDual (RegularEdge _ _ dual) = dual

-- | Exact regular subdivision, opaque so incidence and visibility cannot
-- disagree.
data RegularTriangulation label = RegularTriangulation
  { storedRegularDispositions :: !(Map label (RegularSiteDisposition label))
  , storedRegularFaces :: ![RegularFace label]
  , storedRegularEdges :: ![RegularEdge label]
  , storedRegularNeighbours :: !(Map label (Set label))
  }
  deriving stock (Eq, Show, Generic)
  deriving anyclass (NFData)

data RegularTriangulationReceipt = RegularTriangulationReceipt
  { regularTriangulationInputSites :: !Int
  , regularTriangulationRepresentativeSites :: !Int
  , regularTriangulationVisibleSites :: !Int
  , regularTriangulationLowerDimensionalSites :: !Int
  , regularTriangulationHiddenSites :: !Int
  , regularTriangulationCoincidentSites :: !Int
  , regularTriangulationFaces :: !Int
  , regularTriangulationEdges :: !Int
  , regularTriangulationPeakHullFaces :: !Int
  }
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | Total labelled result, opaque so callers cannot omit a submitted label.
newtype BoundedPowerDiagram label =
  BoundedPowerDiagram (Map label (PowerCellDisposition label))
  deriving stock (Eq, Show, Generic)
  deriving anyclass (NFData)

data PowerDiagramError label
  = PowerSitePositionInvalid !label !PointValidationError
  | DuplicatePowerSiteLabel !label
  | PowerDomainInvalid !ExactHalfPlaneError
  | PowerBisectorInvalid !label !label !ExactHalfPlaneError
  | PowerRegularTopologyFailed !(RegularTopologyError label)
  | PowerCellClipFailed !label !ExactClipError
  | PowerDiagramOutsideLabelCollides !label
  deriving stock (Eq, Show, Generic)
  deriving anyclass (NFData)

data PowerDiagramReceipt = PowerDiagramReceipt
  { powerDiagramDomainVertices :: !Int
  , powerDiagramSubmittedSiteConstraints :: !Int
  , powerDiagramActiveBoundaries :: !Int
  , powerDiagramBoundaryCompatibilityChecks :: !Int
  , powerDiagramExactIntersections :: !Int
  , powerDiagramPublishedCells :: !Int
  , powerDiagramLowerDimensionalCells :: !Int
  , powerDiagramEmptyCells :: !Int
  , powerDiagramCoincidentEquivalentCells :: !Int
  , powerDiagramCoincidentDominatedCells :: !Int
  , powerDiagramRegularFaces :: !Int
  , powerDiagramRegularEdges :: !Int
  , powerDiagramOracleCells :: !Int
  , powerDiagramMaximumCellConstraints :: !Int
  , powerDiagramMaximumInputBits :: !Int
  , powerDiagramMaximumAffineCoefficientBits :: !Int
  , powerDiagramPeakIntermediateCoordinateBits :: !Int
  , powerDiagramFinalCoordinateBits :: !Int
  , powerDiagramFinalDenominatorBits :: !Int
  }
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | Submitted sites, derived from the exhaustive disposition partition.
powerDiagramInputSites :: PowerDiagramReceipt -> Int
powerDiagramInputSites receipt =
  powerDiagramPublishedCells receipt
    + powerDiagramLowerDimensionalCells receipt
    + powerDiagramEmptyCells receipt
    + powerDiagramCoincidentEquivalentCells receipt
    + powerDiagramCoincidentDominatedCells receipt
{-# INLINE powerDiagramInputSites #-}

-- | Peak exact-coordinate width beyond the widest admitted input.
powerDiagramPeakIntermediateBitGrowth :: PowerDiagramReceipt -> Int
powerDiagramPeakIntermediateBitGrowth receipt =
  max
    0
    ( powerDiagramPeakIntermediateCoordinateBits receipt
        - powerDiagramMaximumInputBits receipt
    )
{-# INLINE powerDiagramPeakIntermediateBitGrowth #-}

-- | Published coordinate width beyond the widest admitted input.
powerDiagramFinalCoordinateBitGrowth :: PowerDiagramReceipt -> Int
powerDiagramFinalCoordinateBitGrowth receipt =
  max
    0
    ( powerDiagramFinalCoordinateBits receipt
        - powerDiagramMaximumInputBits receipt
    )
{-# INLINE powerDiagramFinalCoordinateBitGrowth #-}

-- | One exact affine form @c0 + cx*x + cy*y@.
data AffineForm = AffineForm
  { affineFormConstant :: !ExactRational
  , affineFormXCoefficient :: !ExactRational
  , affineFormYCoefficient :: !ExactRational
  }
  deriving stock (Eq, Ord, Show, Generic)
  deriving anyclass (NFData)

-- | Typed obstructions from exact affine argmax decomposition.  Power-cell
-- construction remains the canonical geometric owner; window restriction is
-- the only additional boundary.
data UpperEnvelopeError label
  = UpperEnvelopeEmptyForms
  | UpperEnvelopeWindowHullFailed !MinkowskiError
  | UpperEnvelopePowerConstructionFailed !(PowerDiagramError label)
  | UpperEnvelopeWindowOverlayFailed !(OverlayError (Maybe label) Bool)
  deriving stock (Eq, Show, Generic)
  deriving anyclass (NFData)

-- | Construct the exact regular subdivision before any bounded clipping.
regularTriangulation
  :: Ord label
  => NonEmpty (PowerSite label)
  -> Either
      (PowerDiagramError label)
      (RegularTriangulation label, RegularTriangulationReceipt)
regularTriangulation submitted = do
  sortedSites <- validateAndSortSites submitted
  let generators = fmap (fst . prepareExactPowerGenerator) sortedSites
  section <- resolvedGeneratorSection generators
  pure (publishRegularSection section)

regularSiteDisposition
  :: Ord label
  => label
  -> RegularTriangulation label
  -> Maybe (RegularSiteDisposition label)
regularSiteDisposition label = Map.lookup label . storedRegularDispositions

regularFaces :: RegularTriangulation label -> [RegularFace label]
regularFaces = storedRegularFaces

regularEdges :: RegularTriangulation label -> [RegularEdge label]
regularEdges = storedRegularEdges

regularNeighbours
  :: Ord label
  => label
  -> RegularTriangulation label
  -> Set label
regularNeighbours label =
  Map.findWithDefault Set.empty label . storedRegularNeighbours

boundedPowerDiagram
  :: Ord label
  => ConvexPolygon
  -> NonEmpty (PowerSite label)
  -> Either (PowerDiagramError label) (BoundedPowerDiagram label, PowerDiagramReceipt)
boundedPowerDiagram domain submitted = do
  sortedSites <- validateAndSortSites submitted
  let preparedGenerators = fmap prepareExactPowerGenerator sortedSites
      generators = fmap fst preparedGenerators
      siteInputBits =
        Foldable.foldl' (\bits preparation -> max bits (snd preparation)) 0 preparedGenerators
  (dispositions, clipReceipt, regularReceipt, maximumCellConstraints) <-
    exactGeneratorDispositionsWith (<>) mempty domain generators
  let receipt =
        aggregateReceipt
          domain
          dispositions
          siteInputBits
          clipReceipt
          regularReceipt
          maximumCellConstraints
  pure (BoundedPowerDiagram dispositions, receipt)

-- | Full-dimensional labelled argmax regions for affine forms
-- @c0 + cx*x + cy*y@ inside an admitted polygonal window. The result is a
-- planar projection: lower-dimensional and empty winners intentionally have
-- no region. Use 'regularTriangulation' on corresponding weighted sites when
-- those dispositions or exact unbounded duals are required. Identical forms
-- choose the least label, independent of map construction order.
upperEnvelope
  :: Ord label
  => PolygonComponent
  -> Map label AffineForm
  -> Either (UpperEnvelopeError label) (PlanarLayer (Maybe label))
upperEnvelope window forms = do
  generators <- affineFormGenerators forms
  domain <-
    first UpperEnvelopeWindowHullFailed
      (convexHullPolygon (exactLoopPoints (polygonOuterLoop window)))
  (dispositions, _, _, _) <-
    first UpperEnvelopePowerConstructionFailed
      (exactGeneratorDispositionsWith discardClipReceipt () domain generators)
  let envelopeLayer = affineDispositionLayer dispositions
  if convexPolygonComponent domain == window
    then Right envelopeLayer
    else restrictEnvelopeToWindow window envelopeLayer

powerCellDisposition
  :: Ord label
  => label
  -> BoundedPowerDiagram label
  -> Maybe (PowerCellDisposition label)
powerCellDisposition label (BoundedPowerDiagram dispositions) =
  Map.lookup label dispositions

powerCellDispositions
  :: BoundedPowerDiagram label
  -> [(label, PowerCellDisposition label)]
powerCellDispositions (BoundedPowerDiagram dispositions) = Map.toAscList dispositions

-- | Publish only full-dimensional cells.  The outside label is a caller-owned
-- view choice and may not collide with any submitted site label.
powerDiagramPlanarLayer
  :: Ord label
  => label
  -> BoundedPowerDiagram label
  -> Either (PowerDiagramError label) (PlanarLayer label)
powerDiagramPlanarLayer outside (BoundedPowerDiagram dispositions)
  | Map.member outside dispositions = Left (PowerDiagramOutsideLabelCollides outside)
  | otherwise =
      Right (publishedPowerLayer outside id dispositions)

affineDispositionLayer
  :: Ord label
  => Map label (PowerCellDisposition label)
  -> PlanarLayer (Maybe label)
affineDispositionLayer = publishedPowerLayer Nothing Just

publishedPowerLayer
  :: Ord publishedLabel
  => publishedLabel
  -> (label -> publishedLabel)
  -> Map label (PowerCellDisposition label)
  -> PlanarLayer publishedLabel
publishedPowerLayer outside publishLabel dispositions =
  planarLayerFromAdmittedComponents
    outside
    [ (publishLabel label, convexPolygonComponent polygon)
    | (label, PublishedPowerCell polygon) <- Map.toAscList dispositions
    ]

restrictEnvelopeToWindow
  :: Ord label
  => PolygonComponent
  -> PlanarLayer (Maybe label)
  -> Either (UpperEnvelopeError label) (PlanarLayer (Maybe label))
restrictEnvelopeToWindow window envelopeLayer = do
  let windowLayer = planarLayerFromAdmittedComponents False [(True, window)]
  clipped <-
    first UpperEnvelopeWindowOverlayFailed
      (overlayLayers envelopeLayer windowLayer)
  pure
    ( planarLayerFromAdmittedComponents
        Nothing
        [ (Just label, component)
        | cell <- Vector.toList (overlayResultCells clipped)
        , overlayCellRight cell
        , Just label <- [overlayCellLeft cell]
        , BoundedOverlayCell component <- [overlayCellGeometry cell]
        ]
    )

validateAndSortSites
  :: Ord label
  => NonEmpty (PowerSite label)
  -> Either (PowerDiagramError label) (NonEmpty (PowerSite label))
validateAndSortSites submitted =
  let sorted = NonEmpty.sortBy (comparing powerSiteLabel) submitted
   in case List.find (uncurry sameLabel) (consecutivePairs (NonEmpty.toList sorted)) of
        Just (duplicate, _) -> Left (DuplicatePowerSiteLabel (powerSiteLabel duplicate))
        Nothing -> Right sorted
 where
  sameLabel :: Eq label => PowerSite label -> PowerSite label -> Bool
  sameLabel left right = powerSiteLabel left == powerSiteLabel right

prepareExactPowerGenerator :: PowerSite label -> (ExactPowerGenerator label, Int)
prepareExactPowerGenerator site =
  let point = exactPointFromQueryPoint (powerSiteQueryPoint site)
      (coordinateX, coordinateY) = exactPointCoordinates point
      weight = powerWeightExact (powerSiteWeight site)
   in ( ExactPowerGenerator
          { exactPowerGeneratorLabel = powerSiteLabel site
          , exactPowerGeneratorXCoefficient = 2 * coordinateX
          , exactPowerGeneratorYCoefficient = 2 * coordinateY
          , exactPowerGeneratorConstant =
              weight - coordinateX * coordinateX - coordinateY * coordinateY
          }
      , max (exactPointBitWidth point) (exactRationalBitWidth weight)
      )

affineFormGenerators
  :: Map label AffineForm
  -> Either (UpperEnvelopeError label) (NonEmpty (ExactPowerGenerator label))
affineFormGenerators forms =
  case Map.minViewWithKey forms of
    Nothing -> Left UpperEnvelopeEmptyForms
    Just ((firstLabel, firstForm), remaining) ->
      Right
        ( affineFormGenerator firstLabel firstForm
            :| fmap (uncurry affineFormGenerator) (Map.toAscList remaining)
        )

affineFormGenerator
  :: label
  -> AffineForm
  -> ExactPowerGenerator label
affineFormGenerator label form =
  ExactPowerGenerator
    { exactPowerGeneratorLabel = label
    , exactPowerGeneratorXCoefficient = affineFormXCoefficient form
    , exactPowerGeneratorYCoefficient = affineFormYCoefficient form
    , exactPowerGeneratorConstant = affineFormConstant form
    }

exactGeneratorDispositionsWith
  :: Ord label
  => (summary -> ExactClipReceipt -> summary)
  -> summary
  -> ConvexPolygon
  -> NonEmpty (ExactPowerGenerator label)
  -> Either
      (PowerDiagramError label)
      ( Map label (PowerCellDisposition label)
      , summary
      , GeneratorRegularReceipt
      , Int
      )
exactGeneratorDispositionsWith summarizeReceipt initialSummary domain generators = do
  retainedDomain <-
    first PowerDomainInvalid (exactRetainedPolygon (convexPolygonPoints domain))
  section <- resolvedGeneratorSection generators
  let topology = resolvedRegularTopology section
  constraints <- prepareRegularConstraintSection topology
  clipped <-
    traverse
      (clipPowerCell retainedDomain section constraints)
      (generatorRegularDispositions topology)
  let (publishedAssociations, summary, maximumCellConstraints) =
        Foldable.foldl'
          (\(associations, accumulatedSummary, peak) (label, disposition, cellReceipt, cellAxes) ->
             let !combinedSummary = summarizeReceipt accumulatedSummary cellReceipt
              in ( (label, disposition) : associations
                 , combinedSummary
                 , max peak cellAxes
                 ))
          ([], initialSummary, 0)
          clipped
      publishedDispositions = Map.fromList publishedAssociations
      dispositions =
        fmap coincidentPowerDisposition (resolvedCoincidentDispositions section)
          <> publishedDispositions
  pure
    ( dispositions
    , summary
    , generatorRegularReceipt topology
    , maximumCellConstraints
    )

discardClipReceipt :: () -> ExactClipReceipt -> ()
discardClipReceipt _ _ = ()

groupGeneratorsBySlope
  :: NonEmpty (ExactPowerGenerator label)
  -> NonEmpty (NonEmpty (ExactPowerGenerator label))
groupGeneratorsBySlope (initial :| remaining) =
  let initialSlope = exactGeneratorSlope initial
      (sameInitialSlope, otherGenerators) =
        List.partition ((== initialSlope) . exactGeneratorSlope) remaining
      otherGroups =
        Map.fromListWith
          (<>)
          [ (exactGeneratorSlope generator, generator :| [])
          | generator <- otherGenerators
          ]
   in (initial :| sameInitialSlope) :| Map.elems otherGroups

exactGeneratorSlope
  :: ExactPowerGenerator label
  -> (ExactRational, ExactRational)
exactGeneratorSlope generator =
  ( exactPowerGeneratorXCoefficient generator
  , exactPowerGeneratorYCoefficient generator
  )

resolveCoincidentGeneratorGroup
  :: Ord label
  => NonEmpty (ExactPowerGenerator label)
  -> (ExactPowerGenerator label, [(label, CoincidentGeneratorDisposition label)])
resolveCoincidentGeneratorGroup generators@(initial :| remaining) =
  let representative = List.foldl' chooseRepresentative initial remaining
      representativeLabel = exactPowerGeneratorLabel representative
   in ( representative
      , [ ( exactPowerGeneratorLabel generator
          , classifyCoincidentGenerator representative generator
          )
        | generator <- NonEmpty.toList generators
        , exactPowerGeneratorLabel generator /= representativeLabel
        ]
      )
 where
  chooseRepresentative
    :: Ord label
    => ExactPowerGenerator label
    -> ExactPowerGenerator label
    -> ExactPowerGenerator label
  chooseRepresentative selected candidate =
    case compare
      (exactPowerGeneratorConstant candidate)
      (exactPowerGeneratorConstant selected) of
      GT -> candidate
      LT -> selected
      EQ ->
        if exactPowerGeneratorLabel candidate < exactPowerGeneratorLabel selected
          then candidate
          else selected

classifyCoincidentGenerator
  :: ExactPowerGenerator label
  -> ExactPowerGenerator label
  -> CoincidentGeneratorDisposition label
classifyCoincidentGenerator representative candidate
  | exactPowerGeneratorConstant candidate == exactPowerGeneratorConstant representative =
      CoincidentGeneratorEquivalentTo (exactPowerGeneratorLabel representative)
  | otherwise = CoincidentGeneratorDominatedBy (exactPowerGeneratorLabel representative)

data CoincidentGeneratorDisposition label
  = CoincidentGeneratorEquivalentTo !label
  | CoincidentGeneratorDominatedBy !label

coincidentPowerDisposition
  :: CoincidentGeneratorDisposition label
  -> PowerCellDisposition label
coincidentPowerDisposition disposition =
  case disposition of
    CoincidentGeneratorEquivalentTo label -> CoincidentEquivalentTo label
    CoincidentGeneratorDominatedBy label -> CoincidentDominatedBy label

coincidentRegularDisposition
  :: CoincidentGeneratorDisposition label
  -> RegularSiteDisposition label
coincidentRegularDisposition disposition =
  case disposition of
    CoincidentGeneratorEquivalentTo label -> RegularSiteCoincidentEquivalentTo label
    CoincidentGeneratorDominatedBy label -> RegularSiteCoincidentDominatedBy label

data ResolvedGeneratorSection label = ResolvedGeneratorSection
  { resolvedCoincidentDispositions :: !(Map label (CoincidentGeneratorDisposition label))
  , resolvedRegularTopology :: !(GeneratorRegularTopology label)
  }

resolvedGeneratorSection
  :: Ord label
  => NonEmpty (ExactPowerGenerator label)
  -> Either (PowerDiagramError label) (ResolvedGeneratorSection label)
resolvedGeneratorSection generators = do
  let resolvedGroups = fmap resolveCoincidentGeneratorGroup (groupGeneratorsBySlope generators)
      representatives = fmap fst resolvedGroups
      coincidentDispositions =
        Map.fromList (concatMap snd (NonEmpty.toList resolvedGroups))
  topology <-
    first PowerRegularTopologyFailed
      (regularGeneratorTopology (DistinctSlopeGenerators representatives))
  pure
    ResolvedGeneratorSection
      { resolvedCoincidentDispositions = coincidentDispositions
      , resolvedRegularTopology = topology
      }

publishRegularSection
  :: Ord label
  => ResolvedGeneratorSection label
  -> (RegularTriangulation label, RegularTriangulationReceipt)
publishRegularSection section =
  let topology = resolvedRegularTopology section
      coincident = resolvedCoincidentDispositions section
      dispositions =
        Map.fromList
          [ ( exactPowerGeneratorLabel generator
            , publishGeneratorDisposition disposition
            )
          | (generator, disposition) <-
              NonEmpty.toList (generatorRegularDispositions topology)
          ]
          <> fmap coincidentRegularDisposition coincident
      faces = fmap publishGeneratorFace (generatorRegularFaces topology)
      edges = fmap publishGeneratorEdge (generatorRegularEdges topology)
      generatorReceipt = generatorRegularReceipt topology
      coincidentCount = Map.size coincident
   in ( RegularTriangulation
          { storedRegularDispositions = dispositions
          , storedRegularFaces = faces
          , storedRegularEdges = edges
          , storedRegularNeighbours = regularNeighbourSection topology
          }
      , RegularTriangulationReceipt
          { regularTriangulationInputSites =
              generatorRegularInputSites generatorReceipt + coincidentCount
          , regularTriangulationRepresentativeSites =
              generatorRegularInputSites generatorReceipt
          , regularTriangulationVisibleSites =
              generatorRegularVisibleSites generatorReceipt
          , regularTriangulationLowerDimensionalSites =
              generatorRegularLowerDimensionalSites generatorReceipt
          , regularTriangulationHiddenSites =
              generatorRegularHiddenSites generatorReceipt
          , regularTriangulationCoincidentSites = coincidentCount
          , regularTriangulationFaces = generatorRegularFaceCount generatorReceipt
          , regularTriangulationEdges = generatorRegularEdgeCount generatorReceipt
          , regularTriangulationPeakHullFaces =
              generatorRegularPeakHullFaces generatorReceipt
          }
      )

publishGeneratorDisposition
  :: RegularGeneratorDisposition
  -> RegularSiteDisposition label
publishGeneratorDisposition disposition =
  case disposition of
    RegularGeneratorVisible -> RegularSiteVisible
    RegularGeneratorLowerDimensional -> RegularSiteLowerDimensional
    RegularGeneratorHidden -> RegularSiteHidden

publishGeneratorFace :: GeneratorRegularFace label -> RegularFace label
publishGeneratorFace face =
  RegularFace
    (generatorRegularFaceFirst face)
    (generatorRegularFaceSecond face)
    (generatorRegularFaceThird face)
    (generatorRegularFaceDualPoint face)

publishGeneratorEdge :: GeneratorRegularEdge label -> RegularEdge label
publishGeneratorEdge edge =
  RegularEdge
    (exactPowerGeneratorLabel (generatorRegularEdgeFirst edge))
    (exactPowerGeneratorLabel (generatorRegularEdgeSecond edge))
    (case generatorRegularEdgeDual edge of
       GeneratorDualSegment segment -> BoundedPowerDual segment
       GeneratorDualRay ray -> UnboundedPowerDual ray
       GeneratorDualLine line -> FullLinePowerDual line
       GeneratorDualCollapsed point -> CollapsedPowerDual point)

regularNeighbourSection
  :: Ord label
  => GeneratorRegularTopology label
  -> Map label (Set label)
regularNeighbourSection =
  List.foldl' insertRegularNeighbourEdge Map.empty . generatorRegularEdges

insertRegularNeighbourEdge
  :: Ord label
  => Map label (Set label)
  -> GeneratorRegularEdge label
  -> Map label (Set label)
insertRegularNeighbourEdge neighbours edge =
  let firstLabel = exactPowerGeneratorLabel (generatorRegularEdgeFirst edge)
      secondLabel = exactPowerGeneratorLabel (generatorRegularEdgeSecond edge)
   in Map.insertWith Set.union firstLabel (Set.singleton secondLabel)
        (Map.insertWith Set.union secondLabel (Set.singleton firstLabel) neighbours)

prepareRadicalAxis
  :: ExactPowerGenerator label
  -> ExactPowerGenerator label
  -> Either (PowerDiagramError label) ExactAffineLine
prepareRadicalAxis firstGenerator secondGenerator =
  first
    ( PowerBisectorInvalid
        (exactPowerGeneratorLabel firstGenerator)
        (exactPowerGeneratorLabel secondGenerator)
    )
    (exactGeneratorAxis firstGenerator secondGenerator)

prepareRegularConstraintSection
  :: Ord label
  => GeneratorRegularTopology label
  -> Either (PowerDiagramError label) (Map label [ExactClosedHalfPlane])
prepareRegularConstraintSection =
  Foldable.foldlM prepareRegularConstraintEdge Map.empty . generatorRegularEdges

prepareRegularConstraintEdge
  :: Ord label
  => Map label [ExactClosedHalfPlane]
  -> GeneratorRegularEdge label
  -> Either (PowerDiagramError label) (Map label [ExactClosedHalfPlane])
prepareRegularConstraintEdge constraints edge = do
  let firstGenerator = generatorRegularEdgeFirst edge
      secondGenerator = generatorRegularEdgeSecond edge
      firstLabel = exactPowerGeneratorLabel firstGenerator
      secondLabel = exactPowerGeneratorLabel secondGenerator
  axis <- prepareRadicalAxis firstGenerator secondGenerator
  pure
    ( Map.insertWith (<>) secondLabel [exactClosedHalfPlane (oppositeExactAffineLine axis)]
        (Map.insertWith (<>) firstLabel [exactClosedHalfPlane axis] constraints)
    )

clipPowerCell
  :: Ord label
  => ExactRetainedPolygon
  -> ResolvedGeneratorSection label
  -> Map label [ExactClosedHalfPlane]
  -> (ExactPowerGenerator label, RegularGeneratorDisposition)
  -> Either
      (PowerDiagramError label)
      (label, PowerCellDisposition label, ExactClipReceipt, Int)
clipPowerCell retainedDomain section constraints (ownerGenerator, disposition) =
  let ownerLabel = exactPowerGeneratorLabel ownerGenerator
      topology = resolvedRegularTopology section
   in case disposition of
    RegularGeneratorHidden ->
      pure (ownerLabel, EmptyPowerCell, mempty, 0)
    RegularGeneratorVisible ->
      finishPowerCell
        retainedDomain
        ownerGenerator
        (Map.findWithDefault [] ownerLabel constraints)
    RegularGeneratorLowerDimensional -> do
      let competitors =
            [ generator
            | (generator, _) <-
                NonEmpty.toList (generatorRegularDispositions topology)
            , exactPowerGeneratorLabel generator /= ownerLabel
            ]
      halfPlanes <- traverse (preparedDirectHalfPlane ownerGenerator) competitors
      finishPowerCell retainedDomain ownerGenerator halfPlanes

preparedDirectHalfPlane
  :: ExactPowerGenerator label
  -> ExactPowerGenerator label
  -> Either (PowerDiagramError label) ExactClosedHalfPlane
preparedDirectHalfPlane owner competitor =
  exactClosedHalfPlane <$> prepareRadicalAxis owner competitor

finishPowerCell
  :: ExactRetainedPolygon
  -> ExactPowerGenerator label
  -> [ExactClosedHalfPlane]
  -> Either
      (PowerDiagramError label)
      (label, PowerCellDisposition label, ExactClipReceipt, Int)
finishPowerCell retainedDomain ownerGenerator halfPlanes = do
  (exactDisposition, receipt) <-
    first (PowerCellClipFailed (exactPowerGeneratorLabel ownerGenerator))
      (exactClipRetainedPolygon retainedDomain halfPlanes)
  pure
    ( exactPowerGeneratorLabel ownerGenerator
    , case exactDisposition of
        ExactClipFullDimensional retained ->
          PublishedPowerCell (convexPolygonFromRetained retained)
        ExactClipLowerDimensional points -> LowerDimensionalPowerCell points
        ExactClipEmpty -> EmptyPowerCell
    , receipt
    , length halfPlanes
    )

aggregateReceipt
  :: ConvexPolygon
  -> Map label (PowerCellDisposition label)
  -> Int
  -> ExactClipReceipt
  -> GeneratorRegularReceipt
  -> Int
  -> PowerDiagramReceipt
aggregateReceipt domain dispositions siteInputBits clipReceipt regularReceipt maximumCellConstraints =
  let inputBits = max siteInputBits (exactClipInputCoordinateBits clipReceipt)
      peakBits = exactClipPeakIntermediateCoordinateBits clipReceipt
      finalBits = exactClipFinalCoordinateBits clipReceipt
      dispositionCounts = countPowerDispositions dispositions
   in PowerDiagramReceipt
        { powerDiagramDomainVertices = NonEmpty.length (convexPolygonPoints domain)
        , powerDiagramSubmittedSiteConstraints = exactClipSubmittedHalfPlanes clipReceipt
        , powerDiagramActiveBoundaries = exactClipActiveBoundaries clipReceipt
        , powerDiagramBoundaryCompatibilityChecks = exactClipBoundaryCompatibilityChecks clipReceipt
        , powerDiagramExactIntersections = exactClipExactIntersections clipReceipt
        , powerDiagramPublishedCells = countedPublishedCells dispositionCounts
        , powerDiagramLowerDimensionalCells = countedLowerDimensionalCells dispositionCounts
        , powerDiagramEmptyCells = countedEmptyCells dispositionCounts
        , powerDiagramCoincidentEquivalentCells = countedCoincidentEquivalentCells dispositionCounts
        , powerDiagramCoincidentDominatedCells = countedCoincidentDominatedCells dispositionCounts
        , powerDiagramRegularFaces = generatorRegularFaceCount regularReceipt
        , powerDiagramRegularEdges = generatorRegularEdgeCount regularReceipt
        , powerDiagramOracleCells = generatorRegularLowerDimensionalSites regularReceipt
        , powerDiagramMaximumCellConstraints = maximumCellConstraints
        , powerDiagramMaximumInputBits = inputBits
        , powerDiagramMaximumAffineCoefficientBits = exactClipMaximumAffineCoefficientBits clipReceipt
        , powerDiagramPeakIntermediateCoordinateBits = peakBits
        , powerDiagramFinalCoordinateBits = finalBits
        , powerDiagramFinalDenominatorBits = exactClipFinalDenominatorBits clipReceipt
        }

data PowerDispositionCounts = PowerDispositionCounts
  { countedPublishedCells :: !Int
  , countedLowerDimensionalCells :: !Int
  , countedEmptyCells :: !Int
  , countedCoincidentEquivalentCells :: !Int
  , countedCoincidentDominatedCells :: !Int
  }

countPowerDispositions
  :: Map label (PowerCellDisposition label)
  -> PowerDispositionCounts
countPowerDispositions =
  Map.foldl'
    (\counts disposition -> case disposition of
        PublishedPowerCell _ ->
          counts {countedPublishedCells = countedPublishedCells counts + 1}
        LowerDimensionalPowerCell _ ->
          counts {countedLowerDimensionalCells = countedLowerDimensionalCells counts + 1}
        EmptyPowerCell ->
          counts {countedEmptyCells = countedEmptyCells counts + 1}
        CoincidentEquivalentTo _ ->
          counts {countedCoincidentEquivalentCells = countedCoincidentEquivalentCells counts + 1}
        CoincidentDominatedBy _ ->
          counts {countedCoincidentDominatedCells = countedCoincidentDominatedCells counts + 1})
    PowerDispositionCounts
      { countedPublishedCells = 0
      , countedLowerDimensionalCells = 0
      , countedEmptyCells = 0
      , countedCoincidentEquivalentCells = 0
      , countedCoincidentDominatedCells = 0
      }