diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -26,4 +26,7 @@ fixed Histogram(2D)PDF sampling 0.1.2.2: - + removed warning causes in Distribution.Discrete/Continuous + added fromVectors, modified fromMatrix in Histogram(2D) + added countInstance, probability to Histogram(2D) + 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.2.1 +Version: 0.1.2.2 License: GPL License-file: LICENSE Copyright: (c) A.V.H. McPhail 2010 @@ -8,7 +8,13 @@ Stability: provisional Homepage: http://code.haskell.org/hmatrix-gsl-stats Synopsis: GSL Statistics interface -Description: Purely functional interface for statistics based on hmatrix and GSL +Description: + Purely functional interface for statistics based on hmatrix and GSL + . + Relieves the burden of alloc/free routines and is otherwise essentially 1:1 + . + Documentation can be found at + Category: Math, Statistics tested-with: GHC ==6.12.1 diff --git a/lib/Numeric/GSL/Distribution/Common.hs b/lib/Numeric/GSL/Distribution/Common.hs --- a/lib/Numeric/GSL/Distribution/Common.hs +++ b/lib/Numeric/GSL/Distribution/Common.hs @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Common --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL common data types for distributions +-- +-- -- ----------------------------------------------------------------------------- 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 @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Continuous --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -11,6 +11,8 @@ -- -- GSL continuous random distribution functions -- +-- +-- ----------------------------------------------------------------------------- module Numeric.GSL.Distribution.Continuous ( @@ -205,12 +207,10 @@ -> Double -- ^ parameter 3 -> Double -- ^ value -> Double -- ^ result -density_3p d f p1 p2 p3 x = unsafePerformIO $ do +density_3p d _ _ _ _ _ = unsafePerformIO $ do case d of LevySkew -> error "Levy skew has no PDF or CDF" - where density_only f' d' p1' p2' p3' x' = if f' /= Density - then error "distribution has no CDF" - else distribution_dist_three_param (fromei f') (fromei d') x' p1' p2' p3' + _ -> error "unknown 3 parameter distribution" foreign import ccall "distribution-aux.h random3" distribution_random_three_param :: CInt -> CInt -> Double -> Double -> Double -> Ptr Double -> IO CInt foreign import ccall "distribution-aux.h random3_dist" distribution_dist_three_param :: CInt -> CInt -> Double -> Double -> Double -> Double -> IO Double @@ -239,7 +239,7 @@ where density_only f' d' p' q' = if f' /= Density then error "distribution has no CDF" else alloca $ \r -> do - app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p vec q "density_mp" + app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p' vec q' "density_mp" r' <- peek r return r' 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 @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Distribution.Discrete --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -11,6 +11,8 @@ -- -- GSL discrete random distribution functions -- +-- +-- ----------------------------------------------------------------------------- module Numeric.GSL.Distribution.Discrete ( @@ -101,10 +103,10 @@ _ -> distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p where pdf_only f' d' p' x' = if f' /= Density then error "no CDF" - else distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p + else distribution_dist_one_param (fromei f') (fromei d') (fromIntegral x') p' no_inverse f' d' p' x' = if (f' == LowInv || f' == UppInv) then error "No inverse CDF" - else distribution_dist_one_param (fromei f) (fromei d) (fromIntegral x) p + else distribution_dist_one_param (fromei f') (fromei d') (fromIntegral x') p' foreign import ccall "distribution-aux.h discrete1" distribution_discrete_one_param :: CInt -> CInt -> Double -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete1_dist" distribution_dist_one_param :: CInt -> CInt -> CUInt -> Double -> IO Double @@ -168,7 +170,7 @@ HyperGeometric -> no_inverse f d p1 p2 p3 x where no_inverse f' d' p1' p2' p3' x' = if (f' == LowInv || f' == UppInv) then error "No inverse CDF" - else distribution_dist_three_param (fromei f') (fromei d') (fromIntegral x') (fromIntegral p1') (fromIntegral p1') (fromIntegral p1') + else distribution_dist_three_param (fromei f') (fromei d') (fromIntegral x') (fromIntegral p1') (fromIntegral p2') (fromIntegral p3') foreign import ccall "distribution-aux.h discrete3" distribution_discrete_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> Ptr CUInt -> IO CInt foreign import ccall "distribution-aux.h discrete3_dist" distribution_dist_three_param :: CInt -> CInt -> CUInt -> CUInt -> CUInt -> CUInt -> IO Double @@ -198,7 +200,7 @@ where density_only f' d' p' q' = if f' /= Density then error "distribution has no CDF" else alloca $ \r -> do - app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p vec (mapVector (\x -> (fromIntegral x) :: CUInt) q) "density_mp" + app2 (distribution_dist_multi_param (fromei f') (fromei d') r) vec p' vec (mapVector (\x -> (fromIntegral x) :: CUInt) q') "density_mp" r' <- peek r return r' diff --git a/lib/Numeric/GSL/Fitting/Linear.hs b/lib/Numeric/GSL/Fitting/Linear.hs --- a/lib/Numeric/GSL/Fitting/Linear.hs +++ b/lib/Numeric/GSL/Fitting/Linear.hs @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Fitting.Linear --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL linear regression functions +-- +-- -- ----------------------------------------------------------------------------- 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 @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Histogram --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -11,24 +11,34 @@ -- -- GSL histogram functions -- +-- +-- ----------------------------------------------------------------------------- module Numeric.GSL.Histogram ( + -- * Creation Histogram , emptyRanges, emptyLimits , fromRanges, fromLimits + -- * Loading , addList, addVector, addListWeighted, addVectorWeighted - , toVectors + -- * Marshalling + , toVectors, fromVectors + -- * Information , getBin, getRange , getMax, getMin, getBins + -- * Querying , find - , count, prob + , count, prob, countInstance, probability , maxVal, maxBin, minVal, minBin + -- * Statistics , mean, stddev, sum , equalBins + -- * Mathematics , add, subtract, multiply, divide, shift, scale + -- * Files , fwriteHistogram, freadHistogram, fprintfHistogram, fscanfHistogram - -- + -- * PDF , HistogramPDF , fromHistogram , sample @@ -95,14 +105,6 @@ -} ----------------------------------------------------------------------------- -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 @@ -111,8 +113,7 @@ get = do b <- get c <- get - return $! addVectorWeighted (emptyRanges b) (middle b) c - where middle = fromList . map (\(x1,x2) -> (x1 + x2)/2) . vectorToTuples + return $! fromVectors b c ----------------------------------------------------------------------------- @@ -174,8 +175,27 @@ return h ----------------------------------------------------------------------------- +{- +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)) --- | extract the ranges and bins +middle = fromList . map (\(x1,x2) -> (x1 + x2)/2) . vectorToTuples +-} +----------------------------------------------------------------------------- + +-- | create a histogram from the ranges and bin weights +fromVectors :: Vector Double -> Vector Double -> Histogram +fromVectors r b = unsafePerformIO $ do + h@(H _ h') <- fromRangesIO r + app2 (\rs r' bs b' -> withForeignPtr h' $ \h'' -> histogram_from_vectors h'' rs r' bs b') vec r vec b "fromVectors" + return h + +foreign import ccall "gsl-histogram.h from_vectors" histogram_from_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt + +-- | extract the ranges and bin weights toVectors :: Histogram -> (Vector Double,Vector Double) -- ^ (ranges,bins) toVectors (H b h) = unsafePerformIO $ do rs <- createVector (b+1) @@ -313,6 +333,14 @@ else return Nothing foreign import ccall "gsl-histogram.h gsl_histogram_find" histogram_find :: HistHandle -> Double -> Ptr CInt -> IO CInt + +-- | find the number of occurences for the input +countInstance :: Histogram -> Double -> Double +countInstance h x = let Just x' = find h x in getBin h x' + +-- | find the probability of the input +probability :: Histogram -> Double -> Double +probability h x = let Just x' = find h x in (getBin h x') / (sum h) -- | find the number of occurences for each element of the input vector count :: Histogram -> Vector Double -> Vector Double 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 @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Histogram2D --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -11,25 +11,35 @@ -- -- GSL 2D histogram functions -- +-- +-- ----------------------------------------------------------------------------- module Numeric.GSL.Histogram2D ( + -- * Creation Histogram2D , emptyRanges, emptyLimits , fromRanges, fromLimits + -- * Loading , addList, addVector, addListWeighted, addVectorWeighted - , toMatrix + -- * Marshalling + , toMatrix, fromMatrix + -- * Information , getBin, getXRange, getYRange , getXMax, getYMax, getXMin, getYMin, getXBins, getYBins , reset + -- * Querying , find - , count, prob + , count, prob, probPaired, countPaired, countInstance, probability , maxVal, maxBin, minVal, minBin + -- * Statistics , xmean, ymean, xstddev, ystddev, covariance, sum , equalBins + -- * Mathematics , add, subtract, multiply, divide, shift, scale + -- * Files , fwriteHistogram2D, freadHistogram2D, fprintfHistogram2D, fscanfHistogram2D - -- + -- * PDF , Histogram2DPDF , fromHistogram2D , sample @@ -50,6 +60,8 @@ import Data.Binary import Foreign hiding(shift) +--import Foreign.Storable +--import Foreign.Ptr --import Foreign.ForeignPtr --import Foreign.Marshal.Alloc(alloca) import Foreign.C.Types(CInt,CChar) @@ -65,6 +77,22 @@ ----------------------------------------------------------------------------- +-- This Storable instance for Complex is binary compatible with C99, C++ and +-- Fortran complex data types. +instance (Storable a) => Storable (a,a) where + sizeOf z = 2 * sizeOf (fst z) + alignment z = alignment (fst z) + peek p = do let q = castPtr p + r <- peek q + i <- peekElemOff q 1 + return (r,i) + poke p (r,i) = do let q = (castPtr p) + poke q r + pokeElemOff q 1 i + + +----------------------------------------------------------------------------- + data Hist2D type Hist2DHandle = Ptr Hist2D -- | A histogram structure @@ -190,13 +218,13 @@ foreign import ccall "histogram-aux.h to_matrix" histogram_to_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt ----------------------------------------------------------------------------- - +{- 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 @@ -204,12 +232,13 @@ -> 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 (\(a,b) d -> (a,b,d)) xy w' +fromMatrix x y w = unsafePerformIO $ do + h@(H _ _ h') <- fromRangesIO x y + app3 (\xs x' ys y' rs cs b -> withForeignPtr h' $ \h'' -> histogram_from_matrix h'' xs x' ys y' rs cs b) vec x vec y mat w "fromMatrix" + return h +foreign import ccall "histogram-aux.h from_matrix" histogram_from_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt + ----------------------------------------------------------------------------- -- | create a copy of a histogram @@ -372,6 +401,14 @@ foreign import ccall "gsl-histogram2d.h gsl_histogram2d_find" histogram2d_find :: Hist2DHandle -> Double -> Double -> Ptr CInt -> Ptr CInt -> IO CInt +-- | find the number of occurences for the input +countInstance :: Histogram2D -> (Double,Double) -> Double +countInstance h x = let Just x' = find h x in getBin h x' + +-- | find the probability of the input +probability :: Histogram2D -> (Double,Double) -> Double +probability h x = let Just x' = find h x in (getBin h x') / (sum h) + -- | find the number of occurences for each element of the input vector count :: Histogram2D -> (Vector Double, Vector Double) -> Vector Double count (H _ _ h) (x,y) = unsafePerformIO $ do @@ -385,6 +422,21 @@ prob :: Histogram2D -> (Vector Double,Vector Double) -> Vector Double prob h z = (count h z) / (scalar $ sum h) +-- | find the number of occurences for each element of the input vector +countPaired :: Histogram2D -> Vector (Double,Double) -> Vector Double +countPaired (H _ _ h) x = unsafePerformIO $ do + r <- createVector $ dim x + app2 (\xs' x' rs' r' -> withForeignPtr h $ \h' -> histogram2d_count_pair h' xs' x' rs' r') vec x vec r "histogram2d_count_pair" + return r + +foreign import ccall "histogram-aux.h hist2d_count_pair" histogram2d_count_pair :: Hist2DHandle -> CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> IO CInt + +-- | find the joint probability of occuring for each element of the input vector pair +probPaired :: Histogram2D -> (Vector (Double,Double)) -> Vector Double +probPaired h z = let l = unzip $ toList z + z' = (fromList $ fst l,fromList $ snd l) + in prob h z' + ----------------------------------------------------------------------------- -- | the maximum value contained in the bins @@ -597,4 +649,15 @@ ----------------------------------------------------------------------------- ----------------------------------------------------------------------------- +----------------------------------------------------------------------------- + +unzipPair :: Vector (Double,Double) -> (Vector Double,Vector Double) +unzipPair v = unsafePerformIO $ do + a <- createVector $ dim v + b <- createVector $ dim v + app3 histogram2d_unzip_double_pair vec v vec a vec b "unzipPair" + return (a,b) + +foreign import ccall "histogram-aux.h unzip_double_pair" histogram2d_unzip_double_pair :: CInt -> Ptr (Double,Double) -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt + ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Permutation.hs b/lib/Numeric/GSL/Permutation.hs --- a/lib/Numeric/GSL/Permutation.hs +++ b/lib/Numeric/GSL/Permutation.hs @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Permutation --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL permutation functions +-- +-- -- ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Sort.hs b/lib/Numeric/GSL/Sort.hs --- a/lib/Numeric/GSL/Sort.hs +++ b/lib/Numeric/GSL/Sort.hs @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Sort --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL sorting functions +-- +-- -- ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/Statistics.hs b/lib/Numeric/GSL/Statistics.hs --- a/lib/Numeric/GSL/Statistics.hs +++ b/lib/Numeric/GSL/Statistics.hs @@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- | -- Module : Numeric.GSL.Statistics --- Copyright : (c) Alexander Vivian Hugh McPhail 2010 +-- Copyright : (c) A. V. H. McPhail 2010 -- License : GPL-style -- -- Maintainer : haskell.vivian.mcphail gmail com @@ -10,6 +10,8 @@ -- Portability : uses ffi -- -- GSL statistics functions +-- +-- -- ----------------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/histogram-aux.c b/lib/Numeric/GSL/histogram-aux.c --- a/lib/Numeric/GSL/histogram-aux.c +++ b/lib/Numeric/GSL/histogram-aux.c @@ -6,6 +6,14 @@ #include +int from_vectors(gsl_histogram * H, int rs, const double* r, int bs, const double* b) +{ + if (rs != bs + 1) return 2000; // BAD_SIZE + int i; + for (i = 0; i < rs - 1; i++) + gsl_histogram_accumulate(H,(r[i]+r[i+1])/2,b[i]); + return 0; +} int to_vectors(gsl_histogram * H, int rs, double* r, int bs, double* b) { @@ -76,6 +84,18 @@ ////////////////////// +int from_matrix(gsl_histogram2d * H, int rxs, const double* rx, int rys, const double* ry, int bx, int by, const double* b) +{ + if (rxs != bx + 1) return 2000; // BAD_SIZE + if (rys != by + 1) return 2000; // BAD_SIZE + + int i, j; + for (i = 0; i < rxs - 1; i++) + for (j = 0; j