fix: numeric correctness (statistics, metrics, parsing)

dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs

@@ -217,7 +217,9 @@ varScatter takeSqrt g nGroups v = runST $ do
             | otherwise = do
                 c <- VUM.unsafeRead cnt k
                 mm <- VUM.unsafeRead m2 k
-                let var = if c < 2 then 0 else mm / fromIntegral (c - 1)
+                -- Sample variance is undefined at n = 1: NaN, matching
+                -- 'computeVariance'.
+                let var = if c < 2 then 0 / 0 else mm / fromIntegral (c - 1)
                 VUM.unsafeWrite out k (if takeSqrt then sqrt var else var)
                 fin (k + 1)
     fin 0

dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs

@@ -290,7 +290,9 @@ varPar takeSqrt vis offs nGroups v caps bounds = do
             | otherwise = do
                 c <- VUM.unsafeRead cnt k
                 mm <- VUM.unsafeRead m2 k
-                let var = if c < 2 then 0 else mm / fromIntegral (c - 1)
+                -- Sample variance is undefined at n = 1: NaN, matching
+                -- 'computeVariance'.
+                let var = if c < 2 then 0 / 0 else mm / fromIntegral (c - 1)
                 VUM.unsafeWrite out k (if takeSqrt then sqrt var else var)
                 fin (k + 1)
     fin 0

dataframe-learn/src/DataFrame/Metrics.hs

@@ -73,15 +73,22 @@ columnOf df e = case interpret @Double df e of
     Right (TColumn c) -> fromRight VU.empty (toVector @Double @VU.Vector c)
     Left err -> throw err
 
-n2 :: VU.Vector Double -> Double
-n2 = fromIntegral . VU.length
+{- | Compared pairs: 'VU.zipWith' truncates to the shorter vector, so every
+mean below divides by this, never by the length of 'truth' alone.
+-}
+nCompared :: VU.Vector Double -> VU.Vector Double -> Double
+nCompared preds truth = fromIntegral (min (VU.length preds) (VU.length truth))
 
 -- | Mean squared error.
 mse :: Metric
 mse preds truth
     | VU.null truth = 0
+    -- No predictions is not a perfect score.
+    | n == 0 = 0 / 0
     | otherwise =
-        VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n2 truth
+        VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n
+  where
+    n = nCompared preds truth
 
 -- | Root mean squared error.
 rmse :: Metric
@@ -91,29 +98,40 @@ rmse preds truth = sqrt (mse preds truth)
 mae :: Metric
 mae preds truth
     | VU.null truth = 0
-    | otherwise = VU.sum (VU.zipWith (\p t -> abs (p - t)) preds truth) / n2 truth
+    | n == 0 = 0 / 0
+    | otherwise = VU.sum (VU.zipWith (\p t -> abs (p - t)) preds truth) / n
+  where
+    n = nCompared preds truth
 
 -- | Coefficient of determination @R²@.
 r2 :: Metric
 r2 preds truth
-    | VU.null truth || ssTot == 0 = 0
+    | VU.null truth = 0
+    | n == 0 = 0 / 0
+    | ssTot == 0 = 0
     | otherwise = 1 - ssRes / ssTot
   where
-    mean = VU.sum truth / n2 truth
-    ssRes = VU.sum (VU.zipWith (\p t -> (t - p) ^ (2 :: Int)) preds truth)
-    ssTot = VU.sum (VU.map (\t -> (t - mean) ^ (2 :: Int)) truth)
+    n = nCompared preds truth
+    truth' = VU.take (min (VU.length preds) (VU.length truth)) truth
+    mean = VU.sum truth' / n
+    ssRes = VU.sum (VU.zipWith (\p t -> (t - p) ^ (2 :: Int)) preds truth')
+    ssTot = VU.sum (VU.map (\t -> (t - mean) ^ (2 :: Int)) truth')
 
 -- | Fraction of exact matches.
 accuracy :: Metric
 accuracy preds truth
     | VU.null truth = 0
+    | n == 0 = 0 / 0
     | otherwise =
-        fromIntegral (VU.length (VU.filter id (VU.zipWith (==) preds truth))) / n2 truth
+        fromIntegral (VU.length (VU.filter id (VU.zipWith (==) preds truth))) / n
+  where
+    n = nCompared preds truth
 
 -- | Binary log loss; probabilities clamped away from @0@/@1@.
 logLoss :: Metric
 logLoss probs truth
     | VU.null truth = 0
+    | n == 0 = 0 / 0
     | otherwise =
         negate
             ( VU.sum
@@ -123,8 +141,9 @@ logLoss probs truth
                     truth
                 )
             )
-            / n2 truth
+            / n
   where
+    n = nCompared probs truth
     clampP p = max 1e-15 (min (1 - 1e-15) p)
 
 -- | Averaging strategy for multiclass precision/recall/F1.

dataframe-learn/tests-internal/Learn/EdgeCases.hs

@@ -28,7 +28,7 @@ import DataFrame.LinearModel
 import DataFrame.LinearSolver (sigmoid)
 import DataFrame.PCA
 
-import DataFrame.Internal.Statistics (correlation', variance')
+import DataFrame.Internal.Statistics (correlation', meanSquaredError, variance')
 
 import Test.HUnit
 
@@ -124,14 +124,13 @@ testVarianceConstant = TestCase $ do
     let v = variance' (VU.replicate 100 (7.0 :: Double))
     assertEqual "variance of constant column is 0" 0 v
 
-{- Variance of fewer than two samples is defined to be 0 (computeVariance guard),
-   not NaN from a /0. -}
+{- Sample variance of one observation is undefined: NaN, so a singleton
+   group can never look as tight as a genuinely constant column. -}
 testVarianceSingleton :: Test
 testVarianceSingleton = TestCase $ do
-    assertEqual
-        "variance of one sample is 0"
-        0
-        (variance' (VU.fromList [3.5 :: Double]))
+    assertBool
+        "variance of one sample is NaN"
+        (isNaN (variance' (VU.fromList [3.5 :: Double])))
 
 {- Correlation of a perfectly linear pair is exactly +1 (and -1 reversed),
    computed stably. y = 2x+1 over a spread of x. -}
@@ -169,6 +168,27 @@ testCorrelationTooFew = TestCase $ do
         Nothing
         (correlation' (VU.fromList [1]) (VU.fromList [2]))
 
+{- meanSquaredError refuses length mismatches and empty inputs rather than
+   averaging over terms it never summed (or indexing out of bounds). -}
+testMeanSquaredErrorGuards :: Test
+testMeanSquaredErrorGuards = TestCase $ do
+    assertEqual
+        "mse of mismatched lengths is Nothing"
+        Nothing
+        (meanSquaredError (VU.fromList [0, 0, 0, 0]) (VU.fromList [2, 2]))
+    assertEqual
+        "mse with the longer prediction does not index out of bounds"
+        Nothing
+        (meanSquaredError (VU.fromList [1]) (VU.fromList [1, 2, 3]))
+    assertEqual
+        "mse of empty inputs is Nothing"
+        Nothing
+        (meanSquaredError VU.empty VU.empty)
+    assertEqual
+        "mse of equal-length inputs is the plain mean"
+        (Just 4.0)
+        (meanSquaredError (VU.fromList [0, 0]) (VU.fromList [2, 2]))
+
 -- ===========================================================================
 -- Category 8: stability inside the model expr layer
 -- ===========================================================================
@@ -425,6 +445,7 @@ tests =
     , testCorrelationPerfect
     , testCorrelationConstantColumnIsNaN
     , testCorrelationTooFew
+    , testMeanSquaredErrorGuards
     , testLogisticProbsExtremeFeatures
     , testOLSOneRow
     , testLogisticSingleClass

dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs

@@ -17,7 +17,9 @@ import DataFrame.Errors (DataFrameException (..))
 mean' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
 mean' samp
     | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = rtf (VU.sum samp) / fromIntegral (VU.length samp)
+    -- Widen per element: summing at 'a' wraps for Int columns.
+    | otherwise =
+        VU.foldl' (\acc x -> acc + rtf x) 0 samp / fromIntegral (VU.length samp)
 {-# INLINE [0] mean' #-}
 
 meanDouble' :: VU.Vector Double -> Double
@@ -29,7 +31,9 @@ meanDouble' samp
 meanInt' :: VU.Vector Int -> Double
 meanInt' samp
     | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = fromIntegral (VU.sum samp) / fromIntegral (VU.length samp)
+    | otherwise =
+        VU.foldl' (\acc x -> acc + fromIntegral x) 0 samp
+            / fromIntegral (VU.length samp)
 {-# INLINE meanInt' #-}
 
 {-# RULES
@@ -54,7 +58,8 @@ median' samp
             then pure (rtf middleElement)
             else do
                 prev <- VUM.read mutableSamp (middleIndex - 1)
-                pure (rtf (middleElement + prev) / 2)
+                -- Widen before adding: 'a' addition wraps for Int.
+                pure ((rtf middleElement + rtf prev) / 2)
 {-# INLINE median' #-}
 
 -- accumulator: count, mean, m2
@@ -73,7 +78,9 @@ varianceStep (VarAcc !n !meanVal !m2) !x =
 
 computeVariance :: VarAcc -> Double
 computeVariance (VarAcc !n _ !m2)
-    | n < 2 = 0 -- or error "variance of <2 samples"
+    | n == 0 = throw $ EmptyDataSetException "variance"
+    -- Sample variance is undefined at n = 1: NaN, not a spurious 0.
+    | n < 2 = 0 / 0
     | otherwise = m2 / fromIntegral (n - 1)
 {-# INLINE computeVariance #-}
 
@@ -106,38 +113,37 @@ skewnessStep (SkewAcc !n !meanVal !m2 !m3) !x' =
 computeSkewness :: SkewAcc -> Double
 computeSkewness (SkewAcc n _ m2 m3)
     | n < 3 = 0 -- or error "skewness of <3 samples"
-    | otherwise = (sqrt (fromIntegral n - 1) * m3) / sqrt (m2 ^ (3 :: Int))
+    -- m2, m3 are raw sums, so population g1 = sqrt n * m3 / m2^(3/2).
+    | otherwise = (sqrt (fromIntegral n) * m3) / sqrt (m2 ^ (3 :: Int))
 {-# INLINE computeSkewness #-}
 
 skewness' :: (VU.Unbox a, Real a, Num a) => VU.Vector a -> Double
 skewness' = computeSkewness . VU.foldl' skewnessStep (SkewAcc 0 0 0 0)
 {-# INLINE skewness' #-}
 
-data CorrelationStats
-    = CorrelationStats
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-
+{- | Centered two-pass form: the one-pass @n*Sxy - Sx*Sy@ form cancels
+catastrophically on low-variance columns and can report |r| > 1.
+-}
 correlation' :: VU.Vector Double -> VU.Vector Double -> Maybe Double
 correlation' xs ys
     | n < 2 = Nothing
     | VU.length xs /= VU.length ys = Nothing
     | otherwise =
         let nf = fromIntegral n
-            initial = CorrelationStats 0 0 0 0 0
-            (CorrelationStats sumX sumY sumXX sumYY sumXY) = VU.ifoldl' step initial xs
-
-            !num = nf * sumXY - sumX * sumY
-            !den = sqrt ((nf * sumXX - sumX * sumX) * (nf * sumYY - sumY * sumY))
-         in Just (num / den)
+            !mx = VU.sum xs / nf
+            !my = VU.sum ys / nf
+            !sxy = VU.sum (VU.zipWith (\x y -> (x - mx) * (y - my)) xs ys)
+            !sxx = VU.sum (VU.map (\x -> (x - mx) * (x - mx)) xs)
+            !syy = VU.sum (VU.map (\y -> (y - my) * (y - my)) ys)
+         in Just (clamp (sxy / sqrt (sxx * syy)))
   where
     n = VU.length xs
-    step (CorrelationStats sx sy sxx syy sxy) i x =
-        let !y = VU.unsafeIndex ys i
-         in CorrelationStats (sx + x) (sy + y) (sxx + x * x) (syy + y * y) (sxy + x * y)
+    -- Cauchy-Schwarz: clamp the last rounding, never the formula. Explicit
+    -- guards so the zero-variance NaN passes through ('max' would eat it).
+    clamp r
+        | r > 1 = 1
+        | r < -1 = -1
+        | otherwise = r
 {-# INLINE correlation' #-}
 
 quantiles' ::
@@ -202,11 +208,14 @@ interQuartileRange' samp =
 {-# INLINE interQuartileRange' #-}
 
 meanSquaredError :: VU.Vector Double -> VU.Vector Double -> Maybe Double
-meanSquaredError target prediction =
-    let
-        squareDiff = VU.ifoldl' (\sq i e -> (e - target VU.! i) ^ (2 :: Int) + sq) 0 prediction
-     in
-        Just $ squareDiff / fromIntegral (max (VU.length target) (VU.length prediction))
+meanSquaredError target prediction
+    | VU.length target /= VU.length prediction = Nothing
+    | VU.null target = Nothing
+    | otherwise =
+        Just
+            ( VU.sum (VU.zipWith (\t p -> (p - t) ^ (2 :: Int)) target prediction)
+                / fromIntegral (VU.length target)
+            )
 {-# INLINE meanSquaredError #-}
 
 mutualInformationBinned ::

dataframe-operations/src/DataFrame/Operations/Inference.hs

@@ -65,15 +65,9 @@ byteStringDateParser "%Y-%m-%d" = parseDateField
 byteStringDateParser fmt = readByteStringDate fmt
 {-# INLINE byteStringDateParser #-}
 
-{- | 'DataFrame.Internal.Parsing.readInt' that rejects overflow instead of
-wrapping. Fields of <= 18 chars cannot overflow and keep the Text-level
-parse; longer (rare) fields take the exact byte-level parser, so an
-overflowing cell demotes\/promotes instead of silently wrapping.
--}
+-- | Alias for 'readInt', which now rejects overflow itself.
 readIntStrict :: T.Text -> Maybe Int
-readIntStrict t
-    | T.length t <= 18 = readInt t
-    | otherwise = parseIntField (TE.encodeUtf8 t)
+readIntStrict = readInt
 {-# INLINE readIntStrict #-}
 
 {- | Candidate-mask priority, reproducing the documented fallback order:

dataframe-operations/src/DataFrame/Operations/Statistics.hs

@@ -416,7 +416,10 @@ summarize df =
 
 -- | Round a @Double@ to Specified Precision
 roundTo :: Int -> Double -> Double
-roundTo n x = fromInteger (round $ x * 10 ^ n) / 10.0 ^^ n
+roundTo n x
+    -- 'round' on NaN yields garbage; keep NaN visible in summaries.
+    | isNaN x = x
+    | otherwise = fromInteger (round $ x * 10 ^ n) / 10.0 ^^ n
 
 toPct2dp :: Double -> String
 toPct2dp x

dataframe-operations/src/DataFrame/Operations/Typing.hs

@@ -217,12 +217,12 @@ handleBoolAssumption isNull cols =
         (handleTextAssumption isNull cols)
 
 {- | Int columns: one fused pass with in-place Int -> Double promotion; a cell
-parsing as neither demotes the column to Text. 'readIntStrict' rejects overflow
+parsing as neither demotes the column to Text. 'readInt' rejects overflow
 so a huge integer promotes to 'Double' rather than wrapping.
 -}
 handleIntAssumption :: (T.Text -> Bool) -> V.Vector T.Text -> Column
 handleIntAssumption isNull cols =
-    case promoteIntColumn (\_ t -> isNull t) readIntStrict readDouble cols of
+    case promoteIntColumn (\_ t -> isNull t) readInt readDouble cols of
         Just col -> col
         Nothing -> handleTextAssumption isNull cols
 

dataframe-parsing/dataframe-parsing.cabal

@@ -43,7 +43,6 @@ library
     build-depends:      base >= 4 && < 5,
                         attoparsec >= 0.12 && < 0.16,
                         bytestring >= 0.11 && < 0.14,
-                        bytestring-lexing >= 0.5 && < 0.7,
                         containers >= 0.6.7 && < 0.10,
                         dataframe-core >= 2.4 && < 2.5,
                         text >= 2.1 && < 3,