diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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 diff --git a/cbits/countmin.c b/cbits/countmin.c --- a/cbits/countmin.c +++ b/cbits/countmin.c @@ -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; } diff --git a/cbits/hll.c b/cbits/hll.c --- a/cbits/hll.c +++ b/cbits/hll.c @@ -3,6 +3,14 @@ #include #include +#if defined(__ARM_NEON) || defined(__ARM_NEON__) +#include +#define HLL_HAS_NEON 1 +#elif defined(__SSE2__) +#include +#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 / clzp; 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; } diff --git a/cbits/kll.c b/cbits/kll.c --- a/cbits/kll.c +++ b/cbits/kll.c @@ -4,6 +4,14 @@ #include #include +#if defined(__ARM_NEON) || defined(__ARM_NEON__) +#include +#define KLL_HAS_NEON 1 +#elif defined(__SSE2__) +#include +#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; } diff --git a/data-sketches-core.cabal b/data-sketches-core.cabal --- a/data-sketches-core.cabal +++ b/data-sketches-core.cabal @@ -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 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 diff --git a/src/DataSketches/Core/Internal/URef.hs b/src/DataSketches/Core/Internal/URef.hs deleted file mode 100644 --- a/src/DataSketches/Core/Internal/URef.hs +++ /dev/null @@ -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 #-} diff --git a/src/DataSketches/Core/Snapshot.hs b/src/DataSketches/Core/Snapshot.hs deleted file mode 100644 --- a/src/DataSketches/Core/Snapshot.hs +++ /dev/null @@ -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) diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal.hs deleted file mode 100644 --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal.hs +++ /dev/null @@ -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 diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs deleted file mode 100644 --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs +++ /dev/null @@ -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) - diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs deleted file mode 100644 --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs +++ /dev/null @@ -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 diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs @@ -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) diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs deleted file mode 100644 --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs +++ /dev/null @@ -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 capacity 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 < 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) diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs deleted file mode 100644 --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs +++ /dev/null @@ -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 {A, B} in the values array --- instead of just a single value, where A and B 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 resolve() 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 arr[] and the search key value v, the algorithms for --- the searching criteria are as follows:

--- ---
  • LT: Find the highest ranked adjacent pair {A, B} such that:
    --- arr[A] < v ≤ arr[B]. The normal return is the index A. ---
    • ---
  • LE: Find the highest ranked adjacent pair {A, B} such that:
    --- arr[A] ≤ v < arr[B]. The normal return is the index A. ---
    • ---
  • EQ: Find the adjacent pair {A, B} such that:
    --- arr[A] ≤ v ≤ arr[B]. The normal return is the index A or B whichever --- equals v, otherwise it returns -1. ---
    • ---
  • GE: Find the lowest ranked adjacent pair {A, B} such that:
    --- arr[A] < v ≤ arr[B]. The normal return is the index B. ---
    • ---
  • GT: Find the lowest ranked adjacent pair {A, B} such that:
    --- arr[A] ≤ v < arr[B]. The normal return is the index B. ---
    • --- -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 #-} diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs --- a/src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs @@ -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