Commits

Bryan O'Sullivan committed a831627

Add splitGen

Comments (0)

Files changed (2)

Statistics/Regression.hs

 import Control.Concurrent.Chan (newChan, readChan, writeChan)
 import Control.DeepSeq (rnf)
 import Control.Monad (forM_, replicateM)
-import Data.Word (Word32)
 import GHC.Conc (getNumCapabilities)
 import Prelude hiding (pred, sum)
 import Statistics.Function as F
 import Statistics.Matrix hiding (map)
 import Statistics.Matrix.Algorithms (qr)
+import Statistics.Resampling (splitGen)
 import Statistics.Resampling.Bootstrap (Estimate(..))
 import Statistics.Sample (mean)
 import Statistics.Sample.Internal (sum)
-import System.Random.MWC (GenIO, initialize, uniformR, uniformVector)
+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
                  -> IO (V.Vector Estimate, Estimate)
 bootstrapRegress gen0 numResamples ci rgrss preds0 resp0 = do
   caps <- getNumCapabilities
-  gens <- fmap (gen0:) . replicateM (caps-1) $
-          initialize =<< (uniformVector gen0 256 :: IO (U.Vector Word32))
+  gens <- splitGen caps gen0
   done <- newChan
   forM_ (zip gens (balance caps numResamples)) $ \(gen,count) -> do
     forkIO $ do

Statistics/Resampling.hs

     , jackknifeStdDev
     , resample
     , estimate
+    , splitGen
     ) where
 
 import Data.Aeson (FromJSON, ToJSON)
 import Control.Concurrent (forkIO, newChan, readChan, writeChan)
-import Control.Monad (forM_, liftM, replicateM_)
-import Control.Monad.Primitive (PrimState)
+import Control.Monad (forM_, liftM, replicateM, replicateM_)
 import Data.Binary (Binary(..))
 import Data.Data (Data, Typeable)
 import Data.Vector.Algorithms.Intro (sort)
 import Statistics.Function (indices)
 import Statistics.Sample (mean, stdDev, variance, varianceUnbiased)
 import Statistics.Types (Estimator(..), Sample)
-import System.Random.MWC (Gen, initialize, uniform, uniformVector)
+import System.Random.MWC (GenIO, initialize, uniform, uniformVector)
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
 -- available CPUs.  At least with GHC 7.0, parallel performance seems
 -- best if the parallel garbage collector is disabled (RTS option
 -- @-qg@).
-resample :: Gen (PrimState IO)
+resample :: GenIO
          -> [Estimator]         -- ^ Estimation functions.
          -> Int                 -- ^ Number of resamples to compute.
          -> Sample              -- ^ Original sample.
 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))