packages feed

tdigest (empty) → 0

raw patch · 10 files changed

+1381/−0 lines, 10 filesdep +Chartdep +Chart-diagramsdep +basebuild-type:Customsetup-changed

Dependencies added: Chart, Chart-diagrams, base, base-compat, binary, bytes, deepseq, directory, doctest, filepath, machines, mwc-random, optparse-applicative, parallel, reducers, semigroups, statistics, tasty, tasty-quickcheck, tdigest, time, vector, vector-algorithms

Files

+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2016 Futurice Oy++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 Oleg Grenrus 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.
+ README.md view
@@ -0,0 +1,32 @@+# tdigest++A new data structure for accurate on-line accumulation of rank-based statistics such as quantiles and trimmed means.++See original paper: ["Computing extremely accurate quantiles using t-digest"](https://github.com/tdunning/t-digest/blob/master/docs/t-digest-paper/histo.pdf) by Ted Dunning and Otmar Ertl++## Synopsis++```hs+λ *Data.TDigest > median (tdigest [1..1000] :: TDigest 3)+Just 499.0090729817737+```++## Benchmarks+++Using 50M exponentially distributed numbers:++- average: **16s**; incorrect approximation of median, mostly to measure prng speed+- sorting using `vector-algorithms`: **33s**; using 1000MB of memory+- sparking t-digest (using some `par`): **53s**+- buffered t-digest: **68s**+- sequential t-digest: **65s**++### Example histogram++```+tdigest-simple -m tdigest -d standard -s 100000 -c 10 -o output.svg -i 34+inkscape --export-png=example.png --export-dpi=80 --export-background-opacity=0 --without-gui example.svg+```++![Example](https://raw.githubusercontent.com/futurice/haskell-tdigest/master/example.png)
+ Setup.lhs view
@@ -0,0 +1,182 @@+\begin{code}+{-# LANGUAGE CPP #-}+{-# LANGUAGE OverloadedStrings #-}+module Main (main) where++#ifndef MIN_VERSION_cabal_doctest+#define MIN_VERSION_cabal_doctest(x,y,z) 0+#endif+++#if MIN_VERSION_cabal_doctest(1,0,0)+import Distribution.Extra.Doctest ( defaultMainWithDoctests )+#else++-- Otherwise we provide a shim++#ifndef MIN_VERSION_Cabal+#define MIN_VERSION_Cabal(x,y,z) 0+#endif+#ifndef MIN_VERSION_directory+#define MIN_VERSION_directory(x,y,z) 0+#endif+#if MIN_VERSION_Cabal(1,24,0)+#define InstalledPackageId UnitId+#endif++import Control.Monad ( when )+import Data.List ( nub )+import Data.String ( fromString )+import Distribution.Package ( InstalledPackageId )+import Distribution.Package ( PackageId, Package (..), packageVersion )+import Distribution.PackageDescription ( PackageDescription(), TestSuite(..) , Library (..), BuildInfo (..))+import Distribution.Simple ( defaultMainWithHooks, UserHooks(..), simpleUserHooks )+import Distribution.Simple.Utils ( rewriteFile, createDirectoryIfMissingVerbose )+import Distribution.Simple.BuildPaths ( autogenModulesDir )+import Distribution.Simple.Setup ( BuildFlags(buildDistPref, buildVerbosity), fromFlag)+import Distribution.Simple.LocalBuildInfo ( withPackageDB, withLibLBI, withTestLBI, LocalBuildInfo(), ComponentLocalBuildInfo(componentPackageDeps), compiler )+import Distribution.Simple.Compiler ( showCompilerId , PackageDB (..))+import Distribution.Text ( display , simpleParse )+import System.FilePath ( (</>) )++#if MIN_VERSION_Cabal(1,25,0)+import Distribution.Simple.BuildPaths ( autogenComponentModulesDir )+#endif++#if MIN_VERSION_directory(1,2,2)+import System.Directory (makeAbsolute)+#else+import System.Directory (getCurrentDirectory)+import System.FilePath (isAbsolute)++makeAbsolute :: FilePath -> IO FilePath+makeAbsolute p | isAbsolute p = return p+               | otherwise    = do+    cwd <- getCurrentDirectory+    return $ cwd </> p+#endif++generateBuildModule :: String -> BuildFlags -> PackageDescription -> LocalBuildInfo -> IO ()+generateBuildModule testsuiteName flags pkg lbi = do+  let verbosity = fromFlag (buildVerbosity flags)+  let distPref = fromFlag (buildDistPref flags)++  -- Package DBs+  let dbStack = withPackageDB lbi ++ [ SpecificPackageDB $ distPref </> "package.conf.inplace" ]+  let dbFlags = "-hide-all-packages" : packageDbArgs dbStack++  withLibLBI pkg lbi $ \lib libcfg -> do+    let libBI = libBuildInfo lib++    -- modules+    let modules = exposedModules lib ++ otherModules libBI+    -- it seems that doctest is happy to take in module names, not actual files!+    let module_sources = modules++    -- We need the directory with library's cabal_macros.h!+#if MIN_VERSION_Cabal(1,25,0)+    let libAutogenDir = autogenComponentModulesDir lbi libcfg+#else+    let libAutogenDir = autogenModulesDir lbi+#endif++    -- Lib sources and includes+    iArgs <- mapM (fmap ("-i"++) . makeAbsolute) $ libAutogenDir : hsSourceDirs libBI+    includeArgs <- mapM (fmap ("-I"++) . makeAbsolute) $ includeDirs libBI++    -- CPP includes, i.e. include cabal_macros.h+    let cppFlags = map ("-optP"++) $+            [ "-include", libAutogenDir ++ "/cabal_macros.h" ]+            ++ cppOptions libBI++    withTestLBI pkg lbi $ \suite suitecfg -> when (testName suite == fromString testsuiteName) $ do++      -- get and create autogen dir+#if MIN_VERSION_Cabal(1,25,0)+      let testAutogenDir = autogenComponentModulesDir lbi suitecfg+#else+      let testAutogenDir = autogenModulesDir lbi+#endif+      createDirectoryIfMissingVerbose verbosity True testAutogenDir++      -- write autogen'd file+      rewriteFile (testAutogenDir </> "Build_doctests.hs") $ unlines+        [ "module Build_doctests where"+        , ""+        -- -package-id etc. flags+        , "pkgs :: [String]"+        , "pkgs = " ++ (show $ formatDeps $ testDeps libcfg suitecfg)+        , ""+        , "flags :: [String]"+        , "flags = " ++ show (iArgs ++ includeArgs ++ dbFlags ++ cppFlags)+        , ""+        , "module_sources :: [String]"+        , "module_sources = " ++ show (map display module_sources)+        ]+  where+    -- we do this check in Setup, as then doctests don't need to depend on Cabal+    isOldCompiler = maybe False id $ do+      a <- simpleParse $ showCompilerId $ compiler lbi+      b <- simpleParse "7.5"+      return $ packageVersion (a :: PackageId) < b++    formatDeps = map formatOne+    formatOne (installedPkgId, pkgId)+      -- The problem is how different cabal executables handle package databases+      -- when doctests depend on the library+      | packageId pkg == pkgId = "-package=" ++ display pkgId+      | otherwise              = "-package-id=" ++ display installedPkgId++    -- From Distribution.Simple.Program.GHC+    packageDbArgs :: [PackageDB] -> [String]+    packageDbArgs | isOldCompiler = packageDbArgsConf+                  | otherwise     = packageDbArgsDb++    -- GHC <7.6 uses '-package-conf' instead of '-package-db'.+    packageDbArgsConf :: [PackageDB] -> [String]+    packageDbArgsConf dbstack = case dbstack of+      (GlobalPackageDB:UserPackageDB:dbs) -> concatMap specific dbs+      (GlobalPackageDB:dbs)               -> ("-no-user-package-conf")+                                           : concatMap specific dbs+      _ -> ierror+      where+        specific (SpecificPackageDB db) = [ "-package-conf=" ++ db ]+        specific _                      = ierror+        ierror = error $ "internal error: unexpected package db stack: "+                      ++ show dbstack++    -- GHC >= 7.6 uses the '-package-db' flag. See+    -- https://ghc.haskell.org/trac/ghc/ticket/5977.+    packageDbArgsDb :: [PackageDB] -> [String]+    -- special cases to make arguments prettier in common scenarios+    packageDbArgsDb dbstack = case dbstack of+      (GlobalPackageDB:UserPackageDB:dbs)+        | all isSpecific dbs              -> concatMap single dbs+      (GlobalPackageDB:dbs)+        | all isSpecific dbs              -> "-no-user-package-db"+                                           : concatMap single dbs+      dbs                                 -> "-clear-package-db"+                                           : concatMap single dbs+     where+       single (SpecificPackageDB db) = [ "-package-db=" ++ db ]+       single GlobalPackageDB        = [ "-global-package-db" ]+       single UserPackageDB          = [ "-user-package-db" ]+       isSpecific (SpecificPackageDB _) = True+       isSpecific _                     = False++testDeps :: ComponentLocalBuildInfo -> ComponentLocalBuildInfo -> [(InstalledPackageId, PackageId)]+testDeps xs ys = nub $ componentPackageDeps xs ++ componentPackageDeps ys++defaultMainWithDoctests :: String -> IO ()+defaultMainWithDoctests testSuiteName = defaultMainWithHooks simpleUserHooks+  { buildHook = \pkg lbi hooks flags -> do+     generateBuildModule testSuiteName flags pkg lbi+     buildHook simpleUserHooks pkg lbi hooks flags+  }++#endif++main :: IO ()+main = defaultMainWithDoctests "doctests"++\end{code}
+ bench/Simple.hs view
@@ -0,0 +1,296 @@+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE RankNTypes          #-}+{-# LANGUAGE ScopedTypeVariables #-}+module Main (main) where++import Prelude ()+import Prelude.Compat+import Control.Monad               (join, replicateM)+import Control.Monad.ST            (runST)+import Control.Parallel.Strategies (parList, rseq, using)+import Data.Foldable               (for_)+import Data.List                   (sort)+import Data.Machine+import Data.Machine.Runner         (runT1)+import Data.Monoid                 ((<>))+import Data.Proxy                  (Proxy (..))+import Data.Time                   (diffUTCTime, getCurrentTime)+import Data.Word                   (Word32)+import GHC.TypeLits                (KnownNat, SomeNat (..), someNatVal)+import Statistics.Distribution     (ContGen (..), density)++import Statistics.Distribution.Exponential (exponential)+import Statistics.Distribution.Gamma       (gammaDistr)+import Statistics.Distribution.Normal      (standard)+import Statistics.Distribution.Uniform     (uniformDistr)++import qualified Data.Vector.Algorithms.Intro as Intro+import qualified Data.Vector.Unboxed          as V+import qualified Data.Vector.Unboxed.Mutable  as VU+import qualified Options.Applicative          as O+import qualified System.Random.MWC            as MWC++import qualified Graphics.Rendering.Chart.Backend.Diagrams as Chart+import           Graphics.Rendering.Chart.Easy             ((&), (.~), (^.))+import qualified Graphics.Rendering.Chart.Easy             as Chart++import Data.TDigest++-------------------------------------------------------------------------------+-- Data+-------------------------------------------------------------------------------++data Method+    = MethodAverage+    | MethodNaive+    | MethodVector+    | MethodTDigest+    | MethodTDigestBuffered+    | MethodTDigestSparking+  deriving (Show)++data Distrib+    = DistribIncr+    | DistribUniform+    | DistribExponent+    | DistribGamma+    | DistribStandard+  deriving (Show)++timed :: Show a => IO a -> IO ()+timed mx = do+    s <- getCurrentTime+    x <- mx+    print x+    e <- getCurrentTime+    print (diffUTCTime e s)++action :: Method -> Distrib -> Int -> Int -> Word32 -> Maybe FilePath -> IO ()+action m d s c iseed fp = do+    print (m, d, s, c)+    let seed = initSeed (V.singleton iseed)+    let dens = case d of+            DistribIncr     -> density $ uniformDistr 0 (fromIntegral s)+            DistribUniform  -> density $ uniformDistr 0 1     -- median around 0.5+            DistribExponent -> density $ exponential $ log 2  -- median around 1.0+            DistribGamma    -> density $ gammaDistr 0.1 0.1   -- median around .0000593391+            DistribStandard -> density standard+    let input = take s $ case d of+            DistribIncr     -> [1 .. fromIntegral s] -- not sure, but end point prevents floating+            DistribUniform  -> randomStream (uniformDistr 0 1) seed     -- median around 0.5+            DistribExponent -> randomStream (exponential $ log 2) seed  -- median around 1.0+            DistribGamma    -> randomStream (gammaDistr 0.1 0.1) seed   -- median around .0000593391+            DistribStandard -> randomStream standard seed+    let method = case m of+          MethodAverage         -> pure . average+          MethodNaive           -> pure . naiveMedian+          MethodVector          -> pure . vectorMedian+          MethodTDigest         -> reifyNat c $ tdigestMachine fp dens+          MethodTDigestBuffered -> reifyNat c $ tdigestBufferedMachine fp dens+          MethodTDigestSparking -> reifyNat c $ tdigestSparkingMachine fp dens+    timed $ method input++reifyNat :: forall x. Int -> (forall n. KnownNat n => Proxy n -> x) -> x+reifyNat n f = case someNatVal (fromIntegral n) of+    Nothing           -> error "Negative m"+    Just (SomeNat cp) -> f cp++actionParser :: O.Parser (IO ())+actionParser = action+    <$> O.option (maybeReader readMethod) (+        O.short 'm' <> O.long "method" <> O.metavar ":method" <> O.value MethodTDigestBuffered)+    <*> O.option (maybeReader readDistrib) (+        O.short 'd' <> O.long "distrib" <> O.metavar ":distrib" <> O.value DistribUniform)+    <*> O.option O.auto (+        O.short 's' <> O.long "size" <> O.metavar ":size" <> O.value 1000000)+    <*> O.option O.auto (+        O.short 'c' <> O.long "compression" <> O.metavar ":comp" <> O.value 20)+    <*> O.option O.auto (+        O.short 'i' <> O.long "seed" <> O.metavar ":seed" <> O.value 42)+    <*> O.optional (O.strOption (+        O.short 'o' <> O.long "output" <> O.metavar ":output.svg"))+  where+    readMethod "average"  = Just MethodAverage+    readMethod "naive"    = Just MethodNaive+    readMethod "vector"   = Just MethodVector+    readMethod "digest"   = Just MethodTDigest+    readMethod "tdigest"  = Just MethodTDigest+    readMethod "buffered" = Just MethodTDigestBuffered+    readMethod "sparking" = Just MethodTDigestSparking+    readMethod _          = Nothing++    readDistrib "incr"     = Just DistribIncr+    readDistrib "uniform"  = Just DistribUniform+    readDistrib "exponent" = Just DistribExponent+    readDistrib "standard" = Just DistribStandard+    readDistrib "gamma"    = Just DistribGamma+    readDistrib _          = Nothing++-- Only on optparse-applicative-0.13+maybeReader :: (String -> Maybe a) -> O.ReadM a+maybeReader f = O.eitherReader $ \x -> maybe (Left x) Right (f x)++main :: IO ()+main = join (O.execParser opts)+  where+    opts = O.info (O.helper <*> actionParser)+        (O.fullDesc <> O.header "tdigest-simple - a small utility to explore tdigest")++-------------------------------------------------------------------------------+-- Methods+-------------------------------------------------------------------------------++average :: [Double] -> Maybe Double+average []     = Nothing+average (x:xs) = Just $ go x 1 xs+  where+    go z _ []       = z+    go z n (y : ys) = go ((z * n + y) / (n + 1)) (n + 1) ys++naiveMedian :: [Double] -> Maybe Double+naiveMedian [] = Nothing+naiveMedian xs = Just $ sort xs !! (length xs `div` 2)++vectorMedian :: [Double] -> Maybe Double+vectorMedian l+    | null l    = Nothing+    | otherwise = runST $ do+        let v = V.fromList l+        mv <- V.thaw v+        Intro.sort mv+        Just <$> VU.unsafeRead mv (VU.length mv `div` 2)++tdigestMachine+    :: forall comp. KnownNat comp+    => Maybe FilePath -> (Double -> Double) -> Proxy comp -> [Double] -> IO (Maybe Double)+tdigestMachine fp dens _ input = do+    mdigest <- fmap validate <$> runT1 machine+    case mdigest of+        Nothing             -> return Nothing+        Just (Left err)     -> fail $ "Validation error: " ++ err+        Just (Right digest) -> do+            printStats fp dens digest+            return $ median digest+  where+    machine :: MachineT IO k (TDigest comp)+    machine+        =  fold (flip insert) mempty+        <~ source input++tdigestBufferedMachine+    :: forall comp. KnownNat comp+    => Maybe FilePath -> (Double -> Double) -> Proxy comp -> [Double] -> IO (Maybe Double)+tdigestBufferedMachine fp dens _ input = do+    mdigest <- fmap validate <$> runT1 machine+    case mdigest of+        Nothing             -> return Nothing+        Just (Left err)     -> fail $ "Validation error: " ++ err+        Just (Right digest) -> do+            printStats fp dens digest+            return $ median digest+  where+    machine :: MachineT IO k (TDigest comp)+    machine+        =  fold mappend mempty+        <~ mapping tdigest+        <~ buffered 10000+        <~ source input++-- Sparking machine doesn't count+tdigestSparkingMachine+    :: forall comp. KnownNat comp+    => Maybe FilePath -> (Double -> Double) -> Proxy comp -> [Double] -> IO (Maybe Double)+tdigestSparkingMachine fp dens _ input = do+    mdigest <- fmap validate <$> runT1 machine+    case mdigest of+        Nothing             -> return Nothing+        Just (Left err)     -> fail $ "Validation error: " ++ err+        Just (Right digest) -> do+            printStats fp dens digest+            return $ median digest+  where+    machine :: MachineT IO k (TDigest comp)+    machine+        =  fold mappend mempty+        <~ sparking+        <~ mapping tdigest+        <~ buffered 10000+        <~ source input++printStats :: Maybe FilePath -> (Double -> Double) -> TDigest comp -> IO ()+printStats mfp dens digest = do+    -- Extra: print quantiles+    putStrLn "quantiles"+    for_ ([0.1,0.2..0.9] ++ [0.95,0.99,0.999,0.9999,0.99999]) $ \q ->+        putStrLn $ show q ++ ":" ++ show (quantile q digest)+    putStrLn "cdf"+    for_ ([0, 0.25, 0.5, 1, 2]) $ \x ->+        putStrLn $ show x ++ ": " ++ show (cdf x digest)+    let mi = minimumValue digest+    let ma = maximumValue digest+    let points = flip map [0,0.01..1] $ \x -> mi + (ma - mi) * x+    for_ mfp $ \fp -> do+        putStrLn $ "Writing to " ++ fp+        Chart.toFile Chart.def fp $ do+            Chart.layout_title Chart..= "Histogram"+            color <- Chart.takeColor+            let lineStyle = Chart.def+                  & Chart.line_color .~ color+            Chart.plot $ pure $ tdigestToPlot lineStyle digest+            Chart.plot $ Chart.line "theoretical" [map (\x -> (x, dens x)) points]++tdigestToPlot :: Chart.LineStyle -> TDigest comp -> Chart.Plot Double Double+tdigestToPlot lineStyle digest = Chart.Plot+    { Chart._plot_render     = renderHistogram+    , Chart._plot_legend     = []+    , Chart._plot_all_points = unzip allPoints+    }+  where+    hist = histogram digest+    allPoints = flip map hist $ \(HistBin mi ma w _) ->+        let x = (ma + mi) / 2+            d = ma - mi+            y = w / d / tw+        in (x, y)+    tw = totalWeight digest++    renderHistogram pmap = do+        let fillColor = Chart.blend 0.5 (Chart.opaque Chart.white) (lineStyle ^. Chart.line_color)+        let fillStyle = Chart.def & Chart.fill_color .~ fillColor+        Chart.withLineStyle lineStyle $ Chart.withFillStyle fillStyle $+            for_ hist $ \(HistBin mi ma w _) -> do+                let d = ma - mi+                    y = w / d / tw+                    path = Chart.rectPath $ Chart.Rect+                        (Chart.mapXY pmap (mi,0))+                        (Chart.mapXY pmap (ma,y))+                Chart.alignFillPath path >>= Chart.fillPath+                Chart.alignStrokePath path >>= Chart.strokePath++-------------------------------------------------------------------------------+-- Machine additions+-------------------------------------------------------------------------------++sparking :: Process a a+sparking+    =  asParts+    <~ mapping (\x -> x `using` parList rseq)+    <~ buffered 10++-------------------------------------------------------------------------------+-- Statistics additions+-------------------------------------------------------------------------------++randomStream :: ContGen d => d -> MWC.Seed -> [Double]+randomStream d = go+  where+    continue (xs, seed) = xs ++ go seed+    go seed = continue $ runST $ do+        g <- MWC.restore seed+        -- Generate first 10000 elements+        xs <- replicateM 10000 (genContVar d g)+        seed' <- MWC.save g+        pure (xs, seed')++initSeed :: V.Vector Word32 -> MWC.Seed+initSeed v = runST $ MWC.initialize v >>= MWC.save
+ src/Data/TDigest.hs view
@@ -0,0 +1,89 @@+{-# LANGUAGE ScopedTypeVariables #-}+-- |+-- A new data structure for accurate on-line accumulation of rank-based+-- statistics such as quantiles and trimmed means.+--                .+-- See original paper: "Computing extremely accurate quantiles using t-digest"+-- by Ted Dunning and Otmar Ertl for more details+-- <https://github.com/tdunning/t-digest/blob/master/docs/t-digest-paper/histo.pdf>.+--+-- === Examples+--+-- >>> quantile 0.99 (tdigest [1..1000] :: TDigest 25)+-- Just 990.499...+--+-- >>> quantile 0.99 (tdigest [1..1000] :: TDigest 3)+-- Just 992.7...+--+-- t-Digest is more precise in tails, especially median is imprecise:+--+-- >>> median (tdigest [1..1000] :: TDigest 25)+-- Just 502.5...+--+module Data.TDigest (+    -- * Construction+    TDigest,+    tdigest,++    -- ** Population+    singleton,+    insert,+    insert',++    -- * Compression+    --+    -- |+    --+    -- >>> let digest = foldl' (flip insert') mempty [0..1000] :: TDigest 10+    -- >>> (size digest, size $ compress digest)+    -- (1001,54)+    --+    -- >>> (quantile 0.1 digest, quantile 0.1 $ compress digest)+    -- (Just 99.6...,Just 90.1...)+    --+    -- /Note:/ when values are inserted in more random order,+    -- t-Digest self-compresses on the fly:+    --+    -- >>> let digest = foldl' (flip insert') mempty (fairshuffle [0..1000]) :: TDigest 10+    -- >>> (size digest, size $ compress digest, size $ forceCompress digest)+    -- (77,77,44)+    --+    -- >>> quantile 0.1 digest+    -- Just 96.9...+    --+    compress,+    forceCompress,++    -- * Statistics+    totalWeight,+    minimumValue,+    maximumValue,+    -- ** Histogram+    histogram,+    HistBin (..),+    -- ** Percentile+    median,+    quantile,+    -- ** CDF+    cdf,+    icdf,++    -- * Debug+    valid,+    validate,+    ) where++import Prelude ()+import Prelude.Compat ()++import Data.TDigest.Internal.Tree+import Data.TDigest.Postprocess++-- $setup+-- >>> :set -XDataKinds+-- >>> import Prelude.Compat+-- >>> import Data.List.Compat (foldl')+--+-- >>> let merge [] ys = []; merge xs [] = xs; merge (x:xs) (y:ys) = x : y : merge xs ys+-- >>> let fairshuffle' xs = uncurry merge (splitAt (length xs `div` 2) xs)+-- >>> let fairshuffle xs = iterate fairshuffle' xs !! 5
+ src/Data/TDigest/Internal/Tree.hs view
@@ -0,0 +1,496 @@+{-# LANGUAGE DataKinds             #-}+{-# LANGUAGE KindSignatures        #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+-- | Internals of 'TDigest'.+--+-- Tree implementation is based on /Adams’ Trees Revisited/ by Milan Straka+-- <http://fox.ucw.cz/papers/bbtree/bbtree.pdf>+module Data.TDigest.Internal.Tree where++import Prelude ()+import Prelude.Compat+import Control.DeepSeq        (NFData (..))+import Control.Monad.ST       (ST, runST)+import Data.Binary            (Binary (..))+import Data.Either            (isRight)+import Data.Foldable          (toList)+import Data.List.Compat       (foldl')+import Data.Ord               (comparing)+import Data.Proxy             (Proxy (..))+import Data.Semigroup         (Semigroup (..))+import Data.Semigroup.Reducer (Reducer (..))+import GHC.TypeLits           (KnownNat, Nat, natVal)++import qualified Data.Vector.Algorithms.Heap as VHeap+import qualified Data.Vector.Unboxed         as VU+import qualified Data.Vector.Unboxed.Mutable as MVU++assert :: Bool -> String -> a -> a+assert False msg _ = error msg+assert True  _   x = x++-------------------------------------------------------------------------------+-- TDigest+-------------------------------------------------------------------------------++-- TODO: make newtypes+type Mean = Double+type Weight = Double+type Centroid = (Mean, Weight)+type Size = Int++-- | 'TDigest' is a tree of centroids.+--+-- @compression@ is a @1/δ@. The greater the value of @compression@ the less+-- likely value merging will happen.+data TDigest (compression :: Nat)+    -- | Tree node+    = Node+        {-# UNPACK #-} !Size     -- size of this tree/centroid+        {-# UNPACK #-} !Mean     -- mean of the centroid+        {-# UNPACK #-} !Weight   -- weight of the centrod+        {-# UNPACK #-} !Weight   -- total weight of the tree+        !(TDigest compression)   -- left subtree+        !(TDigest compression)   -- right subtree+    -- | Empty tree+    | Nil+  deriving (Show)++-- [Note: keep min & max in the tree]+--+-- We tried it, but it seems the alloc/update cost is bigger than+-- re-calculating them on need (it's O(log n) - calculation!)++-- [Note: singleton node]+-- We tried to add one, but haven't seen change in performance++-- [Note: inlining balanceR and balanceL]+-- We probably can squueze some performance by making+-- 'balanceL' and 'balanceR' check arguments only once (like @containers@ do)+-- and not use 'node' function.+-- *But*, the benefit vs. code explosion is not yet worth.++instance KnownNat comp => Semigroup (TDigest comp) where+    (<>) = combineDigest++-- | Both 'cons' and 'snoc' are 'insert'+instance KnownNat comp => Reducer Double (TDigest comp) where+    cons = insert+    snoc = flip insert+    unit = singleton++instance  KnownNat comp => Monoid (TDigest comp) where+    mempty  = emptyTDigest+    mappend = combineDigest++-- | 'TDigest' has only strict fields.+instance NFData (TDigest comp) where+    rnf x = x `seq` ()++-- | 'TDigest' isn't compressed after de-serialisation,+-- but it can be still smaller.+instance KnownNat comp => Binary (TDigest comp) where+    put = put . getCentroids+    get = foldl' (flip insertCentroid) emptyTDigest . lc <$> get+      where+        lc :: [Centroid] -> [Centroid]+        lc = id++getCentroids :: TDigest comp -> [Centroid]+getCentroids = ($ []) . go+  where+    go Nil                = id+    go (Node _ x w _ l r) = go l . ((x,w) : ) . go r++-- | Total count of samples.+--+-- >>> totalWeight (tdigest [1..100] :: TDigest 3)+-- 100.0+--+totalWeight :: TDigest comp -> Double+totalWeight Nil                 = 0+totalWeight (Node _ _ _ tw _ _) = tw++size :: TDigest comp -> Int+size Nil                    = 0+size (Node s _ _ _ _ _) = s++-- | Center of left-most centroid. Note: may be different than min element inserted.+--+-- >>> minimumValue (tdigest [1..100] :: TDigest 3)+-- 1.0+--+minimumValue :: TDigest comp -> Mean+minimumValue = go posInf+  where+    go  acc Nil                    = acc+    go _acc (Node _ x _ _ l _) = go x l++-- | Center of right-most centroid. Note: may be different than max element inserted.+--+-- >>> maximumValue (tdigest [1..100] :: TDigest 3)+-- 99.0+--+maximumValue :: TDigest comp -> Mean+maximumValue = go negInf+  where+    go  acc Nil                    = acc+    go _acc (Node _ x _ _ _ r) = go x r++-------------------------------------------------------------------------------+-- Implementation+-------------------------------------------------------------------------------++emptyTDigest :: TDigest comp+emptyTDigest = Nil++combineDigest+    :: KnownNat comp+    => TDigest comp+    -> TDigest comp+    -> TDigest comp+combineDigest a Nil = a+combineDigest Nil b = b+combineDigest a@(Node n _ _ _ _ _) b@(Node m _ _ _ _ _)+    -- TODO: merge first, then shuffle and insert (part of compress)+    | n < m     = compress $ foldl' (flip insertCentroid) b (getCentroids a)+    | otherwise = compress $ foldl' (flip insertCentroid) a (getCentroids b)++insertCentroid+    :: forall comp. KnownNat comp+    => Centroid+    -> TDigest comp+    -> TDigest comp+insertCentroid (x, w) Nil        = singNode x w+insertCentroid (mean, weight) td = go 0 mean weight False td+  where+    -- New weight of the tree+    n :: Weight+    n = totalWeight td + weight++    -- 1/delta+    compression :: Double+    compression = fromInteger $ natVal (Proxy :: Proxy comp)++    go+        :: Weight        -- weight to the left of this tree+        -> Mean          -- mean to insert+        -> Weight        -- weight to insert+        -> Bool          -- should insert everything.+                         -- if we merged somewhere on top, rest is inserted as is+        -> TDigest comp  -- subtree to insert/merge centroid into+        -> TDigest comp+    go _   newX newW _ Nil                 = singNode newX newW+    go cum newX newW e (Node s x w tw l r) = case compare newX x of+        -- Exact match, insert here+        EQ -> Node s x (w + newW) (tw + newW) l r -- node x (w + newW) l r++        -- there is *no* room to insert into this node+        LT | thr <= w -> balanceL x w (go cum newX newW e l) r+        GT | thr <= w -> balanceR x w l (go (cum + totalWeight l + w) newX newW e r)++        -- otherwise go left ... or later right+        LT | e -> balanceL x w (go cum newX newW e l) r+        LT -> case l of+            -- always create a new node+            Nil -> case mrw of+                Nothing     -> node' s nx nw (tw + newW) Nil r+                Just rw     -> balanceL nx nw (go cum newX rw True Nil) r+            Node _ _ _ _ _ _+                | lmax < newX && abs (newX - x) < abs (newX - lmax) {- && newX < x -} -> case mrw of+                    Nothing -> node' s nx nw (tw + nw - w) l r+                    -- in this two last LT cases, we have to recalculate size+                    Just rw -> balanceL nx nw (go cum newX rw True l) r+                | otherwise -> balanceL x w (go cum newX newW e l) r+              where+                lmax = maximumValue l++        -- ... or right+        GT | e -> balanceR x w l (go (cum + totalWeight l + w) newX newW True r)+        GT -> case r of+            Nil -> case mrw of+                Nothing     -> node' s nx nw (tw + newW) l Nil+                Just rw     -> balanceR nx nw l (go (cum + totalWeight l + nw) newX rw True Nil)+            Node _ _ _ _ _ _+                | rmin > newX && abs (newX - x) < abs (newX - rmin) {- && newX > x -} -> case mrw of+                    Nothing -> node' s nx nw (tw + newW) l r+                    -- in this two last GT cases, we have to recalculate size+                    Just rw -> balanceR nx nw l (go (cum + totalWeight l + nw) newX rw True r)+                | otherwise -> balanceR x w l (go (cum + totalWeight l + w) newX newW e r)+              where+                rmin = minimumValue r+      where+        -- quantile approximation of current node+        q   = (w / 2 + cum) / n++        -- threshold, max size of current node/centroid+        thr = {- traceShowId $ traceShow (n, q) $ -} threshold n q compression++        -- We later use nx, nw and mrw:++        -- max size of current node+        dw :: Weight+        mrw :: Maybe Weight+        (dw, mrw) =+            let diff = assert (thr > w) "threshold should be larger than current node weight"+                     $ w + newW - thr+            in if diff < 0 -- i.e. there is room+                then (newW, Nothing)+                else (thr - w, Just $ diff)++        -- the change of current node+        (nx, nw) = {- traceShowId $ traceShow (newX, newW, x, dw, mrw) $ -} combinedCentroid x w x dw++-- | Constructor which calculates size and total weight.+node :: Mean -> Weight -> TDigest comp -> TDigest comp -> TDigest comp+node x w l r = Node+    (1 + size l + size r)+    x w+    (w + totalWeight l + totalWeight r)+    l r++-- | Balance after right insertion.+balanceR :: Mean -> Weight -> TDigest comp -> TDigest comp -> TDigest comp+balanceR x w l r+    | size l + size r <= 1 = node x w l r+    | size r > balOmega * size l = case r of+        Nil -> error "balanceR: impossible happened"+        (Node _ rx rw _ Nil rr) ->+            -- assert (0 < balAlpha * size rr) "balanceR" $+                -- single left rotation+                node rx rw (node x w l Nil) rr+        (Node _ rx rw _ rl rr)+            | size rl < balAlpha * size rr ->+                -- single left rotation+                node rx rw (node x w l rl) rr+        (Node _ rx rw _ (Node _ rlx rlw _ rll rlr) rr) ->+                -- double left rotation+                node rlx rlw (node x w l rll) (node rx rw rlr rr)+    | otherwise            = node x w l r++-- | Balance after left insertion.+balanceL :: Mean -> Weight -> TDigest comp -> TDigest comp -> TDigest comp+balanceL x w l r+    | size l + size r <= 1 = node x w l r+    | size l > balOmega * size r = case l of+        Nil -> error "balanceL: impossible happened"+        (Node _ lx lw _ ll Nil) ->+            -- assert (0 < balAlpha * size ll) "balanceL" $+                -- single right rotation+                node lx lw ll (node x w Nil r)+        (Node _ lx lw _ ll lr)+            | size lr < balAlpha * size ll ->+                -- single right rotation+                node lx lw ll (node x w lr r)+        (Node _ lx lw _ ll (Node _ lrx lrw _ lrl lrr)) ->+                -- double left rotation+                node lrx lrw (node lx lw ll lrl) (node x w lrr r)+    | otherwise = node x w l r++-- | Alias to 'Node'+node' :: Int -> Mean -> Weight -> Weight -> TDigest comp -> TDigest comp -> TDigest comp+node' = Node++-- | Create singular node.+singNode :: Mean -> Weight -> TDigest comp+singNode x w = Node 1 x w w Nil Nil++-- | Add two weighted means together.+combinedCentroid+    :: Mean -> Weight+    -> Mean -> Weight+    -> Centroid+combinedCentroid x w x' w' =+    ( (x * w + x' * w') / w'' -- this is probably not num. stable+    , w''+    )+  where+    w'' = w + w'++-- | Calculate the threshold, i.e. maximum weight of centroid.+threshold+    :: Double  -- ^ total weight+    -> Double  -- ^ quantile+    -> Double  -- ^ compression (1/δ)+    -> Double+threshold n q compression = 4 * n * q * (1 - q) / compression++-------------------------------------------------------------------------------+-- Compression+-------------------------------------------------------------------------------++-- | Compress 'TDigest'.+--+-- Reinsert the centroids in "better" order (in original paper: in random)+-- so they have opportunity to merge.+--+-- Compression will happen only if size is both:+-- bigger than @'relMaxSize' * comp@ and bigger than 'absMaxSize'.+--+compress :: forall comp. KnownNat comp => TDigest comp -> TDigest comp+compress Nil = Nil+compress td+    | size td > relMaxSize * compression && size td > absMaxSize+        = forceCompress td+    | otherwise+        = td+  where+    compression = fromInteger $ natVal (Proxy :: Proxy comp)++-- | Perform compression, even if current size says it's not necessary.+forceCompress :: forall comp. KnownNat comp => TDigest comp -> TDigest comp+forceCompress Nil = Nil+forceCompress td =+    foldl' (flip insertCentroid) emptyTDigest $ fmap fst $ VU.toList centroids+  where+    -- Centroids are shuffled based on space+    centroids :: VU.Vector (Centroid, Double)+    centroids = runST $ do+        v <- toMVector td+        -- sort by cumulative weight+        VHeap.sortBy (comparing snd) v+        f <- VU.unsafeFreeze v+        pure f++toMVector+    :: forall comp s. KnownNat comp+    => TDigest comp                           -- ^ t-Digest+    -> ST s (VU.MVector s (Centroid, Double)) -- ^ return also a "space left in the centroid" value for "shuffling"+toMVector td = do+    v <- MVU.new (size td)+    (i, cum) <- go v (0 :: Int) (0 :: Double) td+    pure $ assert (i == size td && abs (cum - totalWeight td) < 1e-6) "traversal in toMVector:" v+  where+    go _ i cum Nil                   = pure (i, cum)+    go v i cum (Node _ x w _ l r) = do+        (i', cum') <- go v i cum l+        MVU.unsafeWrite v i' ((x, w), space w cum')+        go v (i' + 1) (cum' + w) r++    n = totalWeight td+    compression = fromInteger $ natVal (Proxy :: Proxy comp)++    space w cum = thr - w+      where+        q     = (w / 2 + cum) / n+        thr   = threshold n q compression++-------------------------------------------------------------------------------+-- Params+-------------------------------------------------------------------------------++-- | Relative size parameter. Hard-coded value: 25.+relMaxSize :: Int+relMaxSize = 25++-- | Absolute size parameter. Hard-coded value: 1000.+absMaxSize :: Int+absMaxSize = 1000++-------------------------------------------------------------------------------+-- Tree balance parameters+-------------------------------------------------------------------------------++balOmega :: Int+balOmega = 3++balAlpha :: Int+balAlpha = 2++-- balDelta = 0++-------------------------------------------------------------------------------+-- Debug+-------------------------------------------------------------------------------++-- | @'isRight' . 'validate'@+valid :: TDigest comp -> Bool+valid = isRight . validate++-- | Check various invariants in the 'TDigest' tree.+validate :: TDigest comp -> Either String (TDigest comp)+validate td+    | not (all sizeValid   centroids) = Left "invalid sizes"+    | not (all weightValid centroids) = Left "invalid weights"+    | not (all orderValid  centroids) = Left "invalid ordering"+    | not (all balanced    centroids) = Left "tree is ill-balanced"+    | otherwise = Right td+  where+    centroids = goc td++    goc Nil = []+    goc n@(Node _ _ _ _ l r) = n : goc l ++ goc r++    sizeValid Nil = True+    sizeValid (Node s _ _ _ l r) = s == size l + size r + 1++    weightValid Nil = True+    weightValid (Node _ _ w tw l r) = eq tw $ w + totalWeight l + totalWeight r++    orderValid Nil = True+    orderValid (Node _ _ _ _ Nil                 Nil)                 = True+    orderValid (Node _ x _ _ (Node _ lx _ _ _ _) Nil)                 = lx < x+    orderValid (Node _ x _ _ Nil                 (Node _ rx _ _ _ _)) = x < rx+    orderValid (Node _ x _ _ (Node _ lx _ _ _ _) (Node _ rx _ _ _ _)) = lx < x && x < rx++    balanced Nil = True+    balanced (Node _ _ _ _ l r) =+        size l <= max 1 (balOmega * size r) &&+        size r <= max 1 (balOmega * size l)++-------------------------------------------------------------------------------+-- Double helpers+-------------------------------------------------------------------------------++eq :: Double -> Double -> Bool+eq a b = abs (a-b) < 1e-6++negInf :: Double+negInf = negate posInf++posInf :: Double+posInf = 1/0++-------------------------------------------------------------------------------+-- Higher level helpers+-------------------------------------------------------------------------------++-- | Insert single value into 'TDigest'.+insert+    :: KnownNat comp+    => Double         -- ^ element+    -> TDigest comp+    -> TDigest comp+insert x = compress . insert' x++-- | Insert single value, don't compress 'TDigest' even if needed.+--+-- For sensibly bounded input, it makes sense to let 'TDigest' grow (it might+-- grow linearly in size), and after that compress it once.+insert'+    :: KnownNat comp+    => Double         -- ^ element+    -> TDigest comp+    -> TDigest comp+insert' x = insertCentroid (x, 1)++-- | Make a 'TDigest' of a single data point.+singleton :: KnownNat comp => Double -> TDigest comp+singleton x = insert x emptyTDigest++-- | Strict 'foldl'' over 'Foldable' structure.+tdigest :: (Foldable f, KnownNat comp) => f Double -> TDigest comp+tdigest = forceCompress . foldl' insertChunk emptyTDigest . chunks . toList+  where+    -- compress after each chunk, forceCompress at the very end.+    insertChunk td xs =+        compress (foldl' (flip insert') td xs)++    chunks [] = []+    chunks xs =+        let (a, b) = splitAt 1000 xs -- 1000 is totally arbitrary.+        in a : chunks b++-- $setup+-- >>> :set -XDataKinds
+ src/Data/TDigest/Postprocess.hs view
@@ -0,0 +1,98 @@+-- | 'TDigest' postprocessing functions.+--+-- These are re-exported from "Data.TDigest" module.+--+module Data.TDigest.Postprocess (+    -- * Histogram+    histogram,+    HistBin (..),+    -- * Quantiles+    median,+    quantile,+    -- * CDF+    cdf,+    icdf,+    ) where++import Prelude ()+import Prelude.Compat+import Data.TDigest.Internal.Tree++-------------------------------------------------------------------------------+-- Histogram+-------------------------------------------------------------------------------++-- | Histogram bin+data HistBin = HistBin+    { hbMin       :: !Double  -- ^ lower bound+    , hbMax       :: !Double  -- ^ upper bound+    , hbWeight    :: !Double  -- ^ weight ("area" of the bar)+    , hbCumWeight :: !Double  -- ^ weight from the right+    }+  deriving (Show)++-- | Calculate histogram based on the 'TDigest'.+histogram :: TDigest comp -> [HistBin]+histogram = iter Nothing 0 . getCentroids+  where+    -- zero+    iter :: Maybe (Mean, Weight) -> Weight -> [(Mean, Weight)] -> [HistBin]+    iter _ _ [] = []+    -- one+    iter Nothing t [(x, w)] = [HistBin x x w t]+    -- first+    iter Nothing t (c1@(x1, w1) : rest@((x2, _) : _))+        = HistBin x1 (mid x1 x2) w1 t : iter (Just c1) (t + w1) rest+    -- middle+    iter (Just (x0, _)) t (c1@(x1, w1) : rest@((x2, _) : _))+        = HistBin (mid x0 x1) (mid x1 x2) w1 t: iter (Just c1) (t + w1) rest+    -- last+    iter (Just (x0, _)) t [(x1, w1)]+        = [HistBin (mid x0 x1) x1 w1 t]++    mid a b = (a + b) / 2++-------------------------------------------------------------------------------+-- Quantile+-------------------------------------------------------------------------------++-- | Median, i.e. @'quantile' 0.5@.+median :: TDigest comp -> Maybe Double+median = quantile 0.5++-- | Calculate quantile of a specific value.+quantile :: Double -> TDigest comp -> Maybe Double+quantile q td =+    iter $ histogram td+  where+    q' = q * totalWeight td++    iter []                          = Nothing+    iter [HistBin a b w t]           = Just $ a + (b - a) * (q' - t) / w+    iter (HistBin a b w t : rest)+        | {- t < q' && -} q' < t + w = Just $ a + (b - a) * (q' - t) / w+        | otherwise                  = iter rest++-- | Alias of 'quantile'.+icdf :: Double -> TDigest comp -> Maybe Double+icdf = quantile++-------------------------------------------------------------------------------+-- CDF - cumulative distribution function+-------------------------------------------------------------------------------++-- | Cumulative distribution function.+--+-- /Note:/ if this is the only thing you need, it's more efficient to count+-- this directly.+cdf :: Double -> TDigest comp -> Double+cdf x td = +    iter $ histogram td+  where+    n = totalWeight td++    iter [] = 1+    iter (HistBin a b w t : rest)+        | x < a     = 0+        | x < b     = (t + w * (x - a) / (b - a)) / n+        | otherwise = iter rest
+ tdigest.cabal view
@@ -0,0 +1,112 @@+name:           tdigest+version:        0+synopsis:       On-line accumulation of rank-based statistics+description:    A new data structure for accurate on-line accumulation of rank-based statistics such as quantiles and trimmed means.+                .+                See original paper: "Computing extremely accurate quantiles using t-digest" by Ted Dunning and Otmar Ertl+                for more details <https://github.com/tdunning/t-digest/blob/master/docs/t-digest-paper/histo.pdf>.+category:       Numeric+homepage:       https://github.com/futurice/haskell-tdigest#readme+bug-reports:    https://github.com/futurice/haskell-tdigest/issues+author:         Oleg Grenrus <oleg.grenrus@iki.fi>+maintainer:     Oleg Grenrus <oleg.grenrus@iki.fi>+license:        BSD3+license-file:   LICENSE+tested-with:    GHC==7.8.4, GHC==7.10.3, GHC==8.0.1, GHC==8.0.2+build-type:     Custom+cabal-version:  >= 1.10++extra-source-files:+    README.md++source-repository head+  type: git+  location: https://github.com/futurice/haskell-tdigest++custom-setup+  setup-depends:+    base          >=4.5 && <5,+    Cabal         >=1.14,+    cabal-doctest >=1 && <1.1++library+  hs-source-dirs:+      src+  ghc-options: -Wall+  build-depends:+      base                >=4.7      && <4.10+    , base-compat         >=0.9.1    && <0.10+    , deepseq             >=1.3.0.2  && <1.5+    , binary              >=0.7.1.0  && <0.9+    , reducers            >=3.12.1   && <3.13+    , semigroups          >=0.18.2   && <0.19+    , vector              >=0.11     && <0.13+    , vector-algorithms   >=0.7.0.1  && <0.8+  exposed-modules:+      Data.TDigest+      Data.TDigest.Postprocess+      Data.TDigest.Internal.Tree+  default-language: Haskell2010+  other-extensions:+      DataKinds+      KindSignatures+      MultiParamTypeClasses+      ScopedTypeVariables++test-suite tdigest-tests+  default-language: Haskell2010+  type:             exitcode-stdio-1.0+  main-is:          Tests.hs+  ghc-options:      -Wall -threaded+  hs-source-dirs:   tests++  build-depends:+    base,+    tdigest,+    base-compat,+    deepseq,+    binary,+    semigroups,+    vector,+    vector-algorithms,+    tasty            >=0.11.0.4 && <0.12,+    tasty-quickcheck >=0.8.4    && <0.9++test-suite doctests+  default-language: Haskell2010+  type:             exitcode-stdio-1.0+  main-is:          doctests.hs+  ghc-options:      -Wall -threaded+  hs-source-dirs:   tests++  build-depends:+    base,+    bytes,+    directory      >=1.0,+    doctest        >=0.11.1 && <0.12,+    filepath       >=1.2++benchmark tdigest-simple+  type: exitcode-stdio-1.0+  main-is: Simple.hs+  hs-source-dirs:+      bench+  ghc-options: -Wall -threaded+  build-depends:+      base+    , tdigest+    , base-compat+    , deepseq+    , binary+    , semigroups+    , vector+    , vector-algorithms+    , Chart                >=1.8.1    && <1.9+    , Chart-diagrams       >=1.8.1    && <1.9+    , machines             >=0.6.1    && <0.7+    , parallel             >=3.2.0.6  && <3.3+    , mwc-random           >=0.13.4.0 && <0.14+    , statistics           >=0.13.3.0 && <0.14+    , time                 >=1.4.2    && <1.8+    , optparse-applicative >=0.12.1.0 && <0.14+  default-language: Haskell2010
+ tests/Tests.hs view
@@ -0,0 +1,21 @@+{-# LANGUAGE DataKinds #-}+module Main (main) where++import Data.TDigest+import Test.Tasty+import Test.Tasty.QuickCheck++main :: IO ()+main = defaultMain tests++tests :: TestTree+tests = testGroup "properties"+    [ testProperty "valid" propValid+    ]++propValid :: [Double] -> Property+propValid ds = case validate td of+    Right _  -> property True+    Left err -> counterexample (err ++ " " ++ show td) (valid td)+  where+    td = tdigest ds :: TDigest 2
+ tests/doctests.hs view
@@ -0,0 +1,25 @@+-----------------------------------------------------------------------------+-- |+-- Module      :  Main (doctests)+-- Copyright   :  (C) 2012-14 Edward Kmett+-- License     :  BSD-style (see the file LICENSE)+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>+-- Stability   :  provisional+-- Portability :  portable+--+-- This module provides doctests for a project based on the actual versions+-- of the packages it was built with. It requires a corresponding Setup.lhs+-- to be added to the project+-----------------------------------------------------------------------------+module Main where++import Build_doctests (flags, pkgs, module_sources)+import Data.Foldable (traverse_)+import Test.DocTest++main :: IO ()+main = do+    traverse_ putStrLn args+    doctest args+  where+    args = flags ++ pkgs ++ module_sources