diff --git a/GSL/Random/Dist.hs b/GSL/Random/Dist.hs new file mode 100644 --- /dev/null +++ b/GSL/Random/Dist.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +-- | +-- Module : GSL.Random.Dist +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module GSL.Random.Dist ( + getPoisson, + getFlat, + getGaussian, + ) where + +import Foreign.C.Types ( CUInt, CDouble ) +import Foreign.ForeignPtr ( withForeignPtr ) +import Foreign.Ptr ( Ptr ) + +import GSL.Random.Gen.Internal ( RNG(..) ) + +-- | @getPoisson r mu@ gets a poisson random variable with mean @mu@. +getPoisson :: RNG -> Double -> IO Int +getPoisson (MkRNG fptr) mu = + let mu' = realToFrac mu + in withForeignPtr fptr $ \ptr -> do + x <- gsl_ran_poisson ptr mu' + return $ (fromInteger . toInteger) x + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_poisson :: Ptr () -> CDouble -> IO CUInt + +-- | @getFlat r a b@ gets a value uniformly chosen in @[a,b)@. +getFlat :: RNG -> Double -> Double -> IO (Double) +getFlat (MkRNG fptr) a b = + let a' = realToFrac a + b' = realToFrac b + in withForeignPtr fptr $ \ptr -> do + x <- gsl_ran_flat ptr a' b' + return $ realToFrac x + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_flat :: Ptr () -> CDouble -> CDouble -> IO CDouble + + +-- | @getGaussian r sigma@ gets a normal random variable with mean +-- @0@ and standard deviation @sigma@. +getGaussian :: RNG -> Double -> IO Double +getGaussian (MkRNG fptr) sigma = + let sigma' = realToFrac sigma + in withForeignPtr fptr $ \ptr -> do + x <- gsl_ran_gaussian ptr sigma' + return $ realToFrac x + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_ran_gaussian :: Ptr () -> CDouble -> IO CDouble + diff --git a/GSL/Random/Gen.hs b/GSL/Random/Gen.hs new file mode 100644 --- /dev/null +++ b/GSL/Random/Gen.hs @@ -0,0 +1,14 @@ +----------------------------------------------------------------------------- +-- | +-- Module : GSL.Random.Gen +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module GSL.Random.Gen ( + module GSL.Random.Gen.Internal + ) where + +import GSL.Random.Gen.Internal hiding ( MkRNG ) diff --git a/GSL/Random/Gen/Internal.hs b/GSL/Random/Gen/Internal.hs new file mode 100644 --- /dev/null +++ b/GSL/Random/Gen/Internal.hs @@ -0,0 +1,245 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +----------------------------------------------------------------------------- +-- | +-- Module : GSL.Random.Gen.Internal +-- Copyright : Copyright (c) , Patrick Perry +-- License : BSD3 +-- Maintainer : Patrick Perry +-- Stability : experimental +-- + +module GSL.Random.Gen.Internal ( + -- * Data types + RNG(..), + RNGType, + + -- * Initializing + newRNG, + setSeed, + + -- * Sampling + getSample, + getUniform, + getUniformPos, + getUniformInt, + + -- * Auxiliary functions + getName, + getMax, + getMin, + getSize, + getState, + setState, + + -- * Copying state + copyRNG, + cloneRNG, + + -- * Algorithms + mt19937, + + ) where + +import Control.Monad ( liftM ) +import Data.Maybe ( fromJust ) +import Data.Word ( Word8, Word64 ) +import Foreign.C.Types ( CULong, CSize, CInt, CDouble ) +import Foreign.C.String ( CString, withCAString, peekCAString ) +import Foreign.ForeignPtr ( ForeignPtr, newForeignPtr, withForeignPtr ) +import Foreign.Marshal.Array ( peekArray, pokeArray, advancePtr ) +import Foreign.Ptr ( Ptr, FunPtr, castPtr, nullPtr ) +import Foreign.Storable ( peek ) +import System.IO.Unsafe ( unsafePerformIO ) + +newtype RNG = MkRNG (ForeignPtr ()) +newtype RNGType = MkRNGType (Ptr ()) + + +--------------------------- Initialization ---------------------------------- + +-- | Allocate a new random number generator of the given type and initialize +-- it with the default seed. +newRNG :: RNGType -> IO RNG +newRNG t = do + ptr <- gsl_rng_alloc t + fptr <- newForeignPtr p_gsl_rng_free ptr + return $! MkRNG fptr + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_alloc :: RNGType -> IO (Ptr ()) + +foreign import ccall unsafe "gsl/gsl_rng.h &gsl_rng_free" + p_gsl_rng_free :: FunPtr (Ptr () -> IO ()) + +-- | Seed the generator with the given value. +setSeed :: RNG -> Word64 -> IO () +setSeed (MkRNG fptr) seed = + let seed' = (fromInteger . toInteger) seed + in withForeignPtr fptr $ flip gsl_rng_set seed' + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_set :: Ptr () -> CULong -> IO () + + +--------------------------- Sampling ---------------------------------------- + +-- | Returns a value uniform in [rngMin, rngMax] +getSample :: RNG -> IO Word64 +getSample (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + gsl_rng_get ptr >>= return . fromInteger . toInteger + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_get :: Ptr () -> IO CULong + +-- | Returns a value uniform on [0,1) +getUniform :: RNG -> IO Double +getUniform (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + gsl_rng_uniform ptr >>= return . realToFrac + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_uniform :: Ptr () -> IO CDouble + +-- | Returns a value uniform on (0,1) +getUniformPos :: RNG -> IO Double +getUniformPos (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + gsl_rng_uniform_pos ptr >>= return . realToFrac + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_uniform_pos :: Ptr () -> IO CDouble + +-- | Returns an integer uniform on [0,n-1]. @n@ must be greater than @0@. +getUniformInt :: RNG -> Int -> IO Int +getUniformInt (MkRNG fptr) n + | n <= 0 = + ioError $ userError $ + "rngUnifInt: expected \"n\" to be greater than 0" ++ + " but got `" ++ show n ++ "' instead." + | otherwise = + let n' = (fromInteger . toInteger) n + in withForeignPtr fptr $ \ptr -> + gsl_rng_uniform_int ptr n' >>= return . fromInteger . toInteger + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_uniform_int :: Ptr () -> CULong -> IO CULong + + +--------------------------- Auxiliary Functions ----------------------------- + +-- | Get the name of the generator. +getName :: RNG -> IO String +getName (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + peekCAString (gsl_rng_name ptr) + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_name :: Ptr () -> CString + +-- | Get the largest value that the generator can return. +getMax :: RNG -> IO Word64 +getMax (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + (return . fromInteger . toInteger) (gsl_rng_max ptr) + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_max :: Ptr () -> CULong + +-- | Get the smallest value that the generator can return. +getMin :: RNG -> IO Word64 +getMin (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + (return . fromInteger . toInteger) (gsl_rng_min ptr) + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_min :: Ptr () -> CULong + +-- | Get the size of the generator state, in bytes. +getSize :: RNG -> IO Word64 +getSize (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> + (return . fromInteger . toInteger) (gsl_rng_size ptr) + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_size :: Ptr () -> CSize + +-- | Get the generator state. +getState :: RNG -> IO [Word8] +getState rng@(MkRNG fptr) = do + n <- liftM (fromInteger . toInteger) (getSize rng) + withForeignPtr fptr $ \ptr -> + peekArray n (gsl_rng_state ptr) + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_rng_state :: Ptr () -> Ptr Word8 + +-- | Set the generator state. The input array should have size equal +-- to @getSize@ of the generator; otherwise, strange things will happen. +setState :: RNG -> [Word8] -> IO () +setState (MkRNG fptr) state = do + withForeignPtr fptr $ \ptr -> + pokeArray (gsl_rng_state ptr) state + +--------------------------- Copying State ----------------------------------- + +-- | @copyRNG dst src@ copies the state from one generator to another. The +-- two generators must have the same type. +copyRNG :: RNG -> RNG -> IO () +copyRNG (MkRNG fdst) (MkRNG fsrc) = + withForeignPtr fdst $ \dst -> + withForeignPtr fsrc $ \src -> + gsl_rng_memcpy dst src + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_memcpy :: Ptr () -> Ptr () -> IO () + +-- | Allocate a new random number generator that is an exact copy of another +-- generator +cloneRNG :: RNG -> IO RNG +cloneRNG (MkRNG fptr) = + withForeignPtr fptr $ \ptr -> do + ptr' <- gsl_rng_clone ptr + fptr' <- newForeignPtr p_gsl_rng_free ptr' + ptr' `seq` return $! MkRNG fptr' + +foreign import ccall unsafe "gsl/gsl_rng.h" + gsl_rng_clone :: Ptr () -> IO (Ptr ()) + + +--------------------------- Algorithms -------------------------------------- + + +mt19937 :: RNGType +mt19937 = + unsafePerformIO $ + withCAString "mt19937" $ \name -> + liftM fromJust $ getRngType name +{-# NOINLINE mt19937 #-} + +getRngType :: CString -> IO (Maybe RNGType) +getRngType name = + go gsl_rng_types_setup + where + go :: Ptr (Ptr ()) -> IO (Maybe RNGType) + go types = do + t <- peek types + if t == nullPtr + then + return Nothing + + else do + name' <- peek (castPtr t) + cmp <- c_strcmp name name' + + if cmp == 0 + then return $ Just (MkRNGType t) + else go $ types `advancePtr` 1 + + + +foreign import ccall unsafe "string.h strcmp" + c_strcmp :: CString -> CString -> IO CInt + +foreign import ccall unsafe "gsl/gsl_randist.h" + gsl_rng_types_setup :: Ptr (Ptr ()) diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) Patrick Perry 2008 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions +are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF +SUCH DAMAGE. diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,3 @@ +You may need to edit the cabal file to get the right libraries linked in. In +particular, you may need to change "cblas" to "gslcblas". Or, if you are using +ATLAS, you may need to add atlas at the end of the library list. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,6 @@ +#!/usr/bin/env runhaskell +> import Distribution.Simple +> import System.Cmd +> +> main = defaultMain +> diff --git a/gsl-random.cabal b/gsl-random.cabal new file mode 100644 --- /dev/null +++ b/gsl-random.cabal @@ -0,0 +1,32 @@ +name: gsl-random +version: 0.1 +homepage: http://stat.stanford.edu/~patperry/code/gsl-random +synopsis: Bindings the the GSL random number generation facilities. +description: + Bindings to the GNU Scientific Library random number generators and random + distributions. + . +category: Math +license: BSD3 +license-file: LICENSE +copyright: (c) 2008. Patrick Perry +author: Patrick Perry +maintainer: Patrick Perry +cabal-version: >= 1.2.0 +build-type: Simple +tested-with: GHC ==6.8.3 + +extra-source-files: README + +library + exposed-modules: GSL.Random.Gen + GSL.Random.Dist + other-modules: GSL.Random.Gen.Internal + + ghc-options: -Wall + extensions: ForeignFunctionInterface + build-depends: base + + extra-libraries: gsl cblas + extra-lib-dirs: /opt/local/lib +