moonlight-planar-1.1.0.0: test/algebra/Moonlight/Planar/ValuationSpec.hs
-- | Exact intrinsic-volume fixtures and common-subdivision inclusion-exclusion.
module Moonlight.Planar.ValuationSpec (tests) where
import Data.Foldable (traverse_)
import Data.Ratio ((%), numerator, denominator)
import Moonlight.Planar.Exact (ExactPoint, exactPoint, exactRationalNumerator, exactRationalDenominator)
import Moonlight.Planar.Internal.ExactRational
( ExactRational
, exactRational
)
import Moonlight.Planar.Overlay
( overlayClosedIntersection
, overlayClosedUnion
, overlayLayers
, overlaySelectedRegion
)
import Moonlight.Planar.Region
( PlanarRegion
, emptyPlanarRegion
, planarRegion
, planarRegionComponents
)
import Moonlight.Planar.Valuation
( CertifiedInterval (..)
, ExactPlanarMoments
, ExactLengthTerm
, PlanarValuations
, ValuationError (ValuationCellSetNotPureRegion)
, cellSetPerimeter
, cellValuations
, exactAreaValue
, exactLengthBounds
, exactLengthExpression
, exactLengthTerms
, exactPlanarFirstX
, exactPlanarFirstY
, exactPlanarMeasure
, exactPlanarSecondXX
, exactPlanarSecondXY
, exactPlanarSecondYY
, eulerCharacteristicValue
, lengthCoefficient
, polygonComponentArea
, polygonComponentMoments
, orientedBoundaryMoments
, regionPerimeter
, regionValuations
, scaleExactPlanarMoments
, lengthRadicand
, valuationArea
, valuationEuler
, valuationIntrinsic1
)
import Moonlight.Planar.OverlayFixtures (singletonLayer)
import Support
( annulusRegion, assertEqual, polygonComponentOf, polygonRegion
, rectangleComponent, rectangleRegion, requireRight
)
tests :: IO ()
tests = do
testRegionGoldenValues
testExactPlanarMoments
testOrientedMomentRationalOracle
testClosedCellInclusionExclusion
testDimensionalCellFixtures
testMetricInvariance
testRadicalNormalForm
putStrLn "valuation: ok"
testExactPlanarMoments :: IO ()
testExactPlanarMoments = do
unit <- rectangleComponent 0 0 1 1
oneHalfValue <- exactValue 1 2
oneThird <- exactValue 1 3
oneQuarter <- exactValue 1 4
let unitMoments = polygonComponentMoments unit
assertEqual
"area is the measure projection"
(exactPlanarMeasure unitMoments)
(exactAreaValue (polygonComponentArea unit))
assertEqual
"unit-square moments"
(1, oneHalfValue, oneHalfValue, oneThird, oneQuarter, oneThird)
(momentDigest unitMoments)
translated <- rectangleComponent 2 (-1) 4 2
assertEqual
"translated-rectangle moments"
(6, 18, 3, 56, 9, 6)
(momentDigest (polygonComponentMoments translated))
annulus <- annulusRegion (0, 0, 3, 3) (1, 1, 2, 2)
seventyFourThirds <- exactValue 74 3
let annulusMoments = foldMap polygonComponentMoments (planarRegionComponents annulus)
assertEqual
"annulus moments respect clockwise holes"
(8, 12, 12, seventyFourThirds, 18, seventyFourThirds)
(momentDigest annulusMoments)
assertEqual
"moment monoid and scalar action"
(momentDigest (unitMoments <> unitMoments))
(momentDigest (scaleExactPlanarMoments 2 unitMoments))
momentDigest
:: ExactPlanarMoments
-> (ExactRational, ExactRational, ExactRational, ExactRational, ExactRational, ExactRational)
momentDigest moments =
( exactPlanarMeasure moments
, exactPlanarFirstX moments
, exactPlanarFirstY moments
, exactPlanarSecondXX moments
, exactPlanarSecondXY moments
, exactPlanarSecondYY moments
)
exactValue :: Integer -> Integer -> IO ExactRational
exactValue numeratorValue denominatorValue =
requireRight "exact valuation fixture" (exactRational numeratorValue denominatorValue)
testRegionGoldenValues :: IO ()
testRegionGoldenValues = do
assertRegionValuations "empty" emptyPlanarRegion 0 0 []
unit <- rectangleRegion 0 0 1 1
assertRegionValuations "unit square" unit 1 1 [(2, 1)]
unitPerimeter <- requireRight "unit-square perimeter" (regionPerimeter unit)
assertLength "unit-square perimeter" [(4, 1)] (exactLengthTerms (exactLengthExpression unitPerimeter))
assertContains "unit-square perimeter bounds" 4 (exactLengthBounds unitPerimeter)
-- Outer 3×3 square: four edges of length √9 = 3, half-perimeter 6. Hole:
-- half-perimeter 2. The unnormalized shape 2·√1 + 2·√9 folds to 8.
annulus <- annulusRegion (0, 0, 3, 3) (1, 1, 2, 2)
assertRegionValuations "annulus" annulus 0 8 [(8, 1)]
lowerLeft <- rectangleComponent 0 0 1 1
upperRight <- rectangleComponent 1 1 2 2
cornerTouch <- requireRight "corner-touch region" (planarRegion [lowerLeft, upperRight])
assertRegionValuations "corner-touch squares" cornerTouch 1 2 [(4, 1)]
right <- rectangleComponent 1 0 2 1
edgeTouch <- requireRight "edge-touch region" (planarRegion [lowerLeft, right])
assertRegionValuations "edge-sharing squares" edgeTouch 1 2 [(3, 1)]
-- The double cell's bottom edge covers the unit cell's top edge from a
-- shared corner: collinear, positive overlap, shared endpoint. The contact
-- is one segment (V − E = 2 − 1), so Euler 1 + 1 − 1. Area 1 + 2. Intrinsic
-- one-volume: half-perimeters 2 + 3 minus the shared unit length, which is
-- 2·√1 + (1·√1 + 1·√4) − 1·√1 before normalization and 4 after.
doubleAbove <- rectangleComponent 0 1 2 2
overlapTouch <-
requireRight "shared-endpoint overlap region" (planarRegion [lowerLeft, doubleAbove])
assertRegionValuations "shared-endpoint collinear overlap" overlapTouch 1 3 [(4, 1)]
testClosedCellInclusionExclusion :: IO ()
testClosedCellInclusionExclusion = do
leftRegion <- rectangleRegion 0 0 1 1
rightRegion <- rectangleRegion 1 0 2 1
leftLayer <- singletonLayer False True leftRegion
rightLayer <- singletonLayer False True rightRegion
overlay <- requireRight "edge-sharing valuation overlay" (overlayLayers leftLayer rightLayer)
left <- requireRight "left closed cells" (overlayClosedUnion id (const False) overlay)
right <- requireRight "right closed cells" (overlayClosedUnion (const False) id overlay)
union <- requireRight "union closed cells" (overlayClosedUnion id id overlay)
intersection <-
requireRight
"intersection closed cells"
(overlayClosedIntersection id id overlay)
leftValues <- requireRight "left valuations" (cellValuations left)
rightValues <- requireRight "right valuations" (cellValuations right)
unionValues <- requireRight "union valuations" (cellValuations union)
sharedEdgeValues <- requireRight "intersection valuations" (cellValuations intersection)
assertEqual
"Euler inclusion-exclusion retains the shared edge"
(1, 1, 1, 1)
( eulerCharacteristicValue (valuationEuler leftValues)
, eulerCharacteristicValue (valuationEuler rightValues)
, eulerCharacteristicValue (valuationEuler unionValues)
, eulerCharacteristicValue (valuationEuler sharedEdgeValues)
)
assertEqual
"area inclusion-exclusion retains zero-dimensional measure"
( 1
, 1
, 2
, 0
)
( areaExact leftValues
, areaExact rightValues
, areaExact unionValues
, areaExact sharedEdgeValues
)
assertLength
"left intrinsic one-volume"
[(2, 1)]
(exactLengthTerms (exactLengthExpression (valuationIntrinsic1 leftValues)))
assertLength
"right intrinsic one-volume"
[(2, 1)]
(exactLengthTerms (exactLengthExpression (valuationIntrinsic1 rightValues)))
assertLength
"union intrinsic one-volume"
[(3, 1)]
(exactLengthTerms (exactLengthExpression (valuationIntrinsic1 unionValues)))
assertLength
"shared-edge intrinsic one-volume"
[(1, 1)]
(exactLengthTerms (exactLengthExpression (valuationIntrinsic1 sharedEdgeValues)))
case cellSetPerimeter intersection of
Left ValuationCellSetNotPureRegion -> pure ()
other ->
fail
( "an isolated selected edge was accepted as a perimeter: "
<> show other
)
published <-
requireRight
"published edge-sharing union"
(overlaySelectedRegion (uncurry (||)) overlay)
publishedValues <- requireRight "published union valuations" (regionValuations published)
assertEqual
"cell and published region valuations agree"
(valuationDigest unionValues)
(valuationDigest publishedValues)
testDimensionalCellFixtures :: IO ()
testDimensionalCellFixtures = do
disjoint <- intersectionValues (0, 0, 1, 1) (2, 0, 3, 1)
assertEqual
"empty intersection valuations"
(0, 0, [])
(exactDigest disjoint)
point <- intersectionValues (0, 0, 1, 1) (1, 1, 2, 2)
assertEqual
"point intersection valuations"
(1, 0, [])
(exactDigest point)
-- The intersection is the 1×2 rectangle (1,0)–(2,2): half-perimeter 3,
-- which is 1·√1 + 1·√4 before normalization.
area <- intersectionValues (0, 0, 2, 2) (1, 0, 3, 2)
assertEqual
"two-dimensional intersection valuations"
(1, 2, [(3, 1)])
(exactDigest area)
testMetricInvariance :: IO ()
testMetricInvariance = do
original <- polygonRegion [(0, 0), (1, 0), (0, 1)]
translated <- polygonRegion [(5, -3), (6, -3), (5, -2)]
quarterTurned <- polygonRegion [(0, 0), (0, 1), (-1, 0)]
originalValues <- requireRight "original triangle valuations" (regionValuations original)
translatedValues <- requireRight "translated triangle valuations" (regionValuations translated)
quarterTurnedValues <- requireRight "quarter-turned triangle valuations" (regionValuations quarterTurned)
assertEqual "translation invariance" (exactDigest originalValues) (exactDigest translatedValues)
assertEqual "quarter-turn invariance" (exactDigest originalValues) (exactDigest quarterTurnedValues)
perimeter <- requireRight "irrational triangle perimeter" (regionPerimeter original)
assertContains
"certified radical perimeter"
(2 + sqrt 2)
(exactLengthBounds perimeter)
testRadicalNormalForm :: IO ()
testRadicalNormalForm = do
-- Diamonds with edge vectors (1, 1) and (2, 2): half-perimeters 2·√2 and
-- 2·√8 = 4·√2. The square factor of the trial prime 2 folds, and the two
-- lengths are one term. Areas 2 + 8.
small <- polygonComponentOf [(0, 1), (1, 0), (2, 1), (1, 2)]
large <- polygonComponentOf [(4, 2), (6, 0), (8, 2), (6, 4)]
diamonds <- requireRight "diamond region" (planarRegion [small, large])
assertRegionValuations "diamonds share one radical" diamonds 2 10 [(6, 2)]
-- Rhombi with edge vectors (536, 201) and (8, 3): squared lengths
-- 327697 = 67²·73 and 73. The prime 67 is beyond the trial primes and
-- 327697 is not a square, so only the class test (327697·73 = 4891²) can
-- fold 2·√327697 into 134·√73; with the second rhombus's 2·√73 the sum is
-- 136·√73. Areas 327697 + 73.
wide <- polygonComponentOf [(0, 0), (536, 201), (335, 737), (-201, 536)]
narrow <- polygonComponentOf [(1000, 0), (1008, 3), (1005, 11), (997, 8)]
rhombi <- requireRight "rhombus region" (planarRegion [wide, narrow])
assertRegionValuations
"radical classes merge beyond the trial primes"
rhombi
2
327770
[(136, 73)]
intersectionValues
:: (Integer, Integer, Integer, Integer)
-> (Integer, Integer, Integer, Integer)
-> IO PlanarValuations
intersectionValues leftBounds rightBounds = do
leftRegion <- uncurryRectangle leftBounds
rightRegion <- uncurryRectangle rightBounds
leftLayer <- singletonLayer False True leftRegion
rightLayer <- singletonLayer False True rightRegion
overlay <- requireRight "dimensional valuation overlay" (overlayLayers leftLayer rightLayer)
selected <-
requireRight
"dimensional closed intersection"
(overlayClosedIntersection id id overlay)
requireRight "dimensional cell valuations" (cellValuations selected)
uncurryRectangle
:: (Integer, Integer, Integer, Integer)
-> IO PlanarRegion
uncurryRectangle (minimumX, minimumY, maximumX, maximumY) =
rectangleRegion minimumX minimumY maximumX maximumY
assertRegionValuations
:: String
-> PlanarRegion
-> Int
-> Integer
-> [(Integer, Integer)]
-> IO ()
assertRegionValuations label region expectedEuler expectedArea expectedLength = do
values <- requireRight label (regionValuations region)
assertEqual
(label <> " Euler")
expectedEuler
(eulerCharacteristicValue (valuationEuler values))
assertEqual
(label <> " area")
(fromInteger expectedArea)
(exactAreaValue (valuationArea values))
assertLength
(label <> " intrinsic one-volume")
expectedLength
(exactLengthTerms (exactLengthExpression (valuationIntrinsic1 values)))
assertLength :: String -> [(Integer, Integer)] -> [ExactLengthTerm] -> IO ()
assertLength label expected actual =
assertEqual
label
[ (fromInteger coefficient, radicand)
| (coefficient, radicand) <- expected
]
[ (lengthCoefficient term, lengthRadicand term)
| term <- actual
]
assertContains :: String -> Double -> CertifiedInterval -> IO ()
assertContains label expected interval =
if intervalLower interval <= expected && expected <= intervalUpper interval
then pure ()
else fail (label <> ": interval does not contain " <> show expected <> ": " <> show interval)
valuationDigest
:: PlanarValuations
-> (Int, ExactRational, CertifiedInterval)
valuationDigest values =
( eulerCharacteristicValue (valuationEuler values)
, exactAreaValue (valuationArea values)
, exactLengthBounds (valuationIntrinsic1 values)
)
exactDigest
:: PlanarValuations
-> (Int, ExactRational, [(ExactRational, Integer)])
exactDigest values =
( eulerCharacteristicValue (valuationEuler values)
, exactAreaValue (valuationArea values)
, [ (lengthCoefficient term, lengthRadicand term)
| term <- exactLengthTerms (exactLengthExpression (valuationIntrinsic1 values))
]
)
areaExact :: PlanarValuations -> ExactRational
areaExact values =
exactAreaValue (valuationArea values)
-- Independent normalized Rational formulas, deliberately not expressed through
-- the production raw-sum carrier or its homogeneous-integer implementation.
rationalEdgeMoments :: ((Rational, Rational), (Rational, Rational)) -> [Rational]
rationalEdgeMoments ((x0, y0), (x1, y1)) =
let cross = x0 * y1 - y0 * x1
in [ cross / 2
, (x0 + x1) * cross / 6
, (y0 + y1) * cross / 6
, (x0 * x0 + x0 * x1 + x1 * x1) * cross / 12
, (2 * x0 * y0 + x0 * y1 + x1 * y0 + 2 * x1 * y1) * cross / 24
, (y0 * y0 + y0 * y1 + y1 * y1) * cross / 12
]
rationalBoundaryMoments :: [((Rational, Rational), (Rational, Rational))] -> [Rational]
rationalBoundaryMoments = foldr (zipWith (+) . rationalEdgeMoments) [0, 0, 0, 0, 0, 0]
rationalMomentValues :: ExactPlanarMoments -> [Rational]
rationalMomentValues moments =
map (\value -> exactRationalNumerator value % exactRationalDenominator value)
[ exactPlanarMeasure moments, exactPlanarFirstX moments, exactPlanarFirstY moments
, exactPlanarSecondXX moments, exactPlanarSecondXY moments, exactPlanarSecondYY moments
]
momentEdgeFamilies :: [(String, [((Rational, Rational), (Rational, Rational))])]
momentEdgeFamilies =
[ ("integer", [((n, 2*n+1), (n+3, 1-n)) | n <- [-4 .. 3]])
, ("common denominator", [((n/64, (2*n+1)/64), ((n+3)/64, (1-n)/64)) | n <- [-4 .. 3]])
, ("mixed coprime denominators", [((n/101, (2*n+1)/103), ((n+3)/107, (1-n)/109)) | n <- [-4 .. 3]])
, ("high bit width", [(((huge+n)/p, (huge-n)/q), ((huge+2*n+1)/r, (n-huge)/s)) | n <- [-4 .. 3]])
, ("zero cross", [((n/101, 2*n/101), ((n+3)/107, 2*(n+3)/107)) | n <- [-4 .. 3]])
]
where
huge = fromInteger (2 ^ (191 :: Int))
p = fromInteger (2 ^ (127 :: Int) - 1)
q = fromInteger (2 ^ (107 :: Int) - 1)
r = fromInteger (2 ^ (89 :: Int) - 1)
s = fromInteger (2 ^ (61 :: Int) - 1)
admitRationalMomentPoint :: (Rational, Rational) -> IO ExactPoint
admitRationalMomentPoint (x, y) =
exactPoint <$> exactValue (numerator x) (denominator x)
<*> exactValue (numerator y) (denominator y)
admitRationalMomentEdge :: ((Rational, Rational), (Rational, Rational)) -> IO (ExactPoint, ExactPoint)
admitRationalMomentEdge (from, to) =
(,) <$> admitRationalMomentPoint from <*> admitRationalMomentPoint to
testOrientedMomentRationalOracle :: IO ()
testOrientedMomentRationalOracle = do
traverse_ checkFamily momentEdgeFamilies
putStrLn "oriented moments: 40 Rational edge oracles, reversal, subdivision and gluing passed"
where
checkFamily :: (String, [((Rational, Rational), (Rational, Rational))]) -> IO ()
checkFamily (label, edges) = do
admitted <- traverse admitRationalMomentEdge edges
let actual = orientedBoundaryMoments admitted
assertEqual (label <> " six independent Rational formulas")
(rationalBoundaryMoments edges) (rationalMomentValues actual)
traverse_
(\(edge, exactEdge) ->
assertEqual (label <> " individual edge")
(rationalEdgeMoments edge)
(rationalMomentValues (orientedBoundaryMoments [exactEdge])))
(zip edges admitted)
assertEqual (label <> " orientation reversal")
(map negate (rationalMomentValues actual))
(rationalMomentValues (orientedBoundaryMoments (map (\(from, to) -> (to, from)) (reverse admitted))))
subdivided <- traverse admitRationalMomentEdge
(concatMap (\(from@(x0,y0), to@(x1,y1)) ->
let middle = ((x0+x1)/2, (y0+y1)/2)
in [(from,middle),(middle,to)]) edges)
assertEqual (label <> " exact midpoint subdivision") actual (orientedBoundaryMoments subdivided)
let (prefix, suffix) = splitAt 3 admitted
assertEqual (label <> " finite gluing") actual
(orientedBoundaryMoments prefix <> orientedBoundaryMoments suffix)
assertEqual (label <> " empty boundary identity") actual
(orientedBoundaryMoments [] <> actual)