diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/benchmark/Main.hs b/benchmark/Main.hs --- a/benchmark/Main.hs +++ b/benchmark/Main.hs @@ -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 ] ] diff --git a/src/Streamly/Statistics.hs b/src/Streamly/Statistics.hs --- a/src/Streamly/Statistics.hs +++ b/src/Streamly/Statistics.hs @@ -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 diff --git a/src/Streamly/Statistics/Scanl.hs b/src/Streamly/Statistics/Scanl.hs new file mode 100644 --- /dev/null +++ b/src/Streamly/Statistics/Scanl.hs @@ -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 diff --git a/streamly-statistics.cabal b/streamly-statistics.cabal --- a/streamly-statistics.cabal +++ b/streamly-statistics.cabal @@ -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 diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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"