data-sketches-core 0.2.0.1 → 0.3.0.0

raw patch · 14 files changed

+231/−1251 lines, 14 filesdep −mwc-randomdep −vector-algorithmsPVP ok

version bump matches the API change (PVP)

Dependencies removed: mwc-random, vector-algorithms

API changes (from Hackage documentation)

Files

• ChangeLog.md +45/−0
• cbits/countmin.c +3/−1
• cbits/hll.c +120/−7
• cbits/kll.c +50/−2
• data-sketches-core.cabal +1/−12
• src/DataSketches/Core/Internal/URef.hs +0/−73
• src/DataSketches/Core/Snapshot.hs +0/−7
• src/DataSketches/Quantiles/RelativeErrorQuantile/Internal.hs +0/−184
• src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs +0/−189
• src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs +0/−285
• src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs +0/−11
• src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs +0/−340
• src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs +0/−121
• src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs +12/−19

@@ -1,5 +1,50 @@ # Changelog for data-sketches-core +## 0.3.0.0 — 2026-04-08++### Breaking changes++- **Removed old pure-Haskell REQ internals**: The following modules have been+  removed. All sketch operations now go through the C backend (`CInternal`),+  which replaced these modules in 0.2.0.0 but left the dead code in the+  package. Downstream consumers should use the public `data-sketches` API.++  - `DataSketches.Core.Internal.URef`+  - `DataSketches.Core.Snapshot`+  - `DataSketches.Quantiles.RelativeErrorQuantile.Internal`+  - `DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary`+  - `DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor`+  - `DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer`+  - `DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch`++- **`Criterion` type removed from `Types` module**: The `Criterion` type and+  its `InequalitySearch` class instance were only used by the deleted Haskell+  internals. Criterion-based queries are handled entirely in C.++- **`DoubleIsNonFiniteException` moved to `Types` module**: Previously exported+  from `Internal.DoubleBuffer`, now exported from+  `DataSketches.Quantiles.RelativeErrorQuantile.Types`.++### Performance++- **HLL estimate: 2.6x faster** — replaced `ldexp(1.0, -val)` (libm call per+  register) with IEEE 754 bit manipulation (`pow2_neg`). Eliminates ~4096+  function calls per estimate at p=12.++- **NEON/SSE2 intrinsics for HLL and KLL** — explicit SIMD paths for+  `hll_c_merge` (16-wide `uint8` max via `vmaxq_u8`/`_mm_max_epu8`),+  `hll_c_estimate` zero-counting (16-wide via `vceqq_u8`/`_mm_cmpeq_epi8`),+  and `kll_rank` inner loop (2-4 wide `double` comparison via+  `vcltq_f64`/`_mm_cmplt_pd`). Scalar fallback on other architectures.++- Added `__restrict__` qualifiers to `cms_c_merge` to guarantee+  auto-vectorization of the bulk `uint64` addition.++### Dependency changes++- Dropped `mwc-random` and `vector-algorithms` dependencies (only needed by+  the removed Haskell internals).+ ## 0.2.0.1  ### Bug fixes

@@ -71,8 +71,10 @@  void cms_c_merge(cms_sketch_t *dst, const cms_sketch_t *src) {     int n = dst->rows * dst->cols;+    uint64_t *__restrict__ dt = dst->table;+    const uint64_t *__restrict__ st = src->table;     for (int i = 0; i < n; i++)-        dst->table[i] += src->table[i];+        dt[i] += st[i];     dst->total_n += src->total_n; } 

@@ -3,6 +3,14 @@ #include <string.h> #include <math.h> +#if defined(__ARM_NEON) || defined(__ARM_NEON__)+#include <arm_neon.h>+#define HLL_HAS_NEON 1+#elif defined(__SSE2__)+#include <emmintrin.h>+#define HLL_HAS_SSE2 1+#endif+ static inline uint64_t hll_murmur_mix64(uint64_t h) {     h ^= h >> 33; h *= 0xFF51AFD7ED558CCDULL;     h ^= h >> 33; h *= 0xC4CEB9FE1A85EC53ULL;@@ -31,9 +39,6 @@ void hll_c_insert(hll_sketch_t *sk, uint64_t item) {     uint64_t hash = hll_murmur_mix64(item);     int reg_idx = (int)(hash & (uint64_t)(sk->m - 1));-    /* Branchless rho: place a sentinel bit at position (64-p) so-       __builtin_ctzll always terminates within the valid range.-       Eliminates two branches from the old w==0 / clz<bits checks. */     int bits = 64 - sk->p;     uint64_t w = (hash >> sk->p) | (1ULL << bits);     uint8_t rho = (uint8_t)(__builtin_ctzll(w) + 1);@@ -46,18 +51,100 @@         hll_c_insert(sk, items[i]); } +/* 2^(-val) via IEEE 754 bit manipulation.+   For val in [0, 1022]: exponent = 1023 - val, mantissa = 0.+   For val >= 1023 (shouldn't happen with HLL registers ≤ 64): returns 0.0. */+static inline double pow2_neg(int val) {+    if (__builtin_expect(val > 1022, 0)) return 0.0;+    union { uint64_t u; double d; } bits;+    bits.u = (uint64_t)(1023 - val) << 52;+    return bits.d;+}+ double hll_c_estimate(const hll_sketch_t *sk) {     int m = sk->m;     double mf = (double)m;     double harmonic_sum = 0.0;     int zero_count = 0;+    const uint8_t *regs = sk->registers; -    /* Auto-vectorizable: straight-line accumulation, no data-dependent branches */+#if HLL_HAS_NEON+    /* NEON: count zeros 16 registers at a time */+    {+        uint8x16_t zero_vec = vdupq_n_u8(0);+        uint8x16_t zcount16 = vdupq_n_u8(0);+        int i = 0;+        int chunks = m & ~15;+        int batch = 0;+        for (; i < chunks; i += 16) {+            uint8x16_t data = vld1q_u8(regs + i);+            uint8x16_t eq = vceqq_u8(data, zero_vec);+            /* eq lanes are 0xFF where zero, 0x00 otherwise.+               Subtracting 0xFF is adding 1 in unsigned wrapping. */+            zcount16 = vsubq_u8(zcount16, eq);+            batch++;+            /* Flush to avoid uint8 overflow (max 255 accumulated) */+            if (batch == 255) {+                uint16x8_t sum16 = vpaddlq_u8(zcount16);+                uint32x4_t sum32 = vpaddlq_u16(sum16);+                uint64x2_t sum64 = vpaddlq_u32(sum32);+                zero_count += (int)(vgetq_lane_u64(sum64, 0) + vgetq_lane_u64(sum64, 1));+                zcount16 = vdupq_n_u8(0);+                batch = 0;+            }+        }+        /* Flush remaining accumulated zeros */+        {+            uint16x8_t sum16 = vpaddlq_u8(zcount16);+            uint32x4_t sum32 = vpaddlq_u16(sum16);+            uint64x2_t sum64 = vpaddlq_u32(sum32);+            zero_count += (int)(vgetq_lane_u64(sum64, 0) + vgetq_lane_u64(sum64, 1));+        }+        /* Scalar tail */+        for (; i < m; i++)+            zero_count += (regs[i] == 0);+    }+    /* Harmonic sum (scalar with pow2_neg — hard to vectorize the+       uint8→double widening chain profitably) */+    for (int i = 0; i < m; i++)+        harmonic_sum += pow2_neg((int)regs[i]);+#elif HLL_HAS_SSE2+    /* SSE2: count zeros 16 registers at a time */+    {+        __m128i zero_vec = _mm_setzero_si128();+        __m128i zcount16 = _mm_setzero_si128();+        int i = 0;+        int chunks = m & ~15;+        int batch = 0;+        for (; i < chunks; i += 16) {+            __m128i data = _mm_loadu_si128((const __m128i *)(regs + i));+            __m128i eq = _mm_cmpeq_epi8(data, zero_vec);+            zcount16 = _mm_sub_epi8(zcount16, eq);+            batch++;+            if (batch == 255) {+                /* Horizontal sum via SAD against zero */+                __m128i sad = _mm_sad_epu8(zcount16, _mm_setzero_si128());+                zero_count += _mm_extract_epi16(sad, 0) + _mm_extract_epi16(sad, 4);+                zcount16 = _mm_setzero_si128();+                batch = 0;+            }+        }+        {+            __m128i sad = _mm_sad_epu8(zcount16, _mm_setzero_si128());+            zero_count += _mm_extract_epi16(sad, 0) + _mm_extract_epi16(sad, 4);+        }+        for (; i < m; i++)+            zero_count += (regs[i] == 0);+    }+    for (int i = 0; i < m; i++)+        harmonic_sum += pow2_neg((int)regs[i]);+#else     for (int i = 0; i < m; i++) {-        int val = (int)sk->registers[i];-        harmonic_sum += ldexp(1.0, -val);+        int val = (int)regs[i];+        harmonic_sum += pow2_neg(val);         zero_count += (val == 0);     }+#endif      double alpha;     if      (m == 16) alpha = 0.673;@@ -76,15 +163,41 @@         return raw; } -/* Branchless register-wise max for auto-vectorization */ void hll_c_merge(hll_sketch_t *dst, const hll_sketch_t *src) {     int m = dst->m;     uint8_t *__restrict__ d = dst->registers;     const uint8_t *__restrict__ s = src->registers;++#if HLL_HAS_NEON+    int i = 0;+    int chunks = m & ~15;+    for (; i < chunks; i += 16) {+        uint8x16_t dv = vld1q_u8(d + i);+        uint8x16_t sv = vld1q_u8(s + i);+        vst1q_u8(d + i, vmaxq_u8(dv, sv));+    }+    for (; i < m; i++) {+        uint8_t sv = s[i], dv = d[i];+        d[i] = sv > dv ? sv : dv;+    }+#elif HLL_HAS_SSE2+    int i = 0;+    int chunks = m & ~15;+    for (; i < chunks; i += 16) {+        __m128i dv = _mm_loadu_si128((const __m128i *)(d + i));+        __m128i sv = _mm_loadu_si128((const __m128i *)(s + i));+        _mm_storeu_si128((__m128i *)(d + i), _mm_max_epu8(dv, sv));+    }+    for (; i < m; i++) {+        uint8_t sv = s[i], dv = d[i];+        d[i] = sv > dv ? sv : dv;+    }+#else     for (int i = 0; i < m; i++) {         uint8_t sv = s[i], dv = d[i];         d[i] = sv > dv ? sv : dv;     }+#endif }  int hll_c_precision(const hll_sketch_t *sk) { return sk->p; }

@@ -4,6 +4,14 @@ #include <math.h> #include <float.h> +#if defined(__ARM_NEON) || defined(__ARM_NEON__)+#include <arm_neon.h>+#define KLL_HAS_NEON 1+#elif defined(__SSE2__)+#include <emmintrin.h>+#define KLL_HAS_SSE2 1+#endif+ static inline uint64_t rotl64(uint64_t x, int k) {     return (x << k) | (x >> (64 - k)); }@@ -254,8 +262,48 @@     for (int h = 0; h < sk->num_levels; h++) {         uint64_t weight = 1ULL << h;         int lo = sk->levels[h], hi = sk->levels[h + 1];-        for (int i = lo; i < hi; i++)-            count_below += weight * (sk->items[i] < value);+        int n = hi - lo;+        const double *items = sk->items + lo;+#if KLL_HAS_NEON+        {+            float64x2_t val_vec = vdupq_n_f64(value);+            uint64x2_t acc = vdupq_n_u64(0);+            int i = 0;+            int chunks = n & ~3;+            for (; i < chunks; i += 4) {+                /* vcltq_f64 returns all-ones (0xFFFF...) per lane where true.+                   All-ones as uint64 = UINT64_MAX; negate to get 1. */+                uint64x2_t cmp0 = vcltq_f64(vld1q_f64(items + i), val_vec);+                uint64x2_t cmp1 = vcltq_f64(vld1q_f64(items + i + 2), val_vec);+                acc = vsubq_u64(acc, cmp0);+                acc = vsubq_u64(acc, cmp1);+            }+            count_below += weight * (vgetq_lane_u64(acc, 0) + vgetq_lane_u64(acc, 1));+            for (; i < n; i++)+                count_below += weight * (items[i] < value);+        }+#elif KLL_HAS_SSE2+        {+            __m128d val_vec = _mm_set1_pd(value);+            uint64_t local = 0;+            int i = 0;+            int chunks = n & ~1;+            for (; i < chunks; i += 2) {+                __m128d data = _mm_loadu_pd(items + i);+                __m128d cmp = _mm_cmplt_pd(data, val_vec);+                /* Each lane is all-ones (-1 as int64) or zero */+                local -= (uint64_t)_mm_cvtsi128_si64(_mm_castpd_si128(cmp));+                local -= (uint64_t)_mm_cvtsi128_si64(+                    _mm_srli_si128(_mm_castpd_si128(cmp), 8));+            }+            count_below += weight * local;+            for (; i < n; i++)+                count_below += weight * (items[i] < value);+        }+#else+        for (int i = 0; i < n; i++)+            count_below += weight * (items[i] < value);+#endif     }     return (double)count_below / (double)sk->total_n; }

@@ -5,7 +5,7 @@ -- see: https://github.com/sol/hpack  name:           data-sketches-core-version:        0.2.0.1+version:        0.3.0.0 description:    Please see the README on GitHub at <https://github.com/iand675/datasketches-haskell#readme> homepage:       https://github.com/iand675/datasketches-haskell#readme bug-reports:    https://github.com/iand675/datasketches-haskell/issues@@ -27,14 +27,7 @@  library   exposed-modules:-      DataSketches.Core.Internal.URef-      DataSketches.Core.Snapshot-      DataSketches.Quantiles.RelativeErrorQuantile.Internal-      DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary-      DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor       DataSketches.Quantiles.RelativeErrorQuantile.Internal.Constants-      DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer-      DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch       DataSketches.Quantiles.RelativeErrorQuantile.Types       DataSketches.Quantiles.KLL.Internal       DataSketches.Frequencies.CountMin.Internal@@ -68,10 +61,8 @@       base >=4.7 && <5     , deepseq     , ghc-prim-    , mwc-random     , primitive     , vector-    , vector-algorithms   default-language: Haskell2010  test-suite data-sketches-core-test@@ -95,8 +86,6 @@     , data-sketches-core     , deepseq     , ghc-prim-    , mwc-random     , primitive     , vector-    , vector-algorithms   default-language: Haskell2010

@@ -1,73 +0,0 @@-{-# LANGUAGE MagicHash #-}-{-# LANGUAGE UnboxedTuples #-}-module DataSketches.Core.Internal.URef-  ( URef-  , IOURef-  , newURef-  , readURef-  , writeURef-  , modifyURef-  -- * Packed mutable byte arrays for multiple fields-  , MutableFields-  , newMutableFields-  , readField-  , writeField-  , modifyField-  ) where--import Control.Monad.Primitive-import Data.Primitive.ByteArray-import qualified Data.Vector.Unboxed.Mutable as MUVector-import Data.Vector.Unboxed (Unbox)-import Data.Primitive (Prim, readByteArray, writeByteArray, sizeOf)---- | An unboxed reference. Stores a single value in a 'MutableByteArray'.-newtype URef s a = URef (MUVector.MVector s a)--type IOURef = URef (PrimState IO)--newURef :: (PrimMonad m, Unbox a) => a -> m (URef (PrimState m) a)-newURef a = fmap URef (MUVector.replicate 1 a)-{-# INLINE newURef #-}--readURef :: (PrimMonad m, Unbox a) => URef (PrimState m) a -> m a-readURef (URef v) = MUVector.unsafeRead v 0-{-# INLINE readURef #-}--writeURef :: (PrimMonad m, Unbox a) => URef (PrimState m) a -> a -> m ()-writeURef (URef v) = MUVector.unsafeWrite v 0-{-# INLINE writeURef #-}--modifyURef :: (PrimMonad m, Unbox a) => URef (PrimState m) a -> (a -> a) -> m ()-modifyURef (URef v) f = do-  !x <- MUVector.unsafeRead v 0-  MUVector.unsafeWrite v 0 $! f x-{-# INLINE modifyURef #-}---- | A single 'MutableByteArray' that packs multiple typed fields at byte offsets.--- Use 'readField' and 'writeField' with the byte offset of each field.--- Callers are responsible for computing non-overlapping offsets from 'Data.Primitive.sizeOf'.-newtype MutableFields s = MutableFields (MutableByteArray s)---- | Allocate a packed mutable fields block of the given total size in bytes.-newMutableFields :: PrimMonad m => Int -> m (MutableFields (PrimState m))-newMutableFields size = MutableFields <$> newByteArray size-{-# INLINE newMutableFields #-}---- | Read a 'Prim' value at the given ELEMENT index (not byte offset).--- The index is in units of @sizeOf a@.-readField :: (PrimMonad m, Prim a) => MutableFields (PrimState m) -> Int -> m a-readField (MutableFields mba) = readByteArray mba-{-# INLINE readField #-}---- | Write a 'Prim' value at the given ELEMENT index.-writeField :: (PrimMonad m, Prim a) => MutableFields (PrimState m) -> Int -> a -> m ()-writeField (MutableFields mba) = writeByteArray mba-{-# INLINE writeField #-}---- | Strict read-modify-write at the given ELEMENT index.-modifyField :: (PrimMonad m, Prim a) => MutableFields (PrimState m) -> Int -> (a -> a) -> m ()-modifyField mf ix f = do-  !x <- readField mf ix-  writeField mf ix $! f x-{-# INLINE modifyField #-}

@@ -1,7 +0,0 @@-module DataSketches.Core.Snapshot where--import Control.Monad.Primitive--class TakeSnapshot a where-  type Snapshot a-  takeSnapshot :: PrimMonad m => a (PrimState m) -> m (Snapshot a)

@@ -1,184 +0,0 @@-{-# LANGUAGE DeriveGeneric #-}-module DataSketches.Quantiles.RelativeErrorQuantile.Internal where--import Data.Primitive (MutVar, readMutVar)-import Data.Word-import qualified Data.Vector as Vector-import GHC.Generics-import System.Random.MWC (Gen)--import DataSketches.Core.Snapshot-import DataSketches.Quantiles.RelativeErrorQuantile.Types-import DataSketches.Core.Internal.URef (MutableFields, readField, writeField, modifyField)-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary (ReqAuxiliary)-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor (ReqCompactor)-import Control.DeepSeq (NFData, rnf)-import Control.Monad.Primitive (PrimMonad (PrimState))-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor as Compactor-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer as DoubleBuffer-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary as Auxiliary-import Control.Exception (Exception)---- | Mutable sketch state packed into a single 'MutableByteArray'.------ All 6 scalar fields occupy 48 bytes (6 × 8-byte slots), fitting entirely--- within one 64-byte x86-64 cache line. On the insert hot path, totalN,--- retainedItems, maxNominalCapacitiesSize, minValue, maxValue, and sumValue--- are all read/written — a single cache line fetch covers all of them.------ Layout (element index → field):---   0: totalN                 (Int, representing Word64)---   1: minValue               (Double)---   2: maxValue               (Double)---   3: sumValue               (Double)---   4: retainedItems          (Int)---   5: maxNominalCapacitiesSize (Int)-data ReqSketch s = ReqSketch-  { k :: !Word32-  , rankAccuracySetting :: !RankAccuracy-  , criterion :: !Criterion-  , sketchRng :: {-# UNPACK #-} !(Gen s)-  , sketchFields :: {-# UNPACK #-} !(MutableFields s)-  , aux :: {-# UNPACK #-} !(MutVar s (Maybe ReqAuxiliary))-  , compactors :: {-# UNPACK #-} !(MutVar s (Vector.Vector (ReqCompactor s)))-  } deriving (Generic)--instance NFData (ReqSketch s) where-  rnf !_ = ()---- Field indices-fTotalN, fMinValue, fMaxValue, fSumValue, fRetainedItems, fMaxNomCapSize :: Int-fTotalN         = 0-fMinValue       = 1-fMaxValue       = 2-fSumValue       = 3-fRetainedItems  = 4-fMaxNomCapSize  = 5--sketchFieldBytes :: Int-sketchFieldBytes = 6 * 8---- Typed accessors--getTotalN :: PrimMonad m => ReqSketch (PrimState m) -> m Word64-getTotalN sk = do-  v <- readField (sketchFields sk) fTotalN-  pure $! fromIntegral (v :: Int)-{-# INLINE getTotalN #-}--setTotalN :: PrimMonad m => ReqSketch (PrimState m) -> Word64 -> m ()-setTotalN sk v = writeField (sketchFields sk) fTotalN (fromIntegral v :: Int)-{-# INLINE setTotalN #-}--modifyTotalN :: PrimMonad m => ReqSketch (PrimState m) -> (Word64 -> Word64) -> m ()-modifyTotalN sk f = do-  !v <- getTotalN sk-  setTotalN sk $! f v-{-# INLINE modifyTotalN #-}--getMinValue :: PrimMonad m => ReqSketch (PrimState m) -> m Double-getMinValue sk = readField (sketchFields sk) fMinValue-{-# INLINE getMinValue #-}--setMinValue :: PrimMonad m => ReqSketch (PrimState m) -> Double -> m ()-setMinValue sk = writeField (sketchFields sk) fMinValue-{-# INLINE setMinValue #-}--getMaxValue :: PrimMonad m => ReqSketch (PrimState m) -> m Double-getMaxValue sk = readField (sketchFields sk) fMaxValue-{-# INLINE getMaxValue #-}--setMaxValue :: PrimMonad m => ReqSketch (PrimState m) -> Double -> m ()-setMaxValue sk = writeField (sketchFields sk) fMaxValue-{-# INLINE setMaxValue #-}--getSumValue :: PrimMonad m => ReqSketch (PrimState m) -> m Double-getSumValue sk = readField (sketchFields sk) fSumValue-{-# INLINE getSumValue #-}--modifySumValue :: PrimMonad m => ReqSketch (PrimState m) -> (Double -> Double) -> m ()-modifySumValue sk f = modifyField (sketchFields sk) fSumValue f-{-# INLINE modifySumValue #-}--getRetainedItems :: PrimMonad m => ReqSketch (PrimState m) -> m Int-getRetainedItems sk = readField (sketchFields sk) fRetainedItems-{-# INLINE getRetainedItems #-}--setRetainedItems :: PrimMonad m => ReqSketch (PrimState m) -> Int -> m ()-setRetainedItems sk = writeField (sketchFields sk) fRetainedItems-{-# INLINE setRetainedItems #-}--modifyRetainedItems :: PrimMonad m => ReqSketch (PrimState m) -> (Int -> Int) -> m ()-modifyRetainedItems sk = modifyField (sketchFields sk) fRetainedItems-{-# INLINE modifyRetainedItems #-}--getMaxNomCapSize :: PrimMonad m => ReqSketch (PrimState m) -> m Int-getMaxNomCapSize sk = readField (sketchFields sk) fMaxNomCapSize-{-# INLINE getMaxNomCapSize #-}--setMaxNomCapSize :: PrimMonad m => ReqSketch (PrimState m) -> Int -> m ()-setMaxNomCapSize sk = writeField (sketchFields sk) fMaxNomCapSize-{-# INLINE setMaxNomCapSize #-}--modifyMaxNomCapSize :: PrimMonad m => ReqSketch (PrimState m) -> (Int -> Int) -> m ()-modifyMaxNomCapSize sk = modifyField (sketchFields sk) fMaxNomCapSize-{-# INLINE modifyMaxNomCapSize #-}--data ReqSketchSnapshot = ReqSketchSnapshot-    { snapshotRankAccuracySetting :: !RankAccuracy-    , snapshotCriterion :: !Criterion-    , snapshotTotalN :: !Word64-    , snapshotMinValue :: !Double-    , snapshotMaxValue :: !Double-    , snapshotRetainedItems :: !Int-    , snapshotMaxNominalCapacitiesSize :: !Int-    , snapshotCompactors :: !(Vector.Vector (Snapshot ReqCompactor))-    } deriving Show--instance TakeSnapshot ReqSketch where-  type Snapshot ReqSketch = ReqSketchSnapshot-  takeSnapshot sk = do-    tn <- getTotalN sk-    mn <- getMinValue sk-    mx <- getMaxValue sk-    ri <- getRetainedItems sk-    mc <- getMaxNomCapSize sk-    cs <- readMutVar (compactors sk) >>= mapM takeSnapshot-    pure $ ReqSketchSnapshot (rankAccuracySetting sk) (criterion sk) tn mn mx ri mc cs--getCompactors :: PrimMonad m => ReqSketch (PrimState m) -> m (Vector.Vector (ReqCompactor (PrimState m)))-getCompactors = readMutVar . compactors-{-# INLINE getCompactors #-}--computeTotalRetainedItems :: PrimMonad m => ReqSketch (PrimState m) -> m Int-computeTotalRetainedItems this = do-  cs <- getCompactors this-  Vector.foldM countBuffer 0 cs-  where-    countBuffer acc compactor = do-      buff <- Compactor.getBuffer compactor-      buffSize <- DoubleBuffer.getCount buff-      pure $! buffSize + acc--retainedItemCount :: PrimMonad m => ReqSketch (PrimState m) -> m Int-retainedItemCount = getRetainedItems-{-# INLINE retainedItemCount #-}--count :: PrimMonad m => ReqSketch (PrimState m) -> m Word64-count = getTotalN-{-# INLINE count #-}--mkAuxiliaryFromReqSketch :: PrimMonad m => ReqSketch (PrimState m) -> m ReqAuxiliary-mkAuxiliaryFromReqSketch this = do-  total <- count this-  ri <- retainedItemCount this-  cs <- getCompactors this-  Auxiliary.mkAuxiliary (rankAccuracySetting this) total ri cs--data CumulativeDistributionInvariants-  = CumulativeDistributionInvariantsSplitsAreEmpty-  | CumulativeDistributionInvariantsSplitsAreNotFinite-  | CumulativeDistributionInvariantsSplitsAreNotUniqueAndMontonicallyIncreasing-  deriving (Show, Eq)--instance Exception CumulativeDistributionInvariants

@@ -1,189 +0,0 @@-module DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary-  ( ReqAuxiliary(..)-  , MReqAuxiliary (..)-  , mkAuxiliary-  , getQuantile-  -- | Really extra private, just needed for tests-  , mergeSortIn-  ) where--import GHC.TypeLits-import Control.Monad (when)-import Control.Monad.Primitive-import Data.Bits (shiftL)-import Data.Word-import Data.Primitive.MutVar-import Data.Vector.Algorithms.Search-import qualified Data.Vector as Vector-import qualified Data.Vector.Unboxed.Mutable as MUVector-import DataSketches.Quantiles.RelativeErrorQuantile.Types-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor (ReqCompactor)-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor as Compactor-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer (DoubleBuffer)-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer as DoubleBuffer-import qualified Data.Vector.Unboxed as U-import Control.Monad.ST-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch (find)-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch as IS-import Debug.Trace-import qualified Data.Vector.Generic.Mutable as MG--data ReqAuxiliary = ReqAuxiliary-  { raWeightedItems :: {-# UNPACK #-} !(U.Vector (Double, Word64))-  , raHighRankAccuracy :: !RankAccuracy-  , raSize :: {-# UNPACK #-} !Word64-  }-  deriving (Show, Eq)--data MReqAuxiliary s = MReqAuxiliary-  { mraWeightedItems :: {-# UNPACK #-} !(MutVar s (MUVector.MVector s (Double, Word64)))-  , mraHighRankAccuracy :: !RankAccuracy-  , mraSize :: {-# UNPACK #-} !Word64-  }--mkAuxiliary :: (PrimMonad m, s ~ PrimState m) => RankAccuracy -> Word64 -> Int -> Vector.Vector (ReqCompactor s) -> m ReqAuxiliary-mkAuxiliary rankAccuracy totalN retainedItems compactors = do-  items <- newMutVar =<< MUVector.replicate retainedItems (0, 0)-  let this = MReqAuxiliary-        { mraWeightedItems = items-        , mraHighRankAccuracy = rankAccuracy-        , mraSize = totalN-        }-  Vector.foldM_ (mergeBuffers this) 0 compactors-  createCumulativeWeights this-  dedup this-  items' <- U.unsafeFreeze =<< readMutVar items-  pure ReqAuxiliary-    { raWeightedItems = items'-    , raHighRankAccuracy = rankAccuracy-    , raSize = totalN-    }-  where-    mergeBuffers this auxCount compactor = do-      buff <- Compactor.getBuffer compactor-      buffSize <-  DoubleBuffer.getCount buff-      let lgWeight = Compactor.getLgWeight compactor-          weight = 1 `shiftL` fromIntegral lgWeight-      mergeSortIn this buff weight auxCount-      pure $ auxCount + buffSize--getWeightedItems :: PrimMonad m => MReqAuxiliary (PrimState m) -> m (MUVector.MVector (PrimState m) (Double, Word64))-getWeightedItems = readMutVar . mraWeightedItems--getItems :: PrimMonad m => MReqAuxiliary (PrimState m) -> m (MUVector.MVector (PrimState m) Double)-getItems = fmap (fst . MUVector.unzip) . getWeightedItems--getWeights :: PrimMonad m => MReqAuxiliary (PrimState m) -> m (MUVector.MVector (PrimState m) Word64)-getWeights = fmap (snd . MUVector.unzip) . getWeightedItems--getQuantile :: ReqAuxiliary -> Double -> Criterion  -> Double-getQuantile this normalRank ltEq = fst (weightedItems U.! ix)-  where-    ix = if searchResult == U.length weightedItems-      then searchResult - 1-      else searchResult-    searchResult = runST $ do-      v <- U.unsafeThaw $ snd $ U.unzip weightedItems-      let search = case ltEq of-            (:<) -> find (IS.:>)-            (:<=) -> find (IS.:>=)-      search v 0 (weightsSize - 1) rank-    weightedItems = raWeightedItems this-    weightsSize = U.length weightedItems-    rank = floor (normalRank * fromIntegral (raSize this))--createCumulativeWeights :: PrimMonad m => MReqAuxiliary (PrimState m) -> m ()-createCumulativeWeights this = do-  weights <- getWeights this-  let size = MUVector.length weights-  let accumulateM i weight = do-        when (i > 0) $ do-          prevWeight <- MUVector.read weights (i - 1)-          MUVector.unsafeWrite weights i (weight + prevWeight)-  forI_ weights (\i -> MUVector.read weights i >>= \x -> accumulateM i x)-  lastWeight <- MUVector.read weights (size - 1)-  when (lastWeight /= mraSize this) $ do-    error "invariant violated: lastWeight does not equal raSize"-  where-    forI_ :: (Monad m, MG.MVector v a) => v (PrimState m) a -> (Int -> m b) -> m ()-    {-# INLINE forI_ #-}-    forI_ v f = loop 0-      where-        loop i -          | i >= n    = return ()-          | otherwise = f i >> loop (i + 1)-        n = MG.length v---dedup :: PrimMonad m => MReqAuxiliary (PrimState m) -> m ()-dedup this = do-  weightedItems <- getWeightedItems this-  let size = MUVector.length weightedItems-  weightedItemsB <- MUVector.replicate size (0, 0)-  bi <- doDedup weightedItems size weightedItemsB 0 0-  writeMutVar (mraWeightedItems this) $ MUVector.slice 0 bi weightedItemsB-  where-    doDedup weightedItems itemsSize weightedItemsB = go-      where -        go !i !bi-          | i >= itemsSize = pure bi-          | otherwise = do-            let j = i + 1-                hidup = j-                countDups !j !hidup = if j < itemsSize -                  then do-                    (itemI, _) <- MUVector.read weightedItems i-                    (itemJ, _) <- MUVector.read weightedItems j-                    if itemI == itemJ-                      then countDups (j + 1) j-                      else pure (j, hidup)-                  else pure (j, hidup)-            (j', hidup') <- countDups j hidup-            if j' - i == 1 -- no dups-              then do-                (item, weight) <- MUVector.read weightedItems i-                MUVector.unsafeWrite weightedItemsB bi (item, weight)-                go (i + 1) (bi + 1)-              else do-                (item, weight) <- MUVector.read weightedItems hidup'-                MUVector.unsafeWrite weightedItemsB bi (item, weight)-                go j' (bi + 1)--mergeSortIn :: PrimMonad m => MReqAuxiliary (PrimState m) -> DoubleBuffer (PrimState m) -> Word64 -> Int -> m ()-mergeSortIn this bufIn defaultWeight auxCount = do-  DoubleBuffer.sort bufIn-  weightedItems <- getWeightedItems this-  otherItems <- DoubleBuffer.getVector bufIn-  otherBuffSize <- DoubleBuffer.getCount bufIn-  otherBuffCapacity <- DoubleBuffer.getCapacity bufIn-  let totalSize = otherBuffSize + auxCount - 1-      height = case mraHighRankAccuracy this of-        HighRanksAreAccurate -> otherBuffCapacity - 1-        LowRanksAreAccurate -> otherBuffSize - 1-  merge totalSize weightedItems otherItems (auxCount - 1) (otherBuffSize - 1) height -  where-    merge totalSize weightedItems otherItems = go totalSize-      where-        go !k !i !j !h -          | k < 0 = pure ()-          | i >= 0 && j >= 0 = do-            (item, weight) <- MUVector.read weightedItems i-            otherItem <- MUVector.read otherItems h-            if item >= otherItem-               then do-                 MUVector.unsafeWrite weightedItems k (item, weight)-                 continue (i - 1) j h-               else do-                 MUVector.unsafeWrite weightedItems k (otherItem, defaultWeight)-                 continue i (j - 1) (h - 1)-          | i >= 0 = do-            MUVector.read weightedItems i >>= MUVector.write weightedItems k-            continue (i - 1) j h-          | j >= 0 = do-            otherItem <- MUVector.read otherItems h-            MUVector.unsafeWrite weightedItems k (otherItem, defaultWeight)-            continue i (j - 1) (h - 1)-          | otherwise = pure ()-          where-            continue = go (k - 1)-

@@ -1,285 +0,0 @@-module DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor-  ( ReqCompactor-  , mkReqCompactor-  , CompactorReturn (..)-  , compact-  , getBuffer-  , getCoin-  , getLgWeight-  , getNominalCapacity-  , getNumSections-  , merge-  , nearestEven-  ) where--import Data.Bits ((.&.), (.|.), complement, countTrailingZeros, shiftL, shiftR)-import Data.Primitive.MutVar-import Data.Word-import DataSketches.Quantiles.RelativeErrorQuantile.Types-import System.Random.MWC (Variate(uniform), Gen)-import Control.Exception (assert)-import Control.Monad (when)-import Control.Monad.Primitive-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Constants-import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer-import DataSketches.Core.Internal.URef (MutableFields, newMutableFields, readField, writeField, modifyField)-import DataSketches.Core.Snapshot--data CompactorReturn s = CompactorReturn-  { crDeltaRetItems :: {-# UNPACK #-} !Int-  , crDeltaNominalSize :: {-# UNPACK #-} !Int-  , crDoubleBuffer :: {-# UNPACK #-} !(DoubleBuffer s)-  }---- | Mutable compactor state packed into a single 'MutableByteArray'.------ All scalar fields are stored at 8-byte-aligned offsets in one contiguous--- allocation. On x86-64 (64-byte cache lines), the entire 40-byte block--- fits in a single cache line, so reading state + numSections + sectionSize--- during compact (the hot path) never causes a second cache miss.------ Layout (byte offset → field):---   0:  state          (Word64) — compaction counter, read/written every compact---   8:  sectionSizeFlt (Double) — current section size as float---  16:  sectionSize    (Int)    — current section size---  24:  numSections    (Int)    — current number of sections---  32:  lastFlip       (Int)    — last coin flip result (0 or 1)-data ReqCompactor s = ReqCompactor-  { rcRankAccuracy :: !RankAccuracy-  , rcLgWeight :: {-# UNPACK #-} !Word8-  , rcRng :: {-# UNPACK #-} !(Gen s)-  , rcFields :: {-# UNPACK #-} !(MutableFields s)-  , rcBuffer :: {-# UNPACK #-} !(MutVar s (DoubleBuffer s))-  }---- Field indices (element index into the Int/Double-sized slots of MutableByteArray)--- All stored as 8-byte values for uniform alignment.-fState, fSectionSizeFlt, fSectionSize, fNumSections, fLastFlip :: Int-fState          = 0  -- Word64, but stored via Int-sized readField/writeField-fSectionSizeFlt = 1  -- Double-fSectionSize    = 2  -- Int-fNumSections    = 3  -- Int-fLastFlip       = 4  -- Int (0 or 1)--fieldBytes :: Int-fieldBytes = 5 * 8  -- 40 bytes total, fits in one 64-byte cache line---- Typed accessors. These read/write the MutableByteArray at the right element--- index, relying on Prim instances for Int, Double, Word64.--getStateField :: PrimMonad m => ReqCompactor (PrimState m) -> m Word64-getStateField rc = do-  v <- readField (rcFields rc) fState-  pure $! fromIntegral (v :: Int)-{-# INLINE getStateField #-}--setStateField :: PrimMonad m => ReqCompactor (PrimState m) -> Word64 -> m ()-setStateField rc v = writeField (rcFields rc) fState (fromIntegral v :: Int)-{-# INLINE setStateField #-}--modifyStateField :: PrimMonad m => ReqCompactor (PrimState m) -> (Word64 -> Word64) -> m ()-modifyStateField rc f = do-  !v <- getStateField rc-  setStateField rc $! f v-{-# INLINE modifyStateField #-}--getSectionSizeFltField :: PrimMonad m => ReqCompactor (PrimState m) -> m Double-getSectionSizeFltField rc = readField (rcFields rc) fSectionSizeFlt-{-# INLINE getSectionSizeFltField #-}--setSectionSizeFltField :: PrimMonad m => ReqCompactor (PrimState m) -> Double -> m ()-setSectionSizeFltField rc = writeField (rcFields rc) fSectionSizeFlt-{-# INLINE setSectionSizeFltField #-}--getSectionSizeField :: PrimMonad m => ReqCompactor (PrimState m) -> m Int-getSectionSizeField rc = readField (rcFields rc) fSectionSize-{-# INLINE getSectionSizeField #-}--setSectionSizeField :: PrimMonad m => ReqCompactor (PrimState m) -> Int -> m ()-setSectionSizeField rc = writeField (rcFields rc) fSectionSize-{-# INLINE setSectionSizeField #-}--getNumSectionsField :: PrimMonad m => ReqCompactor (PrimState m) -> m Int-getNumSectionsField rc = readField (rcFields rc) fNumSections-{-# INLINE getNumSectionsField #-}--setNumSectionsField :: PrimMonad m => ReqCompactor (PrimState m) -> Int -> m ()-setNumSectionsField rc = writeField (rcFields rc) fNumSections-{-# INLINE setNumSectionsField #-}--getLastFlipField :: PrimMonad m => ReqCompactor (PrimState m) -> m Bool-getLastFlipField rc = (/= (0 :: Int)) <$> readField (rcFields rc) fLastFlip-{-# INLINE getLastFlipField #-}--setLastFlipField :: PrimMonad m => ReqCompactor (PrimState m) -> Bool -> m ()-setLastFlipField rc b = writeField (rcFields rc) fLastFlip (if b then 1 :: Int else 0)-{-# INLINE setLastFlipField #-}--data ReqCompactorSnapshot = ReqCompactorSnapshot-    { snapshotCompactorRankAccuracy :: !RankAccuracy-    , snapshotCompactorRankAccuracyState :: !Word64-    , snapshotCompactorLastFlip :: !Bool-    , snapshotCompactorSectionSizeFlt :: !Double-    , snapshotCompactorSectionSize :: !Word32-    , snapshotCompactorNumSections :: !Word8-    , snapshotCompactorBuffer :: !(Snapshot DoubleBuffer)-    } deriving (Show)--instance TakeSnapshot ReqCompactor where-  type Snapshot ReqCompactor = ReqCompactorSnapshot-  takeSnapshot rc = do-    st <- getStateField rc-    fl <- getLastFlipField rc-    szf <- getSectionSizeFltField rc-    sz <- getSectionSizeField rc-    ns <- getNumSectionsField rc-    buf <- readMutVar (rcBuffer rc) >>= takeSnapshot-    pure $ ReqCompactorSnapshot-      (rcRankAccuracy rc) st fl szf (fromIntegral sz) (fromIntegral ns) buf--mkReqCompactor-  :: PrimMonad m-  => Gen (PrimState m)-  -> Word8-  -> RankAccuracy-  -> Word32-  -> m (ReqCompactor (PrimState m))-mkReqCompactor g lgWeight rankAccuracy sectionSize = do-  let nominalCapacity = fromIntegral $ nomCapMulti * initNumberOfSections * sectionSize-  buff <- mkBuffer (nominalCapacity * 2) nominalCapacity (rankAccuracy == HighRanksAreAccurate)-  fields <- newMutableFields fieldBytes-  writeField fields fState (0 :: Int)-  writeField fields fSectionSizeFlt (fromIntegral sectionSize :: Double)-  writeField fields fSectionSize (fromIntegral sectionSize :: Int)-  writeField fields fNumSections (fromIntegral initNumberOfSections :: Int)-  writeField fields fLastFlip (0 :: Int)-  bufMv <- newMutVar buff-  pure $ ReqCompactor rankAccuracy lgWeight g fields bufMv--nomCapMult :: Int-nomCapMult = 2--compact :: PrimMonad m => ReqCompactor (PrimState m) -> m (CompactorReturn (PrimState m))-compact this = do-  startBuffSize <- getCount =<< getBuffer this-  startNominalCapacity <- getNominalCapacity this-  numSections <- getNumSectionsField this-  sectionSize <- getSectionSizeField this-  state <- getStateField this-  let trailingOnes = succ $ countTrailingZeros $ complement state-      sectionsToCompact = min trailingOnes numSections-  (compactionStart, compactionEnd) <- computeCompactionRange this sectionsToCompact-  assert (compactionEnd - compactionStart >= 2) $ do-    coin <- if state .&. 1 == 1-      then fmap not $ getLastFlipField this-      else flipCoin this-    setLastFlipField this coin-    buff <- getBuffer this-    promote <- getEvensOrOdds buff compactionStart compactionEnd coin-    trimCount buff $ startBuffSize - (compactionEnd - compactionStart)-    modifyStateField this (+ 1)-    ensureEnoughSections this-    endBuffSize <- getCount buff-    promoteBuffSize <- getCount promote-    endNominalCapacity <- getNominalCapacity this-    pure $ CompactorReturn-      { crDeltaRetItems = endBuffSize - startBuffSize + promoteBuffSize-      , crDeltaNominalSize = endNominalCapacity - startNominalCapacity-      , crDoubleBuffer = promote-      }--getLgWeight :: ReqCompactor s -> Word8-getLgWeight = rcLgWeight--getBuffer :: PrimMonad m => ReqCompactor (PrimState m) -> m (DoubleBuffer (PrimState m))-getBuffer = readMutVar . rcBuffer-{-# INLINE getBuffer #-}--flipCoin :: PrimMonad m => ReqCompactor (PrimState m) -> m Bool-flipCoin = uniform . rcRng-{-# INLINE flipCoin #-}--getCoin :: PrimMonad m => ReqCompactor (PrimState m) -> m Bool-getCoin = getLastFlipField--getNominalCapacity :: PrimMonad m => ReqCompactor (PrimState m) -> m Int-getNominalCapacity compactor = do-  numSections <- getNumSectionsField compactor-  sectionSize <- getSectionSizeField compactor-  pure $! nomCapMult * numSections * sectionSize-{-# INLINE getNominalCapacity #-}--getNumSections :: PrimMonad m => ReqCompactor (PrimState m) -> m Word8-getNumSections rc = fromIntegral <$> getNumSectionsField rc--merge-  :: (PrimMonad m, s ~ PrimState m)-  => ReqCompactor (PrimState m)-  -> ReqCompactor (PrimState m)-  -> m (ReqCompactor s)-merge this otherCompactor = assert (rcLgWeight this == rcLgWeight otherCompactor) $ do-  otherState <- getStateField otherCompactor-  modifyStateField this (.|. otherState)-  ensureMaxSections--  buff <- getBuffer this-  sort buff--  otherBuff <- getBuffer otherCompactor-  sort otherBuff--  otherBuffIsBigger <- (>) <$> getCount otherBuff <*> getCount buff-  if otherBuffIsBigger-     then do-        otherBuff' <- copyBuffer otherBuff-        mergeSortIn otherBuff' buff-        writeMutVar (rcBuffer this) otherBuff'-     else mergeSortIn buff otherBuff-  pure this-  where-    ensureMaxSections = do-      adjusted <- ensureEnoughSections this-      when adjusted ensureMaxSections--ensureEnoughSections-  :: PrimMonad m-  => ReqCompactor (PrimState m)-  -> m Bool-ensureEnoughSections compactor = do-  sectionSizeFlt <- getSectionSizeFltField compactor-  let szf = sectionSizeFlt / sqrt2-      ne = nearestEven szf-  state <- getStateField compactor-  numSections <- getNumSectionsField compactor-  sectionSize <- getSectionSizeField compactor-  if state >= (1 `shiftL` (numSections - 1))-     && sectionSize > minK-     && ne >= minK-     then do-       setSectionSizeFltField compactor szf-       setSectionSizeField compactor ne-       setNumSectionsField compactor (numSections `shiftL` 1)-       buf <- getBuffer compactor-       nomCapacity <- getNominalCapacity compactor-       ensureCapacity buf (2 * nomCapacity)-       pure True-     else pure False--computeCompactionRange-  :: PrimMonad m-  => ReqCompactor (PrimState m)-  -> Int-  -> m (Int, Int)-computeCompactionRange this secsToCompact = do-  buffSize <- getCount =<< getBuffer this-  nominalCapacity <- getNominalCapacity this-  numSections <- getNumSectionsField this-  sectionSize <- getSectionSizeField this-  let nonCompact = (nominalCapacity `div` 2) + (numSections - secsToCompact) * sectionSize-      nonCompact' = if (buffSize - nonCompact) .&. 1 == 1 then nonCompact - 1 else nonCompact-  pure $ case rcRankAccuracy this of-    HighRanksAreAccurate -> (0, buffSize - nonCompact')-    LowRanksAreAccurate -> (nonCompact', buffSize)--nearestEven :: Double -> Int-nearestEven x = round (x / 2) `shiftL` 1

@@ -1,18 +1,7 @@ module DataSketches.Quantiles.RelativeErrorQuantile.Internal.Constants where--- Constants-import Data.Word -sqrt2 :: Double -sqrt2 = sqrt 2- initNumberOfSections :: Num a => a initNumberOfSections = 3--minK :: Num a => a-minK = 4--nomCapMulti :: Num a => a-nomCapMulti = 2  relRseFactor :: Double relRseFactor = sqrt (0.0512 / fromIntegral initNumberOfSections)

@@ -1,340 +0,0 @@-{-# LANGUAGE MagicHash #-}-{-# LANGUAGE UnboxedTuples #-}-module DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer-  ( DoubleBuffer-  , Capacity-  , GrowthIncrement-  , SpaceAtBottom-  , DoubleIsNonFiniteException(..)-  , mkBuffer-  , copyBuffer-  , append-  , ensureCapacity-  , getCountWithCriterion-  , getEvensOrOdds-  , (!) -- getItem-  , growthIncrement-  , spaceAtBottom-  , getCapacity-  , getCount-  , getSpace-  , getVector-  , isEmpty-  , isSorted-  , sort-  , mergeSortIn-  , trimCount-  ) where--import DataSketches.Quantiles.RelativeErrorQuantile.Types-    ( Criterion )-import Control.Monad ( unless, when )-import Control.Monad.Primitive ( PrimMonad(PrimState) )-import Data.Primitive.MutVar-    ( newMutVar, readMutVar, writeMutVar, MutVar )-import qualified Data.Vector.Unboxed as UVector-import qualified Data.Vector.Unboxed.Mutable as MUVector-import DataSketches.Core.Internal.URef-    ( URef, newURef, readURef, writeURef, modifyURef-    , MutableFields, newMutableFields, readField, writeField, modifyField )-import Data.Vector.Algorithms.Intro (sortByBounds)-import GHC.Stack ( HasCallStack )-import System.IO.Unsafe ()-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch as IS-import Control.Exception ( Exception, throw )-import DataSketches.Core.Snapshot ( TakeSnapshot(..) )---- | A special buffer of floats specifically designed to support the ReqCompactor class.------ Mutable scalars (count, sorted flag) are packed into a single--- 'MutableFields' to avoid per-field MutableByteArray# overhead.--- Field layout: index 0 = count (Int), index 1 = sorted (Int, 0 or 1).-data DoubleBuffer s = DoubleBuffer-  { vec :: {-# UNPACK #-} !(MutVar s (MUVector.MVector s Double))-  , dbFields :: {-# UNPACK #-} !(MutableFields s)-  , growthIncrement :: {-# UNPACK #-} !Int-  , spaceAtBottom :: !Bool-  }--dbCountIx, dbSortedIx :: Int-dbCountIx = 0-dbSortedIx = 1--data DoubleBufferSnapshot = DoubleBufferSnapshot-    { dbSnapshotVec :: UVector.Vector Double-    , dbSnapshotCount :: !Int-    , dbSnapshotSorted :: !Bool-    , dbSnapshotGrowthIncrement :: !Int-    , dbSnapshotSpaceAtBottom :: !Bool-    } deriving (Show)--instance TakeSnapshot DoubleBuffer where-  type Snapshot DoubleBuffer = DoubleBufferSnapshot--  takeSnapshot DoubleBuffer{..} = do-    v <- readMutVar vec >>= UVector.freeze-    cnt <- readField dbFields dbCountIx-    srt <- readField dbFields dbSortedIx-    pure $ DoubleBufferSnapshot v cnt (srt /= (0 :: Int)) growthIncrement spaceAtBottom--type Capacity = Int-type GrowthIncrement = Int-type SpaceAtBottom = Bool---- | Constructs an new empty FloatBuffer with an initial capacity specified by--- the <code>capacity</code> argument.-mkBuffer :: PrimMonad m => Capacity -> GrowthIncrement -> SpaceAtBottom -> m (DoubleBuffer (PrimState m))-mkBuffer capacity_ growthIncrement spaceAtBottom = do-  vec <- newMutVar =<< MUVector.new capacity_-  -- Pack count and sorted into 2 Int-sized slots-  dbFields <- newMutableFields (2 * 8)-  writeField dbFields dbCountIx (0 :: Int)-  writeField dbFields dbSortedIx (1 :: Int)-  pure $ DoubleBuffer{..}--copyBuffer :: PrimMonad m => DoubleBuffer (PrimState m) -> m (DoubleBuffer (PrimState m))-copyBuffer buf@DoubleBuffer{..} = do-  vec <- newMutVar =<< MUVector.clone =<< getVector buf-  dbFields <- newMutableFields (2 * 8)-  cnt <- getCount buf-  srt <- isSorted buf-  writeField dbFields dbCountIx cnt-  writeField dbFields dbSortedIx (if srt then 1 :: Int else 0)-  pure $ DoubleBuffer {..}---- | Appends the given item to the active array and increments the active count.--- This will expand the array if necessary.-append :: PrimMonad m => DoubleBuffer (PrimState m) -> Double -> m ()-append buf@DoubleBuffer{..} x = do-  ensureSpace buf 1-  count_ <- getCount buf-  index <- if spaceAtBottom-    then do-      capacity_ <- getCapacity buf-      pure $! capacity_ - count_ - 1-    else pure count_-  writeField dbFields dbCountIx (count_ + 1)-  v <- getVector buf-  MUVector.unsafeWrite v index x-  writeField dbFields dbSortedIx (0 :: Int)-{-# INLINE append #-}---- | Ensures that the capacity of this FloatBuffer is at least newCapacity.--- If newCapacity &lt; capacity(), no action is taken.-ensureSpace :: PrimMonad m => DoubleBuffer (PrimState m) -> Int -> m ()-ensureSpace buf@DoubleBuffer{..} space = do-  count_ <- getCount buf-  capacity_ <- getCapacity buf-  when (count_ + space > capacity_) $ do-    let newCap = count_ + space + growthIncrement-    ensureCapacity buf newCap--getVector :: PrimMonad m => DoubleBuffer (PrimState m) -> m (MUVector.MVector (PrimState m) Double)-getVector = readMutVar . vec-{-# INLINE getVector #-}--getCapacity :: PrimMonad m => DoubleBuffer (PrimState m) -> m Int-getCapacity buf = MUVector.length <$> getVector buf-{-# INLINE getCapacity #-}--ensureCapacity :: PrimMonad m => DoubleBuffer (PrimState m) -> Int -> m ()-ensureCapacity buf@DoubleBuffer{..} newCapacity = do-  capacity_ <- getCapacity buf-  when (newCapacity > capacity_) $ do-    count_ <- getCount buf-    (srcPos, destPos) <- if spaceAtBottom-      then do-        pure (capacity_ - count_, newCapacity - count_)-      else pure (0, 0)-    oldVec <- getVector buf-    newVec <- MUVector.new newCapacity-    MUVector.unsafeCopy-      (MUVector.slice destPos count_ newVec)-      (MUVector.slice srcPos count_ oldVec)-    writeMutVar vec newVec-{-# SCC ensureCapacity #-}--newtype DoubleIsNonFiniteException = DoubleIsNonFiniteException Double-  deriving (Show, Eq)--instance Exception DoubleIsNonFiniteException--getCountWithCriterion :: PrimMonad m => DoubleBuffer (PrimState m) -> Double -> Criterion -> m Int-getCountWithCriterion buf@DoubleBuffer{..} value criterion = do-  when (isNaN value || isInfinite value) $ throw $ DoubleIsNonFiniteException value-  sort buf-  count_ <- getCount buf-  vec <- getVector buf-  (low, high) <- if spaceAtBottom-    then do-      capacity_ <- getCapacity buf-      pure (capacity_ - count_, capacity_ - 1)-    else pure (0, count_ - 1)--  ix <- IS.find criterion vec low high value-  pure $! if ix == MUVector.length vec-    then 0-    else ix - low + 1---- data EvensOrOdds = Evens | Odds--getEvensOrOdds :: PrimMonad m => DoubleBuffer (PrimState m) -> Int -> Int -> Bool -> m (DoubleBuffer (PrimState m))-getEvensOrOdds buf@DoubleBuffer{..} startOffset endOffset odds = do-  (start, end) <- if spaceAtBottom-    then do-      basis <- (-) <$> getCapacity buf <*> getCount buf-      pure (basis + startOffset, basis + endOffset)-    else pure (startOffset, endOffset)-  sort buf-  let range = endOffset - startOffset-  vec <- getVector buf-  out <- MUVector.new (range `div` 2)-  go vec out start 0-  where-    odd = if odds then 1 else 0-    go vec !out !i !j = if j < MUVector.length out-      then do-        MUVector.unsafeWrite out j =<< MUVector.unsafeRead vec (i + odd)-        go vec out (i + 2) (j + 1)-      else do-        dbFields <- newMutableFields (2 * 8)-        writeField dbFields dbCountIx (MUVector.length out)-        writeField dbFields dbSortedIx (1 :: Int)-        vec <- newMutVar out-        pure DoubleBuffer-          { vec = vec-          , dbFields = dbFields-          , growthIncrement = 0-          , spaceAtBottom = spaceAtBottom-          }-{-# SCC getEvensOrOdds #-}---(!) :: PrimMonad m => DoubleBuffer (PrimState m) -> Int -> m Double-(!) buf offset = do-  index <- if spaceAtBottom buf-    then do-      capacity_ <- getCapacity buf-      count_ <- getCount buf-      pure $! capacity_ - count_ + offset-    else pure offset-  vec <- getVector buf-  MUVector.read vec index--getCount :: PrimMonad m => DoubleBuffer (PrimState m) -> m Int-getCount DoubleBuffer{..} = readField dbFields dbCountIx-{-# INLINE getCount #-}--getSpace :: PrimMonad m => DoubleBuffer (PrimState m) -> m Int-getSpace buf@DoubleBuffer{..} = do-  cap <- getCapacity buf-  cnt <- getCount buf-  pure $! cap - cnt-{-# INLINE getSpace #-}--isEmpty :: PrimMonad m => DoubleBuffer (PrimState m) -> m Bool-isEmpty buf = (== 0) <$> getCount buf-{-# INLINE isEmpty #-}--isSorted :: PrimMonad m => DoubleBuffer (PrimState m) -> m Bool-isSorted DoubleBuffer{..} = (/= (0 :: Int)) <$> readField dbFields dbSortedIx-{-# INLINE isSorted #-}---- | Sorts the active region-sort :: PrimMonad m => DoubleBuffer (PrimState m) -> m ()-sort buf@DoubleBuffer{..} = do-  sorted_ <- isSorted buf-  unless sorted_ $ do-    capacity_ <- getCapacity buf-    count_ <- getCount buf-    let (start, end) = if spaceAtBottom-          then (capacity_ - count_, capacity_)-          else (0, count_)-    v <- getVector buf-    sortByBounds compare v start end-    writeField dbFields dbSortedIx (1 :: Int)-{-# INLINE sort #-}---- | Merges the incoming sorted buffer into this sorted buffer.-mergeSortIn :: (PrimMonad m, HasCallStack) => DoubleBuffer (PrimState m) -> DoubleBuffer (PrimState m) -> m ()-mergeSortIn this bufIn = do-  sort this-  sort bufIn--  thatBuf <- getVector bufIn-  bufInLen <- getCount bufIn--  ensureSpace this bufInLen-  count_ <- getCount this-  let totalLength = count_ + bufInLen--  thisBuf <- getVector this--  if spaceAtBottom this-    then do -- scan up, insert at bottom-      capacity_ <- getCapacity this-      bufInCapacity_ <- getCapacity bufIn-      inSs <- takeSnapshot bufIn-      let i = capacity_ - count_-      let j = bufInCapacity_ - bufInLen-      let targetStart = capacity_ - totalLength-      let k = targetStart-      mergeUpwards thisBuf thatBuf capacity_ bufInCapacity_ i j k-    else do -- scan down, insert at top-      let i = count_ - 1-      let j = bufInLen - 1-      let k = totalLength-      mergeDownwards thisBuf thatBuf i j (k - 1)--  modifyField (dbFields this) dbCountIx (+ bufInLen)-  writeField (dbFields this) dbSortedIx (1 :: Int)-  pure ()-  where-    mergeUpwards thisBuf thatBuf capacity_ bufInCapacity_ = go-      where-        go !i !j !k-          -- for loop ended-          | k >= capacity_ = pure ()-          -- both valid-          | i < capacity_ && j < bufInCapacity_ = do-            iVal <- MUVector.read thisBuf i-            jVal <- MUVector.read thatBuf j-            if iVal <= jVal-              then MUVector.unsafeWrite thisBuf k iVal >> go (i + 1) j (k + 1)-              else MUVector.unsafeWrite thisBuf k jVal >> go i (j + 1) (k + 1)-          -- i is valid-          | i < capacity_ = do-            MUVector.unsafeWrite thisBuf k =<< MUVector.read thisBuf i-            go (i + 1) j (k + 1)-          -- j is valid-          | j < bufInCapacity_ = do-            MUVector.unsafeWrite thisBuf k =<< MUVector.read thatBuf j-            go i (j + 1) (k + 1)-          -- neither is valid, break;-          | otherwise = pure ()-    mergeDownwards thisBuf thatBuf !i !j !k-      -- for loop ended-      | k < 0 = pure ()-      -- both valid-      | i >= 0 && j >= 0 = do-        iVal <- MUVector.read thisBuf i-        jVal <- MUVector.read thatBuf j-        if iVal >= jVal-          then do-            MUVector.unsafeWrite thisBuf k iVal >> continue (i - 1) j (k - 1)-          else do-            MUVector.unsafeWrite thisBuf k jVal >> continue i (j - 1) (k - 1)-      | i >= 0 = do-        MUVector.unsafeWrite thisBuf k =<< MUVector.read thisBuf i-        continue (i - 1) j (k - 1)-      | j >= 0 = do-        MUVector.unsafeWrite thisBuf k =<< MUVector.read thatBuf j-        continue i (j - 1) (k - 1)-      -- neither is valid, break;-      | otherwise = pure ()-      where-        continue = mergeDownwards thisBuf thatBuf-{-# SCC mergeSortIn #-}--trimCount :: PrimMonad m => DoubleBuffer (PrimState m) -> Int -> m ()-trimCount DoubleBuffer{..} newCount = modifyField dbFields dbCountIx (\oldCount -> if newCount < oldCount then newCount else oldCount)

@@ -1,121 +0,0 @@-module DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch where--import Control.Monad.Primitive-import Data.Vector.Generic.Mutable (MVector)-import qualified Data.Vector.Generic.Mutable as MV--data (:<) = (:<)-data (:<=) = (:<=)-data (:>) = (:>)-data (:>=) = (:>=)--- JavaDoc copypasta--- --- This provides efficient, unique and unambiguous binary searching for inequality comparison criteria---  for ordered arrays of values that may include duplicate values. The inequality criteria include---  <, >, ==, >=, <=. All the inequality criteria use the same search algorithm.---  (Although == is not an inequality, it is included for convenience.)----  In order to make the searching unique and unambiguous, we modified the traditional binary---  search algorithm to search for adjacent pairs of values <i>{A, B}</i> in the values array---  instead of just a single value, where <i>A</i> and <i>B</i> are the array indicies of two---  adjacent values in the array. For all the search criteria, if the algorithm reaches the ends of---  the search range, the algorithm calls the <i>resolve()</i> method to determine what to---   return to the caller. If the key value cannot be resolved, it returns a -1 to the caller.----  Given an array of values <i>arr[]</i> and the search key value <i>v</i>, the algorithms for---  the searching criteria are as follows:</p>--- ---  <li><b>LT:</b> Find the highest ranked adjacent pair <i>{A, B}</i> such that:<br>---  <i>arr[A] &lt; v &le; arr[B]</i>. The normal return is the index <i>A</i>.---  </li>---  <li><b>LE:</b>  Find the highest ranked adjacent pair <i>{A, B}</i> such that:<br>---  <i>arr[A] &le; v &lt; arr[B]</i>. The normal return is the index <i>A</i>.---  </li>---  <li><b>EQ:</b>  Find the adjacent pair <i>{A, B}</i> such that:<br>---  <i>arr[A] &le; v &le; arr[B]</i>. The normal return is the index <i>A</i> or <i>B</i> whichever---  equals <i>v</i>, otherwise it returns -1.---  </li>---  <li><b>GE:</b>  Find the lowest ranked adjacent pair <i>{A, B}</i> such that:<br>---  <i>arr[A] &lt; v &le; arr[B]</i>. The normal return is the index <i>B</i>.---  </li>---  <li><b>GT:</b>  Find the lowest ranked adjacent pair <i>{A, B}</i> such that:<br>---  <i>arr[A] &le; v &lt; arr[B]</i>. The normal return is the index <i>B</i>.---  </li>---  </ul>-class InequalitySearch s where-  inequalityCompare :: Ord a-    => s-    -> a -    -- ^ V-    -> a-    -- ^ A-    -> a-    -- ^ B-    -> Ordering-    -- ^ 'GT' means we must search higher in the array, 'LT' means we must-    -- search lower in the array, or `EQ`, which means we have found -    -- the correct bounding pair.-  getIndex -    :: (PrimMonad m, MVector v a, Ord a) -    => s -    -> v (PrimState m) a -    -> Int -    -> Int -    -> a -    -> m Int-  resolve -    :: s -    -> Int -- Vector length-    -> (Int, Int) -    -- ^ Final low index, high index (lo, hi)-    -> (Int, Int) -    -- ^ Initial search region (low, high)-    -> Int-    -- ^ A thing--instance InequalitySearch (:<) where-  inequalityCompare _ v a b -    | v <= a = LT-    | b < v = GT-    | otherwise = EQ-  getIndex _ _ a _ _ = pure a-  resolve _ vl (lo, hi) (low, high) = if lo >= high then high else vl--instance InequalitySearch (:<=) where-  inequalityCompare _ v a b-    | v < a = LT-    | b <= v = GT-    | otherwise = EQ-  getIndex _ _ a _ _ = pure a-  resolve _ vl (lo, hi) (low, high) = if lo >= high then high else vl --instance InequalitySearch (:>) where-  inequalityCompare _ v a b-    | v < a = LT-    | b <= v = GT-    | otherwise = EQ-  getIndex _ _ _ b _ = pure b-  resolve _ vl (lo, hi) (low, high) = if hi <= low then low else vl--instance InequalitySearch (:>=) where-  inequalityCompare _ v a b-    | v <= a = LT-    | b < v = GT-    | otherwise = EQ-  getIndex _ _ _ b _ = pure b-  resolve _ vl (lo, hi) (low, high) = if hi <= low then low else vl--find :: (InequalitySearch s, PrimMonad m, MVector v a, Ord a) => s -> v (PrimState m) a -> Int -> Int -> a -> m Int-find strat v low high x = go low (high - 1)-  where-    go lo hi -      | lo <= hi && lo < high = do-          let mid = lo + ((hi - lo) `div` 2)-          midV <- MV.read v mid-          midV' <- MV.read v (mid + 1)-          case inequalityCompare strat x midV midV' of-            LT -> go lo (mid - 1)-            EQ -> getIndex strat v mid (mid + 1) x-            GT -> go (mid + 1) hi-      | otherwise = pure $! resolve strat (MV.length v) (lo, hi) (low, high)-{-# INLINE find #-}

@@ -1,25 +1,18 @@-module DataSketches.Quantiles.RelativeErrorQuantile.Types where-import Control.Monad.Primitive-import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch as IS--data Criterion = (:<) | (:<=)-  deriving (Show, Eq)--instance IS.InequalitySearch Criterion where-  inequalityCompare c = case c of-    (:<) -> IS.inequalityCompare (IS.:<)-    (:<=) -> IS.inequalityCompare(IS.:<=)-  resolve c = case c of-    (:<) -> IS.resolve (IS.:<)-    (:<=) -> IS.resolve (IS.:<=)-  getIndex c = case c of-    (:<) -> IS.getIndex (IS.:<)-    (:<=) -> IS.getIndex (IS.:<=)+module DataSketches.Quantiles.RelativeErrorQuantile.Types+  ( RankAccuracy(..)+  , DoubleIsNonFiniteException(..)+  ) where +import Control.Exception (Exception) -data RankAccuracy -  = HighRanksAreAccurate +data RankAccuracy+  = HighRanksAreAccurate   -- ^ High ranks are prioritized for better accuracy.   | LowRanksAreAccurate   -- ^ Low ranks are prioritized for better accuracy   deriving (Show, Eq)++newtype DoubleIsNonFiniteException = DoubleIsNonFiniteException Double+  deriving (Show, Eq)++instance Exception DoubleIsNonFiniteException