streamly-statistics 0.1.0 → 0.2.0

raw patch · 6 files changed

+1233/−251 lines, 6 filesdep ~QuickCheckdep ~containersdep ~deepseqnew-uploaderPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: QuickCheck, containers, deepseq, random, streamly-core, tasty, tasty-bench

API changes (from Hackage documentation)

Files

• CHANGELOG.md +5/−0
• benchmark/Main.hs +302/−177
• src/Streamly/Statistics.hs +80/−40
• src/Streamly/Statistics/Scanl.hs +787/−0
• streamly-statistics.cabal +22/−15
• test/Main.hs +37/−19

@@ -1,5 +1,10 @@ # Changelog +## 0.2.0 (Sep 2025)++* Add a Scanl module+* Support `streamly-core-0.3.0`+ ## 0.1.0 (Apr 2023)  * Initial version

@@ -2,14 +2,17 @@  import Control.DeepSeq (NFData) import Streamly.Data.Fold (Fold)+import Streamly.Data.Scanl (Scanl) import Streamly.Data.Stream (Stream) import System.Random (randomRIO)  import qualified Streamly.Data.Fold as Fold+import qualified Streamly.Internal.Data.RingArray as Ring+import qualified Streamly.Internal.Data.Scanl as Scanl import qualified Streamly.Data.Stream as Stream import qualified Streamly.Data.Array as Array-import qualified Streamly.Internal.Data.Ring.Unboxed as Ring import qualified Streamly.Statistics as Statistics+import qualified Streamly.Statistics.Scanl as StatScan  import Gauge @@ -47,228 +50,350 @@ benchWithFoldInt :: Int -> String -> Fold IO Int Int -> Benchmark benchWithFoldInt len name f = benchWith source len name f +{-# INLINE benchWithScanSrc #-}+benchWithScanSrc :: (Num a) =>+    (Int -> a -> Stream IO a) -> Int -> String -> Scanl IO a a -> Benchmark+benchWithScanSrc src len name f =+    bench name+        $ nfIO+        $ randomRIO (1, 1 :: Int)+        >>= Stream.fold Fold.drain+            . Stream.postscanl f . src len . fromIntegral+ {-# INLINE benchWithPostscan #-}-benchWithPostscan :: Int -> String -> Fold IO Double Double -> Benchmark+benchWithPostscan :: Int -> String -> Scanl IO Double Double -> Benchmark benchWithPostscan len name f =   bench name $ nfIO $ randomRIO (1, 1) >>=-    Stream.fold Fold.drain . Stream.postscan f . source len+    Stream.fold Fold.drain . Stream.postscanl f . source len  {-# INLINE benchWithResample #-} benchWithResample :: Int -> String -> Benchmark benchWithResample len name = bench name $ nfIO $ do     i <- randomRIO (1, 1)-    arr <- Stream.fold Array.write (source len i :: Stream IO Double)+    arr <- Stream.fold Array.create (source len i :: Stream IO Double)     Stream.fold Fold.drain $ Stream.unfold Statistics.resample arr  {-# INLINE benchWithFoldResamples #-} benchWithFoldResamples :: Int -> String -> Fold IO Double Double -> Benchmark benchWithFoldResamples len name f = bench name $ nfIO $ do     i <- randomRIO (1, 1)-    arr <- Stream.fold Array.write (source len i :: Stream IO Double)+    arr <- Stream.fold Array.create (source len i :: Stream IO Double)     Stream.fold Fold.drain $ Statistics.foldResamples len arr f  {-# INLINE numElements #-} numElements :: Int numElements = 100000 -main :: IO ()-main =-  defaultMain-    [ bgroup-        "fold"-        [ benchWithFold numElements "minimum (window size 100)"-            (Ring.slidingWindow 100 Statistics.minimum)-        , benchWithFold numElements "minimum (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.minimum)-        , benchWith sourceDescendingInt numElements-            "minimum descending (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.minimum)+mkBenchmarks ::+       (Int -> String -> Fold IO Double Double -> Benchmark)+    -> [Benchmark]+mkBenchmarks mkBench =+    [+      mkBench numElements "minimum (window size 100)"+        (Ring.slidingWindow 100 Statistics.minimum)+    , mkBench numElements "minimum (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.minimum)+    , benchWith sourceDescendingInt numElements+        "minimum descending (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.minimum) -        , benchWithFold numElements "maximum (window size 100)"-            (Ring.slidingWindow 100 Statistics.maximum)-        , benchWithFold numElements "maximum (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.maximum)-        , benchWith sourceDescendingInt numElements-            "maximum descending (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.maximum)+    , mkBench numElements "maximum (window size 100)"+        (Ring.slidingWindow 100 Statistics.maximum)+    , mkBench numElements "maximum (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.maximum)+    , benchWith sourceDescendingInt numElements+        "maximum descending (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.maximum) -        , benchWithFold numElements "range (window size 100)"-            (Ring.slidingWindow 100 Statistics.range)-        , benchWithFold numElements "range (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.range)+    , mkBench numElements "range (window size 100)"+        (Ring.slidingWindow 100 Statistics.range)+    , mkBench numElements "range (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.range) -        , benchWithFoldInt numElements "sumInt (window size 100)"-            (Ring.slidingWindow 100 Statistics.sumInt)-        , benchWithFoldInt numElements "sum for Int (window size 100)"-            (Ring.slidingWindow 100 Statistics.sum)+    , mkBench numElements "sum (window size 100)"+        (Ring.slidingWindow 100 Statistics.sum)+    , mkBench numElements "sum (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.sum)+    , mkBench numElements "sum (entire stream)"+        (Statistics.cumulative Statistics.sum)+    , mkBench numElements "sum (Data.Fold)" Fold.sum -        , benchWithFold numElements "sum (window size 100)"-            (Ring.slidingWindow 100 Statistics.sum)-        , benchWithFold numElements "sum (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.sum)-        , benchWithFold numElements "sum (entire stream)"-            (Statistics.cumulative Statistics.sum)-        , benchWithFold numElements "sum (Data.Fold)"-            (Fold.sum)+    , mkBench numElements "mean (window size 100)"+        (Ring.slidingWindow 100 Statistics.mean)+    , mkBench numElements "mean (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.mean)+    , mkBench numElements "mean (entire stream)"+        (Statistics.cumulative Statistics.mean)+    , mkBench numElements "mean (Data.Fold)" Fold.mean -        , benchWithFold numElements "mean (window size 100)"-            (Ring.slidingWindow 100 Statistics.mean)-        , benchWithFold numElements "mean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.mean)-        , benchWithFold numElements "mean (entire stream)"-            (Statistics.cumulative Statistics.mean)-        , benchWithFold numElements "mean (Data.Fold)"-            (Fold.mean)+    , mkBench numElements "welfordMean (window size 100)"+        (Ring.slidingWindow 100 Statistics.welfordMean)+    , mkBench numElements "welfordMean (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.welfordMean)+    , mkBench numElements "welfordMean (entire stream)"+        (Statistics.cumulative Statistics.welfordMean) -        , benchWithFold-            numElements-            "welfordMean (window size 100)"-            (Ring.slidingWindow 100 Statistics.welfordMean)-        , benchWithFold-            numElements-            "welfordMean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.welfordMean)-        , benchWithFold-            numElements-            "welfordMean (entire stream)"-            (Statistics.cumulative Statistics.welfordMean)+    , mkBench numElements "geometricMean (window size 100)"+        (Ring.slidingWindow 100 Statistics.geometricMean)+    , mkBench numElements "geometricMean (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.geometricMean)+    , mkBench numElements "geometricMean (entire stream)"+        (Statistics.cumulative Statistics.geometricMean) -        , benchWithFold numElements "geometricMean (window size 100)"-            (Ring.slidingWindow 100 Statistics.geometricMean)-        , benchWithFold numElements "geometricMean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.geometricMean)-        , benchWithFold numElements "geometricMean (entire stream)"-            (Statistics.cumulative Statistics.geometricMean)+    , mkBench numElements "harmonicMean (window size 100)"+        (Ring.slidingWindow 100 Statistics.harmonicMean)+    , mkBench numElements "harmonicMean (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.harmonicMean)+    , mkBench numElements "harmonicMean (entire stream)"+        (Statistics.cumulative Statistics.harmonicMean) -        , benchWithFold numElements "harmonicMean (window size 100)"-            (Ring.slidingWindow 100 Statistics.harmonicMean)-        , benchWithFold numElements "harmonicMean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.harmonicMean)-        , benchWithFold numElements "harmonicMean (entire stream)"-            (Statistics.cumulative Statistics.harmonicMean)+    , mkBench numElements "quadraticMean (window size 100)"+        (Ring.slidingWindow 100 Statistics.quadraticMean)+    , mkBench numElements "quadraticMean (window size 1000)"+        (Ring.slidingWindow 1000 Statistics.quadraticMean)+    , mkBench numElements "quadraticMean (entire stream)"+        (Statistics.cumulative Statistics.quadraticMean) -        , benchWithFold numElements "quadraticMean (window size 100)"-            (Ring.slidingWindow 100 Statistics.quadraticMean)-        , benchWithFold numElements "quadraticMean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.quadraticMean)-        , benchWithFold numElements "quadraticMean (entire stream)"-            (Statistics.cumulative Statistics.quadraticMean)+    , mkBench numElements "powerSum 2 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerSum 2))+    , mkBench numElements "powerSum 2 (entire stream)"+        (Statistics.cumulative (Statistics.powerSum 2)) -        , benchWithFold numElements "powerSum 2 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerSum 2))-        , benchWithFold numElements "powerSum 2 (entire stream)"-            (Statistics.cumulative (Statistics.powerSum 2))+    , mkBench numElements "rawMoment 2 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerSum 2))+    , mkBench numElements "rawMoment 2 (entire stream)"+        (Statistics.cumulative (Statistics.rawMoment 2)) -        , benchWithFold numElements "rawMoment 2 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerSum 2))-        , benchWithFold numElements "rawMoment 2 (entire stream)"-            (Statistics.cumulative (Statistics.rawMoment 2))+    , mkBench numElements "powerMean 1 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMean 1))+    , mkBench numElements "powerMean 2 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMean 2))+    , mkBench numElements "powerMean 10 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMean 10)) -        , benchWithFold numElements "powerMean 1 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMean 1))-        , benchWithFold numElements "powerMean 2 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMean 2))-        , benchWithFold numElements "powerMean 10 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMean 10))+    , mkBench numElements "powerMeanFrac (-1) (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMeanFrac (-1)))+    , mkBench numElements "powerMeanFrac 1 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMeanFrac 1))+    , mkBench numElements "powerMeanFrac 2 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMeanFrac 2))+    , mkBench numElements "powerMeanFrac 10 (window size 100)"+        (Ring.slidingWindow 100 (Statistics.powerMeanFrac 10)) -        , benchWithFold numElements "powerMeanFrac (-1) (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMeanFrac (-1)))-        , benchWithFold numElements "powerMeanFrac 1 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMeanFrac 1))-        , benchWithFold numElements "powerMeanFrac 2 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMeanFrac 2))-        , benchWithFold numElements "powerMeanFrac 10 (window size 100)"-            (Ring.slidingWindow 100 (Statistics.powerMeanFrac 10))+    , mkBench numElements "ewma (entire stream)"+        (Statistics.ewma 0.5)+    , mkBench numElements "ewmaAfterMean (entire stream)"+         (Statistics.ewmaAfterMean 10 0.5)+    , mkBench numElements "ewmaRampUpSmoothing (entire stream)"+        (Statistics.ewmaRampUpSmoothing 0.5 0.5) -        , benchWithFold numElements "ewma (entire stream)"-            (Statistics.ewma 0.5)-        , benchWithFold numElements "ewmaAfterMean (entire stream)"-            (Statistics.ewmaAfterMean 10 0.5)-        , benchWithFold numElements "ewmaRampUpSmoothing (entire stream)"-            (Statistics.ewmaRampUpSmoothing 0.5 0.5)+    , mkBench numElements "variance (window size 100)"+        (Ring.slidingWindow 100 Statistics.variance)+    , mkBench numElements "variance (entire stream)"+        (Statistics.cumulative Statistics.variance)+    -- , mkBench numElements "variance (Data.Fold)" Fold.variance -        , benchWithFold numElements "variance (window size 100)"-            (Ring.slidingWindow 100 (Statistics.variance))-        , benchWithFold numElements "variance (entire stream)"-            (Statistics.cumulative (Statistics.variance))-        -- , benchWithFold numElements "variance (Data.Fold)"-        --     (Fold.variance)+    , mkBench numElements "sampleVariance (window size 100)"+        (Ring.slidingWindow 100 Statistics.sampleVariance)+    , mkBench numElements "sampleVariance (entire stream)"+        (Statistics.cumulative Statistics.sampleVariance) -        , benchWithFold numElements "sampleVariance (window size 100)"-            (Ring.slidingWindow 100 (Statistics.sampleVariance))-        , benchWithFold numElements "sampleVariance (entire stream)"-            (Statistics.cumulative (Statistics.sampleVariance))+    , mkBench numElements "stdDev (window size 100)"+        (Ring.slidingWindow 100 Statistics.stdDev)+    , mkBench numElements "stdDev (entire stream)"+        (Statistics.cumulative Statistics.stdDev)+    -- , mkBench numElements "stdDev (Data.Fold)" Fold.stdDev -        , benchWithFold numElements "stdDev (window size 100)"-            (Ring.slidingWindow 100 (Statistics.stdDev))-        , benchWithFold numElements "stdDev (entire stream)"-            (Statistics.cumulative (Statistics.stdDev))-        -- , benchWithFold numElements "stdDev (Data.Fold)"-        --     (Fold.stdDev)+    , mkBench numElements "sampleStdDev (window size 100)"+        (Ring.slidingWindow 100 Statistics.sampleStdDev)+    , mkBench numElements "sampleStdDev (entire stream)"+        (Statistics.cumulative Statistics.sampleStdDev) -        , benchWithFold numElements "sampleStdDev (window size 100)"-            (Ring.slidingWindow 100 (Statistics.sampleStdDev))-        , benchWithFold numElements "sampleStdDev (entire stream)"-            (Statistics.cumulative (Statistics.sampleStdDev))+    , mkBench numElements "stdErrMean (window size 100)"+        (Ring.slidingWindow 100 Statistics.stdErrMean)+    , mkBench numElements "stdErrMean (entire stream)"+        (Statistics.cumulative Statistics.stdErrMean) -        , benchWithFold numElements "stdErrMean (window size 100)"-            (Ring.slidingWindow 100 (Statistics.stdErrMean))-        , benchWithFold numElements "stdErrMean (entire stream)"-            (Statistics.cumulative (Statistics.stdErrMean))+-- These benchmarks take a lot of time/memory with fusion-plugin possibly+-- because of the use of Tee.+#ifndef FUSION_PLUGIN+    , mkBench numElements "skewness (window size 100)"+        (Ring.slidingWindow 100 Statistics.skewness)+    , mkBench numElements "skewness (entire stream)"+        (Statistics.cumulative Statistics.skewness) +    , mkBench numElements "kurtosis (window size 100)"+        (Ring.slidingWindow 100 Statistics.kurtosis)+    , mkBench numElements "kurtosis (entire stream)"+        (Statistics.cumulative Statistics.kurtosis)+#endif+    , mkBench numElements "md (window size 100)"+        (Ring.slidingWindowWith 100 Statistics.md)++    ]++mkScans ::+       (Int -> String -> Scanl IO Double Double -> Benchmark)+    -> [Benchmark]+mkScans mkBench =+    [+      mkBench numElements "minimum (window size 100)"+        (Scanl.incrScan 100 StatScan.incrMinimum)+    , mkBench numElements "minimum (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrMinimum)+    , benchWithScanSrc sourceDescendingInt numElements+        "minimum descending (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrMinimum)++    , mkBench numElements "maximum (window size 100)"+        (Scanl.incrScan 100 StatScan.incrMaximum)+    , mkBench numElements "maximum (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrMaximum)+    , benchWithScanSrc sourceDescendingInt numElements+        "maximum descending (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrMaximum)++    , mkBench numElements "range (window size 100)"+        (Scanl.incrScan 100 StatScan.incrRange)+    , mkBench numElements "range (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrRange)++    , mkBench numElements "sum (window size 100)"+        (Scanl.incrScan 100 Scanl.incrSum)+    , mkBench numElements "sum (window size 1000)"+        (Scanl.incrScan 1000 Scanl.incrSum)+    , mkBench numElements "sum (entire stream)"+        (Scanl.cumulativeScan Scanl.incrSum)+    , mkBench numElements "sum (Data.Fold)" Scanl.sum++    , mkBench numElements "mean (window size 100)"+        (Scanl.incrScan 100 Scanl.incrMean)+    , mkBench numElements "mean (window size 1000)"+        (Scanl.incrScan 1000 Scanl.incrMean)+    , mkBench numElements "mean (entire stream)"+        (Scanl.cumulativeScan Scanl.incrMean)+    , mkBench numElements "mean (Data.Fold)" Scanl.mean++    , mkBench numElements "welfordMean (window size 100)"+        (Scanl.incrScan 100 StatScan.incrWelfordMean)+    , mkBench numElements "welfordMean (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrWelfordMean)+    , mkBench numElements "welfordMean (entire stream)"+        (Scanl.cumulativeScan StatScan.incrWelfordMean)++    , mkBench numElements "geometricMean (window size 100)"+        (Scanl.incrScan 100 StatScan.incrGeometricMean)+    , mkBench numElements "geometricMean (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrGeometricMean)+    , mkBench numElements "geometricMean (entire stream)"+        (Scanl.cumulativeScan StatScan.incrGeometricMean)++    , mkBench numElements "harmonicMean (window size 100)"+        (Scanl.incrScan 100 StatScan.incrHarmonicMean)+    , mkBench numElements "harmonicMean (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrHarmonicMean)+    , mkBench numElements "harmonicMean (entire stream)"+        (Scanl.cumulativeScan StatScan.incrHarmonicMean)++    , mkBench numElements "quadraticMean (window size 100)"+        (Scanl.incrScan 100 StatScan.incrQuadraticMean)+    , mkBench numElements "quadraticMean (window size 1000)"+        (Scanl.incrScan 1000 StatScan.incrQuadraticMean)+    , mkBench numElements "quadraticMean (entire stream)"+        (Scanl.cumulativeScan StatScan.incrQuadraticMean)++    , mkBench numElements "powerSum 2 (window size 100)"+        (Scanl.incrScan 100 (Scanl.incrPowerSum 2))+    , mkBench numElements "powerSum 2 (entire stream)"+        (Scanl.cumulativeScan (Scanl.incrPowerSum 2))++    , mkBench numElements "rawMoment 2 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrRawMoment 2))+    , mkBench numElements "rawMoment 2 (entire stream)"+        (Scanl.cumulativeScan (StatScan.incrRawMoment 2))++    , mkBench numElements "powerMean 1 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMean 1))+    , mkBench numElements "powerMean 2 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMean 2))+    , mkBench numElements "powerMean 10 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMean 10))++    , mkBench numElements "powerMeanFrac (-1) (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMeanFrac (-1)))+    , mkBench numElements "powerMeanFrac 1 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMeanFrac 1))+    , mkBench numElements "powerMeanFrac 2 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMeanFrac 2))+    , mkBench numElements "powerMeanFrac 10 (window size 100)"+        (Scanl.incrScan 100 (StatScan.incrPowerMeanFrac 10))++    , mkBench numElements "ewma (entire stream)"+        (StatScan.ewma 0.5)+    , mkBench numElements "ewmaRampUpSmoothing (entire stream)"+        (StatScan.ewmaRampUpSmoothing 0.5 0.5)++    , mkBench numElements "variance (window size 100)"+        (Scanl.incrScan 100 StatScan.incrVariance)+    , mkBench numElements "variance (entire stream)"+        (Scanl.cumulativeScan StatScan.incrVariance)+    -- , mkBench numElements "variance (Data.Fold)" Fold.variance++    , mkBench numElements "sampleVariance (window size 100)"+        (Scanl.incrScan 100 StatScan.incrSampleVariance)+    , mkBench numElements "sampleVariance (entire stream)"+        (Scanl.cumulativeScan StatScan.incrSampleVariance)++    , mkBench numElements "stdDev (window size 100)"+        (Scanl.incrScan 100 StatScan.incrStdDev)+    , mkBench numElements "stdDev (entire stream)"+        (Scanl.cumulativeScan StatScan.incrStdDev)+    -- , mkBench numElements "stdDev (Data.Fold)" Fold.stdDev++    , mkBench numElements "sampleStdDev (window size 100)"+        (Scanl.incrScan 100 StatScan.incrSampleStdDev)+    , mkBench numElements "sampleStdDev (entire stream)"+        (Scanl.cumulativeScan StatScan.incrSampleStdDev)++    , mkBench numElements "stdErrMean (window size 100)"+        (Scanl.incrScan 100 StatScan.incrStdErrMean)+    , mkBench numElements "stdErrMean (entire stream)"+        (Scanl.cumulativeScan StatScan.incrStdErrMean)+ -- These benchmarks take a lot of time/memory with fusion-plugin possibly -- because of the use of Tee. #ifndef FUSION_PLUGIN-        , benchWithFold numElements "skewness (window size 100)"-            (Ring.slidingWindow 100 (Statistics.skewness))-        , benchWithFold numElements "skewness (entire stream)"-            (Statistics.cumulative (Statistics.skewness))+    , mkBench numElements "skewness (window size 100)"+        (Scanl.incrScan 100 StatScan.windowSkewness)+    , mkBench numElements "skewness (entire stream)"+        (Scanl.cumulativeScan StatScan.windowSkewness) -        , benchWithFold numElements "kurtosis (window size 100)"-            (Ring.slidingWindow 100 (Statistics.kurtosis))-        , benchWithFold numElements "kurtosis (entire stream)"-            (Statistics.cumulative (Statistics.kurtosis))+    , mkBench numElements "kurtosis (window size 100)"+        (Scanl.incrScan 100 StatScan.windowKurtosis)+    , mkBench numElements "kurtosis (entire stream)"+        (Scanl.cumulativeScan StatScan.windowKurtosis) #endif-        , benchWithFold numElements "md (window size 100)"-            (Ring.slidingWindowWith 100 Statistics.md)-        ]-    , bgroup-        "scan"-        [ benchWithPostscan numElements "minimum (window size 100)"-            (Ring.slidingWindow 100 Statistics.minimum)-        , benchWithPostscan numElements "minimum (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.minimum)-        , benchWithPostscan numElements "maximum (window size 100)"-            (Ring.slidingWindow 100 Statistics.maximum)-        , benchWithPostscan numElements "maximum (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.maximum)-        , benchWithPostscan numElements "range (window size 100)"-            (Ring.slidingWindow 100 Statistics.range)-        , benchWithPostscan numElements "range (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.range)-        , benchWithPostscan numElements "sum (window size 100)"+    , mkBench numElements "md (window size 100)"+        (Scanl.incrScanWith 100 StatScan.incrMd)++    ]++main :: IO ()+main =+  defaultMain+    [+      bgroup "fold" $ mkBenchmarks benchWithFold+    , bgroup "fold_Int"+        [ benchWithFoldInt numElements "sumInt (window size 100)"+            (Ring.slidingWindow 100 Statistics.sumInt)+        , benchWithFoldInt numElements "sum for Int (window size 100)"             (Ring.slidingWindow 100 Statistics.sum)-        , benchWithPostscan numElements "sum (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.sum)-        , benchWithPostscan numElements "mean (window size 100)"-            (Ring.slidingWindow 100 Statistics.mean)-        , benchWithPostscan numElements "mean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.mean)-        , benchWithPostscan-            numElements-            "welfordMean (window size 100)"-            (Ring.slidingWindow 100 Statistics.welfordMean)-        , benchWithPostscan-            numElements-            "welfordMean (window size 1000)"-            (Ring.slidingWindow 1000 Statistics.welfordMean)-        , benchWithPostscan-            numElements-            "md (window size 100)"-            (Ring.slidingWindowWith 100 Statistics.md)-        -- XXX These benchmarks measure the cost of creating the array as well,-        -- we can do that outside the benchmark.-        , benchWithResample numElements "Resample"+        ]+    , bgroup "scan" $ mkScans benchWithPostscan+    -- XXX These benchmarks measure the cost of creating the array as well,+    -- we can do that outside the benchmark.+    , bgroup "resample"+        [ benchWithResample numElements "Resample"         , benchWithFoldResamples 316 "FoldResamples 316" Fold.mean         ]     ]

@@ -22,6 +22,10 @@  -- Resources: --+-- Related:+-- https://hackage.haskell.org/package/foldl-statistics+-- https://hackage.haskell.org/package/foldl-incremental+-- -- This may be another useful resource for incremental (non-windowed) -- computation: --@@ -80,17 +84,17 @@     -- window folds by keeping the second element of the input tuple as     -- @Nothing@.     ---      Window.lmap+      lmap     , Window.cumulative      -- * Summary Statistics     -- | See https://en.wikipedia.org/wiki/Summary_statistics .      -- ** Sums-    , Window.length-    , Window.sum-    , Window.sumInt-    , Window.powerSum+    , length+    , sum+    , sumInt+    , powerSum      -- ** Location     -- | See https://en.wikipedia.org/wiki/Location_parameter .@@ -183,25 +187,24 @@ import Data.Functor.Identity (runIdentity, Identity) import Data.Map.Strict (Map) import Data.Maybe (fromMaybe)-import Streamly.Data.Array (Array, length, Unbox)+import Streamly.Data.Array (Array, Unbox) import Streamly.Data.Fold (Tee(..)) import Streamly.Data.Stream (Stream)-import Streamly.Internal.Data.Array.Type (unsafeIndexIO)-import Streamly.Internal.Data.Fold.Type (Fold(..), Step(..))-import Streamly.Internal.Data.Stream.StreamD.Step (Step(..))+import Streamly.Internal.Data.Array (unsafeIndexIO)+import Streamly.Internal.Data.Fold (Fold(..), Step(..))+import Streamly.Internal.Data.Stream (Step(..)) import Streamly.Internal.Data.Tuple.Strict (Tuple'(..), Tuple3'(..))-import Streamly.Internal.Data.Unfold.Type (Unfold(..))+import Streamly.Internal.Data.Unfold (Unfold(..)) import System.Random.MWC (createSystemRandom, uniformRM)  import qualified Data.Map.Strict as Map import qualified Deque.Strict as Deque import qualified Streamly.Data.Fold as Fold-import qualified Streamly.Data.Array as Array hiding (read)-import qualified Streamly.Internal.Data.Array as Array (read)+import qualified Streamly.Data.Array as Array import qualified Streamly.Data.MutArray as MA-import qualified Streamly.Internal.Data.Array.Mut as MA-    (getIndexUnsafe, putIndexUnsafe, unsafeSwapIndices)-import qualified Streamly.Internal.Data.Fold.Window as Window+import qualified Streamly.Internal.Data.MutArray as MA+    (unsafeSwapIndices)+import qualified Streamly.Internal.Data.Fold as Window import qualified Streamly.Data.Stream as Stream  import Prelude hiding (length, sum, minimum, maximum)@@ -216,6 +219,32 @@ -- think about deduplication.  -------------------------------------------------------------------------------+-- Re-exports+-------------------------------------------------------------------------------++-- XXX Deprecate these once the streamly functions are released.++-- {-# DEPRECATED lmap "Use Streamly.Data.Fold.windowLmap instead" #-}+lmap :: (c -> a) -> Fold m (a, Maybe a) b -> Fold m (c, Maybe c) b+lmap = Window.windowLmap++-- {-# DEPRECATED length "Use Streamly.Data.Fold.windowLength instead" #-}+length :: (Monad m, Num b) => Fold m (a, Maybe a) b+length = Window.windowLength++-- {-# DEPRECATED sum "Use Streamly.Data.Fold.windowSum instead" #-}+sum :: (Monad m, Num a) => Fold m (a, Maybe a) a+sum = Window.windowSum++-- {-# DEPRECATED sumInt "Use Streamly.Data.Fold.windowSumInt instead" #-}+sumInt :: (Monad m, Integral a) => Fold m (a, Maybe a) a+sumInt = Window.windowSumInt++-- {-# DEPRECATED powerSum "Use Streamly.Data.Fold.windowPowerSum instead" #-}+powerSum :: (Monad m, Num a) => Int -> Fold m (a, Maybe a) a+powerSum = Window.windowPowerSum++------------------------------------------------------------------------------- -- Transforms ------------------------------------------------------------------------------- @@ -305,13 +334,13 @@                     let butterfly i | i >= len  = flight (j + 1) (a + e)                                     | otherwise = do                             let i1 = i + l1-                            xi1 :+ yi1 <- MA.getIndexUnsafe i1 marr+                            xi1 :+ yi1 <- MA.unsafeGetIndex i1 marr                             let !c = cos a                                 !s = sin a                                 d  = (c * xi1 - s * yi1) :+ (s * xi1 + c * yi1)-                            ci <- MA.getIndexUnsafe i marr-                            MA.putIndexUnsafe i1 marr (ci - d)-                            MA.putIndexUnsafe i marr (ci + d)+                            ci <- MA.unsafeGetIndex i marr+                            MA.unsafePutIndex i1 marr (ci - d)+                            MA.unsafePutIndex i marr (ci + d)                             butterfly (i + l2)                     butterfly j             flight 0 0@@ -320,6 +349,10 @@ -- Location ------------------------------------------------------------------------------- +-- XXX prefix window to these folds and implement these using the scans. or+-- just remove these as the corresponding scans can be converted to folds. If+-- we remove these we can just export the scans through this module itself.+ -- Theoretically, we can approximate minimum in a rolling window by using a -- 'powerMean' with sufficiently large negative power. --@@ -329,7 +362,7 @@ -- -- | The minimum element in a rolling window. ----- For smaller window sizes (< 30) Streamly.Data.Fold.Window.minimum performs+-- For smaller window sizes (< 30) Streamly.Internal.Data.Fold.windowMinimum performs -- better.  If you want to compute the minimum of the entire stream Fold.min -- from streamly package would be much faster. --@@ -337,7 +370,7 @@ -- {-# INLINE minimum #-} minimum :: (Monad m, Ord a) => Fold m (a, Maybe a) a-minimum = Fold step initial extract+minimum = Fold step initial extract extract      where @@ -389,7 +422,7 @@ -- -- | The maximum element in a rolling window. ----- For smaller window sizes (< 30) Streamly.Data.Fold.Window.maximum performs+-- For smaller window sizes (< 30) Streamly.Internal.Data.Fold.windowMaximum performs -- better.  If you want to compute the maximum of the entire stream -- Streamly.Data.Fold.maximum from streamly package would be much faster. --@@ -397,7 +430,7 @@ -- {-# INLINE maximum #-} maximum :: (Monad m, Ord a) => Fold m (a, Maybe a) a-maximum = Fold step initial extract+maximum = Fold step initial extract extract      where @@ -469,7 +502,7 @@ -- /Time/: \(\mathcal{O}(n)\) {-# INLINE mean #-} mean :: forall m a. (Monad m, Fractional a) => Fold m (a, Maybe a) a-mean = Window.mean+mean = Window.windowMean  -- | Recompute mean from old mean when an item is removed from the sample. {-# INLINE _meanSubtract #-}@@ -511,7 +544,7 @@ -- /Internal/ {-# INLINE welfordMean #-} welfordMean :: forall m a. (Monad m, Fractional a) => Fold m (a, Maybe a) a-welfordMean = Fold step initial extract+welfordMean = Fold step initial extract extract      where @@ -543,7 +576,7 @@ -- -- \(\mu'_k = \frac{\sum_{i=1}^n x_{i}^k}{n}\) ----- >>> rawMoment k = Fold.teeWith (/) (powerSum p) length+-- >>> rawMoment k = Fold.teeWith (/) (Fold.windowPowerSum p) Fold.windowLength -- -- See https://en.wikipedia.org/wiki/Moment_(mathematics) . --@@ -552,16 +585,17 @@ -- /Time/: \(\mathcal{O}(n)\) {-# INLINE rawMoment #-} rawMoment :: (Monad m, Fractional a) => Int -> Fold m (a, Maybe a) a-rawMoment k = Fold.teeWith (/) (Window.powerSum k) Window.length+rawMoment k = Fold.teeWith (/) (Window.windowPowerSum k) Window.windowLength  -- | Like 'rawMoment' but powers can be negative or fractional. This is -- slower than 'rawMoment' for positive intergal powers. ----- >>> rawMomentFrac p = Fold.teeWith (/) (powerSumFrac p) length+-- >>> rawMomentFrac p = Fold.teeWith (/) (Fold.windowPowerSumFrac p) Fold.windowLength -- {-# INLINE rawMomentFrac #-} rawMomentFrac :: (Monad m, Floating a) => a -> Fold m (a, Maybe a) a-rawMomentFrac k = Fold.teeWith (/) (Window.powerSumFrac k) Window.length+rawMomentFrac k =+    Fold.teeWith (/) (Window.windowPowerSumFrac k) Window.windowLength  -- XXX Overflow can happen when large powers or large numbers are used. We can -- keep a running mean instead of running sum but that won't mitigate the@@ -608,7 +642,9 @@ -- {-# INLINE harmonicMean #-} harmonicMean :: (Monad m, Fractional a) => Fold m (a, Maybe a) a-harmonicMean = Fold.teeWith (/) Window.length (Window.lmap recip Window.sum)+harmonicMean =+    Fold.teeWith (/)+        Window.windowLength (Window.windowLmap recip Window.windowSum)  -- | Geometric mean, defined as: --@@ -625,7 +661,7 @@ -- See https://en.wikipedia.org/wiki/Geometric_mean . {-# INLINE geometricMean #-} geometricMean :: (Monad m, Floating a) => Fold m (a, Maybe a) a-geometricMean = exp <$> Window.lmap log mean+geometricMean = exp <$> Window.windowLmap log mean  -- | The quadratic mean or root mean square (rms) of the numbers -- \(x_1, x_2, \ldots, x_n\) is defined as:@@ -685,8 +721,9 @@      extract (Tuple' x _) = x --- XXX It can perhaps perform better if implemented as a custom fold?---+-- XXX It can perhaps perform better if implemented as a custom fold? We can+-- also enable this to be used as a scan that way.+ -- | @ewma n k@ is like 'ewma' but uses the mean of the first @n@ values and -- then uses that as the initial value for the @ewma@ of the rest of the -- values.@@ -694,6 +731,7 @@ -- This can be used to reduce the effect of volatility of the initial value -- when k is too small. --+-- Note that this cannot be used as a scan. {-# INLINE ewmaAfterMean #-} ewmaAfterMean :: Monad m => Int -> Double -> Fold m Double Double ewmaAfterMean n k =@@ -897,7 +935,8 @@ -- {-# INLINE sampleVariance #-} sampleVariance :: (Monad m, Fractional a) => Fold m (a, Maybe a) a-sampleVariance = Fold.teeWith (\n s2 -> n * s2 / (n - 1)) Window.length variance+sampleVariance =+    Fold.teeWith (\n s2 -> n * s2 / (n - 1)) Window.windowLength variance  -- | Sample standard deviation: --@@ -923,7 +962,8 @@ -- /Time/: \(\mathcal{O}(n)\) {-# INLINE stdErrMean #-} stdErrMean :: (Monad m, Floating a) => Fold m (a, Maybe a) a-stdErrMean = Fold.teeWith (\sd n -> sd / sqrt n) sampleStdDev Window.length+stdErrMean =+    Fold.teeWith (\sd n -> sd / sqrt n) sampleStdDev Window.windowLength  ------------------------------------------------------------------------------- -- Resampling@@ -943,7 +983,7 @@ {-# INLINE jackKnifeMean #-} jackKnifeMean :: (Monad m, Fractional a, Unbox a) => Array a -> Stream m a jackKnifeMean arr = do-    let len = fromIntegral (length arr - 1)+    let len = fromIntegral (Array.length arr - 1)         s = foldArray Fold.sum arr      in fmap (\b -> (s - b) / len) $ Array.read arr @@ -955,10 +995,10 @@ jackKnifeVariance :: (Monad m, Fractional a, Unbox a) =>     Array a -> Stream m a jackKnifeVariance arr = do-    let len = fromIntegral $ length arr - 1+    let len = fromIntegral $ Array.length arr - 1         foldSums (s, s2) x = (s + x, s2 + x ^ (2 :: Int))-        (sum, sum2) = foldArray (Fold.foldl' foldSums (0.0, 0.0)) arr-        var x = (sum2 - x ^ (2 :: Int)) / len -  ((sum - x) / len) ^ (2::Int)+        (sum1, sum2) = foldArray (Fold.foldl' foldSums (0.0, 0.0)) arr+        var x = (sum2 - x ^ (2 :: Int)) / len -  ((sum1 - x) / len) ^ (2::Int)      in fmap var $ Array.read arr  -- | Standard deviation computed from 'jackKnifeVariance'.@@ -981,7 +1021,7 @@      inject arr = liftIO $ do         g <- createSystemRandom-        return $ (g, arr, length arr, 0)+        return $ (g, arr, Array.length arr, 0)      chooseOne g arr len = do         i <- uniformRM (0, len - 1) g

@@ -0,0 +1,787 @@+-- |+-- Module      : Streamly.Statistics.Scanl+-- Copyright   : (c) 2024 Composewell Technologies+-- License     : Apache-2.0+-- Maintainer  : streamly@composewell.com+-- Stability   : experimental+-- Portability : GHC+--+-- See "Streamly.Statistics" for general information. This module provides+-- scans instead of folds.++{-# LANGUAGE ScopedTypeVariables #-}+module Streamly.Statistics.Scanl+    (+    -- * Incremental Scans+    -- | Scans of type @Scanl m (a, Maybe a) b@ are incremental sliding window+    -- scans. An input of type @(a, Nothing)@ indicates that the input element+    -- @a@ is being inserted in the window without ejecting an old value+    -- increasing the window size by 1. An input of type @(a, Just a)@+    -- indicates that the first element is being inserted in the window and the+    -- second element is being removed from the window, the window size remains+    -- the same. The window size can only increase and never decrease.+    --+    -- You can compute the statistics over the entire stream using sliding+    -- window folds by keeping the second element of the input tuple as+    -- @Nothing@.+    --+    -- Also see "Streamly.Data.Scanl" for some basic window scans.++    -- * Summary Statistics+    -- | See https://en.wikipedia.org/wiki/Summary_statistics .++    -- ** Location+    -- | See https://en.wikipedia.org/wiki/Location_parameter .+    --+    -- See https://en.wikipedia.org/wiki/Central_tendency .+      incrMinimum+    , incrMaximum+    , incrRawMoment+    , incrRawMomentFrac++    -- Pythagorean means (https://en.wikipedia.org/wiki/Pythagorean_means)+    , incrWelfordMean+    , incrGeometricMean+    , incrHarmonicMean++    , incrQuadraticMean++    -- Generalized mean+    , incrPowerMean+    , incrPowerMeanFrac++    -- ** Weighted Means+    -- | Exponential Smoothing.+    , ewma+    , ewmaRampUpSmoothing+    , incrEwma++    -- ** Spread+    -- | Second order central moment is a statistical measure of dispersion.+    -- The \(k\)th moment about the mean (or \(k\)th central moment) is defined+    -- as:+    --+    -- \(\mu_k = \frac{1}{n}\sum_{i=1}^n {(x_{i}-\mu)}^k\)+    --+    -- See https://mathworld.wolfram.com/CentralMoment.html .+    --+    -- See https://en.wikipedia.org/wiki/Statistical_dispersion .+    , incrRange+    , incrMd+    , incrVariance+    , incrStdDev++    -- ** Shape+    -- | Third and fourth order central moments are a measure of shape.+    --+    -- See https://en.wikipedia.org/wiki/Shape_parameter .+    --+    -- See https://en.wikipedia.org/wiki/Standardized_moment .+    , incrSkewness+    , incrKurtosis++    -- XXX Move to Statistics.Sample or Statistics.Estimation module?+    -- ** Estimation+    , incrSampleVariance+    , incrSampleStdDev+    , incrStdErrMean++    -- ** Probability Distribution+    , incrFrequency+    )+where++import Control.Exception (assert)+import Control.Monad (when)+import Control.Monad.IO.Class (MonadIO(..))+import Data.Function ((&))+import Data.Map.Strict (Map)+import Data.Maybe (fromMaybe)+import Streamly.Internal.Data.Fold (Step(..))+import Streamly.Internal.Data.Scanl (Scanl(..), Incr(..))+import Streamly.Internal.Data.Tuple.Strict (Tuple'(..), Tuple3'(..))++import qualified Data.Map.Strict as Map+import qualified Deque.Strict as Deque+import qualified Streamly.Data.Fold as Fold+import qualified Streamly.Internal.Data.RingArray as Ring+import qualified Streamly.Internal.Data.Scanl as Scanl+import qualified Streamly.Data.Stream as Stream++import Prelude hiding (length, sum, minimum, maximum)++-- TODO: Overflow checks. Would be good if we can directly replace the+-- operations with overflow checked operations.+--+-- See https://hackage.haskell.org/package/safe-numeric+-- See https://hackage.haskell.org/package/safeint+--+-- TODO We have many of these functions in Streamly.Data.Fold as well. Need to+-- think about deduplication.++-------------------------------------------------------------------------------+-- Location+-------------------------------------------------------------------------------++-- Theoretically, we can approximate minimum in a rolling window by using a+-- 'powerMean' with sufficiently large negative power.+--+-- XXX If we need to know the minimum in the window only once in a while then+-- we can use linear search when it is extracted and not pay the cost all the+-- time.++-- | The minimum element in a rolling window.+--+-- For smaller window sizes (< 30) Streamly.Internal.Data.Fold.windowMinimum performs+-- better.  If you want to compute the minimum of the entire stream Fold.min+-- from streamly package would be much faster.+--+-- /Time/: \(\mathcal{O}(n*w)\) where \(w\) is the window size.+--+{-# INLINE incrMinimum #-}+incrMinimum :: (Monad m, Ord a) => Scanl m (Incr a) a+incrMinimum = Scanl step initial extract extract++    where++    initial =+        return+            $ Partial+            $ Tuple3' (0 :: Int) (0 :: Int) (mempty :: Deque.Deque (Int, a))++    step (Tuple3' i w q) (Insert a) =+                return+                    $ Partial+                    $ Tuple3'+                        (i + 1)+                        (w + 1)+                        (headCheck i q (w + 1) & dqloop (i, a))++    step (Tuple3' i w q) (Replace _ new) =+        return+            $ Partial+            $ Tuple3' (i + 1) w (headCheck i q w & dqloop (i, new))++    {-# INLINE headCheck #-}+    headCheck i q w =+        case Deque.uncons q of+            Nothing -> q+            Just (ia', q') ->+                if fst ia' <= i - w+                then q'+                else q++    dqloop ia q =+        case Deque.unsnoc q of+            Nothing -> Deque.snoc ia q+            -- XXX This can be improved for the case of `=`+            Just (ia', q') ->+                if snd ia <= snd ia'+                then dqloop ia q'+                else Deque.snoc ia q++    extract (Tuple3' _ _ q) =+        return+            $ snd+            $ fromMaybe (0, error "minimum: Empty stream")+            $ Deque.head q++-- Theoretically, we can approximate maximum in a rolling window by using a+-- 'powerMean' with sufficiently large positive power.++-- | The maximum element in a rolling window.+--+-- For smaller window sizes (< 30) Streamly.Internal.Data.Fold.windowMaximum+-- performs better.  If you want to compute the maximum of the entire stream+-- Streamly.Data.Fold.maximum from streamly package would be much faster.+--+-- /Time/: \(\mathcal{O}(n*w)\) where \(w\) is the window size.+--+{-# INLINE incrMaximum #-}+incrMaximum :: (Monad m, Ord a) => Scanl m (Incr a) a+incrMaximum = Scanl step initial extract extract++    where++    initial =+        return+            $ Partial+            $ Tuple3' (0 :: Int) (0 :: Int) (mempty :: Deque.Deque (Int, a))++    step (Tuple3' i w q) (Insert a) =+        return+            $ Partial+            $ Tuple3'+                (i + 1)+                (w + 1)+                (headCheck i q (w + 1) & dqloop (i, a))++    step (Tuple3' i w q) (Replace _ new) =+        return+            $ Partial+            $ Tuple3' (i + 1) w (headCheck i q w & dqloop (i, new))++    {-# INLINE headCheck #-}+    headCheck i q w =+        case Deque.uncons q of+            Nothing -> q+            Just (ia', q') ->+                if fst ia' <= i - w+                then q'+                else q++    dqloop ia q =+        case Deque.unsnoc q of+            Nothing -> Deque.snoc ia q+            -- XXX This can be improved for the case of `=`+            Just (ia', q') ->+                if snd ia >= snd ia'+                then dqloop ia q'+                else Deque.snoc ia q++    extract (Tuple3' _ _ q) =+        return+            $ snd+            $ fromMaybe (0, error "maximum: Empty stream")+            $ Deque.head q++-------------------------------------------------------------------------------+-- Mean+-------------------------------------------------------------------------------++-- | Recompute mean from old mean when an item is added to the sample.+{-# INLINE meanAdd #-}+meanAdd :: Fractional a => Int -> a -> a -> a+meanAdd n oldMean newItem =+    let delta = (newItem - oldMean) / fromIntegral (n + 1)+     in oldMean + delta++-- We do not carry rounding errors, therefore, this would be less numerically+-- stable than the kbn mean.++-- | Recompute mean from old mean when an item in the sample is replaced.+{-# INLINE meanReplace #-}+meanReplace :: Fractional a => Int -> a -> a -> a -> a+meanReplace n oldMean oldItem newItem =+    let n1 = fromIntegral n+        -- Compute two deltas instead of a single (newItem - oldItem) because+        -- the latter would be too small causing rounding errors.+        delta1 = (newItem - oldMean) / n1+        delta2 = (oldItem - oldMean) / n1+     in (oldMean + delta1) - delta2++-- | Same as 'mean' but uses Welford's algorithm to compute the mean+-- incrementally.+--+-- It maintains a running mean instead of a running sum and adjusts the mean+-- based on a new value.  This is slower than 'mean' because of using the+-- division operation on each step and it is numerically unstable (as of now).+-- The advantage over 'mean' could be no overflow if the numbers are large,+-- because we do not maintain a sum, but that is a highly unlikely corner case.+--+-- /Internal/+{-# INLINE incrWelfordMean #-}+incrWelfordMean :: forall m a. (Monad m, Fractional a) => Scanl m (Incr a) a+incrWelfordMean = Scanl step initial extract extract++    where++    initial =+        return+            $ Partial+            $ Tuple'+                (0 :: a)   -- running mean+                (0 :: Int) -- count of items in the window++    step (Tuple' oldMean w) (Insert new) =+        return $ Partial $ Tuple' (meanAdd w oldMean new) (w + 1)++    step (Tuple' oldMean w) (Replace old new) =+        return $ Partial $ Tuple' (meanReplace w oldMean old new) w++    extract (Tuple' x _) = return x++-------------------------------------------------------------------------------+-- Moments+-------------------------------------------------------------------------------++-- XXX We may have chances of overflow if the powers are high or the numbers+-- are large. A limited mitigation could be to use welford style avg+-- computation. Do we need an overflow detection?++-- | Raw moment is the moment about 0. The \(k\)th raw moment is defined as:+--+-- \(\mu'_k = \frac{\sum_{i=1}^n x_{i}^k}{n}\)+--+-- >>> rawMoment k = Fold.teeWith (/) (Fold.windowPowerSum p) Fold.windowLength+--+-- See https://en.wikipedia.org/wiki/Moment_(mathematics) .+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE incrRawMoment #-}+incrRawMoment :: (Monad m, Fractional a) => Int -> Scanl m (Incr a) a+incrRawMoment k =+    Scanl.teeWith (/) (Scanl.incrPowerSum k) Scanl.incrCount++-- | Like 'rawMoment' but powers can be negative or fractional. This is+-- slower than 'rawMoment' for positive intergal powers.+--+-- >>> rawMomentFrac p = Fold.teeWith (/) (Fold.windowPowerSumFrac p) Fold.windowLength+--+{-# INLINE incrRawMomentFrac #-}+incrRawMomentFrac :: (Monad m, Floating a) => a -> Scanl m (Incr a) a+incrRawMomentFrac k =+    Scanl.teeWith (/) (Scanl.incrPowerSumFrac k) Scanl.incrCount++-- XXX Overflow can happen when large powers or large numbers are used. We can+-- keep a running mean instead of running sum but that won't mitigate the+-- overflow possibility by much. The overflow can still happen when computing+-- the mean incrementally.++-- | The \(k\)th power mean of numbers \(x_1, x_2, \ldots, x_n\) is:+--+-- \(M_k = \left( \frac{1}{n} \sum_{i=1}^n x_i^k \right)^{\frac{1}{k}}\)+--+-- \(powerMean(k) = (rawMoment(k))^\frac{1}{k}\)+--+-- >>> powerMean k = (** (1 / fromIntegral k)) <$> rawMoment k+--+-- All other means can be expressed in terms of power mean. It is also known as+-- the generalized mean.+--+-- See https://en.wikipedia.org/wiki/Generalized_mean+--+{-# INLINE incrPowerMean #-}+incrPowerMean :: (Monad m, Floating a) => Int -> Scanl m (Incr a) a+incrPowerMean k = (** (1 / fromIntegral k)) <$> incrRawMoment k++-- | Like 'powerMean' but powers can be negative or fractional. This is+-- slower than 'powerMean' for positive intergal powers.+--+-- >>> powerMeanFrac k = (** (1 / k)) <$> rawMomentFrac k+--+{-# INLINE incrPowerMeanFrac #-}+incrPowerMeanFrac :: (Monad m, Floating a) => a -> Scanl m (Incr a) a+incrPowerMeanFrac k = (** (1 / k)) <$> incrRawMomentFrac k++-- | The harmonic mean of the positive numbers \(x_1, x_2, \ldots, x_n\) is+-- defined as:+--+-- \(HM = \frac{n}{\frac1{x_1} + \frac1{x_2} + \cdots + \frac1{x_n}}\)+--+-- \(HM = \left(\frac{\sum\limits_{i=1}^n x_i^{-1}}{n}\right)^{-1}\)+--+-- >>> harmonicMean = Fold.teeWith (/) length (lmap recip sum)+-- >>> harmonicMean = powerMeanFrac (-1)+--+-- See https://en.wikipedia.org/wiki/Harmonic_mean .+--+{-# INLINE incrHarmonicMean #-}+incrHarmonicMean :: (Monad m, Fractional a) => Scanl m (Incr a) a+incrHarmonicMean =+    Scanl.teeWith (/)+        Scanl.incrCount (Scanl.lmap (fmap recip) Scanl.incrSum)++-- | Geometric mean, defined as:+--+-- \(GM = \sqrt[n]{x_1 x_2 \cdots x_n}\)+--+-- \(GM = \left(\prod_{i=1}^n x_i\right)^\frac{1}{n}\)+--+-- or, equivalently, as the arithmetic mean in log space:+--+-- \(GM = e ^{{\frac{\sum_{i=1}^{n}\ln a_i}{n}}}\)+--+-- >>> geometricMean = exp <$> lmap log mean+--+-- See https://en.wikipedia.org/wiki/Geometric_mean .+{-# INLINE incrGeometricMean #-}+incrGeometricMean :: (Monad m, Floating a) => Scanl m (Incr a) a+incrGeometricMean = exp <$> Scanl.lmap (fmap log) Scanl.incrMean++-- | The quadratic mean or root mean square (rms) of the numbers+-- \(x_1, x_2, \ldots, x_n\) is defined as:+--+-- \(RMS = \sqrt{ \frac{1}{n} \left( x_1^2 + x_2^2 + \cdots + x_n^2 \right) }.\)+--+-- >>> quadraticMean = powerMean 2+--+-- See https://en.wikipedia.org/wiki/Root_mean_square .+--+{-# INLINE incrQuadraticMean #-}+incrQuadraticMean :: (Monad m, Floating a) => Scanl m (Incr a) a+incrQuadraticMean = incrPowerMean 2++-------------------------------------------------------------------------------+-- Weighted Means+-------------------------------------------------------------------------------++-- XXX Is this numerically stable? We can use the kbn summation here.+-- XXX change the signature to use the newer value first.++-- | ewmaStep smoothing-factor old-value new-value+--+-- >>> ewmaStep k x0 x1 = k * x1 + (1 - k) * x0+--+{-# INLINE ewmaStep #-}+ewmaStep :: Double -> Double -> Double -> Double+ewmaStep k x0 x1 = (1 - k) * x0 + k * x1++-- XXX Compute this in a sliding window?++-- | Exponentially weighted moving average is given by @ewma w@ where @w@ is+-- the smoothing factor varying between 0 and 1. To compute the moving average+-- the newest input is given a weight of @w@ and the running ewma is given a+-- weight of @1 - w@.+--+-- For an empty stream this API returns 0. For a non-empty stream the first+-- value in the stream is used as the initial ewma.+--+-- The higher the smoothing factor the more weight is given to the new value.+-- Consider some interesting special cases, when @w@ is 1 the ewma is always+-- the newest value, when @w@ is 0 then the ewma is always the oldest value. If+-- @w@ is @0.5@ then the new inputs get exponentially weighted by weights+-- @1\/2, 1\/4, 1\/8, ...@, see below for details.+--+-- Mathematically, we define ewma \(s_n\), of \(n\) values, \(x_1,\ldots,x_n\),+-- recursively as follows:+--+-- \(\begin{align} s_0& = x_0, \quad \text{initial value}\\ s_n & = \alpha x_{n} + (1-\alpha)s_{n-1},\quad n>0 \end{align}\)+--+-- If we expand the recursive term it reveals an exponential series:+--+-- \(s_n = \alpha \left[x_n + (1-\alpha)x_{n-1} + (1-\alpha)^2 x_{n-2} + \cdots + (1-\alpha)^{n-1} x_1 \right] + (1-\alpha)^n x_0\)+--+-- where \(\alpha\), the smoothing factor, is in the range \(0 <\alpha < 1\).+-- More the value of \(\alpha\), the more weight is given to newer values.  As+-- a special case if it is 0 then the weighted sum would always be the same as+-- the oldest value, if it is 1 then the sum would always be the same as the+-- newest value.+--+-- See https://en.wikipedia.org/wiki/Moving_average+--+-- See https://en.wikipedia.org/wiki/Exponential_smoothing+--+{-# INLINE ewma #-}+ewma :: Monad m => Double -> Scanl m Double Double+ewma k = extract <$> Scanl.mkScanl step (Tuple' 0 1)++    where++    step (Tuple' x0 k1) x = Tuple' (ewmaStep k1 x0 x) k++    extract (Tuple' x _) = x++-- | @ewma n k@ is like 'ewma' but the smoothing factor used is itself+-- exponentially smoothened starting from @1@ to the final value @k@ using @n@+-- as its smoothing factor. In other words, the smoothing factor is derived by+-- smoothing the series @1,k,k,k,...@ using @n@ as the smoothing factor. As+-- time goes on, the smoothing factor gets closer to k.+--+{-# INLINE ewmaRampUpSmoothing #-}+ewmaRampUpSmoothing :: Monad m => Double -> Double -> Scanl m Double Double+ewmaRampUpSmoothing n k1 = extract <$> Scanl.mkScanl step initial++    where++    initial = Tuple' 0 1++    step (Tuple' x0 k0) x1 =+        let x = ewmaStep k0 x0 x1+            k = ewmaStep n k0 k1+        in Tuple' x k++    extract (Tuple' x _) = x++data Ewma = EwmaInit | EwmaGo !Double !Double++-- | @incrEwma w@ computes the ewma of the elements incrementally in a window+-- using @w@ as the smoothing factor.+--+-- ewma can be calculated incrementally as follows. If \(x_i\) is the value+-- entering the window, n is the size of window, \(x_1\) is the second last value+-- in the old window and \(x_0\) is the value exiting the window:+--+-- \(s_{new} = \alpha x_i + (1-\alpha)s_{old} + (1-\alpha)^{n} (x_1 - x_0\)\)+--+-- /Unimplemented/+--+{-# INLINE incrEwma #-}+incrEwma :: MonadIO m =>+    Double -> Scanl m (Incr Double, Ring.RingArray Double) Double+incrEwma w = Scanl step initial extract extract++    where++    initial = do+        when (w < 0 || w > 1)+            $ error "incrEwma: weight must be >= 0 and <= 1"+        return $ Partial EwmaInit++    step EwmaInit (Insert x, rng) = do+        let len = Ring.length rng+        assert (len == 0) (return ())+        return $ Partial (EwmaGo x 1)++    step EwmaInit _ = error "incrEwma: the first window operation must be Insert"++    step (EwmaGo s k) (Insert x, _) = do+        let s1 = w * x + (1 - w) * s+        return $ Partial (EwmaGo s1 (k * (1 - w)))++    step (EwmaGo s k) (Replace x0 x, rng) = do+        x1 <- Ring.unsafeGetIndex 0 rng+        let s1 = w * x + (1 - w) * s + k * (x1 - x0)+        return $ Partial (EwmaGo s1 k)++    extract EwmaInit = return 0+    extract (EwmaGo x _) = return x++-------------------------------------------------------------------------------+-- Spread/Dispersion+-------------------------------------------------------------------------------++-- | The difference between the maximum and minimum elements of a rolling window.+--+-- >>> range = Fold.teeWith (-) maximum minimum+--+-- If you want to compute the range of the entire stream @Fold.teeWith (-)+-- Fold.maximum Fold.minimum@ from the streamly package would be much faster.+--+-- /Space/: \(\mathcal{O}(n)\) where @n@ is the window size.+--+-- /Time/: \(\mathcal{O}(n*w)\) where \(w\) is the window size.+--+{-# INLINE incrRange #-}+incrRange :: (Monad m, Num a, Ord a) => Scanl m (Incr a) a+incrRange = Scanl.teeWith (-) incrMaximum incrMinimum++-- | @md n@ computes the mean absolute deviation (or mean deviation) in a+-- sliding window of last @n@ elements in the stream.+--+-- The input of the scan is (incr, ring), where incr is the incremental window+-- operation and ring is the contents of the entire window in a ring array.+--+-- The mean absolute deviation of the numbers \(x_1, x_2, \ldots, x_n\) is:+--+-- \(MD = \frac{1}{n}\sum_{i=1}^n |x_i-\mu|\)+--+-- Note: It is expensive to compute MD in a sliding window. We need to+-- maintain a ring buffer of last n elements and maintain a running mean, when+-- the result is extracted we need to compute the difference of all elements+-- from the mean and get the average. Using standard deviation may be+-- computationally cheaper.+--+-- See https://en.wikipedia.org/wiki/Average_absolute_deviation .+--+-- /Pre-release/+{-# INLINE incrMd #-}+incrMd ::  MonadIO m =>+    Scanl m (Incr Double, Ring.RingArray Double) Double+incrMd =+    Scanl.rmapM computeMD+        $ Scanl.tee+            (Scanl.lmap fst Scanl.incrMean) (Scanl.lmap snd Scanl.latest)++    where++    computeMD (mn, mRng) =+        case mRng of+            Just rng -> do+                Stream.fold Fold.mean+                    $ fmap (\a -> abs (mn - a))+                    $ Ring.read rng+            Nothing -> return 0.0++-- | The variance \(\sigma^2\) of a population of \(n\) equally likely values+-- is defined as the average of the squares of deviations from the mean+-- \(\mu\). In other words, second moment about the mean:+--+-- \(\sigma^2 = \frac{1}{n}\sum_{i=1}^n {(x_{i}-\mu)}^2\)+--+-- \(\sigma^2 = rawMoment(2) - \mu^2\)+--+-- \(\mu_2 = -(\mu'_1)^2 + \mu'_2\)+--+-- Note that the variance would be biased if applied to estimate the population+-- variance from a sample of the population. See 'sampleVariance'.+--+-- See https://en.wikipedia.org/wiki/Variance.+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE incrVariance #-}+incrVariance :: (Monad m, Fractional a) => Scanl m (Incr a) a+incrVariance =+    Scanl.teeWith+        (\p2 m -> p2 - m ^ (2 :: Int)) (incrRawMoment 2) Scanl.incrMean++-- | Standard deviation \(\sigma\) is the square root of 'variance'.+--+-- This is the population standard deviation or uncorrected sample standard+-- deviation.+--+-- >>> stdDev = sqrt <$> variance+--+-- See https://en.wikipedia.org/wiki/Standard_deviation .+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE incrStdDev #-}+incrStdDev :: (Monad m, Floating a) => Scanl m (Incr a) a+incrStdDev = sqrt <$> incrVariance++-- XXX Need a tee3 operation for better performance.++-- | Skewness \(\gamma\) is the standardized third central moment defined as:+--+-- \(\tilde{\mu}_3 = \frac{\mu_3}{\sigma^3}\)+--+-- The third central moment can be computed in terms of raw moments:+--+-- \(\mu_3 = 2(\mu'_1)^3 - 3\mu'_1\mu'_2 + \mu'_3\)+--+-- Substituting \(\mu'_1 = \mu\), and \(\mu'_2 = \mu^2 + \sigma^2\):+--+-- \(\mu_3 = -\mu^3 - 3\mu\sigma^2 + \mu'_3\)+--+-- Skewness is a measure of symmetry of the probability distribution. It is 0+-- for a symmetric distribution, negative for a distribution that is skewed+-- towards left, positive for a distribution skewed towards right.+--+-- For a normal like distribution the median can be found around+-- \(\mu - \frac{\gamma\sigma}{6}\) and the mode can be found around+-- \(\mu - \frac{\gamma \sigma}{2}\).+--+-- See https://en.wikipedia.org/wiki/Skewness .+--+{-# INLINE incrSkewness #-}+incrSkewness :: (Monad m, Floating a) => Scanl m (Incr a) a+incrSkewness =+          (\rm3 sd mu ->+            rm3 / sd ^ (3 :: Int) - 3 * (mu / sd) - (mu / sd) ^ (3 :: Int)+          )+        <$> incrRawMoment 3+        <*> incrStdDev+        <*> Scanl.incrMean++-- XXX We can compute the 2nd, 3rd, 4th raw moments by repeatedly multiplying+-- instead of computing the powers every time.+-- XXX Need a tee4 operation for better performance.++-- | Kurtosis \(\kappa\) is the standardized fourth central moment, defined as:+--+-- \(\tilde{\mu}_4 = \frac{\mu_4}{\sigma^4}\)+--+-- The fourth central moment can be computed in terms of raw moments:+--+-- \(\mu_4 = -3(\mu'_1)^4 + 6(\mu'_1)^2\mu'_2 - 4\mu'_1\mu'_3\ + \mu'_4\)+--+-- Substituting \(\mu'_1 = \mu\), and \(\mu'_2 = \mu^2 + \sigma^2\):+--+-- \(\mu_4 = 3\mu^4 + 6\mu^2\sigma^2 - 4\mu\mu'_3 + \mu'_4\)+--+-- It is always non-negative. It is 0 for a point distribution, low for light+-- tailed (platykurtic) distributions and high for heavy tailed (leptokurtic)+-- distributions.+--+-- \(\kappa >= \gamma^2 + 1\)+--+-- For a normal distribution \(\kappa = 3\sigma^4\).+--+-- See https://en.wikipedia.org/wiki/Kurtosis .+--+{-# INLINE incrKurtosis #-}+incrKurtosis :: (Monad m, Floating a) => Scanl m (Incr a) a+incrKurtosis =+          (\rm4 rm3 sd mu ->+             ( 3 * mu ^ (4 :: Int)+            + 6 * mu ^ (2 :: Int) * sd ^ (2 :: Int)+            - 4 * mu * rm3+            + rm4) / (sd ^ (4 :: Int))+          )+        <$> incrRawMoment 4+        <*> incrRawMoment 3+        <*> incrStdDev+        <*> Scanl.incrMean++-------------------------------------------------------------------------------+-- Estimation+-------------------------------------------------------------------------------++-- | Unbiased sample variance i.e. the variance of a sample corrected to+-- better estimate the variance of the population, defined as:+--+-- \(s^2 = \frac{1}{n - 1}\sum_{i=1}^n {(x_{i}-\mu)}^2\)+--+-- \(s^2 = \frac{n}{n - 1} \times \sigma^2\).+--+-- See https://en.wikipedia.org/wiki/Bessel%27s_correction.+--+{-# INLINE incrSampleVariance #-}+incrSampleVariance :: (Monad m, Fractional a) => Scanl m (Incr a) a+incrSampleVariance =+    Scanl.teeWith (\n s2 -> n * s2 / (n - 1)) Scanl.incrCount incrVariance++-- | Sample standard deviation:+--+-- \(s = \sqrt{sampleVariance}\)+--+-- >>> sampleStdDev = sqrt <$> sampleVariance+--+-- See https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation+-- .+--+{-# INLINE incrSampleStdDev #-}+incrSampleStdDev :: (Monad m, Floating a) => Scanl m (Incr a) a+incrSampleStdDev = sqrt <$> incrSampleVariance++-- | Standard error of the sample mean (SEM), defined as:+--+-- \( SEM = \frac{sampleStdDev}{\sqrt{n}} \)+--+-- See https://en.wikipedia.org/wiki/Standard_error .+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE incrStdErrMean #-}+incrStdErrMean :: (Monad m, Floating a) => Scanl m (Incr a) a+incrStdErrMean =+    Scanl.teeWith (\sd n -> sd / sqrt n) incrSampleStdDev Scanl.incrCount++-------------------------------------------------------------------------------+-- Probability Distribution+-------------------------------------------------------------------------------++-- XXX We can use a Windowed classifyWith operation, that will allow us to+-- express windowed frequency, mode, histograms etc idiomatically.++-- | Count the frequency of elements in a sliding window.+--+-- >>> input = Stream.fromList [1,1,3,4,4::Int]+-- >>> f = Ring.slidingWindow 4 Statistics.frequency+-- >>> Stream.fold f input+-- fromList [(1,1),(3,1),(4,2)]+--+{-# INLINE incrFrequency #-}+incrFrequency :: (Monad m, Ord a) => Scanl m (Incr a) (Map a Int)+incrFrequency = Scanl.mkScanl step Map.empty++    where++    decrement v =+        if v == 1+        then Nothing+        else Just (v - 1)++    step refCountMap (Insert new) =+        Map.insertWith (+) new 1 refCountMap++    step refCountMap (Replace old new) =+        let m1 = Map.insertWith (+) new 1 refCountMap+         in Map.update decrement old m1

@@ -1,6 +1,6 @@ cabal-version:       2.4 name:                streamly-statistics-version:             0.1.0+version:             0.2.0 synopsis:     Statistical measures for finite or infinite data streams. description:@@ -25,14 +25,18 @@ license:             Apache-2.0 license-file:        LICENSE tested-with:-      GHC==8.10.7-    , GHC==9.0.2-    , GHC==9.2.2+      GHC==8.8.4+    , GHC==8.10.7+    , GHC==9.2.7     , GHC==9.4.4+    , GHC==9.6.3+    , GHC==9.8.1+    , GHC==9.10.1+    , GHC==9.12.1 author: Composewell Technologies maintainer: streamly@composewell.com copyright: 2019 Composewell Technologies-category: Streamly, Statistics+category: Streamly, Statistics, Streaming  extra-source-files:     CHANGELOG.md@@ -93,6 +97,7 @@                  -Wincomplete-uni-patterns                  -Wredundant-constraints                  -Wnoncanonical-monad-instances+                 -Wno-deprecations                  -Rghc-timing  common optimization-options@@ -107,10 +112,12 @@  library     import: ghc-options-    exposed-modules:     Streamly.Statistics+    exposed-modules:+        Streamly.Statistics+      , Streamly.Statistics.Scanl     build-depends:       base     >= 4.9 && < 5-                       , streamly-core == 0.1.0-                       , containers  >= 0.5   && < 0.7+                       , streamly-core >= 0.3.0 && < 0.4.0+                       , containers  >= 0.5   && < 0.8                        , random >= 1.2 && < 1.3                        , mwc-random >= 0.15 && < 0.16                        , deque      >= 0.4.4 && < 0.4.5@@ -122,13 +129,13 @@     hs-source-dirs:     test     main-is:            Main.hs     build-depends:      streamly-statistics-                      , streamly-core == 0.1.0+                      , streamly-core >= 0.3.0                       , base           >= 4.9   && < 5-                      , QuickCheck     >= 2.10  && < 2.15+                      , QuickCheck     >= 2.10  && < 2.16                       , hspec          >= 2.0   && < 3                       , hspec-core     >= 2.0   && < 3                       , random         >= 1.0.0 && < 2-                      , containers     >= 0.5   && < 0.7+                      , containers     >= 0.5   && < 0.8                       -- XXX Should remove these dependencies                       , vector         >= 0.11  && < 0.14                       , statistics     >= 0.15  && < 0.17@@ -140,12 +147,12 @@     hs-source-dirs:   benchmark     main-is:          Main.hs     build-depends:      streamly-statistics-                      , streamly-core == 0.1.0+                      , streamly-core >= 0.3.0                       , base           >= 4.9   && < 5                       , random         >= 1.0.0 && < 2-                      , deepseq        >= 1.4.1 && < 1.5-                      , tasty-bench >= 0.3 && < 0.4-                      , tasty     >= 1.4.1 && < 1.5+                      , deepseq        >= 1.4.1 && < 1.6+                      , tasty-bench >= 0.3 && < 0.5+                      , tasty     >= 1.4.1 && < 1.6     mixins: tasty-bench       (Test.Tasty.Bench as Gauge       ,Test.Tasty.Bench as Gauge.Main

@@ -5,6 +5,7 @@ import Data.Functor.Classes (liftEq2) import Streamly.Data.Array (Unbox) import Streamly.Data.Stream (Stream)+import Streamly.Internal.Data.Scanl (Incr(..)) import Test.Hspec.Core.Spec (SpecM) import Test.Hspec.QuickCheck (prop) import Test.QuickCheck@@ -19,9 +20,11 @@ import qualified Streamly.Data.Array as Array import qualified Streamly.Data.Fold as Fold import qualified Streamly.Data.MutArray as MA-import qualified Streamly.Internal.Data.Ring.Unboxed as Ring+import qualified Streamly.Internal.Data.RingArray as Ring import qualified Streamly.Data.Stream as Stream import qualified Streamly.Data.Stream as S+import qualified Streamly.Internal.Data.Scanl as Scanl+import qualified Streamly.Statistics.Scanl as Stat  import Prelude hiding (sum, maximum, minimum) @@ -200,19 +203,30 @@ main :: IO () main = hspec $ do     describe "Numerical stability while streaming" $ do-        let numElem = 80000-            winSize = 800-            testCaseChunk = [9007199254740992, 1, 1.0 :: Double,-                                9007199254740992, 1, 1, 1, 9007199254740992]-            testCase = take numElem $ cycle testCaseChunk+        let winSize = 800+            numElem = winSize * 100+            segment =+                [ 9007199254740992+                , 1+                , 1.0 :: Double+                , 9007199254740992+                , 1+                , 1+                , 1+                , 9007199254740992+                ]+            input = take numElem $ cycle segment             deviationLimit = 1             testFunc f = do-                let c = S.fromList testCase+                let c = S.fromList input                 a <- runIO $ S.fold (Ring.slidingWindow winSize f) c-                b <- runIO $ S.fold f $ S.drop (numElem - winSize)-                        $ fmap (, Nothing) c+                b <- runIO $ S.fold f $ S.take winSize $ fmap (, Nothing) c                 let c1 = a - b-                it ("should not deviate more than " ++ show deviationLimit)+                it ("deviation " ++ show c1 ++ " should not be more than "+                    ++ show deviationLimit+                    ++ " one window fold = " ++ show b+                    ++ " rolling window fold = " ++ show a+                    )                     $ c1 >= -1 * deviationLimit && c1 <= deviationLimit          describe "Sum" $ testFunc sum@@ -226,9 +240,13 @@              testFunc tc f sI sW = do                 let c = S.fromList tc-                a <- runIO $ S.fold Fold.toList $ S.postscan f $ fmap (, Nothing) c-                b <- runIO $ S.fold Fold.toList $ S.postscan-                        (Ring.slidingWindow winSize f) c+                a <- runIO+                        $ S.fold Fold.toList+                        $ S.postscanl f+                        $ fmap Insert c+                b <- runIO+                        $ S.fold Fold.toList+                        $ S.postscanl (Scanl.incrScan winSize f) c                 it "Infinite" $ a  == sI                 it ("Finite " ++ show winSize) $ b == sW @@ -254,27 +272,27 @@         describe "minimum" $ do             let scanInf = [31, 31, 31, 26, 26, 26, 26] :: [Double]                 scanWin = [31, 31, 31, 26, 26, 26, 53] :: [Double]-            testFunc testCase1 minimum scanInf scanWin+            testFunc testCase1 Stat.incrMinimum scanInf scanWin         describe "maximum" $ do             let scanInf = [31, 41, 59, 59, 59, 59, 97] :: [Double]                 scanWin = [31, 41, 59, 59, 59, 58, 97] :: [Double]-            testFunc testCase1 maximum scanInf scanWin+            testFunc testCase1 Stat.incrMaximum scanInf scanWin         describe "range" $ do             let scanInf = [0, 10, 28, 33, 33, 33, 71] :: [Double]                 scanWin = [0, 10, 28, 33, 33, 32, 44] :: [Double]-            testFunc testCase1 range scanInf scanWin+            testFunc testCase1 Stat.incrRange scanInf scanWin         describe "sum" $ do             let scanInf = [1, 2, 3, 4, 5, 12] :: [Double]                 scanWin = [1, 2, 3, 3, 3, 9] :: [Double]-            testFunc testCase2 sum scanInf scanWin+            testFunc testCase2 Scanl.incrSum scanInf scanWin         describe "mean" $ do             let scanInf = [1, 1, 1, 1, 1, 2] :: [Double]                 scanWin = [1, 1, 1, 1, 1, 3] :: [Double]-            testFunc testCase2 mean scanInf scanWin+            testFunc testCase2 Scanl.incrMean scanInf scanWin         describe "welfordMean" $ do             let scanInf = [1, 1, 1, 1, 1, 2] :: [Double]                 scanWin = [1, 1, 1, 1, 1, 3] :: [Double]-            testFunc testCase2 welfordMean scanInf scanWin+            testFunc testCase2 Stat.incrWelfordMean scanInf scanWin          -- Probability Distribution         describe "frequency"