packages feed

gauge (empty) → 0.1.0

raw patch · 44 files changed

+5165/−0 lines, 44 filesdep +HUnitdep +QuickCheckdep +ansi-wl-pprintsetup-changed

Dependencies added: HUnit, QuickCheck, ansi-wl-pprint, base, basement, bytestring, code-page, containers, deepseq, directory, foundation, gauge, math-functions, mwc-random, optparse-applicative, statistics, tasty, tasty-hunit, tasty-quickcheck, vector

Files

+ Gauge.hs view
@@ -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
+ Gauge/Analysis.hs view
@@ -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"
+ Gauge/IO/Printf.hs view
@@ -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
+ Gauge/Internal.hs view
@@ -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 = "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n" +++          printf "<testsuite name=\"Gauge benchmarks\" tests=\"%d\">\n"+          (length rs) +++          concatMap single rs +++          "</testsuite>\n"+    single Report{..} = printf "  <testcase name=\"%s\" time=\"%f\" />\n"+               (attrEsc reportName) (estPoint $ anMean $ reportAnalysis)+    attrEsc = concatMap esc+      where+        esc '\'' = "&apos;"+        esc '"'  = "&quot;"+        esc '<'  = "&lt;"+        esc '>'  = "&gt;"+        esc '&'  = "&amp;"+        esc c    = [c]+-}
+ Gauge/Main.hs view
@@ -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!
+ Gauge/Main/Options.hs view
@@ -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)
+ Gauge/Measurement.hs view
@@ -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
+ Gauge/Monad.hs view
@@ -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
+ Gauge/Monad/ExceptT.hs view
@@ -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))+
+ Gauge/Monad/Internal.hs view
@@ -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)
+ Gauge/Types.hs view
@@ -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&#0178; 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
+ LICENSE view
@@ -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.
+ README.markdown view
@@ -0,0 +1,1 @@+# Gauge: a clone of criterion
+ Setup.lhs view
@@ -0,0 +1,3 @@+#!/usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ cbits/cycles.c view
@@ -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 <unistd.h>+#include <asm-generic/unistd.h>+#include <linux/perf_event.h>++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
+ cbits/time-osx.c view
@@ -0,0 +1,35 @@+#include <mach/mach.h>+#include <mach/mach_time.h>++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));+}
+ cbits/time-posix.c view
@@ -0,0 +1,24 @@+#include <time.h>++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;+}
+ cbits/time-windows.c view
@@ -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 <windows.h>++#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;+}
+ changelog.md view
@@ -0,0 +1,4 @@+# 0.1.0++* remove bunch of dependencies+* initial import of criterion-1.2.2.0
+ gauge.cabal view
@@ -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 <bos@serpentine.com>+maintainer:     Vincent Hanquez <vincent@snarc.org>+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
+ statistics/Statistics/Distribution.hs view
@@ -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/&#8804;/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+/&#948;/x/) equal to /density(x)/&#8901;&#948;/x/+    density :: d -> Double -> Double+    density d = exp . logDensity d++    -- | Inverse of the cumulative distribution function. The value+    -- /x/ for which P(/X/&#8804;/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
+ statistics/Statistics/Distribution/Normal.hs view
@@ -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
+ statistics/Statistics/Function.hs view
@@ -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 #-}
+ statistics/Statistics/Internal.hs view
@@ -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)
+ statistics/Statistics/Math/RootFinding.hs view
@@ -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&#8211;980.
+ statistics/Statistics/Matrix.hs view
@@ -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
+ statistics/Statistics/Matrix/Algorithms.hs view
@@ -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)
+ statistics/Statistics/Matrix/Mutable.hs view
@@ -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 #-}
+ statistics/Statistics/Matrix/Types.hs view
@@ -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
+ statistics/Statistics/Quantile.hs view
@@ -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/.+--   <http://mathworld.wolfram.com/Quantile.html>+--+-- * Hyndman, R.J.; Fan, Y. (1996) Sample quantiles in statistical+--   packages. /American Statistician/+--   50(4):361&#8211;365. <http://www.jstor.org/stable/2684934>
+ statistics/Statistics/Regression.hs view
@@ -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&#0178;/, 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&#0178;/, 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
+ statistics/Statistics/Resampling.hs view
@@ -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))
+ statistics/Statistics/Resampling/Bootstrap.hs view
@@ -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 (<pt) $ resample+        accel = sumCubes / (6 * (sumSquares ** 1.5))+          where (sumSquares :< sumCubes) = U.foldl' f (0 :< 0) jack+                f (s :< c) j = s + d2 :< c + d2 * d+                    where d  = jackMean - j+                          d2 = d * d+                jackMean     = mean jack+        jack  = jackknife est sample++-- $references+--+-- * Davison, A.C; Hinkley, D.V. (1997) Bootstrap methods and their+--   application. <http://statwww.epfl.ch/davison/BMA/>
+ statistics/Statistics/Sample.hs view
@@ -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&#8212;and hence the standard deviation&#8212;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. <ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>+--+-- * 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&#8211;420. <http://www.jstor.org/stable/1266577>+--+-- * West, D.H.D. (1979) Updating mean and variance estimates: an+--   improved method. /Communications of the ACM/+--   22(9):532&#8211;535. <http://doi.acm.org/10.1145/359146.359153>
+ statistics/Statistics/Sample/Histogram.hs view
@@ -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_ #-}+
+ statistics/Statistics/Sample/Internal.hs view
@@ -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 #-}
+ statistics/Statistics/Sample/KernelDensity.hs view
@@ -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&#8211;2957. <http://arxiv.org/pdf/1011.2602>
+ statistics/Statistics/Transform.hs view
@@ -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
+ statistics/Statistics/Types.hs view
@@ -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)+
+ statistics/Statistics/Types/Internal.hs view
@@ -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+
+ tests/Cleanup.hs view
@@ -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+    ]
+ tests/Properties.hs view
@@ -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+  ]
+ tests/Sanity.hs view
@@ -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
+ tests/Tests.hs view
@@ -0,0 +1,9 @@+module Main (main) where++import Properties+import Test.Tasty (defaultMain, testGroup)++main :: IO ()+main = defaultMain $ testGroup "Tests"+       [ Properties.tests+       ]