randomVector

3 files changed, 42 insertions, 2 deletions

diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 21f51e5..836f11a 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -21,7 +21,8 @@ module Data.Packed.Vector (
     vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,
     mapVector, zipVector,
     fscanfVector, fprintfVector, freadVector, fwriteVector,
-    foldLoop, foldVector, foldVectorG
+    foldLoop, foldVector, foldVectorG,
+    RandDist(..), randomVector
 ) where
 
 import Data.Packed.Internal
diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs
index 0366744..2b21de2 100644
--- a/lib/Numeric/GSL/Vector.hs
+++ b/lib/Numeric/GSL/Vector.hs
@@ -17,7 +17,8 @@ module Numeric.GSL.Vector (
     FunCodeS(..), toScalarR,
     FunCodeV(..), vectorMapR, vectorMapC,
     FunCodeSV(..), vectorMapValR, vectorMapValC,
-    FunCodeVV(..), vectorZipR, vectorZipC
+    FunCodeVV(..), vectorZipR, vectorZipC,
+    RandDist(..), randomVector
 ) where
 
 import Data.Packed.Internal.Common
@@ -146,3 +147,21 @@ vectorZipC :: FunCodeVV -> Vector (Complex Double) -> Vector (Complex Double) ->
 vectorZipC = vectorZipAux c_vectorZipC
 
 foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV
+
+-----------------------------------------------------------------------
+
+data RandDist = Uniform  -- ^ uniform distribution in [0,1)
+              | Gaussian -- ^ normal distribution with mean zero and standard deviation one
+              deriving Enum
+
+-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
+randomVector :: Int      -- ^ seed
+             -> RandDist -- ^ distribution
+             -> Int      -- ^ vector size
+             -> Vector Double
+randomVector seed dist n = unsafePerformIO $ do
+    r <- createVector n
+    app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
+    return r
+
+foreign import ccall safe "random_vector" c_random_vector :: CInt -> CInt -> TV
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index d129aeb..0c71ca1 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -20,6 +20,8 @@
 #include <gsl/gsl_multimin.h>
 #include <gsl/gsl_multiroots.h>
 #include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -777,3 +779,21 @@ int rootj(int method, int f(int, double*, int, double*),
     gsl_multiroot_fdfsolver_free(s);
     OK
 }
+
+//////////////////////////////////////////////////////
+
+#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
+
+int random_vector(int seed, int code, RVEC(r)) {
+    DEBUGMSG("random_vector")
+    static gsl_rng * gen = NULL;
+    if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
+    gsl_rng_set (gen, seed);
+    int k;
+    switch (code) {
+        RAN(0,gsl_rng_uniform)
+        RAN(1,gsl_ran_ugaussian)
+        default: ERROR(BAD_CODE);
+    }
+}
+#undef RAN