diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,3 @@ +# Changelog for data-sketches-core + +## Unreleased changes diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Author name here (c) 2021 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Author name here nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,1 @@ +# data-sketches-core diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/data-sketches-core.cabal b/data-sketches-core.cabal new file mode 100644 --- /dev/null +++ b/data-sketches-core.cabal @@ -0,0 +1,84 @@ +cabal-version: 1.12 + +-- This file has been generated from package.yaml by hpack version 0.34.4. +-- +-- see: https://github.com/sol/hpack + +name: data-sketches-core +version: 0.1.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 +author: Ian Duncan +maintainer: ian@iankduncan.com +copyright: 2021 Ian Duncan, Rob Bassi, Mercury Technologies +license: BSD3 +license-file: LICENSE +build-type: Simple +extra-source-files: + README.md + ChangeLog.md + +source-repository head + type: git + location: https://github.com/iand675/datasketches-haskell + +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 + other-modules: + Paths_data_sketches_core + hs-source-dirs: + src + default-extensions: + BangPatterns + FlexibleInstances + RecordWildCards + ScopedTypeVariables + StandaloneDeriving + TypeFamilies + TypeOperators + build-depends: + base >=4.7 && <5 + , deepseq + , ghc-prim + , mwc-random + , primitive + , vector + , vector-algorithms + default-language: Haskell2010 + +test-suite data-sketches-core-test + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: + Paths_data_sketches_core + hs-source-dirs: + test + default-extensions: + BangPatterns + FlexibleInstances + RecordWildCards + ScopedTypeVariables + StandaloneDeriving + TypeFamilies + TypeOperators + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.7 && <5 + , 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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Core/Internal/URef.hs @@ -0,0 +1,49 @@ +module DataSketches.Core.Internal.URef where + +import Control.Monad.Primitive +import qualified Data.Vector.Unboxed.Mutable as MUVector +import Data.Vector.Unboxed (Unbox) + +-- | An unboxed reference. This works like an 'IORef', but the data is +-- stored in a bytearray instead of a heap object, avoiding +-- significant allocation overhead in some cases. For a concrete +-- example, see this Stack Overflow question: +-- . +-- +-- The first parameter is the state token type, the same as would be +-- used for the 'ST' monad. If you're using an 'IO'-based monad, you +-- can use the convenience 'IOURef' type synonym instead. +-- +-- @since 0.0.2.0 +newtype URef s a = URef (MUVector.MVector s a) + +-- | Helpful type synonym for using a 'URef' from an 'IO'-based stack. +-- +-- @since 0.0.2.0 +type IOURef = URef (PrimState IO) + +-- | Create a new 'URef' +-- +-- @since 0.0.2.0 +newURef :: (PrimMonad m, Unbox a) => a -> m (URef (PrimState m) a) +newURef a = fmap URef (MUVector.replicate 1 a) + +-- | Read the value in a 'URef' +-- +-- @since 0.0.2.0 +readURef :: (PrimMonad m, Unbox a) => URef (PrimState m) a -> m a +readURef (URef v) = MUVector.read v 0 + +-- | Write a value into a 'URef'. Note that this action is strict, and +-- will force evalution of the value. +-- +-- @since 0.0.2.0 +writeURef :: (PrimMonad m, Unbox a) => URef (PrimState m) a -> a -> m () +writeURef (URef v) = MUVector.unsafeWrite v 0 + +-- | Modify a value in a 'URef'. Note that this action is strict, and +-- will force evaluation of the result value. +-- +-- @since 0.0.2.0 +modifyURef :: (PrimMonad m, Unbox a) => URef (PrimState m) a -> (a -> a) -> m () +modifyURef u f = readURef u >>= writeURef u . f diff --git a/src/DataSketches/Core/Snapshot.hs b/src/DataSketches/Core/Snapshot.hs new file mode 100644 --- /dev/null +++ b/src/DataSketches/Core/Snapshot.hs @@ -0,0 +1,7 @@ +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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal.hs @@ -0,0 +1,130 @@ +{-# 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 (URef, readURef) +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) + + +{- | +This Relative Error Quantiles Sketch is the Haskell implementation based on the paper +"Relative Error Streaming Quantiles", https://arxiv.org/abs/2004.01668, and loosely derived from +a Python prototype written by Pavel Vesely, ported from the Java equivalent. + +This implementation differs from the algorithm described in the paper in the following: + +The algorithm requires no upper bound on the stream length. +Instead, each relative-compactor counts the number of compaction operations performed +so far (via variable state). Initially, the relative-compactor starts with INIT_NUMBER_OF_SECTIONS. +Each time the number of compactions (variable state) exceeds 2^{numSections - 1}, we double +numSections. Note that after merging the sketch with another one variable state may not correspond +to the number of compactions performed at a particular level, however, since the state variable +never exceeds the number of compactions, the guarantees of the sketch remain valid. + +The size of each section (variable k and sectionSize in the code and parameter k in +the paper) is initialized with a value set by the user via variable k. +When the number of sections doubles, we decrease sectionSize by a factor of sqrt(2). +This is applied at each level separately. Thus, when we double the number of sections, the +nominal compactor size increases by a factor of approx. sqrt(2) (+/- rounding). + +The merge operation here does not perform "special compactions", which are used in the paper +to allow for a tight mathematical analysis of the sketch. + +This implementation provides a number of capabilities not discussed in the paper or provided +in the Python prototype. + +The Python prototype only implemented high accuracy for low ranks. This implementation +provides the user with the ability to choose either high rank accuracy or low rank accuracy at +the time of sketch construction. + +- The Python prototype only implemented a comparison criterion of "<". This implementation +allows the user to switch back and forth between the "<=" criterion and the "<=" criterion. +-} +data ReqSketch s = ReqSketch + { k :: !Word32 + , rankAccuracySetting :: !RankAccuracy + , criterion :: !Criterion + , sketchRng :: {-# UNPACK #-} !(Gen s) + , totalN :: {-# UNPACK #-} !(URef s Word64) + , minValue :: {-# UNPACK #-} !(URef s Double) + , maxValue :: {-# UNPACK #-} !(URef s Double) + , sumValue :: {-# UNPACK #-} !(URef s Double) + , retainedItems :: {-# UNPACK #-} !(URef s Int) + , maxNominalCapacitiesSize :: {-# UNPACK #-} !(URef s Int) + , aux :: {-# UNPACK #-} !(MutVar s (Maybe ReqAuxiliary)) + , compactors :: {-# UNPACK #-} !(MutVar s (Vector.Vector (ReqCompactor s))) + } deriving (Generic) + +instance NFData (ReqSketch s) where + rnf !rs = () + +data ReqSketchSnapshot = ReqSketchSnapshot + { snapshotRankAccuracySetting :: !RankAccuracy + , snapshotCriterion :: !Criterion + , snapshotTotalN :: !Word64 + , snapshotMinValue :: !Double + , snapshotMaxValue :: !Double + , snapshotRetainedItems :: !Int + , snapshotMaxNominalCapacitiesSize :: !Int + -- , aux :: !(MutVar s (Maybe ())) + , snapshotCompactors :: !(Vector.Vector (Snapshot ReqCompactor)) + } deriving Show + +instance TakeSnapshot ReqSketch where + type Snapshot ReqSketch = ReqSketchSnapshot + takeSnapshot ReqSketch{..} = ReqSketchSnapshot rankAccuracySetting criterion + <$> readURef totalN + <*> readURef minValue + <*> readURef maxValue + <*> readURef retainedItems + <*> readURef maxNominalCapacitiesSize + <*> (readMutVar compactors >>= mapM takeSnapshot) + +getCompactors :: PrimMonad m => ReqSketch (PrimState m) -> m (Vector.Vector (ReqCompactor (PrimState m))) +getCompactors = readMutVar . compactors + +computeTotalRetainedItems :: PrimMonad m => ReqSketch (PrimState m) -> m Int +computeTotalRetainedItems this = do + compactors <- getCompactors this + Vector.foldM countBuffer 0 compactors + 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 = readURef . retainedItems + +-- | Get the total number of items inserted into the sketch +count :: PrimMonad m => ReqSketch (PrimState m) -> m Word64 +count = readURef . totalN + +mkAuxiliaryFromReqSketch :: PrimMonad m => ReqSketch (PrimState m) -> m ReqAuxiliary +mkAuxiliaryFromReqSketch this = do + total <- count this + retainedItems <- retainedItemCount this + compactors <- getCompactors this + Auxiliary.mkAuxiliary (rankAccuracySetting this) total retainedItems compactors + +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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs @@ -0,0 +1,189 @@ +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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs @@ -0,0 +1,238 @@ +module DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor + ( ReqCompactor + , mkReqCompactor + , CompactorReturn (..) + , compact + , getBuffer + , getCoin + , getLgWeight + , getNominalCapacity + , getNumSections + , merge + , nearestEven + ) where + +import GHC.TypeLits +import Data.Bits ((.&.), (.|.), complement, countTrailingZeros, shiftL, shiftR) +import Data.Primitive.MutVar +import Data.Proxy +import Data.Semigroup (Semigroup) +import Data.Word +import DataSketches.Quantiles.RelativeErrorQuantile.Types +import System.Random.MWC (create, 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 +import DataSketches.Core.Snapshot + +data CompactorReturn s = CompactorReturn + { crDeltaRetItems :: {-# UNPACK #-} !Int + , crDeltaNominalSize :: {-# UNPACK #-} !Int + , crDoubleBuffer :: {-# UNPACK #-} !(DoubleBuffer s) + } + +data ReqCompactor s = ReqCompactor + -- Configuration constants + { rcRankAccuracy :: !RankAccuracy + , rcLgWeight :: {-# UNPACK #-} !Word8 + , rcRng :: {-# UNPACK #-} !(Gen s) + -- State + , rcState :: {-# UNPACK #-} !(URef s Word64) + , rcLastFlip :: {-# UNPACK #-} !(URef s Bool) + , rcSectionSizeFlt :: {-# UNPACK #-} !(URef s Double) + , rcSectionSize :: {-# UNPACK #-} !(URef s Word32) + , rcNumSections :: {-# UNPACK #-} !(URef s Word8) + , rcBuffer :: {-# UNPACK #-} !(MutVar s (DoubleBuffer s)) + } + +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 ReqCompactor{..} = ReqCompactorSnapshot rcRankAccuracy + <$> readURef rcState + <*> readURef rcLastFlip + <*> readURef rcSectionSizeFlt + <*> readURef rcSectionSize + <*> readURef rcNumSections + <*> (readMutVar rcBuffer >>= takeSnapshot) + + +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) + ReqCompactor rankAccuracy lgWeight g + <$> newURef 0 + <*> newURef False + <*> newURef (fromIntegral sectionSize) + <*> newURef sectionSize + <*> newURef initNumberOfSections + <*> newMutVar buff + +nomCapMult :: Num a => a +nomCapMult = 2 + +toInt :: Integral a => a -> Int +toInt = fromIntegral + +compact :: (PrimMonad m) => ReqCompactor (PrimState m) -> m (CompactorReturn (PrimState m)) +compact this = do + startBuffSize <- getCount =<< getBuffer this + startNominalCapacity <- getNominalCapacity this + numSections <- readURef $ rcNumSections this + sectionSize <- readURef $ rcSectionSize this + state <- readURef $ rcState this + let trailingOnes = succ $ countTrailingZeros $ complement state + sectionsToCompact = min trailingOnes $ fromIntegral numSections + (compactionStart, compactionEnd) <- computeCompactionRange this sectionsToCompact + -- TODO, this fails in GHCi but not in tests? + assert (compactionEnd - compactionStart >= 2) $ do + coin <- if state .&. 1 == 1 + then fmap not $ readURef $ rcLastFlip this + else flipCoin this + writeURef (rcLastFlip this) coin + buff <- getBuffer this + promote <- getEvensOrOdds buff compactionStart compactionEnd coin + trimCount buff $ startBuffSize - (compactionEnd - compactionStart) + modifyURef (rcState 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 + +flipCoin :: (PrimMonad m) => ReqCompactor (PrimState m) -> m Bool +flipCoin = uniform . rcRng + +getCoin :: PrimMonad m => ReqCompactor (PrimState m) -> m Bool +getCoin = readURef . rcLastFlip + +getNominalCapacity :: PrimMonad m => ReqCompactor (PrimState m) -> m Int +getNominalCapacity compactor = do + numSections <- readURef $ rcNumSections compactor + sectionSize <- readURef $ rcSectionSize compactor + pure $ nomCapMult * toInt numSections * toInt sectionSize + +getNumSections :: PrimMonad m => ReqCompactor (PrimState m) -> m Word8 +getNumSections = readURef . rcNumSections + +getSectionSizeFlt :: PrimMonad m => ReqCompactor (PrimState m) -> m Double +getSectionSizeFlt = readURef . rcSectionSizeFlt + +getState :: PrimMonad m => ReqCompactor (PrimState m) -> m Word64 +getState = readURef . rcState + +-- | Merge the other given compactor into this one. They both must have the +-- same @lgWeight@ +merge + :: (PrimMonad m, s ~ PrimState m) + => ReqCompactor (PrimState m) + -- ^ The compactor to merge into + -> ReqCompactor (PrimState m) + -- ^ The compactor to merge from + -> m (ReqCompactor s) +merge this otherCompactor = assert (rcLgWeight this == rcLgWeight otherCompactor) $ do + otherState <- readURef $ rcState otherCompactor + modifyURef (rcState this) (.|. otherState) + ensureMaxSections + + buff <- getBuffer this + sort buff + + otherBuff <- getBuffer otherCompactor + sort otherBuff + + otherBuffIsBigger <- (>) <$> getCount otherBuff <*> getCount buff + finalBuff <- 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 + -- loop until no adjustments can be made + when adjusted ensureMaxSections + +-- | Adjust the sectionSize and numSections if possible. +ensureEnoughSections + :: PrimMonad m + => ReqCompactor (PrimState m) + -> m Bool + -- ^ 'True' if the SectionSize and NumSections were adjusted. +ensureEnoughSections compactor = do + sectionSizeFlt <- readURef $ rcSectionSizeFlt compactor + let szf = sectionSizeFlt / sqrt2 + ne = nearestEven szf + state <- readURef $ rcState compactor + numSections <- readURef $ rcNumSections compactor + sectionSize <- readURef $ rcSectionSize compactor + if state >= (1 `shiftL` toInt (numSections - 1)) + && sectionSize > minK + && ne >= minK + then do + writeURef (rcSectionSizeFlt compactor) szf + writeURef (rcSectionSize compactor) $ fromIntegral ne + modifyURef (rcNumSections compactor) (`shiftL` 1) + buf <- getBuffer compactor + nomCapacity <- getNominalCapacity compactor + ensureCapacity buf (2 * nomCapacity) + pure True + else pure False + +-- | Computes the start and end indices of the compacted region +computeCompactionRange + :: PrimMonad m + => ReqCompactor (PrimState m) + -> Int + -- ^ secsToCompact the number of contiguous sections to compact + -> m (Int, Int) +-- ^ the start and end indices of the compacted region in compact form +computeCompactionRange this secsToCompact = do + buffSize <- getCount =<< getBuffer this + nominalCapacity <- getNominalCapacity this + numSections <- readURef $ rcNumSections this + sectionSize <- readURef $ rcSectionSize this + let nonCompact = (nominalCapacity `div` 2) + (fromIntegral numSections - secsToCompact) * fromIntegral sectionSize + nonCompact' = if (buffSize - nonCompact) .&. 1 == 1 then nonCompact - 1 else nonCompact + pure $ case rcRankAccuracy this of + HighRanksAreAccurate -> (0, fromIntegral $ buffSize - nonCompact') + LowRanksAreAccurate -> (nonCompact', buffSize) + +-- | Returns the nearest even integer to the given value. Also used by test. +nearestEven + :: Double + -- ^ the given value + -> Int + -- ^ the nearest even integer to the given value. +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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs @@ -0,0 +1,21 @@ +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) + +fixRseFactor :: Double +fixRseFactor = 0.084 diff --git a/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs @@ -0,0 +1,321 @@ +{-# 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 ) +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. +data DoubleBuffer s = DoubleBuffer + { vec :: {-# UNPACK #-} !(MutVar s (MUVector.MVector s Double)) + , count :: {-# UNPACK #-} !(URef s Int) + , sorted :: {-# UNPACK #-} !(URef s Bool) + , growthIncrement :: {-# UNPACK #-} !Int + , spaceAtBottom :: !Bool + } + +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{..} = DoubleBufferSnapshot + <$> (readMutVar vec >>= UVector.freeze) + <*> readURef count + <*> readURef sorted + <*> pure growthIncrement + <*> pure 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_ + count <- newURef 0 + sorted <- newURef True + pure $ DoubleBuffer{..} + +copyBuffer :: PrimMonad m => DoubleBuffer (PrimState m) -> m (DoubleBuffer (PrimState m)) +copyBuffer buf@DoubleBuffer{..} = do + vec <- newMutVar =<< MUVector.clone =<< getVector buf + count <- newURef =<< getCount buf + sorted <- newURef =<< readURef sorted + 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 + index <- if spaceAtBottom + then + (\capacity_ count_ -> capacity_ - count_ - 1) + <$> getCapacity buf + <*> getCount buf + else readURef count + modifyURef count (+ 1) + getVector buf >>= \vec -> MUVector.unsafeWrite vec index x + writeURef sorted False +{-# SCC 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_ <- readURef count + capacity_ <- getCapacity buf + let notEnoughSpace = count_ + space > capacity_ + when notEnoughSpace $ do + let newCap = count_ + space + growthIncrement + ensureCapacity buf newCap + +getVector :: (PrimMonad m, PrimState m ~ s) => DoubleBuffer s -> m (MUVector.MVector s Double) +getVector = readMutVar . vec +{-# INLINE getVector #-} + +getCapacity :: PrimMonad m => DoubleBuffer (PrimState m) -> m Int +getCapacity = fmap MUVector.length . getVector +{-# 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_) + + 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 + count <- newURef (MUVector.length out) + sorted <- newURef True + vec <- newMutVar out + pure DoubleBuffer + { vec = vec + , count = count + , sorted = sorted + , 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 = readURef . count + +getSpace :: PrimMonad m => DoubleBuffer (PrimState m) -> m Int +getSpace buf@DoubleBuffer{..} = (-) <$> getCapacity buf <*> getCount buf + +isEmpty :: PrimMonad m => DoubleBuffer (PrimState m) -> m Bool +isEmpty buf = (== 0) <$> getCount buf + +isSorted :: PrimMonad m => DoubleBuffer (PrimState m) -> m Bool +isSorted = readURef . sorted + +-- | 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_) + vec <- getVector buf + sortByBounds compare vec start end + writeURef sorted True +{-# SCC 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) + + modifyURef (count this) (+ bufInLen) + writeURef (sorted this) True + 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 = modifyURef count (\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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs @@ -0,0 +1,121 @@ +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 new file mode 100644 --- /dev/null +++ b/src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs @@ -0,0 +1,25 @@ +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.:<=) + + +data RankAccuracy + = HighRanksAreAccurate + -- ^ High ranks are prioritized for better accuracy. + | LowRanksAreAccurate + -- ^ Low ranks are prioritized for better accuracy + deriving (Show, Eq) diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,2 @@ +main :: IO () +main = putStrLn "Test suite not yet implemented"