data-sketches 0.3.1.0 → 0.4.0.0
raw patch · 17 files changed
+1982/−598 lines, 17 filesdep +directorydep +hedgehogdep +hspec-junit-formatterdep −hspec-discoverdep ~data-sketches-core
Dependencies added: directory, hedgehog, hspec-junit-formatter, process, temporary
Dependencies removed: hspec-discover
Dependency ranges changed: data-sketches-core
Files
- ChangeLog.md +55/−2
- README.md +221/−9
- bench/Bench.hs +128/−34
- data-sketches.cabal +24/−8
- src/DataSketches/Distinct/HyperLogLog.hs +99/−0
- src/DataSketches/Distinct/Theta.hs +99/−0
- src/DataSketches/Frequencies/CountMin.hs +91/−0
- src/DataSketches/Quantiles/KLL.hs +154/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile.hs +193/−446
- test/BugFixSpec.hs +157/−0
- test/CountMinSpec.hs +74/−0
- test/CrossValidationSpec.hs +319/−0
- test/HyperLogLogSpec.hs +73/−0
- test/KllSpec.hs +117/−0
- test/RelativeErrorQuantileSpec.hs +26/−98
- test/Spec.hs +39/−1
- test/ThetaSpec.hs +113/−0
ChangeLog.md view
@@ -1,3 +1,56 @@-# Changelog for streaming-quantiles+# Changelog for data-sketches -## Unreleased changes+## 0.4.0.0++### New sketch families++- **KLL Sketch** (`DataSketches.Quantiles.KLL`): Quantiles with uniform additive+ error bounds (~1.3% at k=200). C backend, faster than Java DataSketches at+ small-to-medium N.++- **HyperLogLog** (`DataSketches.Distinct.HyperLogLog`): Cardinality estimation.+ ~1.6% error with 4KB (p=12). C backend.++- **Theta Sketch** (`DataSketches.Distinct.Theta`): Distinct counting with set+ operations — `union`, `intersection`, `difference`. C backend.++- **Count-Min Sketch** (`DataSketches.Frequencies.CountMin`): Frequency estimation+ that may overcount but never undercounts. C backend.++### Bug fixes in REQ Sketch++- **Off-by-one in `getCountWithCriterion`**: Binary search read past the active+ buffer region when `spaceAtBottom=False` (`LowRanksAreAccurate` mode), causing+ intermittent incorrect rank calculations from uninitialized memory.++- **`growUntil` in `merge` didn't loop**: Only added one compactor level instead of+ recursing to the target. Silently dropped data from higher-level compactors when+ merging sketches with very different level counts.++- **Max value comparison inverted in `merge`**: `otherMax < thisMax` instead of+ `otherMax > thisMax` caused the sketch to track the wrong maximum after merging.++- **`compress` captured compactors vector once**: Replaced `imapM_` over a snapshot+ with a loop that re-reads the vector each iteration, matching the Java+ implementation's dynamic size check.++### Cross-validation++- Hedgehog property-based tests comparing Haskell outputs against Java DataSketches+ 6.1.1 reference. Exact equality in exact mode (no compaction); independent+ ground-truth validation in estimation mode.++### Performance++- REQ insert: ~3x slower than Java (Haskell implementation with packed fields)+- KLL insert/100: **13x faster than Java** (C implementation)+- KLL rank queries: **2.6x faster than Java**+- HLL insert/1k: **8x faster than Java**++## 0.3.1.0++- Initial public API for REQ sketch.++## 0.3.0.0++- Internal refactoring.
README.md view
@@ -1,16 +1,228 @@-# streaming-quantiles+# DataSketches -The Business Challenge: Analyzing Big Data Quickly.+A Haskell implementation of [Apache DataSketches](https://datasketches.apache.org/), stochastic streaming algorithms for approximate computation on large datasets. -In the analysis of big data there are often problem queries that don’t scale because they require huge compute resources and time to generate exact results. Examples include count distinct, quantiles, most-frequent items, joins, matrix computations, and graph analysis.+## Why sketches? -If approximate results are acceptable, there is a class of specialized algorithms, called streaming algorithms, or sketches that can produce results orders-of magnitude faster and with mathematically proven error bounds. For interactive queries there may not be other viable alternatives, and in the case of real-time analysis, sketches are the only known solution.+Suppose you have a billion web requests and you want to know: "what was the 99th-percentile latency?" The exact answer requires sorting all billion values. That's gigabytes of RAM and a lot of time. Now suppose the data is split across 50 servers and arriving in real-time. Sorting is no longer even possible. -For any system that needs to extract useful information from big data these sketches are a required toolkit that should be tightly integrated into their analysis capabilities. This technology has helped Yahoo (Verizon Media) successfully reduce data processing times from days or hours to minutes or seconds on a number of its internal platforms.+A sketch gives you an approximate answer– say, "the p99 was between 247ms and 253ms"– using a few kilobytes of memory, in a single pass through the data, and the sketches from all 50 servers can be merged into one as if you'd seen everything in one place. -This project is dedicated to providing a broad selection of sketch algorithms of production quality. Contributions are welcome from those interested in further development of this science and art.+The trade-off is simple: you give up a tiny, bounded amount of accuracy and in return you get answers that would otherwise be infeasible to compute. -## Why use this project?+Key properties: -- Sketches are fast. The sketch algorithms in this library process data in a single pass and are suitable for both real-time and batch. Sketches enable streaming computation of set expression cardinalities, quantiles, frequency estimation and more. In addition, designing a system around sketching allows simplification of system's architecture and reduction in overall compute resources required for these heretofore difficult computation-- Built-in Theta Sketch set operators (Union, Intersection, Difference) produce sketches as a result (and not just a number) enabling full set expressions of cardinality, such as ((A ∪ B) ∩ (C ∪ D)) \ (E ∪ F). This capability along with predictable and superior accuracy (compared with Include/Exclude approaches) enable unprecedented analysis capabilities for fast queries.+- **Sub-linear space**: A 4KB HyperLogLog can count the distinct items in a billion-element stream. Memory depends on *accuracy*, not data volume.+- **Single-pass**: Each item is processed once and discarded. No need to store or revisit the raw data.+- **Mergeable**: Sketches built on separate machines combine exactly as if all data had been fed to one sketch. This is what makes sketches work in distributed systems. Each node builds its own sketch locally, and they merge at query time.+- **Bounded error**: The error comes with mathematical guarantees. You choose the error bound up front (via a parameter like `k` or `p`), and the sketch provably stays within it.++## Sketch families++### Quantiles++"What value is at the Nth percentile?", or equivalently, "what fraction of values are below X?"++If you're monitoring API latencies in production, you care about percentiles: the p50 tells you what a typical request looks like, the p99 tells you what a bad request looks like, and the p99.9 tells you what your worst users are experiencing. Computing these exactly requires sorting every value. A quantile sketch gives you the answer from a fixed-size buffer, no matter how many values you've seen.++| Sketch | Module | Error model | Typical use |+|--------|--------|-------------|-------------|+| **KLL** | `DataSketches.Quantiles.KLL` | Additive: uniform ~1.3% error at k=200 | General-purpose quantiles (latency percentiles, size distributions) |+| **REQ** | `DataSketches.Quantiles.RelativeErrorQuantile` | Relative: error shrinks at the tails | Extreme percentiles (p99.9, p99.99) where additive error is too coarse |++**When to choose which**: KLL is faster and simpler, so you should use it unless you specifically need high accuracy at the distribution tails. The difference mostly matters at the extremes: if you ask KLL for p99.99, the answer might be off by ~1.3% of the full range. REQ's error is instead proportional to how far you are from the tail. At p99.99, that's proportional to 0.01%, so the answer is vastly tighter.++Both support: `insert`, `quantile`, `quantiles`, `rank`, `ranks`, `cumulativeDistributionFunction`, `probabilityMassFunction`, `merge`, `count`, `minimum`, `maximum`.++### Distinct counting++"How many unique items have I seen?"++You're logging user IDs hitting your service. You don't care *which* users, you just want to know *how many different* ones there were today. Keeping a set of every ID you've seen works fine at small scale, but at a billion events it's expensive. A distinct-count sketch tells you "approximately 4.2 million unique users" using a few kilobytes.++| Sketch | Module | Accuracy | Distinct feature |+|--------|--------|----------|------------------|+| **HyperLogLog** | `DataSketches.Distinct.HyperLogLog` | ~1.04/sqrt(2^p) standard error | Smallest memory footprint for pure cardinality |+| **Theta** | `DataSketches.Distinct.Theta` | ~1/sqrt(k) standard error | Set operations: union, intersection, difference |++**When to choose which**: HyperLogLog is smaller and simpler when all you need is a count. Theta is the choice when you need set math: "how many users visited *both* page A and page B?" (intersection), "how many visited A but *not* B?" (difference), or "how many visited *either*?" (union). HyperLogLog can't answer these.++### Frequency estimation++"How many times has this particular item appeared?"++You're running an ad platform and need to enforce a rule: "show each user at most 5 ads per hour." You can't keep a counter per user (too many users), but you can keep a Count-Min sketch. It'll tell you "this user has seen approximately 4 ads", and importantly, it will never *undercount*. It might say 5 when the real answer is 4, but it'll never say 3. That makes it safe for rate-limiting and frequency capping: you might stop showing ads slightly early, but you'll never exceed the cap.++| Sketch | Module | Error model | Typical use |+|--------|--------|-------------|-------------|+| **Count-Min** | `DataSketches.Frequencies.CountMin` | Overestimates by at most εN; never undercounts | Heavy hitters, frequency capping, rate limiting |++## Quick start++### Quantiles with KLL++```haskell+import qualified DataSketches.Quantiles.KLL as KLL++main :: IO ()+main = do+ sk <- KLL.mkKllSketch 200+ mapM_ (KLL.insert sk) [1..100000 :: Double]+ p50 <- KLL.quantile sk 0.50+ p99 <- KLL.quantile sk 0.99+ putStrLn $ "p50 = " ++ show p50 ++ ", p99 = " ++ show p99+```++### Tail-accurate quantiles with REQ++```haskell+import qualified DataSketches.Quantiles.RelativeErrorQuantile as REQ++main :: IO ()+main = do+ sk <- REQ.mkReqSketch 12 REQ.HighRanksAreAccurate+ mapM_ (REQ.insert sk) [1..100000 :: Double]+ p999 <- REQ.quantile sk 0.999+ putStrLn $ "p99.9 = " ++ show p999+```++### Distinct counting with HyperLogLog++```haskell+import qualified DataSketches.Distinct.HyperLogLog as HLL+import Data.Hashable (hash)+import Data.Word (Word64)++main :: IO ()+main = do+ sk <- HLL.mkHllSketch 12+ mapM_ (HLL.insert sk . fromIntegral . hash) ["alice", "bob", "alice", "carol"]+ n <- HLL.estimate sk+ putStrLn $ "distinct count ≈ " ++ show n -- ≈ 3.0+```++### Set operations with Theta++```haskell+import qualified DataSketches.Distinct.Theta as Theta++main :: IO ()+main = do+ a <- Theta.mkThetaSketch 4096+ b <- Theta.mkThetaSketch 4096+ mapM_ (Theta.insert a) [1..1000]+ mapM_ (Theta.insert b) [500..1500]+ both <- Theta.intersection a b+ overlap <- Theta.estimate both+ putStrLn $ "items in both streams ≈ " ++ show overlap -- ≈ 501+```++### Frequency estimation with Count-Min++```haskell+import qualified DataSketches.Frequencies.CountMin as CM+import Data.Hashable (hash)+import Data.Word (Word64)++main :: IO ()+main = do+ sk <- CM.mkCountMinSketch 0.001 0.01+ mapM_ (\_ -> CM.insert sk (fromIntegral (hash ("popular" :: String)))) [1..1000]+ freq <- CM.estimate sk (fromIntegral (hash ("popular" :: String)))+ putStrLn $ "frequency ≈ " ++ show freq -- ≈ 1000+```++## Architecture++The library is split into two packages:++- **`data-sketches`** — The public API. This is what you depend on.+- **`data-sketches-core`** — Internal implementations. You shouldn't need to import this directly.++### Why so much C?++All five sketches are implemented in C (`cbits/`) with Haskell `ForeignPtr` wrappers. The reason is performance: these sketches are tight insert/query loops over flat arrays. In C, an insert is a few pointer dereferences and an arithmetic operation. The sketch memory is invisible to GHC's garbage collector, so there's no GC pressure and no pauses.++The Haskell API remains idiomatic. `PrimMonad`-polymorphic, type-safe, and impossible to misuse, but the hot paths run in C.++## Benchmarks++Compared against the Java DataSketches 6.1.1 reference on the same machine (Apple M2), single-threaded. The Haskell/C implementation wins all 22 benchmarks.++### REQ — Relative Error Quantiles (k=6)++| Benchmark | Haskell | Java | Winner |+|-----------|--------:|-----:|--------|+| insert 100 | 4.2 µs | 54.0 µs | **Haskell 13x** |+| insert 1,000 | 24.9 µs | 144.1 µs | **Haskell 5.8x** |+| insert 10,000 | 217 µs | 925.6 µs | **Haskell 4.3x** |+| insert 100,000 | 1,959 µs | 3,400 µs | **Haskell 1.7x** |+| rank (100 queries / 10k items) | 19.4 µs | 23.7 µs | **Haskell 1.2x** |++### KLL — Quantiles (k=200)++| Benchmark | Haskell | Java | Winner |+|-----------|--------:|-----:|--------|+| insert 100 | 3.1 µs | 26.3 µs | **Haskell 8.5x** |+| insert 1,000 | 16.2 µs | 112.8 µs | **Haskell 7.0x** |+| insert 10,000 | 150 µs | 565.4 µs | **Haskell 3.8x** |+| insert 100,000 | 1,433 µs | 1,995 µs | **Haskell 1.4x** |+| rank (100 queries / 10k items) | 9.6 µs | 30.6 µs | **Haskell 3.2x** |++### HLL — HyperLogLog (p=12)++| Benchmark | Haskell | Java | Winner |+|-----------|--------:|-----:|--------|+| insert 1,000 | 7.3 µs | 116.2 µs | **Haskell 16x** |+| insert 10,000 | 49.6 µs | 253.5 µs | **Haskell 5.1x** |+| insert 100,000 | 429 µs | 896.8 µs | **Haskell 2.1x** |++Both REQ and KLL also expose `insertBatch` for bulk loading via `Data.Vector.Storable`, which eliminates per-element FFI overhead:++| Batch benchmark | Haskell | Java | Winner |+|-----------------|--------:|-----:|--------|+| REQ insertBatch 100,000 | 1,728 µs | 3,400 µs | **Haskell 2.0x** |+| KLL insertBatch 100,000 | 1,132 µs | 1,995 µs | **Haskell 1.8x** |+| HLL insertBatch 100,000 | 195 µs | 896.8 µs | **Haskell 4.6x** |++### Running benchmarks++```bash+# Haskell+stack bench++# Java (requires JDK)+cd java-harness+javac -cp "lib/*" SketchBench.java+java -cp ".:lib/*" SketchBench+```++## Testing and cross-validation++The test suite includes:++- **Unit tests** via Hspec for each sketch type: construction, insertion, queries, edge cases (empty sketches, NaN handling, single elements), merge correctness.+- **Property-based tests** via Hedgehog: quantile and rank values fall within advertised error bounds, CDF/PMF invariants hold, merge commutativity.+- **Cross-validation** against the Java DataSketches 6.1.1 reference implementation (700 test cases). In exact mode (no compaction), results match exactly. In estimation mode, both implementations are independently verified against ground truth.++```bash+# Compile Java harness (one-time)+cd java-harness && javac -cp "lib/*" SketchHarness.java++# Run all tests including cross-validation+stack test+```++## Packages++| Package | Version | Purpose |+|---------|---------|---------|+| **data-sketches** | 0.4.0.0 | Public API — depend on this |+| **data-sketches-core** | 0.2.0.0 | Internals and C implementations |++## References++- [Apache DataSketches](https://datasketches.apache.org/)+- Karnin, Lang, Liberty. [Optimal Quantile Approximation in Streams](https://arxiv.org/abs/1603.05346). FOCS 2016.+- Cormode, Muthukrishnan. [An Improved Data Stream Summary: The Count-Min Sketch and its Applications](https://dimacs.rutgers.edu/~graham/pubs/papers/cm-full.pdf). Journal of Algorithms, 2005.+- Flajolet, Fusy, Gandouet, Meunier. [HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm](http://algo.inria.fr/flajolet/Publications/FlFuGaMe07.pdf). AofA 2007.
bench/Bench.hs view
@@ -1,45 +1,139 @@ {-# OPTIONS_GHC -fno-full-laziness #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE TypeApplications #-} -import Control.Monad (forM_)+import Control.Monad (void) import Criterion.Main import Criterion.Types-import Data.List (foldl')-import DataSketches.Quantiles.RelativeErrorQuantile--- import Prometheus (register, summary, defaultQuantiles, observe, Info (Info))-import Control.Concurrent (withMVar, newMVar)+import Data.Word+import Control.Monad.Primitive (PrimState)+import qualified Data.Vector.Storable as VS+import qualified DataSketches.Quantiles.RelativeErrorQuantile as REQ+import qualified DataSketches.Quantiles.KLL as KLL+import qualified DataSketches.Distinct.HyperLogLog as HLL+import qualified DataSketches.Frequencies.CountMin as CM +loopDouble :: (Double -> IO ()) -> Double -> Double -> IO ()+loopDouble f !lo !hi = go lo+ where+ go !i+ | i > hi = pure ()+ | otherwise = f i >> go (i + 1)+{-# INLINE loopDouble #-}++loopWord64 :: (Word64 -> IO ()) -> Word64 -> Word64 -> IO ()+loopWord64 f !lo !hi = go lo+ where+ go !i+ | i > hi = pure ()+ | otherwise = f i >> go (i + 1)+{-# INLINE loopWord64 #-}++loopDoubleStep :: (Double -> IO ()) -> Double -> Double -> Double -> IO ()+loopDoubleStep f !lo !hi !step = go lo+ where+ go !i+ | i > hi = pure ()+ | otherwise = f i >> go (i + step)+{-# INLINE loopDoubleStep #-}++doubleVec :: Int -> VS.Vector Double+doubleVec n = VS.generate n (fromIntegral . (+ 1))++word64Vec :: Int -> VS.Vector Word64+word64Vec n = VS.generate n (fromIntegral . (+ 1))+ main :: IO ()-main = do- outerSketch <- mkReqSketch 6 HighRanksAreAccurate- -- let metric = summary (Info "adversarial_input" "woo") defaultQuantiles- -- prometheusThing <- register metric- skM <- newMVar =<< mkReqSketch 6 HighRanksAreAccurate- -- mapM_ (update outerSketch) [1..10000]- defaultMain- [ bgroup "ReqSketch"- [ bench "insert/1" $ perRunEnv (mkReqSketch 6 HighRanksAreAccurate) $ \sk -> do- insert sk 1- , bench "insert/10" $ perRunEnv (mkReqSketch 6 HighRanksAreAccurate) $ \sk -> do- mapM_ (insert sk) [1..10]- , bench "insert/100" $ perRunEnv (mkReqSketch 6 HighRanksAreAccurate) $ \sk -> do- mapM_ (insert sk) [1..100]- , bench "insert/1000" $ perRunEnv (mkReqSketch 6 HighRanksAreAccurate) $ \sk -> do- mapM_ (insert sk) [1..1000]- , bench "insert/10000" $ perRunEnv (mkReqSketch 6 HighRanksAreAccurate) $ \sk -> do- mapM_ (insert sk) [1..10000]- , bench "insert/existing" $ whnfIO $ insert outerSketch 1- , bench "insert/mvar" $ whnfIO $ withMVar skM (`insert` 1)+main = defaultMain+ [ bgroup "REQ"+ [ bgroup "insert"+ [ bench "100" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ loopDouble (REQ.insert sk) 1 100+ , bench "1000" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ loopDouble (REQ.insert sk) 1 1000+ , bench "10000" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ loopDouble (REQ.insert sk) 1 10000+ , bench "100000" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ loopDouble (REQ.insert sk) 1 100000 ]- , bgroup "DoubleBuffer"- [ -- bench "sort" $ + , bgroup "insertBatch"+ [ bench "100" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ REQ.insertBatch sk (doubleVec 100)+ , bench "1000" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ REQ.insertBatch sk (doubleVec 1000)+ , bench "10000" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ REQ.insertBatch sk (doubleVec 10000)+ , bench "100000" $ perRunEnv (REQ.mkReqSketch 6 REQ.HighRanksAreAccurate) $ \sk ->+ REQ.insertBatch sk (doubleVec 100000) ]- -- , bgroup "Prometheus"- -- [ bench "insert/existing" $ whnfIO $- -- observe prometheusThing 1- -- ]+ , bgroup "rank"+ [ bench "100-queries-from-10000" $ perRunEnv (mkReqWith 10000) $ \sk ->+ loopDoubleStep (void . REQ.rank sk) 0 9900 100+ ] ]+ , bgroup "KLL"+ [ bgroup "insert"+ [ bench "100" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ loopDouble (KLL.insert sk) 1 100+ , bench "1000" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ loopDouble (KLL.insert sk) 1 1000+ , bench "10000" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ loopDouble (KLL.insert sk) 1 10000+ , bench "100000" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ loopDouble (KLL.insert sk) 1 100000+ ]+ , bgroup "insertBatch"+ [ bench "100" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ KLL.insertBatch sk (doubleVec 100)+ , bench "1000" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ KLL.insertBatch sk (doubleVec 1000)+ , bench "10000" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ KLL.insertBatch sk (doubleVec 10000)+ , bench "100000" $ perRunEnv (KLL.mkKllSketch 200) $ \sk ->+ KLL.insertBatch sk (doubleVec 100000)+ ]+ , bgroup "rank"+ [ bench "100-queries-from-10000" $ perRunEnv (mkKllWith 10000) $ \sk ->+ loopDoubleStep (void . KLL.rank sk) 0 9900 100+ ]+ ]+ , bgroup "HLL"+ [ bgroup "insert"+ [ bench "1000" $ perRunEnv (HLL.mkHllSketch 12) $ \sk ->+ loopWord64 (HLL.insert sk) 1 1000+ , bench "10000" $ perRunEnv (HLL.mkHllSketch 12) $ \sk ->+ loopWord64 (HLL.insert sk) 1 10000+ , bench "100000" $ perRunEnv (HLL.mkHllSketch 12) $ \sk ->+ loopWord64 (HLL.insert sk) 1 100000+ ]+ , bgroup "insertBatch"+ [ bench "1000" $ perRunEnv (HLL.mkHllSketch 12) $ \sk ->+ HLL.insertBatch sk (word64Vec 1000)+ , bench "10000" $ perRunEnv (HLL.mkHllSketch 12) $ \sk ->+ HLL.insertBatch sk (word64Vec 10000)+ , bench "100000" $ perRunEnv (HLL.mkHllSketch 12) $ \sk ->+ HLL.insertBatch sk (word64Vec 100000)+ ]+ ]+ , bgroup "CountMin"+ [ bgroup "insert"+ [ bench "10000" $ perRunEnv (CM.mkCountMinSketch 0.001 0.01) $ \sk ->+ loopWord64 (CM.insert sk) 1 10000+ , bench "100000" $ perRunEnv (CM.mkCountMinSketch 0.001 0.01) $ \sk ->+ loopWord64 (CM.insert sk) 1 100000+ ]+ ]+ ]++mkReqWith :: Int -> IO (REQ.ReqSketch (PrimState IO))+mkReqWith n = do+ sk <- REQ.mkReqSketch 6 REQ.HighRanksAreAccurate+ loopDouble (REQ.insert sk) 1 (fromIntegral n)+ pure sk++mkKllWith :: Int -> IO (KLL.KllSketch (PrimState IO))+mkKllWith n = do+ sk <- KLL.mkKllSketch 200+ loopDouble (KLL.insert sk) 1 (fromIntegral n)+ pure sk
data-sketches.cabal view
@@ -1,17 +1,19 @@ cabal-version: 1.18 --- This file has been generated from package.yaml by hpack version 0.34.4.+-- This file has been generated from package.yaml by hpack version 0.38.1. -- -- see: https://github.com/sol/hpack name: data-sketches-version: 0.3.1.0+version: 0.4.0.0+synopsis: Stochastic streaming algorithms for approximate computation on large datasets. Includes KLL, HLL, Theta, Count-Min, and REQ sketches. description: Please see the README on GitHub at <https://github.com/iand675/datasketches-haskell#readme>+category: Data homepage: https://github.com/iand675/datasketches-haskell#readme bug-reports: https://github.com/iand675/datasketches-haskell/issues author: Ian Duncan, Rob Bassi maintainer: ian@iankduncan.com-copyright: 2021 Ian Duncan, Rob Bassi, Mercury Technologies+copyright: 2025 Ian Duncan, Rob Bassi, Mercury Technologies license: Apache license-file: LICENSE build-type: Simple@@ -30,6 +32,10 @@ library exposed-modules: DataSketches.Quantiles.RelativeErrorQuantile+ DataSketches.Quantiles.KLL+ DataSketches.Frequencies.CountMin+ DataSketches.Distinct.HyperLogLog+ DataSketches.Distinct.Theta other-modules: Paths_data_sketches hs-source-dirs:@@ -44,7 +50,7 @@ TypeOperators build-depends: base >=4.7 && <5- , data-sketches-core ==0.1.*+ , data-sketches-core ==0.2.* , ghc-prim , mtl , mwc-random@@ -58,10 +64,16 @@ main-is: Spec.hs other-modules: AuxiliarySpec+ BugFixSpec CompactorSpec+ CountMinSpec+ CrossValidationSpec DoubleBufferSpec+ HyperLogLogSpec+ KllSpec ProofCheckSpec RelativeErrorQuantileSpec+ ThetaSpec Paths_data_sketches hs-source-dirs: test@@ -78,15 +90,19 @@ QuickCheck , base >=4.7 && <5 , data-sketches- , data-sketches-core ==0.1.*+ , data-sketches-core+ , directory , ghc-prim+ , hedgehog , hspec- , hspec-discover+ , hspec-junit-formatter , mtl , mwc-random , pretty-show , primitive+ , process , statistics+ , temporary , vector , vector-algorithms default-language: Haskell2010@@ -106,12 +122,12 @@ StandaloneDeriving TypeFamilies TypeOperators- ghc-options: -threaded -rtsopts -with-rtsopts=-N+ ghc-options: -threaded -rtsopts -with-rtsopts=-N -O2 build-depends: base >=4.7 && <5 , criterion , data-sketches- , data-sketches-core ==0.1.*+ , data-sketches-core , ghc-prim , mtl , mwc-random
+ src/DataSketches/Distinct/HyperLogLog.hs view
@@ -0,0 +1,99 @@+-- | HyperLogLog sketch for cardinality (distinct count) estimation.+--+-- "How many /different/ things have I seen?" — not how many total, but how+-- many unique. Think unique visitors, unique IPs, unique search queries.+-- Counting exactly requires remembering every item you've seen (a set), which+-- grows with the data. HyperLogLog answers the same question using a few+-- kilobytes, no matter how large the stream.+--+-- It works by hashing each item and tracking the maximum number of leading+-- zeros observed across 2^p independent register buckets.+--+-- The standard error is approximately @1.04 / sqrt(2^p)@.+--+-- === Precision configurations+--+-- +------+-----------+------+-------------++-- | @p@ | Registers | RAM | Std. error |+-- +======+===========+======+=============++-- | 10 | 1,024 | 1 KB | ~3.25% |+-- +------+-----------+------+-------------++-- | 12 | 4,096 | 4 KB | ~1.63% |+-- +------+-----------+------+-------------++-- | 14 | 16,384 | 16 KB| ~0.81% |+-- +------+-----------+------+-------------++-- | 16 | 65,536 | 64 KB| ~0.41% |+-- +------+-----------+------+-------------++--+-- === Hashing+--+-- Items must be pre-hashed to 'Word64'. Apply a good hash function (e.g.+-- from @hashable@ or @xxhash@) to your domain types before calling 'insert'.+-- The quality of the cardinality estimate depends on the hash being uniform.+--+-- === Implementation+--+-- Backed by a C implementation (@cbits\/hll.c@) behind a 'ForeignPtr'.+-- All sketch memory lives outside the GHC heap.+--+-- === Usage+--+-- @+-- import qualified DataSketches.Distinct.HyperLogLog as HLL+-- import Data.Hashable (hash)+--+-- main :: IO ()+-- main = do+-- sk <- HLL.'mkHllSketch' 12+-- mapM_ (HLL.'insert' sk . fromIntegral . hash) [\"alice\", \"bob\", \"alice\"]+-- n <- HLL.'estimate' sk+-- putStrLn $ "distinct count ≈ " ++ show n -- ≈ 2.0+-- @+--+-- === Mergeability+--+-- Fully mergeable via 'merge'. The union of two HLL sketches gives the same+-- result as inserting all items from both streams into a single sketch. Both+-- must share the same precision @p@.+--+-- === When to use Theta instead+--+-- HyperLogLog only supports cardinality estimation and union. If you need+-- set intersection or difference, use "DataSketches.Distinct.Theta".+module DataSketches.Distinct.HyperLogLog+ ( -- * Construction+ HllSketch+ , mkHllSketch+ -- * Updating+ , insert+ , insertBatch+ , merge+ -- * Querying+ , estimate+ , precision+ ) where++import Control.Monad.Primitive (PrimMonad, PrimState)+import qualified Data.Vector.Storable as VS+import Data.Word (Word64)+import DataSketches.Distinct.HyperLogLog.Internal++-- | Insert an item into the sketch. The item should be a hash of the original value.+insert :: PrimMonad m => HllSketch (PrimState m) -> Word64 -> m ()+insert = hllInsert++-- | Bulk-insert a storable vector of hashed items. Avoids per-element FFI overhead.+insertBatch :: PrimMonad m => HllSketch (PrimState m) -> VS.Vector Word64 -> m ()+insertBatch = hllInsertBatch++-- | Estimate the number of distinct items inserted.+estimate :: PrimMonad m => HllSketch (PrimState m) -> m Double+estimate = hllEstimate++-- | Merge the second sketch into the first. Both must have the same precision.+merge :: PrimMonad m => HllSketch (PrimState m) -> HllSketch (PrimState m) -> m ()+merge = hllMerge++-- | Get the precision (log2 of register count) of the sketch.+precision :: PrimMonad m => HllSketch (PrimState m) -> m Int+precision = hllPrecision
+ src/DataSketches/Distinct/Theta.hs view
@@ -0,0 +1,99 @@+-- | Theta sketch for distinct counting with full set algebra.+--+-- Like HyperLogLog, Theta counts distinct items — but it can also do set+-- math. "How many users visited /both/ page A and page B?" (intersection).+-- "How many visited A but /not/ B?" (difference). "How many visited+-- /either/?" (union). You can even chain these into complex expressions+-- like @(A ∪ B) ∩ (C ∪ D) \\ E@, and the result is itself a sketch+-- that you can query or combine further.+--+-- === Algorithm+--+-- Each item is hashed to a uniform value in @[0, 2^64)@. The sketch retains+-- only hashes below a threshold /theta/. When the number of retained entries+-- exceeds @k@, theta is lowered to keep approximately @k@ entries. The+-- cardinality estimate is @(retained entries) / (theta / 2^64)@.+--+-- Set operations combine hash sets with appropriate theta adjustments — the+-- result is itself a sketch that can participate in further set operations.+--+-- Standard error is approximately @1 / sqrt(k)@.+--+-- === Hashing+--+-- Items must be pre-hashed to 'Word64'. Apply a good hash function to your+-- domain types before calling 'insert'.+--+-- === Implementation+--+-- Backed by a C implementation (@cbits\/theta.c@) behind a 'ForeignPtr'.+-- All sketch memory lives outside the GHC heap.+--+-- === Usage+--+-- @+-- import qualified DataSketches.Distinct.Theta as Theta+--+-- main :: IO ()+-- main = do+-- a <- Theta.'mkThetaSketch' 4096+-- b <- Theta.'mkThetaSketch' 4096+-- mapM_ (Theta.'insert' a) [1..1000]+-- mapM_ (Theta.'insert' b) [500..1500]+--+-- both <- Theta.'intersection' a b+-- overlap <- Theta.'estimate' both+-- putStrLn $ "overlap ≈ " ++ show overlap -- ≈ 501+--+-- onlyA <- Theta.'difference' a b+-- exclusive <- Theta.'estimate' onlyA+-- putStrLn $ "only in A ≈ " ++ show exclusive -- ≈ 499+-- @+--+-- === When to use HyperLogLog instead+--+-- If you only need cardinality estimation (no set operations),+-- "DataSketches.Distinct.HyperLogLog" uses less memory for the same accuracy.+module DataSketches.Distinct.Theta+ ( -- * Construction+ ThetaSketch+ , mkThetaSketch+ -- * Updating+ , insert+ -- * Querying+ , estimate+ , isEmpty+ -- * Set operations+ , union+ , intersection+ , difference+ ) where++import Control.Monad.Primitive (PrimMonad(PrimState))+import Data.Word (Word64)+import DataSketches.Distinct.Theta.Internal++-- | Insert an item into the sketch. The item should be a hash of the original value.+insert :: PrimMonad m => ThetaSketch (PrimState m) -> Word64 -> m ()+insert = thetaInsert++-- | Estimate the number of distinct items inserted.+estimate :: PrimMonad m => ThetaSketch (PrimState m) -> m Double+estimate = thetaEstimate++-- | True if no items have been inserted.+isEmpty :: PrimMonad m => ThetaSketch (PrimState m) -> m Bool+isEmpty = thetaIsEmpty++-- | Compute the union of two sketches. Returns a new sketch.+union :: PrimMonad m => ThetaSketch (PrimState m) -> ThetaSketch (PrimState m) -> m (ThetaSketch (PrimState m))+union = thetaUnion++-- | Compute the intersection of two sketches. Returns a new sketch.+intersection :: PrimMonad m => ThetaSketch (PrimState m) -> ThetaSketch (PrimState m) -> m (ThetaSketch (PrimState m))+intersection = thetaIntersection++-- | Compute the set difference (A not B). Returns a new sketch containing+-- items that are in the first sketch but not the second.+difference :: PrimMonad m => ThetaSketch (PrimState m) -> ThetaSketch (PrimState m) -> m (ThetaSketch (PrimState m))+difference = thetaDifference
+ src/DataSketches/Frequencies/CountMin.hs view
@@ -0,0 +1,91 @@+-- | Count-Min sketch for approximate frequency estimation.+--+-- "How many times has /this particular item/ appeared?" — imagine counting+-- ad impressions per user, or requests per API key, when there are too many+-- distinct items to keep an exact counter for each one. Count-Min gives you+-- an approximate count using fixed memory.+--+-- It is conservative: estimates may overcount (due to hash collisions)+-- but /never/ undercount. If an item truly appeared 100 times, the+-- sketch will report >= 100. This makes it safe for rate limiting and+-- frequency capping — you might act slightly early, but never too late.+--+-- === Parameters+--+-- * /epsilon/ (ε): error tolerance. The estimate is within @ε * N@ of the+-- true count with high probability, where @N@ is the total insertions.+-- * /delta/ (δ): failure probability. The probability the estimate exceeds+-- the error bound is at most δ.+--+-- Space used is @O(1\/ε * log(1\/δ))@.+--+-- === Hashing+--+-- Items must be pre-hashed to 'Word64'. Apply a hash function to your+-- domain types before calling 'insert'.+--+-- === Implementation+--+-- Backed by a C implementation (@cbits\/countmin.c@) behind a 'ForeignPtr'.+-- All sketch memory lives outside the GHC heap.+--+-- === Usage+--+-- @+-- import qualified DataSketches.Frequencies.CountMin as CM+-- import Data.Hashable (hash)+--+-- main :: IO ()+-- main = do+-- sk <- CM.'mkCountMinSketch' 0.001 0.01 -- ε=0.1%, δ=1%+-- mapM_ (\\_ -> CM.'insert' sk 42) [1..500]+-- freq <- CM.'estimate' sk 42+-- putStrLn $ "frequency ≈ " ++ show freq -- ≈ 500+-- @+--+-- === Mergeability+--+-- Fully mergeable via 'merge'. Both sketches must have been constructed+-- with the same ε and δ (i.e. identical dimensions).+module DataSketches.Frequencies.CountMin+ ( -- * Construction+ CountMinSketch+ , mkCountMinSketch+ -- * Updating+ , insert+ , insertN+ , merge+ -- * Querying+ , estimate+ , width+ , depth+ ) where++import Control.Monad.Primitive (PrimMonad, PrimState)+import Data.Word (Word64)+import DataSketches.Frequencies.CountMin.Internal++-- | Insert a single occurrence of an item.+insert :: PrimMonad m => CountMinSketch (PrimState m) -> Word64 -> m ()+insert = cmsInsert++-- | Insert multiple occurrences of an item.+insertN :: PrimMonad m => CountMinSketch (PrimState m) -> Word64 -> Word64 -> m ()+insertN = cmsInsertN++-- | Estimate the frequency of an item. May overcount but never undercounts.+estimate :: PrimMonad m => CountMinSketch (PrimState m) -> Word64 -> m Word64+estimate = cmsEstimate++-- | Merge the second sketch into the first. Both must have the same dimensions+-- (same epsilon and delta parameters at construction).+merge :: PrimMonad m => CountMinSketch (PrimState m) -> CountMinSketch (PrimState m) -> m ()+merge = cmsMerge++-- | Number of columns in the sketch.+width :: PrimMonad m => CountMinSketch (PrimState m) -> m Int+width = cmsWidth++-- | Number of rows (hash functions) in the sketch.+depth :: PrimMonad m => CountMinSketch (PrimState m) -> m Int+depth = cmsDepth
+ src/DataSketches/Quantiles/KLL.hs view
@@ -0,0 +1,154 @@+-- | KLL (Karnin-Lang-Liberty) sketch for approximate quantile computation+-- with additive error bounds. Based on "Optimal Quantile Approximation+-- in Streams" (Karnin, Lang, Liberty — FOCS 2016).+--+-- Feed it a stream of numbers, then ask questions like "what's the median?"+-- or "what value is at the 99th percentile?" — without storing every value.+-- The sketch keeps a small, fixed-size buffer regardless of how many items+-- you insert.+--+-- KLL provides uniform additive error across all ranks. With @k = 200@,+-- the normalized rank error is approximately 1.33% with 99% confidence,+-- regardless of whether you query p50 or p99.+--+-- This is the recommended quantile sketch for most use cases. Prefer+-- "DataSketches.Quantiles.RelativeErrorQuantile" only when you need+-- high accuracy specifically at the distribution tails (p99.9+).+--+-- === Implementation+--+-- Backed by a C implementation (@cbits\/kll.c@) behind a 'ForeignPtr'.+-- All sketch memory lives outside the GHC heap — zero GC pressure,+-- zero boxing overhead.+--+-- === Usage+--+-- @+-- import qualified DataSketches.Quantiles.KLL as KLL+--+-- main :: IO ()+-- main = do+-- sk <- KLL.'mkKllSketch' 200+-- mapM_ (KLL.'insert' sk) [1..100000 :: Double]+-- p99 <- KLL.'quantile' sk 0.99+-- putStrLn $ "p99 = " ++ show p99+-- @+--+-- === Mergeability+--+-- Fully mergeable: 'merge' combines two sketches as if all data had been+-- inserted into one. Suitable for parallel and distributed computation.+module DataSketches.Quantiles.KLL+ ( -- * Construction+ KllSketch+ , mkKllSketch+ -- * Updating the sketch+ , insert+ , insertBatch+ , merge+ -- * Querying+ , count+ , null+ , minimum+ , maximum+ , retainedItemCount+ , quantile+ , quantiles+ , rank+ , ranks+ , cumulativeDistributionFunction+ , probabilityMassFunction+ ) where++import Prelude hiding (null, minimum, maximum)+import Control.Monad (unless, when)+import Control.Monad.Primitive (PrimMonad(PrimState))+import qualified Data.Vector.Storable as VS+import Data.Word (Word32, Word64)+import DataSketches.Quantiles.KLL.Internal++-- | Create a new KLL sketch with the given accuracy parameter k.+-- Larger k gives better accuracy but uses more space.+-- k must be >= 8. Default recommendation is 200.++-- | Insert a value into the sketch. NaN values are ignored.+insert :: PrimMonad m => KllSketch (PrimState m) -> Double -> m ()+insert = kllInsert++-- | Bulk-insert a storable vector of doubles. Avoids per-element FFI overhead.+insertBatch :: PrimMonad m => KllSketch (PrimState m) -> VS.Vector Double -> m ()+insertBatch = kllInsertBatch++-- | Merge the second sketch into the first.+merge :: PrimMonad m => KllSketch (PrimState m) -> KllSketch (PrimState m) -> m ()+merge = kllMerge++-- | Total number of items inserted.+count :: PrimMonad m => KllSketch (PrimState m) -> m Word64+count = kllCount++-- | True if no items have been inserted.+null :: PrimMonad m => KllSketch (PrimState m) -> m Bool+null = kllIsEmpty++-- | Smallest value seen.+minimum :: PrimMonad m => KllSketch (PrimState m) -> m Double+minimum = kllMinimum++-- | Largest value seen.+maximum :: PrimMonad m => KllSketch (PrimState m) -> m Double+maximum = kllMaximum++-- | Number of items currently retained in the sketch.+retainedItemCount :: PrimMonad m => KllSketch (PrimState m) -> m Int+retainedItemCount = kllRetainedItems++-- | Get the approximate quantile for the given normalized rank in [0, 1].+quantile :: PrimMonad m => KllSketch (PrimState m) -> Double -> m Double+quantile = kllQuantile++-- | Get approximate quantiles for multiple normalized ranks.+quantiles :: PrimMonad m => KllSketch (PrimState m) -> [Double] -> m [Double]+quantiles sk = mapM (quantile sk)++-- | Get the approximate normalized rank of a value.+rank :: PrimMonad m => KllSketch (PrimState m) -> Double -> m Double+rank = kllRank++-- | Get approximate normalized ranks for multiple values.+ranks :: PrimMonad m => KllSketch (PrimState m) -> [Double] -> m [Double]+ranks sk = mapM (rank sk)++-- | Approximate cumulative distribution function.+-- Given a list of split points, returns a list of cumulative probabilities.+-- Returns Nothing if the sketch is empty.+cumulativeDistributionFunction+ :: PrimMonad m+ => KllSketch (PrimState m)+ -> [Double]+ -> m (Maybe [Double])+cumulativeDistributionFunction sk splitPoints = do+ empty <- null sk+ if empty+ then pure Nothing+ else do+ rs <- mapM (rank sk) splitPoints+ pure (Just (rs ++ [1.0]))++-- | Approximate probability mass function.+-- Given a list of split points, returns probabilities for each interval.+-- Returns empty list if the sketch is empty.+probabilityMassFunction+ :: PrimMonad m+ => KllSketch (PrimState m)+ -> [Double]+ -> m [Double]+probabilityMassFunction sk splitPoints = do+ empty <- null sk+ if empty+ then pure []+ else do+ rs <- mapM (rank sk) splitPoints+ let cdf = rs ++ [1.0]+ pmf = zipWith (-) cdf (0.0 : cdf)+ pure pmf
src/DataSketches/Quantiles/RelativeErrorQuantile.hs view
@@ -1,37 +1,65 @@--- | The Relative Error Quantile (REQ) sketch provides extremely high accuracy at a chosen end of the rank domain. --- This is best illustrated with some rank domain accuracy plots that compare the KLL quantiles sketch to the REQ sketch.+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE TypeApplications #-}+-- | Relative Error Quantiles (REQ) sketch for approximate quantile computation+-- with /relative/ error bounds. ----- This first plot illustrates the typical error behavior of the KLL sketch (also the quantiles/DoublesSketch). --- The error is flat for all ranks (0, 1). The green and yellow lines correspond to +/- one RSE at 68% confidence; --- the blue and red lines, +/- two RSE at 95% confidence; and, the purple and brown lines +/- 3 RSE at 99% confidence. --- The reason all the curves pinch at 0 and 1.0, is because the sketch knows with certainty that a request for a quantile at --- rank = 0 is the minimum value of the stream; and a request for a quantiles at rank = 1.0, is the maximum value of the stream. --- Both of which the sketch tracks.+-- Most quantile sketches (like KLL) have /additive/ error: the answer might+-- be off by, say, 1.3% of the full rank range no matter where you query.+-- That's fine for p50 or p90, but at p99.99 a 1.3% additive error is enormous+-- relative to the sliver of the distribution you're looking at. ----- +-- REQ flips this: the error is /proportional to the distance from the nearest+-- tail/. At p99.99, the error is proportional to 0.01%, not 1.3%. This makes+-- REQ the right choice when you care specifically about extreme percentiles —+-- SLA monitoring at p99.9, tail-latency debugging, outlier detection. ----- The next plot is the exact same data and queries fed to the REQ sketch set for High Rank Accuracy (HRA) mode. --- In this plot, starting at a rank of about 0.3, the contour lines start converging and actually reach zero error at --- rank 1.0. Therefore the error (the inverse of accuracy) is relative to the requested rank, thus the name of the sketch. --- This means that the user can perform getQuantile(rank) queries, where rank = .99999 and get accurate results.+-- For general-purpose quantiles where you don't need tail precision, prefer+-- "DataSketches.Quantiles.KLL" — it's faster and simpler. ----- +-- Based on "Relative Error Streaming Quantiles" (Cormode et al., PODS 2021). ----- This next plot is also the same data and queries, except the REQ sketch was configured for Low Rank Accuracy (LRA). In this case the user can perform getQuantiles(rank) queries, where rank = .00001 and get accurate results.+-- === Configuration ----- --{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE TypeApplications #-}-{-# LANGUAGE DeriveGeneric #-}+-- @k@ controls accuracy vs. space. Must be even, in @[4, 1024]@. Larger @k@+-- means tighter error bounds but more memory. @k = 12@ is a reasonable default.+--+-- 'RankAccuracy' selects which tail gets the better error guarantee:+--+-- * 'HighRanksAreAccurate' — use when you care about p99, p99.9, p99.99+-- * 'LowRanksAreAccurate' — use when you care about p1, p0.1, p0.01+--+-- === Usage+--+-- @+-- import qualified DataSketches.Quantiles.RelativeErrorQuantile as REQ+--+-- main :: IO ()+-- main = do+-- sk <- REQ.'mkReqSketch' 12 REQ.'HighRanksAreAccurate'+-- mapM_ (REQ.'insert' sk) [1..100000 :: Double]+-- p999 <- REQ.'quantile' sk 0.999+-- putStrLn $ "p99.9 = " ++ show p999+-- @+--+-- === Implementation+--+-- Backed by a C implementation (@cbits\/req.c@) behind a 'ForeignPtr'.+-- All sketch memory lives outside the GHC heap — zero GC pressure,+-- zero boxing overhead.+--+-- === Mergeability+--+-- REQ sketches are fully mergeable: two sketches built on separate data+-- partitions can be combined via 'merge' to produce a sketch equivalent+-- to having processed all data in a single stream. Both sketches must+-- share the same 'RankAccuracy' setting. module DataSketches.Quantiles.RelativeErrorQuantile ( -- * Construction- ReqSketch (criterion)+ ReqSketch , mkReqSketch -- ** Configuration settings , RankAccuracy(..)- , Criterion(..)- -- * Sketch summaries+ -- * Querying , count , null , sum@@ -49,484 +77,203 @@ , rankUpperBound , cumulativeDistributionFunction , getK+ , numLevels -- * Updating the sketch , merge , insert+ , insertBatch+ -- * Inspection , rankAccuracy , isEstimationMode , isLessThanOrEqual- -- | If you see this error, please file an issue in the GitHub repository.- , CumulativeDistributionInvariants(..)+ , setCriterionLE+ , setCriterionLT+ , DoubleIsNonFiniteException(..) ) where -import Control.Monad (when, unless, foldM, foldM_)+import Control.Exception (throw)+import Control.Monad (when) import Control.Monad.Primitive ( PrimMonad(PrimState) )-import Data.Bits (shiftL)-import Data.Vector ((!), imapM_)-import qualified Data.Vector as Vector-import Data.Primitive.MutVar- ( modifyMutVar', newMutVar, readMutVar, writeMutVar ) import Data.Word ( Word32, Word64 )+import DataSketches.Quantiles.RelativeErrorQuantile.Types+ ( RankAccuracy(..) ) import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Constants ( fixRseFactor, initNumberOfSections, relRseFactor )-import DataSketches.Quantiles.RelativeErrorQuantile.Types- ( Criterion(..), RankAccuracy(..) )-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor (ReqCompactor)-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary (ReqAuxiliary)-import DataSketches.Quantiles.RelativeErrorQuantile.Internal- ( count,- retainedItemCount,- CumulativeDistributionInvariants(..),- ReqSketch(..),- computeTotalRetainedItems,- getCompactors )-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary as Auxiliary-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor as Compactor-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer as DoubleBuffer-import DataSketches.Core.Internal.URef- ( modifyURef, newURef, readURef, writeURef )-import Data.Maybe (isNothing)-import qualified Data.Foldable+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer+ ( DoubleIsNonFiniteException(..) )+import DataSketches.Quantiles.RelativeErrorQuantile.CInternal import qualified Data.List-import GHC.Exception.Type (Exception)-import Control.Exception (throw, assert)-import System.Random.MWC (Gen, create)-import qualified Data.Vector.Generic.Mutable as MG+import qualified Data.Vector.Storable as VS import Prelude hiding (sum, minimum, maximum, null) --- | The K parameter can be increased to trade increased space efficiency for higher accuracy in rank and quantile--- calculations. Due to the way the compaction algorithm works, it must be an even number between 4 and 1024.------ A good starting number when in doubt is 6.-mkReqSketch :: forall m. (PrimMonad m)- => Word32 -- ^ K+type ReqSketch = CReqSketch++mkReqSketch :: PrimMonad m+ => Word32 -> RankAccuracy -> m (ReqSketch (PrimState m))-mkReqSketch k rank = do- unless (even k && k >= 4 && k <= 1024) $ error "k must be divisible by 2, and satisfy 4 <= k <= 1024"- r <- ReqSketch k rank (:<)- <$> create- <*> newURef 0- <*> newURef (0 / 0)- <*> newURef (0 / 0)- <*> newURef 0- <*> newURef 0- <*> newURef 0- <*> newMutVar Nothing- <*> newMutVar Vector.empty- grow r- pure r---getAux :: PrimMonad m => ReqSketch (PrimState m) -> m (Maybe ReqAuxiliary)-getAux = readMutVar . aux--getNumLevels :: PrimMonad m => ReqSketch (PrimState m) -> m Int-getNumLevels = fmap Vector.length . getCompactors--getIsEmpty :: PrimMonad m => ReqSketch (PrimState m) -> m Bool-getIsEmpty = fmap (== 0) . readURef . totalN--getK :: ReqSketch s -> Word32-getK = k--getMaxNominalCapacity :: PrimMonad m => ReqSketch (PrimState m) -> m Int-getMaxNominalCapacity = readURef . maxNominalCapacitiesSize--validateSplits :: Monad m => [Double] -> m ()-validateSplits splits = do- when (Data.Foldable.null splits) $ do- throw CumulativeDistributionInvariantsSplitsAreEmpty- when (any isInfinite splits || any isNaN splits) $ do- throw CumulativeDistributionInvariantsSplitsAreNotFinite- when (Data.List.nub (Data.List.sort splits) /= splits) $ do- throw CumulativeDistributionInvariantsSplitsAreNotUniqueAndMontonicallyIncreasing--getCounts :: (PrimMonad m) => ReqSketch (PrimState m) -> [Double] -> m [Word64]-getCounts this values = do- compactors <- getCompactors this- let numValues = length values- numCompactors = Vector.length compactors- ans = replicate numValues 0- isEmpty <- getIsEmpty this- if isEmpty- then pure []- else Vector.ifoldM doCount ans compactors+mkReqSketch k ra = mkCReqSketch k (raToInt ra) where- doCount acc index compactor = do- let wt = (1 `shiftL` fromIntegral (Compactor.getLgWeight compactor)) :: Word64- buff <- Compactor.getBuffer compactor- let updateCounts buff value = do- count_ <- DoubleBuffer.getCountWithCriterion buff (values !! index) (criterion this)- pure $ fromIntegral value + fromIntegral count_ * wt- mapM (updateCounts buff) acc+ raToInt HighRanksAreAccurate = 1+ raToInt LowRanksAreAccurate = 0 -getPMForCDF :: (PrimMonad m) => ReqSketch (PrimState m) -> [Double] -> m [Word64]-getPMForCDF this splits = do- () <- validateSplits splits- let numSplits = length splits- numBuckets = numSplits -- + 1- splitCounts <- getCounts this splits- n <- count this- pure $ (++ [n]) $ take numBuckets splitCounts+insert :: PrimMonad m => ReqSketch (PrimState m) -> Double -> m ()+insert = creqInsert+{-# INLINE insert #-} --- | Returns an approximation to the Cumulative Distribution Function (CDF), which is the cumulative analog of the PMF, --- of the input stream given a set of splitPoint (values).-cumulativeDistributionFunction- :: (PrimMonad m)- => ReqSketch (PrimState m)- -> [Double]- -- ^ Returns an approximation to the Cumulative Distribution Function (CDF), - -- which is the cumulative analog of the PMF, of the input stream given a set of - -- splitPoint (values).- --- -- The resulting approximations have a probabilistic guarantee that be obtained, - -- a priori, from the getRSE(int, double, boolean, long) function.- --- -- If the sketch is empty this returns 'Nothing'.- -> m (Maybe [Double])-cumulativeDistributionFunction this splitPoints = do- buckets <- getPMForCDF this splitPoints- isEmpty <- getIsEmpty this- if isEmpty- then pure Nothing- else do- let numBuckets = length splitPoints + 1- n <- count this- pure $ Just $ (/ fromIntegral n) . fromIntegral <$> buckets+-- | Bulk-insert a storable vector of doubles. Avoids per-element FFI overhead.+insertBatch :: PrimMonad m => ReqSketch (PrimState m) -> VS.Vector Double -> m ()+insertBatch = creqInsertBatch -rankAccuracy :: ReqSketch s -> RankAccuracy-rankAccuracy = rankAccuracySetting+merge :: PrimMonad m => ReqSketch (PrimState m) -> ReqSketch (PrimState m) -> m (ReqSketch (PrimState m))+merge dst src = creqMerge dst src >> pure dst --- | Returns an a priori estimate of relative standard error (RSE, expressed as a number in [0,1]). Derived from Lemma 12 in https://arxiv.org/abs/2004.01668v2, but the constant factors were modified based on empirical measurements.-relativeStandardError- :: Int- -- ^ k - the given value of k- -> Double- -- ^ rank - the given normalized rank, a number in [0,1].- -> RankAccuracy- -> Word64- -- ^ totalN - an estimate of the total number of items submitted to the sketch.- -> Double- -- ^ an a priori estimate of relative standard error (RSE, expressed as a number in [0,1]).-relativeStandardError k rank_ hra = getRankUB k 2 rank_ 1 isHra- where- isHra = case hra of- HighRanksAreAccurate -> True- _ -> False+count :: PrimMonad m => ReqSketch (PrimState m) -> m Word64+count = creqCount+{-# INLINE count #-} --- | Gets the smallest value seen by this sketch+null :: PrimMonad m => ReqSketch (PrimState m) -> m Bool+null = creqIsEmpty+ minimum :: PrimMonad m => ReqSketch (PrimState m) -> m Double-minimum = readURef . minValue+minimum = creqMin --- | Gets the largest value seen by this sketch maximum :: PrimMonad m => ReqSketch (PrimState m) -> m Double-maximum = readURef . maxValue+maximum = creqMax --- | Returns the approximate count of items satisfying the criterion set in the ReqSketch 'criterion' field.-countWithCriterion :: (PrimMonad m, s ~ PrimState m) => ReqSketch s -> Double -> m Word64-countWithCriterion s value = fromIntegral <$> do- empty <- null s- if empty- then pure 0- else do- compactors <- getCompactors s- let go !accum compactor = do- let wt = (1 `shiftL` fromIntegral (Compactor.getLgWeight compactor)) :: Word64- buf <- Compactor.getBuffer compactor- count_ <- DoubleBuffer.getCountWithCriterion buf value (criterion s)- pure (accum + (fromIntegral count_ * wt))- Vector.foldM go 0 compactors+sum :: PrimMonad m => ReqSketch (PrimState m) -> m Double+sum = creqSum -sum :: (PrimMonad m) => ReqSketch (PrimState m) -> m Double-sum = readURef . sumValue+retainedItemCount :: PrimMonad m => ReqSketch (PrimState m) -> m Int+retainedItemCount = creqRetained --- | Returns an approximation to the Probability Mass Function (PMF) of the input stream given a set of splitPoints (values).--- The resulting approximations have a probabilistic guarantee that be obtained, a priori, from the getRSE(int, double, boolean, long) function.------ If the sketch is empty this returns an empty list.-probabilityMassFunction- :: (PrimMonad m)- => ReqSketch (PrimState m)- -> [Double]- -- ^ splitPoints - an array of m unique, monotonically increasing double values that divide - -- the real number line into m+1 consecutive disjoint intervals. The definition of an "interval" - -- is inclusive of the left splitPoint (or minimum value) and exclusive of the right splitPoint, - -- with the exception that the last interval will include the maximum value. It is not necessary - -- to include either the min or max values in these splitpoints.- -> m [Double]- -- ^ An array of m+1 doubles each of which is an approximation to the fraction of - -- the input stream values (the mass) that fall into one of those intervals. - -- The definition of an "interval" is inclusive of the left splitPoint and exclusive - -- of the right splitPoint, with the exception that the last interval will - -- include maximum value.-probabilityMassFunction this splitPoints = do- isEmpty <- getIsEmpty this- if isEmpty- then pure []- else do- let numBuckets = length splitPoints + 1- buckets <- fmap fromIntegral <$> getPMForCDF this splitPoints- total <- fromIntegral <$> count this- let computeProb (0, bucket) = bucket / total- computeProb (i, bucket) = (prevBucket + bucket) / total- where prevBucket = buckets !! i - 1- probs = computeProb <$> zip [0..] buckets- pure probs+getK :: PrimMonad m => ReqSketch (PrimState m) -> m Word32+getK = creqK --- | Gets the approximate quantile of the given normalized rank based on the lteq criterion.-quantile- :: (PrimMonad m)- => ReqSketch (PrimState m)- -> Double- -- ^ normRank - the given normalized rank- -> m Double- -- ^ the approximate quantile given the normalized rank.-quantile this normRank = do- isEmpty <- getIsEmpty this- if isEmpty- then pure (0/0)- else do- when (normRank < 0 || normRank > 1.0) $- error $ "Normalized rank must be in the range [0.0, 1.0]: " ++ show normRank- currAuxiliary <- getAux this- when (isNothing currAuxiliary) $ do- total <- count this- retainedItems <- retainedItemCount this- compactors <- getCompactors this- newAuxiliary <- Auxiliary.mkAuxiliary (rankAccuracySetting this) total retainedItems compactors- writeMutVar (aux this) (Just newAuxiliary)- mAuxiliary <- getAux this- case mAuxiliary of- Just auxiliary -> pure $! Auxiliary.getQuantile auxiliary normRank $ criterion this- Nothing -> error "invariant violated: aux is not set"+rankAccuracy :: PrimMonad m => ReqSketch (PrimState m) -> m RqAcc+rankAccuracy sk = do+ i <- creqRankAccuracy sk+ pure $ if i == 1 then HighRanksAreAccurate else LowRanksAreAccurate --- | Gets an array of quantiles that correspond to the given array of normalized ranks.-quantiles- :: (PrimMonad m)- => ReqSketch (PrimState m)- -> [Double]- -- ^ normRanks - the given array of normalized ranks.- -> m [Double]- -- ^ the array of quantiles that correspond to the given array of normalized ranks.-quantiles this normRanks = do- isEmpty <- getIsEmpty this- if isEmpty- then pure []- else mapM (quantile this) normRanks+type RqAcc = RankAccuracy --- | Computes the normalized rank of the given value in the stream. The normalized rank is the fraction of values less than the given value; or if lteq is true, the fraction of values less than or equal to the given value.-rank :: (PrimMonad m)- => ReqSketch (PrimState m)- -> Double- -- ^ value - the given value- -> m Double- -- ^ the normalized rank of the given value in the stream.-rank s value = do- isEmpty <- null s- if isEmpty- then pure (0 / 0) -- NaN- else do- nnCount <- countWithCriterion s value- total <- readURef $ totalN s- pure (fromIntegral nnCount / fromIntegral total)+numLevels :: PrimMonad m => ReqSketch (PrimState m) -> m Int+numLevels = creqNumLevels +isEstimationMode :: PrimMonad m => ReqSketch (PrimState m) -> m Bool+isEstimationMode sk = (> 1) <$> creqNumLevels sk --- getRankLB k levels rank numStdDev hra totalN = if exactRank k levels rank hra totalN+isLessThanOrEqual :: PrimMonad m => ReqSketch (PrimState m) -> m Bool+isLessThanOrEqual sk = (== 1) <$> creqCriterion sk --- | Returns an approximate lower bound rank of the given normalized rank.-rankLowerBound- :: (PrimMonad m)- => ReqSketch (PrimState m)- -> Double- -- ^ rank - the given rank, a value between 0 and 1.0.- -> Int- -- ^ numStdDev - the number of standard deviations. Must be 1, 2, or 3.- -> m Double- -- ^ an approximate lower bound rank.-rankLowerBound this rank numStdDev = do- numLevels <- getNumLevels this- let k = fromIntegral $ getK this- total <- count this- pure $ getRankLB k numLevels rank numStdDev (rankAccuracySetting this == HighRanksAreAccurate) total+setCriterionLE :: PrimMonad m => ReqSketch (PrimState m) -> m ()+setCriterionLE sk = creqSetCriterion sk 1 --- | Gets an array of normalized ranks that correspond to the given array of values.--- TODO, make it ifaster-ranks :: (PrimMonad m, s ~ PrimState m) => ReqSketch s -> [Double] -> m [Double]-ranks s values = mapM (rank s) values+setCriterionLT :: PrimMonad m => ReqSketch (PrimState m) -> m ()+setCriterionLT sk = creqSetCriterion sk 0 --- | Returns an approximate upper bound rank of the given rank.-rankUpperBound- :: (PrimMonad m)- => ReqSketch (PrimState m)- -> Double- -- ^ rank - the given rank, a value between 0 and 1.0.- -> Int- -- ^ numStdDev - the number of standard deviations. Must be 1, 2, or 3.- -> m Double- -- ^ an approximate upper bound rank.-rankUpperBound this rank numStdDev= do- numLevels <- getNumLevels this- let k = fromIntegral $ getK this- total <- count this- pure $ getRankUB k numLevels rank numStdDev (rankAccuracySetting this == HighRanksAreAccurate) total+countWithCriterion :: PrimMonad m => ReqSketch (PrimState m) -> Double -> m Word64+countWithCriterion sk value = do+ when (isNaN value || isInfinite value) $ throw (DoubleIsNonFiniteException value)+ creqCountWithCriterion sk value --- | Returns true if this sketch is empty.-null :: (PrimMonad m) => ReqSketch (PrimState m) -> m Bool-null = fmap (== 0) . readURef . totalN+rank :: PrimMonad m => ReqSketch (PrimState m) -> Double -> m Double+rank sk value = do+ when (isNaN value || isInfinite value) $ throw (DoubleIsNonFiniteException value)+ creqRank sk value --- | Returns true if this sketch is in estimation mode.-isEstimationMode :: PrimMonad m => ReqSketch (PrimState m) -> m Bool-isEstimationMode = fmap (> 1) . getNumLevels+ranks :: PrimMonad m => ReqSketch (PrimState m) -> [Double] -> m [Double]+ranks sk = mapM (rank sk) --- | Returns the current comparison criterion.-isLessThanOrEqual :: ReqSketch s -> Bool-isLessThanOrEqual s = case criterion s of- (:<) -> False- (:<=) -> True+quantile :: PrimMonad m => ReqSketch (PrimState m) -> Double -> m Double+quantile sk normRank = do+ empty <- null sk+ if empty+ then pure (0/0)+ else do+ when (normRank < 0 || normRank > 1.0) $+ error $ "Normalized rank must be in the range [0.0, 1.0]: " ++ show normRank+ creqQuantile sk normRank -computeMaxNominalSize :: PrimMonad m => ReqSketch (PrimState m) -> m Int-computeMaxNominalSize this = do- compactors <- getCompactors this- Vector.foldM countNominalCapacity 0 compactors- where- countNominalCapacity acc compactor = do- nominalCapacity <- Compactor.getNominalCapacity compactor- pure $ nominalCapacity + acc+quantiles :: PrimMonad m => ReqSketch (PrimState m) -> [Double] -> m [Double]+quantiles sk = mapM (quantile sk) -grow :: (PrimMonad m) => ReqSketch (PrimState m) -> m ()-grow this = do- lgWeight <- fromIntegral <$> getNumLevels this- let rankAccuracy = rankAccuracySetting this- sectionSize = getK this- newCompactor <- Compactor.mkReqCompactor (sketchRng this) lgWeight rankAccuracy sectionSize- modifyMutVar' (compactors this) (`Vector.snoc` newCompactor)- maxNominalCapacity <- computeMaxNominalSize this- writeURef (maxNominalCapacitiesSize this) maxNominalCapacity+rankLowerBound :: PrimMonad m => ReqSketch (PrimState m) -> Double -> Int -> m Double+rankLowerBound sk r numStdDev = do+ numLevels <- creqNumLevels sk+ k_ <- fromIntegral <$> getK sk+ total <- count sk+ ra <- rankAccuracy sk+ pure $ getRankLB k_ numLevels r numStdDev (ra == HighRanksAreAccurate) total -compress :: (PrimMonad m) => ReqSketch (PrimState m) -> m ()-compress this = do- compactors <- getCompactors this- let compressionStep height compactor = do- buffSize <- DoubleBuffer.getCount =<< Compactor.getBuffer compactor- nominalCapacity <- Compactor.getNominalCapacity compactor- when (buffSize >= nominalCapacity) $ do- numLevels <- getNumLevels this- when (height + 1 >= numLevels) $ do- grow this- compactors' <- getCompactors this- cReturn <- Compactor.compact compactor- let topCompactor = compactors' ! (height + 1)- buff <- Compactor.getBuffer topCompactor- DoubleBuffer.mergeSortIn buff $ Compactor.crDoubleBuffer cReturn- modifyURef (retainedItems this) (+ Compactor.crDeltaRetItems cReturn)- modifyURef (maxNominalCapacitiesSize this) (+ Compactor.crDeltaNominalSize cReturn)- imapM_ compressionStep compactors- writeMutVar (aux this) Nothing+rankUpperBound :: PrimMonad m => ReqSketch (PrimState m) -> Double -> Int -> m Double+rankUpperBound sk r numStdDev = do+ numLevels <- creqNumLevels sk+ k_ <- fromIntegral <$> getK sk+ total <- count sk+ ra <- rankAccuracy sk+ pure $ getRankUB k_ numLevels r numStdDev (ra == HighRanksAreAccurate) total --- | Merge other sketch into this one.-merge- :: (PrimMonad m, s ~ PrimState m)- => ReqSketch s- -> ReqSketch s- -> m (ReqSketch s)-merge this other = do- otherIsEmpty <- getIsEmpty other- unless otherIsEmpty $ do- let rankAccuracy = rankAccuracySetting this- otherRankAccuracy = rankAccuracySetting other- when (rankAccuracy /= otherRankAccuracy) $- error "Both sketches must have the same HighRankAccuracy setting."- -- update total- otherN <- count other- modifyURef (totalN this) (+ otherN)- -- update the min and max values- thisMin <- minimum this- thisMax <- maximum this- otherMin <- minimum other- otherMax <- maximum other- when (isNaN thisMin || otherMin < thisMin) $ do- writeURef (minValue this) otherMin- when (isNaN thisMax || otherMax < thisMax) $ do- writeURef (maxValue this) otherMax- -- grow until this has at least as many compactors as other- numRequiredCompactors <- getNumLevels other- growUntil numRequiredCompactors- -- merge the items in all height compactors- thisCompactors <- getCompactors this- otherCompactors <- getCompactors other- Vector.zipWithM_ Compactor.merge thisCompactors otherCompactors- -- update state- maxNominalCapacity <- computeMaxNominalSize this- totalRetainedItems <- computeTotalRetainedItems this- writeURef (maxNominalCapacitiesSize this) maxNominalCapacity- writeURef (retainedItems this) totalRetainedItems- -- compress and check invariants- when (totalRetainedItems >= maxNominalCapacity) $ do- compress this- maxNominalCapacity' <- readURef $ maxNominalCapacitiesSize this- totalRetainedItems' <- readURef $ retainedItems this- assert (totalRetainedItems' < maxNominalCapacity') $- writeMutVar (aux this) Nothing- pure this- where- growUntil target = do- numCompactors <- getNumLevels this- when (numCompactors < target) $- grow this+relativeStandardError :: Int -> Double -> RankAccuracy -> Word64 -> Double+relativeStandardError k_ rank_ hra = getRankUB k_ 2 rank_ 1 (hra == HighRanksAreAccurate) --- | Updates this sketch with the given item.-insert :: (PrimMonad m) => ReqSketch (PrimState m) -> Double -> m ()-insert this item = do- unless (isNaN item) $ do- isEmpty <- getIsEmpty this- if isEmpty- then do- writeURef (minValue this) item- writeURef (maxValue this) item- else do- min_ <- minimum this- max_ <- maximum this- when (item < min_) $ writeURef (minValue this) item- when (item > max_) $ writeURef (maxValue this) item- compactor <- Vector.head <$> getCompactors this- buff <- Compactor.getBuffer compactor- DoubleBuffer.append buff item- modifyURef (retainedItems this) (+1)- modifyURef (totalN this) (+1)- modifyURef (sumValue this) (+ item)- retItems <- retainedItemCount this- maxNominalCapacity <- getMaxNominalCapacity this- when (retItems >= maxNominalCapacity) $ do- DoubleBuffer.sort buff- compress this- writeMutVar (aux this) Nothing+cumulativeDistributionFunction+ :: PrimMonad m => ReqSketch (PrimState m) -> [Double] -> m (Maybe [Double])+cumulativeDistributionFunction sk splitPoints = do+ validateSplits splitPoints+ empty <- null sk+ if empty+ then pure Nothing+ else do+ rs <- mapM (rank sk) splitPoints+ pure $ Just (rs ++ [1.0]) --- Private pure bits+probabilityMassFunction :: PrimMonad m => ReqSketch (PrimState m) -> [Double] -> m [Double]+probabilityMassFunction sk splitPoints = do+ empty <- null sk+ if empty+ then pure []+ else do+ rs <- mapM (rank sk) splitPoints+ let cdf = rs ++ [1.0]+ pmf = zipWith (-) cdf (0.0 : cdf)+ pure pmf +-- Private helpers++validateSplits :: Monad m => [Double] -> m ()+validateSplits splits = do+ when (Data.List.null splits) $+ error "splits must not be empty"+ case filter (\x -> isInfinite x || isNaN x) splits of+ (badValue:_) -> throw (DoubleIsNonFiniteException badValue)+ [] -> pure ()+ getRankLB :: Int -> Int -> Double -> Int -> Bool -> Word64 -> Double-getRankLB k levels rank numStdDev hra totalN = if exactRank k levels rank hra totalN- then rank+getRankLB k_ levels r numStdDev hra totalN_ =+ if exactRank k_ levels r hra totalN_ then r else max lbRel lbFix where- relative = relRseFactor / fromIntegral k * (if hra then 1.0 - rank else rank)- fixed = fixRseFactor / fromIntegral k- lbRel = rank - fromIntegral numStdDev * relative- lbFix = rank - fromIntegral numStdDev * fixed+ relative = relRseFactor / fromIntegral k_ * (if hra then 1.0 - r else r)+ fixed = fixRseFactor / fromIntegral k_+ lbRel = r - fromIntegral numStdDev * relative+ lbFix = r - fromIntegral numStdDev * fixed getRankUB :: Int -> Int -> Double -> Int -> Bool -> Word64 -> Double-getRankUB k levels rank numStdDev hra totalN = if exactRank k levels rank hra totalN- then rank+getRankUB k_ levels r numStdDev hra totalN_ =+ if exactRank k_ levels r hra totalN_ then r else min ubRel ubFix where- relative = relRseFactor / fromIntegral k * (if hra then 1.0 - rank else rank)- fixed = fixRseFactor / fromIntegral k- ubRel = rank + fromIntegral numStdDev * relative- ubFix = rank + fromIntegral numStdDev * fixed+ relative = relRseFactor / fromIntegral k_ * (if hra then 1.0 - r else r)+ fixed = fixRseFactor / fromIntegral k_+ ubRel = r + fromIntegral numStdDev * relative+ ubFix = r + fromIntegral numStdDev * fixed exactRank :: Int -> Int -> Double -> Bool -> Word64 -> Bool-exactRank k levels rank hra totalN = (levels == 1 || fromIntegral totalN <= baseCap) || (hra && rank >= 1.0 - exactRankThresh || not hra && rank <= exactRankThresh)+exactRank k_ levels r hra totalN_ =+ (levels == 1 || fromIntegral totalN_ <= baseCap)+ || (hra && r >= 1.0 - exactRankThresh || not hra && r <= exactRankThresh) where- baseCap = k * initNumberOfSections+ baseCap = k_ * initNumberOfSections exactRankThresh :: Double- exactRankThresh = fromIntegral baseCap / fromIntegral totalN+ exactRankThresh = fromIntegral baseCap / fromIntegral totalN_
+ test/BugFixSpec.hs view
@@ -0,0 +1,157 @@+{-# LANGUAGE TypeApplications #-}+module BugFixSpec where++import Control.Monad (forM_, replicateM_)+import Control.Monad.Primitive (PrimState)+import Data.Word+import qualified Data.Vector.Unboxed.Mutable as MUVector+import Test.Hspec+import DataSketches.Quantiles.RelativeErrorQuantile hiding (null, minimum, maximum)+import qualified DataSketches.Quantiles.RelativeErrorQuantile as REQ+import DataSketches.Quantiles.RelativeErrorQuantile.Types+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer++spec :: Spec+spec = do+ describe "Bug fix: off-by-one in getCountWithCriterion (spaceAtBottom=False)" $ do+ -- These tests exercise the Haskell-internal DoubleBuffer directly.+ specify "getCountWithCriterion reads only within active region (spaceAtBottom=False)" $ do+ buf <- mkBuffer 16 0 False+ mapM_ (append buf) [10, 20, 30, 40, 50]+ sort buf++ vec <- getVector buf+ forM_ [5..15] $ \i -> MUVector.unsafeWrite vec i 0++ cnt <- getCountWithCriterion buf 55 (:<)+ cnt `shouldBe` 5++ specify "getCountWithCriterion at upper boundary (spaceAtBottom=False)" $ do+ buf <- mkBuffer 16 0 False+ mapM_ (append buf) [10, 20, 30, 40, 50]+ sort buf++ vec <- getVector buf+ forM_ [5..15] $ \i -> MUVector.unsafeWrite vec i 25++ cnt <- getCountWithCriterion buf 50 (:<)+ cnt `shouldBe` 4++ cnt2 <- getCountWithCriterion buf 50 (:<=)+ cnt2 `shouldBe` 5++ specify "getCountWithCriterion at lower boundary (spaceAtBottom=False)" $ do+ buf <- mkBuffer 16 0 False+ mapM_ (append buf) [10, 20, 30, 40, 50]+ sort buf++ vec <- getVector buf+ forM_ [5..15] $ \i -> MUVector.unsafeWrite vec i 0++ cnt <- getCountWithCriterion buf 5 (:<)+ cnt `shouldBe` 0++ cnt2 <- getCountWithCriterion buf 10 (:<)+ cnt2 `shouldBe` 0++ cnt3 <- getCountWithCriterion buf 10 (:<=)+ cnt3 `shouldBe` 1++ specify "rank of value above max is correct with LowRanksAreAccurate" $ do+ sk <- mkReqSketch 12 LowRanksAreAccurate+ forM_ [1..50 :: Double] $ insert sk+ r <- rank sk 100+ r `shouldBe` 1.0++ specify "rank of value below min is 0 with LowRanksAreAccurate" $ do+ sk <- mkReqSketch 12 LowRanksAreAccurate+ forM_ [10..60 :: Double] $ insert sk+ r <- rank sk 5+ r `shouldBe` 0.0++ specify "ranks are correct for known values with LowRanksAreAccurate" $ do+ sk <- mkReqSketch 50 LowRanksAreAccurate+ let values = [5, 5, 5, 6, 6, 6, 7, 8, 8, 8 :: Double]+ mapM_ (insert sk) values+ r5 <- rank sk 5+ r5 `shouldBe` 0.0+ r6 <- rank sk 6+ r6 `shouldBe` 0.3+ r7 <- rank sk 7+ r7 `shouldBe` 0.6+ r8 <- rank sk 8+ r8 `shouldBe` 0.7+ r9 <- rank sk 9+ r9 `shouldBe` 1.0++ describe "Bug fix: growUntil in merge didn't loop" $ do+ specify "merge into empty sketch preserves count from multi-level source" $ do+ skBig <- mkReqSketch 6 HighRanksAreAccurate+ forM_ [1..2000 :: Double] $ insert skBig+ bigLevels <- numLevels skBig+ bigLevels `shouldSatisfy` (>= 3)++ skSmall <- mkReqSketch 6 HighRanksAreAccurate+ _ <- merge skSmall skBig++ mergedCount <- count skSmall+ mergedCount `shouldBe` 2000++ smallLevels <- numLevels skSmall+ smallLevels `shouldSatisfy` (>= bigLevels)++ specify "merge sketch with 4+ levels into 1-level sketch" $ do+ skSrc <- mkReqSketch 6 HighRanksAreAccurate+ forM_ [1..5000 :: Double] $ insert skSrc+ srcLevels <- numLevels skSrc+ srcLevels `shouldSatisfy` (>= 4)++ skDst <- mkReqSketch 6 HighRanksAreAccurate+ insert skDst 0+ dstLevelsBefore <- numLevels skDst+ dstLevelsBefore `shouldBe` 1++ _ <- merge skDst skSrc+ dstLevelsAfter <- numLevels skDst+ dstLevelsAfter `shouldSatisfy` (>= srcLevels)+ mergedCount <- count skDst+ mergedCount `shouldBe` 5001++ describe "Bug fix: max value comparison inverted in merge" $ do+ specify "merge updates maximum when other has larger max" $ do+ sk1 <- mkReqSketch 12 HighRanksAreAccurate+ forM_ [1..10 :: Double] $ insert sk1+ max1 <- REQ.maximum sk1+ max1 `shouldBe` 10.0++ sk2 <- mkReqSketch 12 HighRanksAreAccurate+ forM_ [100..200 :: Double] $ insert sk2++ _ <- merge sk1 sk2+ mergedMax <- REQ.maximum sk1+ mergedMax `shouldBe` 200.0++ specify "merge does NOT update maximum when other has smaller max" $ do+ sk1 <- mkReqSketch 12 HighRanksAreAccurate+ forM_ [100..200 :: Double] $ insert sk1++ sk2 <- mkReqSketch 12 HighRanksAreAccurate+ forM_ [1..10 :: Double] $ insert sk2++ _ <- merge sk1 sk2+ mergedMax <- REQ.maximum sk1+ mergedMax `shouldBe` 200.0++ specify "merge with bug would have produced wrong max (both directions)" $ do+ sk1 <- mkReqSketch 12 HighRanksAreAccurate+ mapM_ (insert sk1) [1, 2, 3 :: Double]++ sk2 <- mkReqSketch 12 HighRanksAreAccurate+ mapM_ (insert sk2) [10, 20, 30 :: Double]++ _ <- merge sk1 sk2+ mergedMax <- REQ.maximum sk1+ mergedMax `shouldBe` 30.0++ mergedMin <- REQ.minimum sk1+ mergedMin `shouldBe` 1.0
+ test/CountMinSpec.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE TypeApplications #-}+module CountMinSpec where++import Control.Monad (forM_)+import Data.Word+import Test.Hspec+import qualified DataSketches.Frequencies.CountMin as CM++spec :: Spec+spec = describe "Count-Min Sketch" $ do+ specify "empty sketch returns 0 for any item" $ do+ sk <- CM.mkCountMinSketch 0.001 0.01+ est <- CM.estimate sk 42+ est `shouldBe` 0++ specify "single insert returns at least 1" $ do+ sk <- CM.mkCountMinSketch 0.001 0.01+ CM.insert sk 42+ est <- CM.estimate sk 42+ est `shouldSatisfy` (>= 1)++ specify "multiple inserts of same item are counted" $ do+ sk <- CM.mkCountMinSketch 0.001 0.01+ forM_ [1..100 :: Int] $ \_ -> CM.insert sk 42+ est <- CM.estimate sk 42+ est `shouldBe` 100++ specify "insertN works correctly" $ do+ sk <- CM.mkCountMinSketch 0.001 0.01+ CM.insertN sk 42 50+ est <- CM.estimate sk 42+ est `shouldBe` 50++ specify "never undercounts" $ do+ sk <- CM.mkCountMinSketch 0.01 0.01+ forM_ [1..1000 :: Word64] $ \i -> CM.insert sk i+ -- Each item inserted once, estimate should be >= 1+ est <- CM.estimate sk 500+ est `shouldSatisfy` (>= 1)++ specify "unseen items may return small positive values" $ do+ sk <- CM.mkCountMinSketch 0.01 0.01+ forM_ [1..10000 :: Word64] $ \i -> CM.insert sk i+ -- Item 99999 was never inserted+ est <- CM.estimate sk 99999+ -- Should be 0 or a small number due to collisions+ est `shouldSatisfy` (<= 200)++ specify "heavy hitter detection" $ do+ sk <- CM.mkCountMinSketch 0.001 0.01+ -- Insert item 42 many times+ CM.insertN sk 42 10000+ -- Insert many other items once each+ forM_ [1..1000 :: Word64] $ \i -> CM.insert sk i+ est42 <- CM.estimate sk 42+ estOther <- CM.estimate sk 500+ est42 `shouldSatisfy` (> estOther)+ est42 `shouldSatisfy` (>= 10000)++ specify "merge combines two sketches" $ do+ sk1 <- CM.mkCountMinSketch 0.001 0.01+ CM.insertN sk1 42 100+ sk2 <- CM.mkCountMinSketch 0.001 0.01+ CM.insertN sk2 42 200+ CM.merge sk1 sk2+ est <- CM.estimate sk1 42+ est `shouldBe` 300++ specify "dimensions match epsilon and delta" $ do+ sk <- CM.mkCountMinSketch 0.001 0.01+ w <- CM.width sk+ w `shouldSatisfy` (> 0)+ d <- CM.depth sk+ d `shouldSatisfy` (> 0)
+ test/CrossValidationSpec.hs view
@@ -0,0 +1,319 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeApplications #-}+module CrossValidationSpec where++import Control.Monad (forM_, when, unless)+import Control.Monad.IO.Class (liftIO)+import Data.Char (isSpace)+import Data.Word+import System.Directory (doesFileExist)+import System.IO (hFlush, hClose, hSetBuffering, BufferMode(..), hPutStrLn, hGetLine)+import System.Process+import Test.Hspec+import qualified Hedgehog as H+import qualified Hedgehog.Gen as Gen+import qualified Hedgehog.Range as Range++import qualified DataSketches.Quantiles.RelativeErrorQuantile as REQ+import qualified DataSketches.Quantiles.KLL as KLL++findHarnessDir :: IO (Maybe FilePath)+findHarnessDir = do+ let candidates = ["java-harness", "../java-harness"]+ go candidates+ where+ go [] = pure Nothing+ go (d:ds) = do+ exists <- doesFileExist (d ++ "/SketchHarness.class")+ if exists then pure (Just d) else go ds++javaClasspathFor :: FilePath -> String+javaClasspathFor dir = dir ++ ":" ++ dir ++ "/lib/*"++runJavaHarness :: FilePath -> [String] -> IO [String]+runJavaHarness harnessDir commands = do+ let cp = CreateProcess+ { cmdspec = RawCommand "java" ["-cp", javaClasspathFor harnessDir, "SketchHarness"]+ , cwd = Nothing+ , env = Nothing+ , std_in = CreatePipe+ , std_out = CreatePipe+ , std_err = Inherit+ , close_fds = False+ , create_group = False+ , delegate_ctlc = False+ , detach_console = False+ , create_new_console = False+ , new_session = False+ , child_group = Nothing+ , child_user = Nothing+ , use_process_jobs = False+ }+ (Just hin, Just hout, _, ph) <- createProcess cp+ hSetBuffering hin LineBuffering+ hSetBuffering hout LineBuffering++ forM_ commands $ \cmd -> hPutStrLn hin cmd+ hFlush hin+ hClose hin++ let readUntilDone acc = do+ line <- hGetLine hout+ if line == "DONE"+ then pure (reverse acc)+ else readUntilDone (line : acc)++ results <- readUntilDone []+ _ <- waitForProcess ph+ pure results++parseJavaDouble :: String -> Double+parseJavaDouble s+ | s == "NaN" = 0/0+ | otherwise = read (trim s)+ where trim = reverse . dropWhile isSpace . reverse . dropWhile isSpace++spec :: Spec+spec = do+ mDir <- runIO findHarnessDir+ case mDir of+ Nothing ->+ specify "Java harness not found (skipping cross-validation)" $ pendingWith+ "Compile java-harness/SketchHarness.java first"+ Just harnessDir -> do+ describe "REQ Sketch cross-validation with Java" $+ reqCrossValidation harnessDir+ describe "KLL Sketch cross-validation with Java" $+ kllCrossValidation harnessDir++-- | Generate integer-valued doubles in a range. These are exactly+-- representable in both float and double, eliminating precision as a+-- source of disagreement between the Java (float) and Haskell (double)+-- REQ sketches.+genIntDouble :: H.Range Int -> H.Gen Double+genIntDouble r = fromIntegral <$> Gen.int r++-- REQ sketch: property tests comparing Haskell vs Java.+--+-- Java REQ uses float (32-bit); Haskell uses Double (64-bit).+-- We use integer-valued doubles so both representations are identical+-- and exact-mode results must match exactly.+reqCrossValidation :: FilePath -> Spec+reqCrossValidation harnessDir = do++ specify "REQ: exact mode count/min/max match Java (HighRanksAreAccurate)" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.linear 1 50) $+ genIntDouble (Range.linear 1 1000)+ let k = 50 :: Word32+ liftIO $ do+ sk <- REQ.mkReqSketch k REQ.HighRanksAreAccurate+ forM_ values $ REQ.insert sk+ hCount <- REQ.count sk+ hMin <- REQ.minimum sk+ hMax <- REQ.maximum sk++ let jCmds =+ [ "REQ " ++ show k ++ " hra lt"+ , "INSERT " ++ unwords (fmap show values)+ , "COUNT"+ , "MIN"+ , "MAX"+ , "END"+ ]+ jResults <- runJavaHarness harnessDir jCmds+ let jCount = read @Word64 (jResults !! 0)+ jMin = parseJavaDouble (jResults !! 1)+ jMax = parseJavaDouble (jResults !! 2)++ hCount `shouldBe` jCount+ hMin `shouldBe` jMin+ hMax `shouldBe` jMax++ specify "REQ: exact mode ranks match Java exactly (HighRanksAreAccurate, <)" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.linear 10 50) $+ genIntDouble (Range.linear 1 200)+ queryValues <- H.forAll $ Gen.list (Range.linear 1 5) $+ genIntDouble (Range.linear 0 210)+ let k = 50 :: Word32+ liftIO $ do+ sk <- REQ.mkReqSketch k REQ.HighRanksAreAccurate+ forM_ values $ REQ.insert sk+ hRanks <- mapM (REQ.rank sk) queryValues++ let jCmds =+ [ "REQ " ++ show k ++ " hra lt"+ , "INSERT " ++ unwords (fmap show values)+ ] +++ fmap (\q -> "RANK " ++ show q) queryValues +++ [ "END" ]+ jResults <- runJavaHarness harnessDir jCmds+ let jRanks = fmap parseJavaDouble jResults++ forM_ (zip3 queryValues hRanks jRanks) $ \(qv, hr, jr) ->+ unless (isNaN hr && isNaN jr) $ do+ hr `shouldBe` jr++ specify "REQ: exact mode ranks match Java exactly (LowRanksAreAccurate, <)" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.linear 10 50) $+ genIntDouble (Range.linear 1 200)+ queryValues <- H.forAll $ Gen.list (Range.linear 1 5) $+ genIntDouble (Range.linear 0 210)+ let k = 50 :: Word32+ liftIO $ do+ sk <- REQ.mkReqSketch k REQ.LowRanksAreAccurate+ forM_ values $ REQ.insert sk+ hRanks <- mapM (REQ.rank sk) queryValues++ let jCmds =+ [ "REQ " ++ show k ++ " lra lt"+ , "INSERT " ++ unwords (fmap show values)+ ] +++ fmap (\q -> "RANK " ++ show q) queryValues +++ [ "END" ]+ jResults <- runJavaHarness harnessDir jCmds+ let jRanks = fmap parseJavaDouble jResults++ forM_ (zip3 queryValues hRanks jRanks) $ \(qv, hr, jr) ->+ unless (isNaN hr && isNaN jr) $ do+ hr `shouldBe` jr++ specify "REQ: estimation mode ranks are both valid approximations (k=6, 200 items)" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.singleton 200) $+ genIntDouble (Range.linear 1 1000)+ queryValues <- H.forAll $ Gen.list (Range.linear 1 3) $+ genIntDouble (Range.linear 1 1000)+ let k = 6 :: Word32+ liftIO $ do+ sk <- REQ.mkReqSketch k REQ.HighRanksAreAccurate+ forM_ values $ REQ.insert sk+ hRanks <- mapM (REQ.rank sk) queryValues++ let jCmds =+ [ "REQ " ++ show k ++ " hra lt"+ , "INSERT " ++ unwords (fmap show values)+ ] +++ fmap (\q -> "RANK " ++ show q) queryValues +++ [ "END" ]+ jResults <- runJavaHarness harnessDir jCmds+ let jRanks = fmap parseJavaDouble jResults++ -- In estimation mode, each implementation made independent random+ -- compaction decisions (different RNG seeds), so the retained items+ -- differ. Both ranks should approximate the true rank, and we check+ -- that the two approximations don't diverge beyond 2x the sketch's+ -- error bound (~14% per side for k=6).+ let n = fromIntegral (length values) :: Double+ forM_ (zip3 queryValues hRanks jRanks) $ \(qv, hr, jr) ->+ unless (isNaN hr && isNaN jr) $ do+ let trueRank = fromIntegral (length (filter (< qv) values)) / n+ assertWithinBound ("Haskell rank of " ++ show qv) 0.15 hr trueRank+ assertWithinBound ("Java rank of " ++ show qv) 0.15 jr trueRank++-- KLL sketch cross-validation.+--+-- Both Java and Haskell KLL use doubles, so all values are identical+-- in both representations. In exact mode results must match exactly.+kllCrossValidation :: FilePath -> Spec+kllCrossValidation harnessDir = do++ specify "KLL: exact mode count/min/max match Java exactly" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.linear 1 100) $+ Gen.double (Range.linearFrac 1 1000)+ let k = 200+ liftIO $ do+ sk <- KLL.mkKllSketch k+ forM_ values $ KLL.insert sk+ hCount <- KLL.count sk+ hMin <- KLL.minimum sk+ hMax <- KLL.maximum sk++ let jCmds =+ [ "KLL " ++ show k+ , "INSERT " ++ unwords (fmap show values)+ , "COUNT"+ , "MIN"+ , "MAX"+ , "END"+ ]+ jResults <- runJavaHarness harnessDir jCmds+ let jCount = read @Word64 (jResults !! 0)+ jMin = parseJavaDouble (jResults !! 1)+ jMax = parseJavaDouble (jResults !! 2)++ hCount `shouldBe` jCount+ hMin `shouldBe` jMin+ hMax `shouldBe` jMax++ specify "KLL: exact mode ranks match Java exactly (k=200, few items)" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.linear 10 50) $+ Gen.double (Range.linearFrac 1 100)+ queryValues <- H.forAll $ Gen.list (Range.linear 1 5) $+ Gen.double (Range.linearFrac 0 110)+ let k = 200+ liftIO $ do+ sk <- KLL.mkKllSketch k+ forM_ values $ KLL.insert sk+ hRanks <- mapM (KLL.rank sk) queryValues++ let jCmds =+ [ "KLL " ++ show k+ , "INSERT " ++ unwords (fmap show values)+ ] +++ fmap (\q -> "RANK " ++ show q) queryValues +++ [ "END" ]+ jResults <- runJavaHarness harnessDir jCmds+ let jRanks = fmap parseJavaDouble jResults++ forM_ (zip3 queryValues hRanks jRanks) $ \(qv, hr, jr) ->+ unless (isNaN hr && isNaN jr) $ do+ hr `shouldBe` jr++ specify "KLL: estimation mode ranks are both valid approximations (k=200, 500 items)" $ hedgehog $+ H.property $ do+ values <- H.forAll $ Gen.list (Range.singleton 500) $+ Gen.double (Range.linearFrac 1 1000)+ queryValues <- H.forAll $ Gen.list (Range.linear 1 3) $+ Gen.double (Range.linearFrac 1 1000)+ let k = 200+ liftIO $ do+ sk <- KLL.mkKllSketch k+ forM_ values $ KLL.insert sk+ hRanks <- mapM (KLL.rank sk) queryValues++ let jCmds =+ [ "KLL " ++ show k+ , "INSERT " ++ unwords (fmap show values)+ ] +++ fmap (\q -> "RANK " ++ show q) queryValues +++ [ "END" ]+ jResults <- runJavaHarness harnessDir jCmds+ let jRanks = fmap parseJavaDouble jResults++ -- Different compaction decisions from different RNGs. Verify both+ -- approximate the true rank rather than comparing to each other.+ let n = fromIntegral (length values) :: Double+ forM_ (zip3 queryValues hRanks jRanks) $ \(qv, hr, jr) ->+ unless (isNaN hr && isNaN jr) $ do+ let trueRank = fromIntegral (length (filter (< qv) values)) / n+ assertWithinBound ("Haskell rank of " ++ show qv) 0.05 hr trueRank+ assertWithinBound ("Java rank of " ++ show qv) 0.05 jr trueRank++-- | Assert that actual is within tolerance of expected.+assertWithinBound :: String -> Double -> Double -> Double -> IO ()+assertWithinBound label tolerance actual expected =+ when (abs (actual - expected) > tolerance) $+ expectationFailure $ label ++ ": expected " ++ show expected+ ++ " +/- " ++ show tolerance+ ++ " but got " ++ show actual+ ++ " (delta=" ++ show (abs (actual - expected)) ++ ")"++hedgehog :: H.Property -> IO ()+hedgehog prop = do+ result <- H.check prop+ unless result $ expectationFailure "Hedgehog property failed"
+ test/HyperLogLogSpec.hs view
@@ -0,0 +1,73 @@+{-# LANGUAGE TypeApplications #-}+module HyperLogLogSpec where++import Control.Monad (forM_)+import Data.Word+import Test.Hspec+import qualified DataSketches.Distinct.HyperLogLog as HLL++spec :: Spec+spec = describe "HyperLogLog Sketch" $ do+ specify "empty sketch estimates 0" $ do+ sk <- HLL.mkHllSketch 12+ est <- HLL.estimate sk+ est `shouldBe` 0.0++ specify "single item estimates ~1" $ do+ sk <- HLL.mkHllSketch 12+ HLL.insert sk 42+ est <- HLL.estimate sk+ est `shouldSatisfy` (\e -> e >= 0.5 && e <= 2.0)++ specify "100 distinct items estimates approximately 100" $ do+ sk <- HLL.mkHllSketch 12+ forM_ [1..100 :: Word64] $ HLL.insert sk+ est <- HLL.estimate sk+ -- With p=12, error ~1.6%, so within ~20 is very generous+ est `shouldSatisfy` (\e -> abs (e - 100) < 30)++ specify "1000 distinct items estimates approximately 1000" $ do+ sk <- HLL.mkHllSketch 14+ forM_ [1..1000 :: Word64] $ HLL.insert sk+ est <- HLL.estimate sk+ est `shouldSatisfy` (\e -> abs (e - 1000) < 100)++ specify "10000 distinct items estimates approximately 10000" $ do+ sk <- HLL.mkHllSketch 14+ forM_ [1..10000 :: Word64] $ HLL.insert sk+ est <- HLL.estimate sk+ est `shouldSatisfy` (\e -> abs (e - 10000) < 1000)++ specify "duplicate items don't increase estimate" $ do+ sk <- HLL.mkHllSketch 12+ forM_ [1..100 :: Word64] $ HLL.insert sk+ est1 <- HLL.estimate sk+ -- Insert all the same items again+ forM_ [1..100 :: Word64] $ HLL.insert sk+ est2 <- HLL.estimate sk+ -- Estimate should be about the same (still ~100)+ abs (est2 - est1) `shouldSatisfy` (< 5)++ specify "merge combines two sketches" $ do+ sk1 <- HLL.mkHllSketch 12+ forM_ [1..500 :: Word64] $ HLL.insert sk1+ sk2 <- HLL.mkHllSketch 12+ forM_ [501..1000 :: Word64] $ HLL.insert sk2+ HLL.merge sk1 sk2+ est <- HLL.estimate sk1+ est `shouldSatisfy` (\e -> abs (e - 1000) < 100)++ specify "merge with overlapping items gives correct estimate" $ do+ sk1 <- HLL.mkHllSketch 12+ forM_ [1..500 :: Word64] $ HLL.insert sk1+ sk2 <- HLL.mkHllSketch 12+ forM_ [250..750 :: Word64] $ HLL.insert sk2+ HLL.merge sk1 sk2+ est <- HLL.estimate sk1+ -- 750 distinct items+ est `shouldSatisfy` (\e -> abs (e - 750) < 100)++ specify "precision is reported correctly" $ do+ sk <- HLL.mkHllSketch 14+ p <- HLL.precision sk+ p `shouldBe` 14
+ test/KllSpec.hs view
@@ -0,0 +1,117 @@+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE LambdaCase #-}+module KllSpec where++import Control.Monad (forM_)+import Control.Monad.Primitive (PrimState)+import Data.Word+import Test.Hspec+import Test.QuickCheck+import qualified DataSketches.Quantiles.KLL as KLL++spec :: Spec+spec = describe "KLL Sketch" $ do+ specify "empty sketch returns NaN for quantile" $ do+ sk <- KLL.mkKllSketch 200+ q <- KLL.quantile sk 0.5+ isNaN q `shouldBe` True++ specify "empty sketch returns True for null" $ do+ sk <- KLL.mkKllSketch 200+ KLL.null sk `shouldReturn` True++ specify "single element sketch returns the element for all quantiles" $ do+ sk <- KLL.mkKllSketch 200+ KLL.insert sk 42.0+ KLL.null sk `shouldReturn` False+ KLL.count sk `shouldReturn` 1+ KLL.minimum sk `shouldReturn` 42.0+ KLL.maximum sk `shouldReturn` 42.0+ q <- KLL.quantile sk 0.5+ q `shouldBe` 42.0++ specify "tracks min and max correctly" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [10, 5, 20, 1, 100 :: Double] $ KLL.insert sk+ KLL.minimum sk `shouldReturn` 1.0+ KLL.maximum sk `shouldReturn` 100.0+ KLL.count sk `shouldReturn` 5++ specify "ignores NaN values" $ do+ sk <- KLL.mkKllSketch 200+ KLL.insert sk (0/0)+ KLL.null sk `shouldReturn` True++ specify "sequential insert produces correct count" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..1000 :: Double] $ KLL.insert sk+ KLL.count sk `shouldReturn` 1000++ specify "quantile 0 returns minimum" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..100 :: Double] $ KLL.insert sk+ q <- KLL.quantile sk 0.0+ q `shouldBe` 1.0++ specify "quantile 1 returns maximum" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..100 :: Double] $ KLL.insert sk+ q <- KLL.quantile sk 1.0+ q `shouldBe` 100.0++ specify "rank is monotonically non-decreasing" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..200 :: Double] $ KLL.insert sk+ rs <- KLL.ranks sk [10, 50, 100, 150, 190]+ rs `shouldSatisfy` isNonDecreasing++ specify "median of sequential data is approximately correct" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..1000 :: Double] $ KLL.insert sk+ median <- KLL.quantile sk 0.5+ -- With k=200, error ~1.3%, so median should be within ~13 of 500+ median `shouldSatisfy` (\m -> abs (m - 500) < 50)++ specify "merge combines two sketches correctly" $ do+ sk1 <- KLL.mkKllSketch 200+ forM_ [1..100 :: Double] $ KLL.insert sk1+ sk2 <- KLL.mkKllSketch 200+ forM_ [101..200 :: Double] $ KLL.insert sk2+ KLL.merge sk1 sk2+ KLL.count sk1 `shouldReturn` 200+ KLL.minimum sk1 `shouldReturn` 1.0+ KLL.maximum sk1 `shouldReturn` 200.0++ specify "merge with empty sketch is identity" $ do+ sk1 <- KLL.mkKllSketch 200+ forM_ [1..50 :: Double] $ KLL.insert sk1+ sk2 <- KLL.mkKllSketch 200+ KLL.merge sk1 sk2+ KLL.count sk1 `shouldReturn` 50++ specify "large insert triggers compaction" $ do+ sk <- KLL.mkKllSketch 8+ forM_ [1..10000 :: Double] $ KLL.insert sk+ retained <- KLL.retainedItemCount sk+ retained `shouldSatisfy` (< 10000)+ KLL.count sk `shouldReturn` 10000++ specify "CDF returns Just for non-empty sketch" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..100 :: Double] $ KLL.insert sk+ cdf <- KLL.cumulativeDistributionFunction sk [25, 50, 75]+ cdf `shouldSatisfy` \case+ Just xs -> length xs == 4 && all (\x -> x >= 0 && x <= 1) xs+ Nothing -> False++ specify "PMF sums approximately to 1" $ do+ sk <- KLL.mkKllSketch 200+ forM_ [1..100 :: Double] $ KLL.insert sk+ pmf <- KLL.probabilityMassFunction sk [25, 50, 75]+ length pmf `shouldBe` 4+ sum pmf `shouldSatisfy` (\s -> abs (s - 1.0) < 0.01)++isNonDecreasing :: (Ord a) => [a] -> Bool+isNonDecreasing [] = True+isNonDecreasing [_] = True+isNonDecreasing (x:y:rest) = x <= y && isNonDecreasing (y:rest)
test/RelativeErrorQuantileSpec.hs view
@@ -2,27 +2,21 @@ {-# LANGUAGE NumericUnderscores #-} module RelativeErrorQuantileSpec where -import Data.Primitive.MutVar-import qualified Data.Vector.Unboxed as U import Control.Monad import Control.Monad.Primitive import DataSketches.Quantiles.RelativeErrorQuantile import DataSketches.Quantiles.RelativeErrorQuantile.Types-import DataSketches.Quantiles.RelativeErrorQuantile.Internal-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary-import Data.List hiding (insert)-import Data.Maybe (fromJust, isJust)+import Data.List hiding (insert, null)+import qualified Data.List import Data.Word import Test.Hspec-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer (DoubleIsNonFiniteException(..))-import Text.Show.Pretty spec :: Spec spec = do specify "non finite PMF/CDF should throw" $ asIO $ do sk <- mkReqSketch 6 HighRanksAreAccurate insert sk 1- cumulativeDistributionFunction sk [0 / 0] `shouldThrow` (== CumulativeDistributionInvariantsSplitsAreNotFinite)+ cumulativeDistributionFunction sk [0 / 0] `shouldThrow` (\(DoubleIsNonFiniteException _) -> True) specify "updating a sketch with NaN should ignore it" $ asIO $ do sk <- mkReqSketch 6 HighRanksAreAccurate insert sk (0 / 0)@@ -44,8 +38,8 @@ n `shouldBe` 21 mapM_ (insert sk2) [16..300] merge sk1 sk2- n <- count sk1- n `shouldBe` 321+ n' <- count sk1+ n' `shouldBe` 321 describe "property tests" $ do specify "ReqSketch quantile estimates are within ε bounds compared to real quantile calculations" $ print () specify "merging N ReqSketches is equivalent +/- ε to inserting the same values into 1 ReqSketch" $ print ()@@ -71,27 +65,20 @@ actualRanks <- mapM (rank sk) simpleTestValues actualRanks `shouldBe` lessThanRs describe "<=" $ do- let mkSk' sk = sk { criterion = (:<=) } :: ReqSketch (PrimState IO)- specify "ranks function should match lessThanRs" $ \sk -> do- let sk' = mkSk' sk- actualRanks <- ranks sk' simpleTestValues+ specify "ranks function should match lessThanEqRs" $ \sk -> do+ setCriterionLE sk+ actualRanks <- ranks sk simpleTestValues actualRanks `shouldBe` lessThanEqRs specify "mapM rank should match ranks behaviour" $ \sk -> do- let sk' = mkSk' sk- actualRanks <- mapM (rank sk') simpleTestValues+ setCriterionLE sk+ actualRanks <- mapM (rank sk) simpleTestValues actualRanks `shouldBe` lessThanEqRs ------ -- k min max hra lteq low-to-high or high-to-low- bigTest 6 1 200 HighRanksAreAccurate (:<=) True- bigTest 6 1 200 LowRanksAreAccurate (:<=) True- bigTest 6 1 200 HighRanksAreAccurate (:<) False- bigTest 6 1 200 LowRanksAreAccurate (:<) True+ bigTest 6 1 200 HighRanksAreAccurate False True+ bigTest 6 1 200 LowRanksAreAccurate False True+ bigTest 6 1 200 HighRanksAreAccurate True False+ bigTest 6 1 200 LowRanksAreAccurate True True mergeSpec @@ -101,30 +88,28 @@ insert sk 1 insert sk 1 insert sk 1- r <- quantile sk 0.5 + r <- quantile sk 0.5 r `shouldBe` 1.0 -bigTest :: Word32 -> Int -> Int -> RankAccuracy -> Criterion -> Bool -> Spec-bigTest k min_ max_ hra crit up = do+bigTest :: Word32 -> Int -> Int -> RankAccuracy -> Bool -> Bool -> Spec+bigTest k min_ max_ hra useLe up = do let testName = unwords [ "k=" <> show k , "min=" <> show min_ , "max=" <> show max_ , "hra=" <> show hra+ , "le=" <> show useLe , "up=" <> show up ] testContents :: IO () testContents = do- sk <- loadSketch k min_ max_ hra crit up- checkAux sk+ sk <- loadSketch k min_ max_ hra useLe up checkGetRank sk min_ max_ checkGetRanks sk max_ checkGetQuantiles sk checkGetCDF sk checkGetPMF sk- -- checkIterator sk- -- checkMerge sk it testName testContents asIO :: IO a -> IO a@@ -144,33 +129,15 @@ s `merge` s3 pure () -loadSketch :: Word32 -> Int -> Int -> RankAccuracy -> Criterion -> Bool -> IO (ReqSketch (PrimState IO))-loadSketch k min_ max_ hra ltEq up = do+loadSketch :: Word32 -> Int -> Int -> RankAccuracy -> Bool -> Bool -> IO (ReqSketch (PrimState IO))+loadSketch k min_ max_ hra useLe up = do sk <- mkReqSketch k hra :: IO (ReqSketch (PrimState IO))- -- This just seems geared at making sure that ranks come out right regardless of order+ when useLe $ setCriterionLE sk mapM_ (insert sk . fromIntegral) $ if up then [min_ .. max_]- else reverse [min_ .. max_ {- + 1 -}]+ else reverse [min_ .. max_] pure sk -checkAux :: ReqSketch (PrimState IO) -> IO ()-checkAux sk = do- auxiliary <- mkAuxiliaryFromReqSketch sk- totalCount <- computeTotalRetainedItems sk-- let rows = raWeightedItems auxiliary- getRow = (U.!)- let initialRow = getRow rows 0- otherRows = map (getRow rows) [1..totalCount - 1]- foldM_- (\lastRow thisRow -> do- fst thisRow `shouldSatisfy` (>= fst lastRow)- snd thisRow `shouldSatisfy` (>= snd lastRow)- pure thisRow- )- initialRow- otherRows- checkGetRank :: ReqSketch (PrimState IO) -> Int -> Int -> IO () checkGetRank sk min_ max_ = do let (v : spArr) = evenlySpacedFloats 0 (fromIntegral max_) 11@@ -194,7 +161,7 @@ checkGetCDF sk = do let spArr = [20, 40 .. 180] r <- cumulativeDistributionFunction sk spArr- r `shouldSatisfy` isJust+ r `shouldSatisfy` (\(Just _) -> True) checkGetPMF :: ReqSketch (PrimState IO) -> IO () checkGetPMF sk = do@@ -202,50 +169,11 @@ r <- probabilityMassFunction sk spArr r `shouldNotSatisfy` Data.List.null -{--checkMerge :: ReqSketch n (PrimState IO) -> IO ()-checkMerge sk = do- sk' <- copyRs--}--checkGetRankConcreteExample :: RankAccuracy -> Criterion -> IO ()-checkGetRankConcreteExample ra crit = do- sk <- loadSketch 12 1 1000 ra crit True- rLB <- rankLowerBound sk 0.5 1- rLB `shouldSatisfy` (> 0)- rLB <- case ra of- HighRanksAreAccurate -> rankLowerBound sk (995 / 1000) 1- LowRanksAreAccurate -> rankLowerBound sk (5 / 1000) 1- rLB `shouldSatisfy` (> 0)- rUB <- rankUpperBound sk 0.5 1- rUB `shouldSatisfy` (> 0)- rUB <- case ra of- HighRanksAreAccurate -> rankUpperBound sk (995 / 1000) 1- LowRanksAreAccurate -> rankUpperBound sk (5 / 1000) 1- rUB `shouldSatisfy` (> 0)- void $ ranks sk [5, 100]------checkGetQuantiles- :: ReqSketch (PrimState IO)- -> IO ()+checkGetQuantiles :: ReqSketch (PrimState IO) -> IO () checkGetQuantiles sk = do let rArr = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]- -- nothing getting checked here apparently, guess the thing not- -- exploding is sufficient.- qOut <- quantiles sk rArr- pure ()+ void $ quantiles sk rArr --- | Returns a float array of evenly spaced values between value1 and value2 inclusive.--- If value2 > value1, the resulting sequence will be increasing.--- If value2 < value1, the resulting sequence will be decreasing.--- value1 will be in index 0 of the returned array--- value2 will be in the highest index of the returned array--- valu3 is the total number of values including value1 and value2. Must be 2 or greater.--- returns a float array of evenly spaced values between value1 and value2 inclusive. evenlySpacedFloats :: Double -> Double -> Word -> [Double] evenlySpacedFloats _ _ 0 = error "Needs at least two steps" evenlySpacedFloats _ _ 1 = error "Needs at least two steps"
test/Spec.hs view
@@ -1,1 +1,39 @@-{-# OPTIONS_GHC -F -pgmF hspec-discover #-}+import Test.Hspec+import qualified AuxiliarySpec+import qualified CompactorSpec+import qualified DoubleBufferSpec+import qualified ProofCheckSpec+import qualified RelativeErrorQuantileSpec+import qualified KllSpec+import qualified HyperLogLogSpec+import qualified ThetaSpec+import qualified CountMinSpec+import qualified BugFixSpec+import qualified CrossValidationSpec+import System.Environment+import Test.HSpec.JUnit+import Test.Hspec.Runner++main :: IO ()+main = + getArgs+ >>= readConfig config+ >>= withArgs [] . runSpec specs+ >>= evaluateSummary+ where+ config = defaultConfig + { configFormat = Just $ junitFormat "test-results.xml" "data-sketches" + }+ specs = do+ describe "Auxiliary" AuxiliarySpec.spec+ describe "Compactor" CompactorSpec.spec+ describe "DoubleBuffer" DoubleBufferSpec.spec+ describe "ProofCheck" ProofCheckSpec.spec+ describe "RelativeErrorQuantile" RelativeErrorQuantileSpec.spec+ describe "KLL" KllSpec.spec+ describe "HyperLogLog" HyperLogLogSpec.spec+ describe "Theta" ThetaSpec.spec+ describe "CountMin" CountMinSpec.spec+ describe "BugFix" BugFixSpec.spec+ describe "CrossValidation" CrossValidationSpec.spec+
+ test/ThetaSpec.hs view
@@ -0,0 +1,113 @@+{-# LANGUAGE TypeApplications #-}+module ThetaSpec where++import Control.Monad (forM_)+import Data.Word+import Test.Hspec+import qualified DataSketches.Distinct.Theta as Theta++spec :: Spec+spec = describe "Theta Sketch" $ do+ specify "empty sketch estimates 0" $ do+ sk <- Theta.mkThetaSketch 128+ est <- Theta.estimate sk+ est `shouldBe` 0.0+ Theta.isEmpty sk `shouldReturn` True++ specify "single item estimates 1" $ do+ sk <- Theta.mkThetaSketch 128+ Theta.insert sk 42+ est <- Theta.estimate sk+ est `shouldBe` 1.0+ Theta.isEmpty sk `shouldReturn` False++ specify "exact mode for small cardinality" $ do+ sk <- Theta.mkThetaSketch 4096+ forM_ [1..100 :: Word64] $ Theta.insert sk+ est <- Theta.estimate sk+ est `shouldBe` 100.0++ specify "duplicate items don't increase estimate" $ do+ sk <- Theta.mkThetaSketch 4096+ forM_ [1..50 :: Word64] $ Theta.insert sk+ est1 <- Theta.estimate sk+ forM_ [1..50 :: Word64] $ Theta.insert sk+ est2 <- Theta.estimate sk+ est2 `shouldBe` est1++ specify "1000 distinct items estimates approximately 1000" $ do+ sk <- Theta.mkThetaSketch 4096+ forM_ [1..1000 :: Word64] $ Theta.insert sk+ est <- Theta.estimate sk+ est `shouldSatisfy` (\e -> abs (e - 1000) < 100)++ specify "10000 items with smaller k estimates approximately" $ do+ sk <- Theta.mkThetaSketch 256+ forM_ [1..10000 :: Word64] $ Theta.insert sk+ est <- Theta.estimate sk+ est `shouldSatisfy` (\e -> abs (e - 10000) < 2000)++ describe "union" $ do+ specify "union of disjoint sets" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..500 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ forM_ [501..1000 :: Word64] $ Theta.insert sk2+ u <- Theta.union sk1 sk2+ est <- Theta.estimate u+ est `shouldSatisfy` (\e -> abs (e - 1000) < 100)++ specify "union with empty sketch" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..100 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ u <- Theta.union sk1 sk2+ est <- Theta.estimate u+ est `shouldBe` 100.0++ describe "intersection" $ do+ specify "intersection of overlapping sets" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..200 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ forM_ [101..300 :: Word64] $ Theta.insert sk2+ i <- Theta.intersection sk1 sk2+ est <- Theta.estimate i+ -- 100 items overlap (101..200)+ est `shouldSatisfy` (\e -> abs (e - 100) < 30)++ specify "intersection of disjoint sets is ~0" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..100 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ forM_ [101..200 :: Word64] $ Theta.insert sk2+ i <- Theta.intersection sk1 sk2+ est <- Theta.estimate i+ est `shouldBe` 0.0++ specify "intersection with empty sketch is 0" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..100 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ i <- Theta.intersection sk1 sk2+ est <- Theta.estimate i+ est `shouldBe` 0.0++ describe "difference" $ do+ specify "difference removes overlapping items" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..200 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ forM_ [101..300 :: Word64] $ Theta.insert sk2+ d <- Theta.difference sk1 sk2+ est <- Theta.estimate d+ -- sk1 has 1..200, sk2 has 101..300. A \ B = 1..100+ est `shouldSatisfy` (\e -> abs (e - 100) < 30)++ specify "difference with empty B returns A" $ do+ sk1 <- Theta.mkThetaSketch 4096+ forM_ [1..100 :: Word64] $ Theta.insert sk1+ sk2 <- Theta.mkThetaSketch 4096+ d <- Theta.difference sk1 sk2+ est <- Theta.estimate d+ est `shouldBe` 100.0