diff --git a/Gauge.hs b/Gauge.hs
new file mode 100644
--- /dev/null
+++ b/Gauge.hs
@@ -0,0 +1,69 @@
+{-# LANGUAGE RecordWildCards #-}
+-- |
+-- Module : Gauge
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Core benchmarking code.
+
+module Gauge
+ (
+ -- * Benchmarkable code
+ Benchmarkable
+ -- * Creating a benchmark suite
+ , Benchmark
+ , env
+ , envWithCleanup
+ , perBatchEnv
+ , perBatchEnvWithCleanup
+ , perRunEnv
+ , perRunEnvWithCleanup
+ , toBenchmarkable
+ , bench
+ , bgroup
+ -- ** Running a benchmark
+ , nf
+ , whnf
+ , nfIO
+ , whnfIO
+ -- * For interactive use
+ , benchmark
+ , benchmarkWith
+ , benchmark'
+ , benchmarkWith'
+ ) where
+
+import Control.Monad (void)
+import Gauge.IO.Printf (note)
+import Gauge.Internal (runAndAnalyseOne)
+import Gauge.Main.Options (defaultConfig)
+import Gauge.Measurement (initializeTime)
+import Gauge.Monad (withConfig)
+import Gauge.Types
+
+-- | Run a benchmark interactively, and analyse its performance.
+benchmark :: Benchmarkable -> IO ()
+benchmark bm = void $ benchmark' bm
+
+-- | Run a benchmark interactively, analyse its performance, and
+-- return the analysis.
+benchmark' :: Benchmarkable -> IO Report
+benchmark' = benchmarkWith' defaultConfig
+
+-- | Run a benchmark interactively, and analyse its performance.
+benchmarkWith :: Config -> Benchmarkable -> IO ()
+benchmarkWith cfg bm = void $ benchmarkWith' cfg bm
+
+-- | Run a benchmark interactively, analyse its performance, and
+-- return the analysis.
+benchmarkWith' :: Config -> Benchmarkable -> IO Report
+benchmarkWith' cfg bm = do
+ initializeTime
+ withConfig cfg $ do
+ _ <- note "benchmarking...\n"
+ Analysed rpt <- runAndAnalyseOne 0 "function" bm
+ return rpt
diff --git a/Gauge/Analysis.hs b/Gauge/Analysis.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Analysis.hs
@@ -0,0 +1,263 @@
+{-# LANGUAGE Trustworthy #-}
+{-# LANGUAGE BangPatterns, DeriveDataTypeable, RecordWildCards #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module : Gauge.Analysis
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Analysis code for benchmarks.
+
+module Gauge.Analysis
+ (
+ Outliers(..)
+ , OutlierEffect(..)
+ , OutlierVariance(..)
+ , SampleAnalysis(..)
+ , analyseSample
+ , scale
+ , analyseMean
+ , countOutliers
+ , classifyOutliers
+ , noteOutliers
+ , outlierVariance
+ , resolveAccessors
+ , validateAccessors
+ , regress
+ ) where
+
+-- Temporary: to support pre-AMP GHC 7.8.4:
+import Data.Monoid
+
+import Control.Arrow (second)
+import Control.Monad (unless, when)
+import Foundation.Monad.Reader
+import Foundation.Monad
+import Gauge.IO.Printf (note, prolix)
+import Gauge.Measurement (secs, threshold)
+import Gauge.Monad (Gauge, getGen, getOverhead)
+import Gauge.Monad.ExceptT
+import Gauge.Types
+import Data.Int (Int64)
+import Data.Maybe (fromJust)
+import Statistics.Function (sort)
+import Statistics.Quantile (weightedAvg, Sorted(..))
+import Statistics.Regression (bootstrapRegress, olsRegress)
+import Statistics.Resampling (Estimator(..),resample)
+import Statistics.Sample (mean)
+import Statistics.Sample.KernelDensity (kde)
+import Statistics.Types (Sample)
+import System.Random.MWC (GenIO)
+import qualified Data.List as List
+import qualified Data.Map as Map
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import qualified Statistics.Resampling.Bootstrap as B
+import qualified Statistics.Types as B
+import Prelude
+
+-- | Classify outliers in a data set, using the boxplot technique.
+classifyOutliers :: Sample -> Outliers
+classifyOutliers sa = U.foldl' ((. outlier) . mappend) mempty ssa
+ where outlier e = Outliers {
+ samplesSeen = 1
+ , lowSevere = if e <= loS && e < hiM then 1 else 0
+ , lowMild = if e > loS && e <= loM then 1 else 0
+ , highMild = if e >= hiM && e < hiS then 1 else 0
+ , highSevere = if e >= hiS && e > loM then 1 else 0
+ }
+ !loS = q1 - (iqr * 3)
+ !loM = q1 - (iqr * 1.5)
+ !hiM = q3 + (iqr * 1.5)
+ !hiS = q3 + (iqr * 3)
+ q1 = weightedAvg 1 4 (Sorted ssa)
+ q3 = weightedAvg 3 4 (Sorted ssa)
+ ssa = sort sa
+ iqr = q3 - q1
+
+-- | Compute the extent to which outliers in the sample data affect
+-- the sample mean and standard deviation.
+outlierVariance
+ :: B.Estimate B.ConfInt Double -- ^ Bootstrap estimate of sample mean.
+ -> B.Estimate B.ConfInt Double -- ^ Bootstrap estimate of sample
+ -- standard deviation.
+ -> Double -- ^ Number of original iterations.
+ -> OutlierVariance
+outlierVariance µ σ a = OutlierVariance effect desc varOutMin
+ where
+ ( effect, desc ) | varOutMin < 0.01 = (Unaffected, "no")
+ | varOutMin < 0.1 = (Slight, "slight")
+ | varOutMin < 0.5 = (Moderate, "moderate")
+ | otherwise = (Severe, "severe")
+ varOutMin = (minBy varOut 1 (minBy cMax 0 µgMin)) / σb2
+ varOut c = (ac / a) * (σb2 - ac * σg2) where ac = a - c
+ σb = B.estPoint σ
+ µa = B.estPoint µ / a
+ µgMin = µa / 2
+ σg = min (µgMin / 4) (σb / sqrt a)
+ σg2 = σg * σg
+ σb2 = σb * σb
+ minBy f q r = min (f q) (f r)
+ cMax x = fromIntegral (floor (-2 * k0 / (k1 + sqrt det)) :: Int)
+ where
+ k1 = σb2 - a * σg2 + ad
+ k0 = -a * ad
+ ad = a * d
+ d = k * k where k = µa - x
+ det = k1 * k1 - 4 * σg2 * k0
+
+-- | Count the total number of outliers in a sample.
+countOutliers :: Outliers -> Int64
+countOutliers (Outliers _ a b c d) = a + b + c + d
+{-# INLINE countOutliers #-}
+
+-- | Display the mean of a 'Sample', and characterise the outliers
+-- present in the sample.
+analyseMean :: Sample
+ -> Int -- ^ Number of iterations used to
+ -- compute the sample.
+ -> Gauge Double
+analyseMean a iters = do
+ let µ = mean a
+ _ <- note "mean is %s (%d iterations)\n" (secs µ) iters
+ noteOutliers . classifyOutliers $ a
+ return µ
+
+-- | Multiply the 'Estimate's in an analysis by the given value, using
+-- 'B.scale'.
+scale :: Double -- ^ Value to multiply by.
+ -> SampleAnalysis -> SampleAnalysis
+scale f s@SampleAnalysis{..} = s {
+ anMean = B.scale f anMean
+ , anStdDev = B.scale f anStdDev
+ }
+
+-- | Perform an analysis of a measurement.
+analyseSample :: Int -- ^ Experiment number.
+ -> String -- ^ Experiment name.
+ -> V.Vector Measured -- ^ Sample data.
+ -> ExceptT String Gauge Report
+analyseSample i name meas = do
+ Config{..} <- ask
+ overhead <- lift getOverhead
+ let ests = [Mean,StdDev]
+ -- The use of filter here throws away very-low-quality
+ -- measurements when bootstrapping the mean and standard
+ -- deviations. Without this, the numbers look nonsensical when
+ -- very brief actions are measured.
+ stime = measure (measTime . rescale) .
+ G.filter ((>= threshold) . measTime) . G.map fixTime .
+ G.tail $ meas
+ fixTime m = m { measTime = measTime m - overhead / 2 }
+ n = G.length meas
+ s = G.length stime
+ _ <- lift $ prolix "bootstrapping with %d of %d samples (%d%%)\n"
+ s n ((s * 100) `quot` n)
+ gen <- lift getGen
+ rs <- mapM (\(ps,r) -> regress gen ps r meas) $
+ ((["iters"],"time"):regressions)
+ resamps <- liftIO $ resample gen ests resamples stime
+ let [estMean,estStdDev] = B.bootstrapBCA confInterval stime resamps
+ ov = outlierVariance estMean estStdDev (fromIntegral n)
+ an = SampleAnalysis {
+ anRegress = rs
+ , anOverhead = overhead
+ , anMean = estMean
+ , anStdDev = estStdDev
+ , anOutlierVar = ov
+ }
+ return Report {
+ reportNumber = i
+ , reportName = name
+ , reportKeys = measureKeys
+ , reportMeasured = meas
+ , reportAnalysis = an
+ , reportOutliers = classifyOutliers stime
+ , reportKDEs = [uncurry (KDE "time") (kde 128 stime)]
+ }
+
+
+-- | Regress the given predictors against the responder.
+--
+-- Errors may be returned under various circumstances, such as invalid
+-- names or lack of needed data.
+--
+-- See 'olsRegress' for details of the regression performed.
+regress :: GenIO
+ -> [String] -- ^ Predictor names.
+ -> String -- ^ Responder name.
+ -> V.Vector Measured
+ -> ExceptT String Gauge Regression
+regress gen predNames respName meas = do
+ when (G.null meas) $
+ mFail "no measurements"
+ accs <- ExceptT . return $ validateAccessors predNames respName
+ let unmeasured = [n | (n, Nothing) <- map (second ($ G.head meas)) accs]
+ unless (null unmeasured) $
+ mFail $ "no data available for " ++ renderNames unmeasured
+ let (r:ps) = map ((`measure` meas) . (fromJust .) . snd) accs
+ Config{..} <- ask
+ (coeffs,r2) <- liftIO $
+ bootstrapRegress gen resamples confInterval olsRegress ps r
+ return Regression {
+ regResponder = respName
+ , regCoeffs = Map.fromList (zip (predNames ++ ["y"]) (G.toList coeffs))
+ , regRSquare = r2
+ }
+
+singleton :: [a] -> Bool
+singleton [_] = True
+singleton _ = False
+
+-- | Given a list of accessor names (see 'measureKeys'), return either
+-- a mapping from accessor name to function or an error message if
+-- any names are wrong.
+resolveAccessors :: [String]
+ -> Either String [(String, Measured -> Maybe Double)]
+resolveAccessors names =
+ case unresolved of
+ [] -> Right [(n, a) | (n, Just (a,_)) <- accessors]
+ _ -> Left $ "unknown metric " ++ renderNames unresolved
+ where
+ unresolved = [n | (n, Nothing) <- accessors]
+ accessors = flip map names $ \n -> (n, Map.lookup n measureAccessors)
+
+-- | Given predictor and responder names, do some basic validation,
+-- then hand back the relevant accessors.
+validateAccessors :: [String] -- ^ Predictor names.
+ -> String -- ^ Responder name.
+ -> Either String [(String, Measured -> Maybe Double)]
+validateAccessors predNames respName = do
+ when (null predNames) $
+ Left "no predictors specified"
+ let names = respName:predNames
+ dups = map head . filter (not . singleton) .
+ List.group . List.sort $ names
+ unless (null dups) $
+ Left $ "duplicated metric " ++ renderNames dups
+ resolveAccessors names
+
+renderNames :: [String] -> String
+renderNames = List.intercalate ", " . map show
+
+-- | Display a report of the 'Outliers' present in a 'Sample'.
+noteOutliers :: Outliers -> Gauge ()
+noteOutliers o = do
+ let frac n = (100::Double) * fromIntegral n / fromIntegral (samplesSeen o)
+ check :: Int64 -> Double -> String -> Gauge ()
+ check k t d = when (frac k > t) $
+ note " %d (%.1g%%) %s\n" k (frac k) d
+ outCount = countOutliers o
+ when (outCount > 0) $ do
+ _ <- note "found %d outliers among %d samples (%.1g%%)\n"
+ outCount (samplesSeen o) (frac outCount)
+ check (lowSevere o) 0 "low severe"
+ check (lowMild o) 1 "low mild"
+ check (highMild o) 1 "high mild"
+ check (highSevere o) 0 "high severe"
diff --git a/Gauge/IO/Printf.hs b/Gauge/IO/Printf.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/IO/Printf.hs
@@ -0,0 +1,90 @@
+-- |
+-- Module : Gauge.IO.Printf
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Input and output actions.
+
+{-# LANGUAGE FlexibleInstances, Rank2Types, TypeSynonymInstances #-}
+module Gauge.IO.Printf
+ (
+ CritHPrintfType
+ , note
+ , printError
+ , prolix
+ ) where
+
+import Control.Monad (when)
+import Foundation.Monad.Reader (ask)
+import Foundation.Monad (liftIO)
+import Gauge.Monad (Gauge)
+import Gauge.Types (Config(verbosity), Verbosity(..))
+import System.IO (Handle, hFlush, stderr, stdout)
+import Text.Printf (PrintfArg)
+import qualified Text.Printf (HPrintfType, hPrintf)
+
+-- First item is the action to print now, given all the arguments
+-- gathered together so far. The second item is the function that
+-- will take a further argument and give back a new PrintfCont.
+data PrintfCont = PrintfCont (IO ()) (forall a . PrintfArg a => a -> PrintfCont)
+
+-- | An internal class that acts like Printf/HPrintf.
+--
+-- The implementation is visible to the rest of the program, but the
+-- details of the class are not.
+class CritHPrintfType a where
+ chPrintfImpl :: (Config -> Bool) -> PrintfCont -> a
+
+
+instance CritHPrintfType (Gauge a) where
+ chPrintfImpl check (PrintfCont final _)
+ = do x <- ask
+ when (check x) (liftIO (final >> hFlush stderr >> hFlush stdout))
+ return undefined
+
+instance CritHPrintfType (IO a) where
+ chPrintfImpl _ (PrintfCont final _)
+ = final >> hFlush stderr >> hFlush stdout >> return undefined
+
+instance (CritHPrintfType r, PrintfArg a) => CritHPrintfType (a -> r) where
+ chPrintfImpl check (PrintfCont _ anotherArg) x
+ = chPrintfImpl check (anotherArg x)
+
+chPrintf :: CritHPrintfType r => (Config -> Bool) -> Handle -> String -> r
+chPrintf shouldPrint h s
+ = chPrintfImpl shouldPrint (make (Text.Printf.hPrintf h s)
+ (Text.Printf.hPrintf h s))
+ where
+ make :: IO () -> (forall a r. (PrintfArg a, Text.Printf.HPrintfType r) =>
+ a -> r) -> PrintfCont
+ make curCall curCall' = PrintfCont curCall (\x -> make (curCall' x)
+ (curCall' x))
+
+{- A demonstration of how to write printf in this style, in case it is
+ever needed
+ in fututre:
+
+cPrintf :: CritHPrintfType r => (Config -> Bool) -> String -> r
+cPrintf shouldPrint s
+ = chPrintfImpl shouldPrint (make (Text.Printf.printf s)
+ (Text.Printf.printf s))
+ where
+ make :: IO () -> (forall a r. (PrintfArg a, Text.Printf.PrintfType r) => a -> r) -> PrintfCont
+ make curCall curCall' = PrintfCont curCall (\x -> make (curCall' x) (curCall' x))
+-}
+
+-- | Print a \"normal\" note.
+note :: (CritHPrintfType r) => String -> r
+note = chPrintf ((> Quiet) . verbosity) stdout
+
+-- | Print verbose output.
+prolix :: (CritHPrintfType r) => String -> r
+prolix = chPrintf ((== Verbose) . verbosity) stdout
+
+-- | Print an error message.
+printError :: (CritHPrintfType r) => String -> r
+printError = chPrintf (const True) stderr
diff --git a/Gauge/Internal.hs b/Gauge/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Internal.hs
@@ -0,0 +1,209 @@
+{-# LANGUAGE BangPatterns, RecordWildCards #-}
+-- |
+-- Module : Gauge.Internal
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Core benchmarking code.
+
+module Gauge.Internal
+ (
+ runAndAnalyse
+ , runAndAnalyseOne
+ , runFixedIters
+ ) where
+
+import Control.DeepSeq (rnf)
+import Control.Exception (evaluate)
+import Control.Monad (foldM, forM_, void, when)
+import Foundation.Monad
+import Foundation.Monad.Reader (ask)
+import Data.Int (Int64)
+import Gauge.Analysis (analyseSample, noteOutliers)
+import Gauge.IO.Printf (note, printError, prolix)
+import Gauge.Measurement (runBenchmark, runBenchmarkable_, secs)
+import Gauge.Monad (Gauge)
+import Gauge.Monad.ExceptT
+import Gauge.Types hiding (measure)
+import qualified Data.Map as Map
+import qualified Data.Vector as V
+import Statistics.Types (Estimate(..),ConfInt(..),confidenceInterval,cl95,confidenceLevel)
+import Text.Printf (printf)
+
+-- | Run a single benchmark.
+runOne :: Int -> String -> Benchmarkable -> Gauge DataRecord
+runOne i desc bm = do
+ Config{..} <- ask
+ (meas,timeTaken) <- liftIO $ runBenchmark bm timeLimit
+ when (timeTaken > timeLimit * 1.25) .
+ void $ prolix "measurement took %s\n" (secs timeTaken)
+ return (Measurement i desc meas)
+
+-- | Analyse a single benchmark.
+analyseOne :: Int -> String -> V.Vector Measured -> Gauge DataRecord
+analyseOne i desc meas = do
+ Config{..} <- ask
+ _ <- prolix "analysing with %d resamples\n" resamples
+ erp <- runExceptT $ analyseSample i desc meas
+ case erp of
+ Left err -> printError "*** Error: %s\n" err
+ Right rpt@Report{..} -> do
+ let SampleAnalysis{..} = reportAnalysis
+ OutlierVariance{..} = anOutlierVar
+ wibble = case ovEffect of
+ Unaffected -> "unaffected" :: String
+ Slight -> "slightly inflated"
+ Moderate -> "moderately inflated"
+ Severe -> "severely inflated"
+ (builtin, others) = splitAt 1 anRegress
+ let r2 n = printf "%.3f R\178" n
+ forM_ builtin $ \Regression{..} ->
+ case Map.lookup "iters" regCoeffs of
+ Nothing -> return ()
+ Just t -> bs secs "time" t >> bs r2 "" regRSquare
+ bs secs "mean" anMean
+ bs secs "std dev" anStdDev
+ forM_ others $ \Regression{..} -> do
+ _ <- bs r2 (regResponder ++ ":") regRSquare
+ forM_ (Map.toList regCoeffs) $ \(prd,val) ->
+ bs (printf "%.3g") (" " ++ prd) val
+ --writeCsv
+ -- (desc,
+ -- estPoint anMean, fst $ confidenceInterval anMean, snd $ confidenceInterval anMean,
+ -- estPoint anStdDev, fst $ confidenceInterval anStdDev, snd $ confidenceInterval anStdDev
+ -- )
+ when (verbosity == Verbose || (ovEffect > Slight && verbosity > Quiet)) $ do
+ when (verbosity == Verbose) $ noteOutliers reportOutliers
+ _ <- note "variance introduced by outliers: %d%% (%s)\n"
+ (round (ovFraction * 100) :: Int) wibble
+ return ()
+ _ <- note "\n"
+ return (Analysed rpt)
+ where bs :: (Double -> String) -> String -> Estimate ConfInt Double -> Gauge ()
+ bs f metric e@Estimate{..} =
+ note "%-20s %-10s (%s .. %s%s)\n" metric
+ (f estPoint) (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e)
+ (let cl = confIntCL estError
+ str | cl == cl95 = ""
+ | otherwise = printf ", ci %.3f" (confidenceLevel cl)
+ in str
+ )
+
+
+-- | Run a single benchmark and analyse its performance.
+runAndAnalyseOne :: Int -> String -> Benchmarkable -> Gauge DataRecord
+runAndAnalyseOne i desc bm = do
+ Measurement _ _ meas <- runOne i desc bm
+ analyseOne i desc meas
+
+-- | Run, and analyse, one or more benchmarks.
+runAndAnalyse :: (String -> Bool) -- ^ A predicate that chooses
+ -- whether to run a benchmark by its
+ -- name.
+ -> Benchmark
+ -> Gauge ()
+runAndAnalyse select bs = do
+ -- The type we write to the file is ReportFileContents, a triple.
+ -- But here we ASSUME that the tuple will become a JSON array.
+ -- This assumption lets us stream the reports to the file incrementally:
+ --liftIO $ hPutStr handle $ "[ \"" ++ headerRoot ++ "\", " ++
+ -- "\"" ++ critVersion ++ "\", [ "
+
+ for select bs $ \idx desc bm -> do
+ _ <- note "benchmarking %s\n" desc
+ Analysed _ <- runAndAnalyseOne idx desc bm
+ return ()
+ --unless (idx == 0) $
+ -- liftIO $ hPutStr handle ", "
+ {-
+ liftIO $ L.hPut handle (Aeson.encode (rpt::Report))
+ -}
+
+ --liftIO $ hPutStr handle " ] ]\n"
+ --liftIO $ hClose handle
+
+ return ()
+{-
+ rpts <- liftIO $ do
+ res <- readJSONReports jsonFile
+ case res of
+ Left err -> error $ "error reading file "++jsonFile++":\n "++show err
+ Right (_,_,rs) ->
+ case mbJsonFile of
+ Just _ -> return rs
+ _ -> removeFile jsonFile >> return rs
+
+ rawReport rpts
+ report rpts
+ json rpts
+ junit rpts
+ -}
+
+
+-- | Run a benchmark without analysing its performance.
+runFixedIters :: Int64 -- ^ Number of loop iterations to run.
+ -> (String -> Bool) -- ^ A predicate that chooses
+ -- whether to run a benchmark by its
+ -- name.
+ -> Benchmark
+ -> Gauge ()
+runFixedIters iters select bs =
+ for select bs $ \_idx desc bm -> do
+ _ <- note "benchmarking %s\n" desc
+ liftIO $ runBenchmarkable_ bm iters
+
+-- | Iterate over benchmarks.
+for :: (String -> Bool)
+ -> Benchmark
+ -> (Int -> String -> Benchmarkable -> Gauge ())
+ -> Gauge ()
+for select bs0 handle = go (0::Int) ("", bs0) >> return ()
+ where
+ go !idx (pfx, Environment mkenv cleanenv mkbench)
+ | shouldRun pfx mkbench = do
+ e <- liftIO $ do
+ ee <- mkenv
+ evaluate (rnf ee)
+ return ee
+ go idx (pfx, mkbench e) `finally` liftIO (cleanenv e)
+ | otherwise = return idx
+ go idx (pfx, Benchmark desc b)
+ | select desc' = do handle idx desc' b; return $! idx + 1
+ | otherwise = return idx
+ where desc' = addPrefix pfx desc
+ go idx (pfx, BenchGroup desc bs) =
+ foldM go idx [(addPrefix pfx desc, b) | b <- bs]
+
+ shouldRun pfx mkbench =
+ any (select . addPrefix pfx) . benchNames . mkbench $
+ error "Gauge.env could not determine the list of your benchmarks since they force the environment (see the documentation for details)"
+
+{-
+-- | Write summary JUnit file (if applicable)
+junit :: [Report] -> Gauge ()
+junit rs
+ = do junitOpt <- asks junitFile
+ case junitOpt of
+ Just fn -> liftIO $ writeFile fn msg
+ Nothing -> return ()
+ where
+ msg = "\n" ++
+ printf "\n"
+ (length rs) ++
+ concatMap single rs ++
+ "\n"
+ single Report{..} = printf " \n"
+ (attrEsc reportName) (estPoint $ anMean $ reportAnalysis)
+ attrEsc = concatMap esc
+ where
+ esc '\'' = "'"
+ esc '"' = """
+ esc '<' = "<"
+ esc '>' = ">"
+ esc '&' = "&"
+ esc c = [c]
+-}
diff --git a/Gauge/Main.hs b/Gauge/Main.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Main.hs
@@ -0,0 +1,259 @@
+{-# LANGUAGE Trustworthy #-}
+
+-- |
+-- Module : Gauge.Main
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Wrappers for compiling and running benchmarks quickly and easily.
+-- See 'defaultMain' below for an example.
+
+module Gauge.Main
+ (
+ -- * How to write benchmarks
+ -- $bench
+
+ -- ** Benchmarking IO actions
+ -- $io
+
+ -- ** Benchmarking pure code
+ -- $pure
+
+ -- ** Fully evaluating a result
+ -- $rnf
+
+ -- * Types
+ Benchmarkable
+ , Benchmark
+ -- * Creating a benchmark suite
+ , env
+ , envWithCleanup
+ , perBatchEnv
+ , perBatchEnvWithCleanup
+ , perRunEnv
+ , perRunEnvWithCleanup
+ , toBenchmarkable
+ , bench
+ , bgroup
+ -- ** Running a benchmark
+ , nf
+ , whnf
+ , nfIO
+ , whnfIO
+ -- * Turning a suite of benchmarks into a program
+ , defaultMain
+ , defaultMainWith
+ , defaultConfig
+ -- * Other useful code
+ , makeMatcher
+ , runMode
+ ) where
+
+import Control.Monad (unless)
+import Foundation.Monad
+import Gauge.IO.Printf (printError)
+import Gauge.Internal (runAndAnalyse, runFixedIters)
+import Gauge.Main.Options (MatchType(..), Mode(..), defaultConfig, describe,
+ versionInfo)
+import Gauge.Measurement (initializeTime)
+import Gauge.Monad (withConfig)
+import Gauge.Types
+import Data.Char (toLower)
+import Data.List (isInfixOf, isPrefixOf, sort)
+import Options.Applicative (execParser)
+import System.Environment (getProgName)
+import System.Exit (ExitCode(..), exitWith)
+-- import System.FilePath.Glob
+import System.IO.CodePage (withCP65001)
+
+-- | An entry point that can be used as a @main@ function.
+--
+-- > import Gauge.Main
+-- >
+-- > fib :: Int -> Int
+-- > fib 0 = 0
+-- > fib 1 = 1
+-- > fib n = fib (n-1) + fib (n-2)
+-- >
+-- > main = defaultMain [
+-- > bgroup "fib" [ bench "10" $ whnf fib 10
+-- > , bench "35" $ whnf fib 35
+-- > , bench "37" $ whnf fib 37
+-- > ]
+-- > ]
+defaultMain :: [Benchmark] -> IO ()
+defaultMain = defaultMainWith defaultConfig
+
+-- | Create a function that can tell if a name given on the command
+-- line matches a benchmark.
+makeMatcher :: MatchType
+ -> [String]
+ -- ^ Command line arguments.
+ -> Either String (String -> Bool)
+makeMatcher matchKind args =
+ case matchKind of
+ Prefix -> Right $ \b -> null args || any (`isPrefixOf` b) args
+ Pattern -> Right $ \b -> null args || any (`isInfixOf` b) args
+ IPattern -> Right $ \b -> null args || any (`isInfixOf` map toLower b) (map (map toLower) args)
+
+selectBenches :: MatchType -> [String] -> Benchmark -> IO (String -> Bool)
+selectBenches matchType benches bsgroup = do
+ toRun <- either parseError return . makeMatcher matchType $ benches
+ unless (null benches || any toRun (benchNames bsgroup)) $
+ parseError "none of the specified names matches a benchmark"
+ return toRun
+
+-- | An entry point that can be used as a @main@ function, with
+-- configurable defaults.
+--
+-- Example:
+--
+-- > import Gauge.Main.Options
+-- > import Gauge.Main
+-- >
+-- > myConfig = defaultConfig {
+-- > -- Do not GC between runs.
+-- > forceGC = False
+-- > }
+-- >
+-- > main = defaultMainWith myConfig [
+-- > bench "fib 30" $ whnf fib 30
+-- > ]
+--
+-- If you save the above example as @\"Fib.hs\"@, you should be able
+-- to compile it as follows:
+--
+-- > ghc -O --make Fib
+--
+-- Run @\"Fib --help\"@ on the command line to get a list of command
+-- line options.
+defaultMainWith :: Config
+ -> [Benchmark]
+ -> IO ()
+defaultMainWith defCfg bs = withCP65001 $ do
+ wat <- execParser (describe defCfg)
+ runMode wat bs
+
+-- | Run a set of 'Benchmark's with the given 'Mode'.
+--
+-- This can be useful if you have a 'Mode' from some other source (e.g. from a
+-- one in your benchmark driver's command-line parser).
+runMode :: Mode -> [Benchmark] -> IO ()
+runMode wat bs =
+ case wat of
+ List -> mapM_ putStrLn . sort . concatMap benchNames $ bs
+ Version -> putStrLn versionInfo
+ RunIters cfg iters matchType benches -> do
+ shouldRun <- selectBenches matchType benches bsgroup
+ withConfig cfg $
+ runFixedIters iters shouldRun bsgroup
+ Run cfg matchType benches -> do
+ shouldRun <- selectBenches matchType benches bsgroup
+ withConfig cfg $ do
+ --writeCsv ("Name","Mean","MeanLB","MeanUB","Stddev","StddevLB",
+ -- "StddevUB")
+ liftIO initializeTime
+ runAndAnalyse shouldRun bsgroup
+ where bsgroup = BenchGroup "" bs
+
+-- | Display an error message from a command line parsing failure, and
+-- exit.
+parseError :: String -> IO a
+parseError msg = do
+ _ <- printError "Error: %s\n" msg
+ _ <- printError "Run \"%s --help\" for usage information\n" =<< getProgName
+ exitWith (ExitFailure 64)
+
+-- $bench
+--
+-- The 'Benchmarkable' type is a container for code that can be
+-- benchmarked. The value inside must run a benchmark the given
+-- number of times. We are most interested in benchmarking two
+-- things:
+--
+-- * 'IO' actions. Any 'IO' action can be benchmarked directly.
+--
+-- * Pure functions. GHC optimises aggressively when compiling with
+-- @-O@, so it is easy to write innocent-looking benchmark code that
+-- doesn't measure the performance of a pure function at all. We
+-- work around this by benchmarking both a function and its final
+-- argument together.
+
+-- $io
+--
+-- Any 'IO' action can be benchmarked easily if its type resembles
+-- this:
+--
+-- @
+-- 'IO' a
+-- @
+
+-- $pure
+--
+-- Because GHC optimises aggressively when compiling with @-O@, it is
+-- potentially easy to write innocent-looking benchmark code that will
+-- only be evaluated once, for which all but the first iteration of
+-- the timing loop will be timing the cost of doing nothing.
+--
+-- To work around this, we provide two functions for benchmarking pure
+-- code.
+--
+-- The first will cause results to be fully evaluated to normal form
+-- (NF):
+--
+-- @
+-- 'nf' :: 'NFData' b => (a -> b) -> a -> 'Benchmarkable'
+-- @
+--
+-- The second will cause results to be evaluated to weak head normal
+-- form (the Haskell default):
+--
+-- @
+-- 'whnf' :: (a -> b) -> a -> 'Benchmarkable'
+-- @
+--
+-- As both of these types suggest, when you want to benchmark a
+-- function, you must supply two values:
+--
+-- * The first element is the function, saturated with all but its
+-- last argument.
+--
+-- * The second element is the last argument to the function.
+--
+-- Here is an example that makes the use of these functions clearer.
+-- Suppose we want to benchmark the following function:
+--
+-- @
+-- firstN :: Int -> [Int]
+-- firstN k = take k [(0::Int)..]
+-- @
+--
+-- So in the easy case, we construct a benchmark as follows:
+--
+-- @
+-- 'nf' firstN 1000
+-- @
+
+-- $rnf
+--
+-- The 'whnf' harness for evaluating a pure function only evaluates
+-- the result to weak head normal form (WHNF). If you need the result
+-- evaluated all the way to normal form, use the 'nf' function to
+-- force its complete evaluation.
+--
+-- Using the @firstN@ example from earlier, to naive eyes it might
+-- /appear/ that the following code ought to benchmark the production
+-- of the first 1000 list elements:
+--
+-- @
+-- 'whnf' firstN 1000
+-- @
+--
+-- Since we are using 'whnf', in this case the result will only be
+-- forced until it reaches WHNF, so what this would /actually/
+-- benchmark is merely how long it takes to produce the first list
+-- element!
diff --git a/Gauge/Main/Options.hs b/Gauge/Main/Options.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Main/Options.hs
@@ -0,0 +1,200 @@
+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric, RecordWildCards #-}
+
+-- |
+-- Module : Gauge.Main.Options
+-- Copyright : (c) 2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Benchmarking command-line configuration.
+
+module Gauge.Main.Options
+ (
+ Mode(..)
+ , MatchType(..)
+ , defaultConfig
+ , describe
+ , versionInfo
+ ) where
+
+-- Temporary: to support pre-AMP GHC 7.8.4:
+import Data.Monoid
+
+import Control.Monad (when)
+import Gauge.Analysis (validateAccessors)
+import Gauge.Types (Config(..), Verbosity(..), measureAccessors,
+ measureKeys)
+import Data.Char (isSpace, toLower)
+import Data.Data (Data, Typeable)
+import Data.Int (Int64)
+import Data.List (isPrefixOf)
+import Data.Version (showVersion)
+import GHC.Generics (Generic)
+import Options.Applicative
+import Options.Applicative.Help (Chunk(..), tabulate)
+import Options.Applicative.Help.Pretty ((.$.))
+import Options.Applicative.Types
+import Paths_gauge (version)
+import Statistics.Types (mkCL,cl95)
+import Text.PrettyPrint.ANSI.Leijen (Doc, text)
+import qualified Data.Map as M
+import Prelude
+
+-- | How to match a benchmark name.
+data MatchType = Prefix
+ -- ^ Match by prefix. For example, a prefix of
+ -- @\"foo\"@ will match @\"foobar\"@.
+ | Pattern
+ -- ^ Match by searching given substring in benchmark
+ -- paths.
+ | IPattern
+ -- ^ Same as 'Pattern', but case insensitive.
+ deriving (Eq, Ord, Bounded, Enum, Read, Show, Typeable, Data,
+ Generic)
+
+-- | Execution mode for a benchmark program.
+data Mode = List
+ -- ^ List all benchmarks.
+ | Version
+ -- ^ Print the version.
+ | RunIters Config Int64 MatchType [String]
+ -- ^ Run the given benchmarks, without collecting or
+ -- analysing performance numbers.
+ | Run Config MatchType [String]
+ -- ^ Run and analyse the given benchmarks.
+ deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+-- | Default benchmarking configuration.
+defaultConfig :: Config
+defaultConfig = Config {
+ confInterval = cl95
+ , forceGC = True
+ , timeLimit = 5
+ , resamples = 1000
+ , regressions = []
+ , rawDataFile = Nothing
+ , reportFile = Nothing
+ , csvFile = Nothing
+ , jsonFile = Nothing
+ , junitFile = Nothing
+ , verbosity = Normal
+ , template = "default"
+ }
+
+-- | Parse a command line.
+parseWith :: Config
+ -- ^ Default configuration to use if options are not
+ -- explicitly specified.
+ -> Parser Mode
+parseWith cfg =
+ (matchNames (Run <$> config cfg)) <|>
+ runIters <|>
+ (List <$ switch (long "list" <> short 'l' <> help "List benchmarks")) <|>
+ (Version <$ switch (long "version" <> help "Show version info"))
+ where
+ runIters = matchNames $
+ RunIters <$> config cfg <*> option auto
+ (long "iters" <> short 'n' <> metavar "ITERS" <>
+ help "Run benchmarks, don't analyse")
+ matchNames wat = wat
+ <*> option match
+ (long "match" <> short 'm' <> metavar "MATCH" <> value Prefix <>
+ help "How to match benchmark names (\"prefix\", \"glob\", \"pattern\", or \"ipattern\")")
+ <*> many (argument str (metavar "NAME..."))
+
+-- | Parse a configuration.
+config :: Config -> Parser Config
+config Config{..} = Config
+ <$> option (mkCL <$> range 0.001 0.999)
+ (long "ci" <> short 'I' <> metavar "CI" <> value confInterval <>
+ help "Confidence interval")
+ <*> (not <$> switch (long "no-gc" <> short 'G' <>
+ help "Do not collect garbage between iterations"))
+ <*> option (range 0.1 86400)
+ (long "time-limit" <> short 'L' <> metavar "SECS" <> value timeLimit <>
+ help "Time limit to run a benchmark")
+ <*> option (range 1 1000000)
+ (long "resamples" <> metavar "COUNT" <> value resamples <>
+ help "Number of bootstrap resamples to perform")
+ <*> many (option regressParams
+ (long "regress" <> metavar "RESP:PRED.." <>
+ help "Regressions to perform"))
+ <*> outputOption rawDataFile (long "raw" <>
+ help "File to write raw data to")
+ <*> outputOption reportFile (long "output" <> short 'o' <>
+ help "File to write report to")
+ <*> outputOption csvFile (long "csv" <>
+ help "File to write CSV summary to")
+ <*> outputOption jsonFile (long "json" <>
+ help "File to write JSON summary to")
+ <*> outputOption junitFile (long "junit" <>
+ help "File to write JUnit summary to")
+ <*> (toEnum <$> option (range 0 2)
+ (long "verbosity" <> short 'v' <> metavar "LEVEL" <>
+ value (fromEnum verbosity) <>
+ help "Verbosity level"))
+ <*> strOption (long "template" <> short 't' <> metavar "FILE" <>
+ value template <>
+ help "Template to use for report")
+
+outputOption :: Maybe String -> Mod OptionFields String -> Parser (Maybe String)
+outputOption file m =
+ optional (strOption (m <> metavar "FILE" <> maybe mempty value file))
+
+range :: (Show a, Read a, Ord a) => a -> a -> ReadM a
+range lo hi = do
+ s <- readerAsk
+ case reads s of
+ [(i, "")]
+ | i >= lo && i <= hi -> return i
+ | otherwise -> readerError $ show i ++ " is outside range " ++
+ show (lo,hi)
+ _ -> readerError $ show s ++ " is not a number"
+
+match :: ReadM MatchType
+match = do
+ m <- readerAsk
+ case map toLower m of
+ mm | mm `isPrefixOf` "pfx" -> return Prefix
+ | mm `isPrefixOf` "prefix" -> return Prefix
+ | mm `isPrefixOf` "pattern" -> return Pattern
+ | mm `isPrefixOf` "ipattern" -> return IPattern
+ | otherwise -> readerError $
+ show m ++ " is not a known match type"
+ ++ "Try \"prefix\", \"pattern\", \"ipattern\"."
+
+regressParams :: ReadM ([String], String)
+regressParams = do
+ m <- readerAsk
+ let repl ',' = ' '
+ repl c = c
+ tidy = reverse . dropWhile isSpace . reverse . dropWhile isSpace
+ (r,ps) = break (==':') m
+ when (null r) $
+ readerError "no responder specified"
+ when (null ps) $
+ readerError "no predictors specified"
+ let ret = (words . map repl . drop 1 $ ps, tidy r)
+ either readerError (const (return ret)) $ uncurry validateAccessors ret
+
+-- | Flesh out a command line parser.
+describe :: Config -> ParserInfo Mode
+describe cfg = info (helper <*> parseWith cfg) $
+ header ("Microbenchmark suite - " <> versionInfo) <>
+ fullDesc <>
+ footerDoc (unChunk regressionHelp)
+
+-- | A string describing the version of this benchmark (really, the
+-- version of gauge that was used to build it).
+versionInfo :: String
+versionInfo = "built with gauge " <> showVersion version
+
+-- We sort not by name, but by likely frequency of use.
+regressionHelp :: Chunk Doc
+regressionHelp =
+ fmap (text "Regression metrics (for use with --regress):" .$.) $
+ tabulate [(text n,text d) | (n,(_,d)) <- map f measureKeys]
+ where f k = (k, measureAccessors M.! k)
diff --git a/Gauge/Measurement.hs b/Gauge/Measurement.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Measurement.hs
@@ -0,0 +1,415 @@
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE Trustworthy #-}
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE BangPatterns, CPP, ForeignFunctionInterface,
+ ScopedTypeVariables #-}
+
+#if MIN_VERSION_base(4,10,0)
+-- Disable deprecation warnings for now until we remove the use of getGCStats
+-- and applyGCStats for good
+{-# OPTIONS_GHC -Wno-deprecations #-}
+#endif
+
+-- |
+-- Module : Gauge.Measurement
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Benchmark measurement code.
+
+module Gauge.Measurement
+ (
+ initializeTime
+ , getTime
+ , getCPUTime
+ , getCycles
+ , getGCStatistics
+ , GCStatistics(..)
+ , secs
+ , measure
+ , runBenchmark
+ , runBenchmarkable
+ , runBenchmarkable_
+ , measured
+ , applyGCStatistics
+ , threshold
+ -- * Deprecated
+ , getGCStats
+ , applyGCStats
+ ) where
+
+import Gauge.Types (Benchmarkable(..), Measured(..))
+import Control.Applicative ((<*))
+import Control.DeepSeq (NFData(rnf))
+import Control.Exception (finally,evaluate)
+import Data.Data (Data, Typeable)
+import Data.Int (Int64)
+import Data.List (unfoldr)
+import Data.Word (Word64)
+import GHC.Generics (Generic)
+import GHC.Stats (GCStats(..))
+#if MIN_VERSION_base(4,10,0)
+import GHC.Stats (RTSStats(..), GCDetails(..))
+#endif
+import System.Mem (performGC)
+import Text.Printf (printf)
+import qualified Control.Exception as Exc
+import qualified Data.Vector as V
+import qualified GHC.Stats as Stats
+
+-- | Statistics about memory usage and the garbage collector. Apart from
+-- 'gcStatsCurrentBytesUsed' and 'gcStatsCurrentBytesSlop' all are cumulative values since
+-- the program started.
+--
+-- 'GCStatistics' is cargo-culted from the 'GCStats' data type that "GHC.Stats"
+-- has. Since 'GCStats' was marked as deprecated and will be removed in GHC 8.4,
+-- we use 'GCStatistics' to provide a backwards-compatible view of GC statistics.
+data GCStatistics = GCStatistics
+ { -- | Total number of bytes allocated
+ gcStatsBytesAllocated :: !Int64
+ -- | Number of garbage collections performed (any generation, major and
+ -- minor)
+ , gcStatsNumGcs :: !Int64
+ -- | Maximum number of live bytes seen so far
+ , gcStatsMaxBytesUsed :: !Int64
+ -- | Number of byte usage samples taken, or equivalently
+ -- the number of major GCs performed.
+ , gcStatsNumByteUsageSamples :: !Int64
+ -- | Sum of all byte usage samples, can be used with
+ -- 'gcStatsNumByteUsageSamples' to calculate averages with
+ -- arbitrary weighting (if you are sampling this record multiple
+ -- times).
+ , gcStatsCumulativeBytesUsed :: !Int64
+ -- | Number of bytes copied during GC
+ , gcStatsBytesCopied :: !Int64
+ -- | Number of live bytes at the end of the last major GC
+ , gcStatsCurrentBytesUsed :: !Int64
+ -- | Current number of bytes lost to slop
+ , gcStatsCurrentBytesSlop :: !Int64
+ -- | Maximum number of bytes lost to slop at any one time so far
+ , gcStatsMaxBytesSlop :: !Int64
+ -- | Maximum number of megabytes allocated
+ , gcStatsPeakMegabytesAllocated :: !Int64
+ -- | CPU time spent running mutator threads. This does not include
+ -- any profiling overhead or initialization.
+ , gcStatsMutatorCpuSeconds :: !Double
+
+ -- | Wall clock time spent running mutator threads. This does not
+ -- include initialization.
+ , gcStatsMutatorWallSeconds :: !Double
+ -- | CPU time spent running GC
+ , gcStatsGcCpuSeconds :: !Double
+ -- | Wall clock time spent running GC
+ , gcStatsGcWallSeconds :: !Double
+ -- | Total CPU time elapsed since program start
+ , gcStatsCpuSeconds :: !Double
+ -- | Total wall clock time elapsed since start
+ , gcStatsWallSeconds :: !Double
+ } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+-- | Try to get GC statistics, bearing in mind that the GHC runtime
+-- will throw an exception if statistics collection was not enabled
+-- using \"@+RTS -T@\".
+{-# DEPRECATED getGCStats
+ ["GCStats has been deprecated in GHC 8.2. As a consequence,",
+ "getGCStats has also been deprecated in favor of getGCStatistics.",
+ "getGCStats will be removed in the next major criterion release."] #-}
+getGCStats :: IO (Maybe GCStats)
+getGCStats =
+ (Just `fmap` Stats.getGCStats) `Exc.catch` \(_::Exc.SomeException) ->
+ return Nothing
+
+-- | Try to get GC statistics, bearing in mind that the GHC runtime
+-- will throw an exception if statistics collection was not enabled
+-- using \"@+RTS -T@\".
+getGCStatistics :: IO (Maybe GCStatistics)
+#if MIN_VERSION_base(4,10,0)
+-- Use RTSStats/GCDetails to gather GC stats
+getGCStatistics = do
+ stats <- Stats.getRTSStats
+ let gcdetails :: Stats.GCDetails
+ gcdetails = gc stats
+
+ nsToSecs :: Int64 -> Double
+ nsToSecs ns = fromIntegral ns * 1.0E-9
+
+ return $ Just GCStatistics {
+ gcStatsBytesAllocated = fromIntegral $ gcdetails_allocated_bytes gcdetails
+ , gcStatsNumGcs = fromIntegral $ gcs stats
+ , gcStatsMaxBytesUsed = fromIntegral $ max_live_bytes stats
+ , gcStatsNumByteUsageSamples = fromIntegral $ major_gcs stats
+ , gcStatsCumulativeBytesUsed = fromIntegral $ cumulative_live_bytes stats
+ , gcStatsBytesCopied = fromIntegral $ gcdetails_copied_bytes gcdetails
+ , gcStatsCurrentBytesUsed = fromIntegral $ gcdetails_live_bytes gcdetails
+ , gcStatsCurrentBytesSlop = fromIntegral $ gcdetails_slop_bytes gcdetails
+ , gcStatsMaxBytesSlop = fromIntegral $ max_slop_bytes stats
+ , gcStatsPeakMegabytesAllocated = fromIntegral (max_mem_in_use_bytes stats) `quot` (1024*1024)
+ , gcStatsMutatorCpuSeconds = nsToSecs $ mutator_cpu_ns stats
+ , gcStatsMutatorWallSeconds = nsToSecs $ mutator_elapsed_ns stats
+ , gcStatsGcCpuSeconds = nsToSecs $ gcdetails_cpu_ns gcdetails
+ , gcStatsGcWallSeconds = nsToSecs $ gcdetails_elapsed_ns gcdetails
+ , gcStatsCpuSeconds = nsToSecs $ cpu_ns stats
+ , gcStatsWallSeconds = nsToSecs $ elapsed_ns stats
+ }
+ `Exc.catch`
+ \(_::Exc.SomeException) -> return Nothing
+#else
+-- Use the old GCStats type to gather GC stats
+getGCStatistics = do
+ stats <- Stats.getGCStats
+ return $ Just GCStatistics {
+ gcStatsBytesAllocated = bytesAllocated stats
+ , gcStatsNumGcs = numGcs stats
+ , gcStatsMaxBytesUsed = maxBytesUsed stats
+ , gcStatsNumByteUsageSamples = numByteUsageSamples stats
+ , gcStatsCumulativeBytesUsed = cumulativeBytesUsed stats
+ , gcStatsBytesCopied = bytesCopied stats
+ , gcStatsCurrentBytesUsed = currentBytesUsed stats
+ , gcStatsCurrentBytesSlop = currentBytesSlop stats
+ , gcStatsMaxBytesSlop = maxBytesSlop stats
+ , gcStatsPeakMegabytesAllocated = peakMegabytesAllocated stats
+ , gcStatsMutatorCpuSeconds = mutatorCpuSeconds stats
+ , gcStatsMutatorWallSeconds = mutatorWallSeconds stats
+ , gcStatsGcCpuSeconds = gcCpuSeconds stats
+ , gcStatsGcWallSeconds = gcWallSeconds stats
+ , gcStatsCpuSeconds = cpuSeconds stats
+ , gcStatsWallSeconds = wallSeconds stats
+ }
+ `Exc.catch`
+ \(_::Exc.SomeException) -> return Nothing
+#endif
+
+-- | Measure the execution of a benchmark a given number of times.
+measure :: Benchmarkable -- ^ Operation to benchmark.
+ -> Int64 -- ^ Number of iterations.
+ -> IO (Measured, Double)
+measure bm iters = runBenchmarkable bm iters addResults $ \act -> do
+ startStats <- getGCStatistics
+ startTime <- getTime
+ startCpuTime <- getCPUTime
+ startCycles <- getCycles
+ act
+ endTime <- getTime
+ endCpuTime <- getCPUTime
+ endCycles <- getCycles
+ endStats <- getGCStatistics
+ let !m = applyGCStatistics endStats startStats $ measured {
+ measTime = max 0 (endTime - startTime)
+ , measCpuTime = max 0 (endCpuTime - startCpuTime)
+ , measCycles = max 0 (fromIntegral (endCycles - startCycles))
+ , measIters = iters
+ }
+ return (m, endTime)
+ where
+ addResults :: (Measured, Double) -> (Measured, Double) -> (Measured, Double)
+ addResults (!m1, !d1) (!m2, !d2) = (m3, d1 + d2)
+ where
+ add f = f m1 + f m2
+
+ m3 = Measured
+ { measTime = add measTime
+ , measCpuTime = add measCpuTime
+ , measCycles = add measCycles
+ , measIters = add measIters
+
+ , measAllocated = add measAllocated
+ , measNumGcs = add measNumGcs
+ , measBytesCopied = add measBytesCopied
+ , measMutatorWallSeconds = add measMutatorWallSeconds
+ , measMutatorCpuSeconds = add measMutatorCpuSeconds
+ , measGcWallSeconds = add measGcWallSeconds
+ , measGcCpuSeconds = add measGcCpuSeconds
+ }
+{-# INLINE measure #-}
+
+-- | The amount of time a benchmark must run for in order for us to
+-- have some trust in the raw measurement.
+--
+-- We set this threshold so that we can generate enough data to later
+-- perform meaningful statistical analyses.
+--
+-- The threshold is 30 milliseconds. One use of 'runBenchmark' must
+-- accumulate more than 300 milliseconds of total measurements above
+-- this threshold before it will finish.
+threshold :: Double
+threshold = 0.03
+{-# INLINE threshold #-}
+
+runBenchmarkable :: Benchmarkable -> Int64 -> (a -> a -> a) -> (IO () -> IO a) -> IO a
+runBenchmarkable Benchmarkable{..} i comb f
+ | perRun = work >>= go (i - 1)
+ | otherwise = work
+ where
+ go 0 result = return result
+ go !n !result = work >>= go (n - 1) . comb result
+
+ count | perRun = 1
+ | otherwise = i
+
+ work = do
+ env <- allocEnv count
+ let clean = cleanEnv count env
+ run = runRepeatedly env count
+
+ clean `seq` run `seq` evaluate $ rnf env
+
+ performGC
+ f run `finally` clean <* performGC
+ {-# INLINE work #-}
+{-# INLINE runBenchmarkable #-}
+
+runBenchmarkable_ :: Benchmarkable -> Int64 -> IO ()
+runBenchmarkable_ bm i = runBenchmarkable bm i (\() () -> ()) id
+{-# INLINE runBenchmarkable_ #-}
+
+-- | Run a single benchmark, and return measurements collected while
+-- executing it, along with the amount of time the measurement process
+-- took.
+runBenchmark :: Benchmarkable
+ -> Double
+ -- ^ Lower bound on how long the benchmarking process
+ -- should take. In practice, this time limit may be
+ -- exceeded in order to generate enough data to perform
+ -- meaningful statistical analyses.
+ -> IO (V.Vector Measured, Double)
+runBenchmark bm timeLimit = do
+ runBenchmarkable_ bm 1
+ start <- performGC >> getTime
+ let loop [] !_ !_ _ = error "unpossible!"
+ loop (iters:niters) prev count acc = do
+ (m, endTime) <- measure bm iters
+ let overThresh = max 0 (measTime m - threshold) + prev
+ -- We try to honour the time limit, but we also have more
+ -- important constraints:
+ --
+ -- We must generate enough data that bootstrapping won't
+ -- simply crash.
+ --
+ -- We need to generate enough measurements that have long
+ -- spans of execution to outweigh the (rather high) cost of
+ -- measurement.
+ if endTime - start >= timeLimit &&
+ overThresh > threshold * 10 &&
+ count >= (4 :: Int)
+ then do
+ let !v = V.reverse (V.fromList acc)
+ return (v, endTime - start)
+ else loop niters overThresh (count+1) (m:acc)
+ loop (squish (unfoldr series 1)) 0 0 []
+
+-- Our series starts its growth very slowly when we begin at 1, so we
+-- eliminate repeated values.
+squish :: (Eq a) => [a] -> [a]
+squish ys = foldr go [] ys
+ where go x xs = x : dropWhile (==x) xs
+
+series :: Double -> Maybe (Int64, Double)
+series k = Just (truncate l, l)
+ where l = k * 1.05
+
+-- | An empty structure.
+measured :: Measured
+measured = Measured {
+ measTime = 0
+ , measCpuTime = 0
+ , measCycles = 0
+ , measIters = 0
+
+ , measAllocated = minBound
+ , measNumGcs = minBound
+ , measBytesCopied = minBound
+ , measMutatorWallSeconds = bad
+ , measMutatorCpuSeconds = bad
+ , measGcWallSeconds = bad
+ , measGcCpuSeconds = bad
+ } where bad = -1/0
+
+-- | Apply the difference between two sets of GC statistics to a
+-- measurement.
+{-# DEPRECATED applyGCStats
+ ["GCStats has been deprecated in GHC 8.2. As a consequence,",
+ "applyGCStats has also been deprecated in favor of applyGCStatistics.",
+ "applyGCStats will be removed in the next major criterion release."] #-}
+applyGCStats :: Maybe GCStats
+ -- ^ Statistics gathered at the __end__ of a run.
+ -> Maybe GCStats
+ -- ^ Statistics gathered at the __beginning__ of a run.
+ -> Measured
+ -- ^ Value to \"modify\".
+ -> Measured
+applyGCStats (Just end) (Just start) m = m {
+ measAllocated = diff bytesAllocated
+ , measNumGcs = diff numGcs
+ , measBytesCopied = diff bytesCopied
+ , measMutatorWallSeconds = diff mutatorWallSeconds
+ , measMutatorCpuSeconds = diff mutatorCpuSeconds
+ , measGcWallSeconds = diff gcWallSeconds
+ , measGcCpuSeconds = diff gcCpuSeconds
+ } where diff f = f end - f start
+applyGCStats _ _ m = m
+
+-- | Apply the difference between two sets of GC statistics to a
+-- measurement.
+applyGCStatistics :: Maybe GCStatistics
+ -- ^ Statistics gathered at the __end__ of a run.
+ -> Maybe GCStatistics
+ -- ^ Statistics gathered at the __beginning__ of a run.
+ -> Measured
+ -- ^ Value to \"modify\".
+ -> Measured
+applyGCStatistics (Just end) (Just start) m = m {
+ measAllocated = diff gcStatsBytesAllocated
+ , measNumGcs = diff gcStatsNumGcs
+ , measBytesCopied = diff gcStatsBytesCopied
+ , measMutatorWallSeconds = diff gcStatsMutatorWallSeconds
+ , measMutatorCpuSeconds = diff gcStatsMutatorCpuSeconds
+ , measGcWallSeconds = diff gcStatsGcWallSeconds
+ , measGcCpuSeconds = diff gcStatsGcCpuSeconds
+ } where diff f = f end - f start
+applyGCStatistics _ _ m = m
+
+-- | Convert a number of seconds to a string. The string will consist
+-- of four decimal places, followed by a short description of the time
+-- units.
+secs :: Double -> String
+secs k
+ | k < 0 = '-' : secs (-k)
+ | k >= 1 = k `with` "s"
+ | k >= 1e-3 = (k*1e3) `with` "ms"
+#ifdef mingw32_HOST_OS
+ | k >= 1e-6 = (k*1e6) `with` "us"
+#else
+ | k >= 1e-6 = (k*1e6) `with` "μs"
+#endif
+ | k >= 1e-9 = (k*1e9) `with` "ns"
+ | k >= 1e-12 = (k*1e12) `with` "ps"
+ | k >= 1e-15 = (k*1e15) `with` "fs"
+ | k >= 1e-18 = (k*1e18) `with` "as"
+ | otherwise = printf "%g s" k
+ where with (t :: Double) (u :: String)
+ | t >= 1e9 = printf "%.4g %s" t u
+ | t >= 1e3 = printf "%.0f %s" t u
+ | t >= 1e2 = printf "%.1f %s" t u
+ | t >= 1e1 = printf "%.2f %s" t u
+ | otherwise = printf "%.3f %s" t u
+
+-- | Set up time measurement.
+foreign import ccall unsafe "criterion_inittime" initializeTime :: IO ()
+
+-- | Read the CPU cycle counter.
+foreign import ccall unsafe "criterion_rdtsc" getCycles :: IO Word64
+
+-- | Return the current wallclock time, in seconds since some
+-- arbitrary time.
+--
+-- You /must/ call 'initializeTime' once before calling this function!
+foreign import ccall unsafe "criterion_gettime" getTime :: IO Double
+
+-- | Return the amount of elapsed CPU time, combining user and kernel
+-- (system) time into a single measure.
+foreign import ccall unsafe "criterion_getcputime" getCPUTime :: IO Double
diff --git a/Gauge/Monad.hs b/Gauge/Monad.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Monad.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE Trustworthy #-}
+{-# LANGUAGE TypeFamilies #-}
+-- |
+-- Module : Gauge.Monad
+-- Copyright : (c) 2009 Neil Brown
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- The environment in which most criterion code executes.
+module Gauge.Monad
+ (
+ Gauge
+ , withConfig
+ , getGen
+ , getOverhead
+ ) where
+
+import Foundation.Monad
+import Foundation.Monad.Reader
+import Control.Monad (when)
+import Gauge.Measurement (measure, runBenchmark, secs)
+import Gauge.Monad.Internal (Gauge(..), Crit(..))
+import Gauge.Types hiding (measure)
+import Data.IORef (IORef, newIORef, readIORef, writeIORef)
+import Statistics.Regression (olsRegress)
+import System.Random.MWC (GenIO, createSystemRandom)
+import qualified Data.Vector.Generic as G
+
+-- | Run a 'Gauge' action with the given 'Config'.
+withConfig :: Config -> Gauge a -> IO a
+withConfig cfg (Gauge act) = do
+ g <- newIORef Nothing
+ o <- newIORef Nothing
+ runReaderT act (Crit cfg g o)
+
+-- | Return a random number generator, creating one if necessary.
+--
+-- This is not currently thread-safe, but in a harmless way (we might
+-- call 'createSystemRandom' more than once if multiple threads race).
+getGen :: Gauge GenIO
+getGen = memoise gen createSystemRandom
+
+-- | Return an estimate of the measurement overhead.
+getOverhead :: Gauge Double
+getOverhead = do
+ verbose <- ((== Verbose) . verbosity) <$> ask
+ memoise overhead $ do
+ (meas,_) <- runBenchmark (whnfIO $ measure (whnfIO $ return ()) 1) 1
+ let metric get = G.convert . G.map get $ meas
+ let o = G.head . fst $
+ olsRegress [metric (fromIntegral . measIters)] (metric measTime)
+ when verbose . liftIO $
+ putStrLn $ "measurement overhead " ++ secs o
+ return o
+
+-- | Memoise the result of an 'IO' action.
+--
+-- This is not currently thread-safe, but hopefully in a harmless way.
+-- We might call the given action more than once if multiple threads
+-- race, so our caller's job is to write actions that can be run
+-- multiple times safely.
+memoise :: (Crit -> IORef (Maybe a)) -> IO a -> Gauge a
+memoise ref generate = do
+ r <- Gauge (ref <$> ask)
+ liftIO $ do
+ mv <- readIORef r
+ case mv of
+ Just rv -> return rv
+ Nothing -> do
+ rv <- generate
+ writeIORef r (Just rv)
+ return rv
diff --git a/Gauge/Monad/ExceptT.hs b/Gauge/Monad/ExceptT.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Monad/ExceptT.hs
@@ -0,0 +1,56 @@
+{-# LANGUAGE TypeFamilies #-}
+module Gauge.Monad.ExceptT
+ ( ExceptT(..)
+ , finally
+ -- , try
+ ) where
+
+import Foundation.Monad
+import Foundation.Monad.Reader
+
+newtype ExceptT e m a = ExceptT { runExceptT :: m (Either e a) }
+
+instance (Functor m) => Functor (ExceptT e m) where
+ fmap f = ExceptT . fmap (fmap f) . runExceptT
+
+instance (Functor m, Monad m) => Applicative (ExceptT e m) where
+ pure a = ExceptT $ return (Right a)
+ ExceptT f <*> ExceptT v = ExceptT $ do
+ mf <- f
+ case mf of
+ Left e -> return (Left e)
+ Right k -> do
+ mv <- v
+ case mv of
+ Left e -> return (Left e)
+ Right x -> return (Right (k x))
+
+instance Monad m => MonadFailure (ExceptT e m) where
+ type Failure (ExceptT e m) = e
+ mFail = ExceptT . pure . Left
+
+instance (Monad m) => Monad (ExceptT e m) where
+ return a = ExceptT $ return (Right a)
+ m >>= k = ExceptT $ do
+ a <- runExceptT m
+ case a of
+ Left e -> return (Left e)
+ Right x -> runExceptT (k x)
+ fail = ExceptT . fail
+
+instance MonadReader m => MonadReader (ExceptT e m) where
+ type ReaderContext (ExceptT e m) = ReaderContext m
+ ask = ExceptT (Right <$> ask)
+
+instance MonadTrans (ExceptT e) where
+ lift f = ExceptT (Right <$> f)
+
+instance MonadIO m => MonadIO (ExceptT e m) where
+ liftIO f = ExceptT (Right <$> liftIO f)
+
+finally :: MonadBracket m => m a -> m b -> m a
+finally f g = generalBracket (pure ()) (\() a -> g >> pure a) (\() _ -> g) (const f)
+
+--try :: (MonadCatch, Exception e) => m a -> m (Either e a)
+--try a = catch (a >>= \ v -> return (Right v)) (\e -> return (Left e))
+
diff --git a/Gauge/Monad/Internal.hs b/Gauge/Monad/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Monad/Internal.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# OPTIONS_GHC -funbox-strict-fields #-}
+
+-- |
+-- Module : Gauge.Monad.Internal
+-- Copyright : (c) 2009 Neil Brown
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- The environment in which most criterion code executes.
+module Gauge.Monad.Internal
+ (
+ Gauge(..)
+ , Crit(..)
+ ) where
+
+-- Temporary: to support pre-AMP GHC 7.8.4:
+import Control.Applicative
+
+import Foundation.Monad
+import Foundation.Monad.Reader
+import Gauge.Types (Config)
+import Data.IORef (IORef)
+import System.Random.MWC (GenIO)
+import Prelude
+
+data Crit = Crit {
+ config :: !Config
+ , gen :: !(IORef (Maybe GenIO))
+ , overhead :: !(IORef (Maybe Double))
+ }
+
+-- | The monad in which most criterion code executes.
+newtype Gauge a = Gauge {
+ runGauge :: ReaderT Crit IO a
+ } deriving (Functor, Applicative, Monad, MonadIO, MonadThrow, MonadCatch) -- , MonadBracket)
+
+instance MonadReader Gauge where
+ type ReaderContext Gauge = Config
+ ask = config `fmap` Gauge ask
+
+instance MonadBracket Gauge where
+ generalBracket acq cleanup cleanupExcept innerAction = Gauge $ do
+ c <- ask
+ lift $ generalBracket (runReaderT (runGauge acq) c)
+ (\a b -> runReaderT (runGauge (cleanup a b)) c)
+ (\a exn -> runReaderT (runGauge (cleanupExcept a exn)) c)
+ (\a -> runReaderT (runGauge (innerAction a)) c)
diff --git a/Gauge/Types.hs b/Gauge/Types.hs
new file mode 100644
--- /dev/null
+++ b/Gauge/Types.hs
@@ -0,0 +1,699 @@
+{-# LANGUAGE Trustworthy #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric, GADTs, RecordWildCards #-}
+{-# OPTIONS_GHC -funbox-strict-fields #-}
+
+-- |
+-- Module : Gauge.Types
+-- Copyright : (c) 2009-2014 Bryan O'Sullivan
+--
+-- License : BSD-style
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : GHC
+--
+-- Types for benchmarking.
+--
+-- The core type is 'Benchmarkable', which admits both pure functions
+-- and 'IO' actions.
+--
+-- For a pure function of type @a -> b@, the benchmarking harness
+-- calls this function repeatedly, each time with a different 'Int64'
+-- argument (the number of times to run the function in a loop), and
+-- reduces the result the function returns to weak head normal form.
+--
+-- For an action of type @IO a@, the benchmarking harness calls the
+-- action repeatedly, but does not reduce the result.
+
+module Gauge.Types
+ (
+ -- * Configuration
+ Config(..)
+ , Verbosity(..)
+ -- * Benchmark descriptions
+ , Benchmarkable(..)
+ , Benchmark(..)
+ -- * Measurements
+ , Measured(..)
+ , fromInt
+ , toInt
+ , fromDouble
+ , toDouble
+ , measureAccessors
+ , measureKeys
+ , measure
+ , rescale
+ -- * Benchmark construction
+ , env
+ , envWithCleanup
+ , perBatchEnv
+ , perBatchEnvWithCleanup
+ , perRunEnv
+ , perRunEnvWithCleanup
+ , toBenchmarkable
+ , bench
+ , bgroup
+ , addPrefix
+ , benchNames
+ -- ** Evaluation control
+ , whnf
+ , nf
+ , nfIO
+ , whnfIO
+ -- * Result types
+ , Outliers(..)
+ , OutlierEffect(..)
+ , OutlierVariance(..)
+ , Regression(..)
+ , KDE(..)
+ , Report(..)
+ , SampleAnalysis(..)
+ , DataRecord(..)
+ ) where
+
+-- Temporary: to support pre-AMP GHC 7.8.4:
+import Control.Applicative
+import Data.Monoid
+
+import Control.DeepSeq (NFData(rnf))
+import Control.Exception (evaluate)
+import Data.Data (Data, Typeable)
+import Data.Int (Int64)
+import Data.Map (Map, fromList)
+import GHC.Generics (Generic)
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as U
+import qualified Statistics.Types as St
+import Prelude
+
+-- | Control the amount of information displayed.
+data Verbosity = Quiet
+ | Normal
+ | Verbose
+ deriving (Eq, Ord, Bounded, Enum, Read, Show, Typeable, Data,
+ Generic)
+
+-- | Top-level benchmarking configuration.
+data Config = Config {
+ confInterval :: St.CL Double
+ -- ^ Confidence interval for bootstrap estimation (greater than
+ -- 0, less than 1).
+ , forceGC :: Bool
+ -- ^ /Obsolete, unused/. This option used to force garbage
+ -- collection between every benchmark run, but it no longer has
+ -- an effect (we now unconditionally force garbage collection).
+ -- This option remains solely for backwards API compatibility.
+ , timeLimit :: Double
+ -- ^ Number of seconds to run a single benchmark. (In practice,
+ -- execution time will very slightly exceed this limit.)
+ , resamples :: Int
+ -- ^ Number of resamples to perform when bootstrapping.
+ , regressions :: [([String], String)]
+ -- ^ Regressions to perform.
+ , rawDataFile :: Maybe FilePath
+ -- ^ File to write binary measurement and analysis data to. If
+ -- not specified, this will be a temporary file.
+ , reportFile :: Maybe FilePath
+ -- ^ File to write report output to, with template expanded.
+ , csvFile :: Maybe FilePath
+ -- ^ File to write CSV summary to.
+ , jsonFile :: Maybe FilePath
+ -- ^ File to write JSON-formatted results to.
+ , junitFile :: Maybe FilePath
+ -- ^ File to write JUnit-compatible XML results to.
+ , verbosity :: Verbosity
+ -- ^ Verbosity level to use when running and analysing
+ -- benchmarks.
+ , template :: FilePath
+ -- ^ Template file to use if writing a report.
+ } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+
+-- | A pure function or impure action that can be benchmarked. The
+-- 'Int64' parameter indicates the number of times to run the given
+-- function or action.
+data Benchmarkable = forall a . NFData a =>
+ Benchmarkable
+ { allocEnv :: Int64 -> IO a
+ , cleanEnv :: Int64 -> a -> IO ()
+ , runRepeatedly :: a -> Int64 -> IO ()
+ , perRun :: Bool
+ }
+
+noop :: Monad m => a -> m ()
+noop = const $ return ()
+{-# INLINE noop #-}
+
+-- | Construct a 'Benchmarkable' value from an impure action, where the 'Int64'
+-- parameter indicates the number of times to run the action.
+toBenchmarkable :: (Int64 -> IO ()) -> Benchmarkable
+toBenchmarkable f = Benchmarkable noop (const noop) (const f) False
+{-# INLINE toBenchmarkable #-}
+
+-- | A collection of measurements made while benchmarking.
+--
+-- Measurements related to garbage collection are tagged with __GC__.
+-- They will only be available if a benchmark is run with @\"+RTS
+-- -T\"@.
+--
+-- __Packed storage.__ When GC statistics cannot be collected, GC
+-- values will be set to huge negative values. If a field is labeled
+-- with \"__GC__\" below, use 'fromInt' and 'fromDouble' to safely
+-- convert to \"real\" values.
+data Measured = Measured {
+ measTime :: !Double
+ -- ^ Total wall-clock time elapsed, in seconds.
+ , measCpuTime :: !Double
+ -- ^ Total CPU time elapsed, in seconds. Includes both user and
+ -- kernel (system) time.
+ , measCycles :: !Int64
+ -- ^ Cycles, in unspecified units that may be CPU cycles. (On
+ -- i386 and x86_64, this is measured using the @rdtsc@
+ -- instruction.)
+ , measIters :: !Int64
+ -- ^ Number of loop iterations measured.
+
+ , measAllocated :: !Int64
+ -- ^ __(GC)__ Number of bytes allocated. Access using 'fromInt'.
+ , measNumGcs :: !Int64
+ -- ^ __(GC)__ Number of garbage collections performed. Access
+ -- using 'fromInt'.
+ , measBytesCopied :: !Int64
+ -- ^ __(GC)__ Number of bytes copied during garbage collection.
+ -- Access using 'fromInt'.
+ , measMutatorWallSeconds :: !Double
+ -- ^ __(GC)__ Wall-clock time spent doing real work
+ -- (\"mutation\"), as distinct from garbage collection. Access
+ -- using 'fromDouble'.
+ , measMutatorCpuSeconds :: !Double
+ -- ^ __(GC)__ CPU time spent doing real work (\"mutation\"), as
+ -- distinct from garbage collection. Access using 'fromDouble'.
+ , measGcWallSeconds :: !Double
+ -- ^ __(GC)__ Wall-clock time spent doing garbage collection.
+ -- Access using 'fromDouble'.
+ , measGcCpuSeconds :: !Double
+ -- ^ __(GC)__ CPU time spent doing garbage collection. Access
+ -- using 'fromDouble'.
+ } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+instance NFData Measured where
+ rnf Measured{} = ()
+
+-- THIS MUST REFLECT THE ORDER OF FIELDS IN THE DATA TYPE.
+--
+-- The ordering is used by Javascript code to pick out the correct
+-- index into the vector that represents a Measured value in that
+-- world.
+measureAccessors_ :: [(String, (Measured -> Maybe Double, String))]
+measureAccessors_ = [
+ ("time", (Just . measTime,
+ "wall-clock time"))
+ , ("cpuTime", (Just . measCpuTime,
+ "CPU time"))
+ , ("cycles", (Just . fromIntegral . measCycles,
+ "CPU cycles"))
+ , ("iters", (Just . fromIntegral . measIters,
+ "loop iterations"))
+ , ("allocated", (fmap fromIntegral . fromInt . measAllocated,
+ "(+RTS -T) bytes allocated"))
+ , ("numGcs", (fmap fromIntegral . fromInt . measNumGcs,
+ "(+RTS -T) number of garbage collections"))
+ , ("bytesCopied", (fmap fromIntegral . fromInt . measBytesCopied,
+ "(+RTS -T) number of bytes copied during GC"))
+ , ("mutatorWallSeconds", (fromDouble . measMutatorWallSeconds,
+ "(+RTS -T) wall-clock time for mutator threads"))
+ , ("mutatorCpuSeconds", (fromDouble . measMutatorCpuSeconds,
+ "(+RTS -T) CPU time spent running mutator threads"))
+ , ("gcWallSeconds", (fromDouble . measGcWallSeconds,
+ "(+RTS -T) wall-clock time spent doing GC"))
+ , ("gcCpuSeconds", (fromDouble . measGcCpuSeconds,
+ "(+RTS -T) CPU time spent doing GC"))
+ ]
+
+-- | Field names in a 'Measured' record, in the order in which they
+-- appear.
+measureKeys :: [String]
+measureKeys = map fst measureAccessors_
+
+-- | Field names and accessors for a 'Measured' record.
+measureAccessors :: Map String (Measured -> Maybe Double, String)
+measureAccessors = fromList measureAccessors_
+
+-- | Normalise every measurement as if 'measIters' was 1.
+--
+-- ('measIters' itself is left unaffected.)
+rescale :: Measured -> Measured
+rescale m@Measured{..} = m {
+ measTime = d measTime
+ , measCpuTime = d measCpuTime
+ , measCycles = i measCycles
+ -- skip measIters
+ , measNumGcs = i measNumGcs
+ , measBytesCopied = i measBytesCopied
+ , measMutatorWallSeconds = d measMutatorWallSeconds
+ , measMutatorCpuSeconds = d measMutatorCpuSeconds
+ , measGcWallSeconds = d measGcWallSeconds
+ , measGcCpuSeconds = d measGcCpuSeconds
+ } where
+ d k = maybe k (/ iters) (fromDouble k)
+ i k = maybe k (round . (/ iters)) (fromIntegral <$> fromInt k)
+ iters = fromIntegral measIters :: Double
+
+-- | Convert a (possibly unavailable) GC measurement to a true value.
+-- If the measurement is a huge negative number that corresponds to
+-- \"no data\", this will return 'Nothing'.
+fromInt :: Int64 -> Maybe Int64
+fromInt i | i == minBound = Nothing
+ | otherwise = Just i
+
+-- | Convert from a true value back to the packed representation used
+-- for GC measurements.
+toInt :: Maybe Int64 -> Int64
+toInt Nothing = minBound
+toInt (Just i) = i
+
+-- | Convert a (possibly unavailable) GC measurement to a true value.
+-- If the measurement is a huge negative number that corresponds to
+-- \"no data\", this will return 'Nothing'.
+fromDouble :: Double -> Maybe Double
+fromDouble d | isInfinite d || isNaN d = Nothing
+ | otherwise = Just d
+
+-- | Convert from a true value back to the packed representation used
+-- for GC measurements.
+toDouble :: Maybe Double -> Double
+toDouble Nothing = -1/0
+toDouble (Just d) = d
+
+-- | Apply an argument to a function, and evaluate the result to weak
+-- head normal form (WHNF).
+whnf :: (a -> b) -> a -> Benchmarkable
+whnf = pureFunc id
+{-# INLINE whnf #-}
+
+-- | Apply an argument to a function, and evaluate the result to
+-- normal form (NF).
+nf :: NFData b => (a -> b) -> a -> Benchmarkable
+nf = pureFunc rnf
+{-# INLINE nf #-}
+
+pureFunc :: (b -> c) -> (a -> b) -> a -> Benchmarkable
+pureFunc reduce f0 x0 = toBenchmarkable (go f0 x0)
+ where go f x n
+ | n <= 0 = return ()
+ | otherwise = evaluate (reduce (f x)) >> go f x (n-1)
+{-# INLINE pureFunc #-}
+
+-- | Perform an action, then evaluate its result to normal form.
+-- This is particularly useful for forcing a lazy 'IO' action to be
+-- completely performed.
+nfIO :: NFData a => IO a -> Benchmarkable
+nfIO = toBenchmarkable . impure rnf
+{-# INLINE nfIO #-}
+
+-- | Perform an action, then evaluate its result to weak head normal
+-- form (WHNF). This is useful for forcing an 'IO' action whose result
+-- is an expression to be evaluated down to a more useful value.
+whnfIO :: IO a -> Benchmarkable
+whnfIO = toBenchmarkable . impure id
+{-# INLINE whnfIO #-}
+
+impure :: (a -> b) -> IO a -> Int64 -> IO ()
+impure strategy a = go
+ where go n
+ | n <= 0 = return ()
+ | otherwise = a >>= (evaluate . strategy) >> go (n-1)
+{-# INLINE impure #-}
+
+-- | Specification of a collection of benchmarks and environments. A
+-- benchmark may consist of:
+--
+-- * An environment that creates input data for benchmarks, created
+-- with 'env'.
+--
+-- * A single 'Benchmarkable' item with a name, created with 'bench'.
+--
+-- * A (possibly nested) group of 'Benchmark's, created with 'bgroup'.
+data Benchmark where
+ Environment :: NFData env
+ => IO env -> (env -> IO a) -> (env -> Benchmark) -> Benchmark
+ Benchmark :: String -> Benchmarkable -> Benchmark
+ BenchGroup :: String -> [Benchmark] -> Benchmark
+
+-- | Run a benchmark (or collection of benchmarks) in the given
+-- environment. The purpose of an environment is to lazily create
+-- input data to pass to the functions that will be benchmarked.
+--
+-- A common example of environment data is input that is read from a
+-- file. Another is a large data structure constructed in-place.
+--
+-- __Motivation.__ In earlier versions of criterion, all benchmark
+-- inputs were always created when a program started running. By
+-- deferring the creation of an environment when its associated
+-- benchmarks need the its, we avoid two problems that this strategy
+-- caused:
+--
+-- * Memory pressure distorted the results of unrelated benchmarks.
+-- If one benchmark needed e.g. a gigabyte-sized input, it would
+-- force the garbage collector to do extra work when running some
+-- other benchmark that had no use for that input. Since the data
+-- created by an environment is only available when it is in scope,
+-- it should be garbage collected before other benchmarks are run.
+--
+-- * The time cost of generating all needed inputs could be
+-- significant in cases where no inputs (or just a few) were really
+-- needed. This occurred often, for instance when just one out of a
+-- large suite of benchmarks was run, or when a user would list the
+-- collection of benchmarks without running any.
+--
+-- __Creation.__ An environment is created right before its related
+-- benchmarks are run. The 'IO' action that creates the environment
+-- is run, then the newly created environment is evaluated to normal
+-- form (hence the 'NFData' constraint) before being passed to the
+-- function that receives the environment.
+--
+-- __Complex environments.__ If you need to create an environment that
+-- contains multiple values, simply pack the values into a tuple.
+--
+-- __Lazy pattern matching.__ In situations where a \"real\"
+-- environment is not needed, e.g. if a list of benchmark names is
+-- being generated, @undefined@ will be passed to the function that
+-- receives the environment. This avoids the overhead of generating
+-- an environment that will not actually be used.
+--
+-- The function that receives the environment must use lazy pattern
+-- matching to deconstruct the tuple, as use of strict pattern
+-- matching will cause a crash if @undefined@ is passed in.
+--
+-- __Example.__ This program runs benchmarks in an environment that
+-- contains two values. The first value is the contents of a text
+-- file; the second is a string. Pay attention to the use of a lazy
+-- pattern to deconstruct the tuple in the function that returns the
+-- benchmarks to be run.
+--
+-- > setupEnv = do
+-- > let small = replicate 1000 (1 :: Int)
+-- > big <- map length . words <$> readFile "/usr/dict/words"
+-- > return (small, big)
+-- >
+-- > main = defaultMain [
+-- > -- notice the lazy pattern match here!
+-- > env setupEnv $ \ ~(small,big) -> bgroup "main" [
+-- > bgroup "small" [
+-- > bench "length" $ whnf length small
+-- > , bench "length . filter" $ whnf (length . filter (==1)) small
+-- > ]
+-- > , bgroup "big" [
+-- > bench "length" $ whnf length big
+-- > , bench "length . filter" $ whnf (length . filter (==1)) big
+-- > ]
+-- > ] ]
+--
+-- __Discussion.__ The environment created in the example above is
+-- intentionally /not/ ideal. As Haskell's scoping rules suggest, the
+-- variable @big@ is in scope for the benchmarks that use only
+-- @small@. It would be better to create a separate environment for
+-- @big@, so that it will not be kept alive while the unrelated
+-- benchmarks are being run.
+env :: NFData env =>
+ IO env
+ -- ^ Create the environment. The environment will be evaluated to
+ -- normal form before being passed to the benchmark.
+ -> (env -> Benchmark)
+ -- ^ Take the newly created environment and make it available to
+ -- the given benchmarks.
+ -> Benchmark
+env alloc = Environment alloc noop
+
+-- | Same as `env`, but but allows for an additional callback
+-- to clean up the environment. Resource clean up is exception safe, that is,
+-- it runs even if the 'Benchmark' throws an exception.
+envWithCleanup
+ :: NFData env
+ => IO env
+ -- ^ Create the environment. The environment will be evaluated to
+ -- normal form before being passed to the benchmark.
+ -> (env -> IO a)
+ -- ^ Clean up the created environment.
+ -> (env -> Benchmark)
+ -- ^ Take the newly created environment and make it available to
+ -- the given benchmarks.
+ -> Benchmark
+envWithCleanup = Environment
+
+-- | Create a Benchmarkable where a fresh environment is allocated for every
+-- batch of runs of the benchmarkable.
+--
+-- The environment is evaluated to normal form before the benchmark is run.
+--
+-- When using 'whnf', 'whnfIO', etc. Gauge creates a 'Benchmarkable'
+-- whichs runs a batch of @N@ repeat runs of that expressions. Gauge may
+-- run any number of these batches to get accurate measurements. Environments
+-- created by 'env' and 'envWithCleanup', are shared across all these batches
+-- of runs.
+--
+-- This is fine for simple benchmarks on static input, but when benchmarking
+-- IO operations where these operations can modify (and especially grow) the
+-- environment this means that later batches might have their accuracy effected
+-- due to longer, for example, longer garbage collection pauses.
+--
+-- An example: Suppose we want to benchmark writing to a Chan, if we allocate
+-- the Chan using environment and our benchmark consists of @writeChan env ()@,
+-- the contents and thus size of the Chan will grow with every repeat. If
+-- Gauge runs a 1,000 batches of 1,000 repeats, the result is that the
+-- channel will have 999,000 items in it by the time the last batch is run.
+-- Since GHC GC has to copy the live set for every major GC this means our last
+-- set of writes will suffer a lot of noise of the previous repeats.
+--
+-- By allocating a fresh environment for every batch of runs this function
+-- should eliminate this effect.
+perBatchEnv
+ :: (NFData env, NFData b)
+ => (Int64 -> IO env)
+ -- ^ Create an environment for a batch of N runs. The environment will be
+ -- evaluated to normal form before running.
+ -> (env -> IO b)
+ -- ^ Function returning the IO action that should be benchmarked with the
+ -- newly generated environment.
+ -> Benchmarkable
+perBatchEnv alloc = perBatchEnvWithCleanup alloc (const noop)
+
+-- | Same as `perBatchEnv`, but but allows for an additional callback
+-- to clean up the environment. Resource clean up is exception safe, that is,
+-- it runs even if the 'Benchmark' throws an exception.
+perBatchEnvWithCleanup
+ :: (NFData env, NFData b)
+ => (Int64 -> IO env)
+ -- ^ Create an environment for a batch of N runs. The environment will be
+ -- evaluated to normal form before running.
+ -> (Int64 -> env -> IO ())
+ -- ^ Clean up the created environment.
+ -> (env -> IO b)
+ -- ^ Function returning the IO action that should be benchmarked with the
+ -- newly generated environment.
+ -> Benchmarkable
+perBatchEnvWithCleanup alloc clean work
+ = Benchmarkable alloc clean (impure rnf . work) False
+
+-- | Create a Benchmarkable where a fresh environment is allocated for every
+-- run of the operation to benchmark. This is useful for benchmarking mutable
+-- operations that need a fresh environment, such as sorting a mutable Vector.
+--
+-- As with 'env' and 'perBatchEnv' the environment is evaluated to normal form
+-- before the benchmark is run.
+--
+-- This introduces extra noise and result in reduce accuracy compared to other
+-- Gauge benchmarks. But allows easier benchmarking for mutable operations
+-- than was previously possible.
+perRunEnv
+ :: (NFData env, NFData b)
+ => IO env
+ -- ^ Action that creates the environment for a single run.
+ -> (env -> IO b)
+ -- ^ Function returning the IO action that should be benchmarked with the
+ -- newly genereted environment.
+ -> Benchmarkable
+perRunEnv alloc = perRunEnvWithCleanup alloc noop
+
+-- | Same as `perRunEnv`, but but allows for an additional callback
+-- to clean up the environment. Resource clean up is exception safe, that is,
+-- it runs even if the 'Benchmark' throws an exception.
+perRunEnvWithCleanup
+ :: (NFData env, NFData b)
+ => IO env
+ -- ^ Action that creates the environment for a single run.
+ -> (env -> IO ())
+ -- ^ Clean up the created environment.
+ -> (env -> IO b)
+ -- ^ Function returning the IO action that should be benchmarked with the
+ -- newly genereted environment.
+ -> Benchmarkable
+perRunEnvWithCleanup alloc clean work = bm { perRun = True }
+ where
+ bm = perBatchEnvWithCleanup (const alloc) (const clean) work
+
+-- | Create a single benchmark.
+bench :: String -- ^ A name to identify the benchmark.
+ -> Benchmarkable -- ^ An activity to be benchmarked.
+ -> Benchmark
+bench = Benchmark
+
+-- | Group several benchmarks together under a common name.
+bgroup :: String -- ^ A name to identify the group of benchmarks.
+ -> [Benchmark] -- ^ Benchmarks to group under this name.
+ -> Benchmark
+bgroup = BenchGroup
+
+-- | Add the given prefix to a name. If the prefix is empty, the name
+-- is returned unmodified. Otherwise, the prefix and name are
+-- separated by a @\'\/\'@ character.
+addPrefix :: String -- ^ Prefix.
+ -> String -- ^ Name.
+ -> String
+addPrefix "" desc = desc
+addPrefix pfx desc = pfx ++ '/' : desc
+
+-- | Retrieve the names of all benchmarks. Grouped benchmarks are
+-- prefixed with the name of the group they're in.
+benchNames :: Benchmark -> [String]
+benchNames (Environment _ _ b) = benchNames (b undefined)
+benchNames (Benchmark d _) = [d]
+benchNames (BenchGroup d bs) = map (addPrefix d) . concatMap benchNames $ bs
+
+instance Show Benchmark where
+ show (Environment _ _ b) = "Environment _ _" ++ show (b undefined)
+ show (Benchmark d _) = "Benchmark " ++ show d
+ show (BenchGroup d _) = "BenchGroup " ++ show d
+
+measure :: (U.Unbox a) => (Measured -> a) -> V.Vector Measured -> U.Vector a
+measure f v = U.convert . V.map f $ v
+
+-- | Outliers from sample data, calculated using the boxplot
+-- technique.
+data Outliers = Outliers {
+ samplesSeen :: !Int64
+ , lowSevere :: !Int64
+ -- ^ More than 3 times the interquartile range (IQR) below the
+ -- first quartile.
+ , lowMild :: !Int64
+ -- ^ Between 1.5 and 3 times the IQR below the first quartile.
+ , highMild :: !Int64
+ -- ^ Between 1.5 and 3 times the IQR above the third quartile.
+ , highSevere :: !Int64
+ -- ^ More than 3 times the IQR above the third quartile.
+ } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+instance NFData Outliers
+
+-- | A description of the extent to which outliers in the sample data
+-- affect the sample mean and standard deviation.
+data OutlierEffect = Unaffected -- ^ Less than 1% effect.
+ | Slight -- ^ Between 1% and 10%.
+ | Moderate -- ^ Between 10% and 50%.
+ | Severe -- ^ Above 50% (i.e. measurements
+ -- are useless).
+ deriving (Eq, Ord, Read, Show, Typeable, Data, Generic)
+
+instance NFData OutlierEffect
+
+instance Monoid Outliers where
+ mempty = Outliers 0 0 0 0 0
+ mappend = addOutliers
+
+addOutliers :: Outliers -> Outliers -> Outliers
+addOutliers (Outliers s a b c d) (Outliers t w x y z) =
+ Outliers (s+t) (a+w) (b+x) (c+y) (d+z)
+{-# INLINE addOutliers #-}
+
+-- | Analysis of the extent to which outliers in a sample affect its
+-- standard deviation (and to some extent, its mean).
+data OutlierVariance = OutlierVariance {
+ ovEffect :: OutlierEffect
+ -- ^ Qualitative description of effect.
+ , ovDesc :: String
+ -- ^ Brief textual description of effect.
+ , ovFraction :: Double
+ -- ^ Quantitative description of effect (a fraction between 0 and 1).
+ } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+instance NFData OutlierVariance where
+ rnf OutlierVariance{..} = rnf ovEffect `seq` rnf ovDesc `seq` rnf ovFraction
+
+-- | Results of a linear regression.
+data Regression = Regression {
+ regResponder :: String
+ -- ^ Name of the responding variable.
+ , regCoeffs :: Map String (St.Estimate St.ConfInt Double)
+ -- ^ Map from name to value of predictor coefficients.
+ , regRSquare :: St.Estimate St.ConfInt Double
+ -- ^ R² goodness-of-fit estimate.
+ } deriving (Eq, Read, Show, Typeable, Generic)
+
+instance NFData Regression where
+ rnf Regression{..} =
+ rnf regResponder `seq` rnf regCoeffs `seq` rnf regRSquare
+
+-- | Result of a bootstrap analysis of a non-parametric sample.
+data SampleAnalysis = SampleAnalysis {
+ anRegress :: [Regression]
+ -- ^ Estimates calculated via linear regression.
+ , anOverhead :: Double
+ -- ^ Estimated measurement overhead, in seconds. Estimation is
+ -- performed via linear regression.
+ , anMean :: St.Estimate St.ConfInt Double
+ -- ^ Estimated mean.
+ , anStdDev :: St.Estimate St.ConfInt Double
+ -- ^ Estimated standard deviation.
+ , anOutlierVar :: OutlierVariance
+ -- ^ Description of the effects of outliers on the estimated
+ -- variance.
+ } deriving (Eq, Read, Show, Typeable, Generic)
+
+instance NFData SampleAnalysis where
+ rnf SampleAnalysis{..} =
+ rnf anRegress `seq` rnf anOverhead `seq` rnf anMean `seq`
+ rnf anStdDev `seq` rnf anOutlierVar
+
+-- | Data for a KDE chart of performance.
+data KDE = KDE {
+ kdeType :: String
+ , kdeValues :: U.Vector Double
+ , kdePDF :: U.Vector Double
+ } deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+instance NFData KDE where
+ rnf KDE{..} = rnf kdeType `seq` rnf kdeValues `seq` rnf kdePDF
+
+-- | Report of a sample analysis.
+data Report = Report {
+ reportNumber :: Int
+ -- ^ A simple index indicating that this is the /n/th report.
+ , reportName :: String
+ -- ^ The name of this report.
+ , reportKeys :: [String]
+ -- ^ See 'measureKeys'.
+ , reportMeasured :: V.Vector Measured
+ -- ^ Raw measurements. These are /not/ corrected for the
+ -- estimated measurement overhead that can be found via the
+ -- 'anOverhead' field of 'reportAnalysis'.
+ , reportAnalysis :: SampleAnalysis
+ -- ^ Report analysis.
+ , reportOutliers :: Outliers
+ -- ^ Analysis of outliers.
+ , reportKDEs :: [KDE]
+ -- ^ Data for a KDE of times.
+ } deriving (Eq, Read, Show, Typeable, Generic)
+
+instance NFData Report where
+ rnf Report{..} =
+ rnf reportNumber `seq` rnf reportName `seq` rnf reportKeys `seq`
+ rnf reportMeasured `seq` rnf reportAnalysis `seq` rnf reportOutliers `seq`
+ rnf reportKDEs
+
+data DataRecord = Measurement Int String (V.Vector Measured)
+ | Analysed Report
+ deriving (Eq, Read, Show, Typeable, Generic)
+
+instance NFData DataRecord where
+ rnf (Measurement i n v) = rnf i `seq` rnf n `seq` rnf v
+ rnf (Analysed r) = rnf r
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) 2009, 2010 Bryan O'Sullivan
+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.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.markdown b/README.markdown
new file mode 100644
--- /dev/null
+++ b/README.markdown
@@ -0,0 +1,1 @@
+# Gauge: a clone of criterion
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/cbits/cycles.c b/cbits/cycles.c
new file mode 100644
--- /dev/null
+++ b/cbits/cycles.c
@@ -0,0 +1,57 @@
+#include "Rts.h"
+
+#if x86_64_HOST_ARCH || i386_HOST_ARCH
+
+StgWord64 criterion_rdtsc(void)
+{
+ StgWord32 hi, lo;
+ __asm__ __volatile__ ("rdtsc" : "=a"(lo), "=d"(hi));
+ return ((StgWord64) lo) | (((StgWord64) hi)<<32);
+}
+
+#elif linux_HOST_OS
+
+/*
+ * This should work on all Linux.
+ *
+ * Technique by Austin Seipp found here:
+ *
+ * http://neocontra.blogspot.com/2013/05/user-mode-performance-counters-for.html
+ */
+
+#include
+#include
+#include
+
+static int fddev = -1;
+__attribute__((constructor))
+static void
+init(void)
+{
+ static struct perf_event_attr attr;
+ attr.type = PERF_TYPE_HARDWARE;
+ attr.config = PERF_COUNT_HW_CPU_CYCLES;
+ fddev = syscall (__NR_perf_event_open, &attr, 0, -1, -1, 0);
+}
+
+__attribute__((destructor))
+static void
+fini(void)
+{
+ close(fddev);
+}
+
+StgWord64
+criterion_rdtsc (void)
+{
+ StgWord64 result = 0;
+ if (read (fddev, &result, sizeof(result)) < sizeof(result))
+ return 0;
+ return result;
+}
+
+#else
+
+#error Unsupported OS/architecture/compiler!
+
+#endif
diff --git a/cbits/time-osx.c b/cbits/time-osx.c
new file mode 100644
--- /dev/null
+++ b/cbits/time-osx.c
@@ -0,0 +1,35 @@
+#include
+#include
+
+static mach_timebase_info_data_t timebase_info;
+static double timebase_recip;
+
+void criterion_inittime(void)
+{
+ if (timebase_recip == 0) {
+ mach_timebase_info(&timebase_info);
+ timebase_recip = (timebase_info.denom / timebase_info.numer) / 1e9;
+ }
+}
+
+double criterion_gettime(void)
+{
+ return mach_absolute_time() * timebase_recip;
+}
+
+static double to_double(time_value_t time)
+{
+ return time.seconds + time.microseconds / 1e6;
+}
+
+double criterion_getcputime(void)
+{
+ struct task_thread_times_info thread_info_data;
+ mach_msg_type_number_t thread_info_count = TASK_THREAD_TIMES_INFO_COUNT;
+ kern_return_t kr = task_info(mach_task_self(),
+ TASK_THREAD_TIMES_INFO,
+ (task_info_t) &thread_info_data,
+ &thread_info_count);
+ return (to_double(thread_info_data.user_time) +
+ to_double(thread_info_data.system_time));
+}
diff --git a/cbits/time-posix.c b/cbits/time-posix.c
new file mode 100644
--- /dev/null
+++ b/cbits/time-posix.c
@@ -0,0 +1,24 @@
+#include
+
+void criterion_inittime(void)
+{
+}
+
+double criterion_gettime(void)
+{
+ struct timespec ts;
+
+ clock_gettime(CLOCK_MONOTONIC, &ts);
+
+ return ts.tv_sec + ts.tv_nsec * 1e-9;
+}
+
+
+double criterion_getcputime(void)
+{
+ struct timespec ts;
+
+ clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts);
+
+ return ts.tv_sec + ts.tv_nsec * 1e-9;
+}
diff --git a/cbits/time-windows.c b/cbits/time-windows.c
new file mode 100644
--- /dev/null
+++ b/cbits/time-windows.c
@@ -0,0 +1,80 @@
+/*
+ * Windows has the most amazingly cretinous time measurement APIs you
+ * can possibly imagine.
+ *
+ * Our first possibility is GetSystemTimeAsFileTime, which updates at
+ * roughly 60Hz, and is hence worthless - we'd have to run a
+ * computation for tens or hundreds of seconds to get a trustworthy
+ * number.
+ *
+ * Alternatively, we can use QueryPerformanceCounter, which has
+ * undefined behaviour under almost all interesting circumstances
+ * (e.g. multicore systems, CPU frequency changes). But at least it
+ * increments reasonably often.
+ */
+
+#include
+
+#if 0
+
+void criterion_inittime(void)
+{
+}
+
+double criterion_gettime(void)
+{
+ FILETIME ft;
+ ULARGE_INTEGER li;
+
+ GetSystemTimeAsFileTime(&ft);
+ li.LowPart = ft.dwLowDateTime;
+ li.HighPart = ft.dwHighDateTime;
+
+ return (li.QuadPart - 130000000000000000ull) * 1e-7;
+}
+
+#else
+
+static double freq_recip;
+static LARGE_INTEGER firstClock;
+
+void criterion_inittime(void)
+{
+ LARGE_INTEGER freq;
+
+ if (freq_recip == 0) {
+ QueryPerformanceFrequency(&freq);
+ QueryPerformanceCounter(&firstClock);
+ freq_recip = 1.0 / freq.QuadPart;
+ }
+}
+
+double criterion_gettime(void)
+{
+ LARGE_INTEGER li;
+
+ QueryPerformanceCounter(&li);
+
+ return ((double) (li.QuadPart - firstClock.QuadPart)) * freq_recip;
+}
+
+#endif
+
+static ULONGLONG to_quad_100ns(FILETIME ft)
+{
+ ULARGE_INTEGER li;
+ li.LowPart = ft.dwLowDateTime;
+ li.HighPart = ft.dwHighDateTime;
+ return li.QuadPart;
+}
+
+double criterion_getcputime(void)
+{
+ FILETIME creation, exit, kernel, user;
+ ULONGLONG time;
+
+ GetProcessTimes(GetCurrentProcess(), &creation, &exit, &kernel, &user);
+
+ time = to_quad_100ns(user) + to_quad_100ns(kernel);
+ return time / 1e7;
+}
diff --git a/changelog.md b/changelog.md
new file mode 100644
--- /dev/null
+++ b/changelog.md
@@ -0,0 +1,4 @@
+# 0.1.0
+
+* remove bunch of dependencies
+* initial import of criterion-1.2.2.0
diff --git a/gauge.cabal b/gauge.cabal
new file mode 100644
--- /dev/null
+++ b/gauge.cabal
@@ -0,0 +1,154 @@
+name: gauge
+version: 0.1.0
+synopsis: small framework for performance measurement and analysis
+license: BSD3
+license-file: LICENSE
+author: Bryan O'Sullivan
+maintainer: Vincent Hanquez
+copyright: 2009-2016 Bryan O'Sullivan and others
+category: Development, Performance, Testing, Benchmarking
+homepage: https://github.com/vincenthz/hs-gauge
+bug-reports: https://github.com/vincenthz/hs-gauge/issues
+build-type: Simple
+cabal-version: >= 1.10
+extra-source-files:
+ README.markdown
+ changelog.md
+tested-with:
+ GHC==7.8.4,
+ GHC==7.10.3,
+ GHC==8.0.2,
+ GHC==8.2.1
+
+description:
+ This library provides a powerful but simple way to measure software
+ performance. It provides both a framework for executing and
+ analysing benchmarks and a set of driver functions that makes it
+ easy to build and run benchmarks, and to analyse their results.
+
+library
+ exposed-modules:
+ Gauge
+ Gauge.Main
+ Gauge.Types
+ Gauge.Analysis
+ other-modules:
+ Gauge.IO.Printf
+ Gauge.Internal
+ Gauge.Monad.Internal
+ Gauge.Monad.ExceptT
+ Gauge.Main.Options
+ Gauge.Measurement
+ Gauge.Monad
+ Statistics.Distribution
+ Statistics.Distribution.Normal
+ Statistics.Function
+ Statistics.Internal
+ Statistics.Math.RootFinding
+ Statistics.Matrix
+ Statistics.Matrix.Algorithms
+ Statistics.Matrix.Mutable
+ Statistics.Matrix.Types
+ Statistics.Quantile
+ Statistics.Regression
+ Statistics.Resampling
+ Statistics.Resampling.Bootstrap
+ Statistics.Sample
+ Statistics.Sample.Histogram
+ Statistics.Sample.Internal
+ Statistics.Sample.KernelDensity
+ Statistics.Transform
+ Statistics.Types
+ Statistics.Types.Internal
+
+ hs-source-dirs: . statistics
+
+ c-sources: cbits/cycles.c
+ if os(darwin)
+ c-sources: cbits/time-osx.c
+ else {
+ if os(windows)
+ c-sources: cbits/time-windows.c
+ else
+ c-sources: cbits/time-posix.c
+ }
+
+ other-modules:
+ Paths_gauge
+
+ build-depends:
+ ansi-wl-pprint >= 0.6.7.2,
+ base >= 4.5 && < 5,
+ basement,
+ foundation,
+ code-page,
+ containers,
+ deepseq >= 1.1.0.0,
+ mwc-random >= 0.8.0.3,
+ optparse-applicative >= 0.13,
+ vector >= 0.7.1,
+
+ -- formely statistics dependency that we need
+ math-functions
+
+ default-language: Haskell2010
+ ghc-options: -O2 -Wall -funbox-strict-fields
+
+test-suite sanity
+ type: exitcode-stdio-1.0
+ hs-source-dirs: tests
+ main-is: Sanity.hs
+ default-language: Haskell2010
+ ghc-options: -O2 -Wall -rtsopts
+
+ build-depends:
+ HUnit,
+ base,
+ bytestring,
+ gauge,
+ deepseq,
+ tasty,
+ tasty-hunit
+
+test-suite tests
+ type: exitcode-stdio-1.0
+ hs-source-dirs: tests
+ main-is: Tests.hs
+ default-language: Haskell2010
+ other-modules: Properties
+
+ ghc-options:
+ -Wall -threaded -O0 -rtsopts
+
+ build-depends:
+ QuickCheck >= 2.4,
+ base,
+ gauge,
+ statistics,
+ tasty,
+ tasty-quickcheck,
+ vector
+
+test-suite cleanup
+ type: exitcode-stdio-1.0
+ hs-source-dirs: tests
+ default-language: Haskell2010
+ main-is: Cleanup.hs
+
+ ghc-options:
+ -Wall -threaded -O0 -rtsopts
+
+ build-depends:
+ HUnit,
+ base,
+ bytestring,
+ gauge,
+ deepseq,
+ directory,
+ foundation,
+ tasty,
+ tasty-hunit
+
+source-repository head
+ type: git
+ location: https://github.com/vincenthz/hs-gauge
diff --git a/statistics/Statistics/Distribution.hs b/statistics/Statistics/Distribution.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Distribution.hs
@@ -0,0 +1,70 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}
+-- |
+-- Module : Statistics.Distribution
+-- Copyright : (c) 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Type classes for probability distributions
+
+module Statistics.Distribution
+ (
+ -- * Type classes
+ Distribution(..)
+ , ContDistr(..)
+ ) where
+
+import Prelude hiding (sum)
+
+-- | Type class common to all distributions. Only c.d.f. could be
+-- defined for both discrete and continuous distributions.
+class Distribution d where
+ -- | Cumulative distribution function. The probability that a
+ -- random variable /X/ is less or equal than /x/,
+ -- i.e. P(/X/≤/x/). Cumulative should be defined for
+ -- infinities as well:
+ --
+ -- > cumulative d +∞ = 1
+ -- > cumulative d -∞ = 0
+ cumulative :: d -> Double -> Double
+
+ -- | One's complement of cumulative distibution:
+ --
+ -- > complCumulative d x = 1 - cumulative d x
+ --
+ -- It's useful when one is interested in P(/X/>/x/) and
+ -- expression on the right side begin to lose precision. This
+ -- function have default implementation but implementors are
+ -- encouraged to provide more precise implementation.
+ complCumulative :: d -> Double -> Double
+ complCumulative d x = 1 - cumulative d x
+
+-- | Continuous probability distributuion.
+--
+-- Minimal complete definition is 'quantile' and either 'density' or
+-- 'logDensity'.
+class Distribution d => ContDistr d where
+ -- | Probability density function. Probability that random
+ -- variable /X/ lies in the infinitesimal interval
+ -- [/x/,/x+/δ/x/) equal to /density(x)/⋅δ/x/
+ density :: d -> Double -> Double
+ density d = exp . logDensity d
+
+ -- | Inverse of the cumulative distribution function. The value
+ -- /x/ for which P(/X/≤/x/) = /p/. If probability is outside
+ -- of [0,1] range function should call 'error'
+ quantile :: d -> Double -> Double
+
+ -- | 1-complement of @quantile@:
+ --
+ -- > complQuantile x ≡ quantile (1 - x)
+ complQuantile :: d -> Double -> Double
+ complQuantile d x = quantile d (1 - x)
+
+ -- | Natural logarithm of density.
+ logDensity :: d -> Double -> Double
+ logDensity d = log . density d
diff --git a/statistics/Statistics/Distribution/Normal.hs b/statistics/Statistics/Distribution/Normal.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Distribution/Normal.hs
@@ -0,0 +1,110 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
+-- |
+-- Module : Statistics.Distribution.Normal
+-- Copyright : (c) 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- The normal distribution. This is a continuous probability
+-- distribution that describes data that cluster around a mean.
+
+module Statistics.Distribution.Normal
+ (
+ NormalDistribution
+ -- * Constructors
+ -- , normalDistr
+ --, normalDistrE
+ , standard
+ ) where
+
+import Data.Data (Data, Typeable)
+import GHC.Generics (Generic)
+import Numeric.MathFunctions.Constants (m_sqrt_2, m_sqrt_2_pi)
+import Numeric.SpecFunctions (erfc, invErfc)
+
+import qualified Statistics.Distribution as D
+import Statistics.Internal
+
+
+-- | The normal distribution.
+data NormalDistribution = ND {
+ mean :: {-# UNPACK #-} !Double
+ , stdDev :: {-# UNPACK #-} !Double
+ , ndPdfDenom :: {-# UNPACK #-} !Double
+ , ndCdfDenom :: {-# UNPACK #-} !Double
+ } deriving (Eq, Typeable, Data, Generic)
+
+instance Show NormalDistribution where
+ showsPrec i (ND m s _ _) = defaultShow2 "normalDistr" m s i
+instance Read NormalDistribution where
+ readPrec = defaultReadPrecM2 "normalDistr" normalDistrE
+
+instance D.Distribution NormalDistribution where
+ cumulative = cumulative
+ complCumulative = complCumulative
+
+instance D.ContDistr NormalDistribution where
+ logDensity = logDensity
+ quantile = quantile
+ complQuantile = complQuantile
+
+-- | Standard normal distribution with mean equal to 0 and variance equal to 1
+standard :: NormalDistribution
+standard = ND { mean = 0.0
+ , stdDev = 1.0
+ , ndPdfDenom = log m_sqrt_2_pi
+ , ndCdfDenom = m_sqrt_2
+ }
+
+-- | Create normal distribution from parameters.
+--
+-- IMPORTANT: prior to 0.10 release second parameter was variance not
+-- standard deviation.
+normalDistrE :: Double -- ^ Mean of distribution
+ -> Double -- ^ Standard deviation of distribution
+ -> Maybe NormalDistribution
+normalDistrE m sd
+ | sd > 0 = Just ND { mean = m
+ , stdDev = sd
+ , ndPdfDenom = log $ m_sqrt_2_pi * sd
+ , ndCdfDenom = m_sqrt_2 * sd
+ }
+ | otherwise = Nothing
+
+logDensity :: NormalDistribution -> Double -> Double
+logDensity d x = (-xm * xm / (2 * sd * sd)) - ndPdfDenom d
+ where xm = x - mean d
+ sd = stdDev d
+
+cumulative :: NormalDistribution -> Double -> Double
+cumulative d x = erfc ((mean d - x) / ndCdfDenom d) / 2
+
+complCumulative :: NormalDistribution -> Double -> Double
+complCumulative d x = erfc ((x - mean d) / ndCdfDenom d) / 2
+
+quantile :: NormalDistribution -> Double -> Double
+quantile d p
+ | p == 0 = -inf
+ | p == 1 = inf
+ | p == 0.5 = mean d
+ | p > 0 && p < 1 = x * ndCdfDenom d + mean d
+ | otherwise =
+ error $ "Statistics.Distribution.Normal.quantile: p must be in [0,1] range. Got: "++show p
+ where x = - invErfc (2 * p)
+ inf = 1/0
+
+complQuantile :: NormalDistribution -> Double -> Double
+complQuantile d p
+ | p == 0 = inf
+ | p == 1 = -inf
+ | p == 0.5 = mean d
+ | p > 0 && p < 1 = x * ndCdfDenom d + mean d
+ | otherwise =
+ error $ "Statistics.Distribution.Normal.complQuantile: p must be in [0,1] range. Got: "++show p
+ where x = invErfc (2 * p)
+ inf = 1/0
diff --git a/statistics/Statistics/Function.hs b/statistics/Statistics/Function.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Function.hs
@@ -0,0 +1,204 @@
+{-# LANGUAGE BangPatterns, CPP, FlexibleContexts, Rank2Types #-}
+{-# LANGUAGE TypeFamilies #-}
+#if __GLASGOW_HASKELL__ >= 704
+{-# OPTIONS_GHC -fsimpl-tick-factor=200 #-}
+#endif
+
+-- |
+-- Module : Statistics.Function
+-- Copyright : (c) 2009, 2010, 2011 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Useful functions.
+
+module Statistics.Function
+ (
+ -- * Scanning
+ minMax
+ -- * Sorting
+ , sort
+ , inplaceSortIO
+ -- * Indexing
+ , indices
+ -- * Bit twiddling
+ , nextHighestPowerOfTwo
+ -- * Comparison
+ , within
+ -- * Arithmetic
+ , square
+ -- * Vectors
+ , unsafeModify
+ -- * Combinators
+ , for
+ , rfor
+ ) where
+
+#include "MachDeps.h"
+
+import Control.Monad.ST (ST)
+import Data.Bits ((.|.), shiftR)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as M
+import Numeric.MathFunctions.Comparison (within)
+import Basement.Monad
+
+-- | Sort a vector.
+sort :: U.Vector Double -> U.Vector Double
+sort = G.modify inplaceSortST
+{-# NOINLINE sort #-}
+
+inplaceSortST :: M.MVector s Double
+ -> ST s ()
+inplaceSortST mvec = qsort 0 (M.length mvec-1)
+ where
+ qsort lo hi
+ | lo >= hi = pure ()
+ | otherwise = do
+ p <- partition lo hi
+ qsort lo (pred p)
+ qsort (p+1) hi
+ pivotStrategy low high = do
+ let mid = (low + high) `div` 2
+ pivot <- M.unsafeRead mvec mid
+ M.unsafeRead mvec high >>= M.unsafeWrite mvec mid
+ M.unsafeWrite mvec high pivot
+ pure pivot
+ partition lo hi = do
+ pivot <- pivotStrategy lo hi
+ let go iOrig jOrig = do
+ let fw k = do ak <- M.unsafeRead mvec k
+ if compare ak pivot == LT
+ then fw (k+1)
+ else pure (k, ak)
+ (i, ai) <- fw iOrig
+ let bw k | k==i = pure (i, ai)
+ | otherwise = do ak <- M.unsafeRead mvec k
+ if compare ak pivot /= LT
+ then bw (pred k)
+ else pure (k, ak)
+ (j, aj) <- bw jOrig
+ if i < j
+ then do
+ M.unsafeWrite mvec i aj
+ M.unsafeWrite mvec j ai
+ go (i+1) (pred j)
+ else do
+ M.unsafeWrite mvec hi ai
+ M.unsafeWrite mvec i pivot
+ pure i
+ go lo hi
+
+inplaceSortIO :: M.MVector (PrimState IO) Double
+ -> IO ()
+inplaceSortIO mvec = qsort 0 (M.length mvec-1)
+ where
+ qsort lo hi
+ | lo >= hi = pure ()
+ | otherwise = do
+ p <- partition lo hi
+ qsort lo (pred p)
+ qsort (p+1) hi
+ pivotStrategy low high = do
+ let mid = (low + high) `div` 2
+ pivot <- M.unsafeRead mvec mid
+ M.unsafeRead mvec high >>= M.unsafeWrite mvec mid
+ M.unsafeWrite mvec high pivot
+ pure pivot
+ partition lo hi = do
+ pivot <- pivotStrategy lo hi
+ let go iOrig jOrig = do
+ let fw k = do ak <- M.unsafeRead mvec k
+ if compare ak pivot == LT
+ then fw (k+1)
+ else pure (k, ak)
+ (i, ai) <- fw iOrig
+ let bw k | k==i = pure (i, ai)
+ | otherwise = do ak <- M.unsafeRead mvec k
+ if compare ak pivot /= LT
+ then bw (pred k)
+ else pure (k, ak)
+ (j, aj) <- bw jOrig
+ if i < j
+ then do
+ M.unsafeWrite mvec i aj
+ M.unsafeWrite mvec j ai
+ go (i+1) (pred j)
+ else do
+ M.unsafeWrite mvec hi ai
+ M.unsafeWrite mvec i pivot
+ pure i
+ go lo hi
+
+-- | Return the indices of a vector.
+indices :: (G.Vector v a, G.Vector v Int) => v a -> v Int
+indices a = G.enumFromTo 0 (G.length a - 1)
+{-# INLINE indices #-}
+
+data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double
+
+-- | Compute the minimum and maximum of a vector in one pass.
+minMax :: (G.Vector v Double) => v Double -> (Double, Double)
+minMax = fini . G.foldl' go (MM (1/0) (-1/0))
+ where
+ go (MM lo hi) k = MM (min lo k) (max hi k)
+ fini (MM lo hi) = (lo, hi)
+{-# INLINE minMax #-}
+
+-- | Efficiently compute the next highest power of two for a
+-- non-negative integer. If the given value is already a power of
+-- two, it is returned unchanged. If negative, zero is returned.
+nextHighestPowerOfTwo :: Int -> Int
+nextHighestPowerOfTwo n
+#if WORD_SIZE_IN_BITS == 64
+ = 1 + _i32
+#else
+ = 1 + i16
+#endif
+ where
+ i0 = n - 1
+ i1 = i0 .|. i0 `shiftR` 1
+ i2 = i1 .|. i1 `shiftR` 2
+ i4 = i2 .|. i2 `shiftR` 4
+ i8 = i4 .|. i4 `shiftR` 8
+ i16 = i8 .|. i8 `shiftR` 16
+ _i32 = i16 .|. i16 `shiftR` 32
+-- It could be implemented as
+--
+-- > nextHighestPowerOfTwo n = 1 + foldl' go (n-1) [1, 2, 4, 8, 16, 32]
+-- where go m i = m .|. m `shiftR` i
+--
+-- But GHC do not inline foldl (probably because it's recursive) and
+-- as result function walks list of boxed ints. Hand rolled version
+-- uses unboxed arithmetic.
+
+-- | Multiply a number by itself.
+square :: Double -> Double
+square x = x * x
+
+-- | Simple for loop. Counts from /start/ to /end/-1.
+for :: Monad m => Int -> Int -> (Int -> m ()) -> m ()
+for n0 !n f = loop n0
+ where
+ loop i | i == n = return ()
+ | otherwise = f i >> loop (i+1)
+{-# INLINE for #-}
+
+-- | Simple reverse-for loop. Counts from /start/-1 to /end/ (which
+-- must be less than /start/).
+rfor :: Monad m => Int -> Int -> (Int -> m ()) -> m ()
+rfor n0 !n f = loop n0
+ where
+ loop i | i == n = return ()
+ | otherwise = let i' = i-1 in f i' >> loop i'
+{-# INLINE rfor #-}
+
+unsafeModify :: M.MVector s Double -> Int -> (Double -> Double) -> ST s ()
+unsafeModify v i f = do
+ k <- M.unsafeRead v i
+ M.unsafeWrite v i (f k)
+{-# INLINE unsafeModify #-}
diff --git a/statistics/Statistics/Internal.hs b/statistics/Statistics/Internal.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Internal.hs
@@ -0,0 +1,71 @@
+-- |
+-- Module : Statistics.Internal
+-- Copyright : (c) 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+--
+module Statistics.Internal (
+ -- * Default definitions for Show
+ defaultShow1
+ , defaultShow2
+ -- * Default definitions for Read
+ , defaultReadPrecM1
+ , defaultReadPrecM2
+ -- * Reexports
+ , Show(..)
+ , Read(..)
+ ) where
+
+import Control.Applicative
+import Control.Monad
+import Text.Read
+
+
+
+----------------------------------------------------------------
+-- Default show implementations
+----------------------------------------------------------------
+
+defaultShow1 :: (Show a) => String -> a -> Int -> ShowS
+defaultShow1 con a n
+ = showParen (n >= 11)
+ ( showString con
+ . showChar ' '
+ . showsPrec 11 a
+ )
+
+defaultShow2 :: (Show a, Show b) => String -> a -> b -> Int -> ShowS
+defaultShow2 con a b n
+ = showParen (n >= 11)
+ ( showString con
+ . showChar ' '
+ . showsPrec 11 a
+ . showChar ' '
+ . showsPrec 11 b
+ )
+
+----------------------------------------------------------------
+-- Default read implementations
+----------------------------------------------------------------
+
+defaultReadPrecM1 :: (Read a) => String -> (a -> Maybe r) -> ReadPrec r
+defaultReadPrecM1 con f = parens $ prec 10 $ do
+ expect con
+ a <- readPrec
+ maybe empty return $ f a
+
+defaultReadPrecM2 :: (Read a, Read b) => String -> (a -> b -> Maybe r) -> ReadPrec r
+defaultReadPrecM2 con f = parens $ prec 10 $ do
+ expect con
+ a <- readPrec
+ b <- readPrec
+ maybe empty return $ f a b
+
+expect :: String -> ReadPrec ()
+expect str = do
+ Ident s <- lexP
+ guard (s == str)
diff --git a/statistics/Statistics/Math/RootFinding.hs b/statistics/Statistics/Math/RootFinding.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Math/RootFinding.hs
@@ -0,0 +1,127 @@
+{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric #-}
+
+-- |
+-- Module : Statistics.Math.RootFinding
+-- Copyright : (c) 2011 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Haskell functions for finding the roots of mathematical functions.
+
+module Statistics.Math.RootFinding
+ (
+ Root(..)
+ , fromRoot
+ , ridders
+ -- * References
+ -- $references
+ ) where
+
+import Control.Applicative (Alternative(..), Applicative(..))
+import Control.Monad (MonadPlus(..), ap)
+import Data.Data (Data, Typeable)
+import GHC.Generics (Generic)
+import Numeric.MathFunctions.Comparison (within)
+
+
+-- | The result of searching for a root of a mathematical function.
+data Root a = NotBracketed
+ -- ^ The function does not have opposite signs when
+ -- evaluated at the lower and upper bounds of the search.
+ | SearchFailed
+ -- ^ The search failed to converge to within the given
+ -- error tolerance after the given number of iterations.
+ | Root a
+ -- ^ A root was successfully found.
+ deriving (Eq, Read, Show, Typeable, Data, Generic)
+
+instance Functor Root where
+ fmap _ NotBracketed = NotBracketed
+ fmap _ SearchFailed = SearchFailed
+ fmap f (Root a) = Root (f a)
+
+instance Monad Root where
+ NotBracketed >>= _ = NotBracketed
+ SearchFailed >>= _ = SearchFailed
+ Root a >>= m = m a
+
+ return = Root
+
+instance MonadPlus Root where
+ mzero = SearchFailed
+
+ r@(Root _) `mplus` _ = r
+ _ `mplus` p = p
+
+instance Applicative Root where
+ pure = Root
+ (<*>) = ap
+
+instance Alternative Root where
+ empty = SearchFailed
+
+ r@(Root _) <|> _ = r
+ _ <|> p = p
+
+-- | Returns either the result of a search for a root, or the default
+-- value if the search failed.
+fromRoot :: a -- ^ Default value.
+ -> Root a -- ^ Result of search for a root.
+ -> a
+fromRoot _ (Root a) = a
+fromRoot a _ = a
+
+
+-- | Use the method of Ridders to compute a root of a function.
+--
+-- The function must have opposite signs when evaluated at the lower
+-- and upper bounds of the search (i.e. the root must be bracketed).
+ridders :: Double -- ^ Absolute error tolerance.
+ -> (Double,Double) -- ^ Lower and upper bounds for the search.
+ -> (Double -> Double) -- ^ Function to find the roots of.
+ -> Root Double
+ridders tol (lo,hi) f
+ | flo == 0 = Root lo
+ | fhi == 0 = Root hi
+ | flo*fhi > 0 = NotBracketed -- root is not bracketed
+ | otherwise = go lo flo hi fhi 0
+ where
+ go !a !fa !b !fb !i
+ -- Root is bracketed within 1 ulp. No improvement could be made
+ | within 1 a b = Root a
+ -- Root is found. Check that f(m) == 0 is nessesary to ensure
+ -- that root is never passed to 'go'
+ | fm == 0 = Root m
+ | fn == 0 = Root n
+ | d < tol = Root n
+ -- Too many iterations performed. Fail
+ | i >= (100 :: Int) = SearchFailed
+ -- Ridder's approximation coincide with one of old
+ -- bounds. Revert to bisection
+ | n == a || n == b = case () of
+ _| fm*fa < 0 -> go a fa m fm (i+1)
+ | otherwise -> go m fm b fb (i+1)
+ -- Proceed as usual
+ | fn*fm < 0 = go n fn m fm (i+1)
+ | fn*fa < 0 = go a fa n fn (i+1)
+ | otherwise = go n fn b fb (i+1)
+ where
+ d = abs (b - a)
+ dm = (b - a) * 0.5
+ !m = a + dm
+ !fm = f m
+ !dn = signum (fb - fa) * dm * fm / sqrt(fm*fm - fa*fb)
+ !n = m - signum dn * min (abs dn) (abs dm - 0.5 * tol)
+ !fn = f n
+ !flo = f lo
+ !fhi = f hi
+
+
+-- $references
+--
+-- * Ridders, C.F.J. (1979) A new algorithm for computing a single
+-- root of a real continuous function.
+-- /IEEE Transactions on Circuits and Systems/ 26:979–980.
diff --git a/statistics/Statistics/Matrix.hs b/statistics/Statistics/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Matrix.hs
@@ -0,0 +1,95 @@
+{-# LANGUAGE PatternGuards #-}
+-- |
+-- Module : Statistics.Matrix
+-- Copyright : 2011 Aleksey Khudyakov, 2014 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Basic matrix operations.
+--
+-- There isn't a widely used matrix package for Haskell yet, so
+-- we implement the necessary minimum here.
+
+module Statistics.Matrix
+ ( -- * Data types
+ Matrix(..)
+ , Vector
+ -- * Conversion from/to lists/vectors
+ , fromVector
+ , dimension
+ -- , center
+ , multiplyV
+ , transpose
+ , norm
+ , column
+ -- , row
+ , for
+ , unsafeIndex
+ ) where
+
+import Prelude hiding (exponent, map, sum)
+import qualified Data.Vector.Unboxed as U
+
+import Statistics.Function (for, square)
+import Statistics.Matrix.Types
+import Statistics.Sample.Internal (sum)
+
+
+----------------------------------------------------------------
+-- Conversion to/from vectors/lists
+----------------------------------------------------------------
+
+-- | Convert from a row-major vector.
+fromVector :: Int -- ^ Number of rows.
+ -> Int -- ^ Number of columns.
+ -> U.Vector Double -- ^ Flat list of values, in row-major order.
+ -> Matrix
+fromVector r c v
+ | r*c /= len = error "input size mismatch"
+ | otherwise = Matrix r c 0 v
+ where len = U.length v
+
+----------------------------------------------------------------
+-- Other
+----------------------------------------------------------------
+
+-- | Return the dimensions of this matrix, as a (row,column) pair.
+dimension :: Matrix -> (Int, Int)
+dimension (Matrix r c _ _) = (r, c)
+
+-- | Matrix-vector multiplication.
+multiplyV :: Matrix -> Vector -> Vector
+multiplyV m v
+ | cols m == c = U.generate (rows m) (sum . U.zipWith (*) v . row m)
+ | otherwise = error $ "matrix/vector unconformable " ++ show (cols m,c)
+ where c = U.length v
+
+-- | Calculate the Euclidean norm of a vector.
+norm :: Vector -> Double
+norm = sqrt . sum . U.map square
+
+-- | Return the given column.
+column :: Matrix -> Int -> Vector
+column (Matrix r c _ v) i = U.backpermute v $ U.enumFromStepN i c r
+{-# INLINE column #-}
+
+-- | Return the given row.
+row :: Matrix -> Int -> Vector
+row (Matrix _ c _ v) i = U.slice (c*i) c v
+
+unsafeIndex :: Matrix
+ -> Int -- ^ Row.
+ -> Int -- ^ Column.
+ -> Double
+unsafeIndex = unsafeBounds U.unsafeIndex
+
+-- | Given row and column numbers, calculate the offset into the flat
+-- row-major vector, without checking.
+unsafeBounds :: (Vector -> Int -> r) -> Matrix -> Int -> Int -> r
+unsafeBounds k (Matrix _ cs _ v) r c = k v $! r * cs + c
+{-# INLINE unsafeBounds #-}
+
+
+transpose :: Matrix -> Matrix
+transpose m@(Matrix r0 c0 e _) = Matrix c0 r0 e . U.generate (r0*c0) $ \i ->
+ let (r,c) = i `quotRem` r0
+ in unsafeIndex m c r
diff --git a/statistics/Statistics/Matrix/Algorithms.hs b/statistics/Statistics/Matrix/Algorithms.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Matrix/Algorithms.hs
@@ -0,0 +1,42 @@
+-- |
+-- Module : Statistics.Matrix.Algorithms
+-- Copyright : 2014 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Useful matrix functions.
+
+module Statistics.Matrix.Algorithms
+ (
+ qr
+ ) where
+
+import Control.Applicative ((<$>), (<*>))
+import Control.Monad.ST (ST, runST)
+import Prelude hiding (sum, replicate)
+import Statistics.Matrix (Matrix, column, dimension, for, norm)
+import qualified Statistics.Matrix.Mutable as M
+import Statistics.Sample.Internal (sum)
+import qualified Data.Vector.Unboxed as U
+
+-- | /O(r*c)/ Compute the QR decomposition of a matrix.
+-- The result returned is the matrices (/q/,/r/).
+qr :: Matrix -> (Matrix, Matrix)
+qr mat = runST $ do
+ let (m,n) = dimension mat
+ r <- M.replicate n n 0
+ a <- M.thaw mat
+ for 0 n $ \j -> do
+ cn <- M.immutably a $ \aa -> norm (column aa j)
+ M.unsafeWrite r j j cn
+ for 0 m $ \i -> M.unsafeModify a i j (/ cn)
+ for (j+1) n $ \jj -> do
+ p <- innerProduct a j jj
+ M.unsafeWrite r j jj p
+ for 0 m $ \i -> do
+ aij <- M.unsafeRead a i j
+ M.unsafeModify a i jj $ subtract (p * aij)
+ (,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r
+
+innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double
+innerProduct mmat j k = M.immutably mmat $ \mat ->
+ sum $ U.zipWith (*) (column mat j) (column mat k)
diff --git a/statistics/Statistics/Matrix/Mutable.hs b/statistics/Statistics/Matrix/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Matrix/Mutable.hs
@@ -0,0 +1,63 @@
+-- |
+-- Module : Statistics.Matrix.Mutable
+-- Copyright : (c) 2014 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Basic mutable matrix operations.
+
+module Statistics.Matrix.Mutable
+ (
+ MMatrix(..)
+ , MVector
+ , replicate
+ , thaw
+ , unsafeFreeze
+ , unsafeRead
+ , unsafeWrite
+ , unsafeModify
+ , immutably
+ ) where
+
+import Control.Applicative ((<$>))
+import Control.DeepSeq (NFData(..))
+import Control.Monad.ST (ST)
+import Statistics.Matrix.Types (Matrix(..), MMatrix(..), MVector)
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as M
+import Prelude hiding (replicate)
+
+replicate :: Int -> Int -> Double -> ST s (MMatrix s)
+replicate r c k = MMatrix r c 0 <$> M.replicate (r*c) k
+
+thaw :: Matrix -> ST s (MMatrix s)
+thaw (Matrix r c e v) = MMatrix r c e <$> U.thaw v
+
+unsafeFreeze :: MMatrix s -> ST s Matrix
+unsafeFreeze (MMatrix r c e mv) = Matrix r c e <$> U.unsafeFreeze mv
+
+unsafeRead :: MMatrix s -> Int -> Int -> ST s Double
+unsafeRead mat r c = unsafeBounds mat r c M.unsafeRead
+{-# INLINE unsafeRead #-}
+
+unsafeWrite :: MMatrix s -> Int -> Int -> Double -> ST s ()
+unsafeWrite mat row col k = unsafeBounds mat row col $ \v i ->
+ M.unsafeWrite v i k
+{-# INLINE unsafeWrite #-}
+
+unsafeModify :: MMatrix s -> Int -> Int -> (Double -> Double) -> ST s ()
+unsafeModify mat row col f = unsafeBounds mat row col $ \v i -> do
+ k <- M.unsafeRead v i
+ M.unsafeWrite v i (f k)
+{-# INLINE unsafeModify #-}
+
+-- | Given row and column numbers, calculate the offset into the flat
+-- row-major vector, without checking.
+unsafeBounds :: MMatrix s -> Int -> Int -> (MVector s -> Int -> r) -> r
+unsafeBounds (MMatrix _ cs _ mv) r c k = k mv $! r * cs + c
+{-# INLINE unsafeBounds #-}
+
+immutably :: NFData a => MMatrix s -> (Matrix -> a) -> ST s a
+immutably mmat f = do
+ k <- f <$> unsafeFreeze mmat
+ rnf k `seq` return k
+{-# INLINE immutably #-}
diff --git a/statistics/Statistics/Matrix/Types.hs b/statistics/Statistics/Matrix/Types.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Matrix/Types.hs
@@ -0,0 +1,63 @@
+-- |
+-- Module : Statistics.Matrix.Types
+-- Copyright : 2014 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Basic matrix operations.
+--
+-- There isn't a widely used matrix package for Haskell yet, so
+-- we implement the necessary minimum here.
+
+module Statistics.Matrix.Types
+ (
+ Vector
+ , MVector
+ , Matrix(..)
+ , MMatrix(..)
+ ) where
+
+import Data.Char (isSpace)
+import Numeric (showFFloat)
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as M
+
+type Vector = U.Vector Double
+type MVector s = M.MVector s Double
+
+-- | Two-dimensional matrix, stored in row-major order.
+data Matrix = Matrix {
+ rows :: {-# UNPACK #-} !Int -- ^ Rows of matrix.
+ , cols :: {-# UNPACK #-} !Int -- ^ Columns of matrix.
+ , exponent :: {-# UNPACK #-} !Int
+ -- ^ In order to avoid overflows during matrix multiplication, a
+ -- large exponent is stored separately.
+ , _vector :: !Vector -- ^ Matrix data.
+ } deriving (Eq)
+
+-- | Two-dimensional mutable matrix, stored in row-major order.
+data MMatrix s = MMatrix
+ {-# UNPACK #-} !Int
+ {-# UNPACK #-} !Int
+ {-# UNPACK #-} !Int
+ !(MVector s)
+
+-- The Show instance is useful only for debugging.
+instance Show Matrix where
+ show = debug
+
+debug :: Matrix -> String
+debug (Matrix r c _ vs) = unlines $ zipWith (++) (hdr0 : repeat hdr) rrows
+ where
+ rrows = map (cleanEnd . unwords) . split $ zipWith (++) ldone tdone
+ hdr0 = show (r,c) ++ " "
+ hdr = replicate (length hdr0) ' '
+ pad plus k xs = replicate (k - length xs) ' ' `plus` xs
+ ldone = map (pad (++) (longest lstr)) lstr
+ tdone = map (pad (flip (++)) (longest tstr)) tstr
+ (lstr, tstr) = unzip . map (break (=='.') . render) . U.toList $ vs
+ longest = maximum . map length
+ render k = reverse . dropWhile (=='.') . dropWhile (=='0') . reverse .
+ showFFloat (Just 4) k $ ""
+ split [] = []
+ split xs = i : split rest where (i, rest) = splitAt c xs
+ cleanEnd = reverse . dropWhile isSpace . reverse
diff --git a/statistics/Statistics/Quantile.hs b/statistics/Statistics/Quantile.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Quantile.hs
@@ -0,0 +1,81 @@
+{-# LANGUAGE FlexibleContexts #-}
+-- |
+-- Module : Statistics.Quantile
+-- Copyright : (c) 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Functions for approximating quantiles, i.e. points taken at regular
+-- intervals from the cumulative distribution function of a random
+-- variable.
+--
+-- The number of quantiles is described below by the variable /q/, so
+-- with /q/=4, a 4-quantile (also known as a /quartile/) has 4
+-- intervals, and contains 5 points. The parameter /k/ describes the
+-- desired point, where 0 ≤ /k/ ≤ /q/.
+
+module Statistics.Quantile
+ (
+
+ -- * Quantile estimation functions
+ weightedAvg
+ , Sorted(..)
+ -- * References
+ -- $references
+ ) where
+
+import Data.Vector.Generic ((!))
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+newtype Sorted x = Sorted x
+
+-- | O(/n/ log /n/). Estimate the /k/th /q/-quantile of a sample,
+-- using the weighted average method.
+--
+-- The following properties should hold:
+-- * the length of the input is greater than @0@
+-- * the input does not contain @NaN@
+-- * k ≥ 0 and k ≤ q
+--
+-- otherwise an error will be thrown.
+weightedAvg :: G.Vector v Double =>
+ Int -- ^ /k/, the desired quantile.
+ -> Int -- ^ /q/, the number of quantiles.
+ -> Sorted (v Double) -- ^ /x/, the sample data.
+ -> Double
+weightedAvg k q (Sorted x)
+ | G.any isNaN x = modErr "weightedAvg" "Sample contains NaNs"
+ | n == 0 = modErr "weightedAvg" "Sample is empty"
+ | n == 1 = G.head x
+ | q < 2 = modErr "weightedAvg" "At least 2 quantiles is needed"
+ | k == q = G.maximum x
+ | k >= 0 || k < q = xj + g * (xj1 - xj)
+ | otherwise = modErr "weightedAvg" "Wrong quantile number"
+ where
+ j = floor idx
+ idx = fromIntegral (n - 1) * fromIntegral k / fromIntegral q
+ g = idx - fromIntegral j
+ xj = x ! j
+ xj1 = x ! (j+1)
+ n = G.length x
+{-# SPECIALIZE weightedAvg :: Int -> Int -> Sorted (U.Vector Double) -> Double #-}
+{-# SPECIALIZE weightedAvg :: Int -> Int -> Sorted (V.Vector Double) -> Double #-}
+
+modErr :: String -> String -> a
+modErr f err = error $ "Statistics.Quantile." ++ f ++ ": " ++ err
+
+
+
+-- $references
+--
+-- * Weisstein, E.W. Quantile. /MathWorld/.
+--
+--
+-- * Hyndman, R.J.; Fan, Y. (1996) Sample quantiles in statistical
+-- packages. /American Statistician/
+-- 50(4):361–365.
diff --git a/statistics/Statistics/Regression.hs b/statistics/Statistics/Regression.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Regression.hs
@@ -0,0 +1,152 @@
+-- |
+-- Module : Statistics.Regression
+-- Copyright : 2014 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Functions for regression analysis.
+
+module Statistics.Regression
+ (
+ olsRegress
+ , bootstrapRegress
+ ) where
+
+import Control.Applicative ((<$>))
+import Control.Concurrent (forkIO)
+import Control.Concurrent.Chan (newChan, readChan, writeChan)
+import Control.DeepSeq (rnf)
+import Control.Monad (forM_, replicateM)
+import GHC.Conc (getNumCapabilities)
+import Prelude hiding (pred, sum)
+import Statistics.Function as F
+import Statistics.Matrix
+import Statistics.Matrix.Algorithms (qr)
+import Statistics.Resampling (splitGen)
+import Statistics.Types (Estimate(..),ConfInt,CL,estimateFromInterval,significanceLevel)
+import Statistics.Sample (mean)
+import Statistics.Sample.Internal (sum)
+import System.Random.MWC (GenIO, uniformR)
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as M
+
+-- | Perform an ordinary least-squares regression on a set of
+-- predictors, and calculate the goodness-of-fit of the regression.
+--
+-- The returned pair consists of:
+--
+-- * A vector of regression coefficients. This vector has /one more/
+-- element than the list of predictors; the last element is the
+-- /y/-intercept value.
+--
+-- * /R²/, the coefficient of determination (see 'rSquare' for
+-- details).
+olsRegress :: [Vector]
+ -- ^ Non-empty list of predictor vectors. Must all have
+ -- the same length. These will become the columns of
+ -- the matrix /A/ solved by 'ols'.
+ -> Vector
+ -- ^ Responder vector. Must have the same length as the
+ -- predictor vectors.
+ -> (Vector, Double)
+olsRegress preds@(_:_) resps
+ | any (/=n) ls = error $ "predictor vector length mismatch " ++
+ show lss
+ | G.length resps /= n = error $ "responder/predictor length mismatch " ++
+ show (G.length resps, n)
+ | otherwise = (coeffs, rSquare mxpreds resps coeffs)
+ where
+ coeffs = ols mxpreds resps
+ mxpreds = transpose .
+ fromVector (length lss + 1) n .
+ G.concat $ preds ++ [G.replicate n 1]
+ lss@(n:ls) = map G.length preds
+olsRegress _ _ = error "no predictors given"
+
+-- | Compute the ordinary least-squares solution to /A x = b/.
+ols :: Matrix -- ^ /A/ has at least as many rows as columns.
+ -> Vector -- ^ /b/ has the same length as columns in /A/.
+ -> Vector
+ols a b
+ | rs < cs = error $ "fewer rows than columns " ++ show d
+ | otherwise = solve r (transpose q `multiplyV` b)
+ where
+ d@(rs,cs) = dimension a
+ (q,r) = qr a
+
+-- | Solve the equation /R x = b/.
+solve :: Matrix -- ^ /R/ is an upper-triangular square matrix.
+ -> Vector -- ^ /b/ is of the same length as rows\/columns in /R/.
+ -> Vector
+solve r b
+ | n /= l = error $ "row/vector mismatch " ++ show (n,l)
+ | otherwise = U.create $ do
+ s <- U.thaw b
+ rfor n 0 $ \i -> do
+ si <- (/ unsafeIndex r i i) <$> M.unsafeRead s i
+ M.unsafeWrite s i si
+ for 0 i $ \j -> F.unsafeModify s j $ subtract (unsafeIndex r j i * si)
+ return s
+ where n = rows r
+ l = U.length b
+
+-- | Compute /R²/, the coefficient of determination that
+-- indicates goodness-of-fit of a regression.
+--
+-- This value will be 1 if the predictors fit perfectly, dropping to 0
+-- if they have no explanatory power.
+rSquare :: Matrix -- ^ Predictors (regressors).
+ -> Vector -- ^ Responders.
+ -> Vector -- ^ Regression coefficients.
+ -> Double
+rSquare pred resp coeff = 1 - r / t
+ where
+ r = sum $ flip U.imap resp $ \i x -> square (x - p i)
+ t = sum $ flip U.map resp $ \x -> square (x - mean resp)
+ p i = sum . flip U.imap coeff $ \j -> (* unsafeIndex pred i j)
+
+-- | Bootstrap a regression function. Returns both the results of the
+-- regression and the requested confidence interval values.
+bootstrapRegress
+ :: GenIO
+ -> Int -- ^ Number of resamples to compute.
+ -> CL Double -- ^ Confidence level.
+ -> ([Vector] -> Vector -> (Vector, Double))
+ -- ^ Regression function.
+ -> [Vector] -- ^ Predictor vectors.
+ -> Vector -- ^ Responder vector.
+ -> IO (V.Vector (Estimate ConfInt Double), Estimate ConfInt Double)
+bootstrapRegress gen0 numResamples cl rgrss preds0 resp0
+ | numResamples < 1 = error $ "bootstrapRegress: number of resamples " ++
+ "must be positive"
+ | otherwise = do
+ caps <- getNumCapabilities
+ gens <- splitGen caps gen0
+ done <- newChan
+ forM_ (zip gens (balance caps numResamples)) $ \(gen,count) -> forkIO $ do
+ v <- V.replicateM count $ do
+ let n = U.length resp0
+ ixs <- U.replicateM n $ uniformR (0,n-1) gen
+ let resp = U.backpermute resp0 ixs
+ preds = map (flip U.backpermute ixs) preds0
+ return $ rgrss preds resp
+ rnf v `seq` writeChan done v
+ (coeffsv, r2v) <- (G.unzip . V.concat) <$> replicateM caps (readChan done)
+ let coeffs = flip G.imap (G.convert coeffss) $ \i x ->
+ est x . U.generate numResamples $ \k -> (coeffsv G.! k) G.! i
+ r2 = est r2s (G.convert r2v)
+ (coeffss, r2s) = rgrss preds0 resp0
+ est s v = estimateFromInterval s (w G.! lo, w G.! hi) cl
+ where w = F.sort v
+ lo = round c
+ hi = truncate (n - c)
+ n = fromIntegral numResamples
+ c = n * (significanceLevel cl / 2)
+ return (coeffs, r2)
+
+-- | Balance units of work across workers.
+balance :: Int -> Int -> [Int]
+balance numSlices numItems = zipWith (+) (replicate numSlices q)
+ (replicate r 1 ++ repeat 0)
+ where (q,r) = numItems `quotRem` numSlices
diff --git a/statistics/Statistics/Resampling.hs b/statistics/Statistics/Resampling.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Resampling.hs
@@ -0,0 +1,219 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE DeriveFoldable #-}
+{-# LANGUAGE DeriveTraversable #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric, FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module : Statistics.Resampling
+-- Copyright : (c) 2009, 2010 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Resampling statistics.
+
+module Statistics.Resampling
+ ( -- * Data types
+ Bootstrap(..)
+ , Estimator(..)
+ , resample
+ -- * Jackknife
+ , jackknife
+ -- * Helper functions
+ , splitGen
+ ) where
+
+import Control.Concurrent (forkIO, newChan, readChan, writeChan)
+import Control.Monad
+import Basement.Monad (PrimMonad(..))
+import Data.Data (Data, Typeable)
+import Data.Vector.Generic (unsafeFreeze)
+import Data.Word (Word32)
+import qualified Data.Foldable as T
+import qualified Data.Traversable as T
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as MU
+
+import GHC.Conc (numCapabilities)
+import GHC.Generics (Generic)
+import Numeric.Sum (Summation(..), kbn)
+import Statistics.Function (indices, inplaceSortIO)
+import Statistics.Sample (mean, stdDev, variance, varianceUnbiased)
+import Statistics.Types (Sample)
+import System.Random.MWC (Gen, GenIO, initialize, uniformR, uniformVector)
+
+
+----------------------------------------------------------------
+-- Data types
+----------------------------------------------------------------
+
+data Bootstrap v a = Bootstrap
+ { fullSample :: !a
+ , resamples :: v a
+ }
+ deriving (Eq, Read, Show , Generic, Functor, T.Foldable, T.Traversable
+#if __GLASGOW_HASKELL__ >= 708
+ , Typeable, Data
+#endif
+ )
+
+-- | An estimator of a property of a sample, such as its 'mean'.
+--
+-- The use of an algebraic data type here allows functions such as
+-- 'jackknife' and 'bootstrapBCA' to use more efficient algorithms
+-- when possible.
+data Estimator = Mean
+ | Variance
+ | VarianceUnbiased
+ | StdDev
+ | Function (Sample -> Double)
+
+-- | Run an 'Estimator' over a sample.
+estimate :: Estimator -> Sample -> Double
+estimate Mean = mean
+estimate Variance = variance
+estimate VarianceUnbiased = varianceUnbiased
+estimate StdDev = stdDev
+estimate (Function est) = est
+
+
+----------------------------------------------------------------
+-- Resampling
+----------------------------------------------------------------
+
+-- | /O(e*r*s)/ Resample a data set repeatedly, with replacement,
+-- computing each estimate over the resampled data.
+--
+-- This function is expensive; it has to do work proportional to
+-- /e*r*s/, where /e/ is the number of estimation functions, /r/ is
+-- the number of resamples to compute, and /s/ is the number of
+-- original samples.
+--
+-- To improve performance, this function will make use of all
+-- available CPUs. At least with GHC 7.0, parallel performance seems
+-- best if the parallel garbage collector is disabled (RTS option
+-- @-qg@).
+resample :: GenIO
+ -> [Estimator] -- ^ Estimation functions.
+ -> Int -- ^ Number of resamples to compute.
+ -> U.Vector Double -- ^ Original sample.
+ -> IO [(Estimator, Bootstrap U.Vector Double)]
+resample gen ests numResamples samples = do
+ let ixs = scanl (+) 0 $
+ zipWith (+) (replicate numCapabilities q)
+ (replicate r 1 ++ repeat 0)
+ where (q,r) = numResamples `quotRem` numCapabilities
+ results <- mapM (const (MU.new numResamples)) ests
+ done <- newChan
+ gens <- splitGen numCapabilities gen
+ forM_ (zip3 ixs (tail ixs) gens) $ \ (start,!end,gen') ->
+ forkIO $ do
+ let loop k ers | k >= end = writeChan done ()
+ | otherwise = do
+ re <- resampleVector gen' samples
+ forM_ ers $ \(est,arr) ->
+ MU.write arr k . est $ re
+ loop (k+1) ers
+ loop start (zip ests' results)
+ replicateM_ numCapabilities $ readChan done
+ mapM_ inplaceSortIO results
+ -- Build resamples
+ res <- mapM unsafeFreeze results
+ return $ zip ests
+ $ zipWith Bootstrap [estimate e samples | e <- ests]
+ res
+ where
+ ests' = map estimate ests
+
+-- | Create vector using resamples
+resampleVector :: G.Vector v a
+ => Gen (PrimState IO) -> v a -> IO (v a)
+resampleVector gen v
+ = G.replicateM n $ do i <- uniformR (0,n-1) gen
+ return $! G.unsafeIndex v i
+ where
+ n = G.length v
+
+----------------------------------------------------------------
+-- Jackknife
+----------------------------------------------------------------
+
+-- | /O(n) or O(n^2)/ Compute a statistical estimate repeatedly over a
+-- sample, each time omitting a successive element.
+jackknife :: Estimator -> Sample -> U.Vector Double
+jackknife Mean sample = jackknifeMean sample
+jackknife Variance sample = jackknifeVariance sample
+jackknife VarianceUnbiased sample = jackknifeVarianceUnb sample
+jackknife StdDev sample = jackknifeStdDev sample
+jackknife (Function est) sample
+ | G.length sample == 1 = singletonErr "jackknife"
+ | otherwise = U.map f . indices $ sample
+ where f i = est (dropAt i sample)
+
+-- | /O(n)/ Compute the jackknife mean of a sample.
+jackknifeMean :: Sample -> U.Vector Double
+jackknifeMean samp
+ | len == 1 = singletonErr "jackknifeMean"
+ | otherwise = G.map (/l) $ G.zipWith (+) (pfxSumL samp) (pfxSumR samp)
+ where
+ l = fromIntegral (len - 1)
+ len = G.length samp
+
+-- | /O(n)/ Compute the jackknife variance of a sample with a
+-- correction factor @c@, so we can get either the regular or
+-- \"unbiased\" variance.
+jackknifeVariance_ :: Double -> Sample -> U.Vector Double
+jackknifeVariance_ c samp
+ | len == 1 = singletonErr "jackknifeVariance"
+ | otherwise = G.zipWith4 go als ars bls brs
+ where
+ als = pfxSumL . G.map goa $ samp
+ ars = pfxSumR . G.map goa $ samp
+ goa x = v * v where v = x - m
+ bls = pfxSumL . G.map (subtract m) $ samp
+ brs = pfxSumR . G.map (subtract m) $ samp
+ m = mean samp
+ n = fromIntegral len
+ go al ar bl br = (al + ar - (b * b) / q) / (q - c)
+ where b = bl + br
+ q = n - 1
+ len = G.length samp
+
+-- | /O(n)/ Compute the unbiased jackknife variance of a sample.
+jackknifeVarianceUnb :: Sample -> U.Vector Double
+jackknifeVarianceUnb = jackknifeVariance_ 1
+
+-- | /O(n)/ Compute the jackknife variance of a sample.
+jackknifeVariance :: Sample -> U.Vector Double
+jackknifeVariance = jackknifeVariance_ 0
+
+-- | /O(n)/ Compute the jackknife standard deviation of a sample.
+jackknifeStdDev :: Sample -> U.Vector Double
+jackknifeStdDev = G.map sqrt . jackknifeVarianceUnb
+
+pfxSumL :: U.Vector Double -> U.Vector Double
+pfxSumL = G.map kbn . G.scanl add zero
+
+pfxSumR :: U.Vector Double -> U.Vector Double
+pfxSumR = G.tail . G.map kbn . G.scanr (flip add) zero
+
+-- | Drop the /k/th element of a vector.
+dropAt :: U.Unbox e => Int -> U.Vector e -> U.Vector e
+dropAt n v = U.slice 0 n v U.++ U.slice (n+1) (U.length v - n - 1) v
+
+singletonErr :: String -> a
+singletonErr func = error $
+ "Statistics.Resampling." ++ func ++ ": singleton input"
+
+-- | Split a generator into several that can run independently.
+splitGen :: Int -> GenIO -> IO [GenIO]
+splitGen n gen
+ | n <= 0 = return []
+ | otherwise =
+ fmap (gen:) . replicateM (n-1) $
+ initialize =<< (uniformVector gen 256 :: IO (U.Vector Word32))
diff --git a/statistics/Statistics/Resampling/Bootstrap.hs b/statistics/Statistics/Resampling/Bootstrap.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Resampling/Bootstrap.hs
@@ -0,0 +1,80 @@
+-- |
+-- Module : Statistics.Resampling.Bootstrap
+-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- The bootstrap method for statistical inference.
+
+module Statistics.Resampling.Bootstrap
+ ( bootstrapBCA
+ -- * References
+ -- $references
+ ) where
+
+import Data.Vector.Generic ((!))
+import qualified Data.Vector.Unboxed as U
+
+import Statistics.Distribution (cumulative, quantile)
+import Statistics.Distribution.Normal
+import Statistics.Resampling (Bootstrap(..), jackknife)
+import Statistics.Sample (mean)
+import Statistics.Types (Sample, CL, Estimate, ConfInt, estimateFromInterval,
+ estimateFromErr, CL, significanceLevel)
+import qualified Statistics.Resampling as R
+
+
+data T = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double
+infixl 2 :<
+
+-- | Bias-corrected accelerated (BCA) bootstrap. This adjusts for both
+-- bias and skewness in the resampled distribution.
+--
+-- BCA algorithm is described in ch. 5 of Davison, Hinkley "Confidence
+-- intervals" in section 5.3 "Percentile method"
+bootstrapBCA
+ :: CL Double -- ^ Confidence level
+ -> Sample -- ^ Full data sample
+ -> [(R.Estimator, Bootstrap U.Vector Double)]
+ -- ^ Estimates obtained from resampled data and estimator used for
+ -- this.
+ -> [Estimate ConfInt Double]
+bootstrapBCA confidenceLevel sample resampledData
+ = map e resampledData
+ where
+ e (est, Bootstrap pt resample)
+ | U.length sample == 1 || isInfinite bias =
+ estimateFromErr pt (0,0) confidenceLevel
+ | otherwise =
+ estimateFromInterval pt (resample ! lo, resample ! hi) confidenceLevel
+ where
+ -- Quantile estimates for given CL
+ lo = max (cumn a1) 0
+ where a1 = bias + b1 / (1 - accel * b1)
+ b1 = bias + z1
+ hi = min (cumn a2) (ni - 1)
+ where a2 = bias + b2 / (1 - accel * b2)
+ b2 = bias - z1
+ -- Number of resamples
+ ni = U.length resample
+ n = fromIntegral ni
+ -- Corrections
+ z1 = quantile standard (significanceLevel confidenceLevel / 2)
+ cumn = round . (*n) . cumulative standard
+ bias = quantile standard (probN / n)
+ where probN = fromIntegral . U.length . U.filter (
diff --git a/statistics/Statistics/Sample.hs b/statistics/Statistics/Sample.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Sample.hs
@@ -0,0 +1,117 @@
+{-# LANGUAGE FlexibleContexts #-}
+-- |
+-- Module : Statistics.Sample
+-- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Commonly used sample statistics, also known as descriptive
+-- statistics.
+
+module Statistics.Sample
+ (
+ -- * Statistics of location
+ mean
+
+ -- ** Two-pass functions (numerically robust)
+ -- $robust
+ , variance
+ , varianceUnbiased
+ , stdDev
+
+ -- * References
+ -- $references
+ ) where
+
+import Statistics.Sample.Internal (robustSumVar, sum)
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+-- Operator ^ will be overriden
+import Prelude hiding ((^), sum)
+
+-- | /O(n)/ Arithmetic mean. This uses Kahan-Babuška-Neumaier
+-- summation, so is more accurate than 'welfordMean' unless the input
+-- values are very large.
+mean :: (G.Vector v Double) => v Double -> Double
+mean xs = sum xs / fromIntegral (G.length xs)
+{-# SPECIALIZE mean :: U.Vector Double -> Double #-}
+{-# SPECIALIZE mean :: V.Vector Double -> Double #-}
+
+-- $variance
+--
+-- The variance—and hence the standard deviation—of a
+-- sample of fewer than two elements are both defined to be zero.
+
+-- $robust
+--
+-- These functions use the compensated summation algorithm of Chan et
+-- al. for numerical robustness, but require two passes over the
+-- sample data as a result.
+--
+-- Because of the need for two passes, these functions are /not/
+-- subject to stream fusion.
+
+-- | Maximum likelihood estimate of a sample's variance. Also known
+-- as the population variance, where the denominator is /n/.
+variance :: (G.Vector v Double) => v Double -> Double
+variance samp
+ | n > 1 = robustSumVar (mean samp) samp / fromIntegral n
+ | otherwise = 0
+ where
+ n = G.length samp
+{-# SPECIALIZE variance :: U.Vector Double -> Double #-}
+{-# SPECIALIZE variance :: V.Vector Double -> Double #-}
+
+
+-- | Unbiased estimate of a sample's variance. Also known as the
+-- sample variance, where the denominator is /n/-1.
+varianceUnbiased :: (G.Vector v Double) => v Double -> Double
+varianceUnbiased samp
+ | n > 1 = robustSumVar (mean samp) samp / fromIntegral (n-1)
+ | otherwise = 0
+ where
+ n = G.length samp
+{-# SPECIALIZE varianceUnbiased :: U.Vector Double -> Double #-}
+{-# SPECIALIZE varianceUnbiased :: V.Vector Double -> Double #-}
+
+-- | Standard deviation. This is simply the square root of the
+-- unbiased estimate of the variance.
+stdDev :: (G.Vector v Double) => v Double -> Double
+stdDev = sqrt . varianceUnbiased
+{-# SPECIALIZE stdDev :: U.Vector Double -> Double #-}
+{-# SPECIALIZE stdDev :: V.Vector Double -> Double #-}
+
+-- $cancellation
+--
+-- The functions prefixed with the name @fast@ below perform a single
+-- pass over the sample data using Knuth's algorithm. They usually
+-- work well, but see below for caveats. These functions are subject
+-- to array fusion.
+--
+-- /Note/: in cases where most sample data is close to the sample's
+-- mean, Knuth's algorithm gives inaccurate results due to
+-- catastrophic cancellation.
+
+-- $references
+--
+-- * Chan, T. F.; Golub, G.H.; LeVeque, R.J. (1979) Updating formulae
+-- and a pairwise algorithm for computing sample
+-- variances. Technical Report STAN-CS-79-773, Department of
+-- Computer Science, Stanford
+-- University.
+--
+-- * Knuth, D.E. (1998) The art of computer programming, volume 2:
+-- seminumerical algorithms, 3rd ed., p. 232.
+--
+-- * Welford, B.P. (1962) Note on a method for calculating corrected
+-- sums of squares and products. /Technometrics/
+-- 4(3):419–420.
+--
+-- * West, D.H.D. (1979) Updating mean and variance estimates: an
+-- improved method. /Communications of the ACM/
+-- 22(9):532–535.
diff --git a/statistics/Statistics/Sample/Histogram.hs b/statistics/Statistics/Sample/Histogram.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Sample/Histogram.hs
@@ -0,0 +1,56 @@
+{-# LANGUAGE FlexibleContexts, BangPatterns #-}
+
+-- |
+-- Module : Statistics.Sample.Histogram
+-- Copyright : (c) 2011 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Functions for computing histograms of sample data.
+
+module Statistics.Sample.Histogram
+ (
+ -- * Building blocks
+ histogram_
+ ) where
+
+import Numeric.MathFunctions.Constants (m_epsilon)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as GM
+
+-- | /O(n)/ Compute a histogram over a data set.
+--
+-- Interval (bin) sizes are uniform, based on the supplied upper
+-- and lower bounds.
+histogram_ :: (Num b, RealFrac a, G.Vector v0 a, G.Vector v1 b) =>
+ Int
+ -- ^ Number of bins. This value must be positive. A zero
+ -- or negative value will cause an error.
+ -> a
+ -- ^ Lower bound on interval range. Sample data less than
+ -- this will cause an error.
+ -> a
+ -- ^ Upper bound on interval range. This value must not be
+ -- less than the lower bound. Sample data that falls above
+ -- the upper bound will cause an error.
+ -> v0 a
+ -- ^ Sample data.
+ -> v1 b
+histogram_ numBins lo hi xs0 = G.create (GM.replicate numBins 0 >>= bin xs0)
+ where
+ bin xs bins = go 0
+ where
+ go i | i >= len = return bins
+ | otherwise = do
+ let x = xs `G.unsafeIndex` i
+ b = truncate $ (x - lo) / d
+ write' bins b . (+1) =<< GM.read bins b
+ go (i+1)
+ write' bins b !e = GM.write bins b e
+ len = G.length xs
+ d = ((hi - lo) * (1 + realToFrac m_epsilon)) / fromIntegral numBins
+{-# INLINE histogram_ #-}
+
diff --git a/statistics/Statistics/Sample/Internal.hs b/statistics/Statistics/Sample/Internal.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Sample/Internal.hs
@@ -0,0 +1,30 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module : Statistics.Sample.Internal
+-- Copyright : (c) 2013 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Internal functions for computing over samples.
+module Statistics.Sample.Internal
+ (
+ robustSumVar
+ , sum
+ ) where
+
+import Numeric.Sum (kbn, sumVector)
+import Prelude hiding (sum)
+import Statistics.Function (square)
+import qualified Data.Vector.Generic as G
+
+robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double
+robustSumVar m = sum . G.map (square . subtract m)
+{-# INLINE robustSumVar #-}
+
+sum :: (G.Vector v Double) => v Double -> Double
+sum = sumVector kbn
+{-# INLINE sum #-}
diff --git a/statistics/Statistics/Sample/KernelDensity.hs b/statistics/Statistics/Sample/KernelDensity.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Sample/KernelDensity.hs
@@ -0,0 +1,123 @@
+{-# LANGUAGE BangPatterns, FlexibleContexts, UnboxedTuples #-}
+-- |
+-- Module : Statistics.Sample.KernelDensity
+-- Copyright : (c) 2011 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Kernel density estimation. This module provides a fast, robust,
+-- non-parametric way to estimate the probability density function of
+-- a sample.
+--
+-- This estimator does not use the commonly employed \"Gaussian rule
+-- of thumb\". As a result, it outperforms many plug-in methods on
+-- multimodal samples with widely separated modes.
+
+module Statistics.Sample.KernelDensity
+ (
+ -- * Estimation functions
+ kde
+ -- , kde_
+ -- * References
+ -- $references
+ ) where
+
+import Numeric.MathFunctions.Constants (m_sqrt_2_pi)
+import Prelude hiding (const, min, max, sum)
+import Statistics.Function (minMax, nextHighestPowerOfTwo)
+import Statistics.Math.RootFinding (fromRoot, ridders)
+import Statistics.Sample.Histogram (histogram_)
+import Statistics.Sample.Internal (sum)
+import Statistics.Transform (CD, dct, idct)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector as V
+
+
+-- | Gaussian kernel density estimator for one-dimensional data, using
+-- the method of Botev et al.
+--
+-- The result is a pair of vectors, containing:
+--
+-- * The coordinates of each mesh point. The mesh interval is chosen
+-- to be 20% larger than the range of the sample. (To specify the
+-- mesh interval, use 'kde_'.)
+--
+-- * Density estimates at each mesh point.
+kde :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
+ => Int
+ -- ^ The number of mesh points to use in the uniform discretization
+ -- of the interval @(min,max)@. If this value is not a power of
+ -- two, then it is rounded up to the next power of two.
+ -> v Double -> (v Double, v Double)
+kde n0 xs = kde_ n0 (lo - range / 10) (hi + range / 10) xs
+ where
+ (lo,hi) = minMax xs
+ range | G.length xs <= 1 = 1 -- Unreasonable guess
+ | lo == hi = 1 -- All elements are equal
+ | otherwise = hi - lo
+{-# INLINABLE kde #-}
+{-# SPECIAlIZE kde :: Int -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}
+{-# SPECIAlIZE kde :: Int -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}
+
+
+-- | Gaussian kernel density estimator for one-dimensional data, using
+-- the method of Botev et al.
+--
+-- The result is a pair of vectors, containing:
+--
+-- * The coordinates of each mesh point.
+--
+-- * Density estimates at each mesh point.
+kde_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int)
+ => Int
+ -- ^ The number of mesh points to use in the uniform discretization
+ -- of the interval @(min,max)@. If this value is not a power of
+ -- two, then it is rounded up to the next power of two.
+ -> Double
+ -- ^ Lower bound (@min@) of the mesh range.
+ -> Double
+ -- ^ Upper bound (@max@) of the mesh range.
+ -> v Double
+ -> (v Double, v Double)
+kde_ n0 min max xs
+ | G.null xs = error "Statistics.KernelDensity.kde: empty sample"
+ | n0 <= 1 = error "Statistics.KernelDensity.kde: invalid number of points"
+ | otherwise = (mesh, density)
+ where
+ mesh = G.generate ni $ \z -> min + (d * fromIntegral z)
+ where d = r / (n-1)
+ density = G.map (/(2 * r)) . idct $ G.zipWith f a (G.enumFromTo 0 (n-1))
+ where f b z = b * exp (sqr z * sqr pi * t_star * (-0.5))
+ !n = fromIntegral ni
+ !ni = nextHighestPowerOfTwo n0
+ !r = max - min
+ a = dct . G.map (/ sum h) $ h
+ where h = G.map (/ len) $ histogram_ ni min max xs
+ !len = fromIntegral (G.length xs)
+ !t_star = fromRoot (0.28 * len ** (-0.4)) . ridders 1e-14 (0,0.1) $ \x ->
+ x - (len * (2 * sqrt pi) * go 6 (f 7 x)) ** (-0.4)
+ where
+ f q t = 2 * pi ** (q*2) * sum (G.zipWith g iv a2v)
+ where g i a2 = i ** q * a2 * exp ((-i) * sqr pi * t)
+ a2v = G.map (sqr . (*0.5)) $ G.tail a
+ iv = G.map sqr $ G.enumFromTo 1 (n-1)
+ go s !h | s == 1 = h
+ | otherwise = go (s-1) (f s time)
+ where time = (2 * const * k0 / len / h) ** (2 / (3 + 2 * s))
+ const = (1 + 0.5 ** (s+0.5)) / 3
+ k0 = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi
+ sqr x = x * x
+{-# INLINABLE kde_ #-}
+{-# SPECIAlIZE kde_ :: Int -> Double -> Double -> U.Vector Double -> (U.Vector Double, U.Vector Double) #-}
+{-# SPECIAlIZE kde_ :: Int -> Double -> Double -> V.Vector Double -> (V.Vector Double, V.Vector Double) #-}
+
+
+-- $references
+--
+-- Botev. Z.I., Grotowski J.F., Kroese D.P. (2010). Kernel density
+-- estimation via diffusion. /Annals of Statistics/
+-- 38(5):2916–2957.
diff --git a/statistics/Statistics/Transform.hs b/statistics/Statistics/Transform.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Transform.hs
@@ -0,0 +1,155 @@
+{-# LANGUAGE BangPatterns, FlexibleContexts #-}
+-- |
+-- Module : Statistics.Transform
+-- Copyright : (c) 2011 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Fourier-related transformations of mathematical functions.
+--
+-- These functions are written for simplicity and correctness, not
+-- speed. If you need a fast FFT implementation for your application,
+-- you should strongly consider using a library of FFTW bindings
+-- instead.
+
+module Statistics.Transform
+ (
+ -- * Type synonyms
+ CD
+ -- * Discrete cosine transform
+ , dct
+ , idct
+ ) where
+
+import Control.Monad (when)
+import Control.Monad.ST (ST)
+import Data.Bits (shiftL, shiftR)
+import Data.Complex (Complex(..), conjugate, realPart)
+import Numeric.SpecFunctions (log2)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as M
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector as V
+
+type CD = Complex Double
+
+-- | Discrete cosine transform (DCT-II).
+dct :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v Double -> v Double
+dct = dctWorker . G.map (:+0)
+{-# INLINABLE dct #-}
+{-# SPECIAlIZE dct :: U.Vector Double -> U.Vector Double #-}
+{-# SPECIAlIZE dct :: V.Vector Double -> V.Vector Double #-}
+
+dctWorker :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v CD -> v Double
+{-# INLINE dctWorker #-}
+dctWorker xs
+ -- length 1 is special cased because shuffle algorithms fail for it.
+ | G.length xs == 1 = G.map ((2*) . realPart) xs
+ | vectorOK xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)
+ | otherwise = error "Statistics.Transform.dct: bad vector length"
+ where
+ interleaved = G.backpermute xs $ G.enumFromThenTo 0 2 (len-2) G.++
+ G.enumFromThenTo (len-1) (len-3) 1
+ weights = G.cons 2 . G.generate (len-1) $ \x ->
+ 2 * exp ((0:+(-1))*fi (x+1)*pi/(2*n))
+ where n = fi len
+ len = G.length xs
+
+
+
+-- | Inverse discrete cosine transform (DCT-III). It's inverse of
+-- 'dct' only up to scale parameter:
+--
+-- > (idct . dct) x = (* length x)
+idct :: (G.Vector v CD, G.Vector v Double) => v Double -> v Double
+idct = idctWorker . G.map (:+0)
+{-# INLINABLE idct #-}
+{-# SPECIAlIZE idct :: U.Vector Double -> U.Vector Double #-}
+{-# SPECIAlIZE idct :: V.Vector Double -> V.Vector Double #-}
+
+idctWorker :: (G.Vector v CD, G.Vector v Double) => v CD -> v Double
+{-# INLINE idctWorker #-}
+idctWorker xs
+ | vectorOK xs = G.generate len interleave
+ | otherwise = error "Statistics.Transform.dct: bad vector length"
+ where
+ interleave z | even z = vals `G.unsafeIndex` halve z
+ | otherwise = vals `G.unsafeIndex` (len - halve z - 1)
+ vals = G.map realPart . ifft $ G.zipWith (*) weights xs
+ weights
+ = G.cons n
+ $ G.generate (len - 1) $ \x -> 2 * n * exp ((0:+1) * fi (x+1) * pi/(2*n))
+ where n = fi len
+ len = G.length xs
+
+
+
+-- | Inverse fast Fourier transform.
+ifft :: G.Vector v CD => v CD -> v CD
+ifft xs
+ | vectorOK xs = G.map ((/fi (G.length xs)) . conjugate) . fft . G.map conjugate $ xs
+ | otherwise = error "Statistics.Transform.ifft: bad vector length"
+{-# INLINABLE ifft #-}
+{-# SPECIAlIZE ifft :: U.Vector CD -> U.Vector CD #-}
+{-# SPECIAlIZE ifft :: V.Vector CD -> V.Vector CD #-}
+
+-- | Radix-2 decimation-in-time fast Fourier transform.
+fft :: G.Vector v CD => v CD -> v CD
+fft v | vectorOK v = G.create $ do mv <- G.thaw v
+ mfft mv
+ return mv
+ | otherwise = error "Statistics.Transform.fft: bad vector length"
+{-# INLINABLE fft #-}
+{-# SPECIAlIZE fft :: U.Vector CD -> U.Vector CD #-}
+{-# SPECIAlIZE fft :: V.Vector CD -> V.Vector CD #-}
+
+-- Vector length must be power of two. It's not checked
+mfft :: (M.MVector v CD) => v s CD -> ST s ()
+{-# INLINE mfft #-}
+mfft vec = bitReverse 0 0
+ where
+ bitReverse i j | i == len-1 = stage 0 1
+ | otherwise = do
+ when (i < j) $ M.swap vec i j
+ let inner k l | k <= l = inner (k `shiftR` 1) (l-k)
+ | otherwise = bitReverse (i+1) (l+k)
+ inner (len `shiftR` 1) j
+ stage l !l1 | l == m = return ()
+ | otherwise = do
+ let !l2 = l1 `shiftL` 1
+ !e = -6.283185307179586/fromIntegral l2
+ flight j !a | j == l1 = stage (l+1) l2
+ | otherwise = do
+ let butterfly i | i >= len = flight (j+1) (a+e)
+ | otherwise = do
+ let i1 = i + l1
+ xi1 :+ yi1 <- M.read vec i1
+ let !c = cos a
+ !s = sin a
+ d = (c*xi1 - s*yi1) :+ (s*xi1 + c*yi1)
+ ci <- M.read vec i
+ M.write vec i1 (ci - d)
+ M.write vec i (ci + d)
+ butterfly (i+l2)
+ butterfly j
+ flight 0 0
+ len = M.length vec
+ m = log2 len
+
+
+----------------------------------------------------------------
+-- Helpers
+----------------------------------------------------------------
+
+fi :: Int -> CD
+fi = fromIntegral
+
+halve :: Int -> Int
+halve = (`shiftR` 1)
+
+vectorOK :: G.Vector v a => v a -> Bool
+{-# INLINE vectorOK #-}
+vectorOK v = (1 `shiftL` log2 n) == n where n = G.length v
diff --git a/statistics/Statistics/Types.hs b/statistics/Statistics/Types.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Types.hs
@@ -0,0 +1,283 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
+-- |
+-- Module : Statistics.Types
+-- Copyright : (c) 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Data types common used in statistics
+module Statistics.Types
+ ( -- * Confidence level
+ CL
+ -- ** Accessors
+ , confidenceLevel
+ , significanceLevel
+ -- ** Constructors
+ , mkCL
+ -- ** Constants and conversion to nσ
+ , cl95
+ -- * Estimates and upper/lower limits
+ , Estimate(..)
+ -- , NormalErr(..)
+ , ConfInt(..)
+ -- ** Constructors
+ -- , estimateNormErr
+ , estimateFromInterval
+ , estimateFromErr
+ -- ** Accessors
+ , confidenceInterval
+ , Scale(..)
+ -- * Other
+ , Sample
+ ) where
+
+import Control.DeepSeq (NFData(..))
+import Data.Data (Data,Typeable)
+import Data.Maybe (fromMaybe)
+import GHC.Generics (Generic)
+
+#if __GLASGOW_HASKELL__ == 704
+import qualified Data.Vector.Generic
+import qualified Data.Vector.Generic.Mutable
+#endif
+
+import Statistics.Internal
+import Statistics.Types.Internal
+
+
+----------------------------------------------------------------
+-- Data type for confidence level
+----------------------------------------------------------------
+
+-- |
+-- Confidence level. In context of confidence intervals it's
+-- probability of said interval covering true value of measured
+-- value. In context of statistical tests it's @1-α@ where α is
+-- significance of test.
+--
+-- Since confidence level are usually close to 1 they are stored as
+-- @1-CL@ internally. There are two smart constructors for @CL@:
+-- 'mkCL' and 'mkCLFromSignificance' (and corresponding variant
+-- returning @Maybe@). First creates @CL@ from confidence level and
+-- second from @1 - CL@ or significance level.
+--
+-- >>> cl95
+-- mkCLFromSignificance 0.05
+--
+-- Prior to 0.14 confidence levels were passed to function as plain
+-- @Doubles@. Use 'mkCL' to convert them to @CL@.
+newtype CL a = CL a
+ deriving (Eq, Typeable, Data, Generic)
+
+instance Show a => Show (CL a) where
+ showsPrec n (CL p) = defaultShow1 "mkCLFromSignificance" p n
+instance (Num a, Ord a, Read a) => Read (CL a) where
+ readPrec = defaultReadPrecM1 "mkCLFromSignificance" mkCLFromSignificanceE
+
+instance NFData a => NFData (CL a) where
+ rnf (CL a) = rnf a
+
+-- |
+-- >>> cl95 > cl90
+-- True
+instance Ord a => Ord (CL a) where
+ CL a < CL b = a > b
+ CL a <= CL b = a >= b
+ CL a > CL b = a < b
+ CL a >= CL b = a <= b
+ max (CL a) (CL b) = CL (min a b)
+ min (CL a) (CL b) = CL (max a b)
+
+
+-- | Create confidence level from probability β or probability
+-- confidence interval contain true value of estimate. Will throw
+-- exception if parameter is out of [0,1] range
+--
+-- >>> mkCL 0.95 -- same as cl95
+-- mkCLFromSignificance 0.05
+mkCL :: (Ord a, Num a) => a -> CL a
+mkCL
+ = fromMaybe (error "Statistics.Types.mkCL: probability is out if [0,1] range")
+ . mkCLE
+
+-- | Same as 'mkCL' but returns @Nothing@ instead of error if
+-- parameter is out of [0,1] range
+--
+-- >>> mkCLE 0.95 -- same as cl95
+-- Just (mkCLFromSignificance 0.05)
+mkCLE :: (Ord a, Num a) => a -> Maybe (CL a)
+mkCLE p
+ | p >= 0 && p <= 1 = Just $ CL (1 - p)
+ | otherwise = Nothing
+
+-- | Same as 'mkCLFromSignificance' but returns @Nothing@ instead of error if
+-- parameter is out of [0,1] range
+--
+-- >>> mkCLFromSignificanceE 0.05 -- same as cl95
+-- Just (mkCLFromSignificance 0.05)
+mkCLFromSignificanceE :: (Ord a, Num a) => a -> Maybe (CL a)
+mkCLFromSignificanceE p
+ | p >= 0 && p <= 1 = Just $ CL p
+ | otherwise = Nothing
+
+-- | Get confidence level. This function is subject to rounding
+-- errors. If @1 - CL@ is needed use 'significanceLevel' instead
+confidenceLevel :: (Num a) => CL a -> a
+confidenceLevel (CL p) = 1 - p
+
+-- | Get significance level.
+significanceLevel :: CL a -> a
+significanceLevel (CL p) = p
+
+
+
+-- | 95% confidence level
+cl95 :: Fractional a => CL a
+cl95 = CL 0.05
+
+----------------------------------------------------------------
+-- Data type for p-value
+----------------------------------------------------------------
+
+-- | Newtype wrapper for p-value.
+newtype PValue a = PValue a
+ deriving (Eq,Ord, Typeable, Data, Generic)
+
+instance Show a => Show (PValue a) where
+ showsPrec n (PValue p) = defaultShow1 "mkPValue" p n
+instance (Num a, Ord a, Read a) => Read (PValue a) where
+ readPrec = defaultReadPrecM1 "mkPValue" mkPValueE
+
+instance NFData a => NFData (PValue a) where
+ rnf (PValue a) = rnf a
+
+
+-- | Construct PValue. Returns @Nothing@ if argument is out of [0,1] range.
+mkPValueE :: (Ord a, Num a) => a -> Maybe (PValue a)
+mkPValueE p
+ | p >= 0 && p <= 1 = Just $ PValue p
+ | otherwise = Nothing
+
+----------------------------------------------------------------
+-- Point estimates
+----------------------------------------------------------------
+
+-- |
+-- A point estimate and its confidence interval. It's parametrized by
+-- both error type @e@ and value type @a@. This module provides two
+-- types of error: 'NormalErr' for normally distributed errors and
+-- 'ConfInt' for error with normal distribution. See their
+-- documentation for more details.
+--
+-- For example @144 ± 5@ (assuming normality) could be expressed as
+--
+-- > Estimate { estPoint = 144
+-- > , estError = NormalErr 5
+-- > }
+--
+-- Or if we want to express @144 + 6 - 4@ at CL95 we could write:
+--
+-- > Estimate { estPoint = 144
+-- > , estError = ConfInt
+-- > { confIntLDX = 4
+-- > , confIntUDX = 6
+-- > , confIntCL = cl95
+-- > }
+--
+-- Prior to statistics 0.14 @Estimate@ data type used following definition:
+--
+-- > data Estimate = Estimate {
+-- > estPoint :: {-# UNPACK #-} !Double
+-- > , estLowerBound :: {-# UNPACK #-} !Double
+-- > , estUpperBound :: {-# UNPACK #-} !Double
+-- > , estConfidenceLevel :: {-# UNPACK #-} !Double
+-- > }
+--
+-- Now type @Estimate ConfInt Double@ should be used instead. Function
+-- 'estimateFromInterval' allow to easily construct estimate from same inputs.
+data Estimate e a = Estimate
+ { estPoint :: !a
+ -- ^ Point estimate.
+ , estError :: !(e a)
+ -- ^ Confidence interval for estimate.
+ } deriving (Eq, Read, Show, Generic
+#if __GLASGOW_HASKELL__ >= 708
+ , Typeable, Data
+#endif
+ )
+
+instance (NFData (e a), NFData a) => NFData (Estimate e a) where
+ rnf (Estimate x dx) = rnf x `seq` rnf dx
+
+
+-- | Confidence interval. It assumes that confidence interval forms
+-- single interval and isn't set of disjoint intervals.
+data ConfInt a = ConfInt
+ { confIntLDX :: !a
+ -- ^ Lower error estimate, or distance between point estimate and
+ -- lower bound of confidence interval.
+ , confIntUDX :: !a
+ -- ^ Upper error estimate, or distance between point estimate and
+ -- upper bound of confidence interval.
+ , confIntCL :: !(CL Double)
+ -- ^ Confidence level corresponding to given confidence interval.
+ }
+ deriving (Read,Show,Eq,Typeable,Data,Generic)
+
+instance NFData a => NFData (ConfInt a) where
+ rnf (ConfInt x y _) = rnf x `seq` rnf y
+
+
+
+----------------------------------------
+-- Constructors
+
+-- | Create estimate with asymmetric error.
+estimateFromErr
+ :: a -- ^ Central estimate
+ -> (a,a) -- ^ Lower and upper errors. Both should be
+ -- positive but it's not checked.
+ -> CL Double -- ^ Confidence level for interval
+ -> Estimate ConfInt a
+estimateFromErr x (ldx,udx) cl = Estimate x (ConfInt ldx udx cl)
+
+-- | Create estimate with asymmetric error.
+estimateFromInterval
+ :: Num a
+ => a -- ^ Point estimate. Should lie within
+ -- interval but it's not checked.
+ -> (a,a) -- ^ Lower and upper bounds of interval
+ -> CL Double -- ^ Confidence level for interval
+ -> Estimate ConfInt a
+estimateFromInterval x (lx,ux) cl
+ = Estimate x (ConfInt (x-lx) (ux-x) cl)
+
+
+----------------------------------------
+-- Accessors
+
+-- | Get confidence interval
+confidenceInterval :: Num a => Estimate ConfInt a -> (a,a)
+confidenceInterval (Estimate x (ConfInt ldx udx _))
+ = (x - ldx, x + udx)
+
+
+-- | Data types which could be multiplied by constant.
+class Scale e where
+ scale :: (Ord a, Num a) => a -> e a -> e a
+
+instance Scale ConfInt where
+ scale a (ConfInt l u cl) | a >= 0 = ConfInt (a*l) (a*u) cl
+ | otherwise = ConfInt (-a*u) (-a*l) cl
+
+instance Scale e => Scale (Estimate e) where
+ scale a (Estimate x dx) = Estimate (a*x) (scale a dx)
+
diff --git a/statistics/Statistics/Types/Internal.hs b/statistics/Statistics/Types/Internal.hs
new file mode 100644
--- /dev/null
+++ b/statistics/Statistics/Types/Internal.hs
@@ -0,0 +1,24 @@
+-- |
+-- Module : Statistics.Types.Internal
+-- Copyright : (c) 2009 Bryan O'Sullivan
+-- License : BSD3
+--
+-- Maintainer : bos@serpentine.com
+-- Stability : experimental
+-- Portability : portable
+--
+-- Types for working with statistics.
+module Statistics.Types.Internal where
+
+
+import qualified Data.Vector.Unboxed as U (Vector)
+
+-- | Sample data.
+type Sample = U.Vector Double
+
+-- | Sample with weights. First element of sample is data, second is weight
+--type WeightedSample = U.Vector (Double,Double)
+
+-- | Weights for affecting the importance of elements of a sample.
+--type Weights = U.Vector Double
+
diff --git a/tests/Cleanup.hs b/tests/Cleanup.hs
new file mode 100644
--- /dev/null
+++ b/tests/Cleanup.hs
@@ -0,0 +1,113 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+import Gauge.Main (Benchmark, bench, nfIO)
+import Gauge.Types (Config(..), Verbosity(Quiet))
+import Control.Applicative (pure)
+import Control.DeepSeq (NFData(..))
+import Control.Exception (Exception, try)
+import Control.Monad (when)
+import Foundation.Monad
+import Data.Typeable (Typeable)
+import System.Directory (doesFileExist, removeFile)
+import System.Environment (withArgs)
+import System.IO ( Handle, IOMode(ReadWriteMode), SeekMode(AbsoluteSeek)
+ , hClose, hFileSize, hSeek, openFile)
+import Test.Tasty (TestTree, defaultMain, testGroup)
+import Test.Tasty.HUnit (testCase)
+import Test.HUnit (assertFailure)
+import qualified Gauge.Main as C
+import Data.ByteString (ByteString)
+import qualified Data.ByteString as BS
+
+instance NFData Handle where
+ rnf !_ = ()
+
+data CheckResult = ShouldThrow | WrongData deriving (Show, Typeable)
+
+instance Exception CheckResult
+
+type BenchmarkWithFile =
+ String -> IO Handle -> (Handle -> IO ()) -> (Handle -> IO ()) -> Benchmark
+
+perRun :: BenchmarkWithFile
+perRun name alloc clean work =
+ bench name $ C.perRunEnvWithCleanup alloc clean work
+
+perBatch :: BenchmarkWithFile
+perBatch name alloc clean work =
+ bench name $ C.perBatchEnvWithCleanup (const alloc) (const clean) work
+
+envWithCleanup :: BenchmarkWithFile
+envWithCleanup name alloc clean work =
+ C.envWithCleanup alloc clean $ bench name . nfIO . work
+
+testCleanup :: Bool -> String -> BenchmarkWithFile -> TestTree
+testCleanup shouldFail name withEnvClean = testCase name $ do
+ existsBefore <- doesFileExist testFile
+ when existsBefore $ failTest "Input file already exists"
+
+ result <- runTest . withEnvClean name alloc clean $ \hnd -> do
+ result <- hFileSize hnd >>= BS.hGet hnd . fromIntegral
+ resetHandle hnd
+ when (result /= testData) $ throw WrongData
+ when shouldFail $ throw ShouldThrow
+
+ case result of
+ Left WrongData -> failTest "Incorrect result read from file"
+ Left ShouldThrow -> return ()
+ Right _ | shouldFail -> failTest "Failed to throw exception"
+ | otherwise -> return ()
+
+ existsAfter <- doesFileExist testFile
+ when existsAfter $ do
+ removeFile testFile
+ failTest "Failed to delete file"
+ where
+ testFile :: String
+ testFile = "tmp"
+
+ testData :: ByteString
+ testData = "blah"
+
+ runTest :: Benchmark -> IO (Either CheckResult ())
+ runTest = withArgs (["-n","1"]) . try . C.defaultMainWith config . pure
+ where
+ config = C.defaultConfig { verbosity = Quiet , timeLimit = 1 }
+
+ failTest :: String -> IO ()
+ failTest s = assertFailure $ s ++ " in test: " ++ name ++ "!"
+
+ resetHandle :: Handle -> IO ()
+ resetHandle hnd = hSeek hnd AbsoluteSeek 0
+
+ alloc :: IO Handle
+ alloc = do
+ hnd <- openFile testFile ReadWriteMode
+ BS.hPut hnd testData
+ resetHandle hnd
+ return hnd
+
+ clean :: Handle -> IO ()
+ clean hnd = do
+ hClose hnd
+ removeFile testFile
+
+testSuccess :: String -> BenchmarkWithFile -> TestTree
+testSuccess = testCleanup False
+
+testFailure :: String -> BenchmarkWithFile -> TestTree
+testFailure = testCleanup True
+
+main :: IO ()
+main = defaultMain $ testGroup "cleanup"
+ [ testSuccess "perRun Success" perRun
+ , testFailure "perRun Failure" perRun
+ , testSuccess "perBatch Success" perBatch
+ , testFailure "perBatch Failure" perBatch
+ , testSuccess "envWithCleanup Success" envWithCleanup
+ , testFailure "envWithCleanup Failure" envWithCleanup
+ ]
diff --git a/tests/Properties.hs b/tests/Properties.hs
new file mode 100644
--- /dev/null
+++ b/tests/Properties.hs
@@ -0,0 +1,42 @@
+{-# LANGUAGE CPP #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module Properties (tests) where
+
+import Control.Applicative as A ((<$>))
+import Gauge.Analysis
+import Statistics.Types (Sample)
+import Test.Tasty (TestTree, testGroup)
+import Test.Tasty.QuickCheck (testProperty)
+import Test.QuickCheck
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+#if __GLASGOW_HASKELL__ >= 704
+import Data.Monoid ((<>))
+#else
+import Data.Monoid
+
+(<>) :: Monoid m => m -> m -> m
+(<>) = mappend
+infixr 6 <>
+#endif
+
+instance (Arbitrary a, U.Unbox a) => Arbitrary (U.Vector a) where
+ arbitrary = U.fromList A.<$> arbitrary
+ shrink = map U.fromList . shrink . U.toList
+
+outlier_bucketing :: Double -> Sample -> Bool
+outlier_bucketing y ys =
+ countOutliers (classifyOutliers xs) <= fromIntegral (G.length xs)
+ where xs = U.cons y ys
+
+outlier_bucketing_weighted :: Double -> Sample -> Bool
+outlier_bucketing_weighted x xs =
+ outlier_bucketing x (xs <> G.replicate (G.length xs * 10) 0)
+
+tests :: TestTree
+tests = testGroup "Properties" [
+ testProperty "outlier_bucketing" outlier_bucketing
+ , testProperty "outlier_bucketing_weighted" outlier_bucketing_weighted
+ ]
diff --git a/tests/Sanity.hs b/tests/Sanity.hs
new file mode 100644
--- /dev/null
+++ b/tests/Sanity.hs
@@ -0,0 +1,65 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+import Gauge.Main (bench, bgroup, env, whnf)
+import System.Environment (getEnv, withArgs)
+import System.Timeout (timeout)
+import Test.Tasty (defaultMain)
+import Test.Tasty.HUnit (testCase)
+import Test.HUnit (Assertion, assertFailure)
+import qualified Gauge.Main as C
+import qualified Control.Exception as E
+import qualified Data.ByteString as B
+
+#if !MIN_VERSION_bytestring(0,10,0)
+import Control.DeepSeq (NFData (..))
+#endif
+
+fib :: Int -> Int
+fib = sum . go
+ where go 0 = [0]
+ go 1 = [1]
+ go n = go (n-1) ++ go (n-2)
+
+-- Additional arguments to include along with the ARGS environment variable.
+extraArgs :: [String]
+extraArgs = [ "--raw=sanity.dat", "--json=sanity.json", "--csv=sanity.csv"
+ , "--output=sanity.html", "--junit=sanity.junit" ]
+
+sanity :: Assertion
+sanity = do
+ args <- getArgEnv
+ withArgs (extraArgs ++ args) $ do
+ let tooLong = 30
+ wat <- timeout (tooLong * 1000000) $
+ C.defaultMain [
+ bgroup "fib" [
+ bench "fib 10" $ whnf fib 10
+ , bench "fib 22" $ whnf fib 22
+ ]
+ , env (return (replicate 1024 0)) $ \xs ->
+ bgroup "length . filter" [
+ bench "string" $ whnf (length . filter (==0)) xs
+ , env (return (B.pack xs)) $ \bs ->
+ bench "bytestring" $ whnf (B.length . B.filter (==0)) bs
+ ]
+ ]
+ case wat of
+ Just () -> return ()
+ Nothing -> assertFailure $ "killed for running longer than " ++
+ show tooLong ++ " seconds!"
+
+main :: IO ()
+main = defaultMain $ testCase "sanity" sanity
+
+-- This is a workaround to in pass arguments that sneak past
+-- test-framework to get to criterion.
+getArgEnv :: IO [String]
+getArgEnv =
+ fmap words (getEnv "ARGS") `E.catch`
+ \(_ :: E.SomeException) -> return []
+
+#if !MIN_VERSION_bytestring(0,10,0)
+instance NFData B.ByteString where
+ rnf bs = bs `seq` ()
+#endif
diff --git a/tests/Tests.hs b/tests/Tests.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests.hs
@@ -0,0 +1,9 @@
+module Main (main) where
+
+import Properties
+import Test.Tasty (defaultMain, testGroup)
+
+main :: IO ()
+main = defaultMain $ testGroup "Tests"
+ [ Properties.tests
+ ]