diff --git a/Gauge/Analysis.hs b/Gauge/Analysis.hs --- a/Gauge/Analysis.hs +++ b/Gauge/Analysis.hs @@ -37,20 +37,18 @@ , benchmarkWith' ) where --- Temporary: to support pre-AMP GHC 7.8.4: -import Data.Monoid - import Control.Applicative import Control.Arrow (second) import Control.DeepSeq (NFData(rnf)) import Control.Monad (forM_, when) -import Gauge.Benchmark (Benchmark (..), Benchmarkable, runBenchmark) +import Gauge.Benchmark (Benchmark (..), BenchmarkAnalysis(..), Benchmarkable, runBenchmark) import Gauge.IO.Printf (note, printError, prolix, rewindClearLine) import Gauge.Main.Options (defaultConfig, Config(..), Verbosity (..), DisplayMode (..)) import Gauge.Measurement (Measured(measTime), secs, rescale, measureKeys, measureAccessors_, validateAccessors, renderNames) import Gauge.Monad (Gauge, askConfig, gaugeIO, Crit(..), askCrit, withConfig) +import Gauge.Format import qualified Gauge.CSV as CSV import Data.Data (Data, Typeable) import Data.Int (Int64) @@ -106,9 +104,8 @@ instance NFData OutlierEffect -instance Monoid Outliers where - mempty = Outliers 0 0 0 0 0 - mappend = addOutliers +outliersEmpty :: Outliers +outliersEmpty = Outliers 0 0 0 0 0 addOutliers :: Outliers -> Outliers -> Outliers addOutliers (Outliers s a b c d) (Outliers t w x y z) = @@ -192,9 +189,10 @@ rnf reportName `seq` rnf reportKeys `seq` rnf reportMeasured `seq` rnf reportAnalysis `seq` rnf reportOutliers `seq` rnf reportKDEs + -- | Classify outliers in a data set, using the boxplot technique. classifyOutliers :: Sample -> Outliers -classifyOutliers sa = U.foldl' ((. outlier) . mappend) mempty ssa +classifyOutliers sa = U.foldl' ((. outlier) . addOutliers) outliersEmpty ssa where outlier e = Outliers { samplesSeen = 1 , lowSevere = if e <= loS && e < hiM then 1 else 0 @@ -412,7 +410,7 @@ case displayMode of StatsTable -> do - _ <- note "%sbenchmarked %s\n" rewindClearLine desc + _ <- note "%sbenchmarked %s%s%s\n" rewindClearLine green desc reset let r2 n = printf "%.3f R\178" n forM_ builtin $ \Regression{..} -> case Map.lookup "iters" regCoeffs of @@ -436,7 +434,7 @@ _ <- note "\n" pure () Condensed -> do - _ <- note "%s%-40s " rewindClearLine desc + _ <- note "%s%s%-40s%s " rewindClearLine green desc reset bsSmall secs "mean" anMean bsSmall secs "( +-" anStdDev _ <- note ")\n" @@ -445,8 +443,9 @@ return 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) + note "%s%-20s%s %s%-10s%s (%s .. %s%s)\n" yellow metric reset + red (f estPoint) reset + (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e) (let cl = confIntCL estError str | cl == cl95 = "" | otherwise = printf ", ci %.3f" (confidenceLevel cl) @@ -454,7 +453,7 @@ ) bsSmall :: (Double -> String) -> String -> Estimate ConfInt Double -> Gauge () bsSmall f metric Estimate{..} = - note "%s %-10s" metric (f estPoint) + note "%s %s%-10s%s" metric red (f estPoint) reset reportStat :: Verbosity -> (Measured -> Maybe Double) @@ -487,9 +486,8 @@ -- | Run a benchmark interactively and analyse its performance. benchmarkWith' :: Config -> Benchmarkable -> IO () -benchmarkWith' cfg bm = - withConfig cfg $ - runBenchmark (const True) (Benchmark "function" bm) analyseBenchmark +benchmarkWith' cfg bm = withConfig cfg $ + runBenchmark (const True) (Benchmark "function" bm) (BenchmarkNormal analyseBenchmark) -- | Run a benchmark interactively and analyse its performanc. benchmark' :: Benchmarkable -> IO () diff --git a/Gauge/Benchmark.hs b/Gauge/Benchmark.hs --- a/Gauge/Benchmark.hs +++ b/Gauge/Benchmark.hs @@ -59,7 +59,7 @@ -- * Running Benchmarks , runBenchmark - , runBenchmarkIters + , BenchmarkAnalysis(..) ) where import Control.Applicative @@ -652,29 +652,25 @@ -- Running benchmarks ------------------------------------------------------------------------------- +-- | The function to run after measurement +data BenchmarkAnalysis = + forall a . BenchmarkNormal (String -> V.Vector Measured -> Gauge a) + | BenchmarkIters Int64 + -- | Run benchmarkables, selected by a given selector function, under a given -- benchmark and analyse the output using the given analysis function. -runBenchmark - :: (String -> Bool) -- ^ Select benchmarks by name. - -> Benchmark - -> (String -> V.Vector Measured -> Gauge a) -- ^ Analysis function. - -> Gauge () +runBenchmark :: (String -> Bool) -- ^ Select benchmarks by name. + -> Benchmark + -> BenchmarkAnalysis -- ^ Analysis function + -> Gauge () runBenchmark selector bs analyse = - for selector bs $ \_idx desc bm -> - runBenchmarkable desc bm >>= analyse desc >>= \_ -> return () - --- XXX For consistency, this should also use a separate process when --- --measure-with is specified. --- | Run a benchmark without analysing its performance. -runBenchmarkIters - :: (String -> Bool) -- ^ Select benchmarks by name. - -> Benchmark - -> Int64 -- ^ Number of iterations to run. - -> Gauge () -runBenchmarkIters selector bs iters = - for selector bs $ \_idx desc bm -> do - _ <- note "benchmarking %s\r" desc - gaugeIO $ iterateBenchmarkable_ bm iters + for selector bs $ \_idx desc bm -> + case analyse of + BenchmarkNormal step -> + runBenchmarkable desc bm >>= step desc >>= \_ -> return () + BenchmarkIters iters -> do + _ <- note "benchmarking %s\r" desc + gaugeIO $ iterateBenchmarkable_ bm iters -- | Iterate over benchmarks. for :: (String -> Bool) diff --git a/Gauge/Format.hs b/Gauge/Format.hs new file mode 100644 --- /dev/null +++ b/Gauge/Format.hs @@ -0,0 +1,119 @@ +-- | +-- Module : Gauge.Format +-- Copyright : (c) 2017 Vincent Hanquez +-- +-- Formatting helpers +-- +-- shame there's no leftPad package to use. /s +-- +module Gauge.Format + ( printNanoseconds + , printSubNanoseconds + , tableMarkdown + , reset + , green + , red + , yellow + ) where + +import Gauge.Time +import Data.List +import Data.Word +import Text.Printf +import qualified Basement.Terminal.ANSI as ANSI +import Basement.Bounded (zn64) +import GHC.Exts (toList) + +-- | Print a NanoSeconds quantity with a human friendly format +-- that make it easy to compare different values +-- +-- Given a separator Char of '_': +-- +-- 0 -> " 0" +-- 1000 -> " 1_000" +-- 1234567 -> " 1_234_567" +-- 10200300400 -> "10_200_300_400" +-- +-- Note that the seconds parameters is aligned considered +-- maximum of 2 characters (i.e. 99 seconds). +-- +printNanoseconds :: Maybe Char -> NanoSeconds -> String +printNanoseconds thousandSeparator (NanoSeconds absNs) = + case divSub1000 0 absNs of + [ns] -> padLeft maxLength $ printSpace ns + [ns,us] -> padLeft maxLength $ addSeparators1000 [printSpace us,print3 ns] + [ns,us,ms] -> padLeft maxLength $ addSeparators1000 [printSpace ms,print3 us,print3 ns] + [ns,us,ms,s] -> padLeft maxLength $ addSeparators1000 [printSpace s,print3 ms,print3 us,print3 ns] + _ -> error "printNanoSeconds: internal error: invalid format" + where + maxLength = 3 + 3 + 3 + 2 + (sepLength * 3) + + (addSeparators1000, sepLength) = + case thousandSeparator of + Nothing -> (concat, 0) + Just c -> (intercalate [c], 1) + + printSpace :: Word64 -> String + printSpace n = printf "%3d" n + print3 :: Word64 -> String + print3 n = printf "%03d" n + + divSub1000 :: Int -> Word64 -> [Word64] + divSub1000 n i + | n == 3 = [i] + | otherwise = + let (d,m) = i `divMod` 1000 + in if d == 0 then [m] else m : divSub1000 (n+1) d + +printSubNanoseconds :: Maybe Char -> PicoSeconds100 -> String +printSubNanoseconds ts p = + printNanoseconds ts ns ++ "." ++ show fragment + where + (ns, fragment) = picosecondsToNanoSeconds p + + +-- | Produce a table in markdown +-- +-- This is handy when wanting to copy paste to a markdown flavor destination. +tableMarkdown :: String -- ^ top left corner label + -> [String] -- ^ columns labels + -> [[String]] -- ^ a list of row labels followed by content rows + -> String -- ^ the resulting string +tableMarkdown name cols rows = + let hdr = "| " ++ intercalate " | " (padList (name : cols)) ++ " |\n" + sep = "|-" ++ intercalate "-|-" (map (map (const '-')) (padList (name : cols))) ++ "-|\n" + in hdr ++ sep ++ concatMap printRow (map padList rows) + where + printRow :: [String] -> String + printRow l = "| " ++ intercalate " | " l ++ " |\n" + + getColN n = map (flip (!!) n) rows + + sizeCols :: [Int] + sizeCols = map (\(i, c) -> maximum $ map length (c : getColN i)) $ zip [0..] (name : cols) + + padList l = zipWith padCenter sizeCols l + +padLeft :: Int -> String -> String +padLeft sz s + | sz <= len = s + | otherwise = replicate leftPad ' ' ++ s + where + len = length s + leftPad = (sz - len) + +padCenter :: Int -> String -> String +padCenter sz s + | sz <= len = s + | otherwise = replicate leftPad ' ' ++ s ++ replicate rightPad ' ' + where + len = length s + (leftPad, r) = (sz - len) `divMod` 2 + rightPad = leftPad + r + +-- | reset, green, red, yellow ANSI escape +reset, green, red, yellow :: String +reset = toList ANSI.sgrReset +green = toList $ ANSI.sgrForeground (zn64 2) True +red = toList $ ANSI.sgrForeground (zn64 1) True +yellow = toList $ ANSI.sgrForeground (zn64 3) True diff --git a/Gauge/Main.hs b/Gauge/Main.hs --- a/Gauge/Main.hs +++ b/Gauge/Main.hs @@ -23,6 +23,7 @@ -- * Running Benchmarks Interactively , benchmark , benchmarkWith + , module Gauge.Benchmark ) where import Control.Applicative @@ -109,7 +110,7 @@ benchmarkWith :: Config -> Benchmarkable -> IO () benchmarkWith cfg bm = withConfig cfg $ - runBenchmark (const True) (Benchmark "function" bm) quickAnalyse + runBenchmark (const True) (Benchmark "function" bm) (BenchmarkNormal quickAnalyse) -- | Run a benchmark interactively with default config, and analyse its -- performance. @@ -160,35 +161,36 @@ -- one in your benchmark driver's command-line parser). runMode :: Mode -> Config -> [String] -> [Benchmark] -> IO () runMode wat cfg benches bs = - -- TBD: This has become messy. We use mode as well as cfg options for the - -- same purpose It is possible to specify multiple exclusive options. We - -- need to handle the exclusive options in a better way. - case wat of - List -> mapM_ putStrLn . sort . concatMap benchNames $ bs - Version -> putStrLn versionInfo - Help -> putStrLn describe - DefaultMode -> - case measureOnly cfg of - Just outfile -> runWithConfig runBenchmark (\_ r -> - gaugeIO (writeFile outfile (show r))) - Nothing -> - case iters cfg of - Just nbIters -> runWithConfig runBenchmarkIters nbIters - Nothing -> - case quickMode cfg of - True -> runWithConfig runBenchmark quickAnalyse - False -> do + -- TBD: This has become messy. We use mode as well as cfg options for the + -- same purpose It is possible to specify multiple exclusive options. We + -- need to handle the exclusive options in a better way. + case wat of + List -> mapM_ putStrLn . sort . concatMap benchNames $ bs + Version -> putStrLn versionInfo + Help -> putStrLn describe + DefaultMode -> runDefault + where + runDefault = do + CSV.write (csvRawFile cfg) $ CSV.Row $ map (CSV.string . fst) measureAccessors_ + CSV.write (csvFile cfg) $ CSV.Row $ map CSV.string + ["Name", "Mean","MeanLB","MeanUB","Stddev","StddevLB","StddevUB"] + + hSetBuffering stdout NoBuffering + selector <- selectBenches (match cfg) benches bsgroup + + -- if compiled without analysis step, then default to quickmode #ifdef HAVE_ANALYSIS - CSV.write (csvRawFile cfg) $ CSV.Row $ map CSV.string - (map fst measureAccessors_) - CSV.write (csvFile cfg) $ CSV.Row $ map CSV.string - ["Name", "Mean","MeanLB","MeanUB","Stddev","StddevLB","StddevUB"] - runWithConfig runBenchmark analyseBenchmark + let compiledAnalyseStep = analyseBenchmark #else - runWithConfig runBenchmark quickAnalyse + let compiledAnalyseStep = quickAnalyse #endif - where bsgroup = BenchGroup "" bs - runWithConfig f arg = do - hSetBuffering stdout NoBuffering - selector <- selectBenches (match cfg) benches bsgroup - withConfig cfg $ f selector bsgroup arg + + let mode = case (measureOnly cfg, iters cfg, quickMode cfg) of + (Just outfile, _ , _ ) -> BenchmarkNormal $ \_ r -> gaugeIO (writeFile outfile (show r)) + (Nothing , Just nbIters, _ ) -> BenchmarkIters nbIters + (Nothing , Nothing , True) -> BenchmarkNormal quickAnalyse + (Nothing , Nothing , False) -> BenchmarkNormal compiledAnalyseStep + + withConfig cfg $ runBenchmark selector bsgroup mode + + bsgroup = BenchGroup "" bs diff --git a/Gauge/Optional.hs b/Gauge/Optional.hs --- a/Gauge/Optional.hs +++ b/Gauge/Optional.hs @@ -1,4 +1,6 @@ {-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ConstraintKinds #-} -- | -- Module : Gauge.Optional -- Copyright : (c) 2017-2018 Vincent Hanquez @@ -29,6 +31,7 @@ import Data.Word import Data.Data import GHC.Generics +import Basement.Compat.CallStack -- | A type representing a sum-type free Maybe a -- where a specific tag represent Nothing @@ -50,7 +53,7 @@ isOptionalTag d = isInfinite d || isNaN d -- | Create an optional value from a -toOptional :: OptionalTag a => a -> Optional a +toOptional :: (HasCallStack, OptionalTag a) => a -> Optional a toOptional v | isOptionalTag v = error "Creating an optional valid value using the optional tag" | otherwise = Optional v @@ -70,7 +73,7 @@ | otherwise = Just v {-# INLINE toMaybe #-} -fromMaybe :: OptionalTag a => Maybe a -> Optional a +fromMaybe :: (HasCallStack, OptionalTag a) => Maybe a -> Optional a fromMaybe Nothing = Optional optionalTag fromMaybe (Just v) | isOptionalTag v = error "fromMaybe: creating an optional value using the optional tag" @@ -82,7 +85,7 @@ | otherwise = Optional (f v) {-# INLINE map #-} -both :: OptionalTag a => (a -> a -> a) -> Optional a -> Optional a -> Optional a +both :: (HasCallStack, OptionalTag a) => (a -> a -> a) -> Optional a -> Optional a -> Optional a both f o1 o2 | isOmitted o1 = o2 | isOmitted o2 = o1 diff --git a/README.markdown b/README.markdown --- a/README.markdown +++ b/README.markdown @@ -36,7 +36,7 @@ Number of total dependencies (direct & indirect): -* gauge: 18 dependencies +* gauge: 12 dependencies * criterion: 63 dependencies Dependencies removed: @@ -56,30 +56,39 @@ * bytestring 0.10.8.1 * cassava 0.4.5.1 * cereal 0.5.4.0 +* code-page 0.1.3 +* containers 0.5.7.1 * criterion 1.2.2.0 * directory 1.3.0.0 * dlist 0.8.0.3 * erf 2.0.0.0 * exceptions 0.8.3 * filepath 1.4.1.1 +* ghc-boot-th 8.0.2 * hashable 1.2.6.1 +* integer-gmp 1.0.0.1 * integer-logarithms 1.0.2 * js-flot 0.8.3 * js-jquery 3.2.1 +* math-functions 0.2.1.0 * microstache 1.0.1.1 * monad-par 0.3.4.8 * monad-par-extras 0.3.3 * mtl 2.2.1 +* mwc-random 0.13.6.0 * optparse-applicative 0.13.2.0 * parallel 3.2.1.1 * parsec 3.1.11 +* pretty 1.1.3.3 * process 1.4.3.0 * random 1.1 * scientific 0.3.5.2 * statistics 0.14.0.2 * stm 2.4.4.1 * tagged 0.8.5 +* template-haskell 2.11.1.0 * text 1.2.2.2 +* time 1.6.0.1 * time-locale-compat 0.1.1.3 * transformers-compat 0.5.1.4 * unix 2.7.2.1 @@ -87,7 +96,8 @@ * uuid-types 1.0.3 * vector-algorithms 0.7.0.1 * vector-binary-instances 0.2.3.5 -* code-page 0.1.3 +* vector-th-unbox 0.2.1.6 + Criterion graph of dependencies: diff --git a/cbits/cycles.h b/cbits/cycles.h --- a/cbits/cycles.h +++ b/cbits/cycles.h @@ -2,6 +2,7 @@ #define CYCLES_H #include +#include "Rts.h" #if x86_64_HOST_ARCH || i386_HOST_ARCH diff --git a/cbits/time-posix.c b/cbits/time-posix.c --- a/cbits/time-posix.c +++ b/cbits/time-posix.c @@ -2,19 +2,30 @@ #include #include -#include -#include #include "gauge-time.h" +// #define USE_PERF_EVENT + +#ifdef USE_PERF_EVENT +#include +#include +#else +#include +#endif + +#ifdef USE_PERF_EVENT static int gauge_rdtsc_fddev = -1; +#endif void gauge_inittime(void) { +#ifdef USE_PERF_EVENT static struct perf_event_attr attr; attr.type = PERF_TYPE_HARDWARE; attr.config = PERF_COUNT_HW_CPU_CYCLES; gauge_rdtsc_fddev = syscall (__NR_perf_event_open, &attr, 0, -1, -1, 0); +#endif } #define timespec_to_uint64(x) ( \ @@ -25,14 +36,18 @@ void gauge_record(struct gauge_time *tr) { struct timespec ts, ts2; - uint64_t res; + uint64_t res = 0; clock_gettime(CLOCK_MONOTONIC, &ts); clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts2); tr->clock_nanosecs = timespec_to_uint64(ts); tr->cpu_nanosecs = timespec_to_uint64(ts2); +#ifdef USE_PERF_EVENT tr->rdtsc = (read (gauge_rdtsc_fddev, &res, sizeof(res)) < sizeof(res)) ? 0 : res; +#else + tr->rdtsc = instruction_rdtsc(); +#endif } double gauge_gettime(void) diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,15 @@ +# 0.2.1 + +* Inline math-functions & mwc-random: + * Remove most functions, instances and types, that are unnecessary for gauge + * Remove unsafe seeding with partial seed (unused in gauge anyway) + * Remove vector-th-unbox dependency (transitively template-haskell, pretty, ghc-boot-th) + * Remove time dependency +* Re-add Gauge.Benchmark to Gauge.Main to keep the transition between criterion and gauge easy +* Fix cycles reporting on linux, osx and windows +* Add some extra callstack for reporting on partial function +* Fix compilation with Semigroup => Monoid (compilation on 8.4). still unsupported +* Add some color on terminal output # 0.2.0 diff --git a/gauge.cabal b/gauge.cabal --- a/gauge.cabal +++ b/gauge.cabal @@ -1,5 +1,5 @@ name: gauge -version: 0.2.0 +version: 0.2.1 synopsis: small framework for performance measurement and analysis license: BSD3 license-file: LICENSE @@ -19,7 +19,7 @@ GHC==7.8.4, GHC==7.10.3, GHC==8.0.2, - GHC==8.2.1 + GHC==8.2.2 description: This library provides a powerful but simple way to measure software @@ -46,6 +46,7 @@ Gauge.Time Gauge.Optional Gauge.CSV + Gauge.Format Gauge.Source.RUsage Gauge.Source.GC @@ -75,9 +76,16 @@ Statistics.Transform Statistics.Types Statistics.Types.Internal + System.Random.MWC + Numeric.MathFunctions.Comparison + Numeric.MathFunctions.Constants + Numeric.SpecFunctions + Numeric.SpecFunctions.Internal + Numeric.Sum - hs-source-dirs: . statistics + hs-source-dirs: . statistics mwc-random math-functions + include-Dirs: cbits c-sources: cbits/cycles.c if os(darwin) @@ -96,16 +104,10 @@ base >= 4.7 && < 5, basement >= 0.0.4, deepseq >= 1.1.0.0, - mwc-random >= 0.8.0.3, vector >= 0.7.1, process, directory - -- formely statistics dependency that we need - if flag(analysis) - build-depends: - math-functions >= 0.1.7 - default-language: Haskell2010 ghc-options: -O2 -Wall -funbox-strict-fields if flag(analysis) @@ -123,7 +125,6 @@ base, bytestring, gauge, - deepseq, tasty, tasty-hunit @@ -140,7 +141,6 @@ base, bytestring, gauge, - deepseq, tasty, tasty-hunit @@ -157,7 +157,6 @@ base, bytestring, gauge, - deepseq, tasty, tasty-hunit @@ -174,7 +173,6 @@ base, bytestring, gauge, - deepseq, tasty, tasty-hunit diff --git a/math-functions/Numeric/MathFunctions/Comparison.hs b/math-functions/Numeric/MathFunctions/Comparison.hs new file mode 100644 --- /dev/null +++ b/math-functions/Numeric/MathFunctions/Comparison.hs @@ -0,0 +1,57 @@ +-- | +-- Module : Numeric.MathFunctions.Comparison +-- Copyright : (c) 2011 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Functions for approximate comparison of floating point numbers. +-- +-- Approximate floating point comparison, based on Bruce Dawson's +-- \"Comparing floating point numbers\": +-- +module Numeric.MathFunctions.Comparison + ( within + ) where + +import Basement.Floating (doubleToWord) +import Data.Word (Word64) + +-- | +-- Measure distance between two @Double@s in ULPs (units of least +-- precision). Note that it's different from @abs (ulpDelta a b)@ +-- since it returns correct result even when 'ulpDelta' overflows. +ulpDistance :: Double + -> Double + -> Word64 +ulpDistance a b = + -- IEEE754 floats use most significant bit as sign bit (not + -- 2-complement) and we need to rearrange representations of float + -- number so that they could be compared lexicographically as + -- Word64. + let big = 0x8000000000000000 + order i | i < big = i + big + | otherwise = maxBound - i + ai = order ai0 + bi = order bi0 + d | ai > bi = ai - bi + | otherwise = bi - ai + in d + where + ai0 = doubleToWord a + bi0 = doubleToWord b + + +-- | Compare two 'Double' values for approximate equality, using +-- Dawson's method. +-- +-- The required accuracy is specified in ULPs (units of least +-- precision). If the two numbers differ by the given number of ULPs +-- or less, this function returns @True@. +within :: Int -- ^ Number of ULPs of accuracy desired. + -> Double -> Double -> Bool +within ulps a b + | ulps < 0 = False + | otherwise = ulpDistance a b <= fromIntegral ulps diff --git a/math-functions/Numeric/MathFunctions/Constants.hs b/math-functions/Numeric/MathFunctions/Constants.hs new file mode 100644 --- /dev/null +++ b/math-functions/Numeric/MathFunctions/Constants.hs @@ -0,0 +1,115 @@ +-- | +-- Module : Numeric.MathFunctions.Constants +-- Copyright : (c) 2009, 2011 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Constant values common to much numeric code. + +module Numeric.MathFunctions.Constants + ( + -- * IEE754 constants + m_epsilon + , m_huge + , m_tiny + , m_max_exp + , m_pos_inf + , m_neg_inf + , m_NaN + , m_max_log + , m_min_log + -- * Mathematical constants + , m_1_sqrt_2 + , m_2_sqrt_pi + , m_ln_sqrt_2_pi + , m_sqrt_2 + , m_sqrt_2_pi + , m_eulerMascheroni + ) where + +---------------------------------------------------------------- +-- IEE754 constants +---------------------------------------------------------------- + +-- | Largest representable finite value. +m_huge :: Double +m_huge = 1.7976931348623157e308 +{-# INLINE m_huge #-} + +-- | The smallest representable positive normalized value. +m_tiny :: Double +m_tiny = 2.2250738585072014e-308 +{-# INLINE m_tiny #-} + +-- | The largest 'Int' /x/ such that 2**(/x/-1) is approximately +-- representable as a 'Double'. +m_max_exp :: Int +m_max_exp = 1024 + +-- | Positive infinity. +m_pos_inf :: Double +m_pos_inf = 1/0 +{-# INLINE m_pos_inf #-} + +-- | Negative infinity. +m_neg_inf :: Double +m_neg_inf = -1/0 +{-# INLINE m_neg_inf #-} + +-- | Not a number. +m_NaN :: Double +m_NaN = 0/0 +{-# INLINE m_NaN #-} + +-- | Maximum possible finite value of @log x@ +m_max_log :: Double +m_max_log = 709.782712893384 +{-# INLINE m_max_log #-} + +-- | Logarithm of smallest normalized double ('m_tiny') +m_min_log :: Double +m_min_log = -708.3964185322641 +{-# INLINE m_min_log #-} + + +---------------------------------------------------------------- +-- Mathematical constants +---------------------------------------------------------------- + +-- | @sqrt 2@ +m_sqrt_2 :: Double +m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766 +{-# INLINE m_sqrt_2 #-} + +-- | @sqrt (2 * pi)@ +m_sqrt_2_pi :: Double +m_sqrt_2_pi = 2.5066282746310005024157652848110452530069867406099383166299 +{-# INLINE m_sqrt_2_pi #-} + +-- | @2 / sqrt pi@ +m_2_sqrt_pi :: Double +m_2_sqrt_pi = 1.1283791670955125738961589031215451716881012586579977136881 +{-# INLINE m_2_sqrt_pi #-} + +-- | @1 / sqrt 2@ +m_1_sqrt_2 :: Double +m_1_sqrt_2 = 0.7071067811865475244008443621048490392848359376884740365883 +{-# INLINE m_1_sqrt_2 #-} + +-- | The smallest 'Double' ε such that 1 + ε ≠ 1. +m_epsilon :: Double +m_epsilon = encodeFloat (signif+1) expo - 1.0 + where (signif,expo) = decodeFloat (1.0::Double) + +-- | @log(sqrt((2*pi))@ +m_ln_sqrt_2_pi :: Double +m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171 +{-# INLINE m_ln_sqrt_2_pi #-} + +-- | Euler–Mascheroni constant (γ = 0.57721...) +m_eulerMascheroni :: Double +m_eulerMascheroni = 0.5772156649015328606065121 +{-# INLINE m_eulerMascheroni #-} diff --git a/math-functions/Numeric/SpecFunctions.hs b/math-functions/Numeric/SpecFunctions.hs new file mode 100644 --- /dev/null +++ b/math-functions/Numeric/SpecFunctions.hs @@ -0,0 +1,3 @@ +module Numeric.SpecFunctions ( module X ) where + +import Numeric.SpecFunctions.Internal as X diff --git a/math-functions/Numeric/SpecFunctions/Internal.hs b/math-functions/Numeric/SpecFunctions/Internal.hs new file mode 100644 --- /dev/null +++ b/math-functions/Numeric/SpecFunctions/Internal.hs @@ -0,0 +1,116 @@ +{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-} +-- | +-- Module : Numeric.SpecFunctions.Internal +-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Internal module with implementation of special functions. +module Numeric.SpecFunctions.Internal + ( erf + , erfc + , invErf + , invErfc + , log2 + ) where + +import Data.Bits ((.&.), (.|.), shiftR) +import Data.Word (Word64) +import qualified Data.Vector.Unboxed as U + +import Numeric.MathFunctions.Constants + +---------------------------------------------------------------- +-- Error function +---------------------------------------------------------------- + +-- | Error function. +-- +-- \[ +-- \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t^2) dt +-- \] +-- +-- Function limits are: +-- +-- \[ +-- \begin{aligned} +-- &\operatorname{erf}(-\infty) &=& -1 \\ +-- &\operatorname{erf}(0) &=& \phantom{-}\,0 \\ +-- &\operatorname{erf}(+\infty) &=& \phantom{-}\,1 \\ +-- \end{aligned} +-- \] +erf :: Double -> Double +{-# INLINE erf #-} +erf = c_erf + +-- | Complementary error function. +-- +-- \[ +-- \operatorname{erfc}(x) = 1 - \operatorname{erf}(x) +-- \] +-- +-- Function limits are: +-- +-- \[ +-- \begin{aligned} +-- &\operatorname{erf}(-\infty) &=&\, 2 \\ +-- &\operatorname{erf}(0) &=&\, 1 \\ +-- &\operatorname{erf}(+\infty) &=&\, 0 \\ +-- \end{aligned} +-- \] +erfc :: Double -> Double +{-# INLINE erfc #-} +erfc = c_erfc + +foreign import ccall "erf" c_erf :: Double -> Double +foreign import ccall "erfc" c_erfc :: Double -> Double + + +-- | Inverse of 'erf'. +invErf :: Double -- ^ /p/ ∈ [-1,1] + -> Double +invErf p = invErfc (1 - p) + +-- | Inverse of 'erfc'. +invErfc :: Double -- ^ /p/ ∈ [0,2] + -> Double +invErfc p + | p == 2 = m_neg_inf + | p == 0 = m_pos_inf + | p >0 && p < 2 = if p <= 1 then r else -r + | otherwise = modErr $ "invErfc: p must be in [0,2] got " ++ show p + where + pp = if p <= 1 then p else 2 - p + t = sqrt $ -2 * log( 0.5 * pp) + -- Initial guess + x0 = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t) + r = loop 0 x0 + -- + loop :: Int -> Double -> Double + loop !j !x + | j >= 2 = x + | otherwise = let err = erfc x - pp + x' = x + err / (1.12837916709551257 * exp(-x * x) - x * err) -- // Halley + in loop (j+1) x' + +-- | /O(log n)/ Compute the logarithm in base 2 of the given value. +log2 :: Int -> Int +log2 v0 + | v0 <= 0 = modErr $ "log2: nonpositive input, got " ++ show v0 + | otherwise = go 5 0 v0 + where + go !i !r !v | i == -1 = r + | v .&. b i /= 0 = let si = U.unsafeIndex sv i + in go (i-1) (r .|. si) (v `shiftR` si) + | otherwise = go (i-1) r v + b = U.unsafeIndex bv + !bv = U.fromList [ 0x02, 0x0c, 0xf0, 0xff00 + , fromIntegral (0xffff0000 :: Word64) + , fromIntegral (0xffffffff00000000 :: Word64)] + !sv = U.fromList [1,2,4,8,16,32] + +modErr :: String -> a +modErr msg = error $ "Numeric.SpecFunctions." ++ msg diff --git a/math-functions/Numeric/Sum.hs b/math-functions/Numeric/Sum.hs new file mode 100644 --- /dev/null +++ b/math-functions/Numeric/Sum.hs @@ -0,0 +1,198 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts, + MultiParamTypeClasses, TypeFamilies #-} +{-# OPTIONS_GHC -fno-warn-name-shadowing #-} +-- | +-- Module : Numeric.Sum +-- Copyright : (c) 2014 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Functions for summing floating point numbers more accurately than +-- the naive 'Prelude.sum' function and its counterparts in the +-- @vector@ package and elsewhere. +-- +-- When used with floating point numbers, in the worst case, the +-- 'Prelude.sum' function accumulates numeric error at a rate +-- proportional to the number of values being summed. The algorithms +-- in this module implement different methods of /compensated +-- summation/, which reduce the accumulation of numeric error so that +-- it either grows much more slowly than the number of inputs +-- (e.g. logarithmically), or remains constant. +module Numeric.Sum ( + -- * Summation type class + Summation(..) + , sumVector + , kbn + ) where + +import Control.DeepSeq (NFData(..)) +import Control.Monad +import Data.Data (Typeable, Data) +import Data.Vector.Generic (Vector(..), foldl') +import qualified Data.Vector.Generic.Mutable as M + +import qualified Data.Foldable as F +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Unboxed as U + +-- | A class for summation of floating point numbers. +class Summation s where + -- | The identity for summation. + zero :: s + + -- | Add a value to a sum. + add :: s -> Double -> s + + -- | Sum a collection of values. + -- + -- Example: + -- @foo = 'sum' 'kbn' [1,2,3]@ + sum :: (F.Foldable f) => (s -> Double) -> f Double -> Double + sum f = f . F.foldl' add zero + {-# INLINE sum #-} + +instance Summation Double where + zero = 0 + add = (+) + +-- | Kahan-Babuška-Neumaier summation. This is a little more +-- computationally costly than plain Kahan summation, but is /always/ +-- at least as accurate. +data KBNSum = KBNSum {-# UNPACK #-} !Double {-# UNPACK #-} !Double + deriving (Eq, Show, Typeable, Data) + +newtype instance U.MVector s KBNSum = MV_KBNSum (U.MVector s (Double,Double)) +newtype instance U.Vector KBNSum = V_KBNSum (U.Vector (Double,Double)) + +instance M.MVector U.MVector KBNSum where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MV_KBNSum v) = M.basicLength v + basicUnsafeSlice i n (MV_KBNSum v) = MV_KBNSum $ M.basicUnsafeSlice i n v + basicOverlaps (MV_KBNSum v1) (MV_KBNSum v2) = M.basicOverlaps v1 v2 + basicUnsafeNew n = MV_KBNSum `liftM` M.basicUnsafeNew n + basicUnsafeReplicate n (KBNSum a b) = MV_KBNSum `liftM` M.basicUnsafeReplicate n (a,b) + basicUnsafeRead (MV_KBNSum v) i = uncurry KBNSum `liftM` M.basicUnsafeRead v i + basicUnsafeWrite (MV_KBNSum v) i (KBNSum a b) = M.basicUnsafeWrite v i (a,b) + basicClear (MV_KBNSum v) = M.basicClear v + basicSet (MV_KBNSum v) (KBNSum a b) = M.basicSet v (a,b) + basicUnsafeCopy (MV_KBNSum v1) (MV_KBNSum v2) = M.basicUnsafeCopy v1 v2 + basicUnsafeMove (MV_KBNSum v1) (MV_KBNSum v2) = M.basicUnsafeMove v1 v2 + basicUnsafeGrow (MV_KBNSum v) n = MV_KBNSum `liftM` M.basicUnsafeGrow v n +#if MIN_VERSION_vector(0,11,0) + {-# INLINE basicInitialize #-} + basicInitialize (MV_KBNSum v) = M.basicInitialize v +#endif + +instance G.Vector U.Vector KBNSum where + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicUnsafeFreeze (MV_KBNSum v) = V_KBNSum `liftM` G.basicUnsafeFreeze v + basicUnsafeThaw (V_KBNSum v) = MV_KBNSum `liftM` G.basicUnsafeThaw v + basicLength (V_KBNSum v) = G.basicLength v + basicUnsafeSlice i n (V_KBNSum v) = V_KBNSum $ G.basicUnsafeSlice i n v + basicUnsafeIndexM (V_KBNSum v) i = uncurry KBNSum `liftM` G.basicUnsafeIndexM v i + basicUnsafeCopy (MV_KBNSum mv) (V_KBNSum v) = G.basicUnsafeCopy mv v + elemseq _ = seq + + +instance U.Unbox KBNSum + +instance Summation KBNSum where + zero = KBNSum 0 0 + add = kbnAdd + +instance NFData KBNSum where + rnf !_ = () + +kbnAdd :: KBNSum -> Double -> KBNSum +kbnAdd (KBNSum sum c) x = KBNSum sum' c' + where c' | abs sum >= abs x = c + ((sum - sum') + x) + | otherwise = c + ((x - sum') + sum) + sum' = sum + x + +-- | Return the result of a Kahan-Babuška-Neumaier sum. +kbn :: KBNSum -> Double +kbn (KBNSum sum c) = sum + c + +-- | /O(n)/ Sum a vector of values. +sumVector :: (Vector v Double, Summation s) => + (s -> Double) -> v Double -> Double +sumVector f = f . foldl' add zero +{-# INLINE sumVector #-} + +-- $usage +-- +-- Most of these summation algorithms are intended to be used via the +-- 'Summation' typeclass interface. Explicit type annotations should +-- not be necessary, as the use of a function such as 'kbn' or 'kb2' +-- to extract the final sum out of a 'Summation' instance gives the +-- compiler enough information to determine the precise type of +-- summation algorithm to use. +-- +-- As an example, here is a (somewhat silly) function that manually +-- computes the sum of elements in a list. +-- +-- @ +-- sillySumList :: [Double] -> Double +-- sillySumList = loop 'zero' +-- where loop s [] = 'kbn' s +-- loop s (x:xs) = 'seq' s' loop s' xs +-- where s' = 'add' s x +-- @ +-- +-- In most instances, you can simply use the much more general 'sum' +-- function instead of writing a summation function by hand. +-- +-- @ +-- -- Avoid ambiguity around which sum function we are using. +-- import Prelude hiding (sum) +-- -- +-- betterSumList :: [Double] -> Double +-- betterSumList xs = 'sum' 'kbn' xs +-- @ + +-- Note well the use of 'seq' in the example above to force the +-- evaluation of intermediate values. If you must write a summation +-- function by hand, and you forget to evaluate the intermediate +-- values, you are likely to incur a space leak. +-- +-- Here is an example of how to compute a prefix sum in which the +-- intermediate values are as accurate as possible. +-- +-- @ +-- prefixSum :: [Double] -> [Double] +-- prefixSum xs = map 'kbn' . 'scanl' 'add' 'zero' $ xs +-- @ + +-- $references +-- +-- * Kahan, W. (1965), Further remarks on reducing truncation +-- errors. /Communications of the ACM/ 8(1):40. +-- +-- * Neumaier, A. (1974), Rundungsfehleranalyse einiger Verfahren zur +-- Summation endlicher Summen. +-- /Zeitschrift für Angewandte Mathematik und Mechanik/ 54:39–51. +-- +-- * Klein, A. (2006), A Generalized +-- Kahan-Babuška-Summation-Algorithm. /Computing/ 76(3):279-293. +-- +-- * Higham, N.J. (1993), The accuracy of floating point +-- summation. /SIAM Journal on Scientific Computing/ 14(4):783–799. diff --git a/mwc-random/System/Random/MWC.hs b/mwc-random/System/Random/MWC.hs new file mode 100644 --- /dev/null +++ b/mwc-random/System/Random/MWC.hs @@ -0,0 +1,428 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts, + MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, + ForeignFunctionInterface #-} +-- | +-- Module : System.Random.MWC +-- Copyright : (c) 2009-2012 Bryan O'Sullivan +-- License : BSD3 +-- +-- Maintainer : bos@serpentine.com +-- Stability : experimental +-- Portability : portable +-- +-- Pseudo-random number generation. This module contains code for +-- generating high quality random numbers that follow a uniform +-- distribution. +-- +-- For non-uniform distributions, see the +-- 'System.Random.MWC.Distributions' module. +-- +-- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) +-- multiply-with-carry generator, which has a period of 2^8222 and +-- fares well in tests of randomness. It is also extremely fast, +-- between 2 and 3 times faster than the Mersenne Twister. +-- +-- The generator state is stored in the 'Gen' data type. It can be +-- created in several ways: +-- +-- 1. Using the 'withSystemRandom' call, which creates a random state. +-- +-- 2. Supply your own seed to 'initialize' function. +-- +-- 3. Finally, 'create' makes a generator from a fixed seed. +-- Generators created in this way aren't really random. +-- +-- For repeatability, the state of the generator can be snapshotted +-- and replayed using the 'save' and 'restore' functions. +-- +-- The simplest use is to generate a vector of uniformly distributed values: +-- +-- @ +-- vs \<- 'withSystemRandom' . 'asGenST' $ \\gen -> 'uniformVector' gen 100 +-- @ +-- +-- These values can be of any type which is an instance of the class +-- 'Variate'. +-- +-- To generate random values on demand, first 'create' a random number +-- generator. +-- +-- @ +-- gen <- 'create' +-- @ +-- +-- Hold onto this generator and use it wherever random values are +-- required (creating a new generator is expensive compared to +-- generating a random number, so you don't want to throw them +-- away). Get a random value using 'uniform' or 'uniformR': +-- +-- @ +-- v <- 'uniform' gen +-- @ +-- +-- @ +-- v <- 'uniformR' (1, 52) gen +-- @ +module System.Random.MWC + ( + -- * Gen: Pseudo-Random Number Generators + Gen + , initialize + , createSystemRandom + , GenIO + , splitGen + + -- * Variates: uniformly distributed values + , Variate(..) + , uniformVector + + ) where + +#if defined(__GLASGOW_HASKELL__) && !defined(__HADDOCK__) +#include "MachDeps.h" +#endif + +import Control.Monad (liftM, replicateM) +import Data.Bits ((.|.), shiftL, shiftR) +import Data.Int (Int8, Int16, Int32, Int64) +import Data.Vector.Generic (Vector) +import Data.Word (Word8, Word16, Word32, Word64) +#if !MIN_VERSION_base(4,8,0) +import Data.Word (Word) +#endif +import Foreign.Marshal.Alloc (allocaBytes) +import Foreign.Marshal.Array (peekArray) +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Unboxed as I +import qualified Data.Vector.Unboxed.Mutable as M +#if defined(mingw32_HOST_OS) +import Foreign.Ptr +import Foreign.C.Types +#else +import System.IO (IOMode(..), hGetBuf, withBinaryFile) +#endif + +import Basement.Monad (PrimMonad(..)) + + +-- | The class of types for which we can generate uniformly +-- distributed random variates. +-- +-- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) +-- multiply-with-carry generator, which has a period of 2^8222 and +-- fares well in tests of randomness. It is also extremely fast, +-- between 2 and 3 times faster than the Mersenne Twister. +-- +-- /Note/: Marsaglia's PRNG is not known to be cryptographically +-- secure, so you should not use it for cryptographic operations. +class Variate a where + -- | Generate a single uniformly distributed random variate. The + -- range of values produced varies by type: + -- + -- * For fixed-width integral types, the type's entire range is + -- used. + -- + -- * For floating point numbers, the range (0,1] is used. Zero is + -- explicitly excluded, to allow variates to be used in + -- statistical calculations that require non-zero values + -- (e.g. uses of the 'log' function). + -- + -- To generate a 'Float' variate with a range of [0,1), subtract + -- 2**(-33). To do the same with 'Double' variates, subtract + -- 2**(-53). + uniform :: Gen -> IO a + -- | Generate single uniformly distributed random variable in a + -- given range. + -- + -- * For integral types inclusive range is used. + -- + -- * For floating point numbers range (a,b] is used if one ignores + -- rounding errors. + uniformR :: (a,a) -> Gen -> IO a + +instance Variate Word32 where + uniform = uniform1 fromIntegral + uniformR a b = uniformRange a b + {-# INLINE uniform #-} + {-# INLINE uniformR #-} + +instance Variate Word64 where + uniform = uniform2 wordsTo64Bit + uniformR a b = uniformRange a b + {-# INLINE uniform #-} + {-# INLINE uniformR #-} + +instance Variate Int where +#if WORD_SIZE_IN_BITS < 64 + uniform = uniform1 fromIntegral +#else + uniform = uniform2 wordsTo64Bit +#endif + uniformR a b = uniformRange a b + {-# INLINE uniform #-} + {-# INLINE uniformR #-} + +instance Variate Word where +#if WORD_SIZE_IN_BITS < 64 + uniform = uniform1 fromIntegral +#else + uniform = uniform2 wordsTo64Bit +#endif + uniformR a b = uniformRange a b + {-# INLINE uniform #-} + {-# INLINE uniformR #-} + +wordsTo64Bit :: (Integral a) => Word32 -> Word32 -> a +wordsTo64Bit x y = + fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64) +{-# INLINE wordsTo64Bit #-} + +-- | State of the pseudo-random number generator. It uses mutable +-- state so same generator shouldn't be used from the different +-- threads simultaneously. +newtype Gen = Gen (M.MVector (PrimState IO) Word32) + +-- | A shorter name for PRNG state in the 'IO' monad. +type GenIO = Gen + +ioff, coff :: Int +ioff = 256 +coff = 257 + +-- | Create a generator for variates using the given seed of 256 elements +-- +-- @gen' <- 'initialize' . 'fromSeed' =<< 'save'@ +initialize :: I.Vector Word32 -> IO Gen +initialize seed + | fini /= 256 = error "mwc seed invalid size" + | otherwise = do + q <- M.unsafeNew 258 + fill q + M.unsafeWrite q ioff 255 + M.unsafeWrite q coff 362436 + return (Gen q) + where + fini = G.length seed + fill q = go 0 where + go i | i == 256 = return () + | otherwise = M.unsafeWrite q i (G.unsafeIndex seed i) >> go (i+1) +{-# INLINE initialize #-} + +-- | Acquire seed from the system entropy source. On Unix machines, +-- this will attempt to use @/dev/urandom@. On Windows, it will internally +-- use @RtlGenRandom@. +acquireSeedSystem :: IO [Word32] +acquireSeedSystem = do +#if !defined(mingw32_HOST_OS) + -- Read 256 random Word32s from /dev/urandom + let nbytes = 1024 + random = "/dev/urandom" + allocaBytes nbytes $ \buf -> do + nread <- withBinaryFile random ReadMode $ + \h -> hGetBuf h buf nbytes + peekArray (nread `div` 4) buf +#else + let nbytes = 1024 + -- Generate 256 random Word32s from RtlGenRandom + allocaBytes nbytes $ \buf -> do + ok <- c_RtlGenRandom buf (fromIntegral nbytes) + if ok then return () else fail "Couldn't use RtlGenRandom" + peekArray (nbytes `div` 4) buf + +-- Note: on 64-bit Windows, the 'stdcall' calling convention +-- isn't supported, so we use 'ccall' instead. +#if defined(i386_HOST_ARCH) +# define WINDOWS_CCONV stdcall +#elif defined(x86_64_HOST_ARCH) +# define WINDOWS_CCONV ccall +#else +# error Unknown mingw32 architecture! +#endif + +-- Note: On Windows, the typical convention would be to use +-- the CryptoGenRandom API in order to generate random data. +-- However, here we use 'SystemFunction036', AKA RtlGenRandom. +-- +-- This is a commonly used API for this purpose; one bonus is +-- that it avoids having to bring in the CryptoAPI library, +-- and completely sidesteps the initialization cost of CryptoAPI. +-- +-- While this function is technically "subject to change" that is +-- extremely unlikely in practice: rand_s in the Microsoft CRT uses +-- this, and they can't change it easily without also breaking +-- backwards compatibility with e.g. statically linked applications. +-- +-- The name 'SystemFunction036' is the actual link-time name; the +-- display name is just for giggles, I guess. +-- +-- See also: +-- - http://blogs.msdn.com/b/michael_howard/archive/2005/01/14/353379.aspx +-- - https://bugzilla.mozilla.org/show_bug.cgi?id=504270 +-- +foreign import WINDOWS_CCONV unsafe "SystemFunction036" + c_RtlGenRandom :: Ptr a -> CULong -> IO Bool +#endif + +-- | Seed a PRNG with data from the system's fast source of pseudo-random +-- numbers. +createSystemRandom :: IO GenIO +createSystemRandom = do + seed <- acquireSeedSystem + initialize (I.fromList seed) + +-- | Compute the next index into the state pool. This is simply +-- addition modulo 256. +nextIndex :: Integral a => a -> Int +nextIndex i = fromIntegral j + where j = fromIntegral (i+1) :: Word8 +{-# INLINE nextIndex #-} + +aa :: Word64 +aa = 1540315826 +{-# INLINE aa #-} + +data DoubleWord32 = DoubleWord32 {-# UNPACK #-} !Word32 {-# UNPACK #-} !Word32 + +uniformWord32 :: Gen -> IO Word32 +uniformWord32 (Gen q) = do + i <- nextIndex `liftM` M.unsafeRead q ioff + c <- fromIntegral `liftM` M.unsafeRead q coff + qi <- fromIntegral `liftM` M.unsafeRead q i + let t = aa * qi + c + c' = fromIntegral (t `shiftR` 32) + x = fromIntegral t + c' + (DoubleWord32 x' c'') | x < c' = DoubleWord32 (x + 1) (c' + 1) + | otherwise = DoubleWord32 x c' + M.unsafeWrite q i x' + M.unsafeWrite q ioff (fromIntegral i) + M.unsafeWrite q coff (fromIntegral c'') + return x' +{-# INLINE uniformWord32 #-} + +uniform1 :: (Word32 -> a) -> Gen -> IO a +uniform1 f gen = do + i <- uniformWord32 gen + return $! f i +{-# INLINE uniform1 #-} + +uniform2 :: (Word32 -> Word32 -> a) -> Gen -> IO a +uniform2 f (Gen q) = do + i <- nextIndex `liftM` M.unsafeRead q ioff + let j = nextIndex i + c <- fromIntegral `liftM` M.unsafeRead q coff + qi <- fromIntegral `liftM` M.unsafeRead q i + qj <- fromIntegral `liftM` M.unsafeRead q j + let t = aa * qi + c + c' = fromIntegral (t `shiftR` 32) + x = fromIntegral t + c' + DoubleWord32 x' c'' | x < c' = DoubleWord32 (x + 1) (c' + 1) + | otherwise = DoubleWord32 x c' + u = aa * qj + fromIntegral c'' + d' = fromIntegral (u `shiftR` 32) + y = fromIntegral u + d' + DoubleWord32 y' d'' | y < d' = DoubleWord32 (y + 1) (d' + 1) + | otherwise = DoubleWord32 y d' + M.unsafeWrite q i x' + M.unsafeWrite q j y' + M.unsafeWrite q ioff (fromIntegral j) + M.unsafeWrite q coff (fromIntegral d'') + return $! f x' y' +{-# INLINE uniform2 #-} + +-- Type family for fixed size integrals. For signed data types it's +-- its unsigned couterpart with same size and for unsigned data types +-- it's same type +type family Unsigned a :: * + +type instance Unsigned Int8 = Word8 +type instance Unsigned Int16 = Word16 +type instance Unsigned Int32 = Word32 +type instance Unsigned Int64 = Word64 + +type instance Unsigned Word8 = Word8 +type instance Unsigned Word16 = Word16 +type instance Unsigned Word32 = Word32 +type instance Unsigned Word64 = Word64 + +-- This is workaround for bug #25. +-- +-- GHC-7.6 has a bug (#8072) which results in calculation of wrong +-- number of buckets in function `uniformRange'. Consequently uniformR +-- generates values in wrong range. +-- +-- Bug only affects 32-bit systems and Int/Word data types. Word32 +-- works just fine. So we set Word32 as unsigned counterpart for Int +-- and Word on 32-bit systems. It's done only for GHC-7.6 because +-- other versions are unaffected by the bug and we expect that GHC may +-- optimise code which uses Word better. +#if (WORD_SIZE_IN_BITS < 64) && (__GLASGOW_HASKELL__ == 706) +type instance Unsigned Int = Word32 +type instance Unsigned Word = Word32 +#else +type instance Unsigned Int = Word +type instance Unsigned Word = Word +#endif + + +-- Subtract two numbers under assumption that x>=y and store result in +-- unsigned data type of same size +sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a +sub x y = fromIntegral x - fromIntegral y +{-# INLINE sub #-} + +add :: (Integral a, Integral (Unsigned a)) => a -> Unsigned a -> a +add m x = m + fromIntegral x +{-# INLINE add #-} + +-- Generate uniformly distributed value in inclusive range. +-- +-- NOTE: This function must be fully applied. Otherwise it won't be +-- inlined, which will cause a severe performance loss. +-- +-- > uniformR = uniformRange -- won't be inlined +-- > uniformR a b = uniformRange a b -- will be inlined +uniformRange :: ( Integral a, Bounded a, Variate a + , Integral (Unsigned a), Bounded (Unsigned a), Variate (Unsigned a)) + => (a,a) -> Gen -> IO a +uniformRange (x1,x2) g + | n == 0 = uniform g -- Abuse overflow in unsigned types + | otherwise = loop + where + -- Allow ranges where x2 Gen -> Int -> IO (v a) +uniformVector gen n = G.replicateM n (uniform gen) +{-# INLINE uniformVector #-} + +-- | 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 + +-- $references +-- +-- * Marsaglia, G. (2003) Seeds for random number generators. +-- /Communications of the ACM/ 46(5):90–93. +-- +-- +-- * Doornik, J.A. (2007) Conversion of high-period random numbers to +-- floating point. +-- /ACM Transactions on Modeling and Computer Simulation/ 17(1). +-- diff --git a/statistics/Statistics/Regression.hs b/statistics/Statistics/Regression.hs --- a/statistics/Statistics/Regression.hs +++ b/statistics/Statistics/Regression.hs @@ -20,11 +20,10 @@ 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 System.Random.MWC (GenIO, uniformR, splitGen) import qualified Data.Vector as V import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U diff --git a/statistics/Statistics/Resampling.hs b/statistics/Statistics/Resampling.hs --- a/statistics/Statistics/Resampling.hs +++ b/statistics/Statistics/Resampling.hs @@ -23,16 +23,12 @@ , 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 @@ -45,7 +41,7 @@ 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) +import System.Random.MWC (Gen, GenIO, uniformR, splitGen) ---------------------------------------------------------------- @@ -131,8 +127,7 @@ ests' = map estimate ests -- | Create vector using resamples -resampleVector :: G.Vector v a - => Gen (PrimState IO) -> v a -> IO (v a) +resampleVector :: G.Vector v a => Gen -> v a -> IO (v a) resampleVector gen v = G.replicateM n $ do i <- uniformR (0,n-1) gen return $! G.unsafeIndex v i @@ -209,11 +204,3 @@ 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))