hstatistics-0.1.0.4: lib/Numeric/GSL/Histogram.hs
{-# OPTIONS_GHC -fglasgow-exts #-}
-----------------------------------------------------------------------------
-- |
-- Module : Numeric.GSL.Histogram
-- Copyright : (c) Alexander Vivian Hugh McPhail 2010
-- License : GPL-style
--
-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com
-- Stability : provisional
-- Portability : uses ffi
--
-- GSL histogram functions
--
-----------------------------------------------------------------------------
module Numeric.GSL.Histogram (
Histogram
, fromRanges, fromLimits
, addList, addVector, addListWeighted, addVectorWeighted
, toVectors
, getBin, getRange
, getMax, getMin, getBins
, find
, maxVal, maxBin, minVal, minBin
, mean, stddev, sum
, equalBins
, add, subtract, multiply, divide, shift, scale
, fwriteHistogram, freadHistogram, fprintfHistogram, fscanfHistogram
--
, HistogramPDF
, fromHistogram
, sample
) where
-----------------------------------------------------------------------------
import Data.Packed.Vector
--import Data.Packed.Matrix hiding(toLists)
import Data.Packed.Development
--import Numeric.LinearAlgebra.Linear
--import Control.Monad
import Foreign hiding(shift)
--import Foreign.ForeignPtr
--import Foreign.Marshal.Alloc(alloca)
import Foreign.C.Types(CInt,CChar)
--import Foreign.C.String(newCString,peekCString)
import Foreign.C.String(newCString)
--import Control.Monad(when)
--import GHC.ForeignPtr (mallocPlainForeignPtrBytes)
--import GHC.Base
--import GHC.IOBase
import Prelude hiding(subtract,sum)
-----------------------------------------------------------------------------
data Hist
type HistHandle = Ptr Hist
-- | A histogram structure
data Histogram = H { hdim :: {-# UNPACK #-} !Int -- ^ number of bins
, hist :: {-# UNPACK #-} !(ForeignPtr Hist) }
data PDF
type PDFHandle = Ptr PDF
-- | A histogram-derived cumulative distribution function (CDF)
data HistogramPDF = P { pdf :: {-# UNPACK #-} !(ForeignPtr PDF)}
-----------------------------------------------------------------------------
instance Eq Histogram where
(==) = equalBins
{-
instance Num Histogram where
(+) = add
(-) = subtract
negate = flip scale (-1.0)
(*) = multiply
signum = error "can't signumm Histogram"
abs = error "can't abs Histogram"
fromInteger x = fromLimits (fromInteger x) (0,1)
instance Fractional Histogram where
fromRational x = fromLimits (round x) (0,fromRational x)
(/) = divide
-}
-----------------------------------------------------------------------------
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 ())
-----------------------------------------------------------------------------
-- | create a histogram with n bins from ranges (x0->x1),(x1->x2)..(xn->xn+1)
fromRangesIO :: Vector Double -> IO Histogram
fromRangesIO v = do
let sz = fromIntegral $ dim v - 1
h <- histogram_new sz
h' <- newForeignPtr histogram_free h
app1 (\d p -> withForeignPtr h' $ \f -> histogram_set_ranges f (fromIntegral d) p) vec v "fromRanges"
return $ H (fromIntegral sz) h'
-- | create a histogram with n bins and lower and upper limits
fromLimitsIO :: Int -> (Double,Double) -> IO Histogram
fromLimitsIO n (l,u) = do
h <- histogram_new (fromIntegral n)
h' <- newForeignPtr histogram_free h
check "set_ranges_uniform" $ withForeignPtr h' (\f -> histogram_set_ranges_uniform f l u)
return $ H n h'
foreign import ccall "gsl-histogram.h gsl_histogram_set_ranges" histogram_set_ranges :: HistHandle -> CInt -> Ptr Double -> IO CInt
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) and increment from a vector
fromRanges :: Vector Double -- ^ the ranges
-> Vector Double -- ^ the data
-> Histogram -- ^ result
fromRanges r d = unsafePerformIO $ do
h <- fromRangesIO r
incrementVectorIO h d
return h
-- | create a histogram with n bins and lower and upper limits and increment from a vector
fromLimits :: Int -- ^ bins
-> (Double,Double) -- ^ lower and upper limits
-> Vector Double -- ^ the data
-> Histogram -- ^ result
fromLimits n r d = unsafePerformIO $ do
h <- fromLimitsIO n r
incrementVectorIO h d
return h
-----------------------------------------------------------------------------
-- | extract the ranges and bins
toVectors :: Histogram -> (Vector Double,Vector Double) -- ^ (ranges,bins)
toVectors (H b h) = unsafePerformIO $ do
rs <- createVector (b+1)
bs <- createVector b
app2 (\s1 p1 s2 p2 -> withForeignPtr h $ \f -> histogram_to_vectors f s1 p1 s2 p2) vec rs vec bs "toVectors"
return (rs,bs)
foreign import ccall "gsl-histogram.h to_vectors" histogram_to_vectors :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | create a copy of a histogram
cloneHistogram :: Histogram -> IO Histogram
cloneHistogram (H b h) = do
h' <- withForeignPtr h histogram_clone
h'' <- newForeignPtr histogram_free h'
return $ H b h''
foreign import ccall "gsl-histogram.h gsl_histogram_clone" histogram_clone :: HistHandle -> IO HistHandle
-----------------------------------------------------------------------------
-- | adds 1.0 to the correct bin for each element of the list
addList :: Histogram -> [Double] -> Histogram
addList h xs = unsafePerformIO $ do
h' <- cloneHistogram h
incrementListIO h' xs
return h'
-- | adds 1.0 to the correct bin for each element of the vector
addVector :: Histogram -> Vector Double -> Histogram
addVector h v = unsafePerformIO $ do
h' <- cloneHistogram h
incrementVectorIO h' v
return h'
-- | adds the appropriate weight for each element of the list
addListWeighted :: Histogram -> [(Double,Double)] -> Histogram
addListWeighted h xs = unsafePerformIO $ do
h' <- cloneHistogram h
accumulateListIO h' xs
return h'
-- | adds the appropriate weight for each element of the list
addVectorWeighted :: Histogram -> Vector Double -> Vector Double -> Histogram
addVectorWeighted h v w = unsafePerformIO $ do
h' <- cloneHistogram h
accumulateVectorIO h' v w
return h'
-- | add 1.0 to the correct bin, fails silently if the value is outside the range
incrementIO :: Histogram -> Double -> IO ()
incrementIO (H _ h) x = withForeignPtr h (\f -> check "increment" $ histogram_increment f x)
-- | add 1.0 to the correct bin for each element of the vector, fails silently if the value is outside the range
incrementListIO :: Histogram -> [Double] -> IO ()
incrementListIO (H _ h) xs = withForeignPtr h (\f -> mapM_ (histogram_increment f) xs)
-- | add 1.0 to the correct bin for each element of the vector, fails silently if the value is outside the range
incrementVectorIO :: Histogram -> Vector Double -> IO ()
incrementVectorIO (H _ h) v = do
app1 (\s p -> withForeignPtr h (\f -> histogram_increment_vector f s p)) vec v "incrementVector"
return ()
-- | adds the weight (second Double) to the bin appropriate for the value (first Double)
accumulateIO :: Histogram -> Double -> Double -> IO ()
accumulateIO (H _ h) x w = withForeignPtr h (\f -> check "accumulate" $ histogram_accumulate f x w)
-- | add the weight (snd) to the correct bin for each (fst) element of the list, fails silently if the value is outside the range
accumulateListIO :: Histogram -> [(Double,Double)] -> IO ()
accumulateListIO (H _ h) xs = do
withForeignPtr h (\f -> mapM_ (\(x,w) -> histogram_accumulate f x w) xs)
return ()
-- | add the weight (second vector) to the correct bin for each element of the first vector, fails silently if the value is outside the range
accumulateVectorIO :: Histogram -> Vector Double -> Vector Double -> IO ()
accumulateVectorIO (H _ h) v w = do
app2 (\s1 p1 s2 p2 -> withForeignPtr h (\f -> histogram_accumulate_vector f s1 p1 s2 p2)) vec v vec w "accumulateVector"
return ()
-- | returns the contents of the i-th bin
getBin :: Histogram -> Int -> Double
getBin (H _ h) b = unsafePerformIO $ do
withForeignPtr h (\f -> histogram_get f (fromIntegral b))
-- | returns the upper and lower limits of the i-th bin
getRange :: Histogram -> Int -> (Double,Double)
getRange (H _ h) b = unsafePerformIO $ do
alloca $ \l ->
alloca $ \u -> do
check "get_range" $ withForeignPtr h (\f -> histogram_get_range f (fromIntegral b) l u)
l' <- peek l
u' <- peek u
return (l',u')
-- | the maximum upper range limit
getMax :: Histogram -> Double
getMax (H _ h) = unsafePerformIO $ withForeignPtr h histogram_max
-- | the minimum lower range limit
getMin :: Histogram -> Double
getMin (H _ h) = unsafePerformIO $ withForeignPtr h histogram_min
-- | the number of bins
getBins :: Histogram -> Int
getBins (H b _) = b
-- | reset all the bins to zero
reset :: Histogram -> IO ()
reset (H _ h) = withForeignPtr h histogram_reset
foreign import ccall "gsl-histogram.h gsl_histogram_increment" histogram_increment :: HistHandle -> Double -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_accumulate" histogram_accumulate :: HistHandle -> Double -> Double -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_get" histogram_get :: HistHandle -> CInt -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_get_range" histogram_get_range :: HistHandle -> CInt -> Ptr Double -> Ptr Double -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_max" histogram_max :: HistHandle -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_min" histogram_min :: HistHandle -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_bins" histogram_bins :: HistHandle -> IO Int
foreign import ccall "gsl-histogram.h gsl_histogram_reset" histogram_reset :: HistHandle -> IO ()
foreign import ccall "histogram-aux.h increment_vector" histogram_increment_vector :: HistHandle -> CInt -> Ptr Double -> IO CInt
foreign import ccall "histogram-aux.h accumulate_vector" histogram_accumulate_vector :: HistHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | find the bin corresponding to the value
find :: Histogram -> Double -> Maybe Int
find (H _ h) x = unsafePerformIO $ do
alloca $ \b -> do
err <- withForeignPtr h (\f -> histogram_find f x b)
if err == 0
then do
b' <- peek b
return $ Just $ fromIntegral b'
else return Nothing
foreign import ccall "gsl-histogram.h gsl_histogram_find" histogram_find :: HistHandle -> Double -> Ptr CInt -> IO CInt
-----------------------------------------------------------------------------
-- | the maximum value contained in the bins
maxVal :: Histogram -> Double
maxVal (H _ h) = unsafePerformIO $ withForeignPtr h histogram_max_val
-- | the index of the bin containing the maximum value
maxBin :: Histogram -> Int
maxBin (H _ h) = unsafePerformIO $ do
i <- withForeignPtr h histogram_max_bin
return $ fromIntegral i
-- | the minimum value contained in the bins
minVal :: Histogram -> Double
minVal (H _ h) = unsafePerformIO $ withForeignPtr h histogram_min_val
-- | the index of the bin containing the minimum value
minBin :: Histogram -> Int
minBin (H _ h) = unsafePerformIO $ do
i <- withForeignPtr h histogram_min_bin
return $ fromIntegral i
-- | the mean of the values, accuracy limited by bin width
mean :: Histogram -> Double
mean (H _ h) = unsafePerformIO $ withForeignPtr h histogram_mean
-- | the standard deviation of the values, accuracy limited by bin width
stddev :: Histogram -> Double
stddev (H _ h) = unsafePerformIO $ withForeignPtr h histogram_sigma
-- | the sum of the values, accuracy limited by bin width
sum :: Histogram -> Double
sum (H _ h) = unsafePerformIO $ withForeignPtr h histogram_sum
foreign import ccall "gsl-histogram.h gsl_histogram_max_val" histogram_max_val :: HistHandle -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_max_bin" histogram_max_bin :: HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_min_val" histogram_min_val :: HistHandle -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_min_bin" histogram_min_bin :: HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_mean" histogram_mean :: HistHandle -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_sigma" histogram_sigma :: HistHandle -> IO Double
foreign import ccall "gsl-histogram.h gsl_histogram_sum" histogram_sum :: HistHandle -> IO Double
-----------------------------------------------------------------------------
-- | returns True of all the individual bin ranges of the two histograms are identical
equalBins :: Histogram -> Histogram -> Bool
equalBins (H _ h1) (H _ h2) = unsafePerformIO $ do
i <- withForeignPtr h1 $ \p1 -> do
withForeignPtr h2 $ \p2 -> histogram_equal_bins p1 p2
if (fromIntegral i) == (1 :: Int)
then return True
else return False
-- | adds the contents of the bins of the second histogram to the first
add :: Histogram -> Histogram -> Histogram
add d (H _ s) = unsafePerformIO $ do
h@(H _ d') <- cloneHistogram d
check "add" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram_add dp sp
return h
-- | subtracts the contents of the bins of the second histogram from the first
subtract :: Histogram -> Histogram -> Histogram
subtract d (H _ s) = unsafePerformIO $ do
h@(H _ d') <- cloneHistogram d
check "subtract" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram_sub dp sp
return h
-- | multiplies the contents of the bins of the second histogram by the first
multiply :: Histogram -> Histogram -> Histogram
multiply d (H _ s) = unsafePerformIO $ do
h@(H _ d') <- cloneHistogram d
check "multiply" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram_mul dp sp
return h
-- | divides the contents of the bins of the first histogram by the second
divide :: Histogram -> Histogram -> Histogram
divide d (H _ s) = unsafePerformIO $ do
h@(H _ d') <- cloneHistogram d
check "divide" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram_div dp sp
return h
-- | multiplies the contents of the bins by a constant
scale :: Histogram -> Double -> Histogram
scale d s = unsafePerformIO $ do
h@(H _ d') <- cloneHistogram d
check "scale" $
withForeignPtr d' $ (\f -> histogram_scale f s)
return h
-- | adds a constant to the contents of the bins
shift :: Histogram -> Double -> Histogram
shift d s = unsafePerformIO $ do
h@(H _ d') <- cloneHistogram d
check "shift" $
withForeignPtr d' $ (\f -> histogram_shift f s)
return h
foreign import ccall "gsl-histogram.h gsl_histogram_equal_bins_p" histogram_equal_bins :: HistHandle -> HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_add" histogram_add :: HistHandle -> HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_sub" histogram_sub :: HistHandle -> HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_mul" histogram_mul :: HistHandle -> HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_div" histogram_div :: HistHandle -> HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_scale" histogram_scale :: HistHandle -> Double -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_shift" histogram_shift :: HistHandle -> Double -> IO CInt
-----------------------------------------------------------------------------
-- | write a histogram in the native binary format (may not be portable)
fwriteHistogram :: FilePath -> Histogram -> IO ()
fwriteHistogram fn (H _ h) = do
cn <- newCString fn
check "fwriteHistogram" $
withForeignPtr h $ histogram_fwrite cn
free cn
-- | read a histogram in the native binary format, number of bins must be known
freadHistogram :: FilePath -> Int -> IO Histogram
freadHistogram fn b = do
h <- histogram_new (fromIntegral b)
h' <- newForeignPtr histogram_free h
cn <- newCString fn
check "freadHistogram" $
withForeignPtr h' $ histogram_fread cn
return $ H b h'
-- | saves the histogram with the given formats (%f,%e,%g) for ranges and bins
-- each line comprises: range[i] range[i+1] bin[i]
fprintfHistogram :: FilePath -> String -> String -> Histogram -> IO ()
fprintfHistogram fn fr fb (H _ h) = do
cn <- newCString fn
cr <- newCString fr
cb <- newCString fb
check "fprintfHistogram" $
withForeignPtr h $ histogram_fprintf cn cr cb
free cn
free cr
free cb
return ()
-- | reads formatted data as written by fprintf, the number of bins must be known in advance
fscanfHistogram :: FilePath -> Int -> IO Histogram
fscanfHistogram fn b = do
h <- histogram_new (fromIntegral b)
h' <- newForeignPtr histogram_free h
cn <- newCString fn
check "fscanfHistogram" $
withForeignPtr h' $ histogram_fscanf cn
return $ H b h'
foreign import ccall "histogram-aux.h hist_fwrite" histogram_fwrite :: Ptr CChar -> HistHandle -> IO CInt
foreign import ccall "histogram-aux.h hist_fread" histogram_fread :: Ptr CChar -> HistHandle -> IO CInt
foreign import ccall "histogram-aux.h hist_fprintf" histogram_fprintf :: Ptr CChar -> Ptr CChar -> Ptr CChar -> HistHandle -> IO CInt
foreign import ccall "histogram-aux.h hist_fscanf" histogram_fscanf :: Ptr CChar -> HistHandle -> IO CInt
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
foreign import ccall "gsl-histogram.h gsl_histogram_pdf_alloc" histogram_pdf_new :: CInt -> IO PDFHandle
foreign import ccall "gsl-histogram.h &gsl_histogram_pdf_free" histogram_pdf_free :: FunPtr (PDFHandle -> IO ())
-----------------------------------------------------------------------------
-- | create a histogram PDF from a histogram
fromHistogram :: Histogram -> HistogramPDF
fromHistogram (H b h) = unsafePerformIO $ do
p <- histogram_pdf_new $ fromIntegral b
p' <- newForeignPtr histogram_pdf_free p
withForeignPtr p' $ \p'' ->
withForeignPtr h $ \h' -> check "pdf_init" $ histogram_pdf_init p'' h'
return $ P p'
-- | given a randomm from the uniform distribution [0,1], draw a random sample from the PDF
sample :: HistogramPDF -> Double
sample (P p) = unsafePerformIO $ withForeignPtr p $ \p' -> histogram_pdf_sample p'
foreign import ccall "gsl-histogram.h gsl_histogram_pdf_init" histogram_pdf_init :: PDFHandle -> HistHandle -> IO CInt
foreign import ccall "gsl-histogram.h gsl_histogram_pdf_sample" histogram_pdf_sample :: PDFHandle -> IO Double
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------