data-sketches (empty) → 0.1.0.0
raw patch · 23 files changed
+2461/−0 lines, 23 filesdep +QuickCheckdep +basedep +cerealsetup-changedbinary-added
Dependencies added: QuickCheck, base, cereal, criterion, data-sketches, deepseq, ghc-prim, hspec, hspec-discover, mtl, mwc-random, pretty-show, prettyprinter, primitive, prometheus-client, statistics, vector, vector-algorithms
Files
- ChangeLog.md +3/−0
- LICENSE +206/−0
- README.md +16/−0
- Setup.hs +2/−0
- bench/Bench.hs +45/−0
- data-sketches.cabal +143/−0
- docs/images/KllErrorK100SL11.png binary
- docs/images/ReqErrorHraK12SL11_LT.png binary
- docs/images/ReqErrorLraK12SL11_LE.png binary
- src/DataSketches/Quantiles/RelativeErrorQuantile.hs +587/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs +177/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs +237/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs +21/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs +319/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs +121/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/URef.hs +49/−0
- src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs +29/−0
- test/AuxiliarySpec.hs +39/−0
- test/CompactorSpec.hs +9/−0
- test/DoubleBufferSpec.hs +116/−0
- test/ProofCheckSpec.hs +80/−0
- test/RelativeErrorQuantileSpec.hs +261/−0
- test/Spec.hs +1/−0
+ ChangeLog.md view
@@ -0,0 +1,3 @@+# Changelog for streaming-quantiles++## Unreleased changes
+ LICENSE view
@@ -0,0 +1,206 @@++ Apache License+ Version 2.0, January 2004+ http://www.apache.org/licenses/++ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION++ 1. Definitions.++ "License" shall mean the terms and conditions for use, reproduction,+ and distribution as defined by Sections 1 through 9 of this document.++ "Licensor" shall mean the copyright owner or entity authorized by+ the copyright owner that is granting the License.++ "Legal Entity" shall mean the union of the acting entity and all+ other entities that control, are controlled by, or are under common+ control with that entity. For the purposes of this definition,+ "control" means (i) the power, direct or indirect, to cause the+ direction or management of such entity, whether by contract or+ otherwise, or (ii) ownership of fifty percent (50%) or more of the+ outstanding shares, or (iii) beneficial ownership of such entity.++ "You" (or "Your") shall mean an individual or Legal Entity+ exercising permissions granted by this License.++ "Source" form shall mean the preferred form for making modifications,+ including but not limited to software source code, documentation+ source, and configuration files.++ "Object" form shall mean any form resulting from mechanical+ transformation or translation of a Source form, including but+ not limited to compiled object code, generated documentation,+ and conversions to other media types.++ "Work" shall mean the work of authorship, whether in Source or+ Object form, made available under the License, as indicated by a+ copyright notice that is included in or attached to the work+ (an example is provided in the Appendix below).++ "Derivative Works" shall mean any work, whether in Source or Object+ form, that is based on (or derived from) the Work and for which the+ editorial revisions, annotations, elaborations, or other modifications+ represent, as a whole, an original work of authorship. For the purposes+ of this License, Derivative Works shall not include works that remain+ separable from, or merely link (or bind by name) to the interfaces of,+ the Work and Derivative Works thereof.++ "Contribution" shall mean any work of authorship, including+ the original version of the Work and any modifications or additions+ to that Work or Derivative Works thereof, that is intentionally+ submitted to Licensor for inclusion in the Work by the copyright owner+ or by an individual or Legal Entity authorized to submit on behalf of+ the copyright owner. For the purposes of this definition, "submitted"+ means any form of electronic, verbal, or written communication sent+ to the Licensor or its representatives, including but not limited to+ communication on electronic mailing lists, source code control systems,+ and issue tracking systems that are managed by, or on behalf of, the+ Licensor for the purpose of discussing and improving the Work, but+ excluding communication that is conspicuously marked or otherwise+ designated in writing by the copyright owner as "Not a Contribution."++ "Contributor" shall mean Licensor and any individual or Legal Entity+ on behalf of whom a Contribution has been received by Licensor and+ subsequently incorporated within the Work.++ 2. Grant of Copyright License. Subject to the terms and conditions of+ this License, each Contributor hereby grants to You a perpetual,+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable+ copyright license to reproduce, prepare Derivative Works of,+ publicly display, publicly perform, sublicense, and distribute the+ Work and such Derivative Works in Source or Object form.++ 3. Grant of Patent License. Subject to the terms and conditions of+ this License, each Contributor hereby grants to You a perpetual,+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable+ (except as stated in this section) patent license to make, have made,+ use, offer to sell, sell, import, and otherwise transfer the Work,+ where such license applies only to those patent claims licensable+ by such Contributor that are necessarily infringed by their+ Contribution(s) alone or by combination of their Contribution(s)+ with the Work to which such Contribution(s) was submitted. If You+ institute patent litigation against any entity (including a+ cross-claim or counterclaim in a lawsuit) alleging that the Work+ or a Contribution incorporated within the Work constitutes direct+ or contributory patent infringement, then any patent licenses+ granted to You under this License for that Work shall terminate+ as of the date such litigation is filed.++ 4. Redistribution. You may reproduce and distribute copies of the+ Work or Derivative Works thereof in any medium, with or without+ modifications, and in Source or Object form, provided that You+ meet the following conditions:++ (a) You must give any other recipients of the Work or+ Derivative Works a copy of this License; and++ (b) You must cause any modified files to carry prominent notices+ stating that You changed the files; and++ (c) You must retain, in the Source form of any Derivative Works+ that You distribute, all copyright, patent, trademark, and+ attribution notices from the Source form of the Work,+ excluding those notices that do not pertain to any part of+ the Derivative Works; and++ (d) If the Work includes a "NOTICE" text file as part of its+ distribution, then any Derivative Works that You distribute must+ include a readable copy of the attribution notices contained+ within such NOTICE file, excluding those notices that do not+ pertain to any part of the Derivative Works, in at least one+ of the following places: within a NOTICE text file distributed+ as part of the Derivative Works; within the Source form or+ documentation, if provided along with the Derivative Works; or,+ within a display generated by the Derivative Works, if and+ wherever such third-party notices normally appear. The contents+ of the NOTICE file are for informational purposes only and+ do not modify the License. You may add Your own attribution+ notices within Derivative Works that You distribute, alongside+ or as an addendum to the NOTICE text from the Work, provided+ that such additional attribution notices cannot be construed+ as modifying the License.++ You may add Your own copyright statement to Your modifications and+ may provide additional or different license terms and conditions+ for use, reproduction, or distribution of Your modifications, or+ for any such Derivative Works as a whole, provided Your use,+ reproduction, and distribution of the Work otherwise complies with+ the conditions stated in this License.++ 5. Submission of Contributions. Unless You explicitly state otherwise,+ any Contribution intentionally submitted for inclusion in the Work+ by You to the Licensor shall be under the terms and conditions of+ this License, without any additional terms or conditions.+ Notwithstanding the above, nothing herein shall supersede or modify+ the terms of any separate license agreement you may have executed+ with Licensor regarding such Contributions.++ 6. Trademarks. This License does not grant permission to use the trade+ names, trademarks, service marks, or product names of the Licensor,+ except as required for reasonable and customary use in describing the+ origin of the Work and reproducing the content of the NOTICE file.++ 7. Disclaimer of Warranty. Unless required by applicable law or+ agreed to in writing, Licensor provides the Work (and each+ Contributor provides its Contributions) on an "AS IS" BASIS,+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or+ implied, including, without limitation, any warranties or conditions+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A+ PARTICULAR PURPOSE. You are solely responsible for determining the+ appropriateness of using or redistributing the Work and assume any+ risks associated with Your exercise of permissions under this License.++ 8. Limitation of Liability. In no event and under no legal theory,+ whether in tort (including negligence), contract, or otherwise,+ unless required by applicable law (such as deliberate and grossly+ negligent acts) or agreed to in writing, shall any Contributor be+ liable to You for damages, including any direct, indirect, special,+ incidental, or consequential damages of any character arising as a+ result of this License or out of the use or inability to use the+ Work (including but not limited to damages for loss of goodwill,+ work stoppage, computer failure or malfunction, or any and all+ other commercial damages or losses), even if such Contributor+ has been advised of the possibility of such damages.++ 9. Accepting Warranty or Additional Liability. While redistributing+ the Work or Derivative Works thereof, You may choose to offer,+ and charge a fee for, acceptance of support, warranty, indemnity,+ or other liability obligations and/or rights consistent with this+ License. However, in accepting such obligations, You may act only+ on Your own behalf and on Your sole responsibility, not on behalf+ of any other Contributor, and only if You agree to indemnify,+ defend, and hold each Contributor harmless for any liability+ incurred by, or claims asserted against, such Contributor by reason+ of your accepting any such warranty or additional liability.++ END OF TERMS AND CONDITIONS++++APPENDIX A: How to apply the Apache License to your work.++ To apply the Apache License to your work, attach the following+ boilerplate notice, with the fields enclosed by brackets "[]"+ replaced with your own identifying information. (Don't include+ the brackets!) The text should be enclosed in the appropriate+ comment syntax for the file format. We also recommend that a+ file or class name and description of purpose be included on the+ same "printed page" as the copyright notice for easier+ identification within third-party archives.+ -------------------------------------------------------------+ Copyright [yyyy] [name of copyright owner]++ Licensed under the Apache License, Version 2.0 (the "License");+ you may not use this file except in compliance with the License.+ You may obtain a copy of the License at++ http://www.apache.org/licenses/LICENSE-2.0++ Unless required by applicable law or agreed to in writing, software+ distributed under the License is distributed on an "AS IS" BASIS,+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.+ See the License for the specific language governing permissions and+ limitations under the License.+ -------------------------------------------------------------+
+ README.md view
@@ -0,0 +1,16 @@+# streaming-quantiles++The Business Challenge: Analyzing Big Data Quickly.++In the analysis of big data there are often problem queries that don’t scale because they require huge compute resources and time to generate exact results. Examples include count distinct, quantiles, most-frequent items, joins, matrix computations, and graph analysis.++If approximate results are acceptable, there is a class of specialized algorithms, called streaming algorithms, or sketches that can produce results orders-of magnitude faster and with mathematically proven error bounds. For interactive queries there may not be other viable alternatives, and in the case of real-time analysis, sketches are the only known solution.++For any system that needs to extract useful information from big data these sketches are a required toolkit that should be tightly integrated into their analysis capabilities. This technology has helped Yahoo (Verizon Media) successfully reduce data processing times from days or hours to minutes or seconds on a number of its internal platforms.++This project is dedicated to providing a broad selection of sketch algorithms of production quality. Contributions are welcome from those interested in further development of this science and art.++## Why use this project?++- Sketches are fast. The sketch algorithms in this library process data in a single pass and are suitable for both real-time and batch. Sketches enable streaming computation of set expression cardinalities, quantiles, frequency estimation and more. In addition, designing a system around sketching allows simplification of system's architecture and reduction in overall compute resources required for these heretofore difficult computation+- Built-in Theta Sketch set operators (Union, Intersection, Difference) produce sketches as a result (and not just a number) enabling full set expressions of cardinality, such as ((A ∪ B) ∩ (C ∪ D)) \ (E ∪ F). This capability along with predictable and superior accuracy (compared with Include/Exclude approaches) enable unprecedented analysis capabilities for fast queries.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ bench/Bench.hs view
@@ -0,0 +1,45 @@+{-# OPTIONS_GHC -fno-full-laziness #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE TypeApplications #-}++import Control.Monad (forM_)+import Criterion.Main+import Criterion.Types+import Data.List (foldl')+import DataSketches.Quantiles.RelativeErrorQuantile+import Prometheus (register, summary, defaultQuantiles, observe, Info (Info))+import Control.Concurrent (withMVar, newMVar)++main :: IO ()+main = do+ outerSketch <- mkReqSketch @6 HighRanksAreAccurate+ let metric = summary (Info "adversarial_input" "woo") defaultQuantiles+ prometheusThing <- register metric+ skM <- newMVar =<< mkReqSketch @6 HighRanksAreAccurate+ -- mapM_ (update outerSketch) [1..10000]+ defaultMain+ [ bgroup "ReqSketch"+ [ bench "insert/1" $ perRunEnv (mkReqSketch @6 HighRanksAreAccurate) $ \sk -> do+ update sk 1+ , bench "insert/10" $ perRunEnv (mkReqSketch @6 HighRanksAreAccurate) $ \sk -> do+ mapM_ (update sk) [1..10]+ , bench "insert/100" $ perRunEnv (mkReqSketch @6 HighRanksAreAccurate) $ \sk -> do+ mapM_ (update sk) [1..100]+ , bench "insert/1000" $ perRunEnv (mkReqSketch @6 HighRanksAreAccurate) $ \sk -> do+ mapM_ (update sk) [1..1000]+ , bench "insert/10000" $ perRunEnv (mkReqSketch @6 HighRanksAreAccurate) $ \sk -> do+ mapM_ (update sk) [1..10000]+ , bench "insert/existing" $ whnfIO $ update outerSketch 1+ , bench "insert/mvar" $ whnfIO $ withMVar skM (`update` 1)+ ]+ , bgroup "DoubleBuffer"+ [ -- bench "sort" $ + ]+ , bgroup "Prometheus"+ [ bench "insert/existing" $ whnfIO $+ observe prometheusThing 1+ ]+ ]
+ data-sketches.cabal view
@@ -0,0 +1,143 @@+cabal-version: 1.18++-- This file has been generated from package.yaml by hpack version 0.34.4.+--+-- see: https://github.com/sol/hpack++name: data-sketches+version: 0.1.0.0+description: Please see the README on GitHub at <https://github.com/iand675/streaming-quantiles#readme>+homepage: https://github.com/iand675/data-sketches#readme+bug-reports: https://github.com/iand675/data-sketches/issues+author: Ian Duncan, Rob Bassi+maintainer: ian@iankduncan.com+copyright: 2021 Ian Duncan, Rob Bassi, Mercury Bank+license: Apache+license-file: LICENSE+build-type: Simple+extra-source-files:+ README.md+ ChangeLog.md+extra-doc-files:+ docs/images/KllErrorK100SL11.png+ docs/images/ReqErrorHraK12SL11_LT.png+ docs/images/ReqErrorLraK12SL11_LE.png++source-repository head+ type: git+ location: https://github.com/iand675/data-sketches++library+ exposed-modules:+ DataSketches.Quantiles.RelativeErrorQuantile+ 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.Internal.URef+ DataSketches.Quantiles.RelativeErrorQuantile.Types+ other-modules:+ Paths_data_sketches+ hs-source-dirs:+ src+ default-extensions:+ BangPatterns+ DataKinds+ FlexibleInstances+ KindSignatures+ RecordWildCards+ ScopedTypeVariables+ StandaloneDeriving+ TypeFamilies+ TypeOperators+ build-depends:+ base >=4.7 && <5+ , cereal+ , deepseq+ , ghc-prim+ , mtl+ , mwc-random+ , pretty-show+ , prettyprinter+ , primitive+ , vector+ , vector-algorithms+ default-language: Haskell2010++test-suite data-sketches-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ AuxiliarySpec+ CompactorSpec+ DoubleBufferSpec+ ProofCheckSpec+ RelativeErrorQuantileSpec+ Paths_data_sketches+ hs-source-dirs:+ test+ default-extensions:+ BangPatterns+ DataKinds+ FlexibleInstances+ KindSignatures+ RecordWildCards+ ScopedTypeVariables+ StandaloneDeriving+ TypeFamilies+ TypeOperators+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ QuickCheck+ , base >=4.7 && <5+ , cereal+ , data-sketches+ , deepseq+ , ghc-prim+ , hspec+ , hspec-discover+ , mtl+ , mwc-random+ , pretty-show+ , prettyprinter+ , primitive+ , statistics+ , vector+ , vector-algorithms+ default-language: Haskell2010++benchmark bench+ type: exitcode-stdio-1.0+ main-is: Bench.hs+ other-modules:+ Paths_data_sketches+ hs-source-dirs:+ bench+ default-extensions:+ BangPatterns+ DataKinds+ FlexibleInstances+ KindSignatures+ RecordWildCards+ ScopedTypeVariables+ StandaloneDeriving+ TypeFamilies+ TypeOperators+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ build-depends:+ base >=4.7 && <5+ , cereal+ , criterion+ , data-sketches+ , deepseq+ , ghc-prim+ , mtl+ , mwc-random+ , pretty-show+ , prettyprinter+ , primitive+ , prometheus-client+ , vector+ , vector-algorithms+ default-language: Haskell2010
+ docs/images/KllErrorK100SL11.png view
binary file changed (absent → 49325 bytes)
+ docs/images/ReqErrorHraK12SL11_LT.png view
binary file changed (absent → 69100 bytes)
+ docs/images/ReqErrorLraK12SL11_LE.png view
binary file changed (absent → 67275 bytes)
+ src/DataSketches/Quantiles/RelativeErrorQuantile.hs view
@@ -0,0 +1,587 @@+-- | The Relative Error Quantile (REQ) sketch provides extremely high accuracy at a chosen end of the rank domain. +-- This is best illustrated with some rank domain accuracy plots that compare the KLL quantiles sketch to the REQ sketch.+--+-- This first plot illustrates the typical error behavior of the KLL sketch (also the quantiles/DoublesSketch). +-- The error is flat for all ranks (0, 1). The green and yellow lines correspond to +/- one RSE at 68% confidence; +-- the blue and red lines, +/- two RSE at 95% confidence; and, the purple and brown lines +/- 3 RSE at 99% confidence. +-- The reason all the curves pinch at 0 and 1.0, is because the sketch knows with certainty that a request for a quantile at +-- rank = 0 is the minimum value of the stream; and a request for a quantiles at rank = 1.0, is the maximum value of the stream. +-- Both of which the sketch tracks.+--+-- +--+-- The next plot is the exact same data and queries fed to the REQ sketch set for High Rank Accuracy (HRA) mode. +-- In this plot, starting at a rank of about 0.3, the contour lines start converging and actually reach zero error at +-- rank 1.0. Therefore the error (the inverse of accuracy) is relative to the requested rank, thus the name of the sketch. +-- This means that the user can perform getQuantile(rank) queries, where rank = .99999 and get accurate results.+--+-- +--+-- This next plot is also the same data and queries, except the REQ sketch was configured for Low Rank Accuracy (LRA). In this case the user can perform getQuantiles(rank) queries, where rank = .00001 and get accurate results.+--+-- ++{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE DeriveGeneric #-}+module DataSketches.Quantiles.RelativeErrorQuantile+ ( ReqSketch (..)+ , ValidK+ , mkReqSketch+ , getN+ , getSum+ , getCompactors+ , getNumLevels+ , getRetainedItems+ , computeTotalRetainedItems+ , cumulativeDistributionFunction+ , RankAccuracy(..)+ , rankAccuracy+ , relativeStandardError+ , maximum+ , count+ , probabilityMassFunction+ , quantile+ , quantiles+ , rank+ , rankLowerBound+ , ranks+ , rankUpperBound+ , null+ , isEstimationMode+ , isLessThanOrEqual+ , merge+ , update+ , CumulativeDistributionInvariants(..)+ -- Test only?+ , mkAuxiliaryFromReqSketch+ ) where++import Control.DeepSeq+import Control.Monad (when, unless, foldM, foldM_)+import Control.Monad.Primitive+import Control.Monad.Trans+import Data.Bits (shiftL)+import Data.Vector ((!))+import qualified Data.Vector as Vector+import Data.Primitive.MutVar+import Data.Proxy+import Data.Word+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Constants+import DataSketches.Quantiles.RelativeErrorQuantile.Types+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor (ReqCompactor)+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary (ReqAuxiliary)+import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary as Auxiliary+import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor as Compactor+import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer as DoubleBuffer+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.URef+import GHC.TypeLits+import Prelude hiding (null, minimum, maximum)+import Data.Maybe (isNothing)+import qualified Data.Foldable+import qualified Data.List+import GHC.Exception.Type (Exception)+import Control.Exception (throw, assert)+import GHC.Generics+import System.Random.MWC (Gen, create)++data ReqSketch (k :: Nat) s = ReqSketch+ { 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 k s) where+ rnf !rs = ()++instance TakeSnapshot (ReqSketch k) where+ data Snapshot (ReqSketch k) = ReqSketchSnapshot+ { snapshotRankAccuracySetting :: !RankAccuracy+ , snapshotCriterion :: !Criterion+ , snapshotTotalN :: !Word64+ , snapshotMinValue :: !Double+ , snapshotMaxValue :: !Double+ , snapshotRetainedItems :: !Int+ , snapshotMaxNominalCapacitiesSize :: !Int+ -- , aux :: !(MutVar s (Maybe ()))+ , snapshotCompactors :: !(Vector.Vector (Snapshot ReqCompactor))+ }+ takeSnapshot ReqSketch{..} = ReqSketchSnapshot rankAccuracySetting criterion+ <$> readURef totalN+ <*> readURef minValue+ <*> readURef maxValue+ <*> readURef retainedItems+ <*> readURef maxNominalCapacitiesSize+ <*> (readMutVar compactors >>= mapM takeSnapshot)++deriving instance Show (Snapshot (ReqSketch k))++type ValidK k = (4 <= k, k <= 1024, (k `Mod` 2) ~ 0)++mkReqSketch :: forall k m. (PrimMonad m, ValidK k, KnownNat k) => RankAccuracy -> m (ReqSketch k (PrimState m))+mkReqSketch rank = do+ r <- ReqSketch rank (:<)+ <$> create+ <*> newURef 0+ <*> newURef (0 / 0)+ <*> newURef (0 / 0)+ <*> newURef 0+ <*> newURef 0+ <*> newURef 0+ <*> newMutVar Nothing+ <*> newMutVar Vector.empty+ grow r+ pure r++mkAuxiliaryFromReqSketch :: PrimMonad m => ReqSketch k (PrimState m) -> m ReqAuxiliary+mkAuxiliaryFromReqSketch this = do+ total <- getN this+ retainedItems <- getRetainedItems this+ compactors <- getCompactors this+ Auxiliary.mkAuxiliary (rankAccuracySetting this) total retainedItems compactors++getAux :: PrimMonad m => ReqSketch k (PrimState m) -> m (Maybe ReqAuxiliary)+getAux = readMutVar . aux++getCompactors :: PrimMonad m => ReqSketch k (PrimState m) -> m (Vector.Vector (ReqCompactor (PrimState m)))+getCompactors = readMutVar . compactors++getNumLevels :: PrimMonad m => ReqSketch k (PrimState m) -> m Int+getNumLevels = fmap Vector.length . getCompactors++getIsEmpty :: PrimMonad m => ReqSketch k (PrimState m) -> m Bool+getIsEmpty = fmap (== 0) . readURef . totalN++getN :: PrimMonad m => ReqSketch k (PrimState m) -> m Word64+getN = readURef . totalN++getK :: forall k s. KnownNat k => ReqSketch k s -> Word32+getK _ = fromIntegral (natVal (Proxy :: Proxy k))++getRetainedItems :: PrimMonad m => ReqSketch k (PrimState m) -> m Int+getRetainedItems = readURef . retainedItems++getMaxNominalCapacity :: PrimMonad m => ReqSketch k (PrimState m) -> m Int+getMaxNominalCapacity = readURef . maxNominalCapacitiesSize++data CumulativeDistributionInvariants+ = CumulativeDistributionInvariantsSplitsAreEmpty+ | CumulativeDistributionInvariantsSplitsAreNotFinite+ | CumulativeDistributionInvariantsSplitsAreNotUniqueAndMontonicallyIncreasing+ deriving (Show, Eq)++instance Exception CumulativeDistributionInvariants++validateSplits :: Monad m => [Double] -> m ()+validateSplits splits = do+ when (Data.Foldable.null splits) $ do+ throw CumulativeDistributionInvariantsSplitsAreEmpty+ when (any isInfinite splits || any isNaN splits) $ do+ throw CumulativeDistributionInvariantsSplitsAreNotFinite+ when (Data.List.nub (Data.List.sort splits) /= splits) $ do+ throw CumulativeDistributionInvariantsSplitsAreNotUniqueAndMontonicallyIncreasing++getCounts :: (PrimMonad m, KnownNat k) => ReqSketch k (PrimState m) -> [Double] -> m [Word64]+getCounts this values = do+ compactors <- getCompactors this+ let numValues = length values+ numCompactors = Vector.length compactors+ ans = replicate numValues 0+ isEmpty <- getIsEmpty this+ if isEmpty+ then pure []+ else Vector.ifoldM doCount ans compactors+ where+ doCount acc index compactor = do+ let wt = (1 `shiftL` fromIntegral (Compactor.getLgWeight compactor)) :: Word64+ buff <- Compactor.getBuffer compactor+ let updateCounts buff value = do+ count_ <- DoubleBuffer.getCountWithCriterion buff (values !! index) (criterion this)+ pure $ fromIntegral value + fromIntegral count_ * wt+ mapM (updateCounts buff) acc++getPMForCDF :: (PrimMonad m, KnownNat k) => ReqSketch k (PrimState m) -> [Double] -> m [Word64]+getPMForCDF this splits = do+ () <- validateSplits splits+ let numSplits = length splits+ numBuckets = numSplits -- + 1+ splitCounts <- getCounts this splits+ n <- getN this+ pure $ (++ [n]) $ take numBuckets splitCounts++-- | Returns an approximation to the Cumulative Distribution Function (CDF), which is the cumulative analog of the PMF, +-- of the input stream given a set of splitPoint (values).+cumulativeDistributionFunction+ :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> [Double]+ -- ^ Returns an approximation to the Cumulative Distribution Function (CDF), + -- which is the cumulative analog of the PMF, of the input stream given a set of + -- splitPoint (values).+ --+ -- The resulting approximations have a probabilistic guarantee that be obtained, + -- a priori, from the getRSE(int, double, boolean, long) function.+ --+ -- If the sketch is empty this returns 'Nothing'.+ -> m (Maybe [Double])+cumulativeDistributionFunction this splitPoints = do+ buckets <- getPMForCDF this splitPoints+ isEmpty <- getIsEmpty this+ if isEmpty+ then pure Nothing+ else do+ let numBuckets = length splitPoints + 1+ n <- getN this+ pure $ Just $ (/ fromIntegral n) . fromIntegral <$> buckets++rankAccuracy :: ReqSketch s k -> RankAccuracy+rankAccuracy = rankAccuracySetting++-- | Returns an a priori estimate of relative standard error (RSE, expressed as a number in [0,1]). Derived from Lemma 12 in https://arxiv.org/abs/2004.01668v2, but the constant factors were modified based on empirical measurements.+relativeStandardError+ :: ReqSketch s k+ -> Int+ -- ^ k - the given value of k+ -> Double+ -- ^ rank - the given normalized rank, a number in [0,1].+ -> RankAccuracy+ -> Int+ -- ^ totalN - an estimate of the total number of items submitted to the sketch.+ -> Double+ -- ^ an a priori estimate of relative standard error (RSE, expressed as a number in [0,1]).+relativeStandardError = undefined++-- | Gets the smallest value seen by this sketch+minimum :: PrimMonad m => ReqSketch k (PrimState m) -> m Double+minimum = readURef . minValue++-- | Gets the largest value seen by this sketch+maximum :: PrimMonad m => ReqSketch k (PrimState m) -> m Double+maximum = readURef . maxValue++count :: (PrimMonad m, s ~ PrimState m, KnownNat k) => ReqSketch k s -> Double -> m Word64+count s value = fromIntegral <$> do+ empty <- null s+ if empty+ then pure 0+ else do+ compactors <- getCompactors s+ let go !accum compactor = do+ let wt = (1 `shiftL` fromIntegral (Compactor.getLgWeight compactor)) :: Word64+ buf <- Compactor.getBuffer compactor+ count_ <- DoubleBuffer.getCountWithCriterion buf value (criterion s)+ pure (accum + (fromIntegral count_ * wt))+ Vector.foldM go 0 compactors++getSum :: (PrimMonad m) => ReqSketch k (PrimState m) -> m Double+getSum = readURef . sumValue++-- | Returns an approximation to the Probability Mass Function (PMF) of the input stream given a set of splitPoints (values).+-- The resulting approximations have a probabilistic guarantee that be obtained, a priori, from the getRSE(int, double, boolean, long) function.+--+-- If the sketch is empty this returns an empty list.+probabilityMassFunction+ :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> [Double]+ -- ^ splitPoints - an array of m unique, monotonically increasing double values that divide + -- the real number line into m+1 consecutive disjoint intervals. The definition of an "interval" + -- is inclusive of the left splitPoint (or minimum value) and exclusive of the right splitPoint, + -- with the exception that the last interval will include the maximum value. It is not necessary + -- to include either the min or max values in these splitpoints.+ -> m [Double]+ -- ^ An array of m+1 doubles each of which is an approximation to the fraction of + -- the input stream values (the mass) that fall into one of those intervals. + -- The definition of an "interval" is inclusive of the left splitPoint and exclusive + -- of the right splitPoint, with the exception that the last interval will + -- include maximum value.+probabilityMassFunction this splitPoints = do+ isEmpty <- getIsEmpty this+ if isEmpty+ then pure []+ else do+ let numBuckets = length splitPoints + 1+ buckets <- fmap fromIntegral <$> getPMForCDF this splitPoints+ total <- fromIntegral <$> getN this+ let computeProb (0, bucket) = bucket / total+ computeProb (i, bucket) = (prevBucket + bucket) / total+ where prevBucket = buckets !! i - 1+ probs = computeProb <$> zip [0..] buckets+ pure probs++-- | Gets the approximate quantile of the given normalized rank based on the lteq criterion.+quantile+ :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> Double+ -- ^ normRank - the given normalized rank+ -> m Double+ -- ^ the approximate quantile given the normalized rank.+quantile this normRank = do+ isEmpty <- getIsEmpty this+ if isEmpty+ then pure (0/0)+ else do+ when (normRank < 0 || normRank > 1.0) $+ error $ "Normalized rank must be in the range [0.0, 1.0]: " ++ show normRank+ currAuxiliary <- getAux this+ when (isNothing currAuxiliary) $ do+ total <- getN this+ retainedItems <- getRetainedItems this+ compactors <- getCompactors this+ newAuxiliary <- Auxiliary.mkAuxiliary (rankAccuracySetting this) total retainedItems compactors+ writeMutVar (aux this) (Just newAuxiliary)+ mAuxiliary <- getAux this+ case mAuxiliary of+ Just auxiliary -> pure $! Auxiliary.getQuantile auxiliary normRank $ criterion this+ Nothing -> error "invariant violated: aux is not set"++-- | Gets an array of quantiles that correspond to the given array of normalized ranks.+quantiles+ :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> [Double]+ -- ^ normRanks - the given array of normalized ranks.+ -> m [Double]+ -- ^ the array of quantiles that correspond to the given array of normalized ranks.+quantiles this normRanks = do+ isEmpty <- getIsEmpty this+ if isEmpty+ then pure []+ else mapM (quantile this) normRanks++-- | Computes the normalized rank of the given value in the stream. The normalized rank is the fraction of values less than the given value; or if lteq is true, the fraction of values less than or equal to the given value.+rank :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> Double+ -- ^ value - the given value+ -> m Double+ -- ^ the normalized rank of the given value in the stream.+rank s value = do+ isEmpty <- null s+ if isEmpty+ then pure (0 / 0) -- NaN+ else do+ nnCount <- count s value+ total <- readURef $ totalN s+ pure (fromIntegral nnCount / fromIntegral total)+++-- getRankLB k levels rank numStdDev hra totalN = if exactRank k levels rank hra totalN++-- | Returns an approximate lower bound rank of the given normalized rank.+rankLowerBound+ :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> Double+ -- ^ rank - the given rank, a value between 0 and 1.0.+ -> Int+ -- ^ numStdDev - the number of standard deviations. Must be 1, 2, or 3.+ -> m Double+ -- ^ an approximate lower bound rank.+rankLowerBound this rank numStdDev = do+ numLevels <- getNumLevels this+ let k = fromIntegral $ getK this+ total <- getN this+ pure $ getRankLB k numLevels rank numStdDev (rankAccuracySetting this == HighRanksAreAccurate) total++-- | Gets an array of normalized ranks that correspond to the given array of values.+-- TODO, make it ifaster+ranks :: (PrimMonad m, s ~ PrimState m, KnownNat k) => ReqSketch k s -> [Double] -> m [Double]+ranks s values = mapM (rank s) values++-- | Returns an approximate upper bound rank of the given rank.+rankUpperBound+ :: (PrimMonad m, KnownNat k)+ => ReqSketch k (PrimState m)+ -> Double+ -- ^ rank - the given rank, a value between 0 and 1.0.+ -> Int+ -- ^ numStdDev - the number of standard deviations. Must be 1, 2, or 3.+ -> m Double+ -- ^ an approximate upper bound rank.+rankUpperBound this rank numStdDev= do+ numLevels <- getNumLevels this+ let k = fromIntegral $ getK this+ total <- getN this+ pure $ getRankUB k numLevels rank numStdDev (rankAccuracySetting this == HighRanksAreAccurate) total++-- | Returns true if this sketch is empty.+null :: (PrimMonad m) => ReqSketch k (PrimState m) -> m Bool+null = fmap (== 0) . readURef . totalN++-- | Returns true if this sketch is in estimation mode.+isEstimationMode :: PrimMonad m => ReqSketch k (PrimState m) -> m Bool+isEstimationMode = fmap (> 1) . getNumLevels++-- | Returns the current comparison criterion.+isLessThanOrEqual :: ReqSketch s k -> Bool+isLessThanOrEqual s = case criterion s of+ (:<) -> False+ (:<=) -> True++computeMaxNominalSize :: PrimMonad m => ReqSketch k (PrimState m) -> m Int+computeMaxNominalSize this = do+ compactors <- getCompactors this+ Vector.foldM countNominalCapacity 0 compactors+ where+ countNominalCapacity acc compactor = do+ nominalCapacity <- Compactor.getNominalCapacity compactor+ pure $ nominalCapacity + acc++computeTotalRetainedItems :: PrimMonad m => ReqSketch k (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++grow :: (PrimMonad m, KnownNat k) => ReqSketch k (PrimState m) -> m ()+grow this = do+ lgWeight <- fromIntegral <$> getNumLevels this+ let rankAccuracy = rankAccuracySetting this+ sectionSize = getK this+ newCompactor <- Compactor.mkReqCompactor (sketchRng this) lgWeight rankAccuracy sectionSize+ modifyMutVar' (compactors this) (`Vector.snoc` newCompactor)+ maxNominalCapacity <- computeMaxNominalSize this+ writeURef (maxNominalCapacitiesSize this) maxNominalCapacity++compress :: (PrimMonad m, KnownNat k) => ReqSketch k (PrimState m) -> m ()+compress this = do+ compactors <- getCompactors this+ Vector.iforM_ compactors $ \height compactor -> do+ buffSize <- DoubleBuffer.getCount =<< Compactor.getBuffer compactor+ nominalCapacity <- Compactor.getNominalCapacity compactor+ when (buffSize >= nominalCapacity) $ do+ numLevels <- getNumLevels this+ when (height + 1 >= numLevels) $ do+ grow this+ compactors' <- getCompactors this+ cReturn <- Compactor.compact compactor+ let topCompactor = compactors' ! (height + 1)+ buff <- Compactor.getBuffer topCompactor+ DoubleBuffer.mergeSortIn buff $ Compactor.crDoubleBuffer cReturn+ modifyURef (retainedItems this) (+ Compactor.crDeltaRetItems cReturn)+ modifyURef (maxNominalCapacitiesSize this) (+ Compactor.crDeltaNominalSize cReturn)+ writeMutVar (aux this) Nothing++-- | Merge other sketch into this one.+merge+ :: (PrimMonad m, s ~ PrimState m, KnownNat k)+ => ReqSketch k s+ -> ReqSketch k s+ -> m (ReqSketch k s)+merge this other = do+ otherIsEmpty <- getIsEmpty other+ unless otherIsEmpty $ do+ let rankAccuracy = rankAccuracySetting this+ otherRankAccuracy = rankAccuracySetting other+ when (rankAccuracy /= otherRankAccuracy) $+ error "Both sketches must have the same HighRankAccuracy setting."+ -- update total+ otherN <- getN other+ modifyURef (totalN this) (+ otherN)+ -- update the min and max values+ thisMin <- minimum this+ thisMax <- maximum this+ otherMin <- minimum other+ otherMax <- maximum other+ when (isNaN thisMin || otherMin < thisMin) $ do+ writeURef (minValue this) otherMin+ when (isNaN thisMax || otherMax < thisMax) $ do+ writeURef (maxValue this) otherMax+ -- grow until this has at least as many compactors as other+ numRequiredCompactors <- getNumLevels other+ growUntil numRequiredCompactors+ -- merge the items in all height compactors+ thisCompactors <- getCompactors this+ otherCompactors <- getCompactors other+ Vector.zipWithM_ Compactor.merge thisCompactors otherCompactors+ -- update state+ maxNominalCapacity <- computeMaxNominalSize this+ totalRetainedItems <- computeTotalRetainedItems this+ writeURef (maxNominalCapacitiesSize this) maxNominalCapacity+ writeURef (retainedItems this) totalRetainedItems+ -- compress and check invariants+ when (totalRetainedItems >= maxNominalCapacity) $ do+ compress this+ maxNominalCapacity' <- readURef $ maxNominalCapacitiesSize this+ totalRetainedItems' <- readURef $ retainedItems this+ assert (totalRetainedItems' < maxNominalCapacity') $ + writeMutVar (aux this) Nothing+ pure this+ where+ growUntil target = do+ numCompactors <- getNumLevels this+ when (numCompactors < target) $+ grow this++-- | Updates this sketch with the given item.+update :: (PrimMonad m, KnownNat k) => ReqSketch k (PrimState m) -> Double -> m ()+update this item = do+ unless (isNaN item) $ do+ isEmpty <- getIsEmpty this+ if isEmpty+ then do+ writeURef (minValue this) item+ writeURef (maxValue this) item+ else do+ min_ <- minimum this+ max_ <- maximum this+ when (item < min_) $ writeURef (minValue this) item+ when (item > max_) $ writeURef (maxValue this) item+ compactor <- Vector.head <$> getCompactors this+ buff <- Compactor.getBuffer compactor+ DoubleBuffer.append buff item+ modifyURef (retainedItems this) (+1)+ modifyURef (totalN this) (+1)+ modifyURef (sumValue this) (+ item)+ retItems <- getRetainedItems this+ maxNominalCapacity <- getMaxNominalCapacity this+ when (retItems >= maxNominalCapacity) $ do+ DoubleBuffer.sort buff+ compress this+ writeMutVar (aux this) Nothing++-- Private pure bits++getRankLB :: Int -> Int -> Double -> Int -> Bool -> Word64 -> Double+getRankLB k levels rank numStdDev hra totalN = if exactRank k levels rank hra totalN+ then rank+ else max lbRel lbFix+ where+ relative = relRseFactor / fromIntegral k * (if hra then 1.0 - rank else rank)+ fixed = fixRseFactor / fromIntegral k+ lbRel = rank - fromIntegral numStdDev * relative+ lbFix = rank - fromIntegral numStdDev * fixed++getRankUB :: Int -> Int -> Double -> Int -> Bool -> Word64 -> Double+getRankUB k levels rank numStdDev hra totalN = if exactRank k levels rank hra totalN+ then rank+ else min ubRel ubFix+ where+ relative = relRseFactor / fromIntegral k * (if hra then 1.0 - rank else rank)+ fixed = fixRseFactor / fromIntegral k+ ubRel = rank + fromIntegral numStdDev * relative+ ubFix = rank + fromIntegral numStdDev * fixed++exactRank :: Int -> Int -> Double -> Bool -> Word64 -> Bool+exactRank k levels rank hra totalN = (levels == 1 || fromIntegral totalN <= baseCap) || (hra && rank >= 1.0 - exactRankThresh || not hra && rank <= exactRankThresh)+ where+ baseCap = k * initNumberOfSections+ exactRankThresh :: Double+ exactRankThresh = fromIntegral baseCap / fromIntegral totalN++getRSE :: Int -> Double -> Bool -> Word64 -> Double+getRSE k rank hra totalN = getRankUB k 2 rank 1 hra totalN
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Auxiliary.hs view
@@ -0,0 +1,177 @@+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++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+ MUVector.iforM_ weights $ \i weight -> do+ when (i > 0) $ do+ prevWeight <- MUVector.read weights (i - 1)+ MUVector.unsafeWrite weights i (weight + prevWeight)+ lastWeight <- MUVector.read weights (size - 1)+ when (lastWeight /= mraSize this) $ do+ error "invariant violated: lastWeight does not equal raSize"++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)+
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Compactor.hs view
@@ -0,0 +1,237 @@+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), Uniform (uniformM), Gen)+import Control.Exception (assert)+import Control.Monad (when)+import Control.Monad.Trans+import Control.Monad.Primitive+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Constants+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.URef++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))+ }++instance TakeSnapshot ReqCompactor where+ data Snapshot ReqCompactor = ReqCompactorSnapshot+ { snapshotCompactorRankAccuracy :: !RankAccuracy+ , snapshotCompactorRankAccuracyState :: !Word64+ , snapshotCompactorLastFlip :: !Bool+ , snapshotCompactorSectionSizeFlt :: !Double+ , snapshotCompactorSectionSize :: !Word32+ , snapshotCompactorNumSections :: !Word8+ , snapshotCompactorBuffer :: !(Snapshot DoubleBuffer)+ }+ takeSnapshot ReqCompactor{..} = ReqCompactorSnapshot rcRankAccuracy+ <$> readURef rcState+ <*> readURef rcLastFlip+ <*> readURef rcSectionSizeFlt+ <*> readURef rcSectionSize+ <*> readURef rcNumSections+ <*> (readMutVar rcBuffer >>= takeSnapshot)++deriving instance Show (Snapshot ReqCompactor)++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
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/Constants.hs view
@@ -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
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/DoubleBuffer.hs view
@@ -0,0 +1,319 @@+{-# 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+import Control.Monad+import Control.Monad.Primitive+import Control.Monad.Reader.Class+import Data.Primitive.MutVar+import qualified Data.Vector.Unboxed as UVector+import qualified Data.Vector.Unboxed.Mutable as MUVector+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.URef+import Data.Vector.Algorithms.Intro (sortByBounds)+import Data.Vector.Algorithms.Search+import GHC.Prim+import GHC.Stack+import System.IO.Unsafe ()+import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.InequalitySearch as IS+import Control.Exception++-- | 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+ }++instance TakeSnapshot DoubleBuffer where+ data Snapshot DoubleBuffer = DoubleBufferSnapshot+ { dbSnapshotVec :: UVector.Vector Double+ , dbSnapshotCount :: !Int+ , dbSnapshotSorted :: !Bool+ , dbSnapshotGrowthIncrement :: !Int+ , dbSnapshotSpaceAtBottom :: !Bool+ }+ takeSnapshot DoubleBuffer{..} = DoubleBufferSnapshot+ <$> (readMutVar vec >>= UVector.freeze)+ <*> readURef count+ <*> readURef sorted+ <*> pure growthIncrement+ <*> pure spaceAtBottom++deriving instance Show (Snapshot DoubleBuffer)++type Capacity = Int+type GrowthIncrement = Int+type SpaceAtBottom = Bool++-- | Constructs an new empty FloatBuffer with an initial capacity specified by+-- the <code>capacity</code> argument.+mkBuffer :: PrimMonad m => Capacity -> GrowthIncrement -> SpaceAtBottom -> m (DoubleBuffer (PrimState m))+mkBuffer capacity_ growthIncrement spaceAtBottom = do+ vec <- newMutVar =<< MUVector.new capacity_+ 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)
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/InequalitySearch.hs view
@@ -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 <i>{A, B}</i> in the values array+-- instead of just a single value, where <i>A</i> and <i>B</i> are the array indicies of two+-- adjacent values in the array. For all the search criteria, if the algorithm reaches the ends of+-- the search range, the algorithm calls the <i>resolve()</i> method to determine what to+-- return to the caller. If the key value cannot be resolved, it returns a -1 to the caller.++-- Given an array of values <i>arr[]</i> and the search key value <i>v</i>, the algorithms for+-- the searching criteria are as follows:</p>+-- +-- <li><b>LT:</b> Find the highest ranked adjacent pair <i>{A, B}</i> such that:<br>+-- <i>arr[A] < v ≤ arr[B]</i>. The normal return is the index <i>A</i>.+-- </li>+-- <li><b>LE:</b> Find the highest ranked adjacent pair <i>{A, B}</i> such that:<br>+-- <i>arr[A] ≤ v < arr[B]</i>. The normal return is the index <i>A</i>.+-- </li>+-- <li><b>EQ:</b> Find the adjacent pair <i>{A, B}</i> such that:<br>+-- <i>arr[A] ≤ v ≤ arr[B]</i>. The normal return is the index <i>A</i> or <i>B</i> whichever+-- equals <i>v</i>, otherwise it returns -1.+-- </li>+-- <li><b>GE:</b> Find the lowest ranked adjacent pair <i>{A, B}</i> such that:<br>+-- <i>arr[A] < v ≤ arr[B]</i>. The normal return is the index <i>B</i>.+-- </li>+-- <li><b>GT:</b> Find the lowest ranked adjacent pair <i>{A, B}</i> such that:<br>+-- <i>arr[A] ≤ v < arr[B]</i>. The normal return is the index <i>B</i>.+-- </li>+-- </ul>+class InequalitySearch s where+ inequalityCompare :: Ord a+ => s+ -> a + -- ^ V+ -> a+ -- ^ A+ -> a+ -- ^ B+ -> Ordering+ -- ^ 'GT' means we must search higher in the array, 'LT' means we must+ -- search lower in the array, or `EQ`, which means we have found + -- the correct bounding pair.+ getIndex + :: (PrimMonad m, MVector v a, Ord a) + => s + -> v (PrimState m) a + -> Int + -> Int + -> a + -> m Int+ resolve + :: s + -> Int -- Vector length+ -> (Int, Int) + -- ^ Final low index, high index (lo, hi)+ -> (Int, Int) + -- ^ Initial search region (low, high)+ -> Int+ -- ^ A thing++instance InequalitySearch (:<) where+ inequalityCompare _ v a b + | v <= a = LT+ | b < v = GT+ | otherwise = EQ+ getIndex _ _ a _ _ = pure a+ resolve _ vl (lo, hi) (low, high) = if lo >= high then high else vl++instance InequalitySearch (:<=) where+ inequalityCompare _ v a b+ | v < a = LT+ | b <= v = GT+ | otherwise = EQ+ getIndex _ _ a _ _ = pure a+ resolve _ vl (lo, hi) (low, high) = if lo >= high then high else vl ++instance InequalitySearch (:>) where+ inequalityCompare _ v a b+ | v < a = LT+ | b <= v = GT+ | otherwise = EQ+ getIndex _ _ _ b _ = pure b+ resolve _ vl (lo, hi) (low, high) = if hi <= low then low else vl++instance InequalitySearch (:>=) where+ inequalityCompare _ v a b+ | v <= a = LT+ | b < v = GT+ | otherwise = EQ+ getIndex _ _ _ b _ = pure b+ resolve _ vl (lo, hi) (low, high) = if hi <= low then low else vl++find :: (InequalitySearch s, PrimMonad m, MVector v a, Ord a) => s -> v (PrimState m) a -> Int -> Int -> a -> m Int+find strat v low high x = go low (high - 1)+ where+ go lo hi + | lo <= hi && lo < high = do+ let mid = lo + ((hi - lo) `div` 2)+ midV <- MV.read v mid+ midV' <- MV.read v (mid + 1)+ case inequalityCompare strat x midV midV' of+ LT -> go lo (mid - 1)+ EQ -> getIndex strat v mid (mid + 1) x+ GT -> go (mid + 1) hi+ | otherwise = pure $! resolve strat (MV.length v) (lo, hi) (low, high)+{-# INLINE find #-}
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Internal/URef.hs view
@@ -0,0 +1,49 @@+module DataSketches.Quantiles.RelativeErrorQuantile.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:+-- <https://stackoverflow.com/questions/27261813/why-is-my-little-stref-int-require-allocating-gigabytes>.+--+-- 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
+ src/DataSketches/Quantiles/RelativeErrorQuantile/Types.hs view
@@ -0,0 +1,29 @@+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)++class TakeSnapshot a where+ data Snapshot a+ takeSnapshot :: PrimMonad m => a (PrimState m) -> m (Snapshot a)
+ test/AuxiliarySpec.hs view
@@ -0,0 +1,39 @@+module AuxiliarySpec where++import Data.Primitive.MutVar+import qualified Data.Vector.Unboxed.Mutable as MUVector+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary+import qualified DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary as Aux+import DataSketches.Quantiles.RelativeErrorQuantile.Types+import Test.Hspec+import qualified Data.Vector.Unboxed as U+import qualified Data.List+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer+++spec :: Spec+spec = do+ mapM_ checkMergeSortIn [HighRanksAreAccurate, LowRanksAreAccurate]++checkMergeSortIn :: RankAccuracy -> Spec+checkMergeSortIn ra = specify ("mergeSortIn works. hra=" ++ show ra) $ do+ let hraBool = case ra of+ HighRanksAreAccurate -> True+ LowRanksAreAccurate -> False+ let oddItems = [1,3..11] -- 6 items+ let evenItems = [2,4..12]++ buf1 <- mkBuffer 25 0 hraBool+ mapM_ (append buf1) oddItems+ + buf2 <- mkBuffer 25 0 hraBool+ mapM_ (append buf2) evenItems++ let n = 12+ weightedItems <- newMutVar =<< MUVector.new 25 + let aux = MReqAuxiliary weightedItems ra n+ Aux.mergeSortIn aux buf1 1 0+ Aux.mergeSortIn aux buf2 2 6+ (items, _) <- fmap U.unzip . U.freeze =<< readMutVar (mraWeightedItems aux)+ let itemsToCheck = U.slice 0 (fromIntegral n) items+ itemsToCheck `shouldBe` U.fromList (Data.List.sort (oddItems ++ evenItems))
+ test/CompactorSpec.hs view
@@ -0,0 +1,9 @@+module CompactorSpec where++import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Compactor+import Test.Hspec++spec :: Spec+spec = do+ specify "nearestEven behaves as expected" $ do+ nearestEven (-0.9) `shouldBe` 0
+ test/DoubleBufferSpec.hs view
@@ -0,0 +1,116 @@+module DoubleBufferSpec where++import Control.Monad+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer+import DataSketches.Quantiles.RelativeErrorQuantile.Types+import qualified Data.Vector.Unboxed as UV+import Test.Hspec+import qualified Data.List++spec :: Spec+spec = do+ describe "Merge Sort" $ do+ let mergeSortTestImpl direction = do+ buf <- mkBuffer 16 0 direction+ forM [0..31] (append buf)+ buf2 <- copyBuffer buf+ mergeSortIn buf buf2+ v <- getVector buf+ UV.freeze v+ it "works upwards" $ do+ v <- mergeSortTestImpl False+ v `shouldBe` UV.fromList (Data.List.sort ([0..31] ++ [0..31]))+ it "works downwards" $ do+ v <- mergeSortTestImpl True+ v `shouldBe` UV.fromList (Data.List.sort ([0..31] ++ [0..31]))+ describe "getEvensOrOdds" $ do+ checkGetEvensOrOdds False False+ checkGetEvensOrOdds False True+ checkGetEvensOrOdds True False+ checkGetEvensOrOdds True True+ describe "checkAppendAndSpaceTop" $ do+ checkAppendAndSpace True+ checkAppendAndSpace False+ describe "checkEnsureCapacity" $ do+ mapM_ checkEnsureCapacity [True, False]+ describe "checkCountLessThanImpl" $ do+ mapM_ checkCountLessThanImpl [True, False]++checkGetEvensOrOdds :: Bool -> Bool -> Spec+checkGetEvensOrOdds odds spaceAtBottom = + it ("works for odds=" ++ show odds ++ " spaceAtBottom=" ++ show spaceAtBottom) $ do+ buf <- mkBuffer cap 0 spaceAtBottom+ forM_ [0..cap `div` 2] (append buf . fromIntegral)+ out <- getEvensOrOdds buf 0 (cap `div` 2) odds+ v <- UV.freeze =<< getVector out+ v `shouldSatisfy` UV.all (\x -> if odds then floor x `mod` 2 == 1 else floor x `mod` 2 == 0)+ where+ cap = 16++checkAppendAndSpace :: Bool -> Spec+checkAppendAndSpace spaceAtBottom = it ("works for spaceAtBottom=" ++ show spaceAtBottom) $ do+ buf <- mkBuffer 2 2 spaceAtBottom+ do+ count_ <- getCount buf+ count_ `shouldBe` 0+ capacity_ <- getCapacity buf + capacity_ `shouldBe` 2+ space_ <- getSpace buf+ space_ `shouldBe` 2+ append buf 1+ do+ count_ <- getCount buf+ count_ `shouldBe` 1+ capacity_ <- getCapacity buf + capacity_ `shouldBe` 2+ space_ <- getSpace buf+ space_ `shouldBe` 1+ append buf 2+ do+ count_ <- getCount buf+ count_ `shouldBe` 2+ capacity_ <- getCapacity buf + capacity_ `shouldBe` 2+ space_ <- getSpace buf+ space_ `shouldBe` 0+ append buf 3+ do+ count_ <- getCount buf+ count_ `shouldBe` 3+ capacity_ <- getCapacity buf + capacity_ `shouldBe` 5+ space_ <- getSpace buf+ space_ `shouldBe` 2++checkEnsureCapacity :: Bool -> Spec+checkEnsureCapacity spaceAtBottom = it ("works for spaceAtBottom=" ++ show spaceAtBottom) $ do+ buf <- mkBuffer 4 2 spaceAtBottom+ append buf 2+ append buf 1+ append buf 3+ ensureCapacity buf 8+ sort buf+ x1 <- buf ! 0+ x1 `shouldBe` 1+ x2 <- buf ! 1+ x2 `shouldBe` 2+ x3 <- buf ! 2+ x3 `shouldBe` 3++checkCountLessThanImpl :: Bool -> Spec+checkCountLessThanImpl spaceAtBottom = it ("works for spaceAtBottom=" ++ show spaceAtBottom) $ do + buf <- mkBuffer 7 0 spaceAtBottom+ mapM_ (append buf) [1..7]+ sort buf+ buf2 <- mkBuffer 7 0 spaceAtBottom+ mergeSortIn buf2 buf+ do+ count_ <- getCountWithCriterion buf2 4 (:<)+ count_ `shouldBe` 3+ mergeSortIn buf2 buf+ do+ count_ <- getCountWithCriterion buf2 4 (:<)+ count_ `shouldBe` 6+ do+ count_ <- getCount buf2+ count_ `shouldBe` 14
+ test/ProofCheckSpec.hs view
@@ -0,0 +1,80 @@+{-# LANGUAGE TypeApplications #-}+module ProofCheckSpec where++import Statistics.Quantile+import Test.Hspec+import qualified Data.Vector as V+import DataSketches.Quantiles.RelativeErrorQuantile (mkReqSketch, RankAccuracy (HighRanksAreAccurate), update)+import qualified DataSketches.Quantiles.RelativeErrorQuantile as SK+import Test.QuickCheck+import Control.Monad.ST+import Control.Monad.Primitive+import GHC.TypeLits+import Control.Monad+import Debug.Trace++sampleData :: V.Vector Double+sampleData = V.fromList [1..200]++{-+compareRealToApproximate :: IO ()+compareRealToApproximate = do+ let realQuantiles = quantiles spss [50, 90, 95, 99] 100 sampleData+ -- print realQuantiles+ sk <- mkReqSketch @6 HighRanksAreAccurate+ mapM_ (update sk) sampleData+ let ranks = [0.01, 0.02 .. 0.99]+ rankInts = map (floor . (* 100)) ranks+ -- print =<< SK.quantiles sk ranks+ upperLowerBounds <- mapM (upperAndLowerBound sk) ranks+ forM_ upperAndLowerBound $ \(actualRank, l, u) -> do+ let actual = quantile spss actualRank 100 sampleData+ assert (actual >= l && actual <= u)+-}++upperAndLowerBound :: (PrimMonad m, KnownNat k) => SK.ReqSketch k (Control.Monad.Primitive.PrimState m)+ -> Double -> m (Double, Double, Double)+upperAndLowerBound sk r = do+ l <- SK.rankLowerBound sk r 3+ u <- SK.rankUpperBound sk r 3+ pure (r, l, u)++spec :: Spec+spec = do+ let sampleInputGen = arbitrary `suchThat` (all (\x -> not (isNaN x) && not (isInfinite x)) . getNonEmpty)++ specify "quantile ranks should fall within advertised upper and lower bounds" $+ property $ forAll sampleInputGen $ \(NonEmpty doubles) -> runST $ do+ let sampleData = V.fromList doubles+ sk <- mkReqSketch @6 HighRanksAreAccurate+ mapM_ (update sk) sampleData+ let ranks = [0.01, 0.02 .. 0.99]+ rankInts = map (floor . (* 100)) ranks+ upperLowerBounds <- mapM (upperAndLowerBound sk) ranks+ ranksInAdvertisedRanges <- forM upperLowerBounds $ \(actualRank, l, u) -> do+ pure (actualRank >= l && actualRank <= u)+ pure $! and ranksInAdvertisedRanges++ {- TODO+ specify "values at quantiles should be close to real quantile" $+ property $ forAll sampleInputGen $ \(NonEmpty doubles) -> runST $ do+ let sampleData = V.fromList doubles+ sk <- mkReqSketch @6 HighRanksAreAccurate+ mapM_ (update sk) sampleData+ let ranks = [0.01, 0.02 .. 0.99]+ rankInts = map (floor . (* 100)) ranks+ upperLowerBounds <- mapM (upperAndLowerBound sk) ranks+ valuesInAdvertisedRanges <- forM upperLowerBounds $ \(actualRank, l, u) -> do+ let actual = quantile spss (floor (actualRank * 100)) 100 sampleData+ estimatedL <- SK.quantile sk l+ estimatedH <- SK.quantile sk u+ traceShow (estimatedL, actual, estimatedH) $ pure+ ((actual > estimatedL || actual `approximatelyEqual` estimatedL) && (actual < estimatedH || actual `approximatelyEqual` estimatedH))+ pure $! and valuesInAdvertisedRanges+ -}++approximatelyEqual :: Double -> Double -> Bool+approximatelyEqual x y = (abs (x - y) / max (abs x) (abs y)) < eEqD++eEqD :: Double+eEqD = 1e-5
+ test/RelativeErrorQuantileSpec.hs view
@@ -0,0 +1,261 @@+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE NumericUnderscores #-}+module RelativeErrorQuantileSpec where++import Data.Primitive.MutVar+import qualified Data.Vector.Unboxed as U+import Control.Monad+import Control.Monad.Primitive+import DataSketches.Quantiles.RelativeErrorQuantile+import DataSketches.Quantiles.RelativeErrorQuantile.Types+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.Auxiliary+import Data.List+import Data.Maybe (fromJust, isJust)+import Data.Proxy+import GHC.TypeLits+import Test.Hspec+import DataSketches.Quantiles.RelativeErrorQuantile.Internal.DoubleBuffer (DoubleIsNonFiniteException(..))+import Text.Show.Pretty++spec :: Spec+spec = do+ specify "non finite PMF/CDF should throw" $ asIO $ do+ sk <- mkReqSketch @6 HighRanksAreAccurate+ update sk 1+ cumulativeDistributionFunction sk [0 / 0] `shouldThrow` (== CumulativeDistributionInvariantsSplitsAreNotFinite)+ specify "updating a sketch with NaN should ignore it" $ asIO $ do+ sk <- mkReqSketch @6 HighRanksAreAccurate+ update sk (0 / 0)+ isEmpty <- DataSketches.Quantiles.RelativeErrorQuantile.null sk+ isEmpty `shouldBe` True+ specify "non finite rank should throw" $ asIO $ do+ let infinity = read "Infinity"::Double+ sk <- mkReqSketch @6 HighRanksAreAccurate+ update sk 1+ (rank sk infinity >>= print) `shouldThrow` (\(DoubleIsNonFiniteException _) -> True)+ specify "big merge doesn't explode" $ asIO $ do+ sk1 <- mkReqSketch @6 HighRanksAreAccurate+ mapM_ (update sk1) [5..10]+ sk2 <- mkReqSketch @6 HighRanksAreAccurate+ merge sk1 sk2+ mapM_ (update sk2) [1..15]+ merge sk1 sk2+ n <- getN sk1+ n `shouldBe` 21+ mapM_ (update sk2) [16..300]+ merge sk1 sk2+ n <- getN sk1+ n `shouldBe` 321+ describe "property tests" $ do+ specify "ReqSketch quantile estimates are within ε bounds compared to real quantile calculations" $ print ()+ specify "merging N ReqSketches is equivalent +/- ε to inserting the same values into 1 ReqSketch" $ print ()++ let simpleTestValues = [5, 5, 5, 6, 6, 6, 7, 8, 8, 8]+ let lessThanRs = [0.0, 0.0, 0.0, 0.3, 0.3, 0.3, 0.6, 0.7, 0.7, 0.7]+ let lessThanEqRs = [0.3, 0.3, 0.3, 0.6, 0.6, 0.6, 0.7, 1.0, 1.0, 1.0]+ let simpleTestSetup = do+ sk <- mkReqSketch @50 HighRanksAreAccurate+ mapM_ (update sk) simpleTestValues+ pure sk+ specify "lots of repeat values should work" $ asIO $ do+ sk <- mkReqSketch @6 HighRanksAreAccurate+ replicateM_ 10_000 (update sk 1 >> getRetainedItems sk)+ print =<< getRetainedItems sk++ describe "simple test" $ before simpleTestSetup $ do+ describe "<" $ do+ specify "ranks function should match lessThanRs" $ \sk -> do+ actualRanks <- ranks sk simpleTestValues+ actualRanks `shouldBe` lessThanRs+ specify "mapM rank should match ranks behaviour" $ \sk -> do+ actualRanks <- mapM (rank sk) simpleTestValues+ actualRanks `shouldBe` lessThanRs+ describe "<=" $ do+ let mkSk' sk = sk { criterion = (:<=) } :: ReqSketch 50 (PrimState IO)+ specify "ranks function should match lessThanRs" $ \sk -> do+ let sk' = mkSk' sk+ actualRanks <- ranks sk' simpleTestValues+ actualRanks `shouldBe` lessThanEqRs+ specify "mapM rank should match ranks behaviour" $ \sk -> do+ let sk' = mkSk' sk+ actualRanks <- mapM (rank sk') simpleTestValues+ actualRanks `shouldBe` lessThanEqRs++++++++ -- k min max hra lteq low-to-high or high-to-low+ bigTest Proxy 1 200 HighRanksAreAccurate (:<=) True+ bigTest Proxy 1 200 LowRanksAreAccurate (:<=) True+ bigTest Proxy 1 200 HighRanksAreAccurate (:<) False+ bigTest Proxy 1 200 LowRanksAreAccurate (:<) True++ mergeSpec++ describe "quantiles" $ do+ it "should be reasonable" $ asIO $ do+ sk <- mkReqSketch @6 HighRanksAreAccurate+ update sk 1+ update sk 1+ update sk 1+ r <- quantile sk 0.5 + r `shouldBe` 1.0+++bigTest :: Proxy 6 -> Int -> Int -> RankAccuracy -> Criterion -> Bool -> Spec+bigTest k min_ max_ hra crit up = do+ let testName = unwords+ [ "k=" <> show (natVal k)+ , "min=" <> show min_+ , "max=" <> show max_+ , "hra=" <> show hra+ , "up=" <> show up+ ]+ testContents :: IO ()+ testContents = do+ sk <- loadSketch k min_ max_ hra crit up+ checkAux sk+ checkGetRank sk min_ max_+ checkGetRanks sk max_+ checkGetQuantiles sk+ checkGetCDF sk+ checkGetPMF sk+ -- checkIterator sk+ -- checkMerge sk+ it testName testContents++asIO :: IO a -> IO a+asIO = id++mergeSpec :: Spec+mergeSpec = specify "merge works" $ asIO $ do+ s1 <- mkReqSketch HighRanksAreAccurate :: IO (ReqSketch 12 (PrimState IO))+ mapM_ (update s1) [0..40]+ s2 <- mkReqSketch HighRanksAreAccurate :: IO (ReqSketch 12 (PrimState IO))+ mapM_ (update s2) [0..40]+ s3 <- mkReqSketch HighRanksAreAccurate :: IO (ReqSketch 12 (PrimState IO))+ mapM_ (update s3) [0..40]+ s <- mkReqSketch HighRanksAreAccurate :: IO (ReqSketch 12 (PrimState IO))+ s `merge` s1+ s `merge` s2+ s `merge` s3+ pure ()++loadSketch :: forall n. (KnownNat n, ValidK n) => Proxy n -> Int -> Int -> RankAccuracy -> Criterion -> Bool -> IO (ReqSketch n (PrimState IO))+loadSketch k min_ max_ hra ltEq up = do+ sk <- mkReqSketch hra :: IO (ReqSketch n (PrimState IO))+ -- This just seems geared at making sure that ranks come out right regardless of order+ mapM_ (update sk . fromIntegral) $ if up+ then [min_ .. max_]+ else reverse [min_ .. max_ {- + 1 -}]+ pure sk++checkAux :: ReqSketch n (PrimState IO) -> IO ()+checkAux sk = do+ auxiliary <- mkAuxiliaryFromReqSketch sk+ totalCount <- computeTotalRetainedItems sk++ let rows = raWeightedItems auxiliary+ getRow = (U.!)+ let initialRow = getRow rows 0+ otherRows = map (getRow rows) [1..totalCount - 1]+ foldM_+ (\lastRow thisRow -> do+ fst thisRow `shouldSatisfy` (>= fst lastRow)+ snd thisRow `shouldSatisfy` (>= snd lastRow)+ pure thisRow+ )+ initialRow+ otherRows++checkGetRank :: KnownNat n => ReqSketch n (PrimState IO) -> Int -> Int -> IO ()+checkGetRank sk min_ max_ = do+ let (v : spArr) = evenlySpacedFloats 0 (fromIntegral max_) 11+ initialRank <- rank sk v+ foldM_+ (\(oldV, oldRank) v -> do+ r <- rank sk v+ v `shouldSatisfy` (>= oldV)+ r `shouldSatisfy` (>= oldRank)+ pure (v, r)+ )+ (v, initialRank)+ spArr++checkGetRanks :: KnownNat n => ReqSketch n (PrimState IO) -> Int -> IO ()+checkGetRanks sk max_ = do+ let sp = evenlySpacedFloats 0 (fromIntegral max_) 11+ void $ ranks sk sp++checkGetCDF :: KnownNat n => ReqSketch n (PrimState IO) -> IO ()+checkGetCDF sk = do+ let spArr = [20, 40 .. 180]+ r <- cumulativeDistributionFunction sk spArr+ r `shouldSatisfy` isJust++checkGetPMF :: KnownNat n => ReqSketch n (PrimState IO) -> IO ()+checkGetPMF sk = do+ let spArr = [20, 40 .. 180]+ r <- probabilityMassFunction sk spArr+ r `shouldNotSatisfy` Data.List.null++{-+checkMerge :: ReqSketch n (PrimState IO) -> IO ()+checkMerge sk = do+ sk' <- copyRs+-}++checkGetRankConcreteExample :: RankAccuracy -> Criterion -> IO ()+checkGetRankConcreteExample ra crit = do+ sk <- loadSketch (Proxy :: Proxy 12) 1 1000 ra crit True+ rLB <- rankLowerBound sk 0.5 1+ rLB `shouldSatisfy` (> 0)+ rLB <- case ra of+ HighRanksAreAccurate -> rankLowerBound sk (995 / 1000) 1+ LowRanksAreAccurate -> rankLowerBound sk (5 / 1000) 1+ rLB `shouldSatisfy` (> 0)+ rUB <- rankUpperBound sk 0.5 1+ rUB `shouldSatisfy` (> 0)+ rUB <- case ra of+ HighRanksAreAccurate -> rankUpperBound sk (995 / 1000) 1+ LowRanksAreAccurate -> rankUpperBound sk (5 / 1000) 1+ rUB `shouldSatisfy` (> 0)+ void $ ranks sk [5, 100]++++++checkGetQuantiles+ :: KnownNat n+ => ReqSketch n (PrimState IO)+ -> IO ()+checkGetQuantiles sk = do+ let rArr = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]+ -- nothing getting checked here apparently, guess the thing not+ -- exploding is sufficient.+ qOut <- quantiles sk rArr+ pure ()++-- | Returns a float array of evenly spaced values between value1 and value2 inclusive.+-- If value2 > value1, the resulting sequence will be increasing.+-- If value2 < value1, the resulting sequence will be decreasing.+-- value1 will be in index 0 of the returned array+-- value2 will be in the highest index of the returned array+-- valu3 is the total number of values including value1 and value2. Must be 2 or greater.+-- returns a float array of evenly spaced values between value1 and value2 inclusive.+evenlySpacedFloats :: Double -> Double -> Word -> [Double]+evenlySpacedFloats _ _ 0 = error "Needs at least two steps"+evenlySpacedFloats _ _ 1 = error "Needs at least two steps"+evenlySpacedFloats value1 value2 steps = unfoldr+ (\ix -> let val = fromIntegral ix * delta + value1 in case value1 `compare` value2 of+ LT -> if val > value2 then Nothing else Just (val, ix + 1)+ EQ -> Nothing+ GT -> if val < value2 then Nothing else Just (val, ix + 1)+ )+ 0+ where+ delta = (value2 - value1) / (fromIntegral steps - 1)
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}