diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -10,3 +10,11 @@ 0.1.1.3: minor improvements + +0.1.1.4: + added emptyRanges, emptyLimits to Histogram(2D) + fixed lack of random0_dist + changed random_mm_dist to random_mp_dist + changed discrete_mm_dist to discrete_mp_dist + added Binary instances of Histogram(2D) + added Histogram2D.fromMatrix diff --git a/hmatrix-gsl-stats.cabal b/hmatrix-gsl-stats.cabal --- a/hmatrix-gsl-stats.cabal +++ b/hmatrix-gsl-stats.cabal @@ -1,5 +1,5 @@ Name: hmatrix-gsl-stats -Version: 0.1.1.3 +Version: 0.1.1.4 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010 @@ -25,8 +25,8 @@ library - Build-Depends: base >= 3 && < 5, - storable-complex, hmatrix >= 0.9.3 + Build-Depends: base >= 3 && < 5, ghc-binary, + storable-complex, hmatrix >= 0.10.0 Extensions: ForeignFunctionInterface diff --git a/lib/Numeric/GSL/Distribution/Continuous.hs b/lib/Numeric/GSL/Distribution/Continuous.hs --- a/lib/Numeric/GSL/Distribution/Continuous.hs +++ b/lib/Numeric/GSL/Distribution/Continuous.hs @@ -244,7 +244,7 @@ return r' foreign import ccall "distribution-aux.h random_mp" distribution_random_multi_param :: CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt -foreign import ccall "distribution-aux.h random_mm_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt +foreign import ccall "distribution-aux.h random_mp_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt ----------------------------------------------------------------------------- @@ -292,7 +292,7 @@ app1 (distribution_spherical_vector (fromIntegral s)) vec r "spherical_vector" return r -foreign import ccall "distribution-aux.h random_vector" distribution_spherical_vector :: CInt -> CInt -> Ptr Double -> IO CInt +foreign import ccall "distribution-aux.h spherical_vector" distribution_spherical_vector :: CInt -> CInt -> Ptr Double -> IO CInt ----------------------------------------------------------------------------- ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Distribution/Discrete.hs b/lib/Numeric/GSL/Distribution/Discrete.hs --- a/lib/Numeric/GSL/Distribution/Discrete.hs +++ b/lib/Numeric/GSL/Distribution/Discrete.hs @@ -203,7 +203,7 @@ return r' foreign import ccall "distribution-aux.h discrete_mp" distribution_discrete_multi_param :: CInt -> CInt -> CUInt -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt -foreign import ccall "distribution-aux.h discrete_mm_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt +foreign import ccall "distribution-aux.h discrete_mp_dist" distribution_dist_multi_param :: CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr CUInt -> IO CInt ----------------------------------------------------------------------------- ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Histogram.hs b/lib/Numeric/GSL/Histogram.hs --- a/lib/Numeric/GSL/Histogram.hs +++ b/lib/Numeric/GSL/Histogram.hs @@ -15,6 +15,7 @@ module Numeric.GSL.Histogram ( Histogram + , emptyRanges, emptyLimits , fromRanges, fromLimits , addList, addVector, addListWeighted, addVectorWeighted , toVectors @@ -44,6 +45,7 @@ import Numeric.LinearAlgebra.Interface() --import Control.Monad +import Data.Binary import Foreign hiding(shift) --import Foreign.ForeignPtr @@ -93,6 +95,27 @@ -} ----------------------------------------------------------------------------- +vectorToTuples = toTuples . toList + where toTuples [] = error "need a minimum of two elements" + toTuples [_] = error "need a minimum of two elements" + toTuples [x1,x2] = [(x1,x2)] + toTuples (x1:x2:xs) = (x1,x2) : (toTuples (x2:xs)) + +----------------------------------------------------------------------------- + +instance Binary Histogram where + put h = do + let (b,c) = toVectors h + put b + put c + get = do + b <- get + c <- get + return $! addVectorWeighted (emptyRanges b) (middle b) c + where middle = fromList . map (\(x1,x2) -> (x1 + x2)/2) . vectorToTuples + +----------------------------------------------------------------------------- + foreign import ccall "gsl-histogram.h gsl_histogram_alloc" histogram_new :: CInt -> IO HistHandle foreign import ccall "gsl-histogram.h &gsl_histogram_free" histogram_free :: FunPtr (HistHandle -> IO ()) @@ -119,6 +142,17 @@ foreign import ccall "gsl-histogram.h gsl_histogram_set_ranges_uniform" histogram_set_ranges_uniform :: HistHandle -> Double -> Double -> IO CInt ----------------------------------------------------------------------------- + +-- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) +emptyRanges :: Vector Double -- ^ the ranges + -> Histogram -- ^ result +emptyRanges x = unsafePerformIO $ fromRangesIO x + +-- | create a histogram with n bins and lower and upper limits +emptyLimits :: Int -- ^ bins + -> (Double,Double) -- ^ lower and upper limits + -> Histogram -- ^ result +emptyLimits n x = unsafePerformIO $ fromLimitsIO n x -- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) and increment from a vector fromRanges :: Vector Double -- ^ the ranges diff --git a/lib/Numeric/GSL/Histogram2D.hs b/lib/Numeric/GSL/Histogram2D.hs --- a/lib/Numeric/GSL/Histogram2D.hs +++ b/lib/Numeric/GSL/Histogram2D.hs @@ -15,6 +15,7 @@ module Numeric.GSL.Histogram2D ( Histogram2D + , emptyRanges, emptyLimits , fromRanges, fromLimits , addList, addVector, addListWeighted, addVectorWeighted , toMatrix @@ -46,6 +47,8 @@ --import Control.Monad --import Control.Monad(when) +import Data.Binary + import Foreign hiding(shift) --import Foreign.ForeignPtr --import Foreign.Marshal.Alloc(alloca) @@ -94,6 +97,20 @@ -} ----------------------------------------------------------------------------- +instance Binary Histogram2D where + put h = do + let (rx,ry,w) = toMatrix h + put rx + put ry + put w + get = do + rx <- get + ry <- get + w <- get + return $! fromMatrix rx ry w + +----------------------------------------------------------------------------- + foreign import ccall "gsl-histogram2d.h gsl_histogram2d_alloc" histogram2d_new :: CInt -> CInt -> IO Hist2DHandle foreign import ccall "gsl-histogram2d.h &gsl_histogram2d_free" histogram2d_free :: FunPtr (Hist2DHandle -> IO ()) @@ -125,6 +142,19 @@ ----------------------------------------------------------------------------- +-- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) +emptyRanges :: Vector Double -- ^ the x ranges + -> Vector Double -- ^ the y ranges + -> Histogram2D -- ^ result +emptyRanges x y = unsafePerformIO $ fromRangesIO x y + +-- | create a histogram with n bins and lower and upper limits +emptyLimits :: Int -> Int -- ^ bins + -> (Double,Double) -- ^ lower and upper limits x + -> (Double,Double) -- ^ lower and upper limits y + -> Histogram2D -- ^ result +emptyLimits nx ny x y = unsafePerformIO $ fromLimitsIO nx ny x y + -- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1) and increment from a vector fromRanges :: Vector Double -- ^ the x ranges -> Vector Double -- ^ the y ranges @@ -161,6 +191,27 @@ ----------------------------------------------------------------------------- +vectorToTuples = toTuples . toList + where toTuples [] = error "need a minimum of two elements" + toTuples [_] = error "need a minimum of two elements" + toTuples [x1,x2] = [(x1,x2)] + toTuples (x1:x2:xs) = (x1,x2) : (toTuples (x2:xs)) + +----------------------------------------------------------------------------- + +-- | create from ranges and bins +fromMatrix :: Vector Double -- ^ x ranges + -> Vector Double -- ^ y ranges + -> Matrix Double -- ^ bins + -> Histogram2D -- ^result +fromMatrix x y w = let x' = map (\(x1,x2) -> (x1 + x2)/2) $ vectorToTuples x + y' = map (\(x1,x2) -> (x1 + x2)/2) $ vectorToTuples y + w' = toList $ flatten w + xy = concat $ map (\j -> zip x' (replicate (length x') j)) y' + in addListWeighted (emptyRanges x y) $ zipWith (\(x,y) d -> (x,y,d)) xy w' + +----------------------------------------------------------------------------- + -- | create a copy of a histogram cloneHistogram2D :: Histogram2D -> IO Histogram2D cloneHistogram2D (H nx ny h) = do @@ -404,9 +455,9 @@ -- | returns True of all the individual bin ranges of the two histograms are identical equalBins :: Histogram2D -> Histogram2D -> Bool equalBins (H _ _ h1) (H _ _ h2) = unsafePerformIO $ do - i <- withForeignPtr h1 $ \p1 -> do + j <- withForeignPtr h1 $ \p1 -> do withForeignPtr h2 $ \p2 -> histogram2d_equal_bins p1 p2 - if (fromIntegral i) == (1 :: Int) + if (fromIntegral j) == (1 :: Int) then return True else return False diff --git a/lib/Numeric/GSL/distribution-aux.c b/lib/Numeric/GSL/distribution-aux.c --- a/lib/Numeric/GSL/distribution-aux.c +++ b/lib/Numeric/GSL/distribution-aux.c @@ -2,6 +2,41 @@ #include #include +int random0(int s, int type, double* r) +{ + const gsl_rng_type * T; + gsl_rng * rng; + + gsl_rng_env_setup(); + T = gsl_rng_default; + rng = gsl_rng_alloc(T); + gsl_rng_set(rng,s); + + switch(type) { + case 0: { (*r) = gsl_ran_landau(rng); break; } + } + + gsl_rng_free(rng); + + return 0; +} + +double random0_pdf(int type, double x) +{ + switch (type) { + case 0: return gsl_ran_landau_pdf(x); + } +} + +double random0_dist(int df, int type, double x) +{ + switch(df) { + case 0: return random0_pdf(type,x); + } +} + +////////////////////////////////////////////////////////////////////// + int random1(int s, int type, double par, double* r) { const gsl_rng_type * T; @@ -440,7 +475,7 @@ ////////////////////////////////////////////////////////////////////// -int random_vector(int s, int rs, double* r) +int spherical_vector(int s, int rs, double* r) { const gsl_rng_type * T; gsl_rng * rng;