hstatistics-0.1.0.2: lib/Numeric/GSL/Histogram2D.hs
{-# OPTIONS_GHC -fglasgow-exts #-}
-----------------------------------------------------------------------------
-- |
-- Module : Numeric.GSL.Histogram2D
-- Copyright : (c) Alexander Vivian Hugh McPhail 2010
-- License : GPL-style
--
-- Maintainer : haskell.vivian.mcphail <at> gmail <dot> com
-- Stability : provisional
-- Portability : uses ffi
--
-- GSL 2D histogram functions
--
-----------------------------------------------------------------------------
module Numeric.GSL.Histogram2D (
Histogram2D
, fromRanges, fromLimits
, addList, addVector, addListWeighted, addVectorWeighted
, toMatrix
, getBin, getXRange, getYRange
, getXMax, getYMax, getXMin, getYMin, getXBins, getYBins
, reset
, find
, maxVal, maxBin, minVal, minBin
, xmean, ymean, xstddev, ystddev, covariance, sum
, equalBins
, add, subtract, multiply, divide, shift, scale
, fwriteHistogram2D, freadHistogram2D, fprintfHistogram2D, fscanfHistogram2D
--
, Histogram2DPDF
, fromHistogram2D
, sample
) where
-----------------------------------------------------------------------------
import Data.Packed.Vector
import Data.Packed.Matrix
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 Control.Monad(when)
--import GHC.ForeignPtr (mallocPlainForeignPtrBytes)
--import GHC.Base
--import GHC.IOBase
import Prelude hiding(subtract,sum)
-----------------------------------------------------------------------------
data Hist2D
type Hist2DHandle = Ptr Hist2D
-- | A histogram structure
data Histogram2D = H { hxdim :: {-# UNPACK #-} !Int -- ^ number of bins
, hydim :: {-# UNPACK #-} !Int -- ^ number of bins
, hist :: {-# UNPACK #-} !(ForeignPtr Hist2D) }
data PDF
type PDFHandle = Ptr PDF
-- | A histogram-derived cumulative distribution function (CDF)
data Histogram2DPDF = P { pdf :: {-# UNPACK #-} !(ForeignPtr PDF)}
-----------------------------------------------------------------------------
instance Eq Histogram2D where
(==) = equalBins
{-
instance Num Histogram2D where
(+) = add
(-) = subtract
negate = flip scale (-1.0)
(*) = multiply
signum = error "can't signumm Histogram2D"
abs = error "can't abs Histogram2D"
fromInteger x = fromLimits (fromInteger x) (0,1)
instance Fractional Histogram2D where
fromRational x = fromLimits (round x) (0,fromRational x)
(/) = divide
-}
-----------------------------------------------------------------------------
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 ())
-----------------------------------------------------------------------------
-- | create a histogram with n bins from ranges (x0->x1),(x1->x2),..,(xn->xn+1) and y0,..,yn+1
fromRangesIO :: Vector Double -> Vector Double -> IO Histogram2D
fromRangesIO v w = do
let sx = fromIntegral $ dim v - 1
let sy = fromIntegral $ dim w - 1
h <- histogram2d_new sx sy
h' <- newForeignPtr histogram2d_free h
app2 (\xs xp ys yp -> withForeignPtr h' (\f -> histogram2d_set_ranges f xp xs yp ys)) vec v vec w "fromRanges"
return $ H (fromIntegral sx) (fromIntegral sy) h'
-- | create a histogram with n bins and lower and upper limits
fromLimitsIO :: Int -> Int -- ^ number of bins
-> (Double,Double) -- ^ xmin, xmax
-> (Double,Double) -- ^ ymin, ymax
-> IO Histogram2D
fromLimitsIO nx ny (lx,ux) (uy,ly) = do
h <- histogram2d_new (fromIntegral nx) (fromIntegral ny)
h' <- newForeignPtr histogram2d_free h
withForeignPtr h' (\f -> histogram2d_set_ranges_uniform f lx ux ly uy)
return $ H nx ny h'
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_set_ranges" histogram2d_set_ranges :: Hist2DHandle -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_set_ranges_uniform" histogram2d_set_ranges_uniform :: Hist2DHandle -> Double -> Double -> 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 x ranges
-> Vector Double -- ^ the y ranges
-> [(Double,Double)] -- ^ the data
-> Histogram2D -- ^ result
fromRanges rx ry d = unsafePerformIO $ do
h <- fromRangesIO rx ry
incrementListIO h d
return h
-- | create a histogram with n bins and lower and upper limits and increment from a vector
fromLimits :: Int -> Int -- ^ bins
-> (Double,Double) -- ^ x lower and upper limits
-> (Double,Double) -- ^ y lower and upper limits
-> [(Double,Double)] -- ^ the data
-> Histogram2D -- ^ result
fromLimits nx ny rx ry d = unsafePerformIO $ do
h <- fromLimitsIO nx ny rx ry
incrementListIO h d
return h
-----------------------------------------------------------------------------
-- | extract the ranges and bins
toMatrix :: Histogram2D -> (Vector Double,Vector Double,Matrix Double) -- ^ (ranges,bins)
toMatrix (H bx by h) = unsafePerformIO $ do
rx <- createVector (bx+1)
ry <- createVector (by+1)
bs <- createMatrix RowMajor bx by
app3 (\s1 p1 s2 p2 sx sy p -> withForeignPtr h $ \f -> histogram_to_matrix f s1 p1 s2 p2 sx sy p) vec rx vec ry mat bs "toMatrix"
return (rx,ry,bs)
foreign import ccall "histogram-aux.h to_matrix" histogram_to_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | create a copy of a histogram
cloneHistogram2D :: Histogram2D -> IO Histogram2D
cloneHistogram2D (H nx ny h) = do
h' <- withForeignPtr h histogram2d_clone
h'' <- newForeignPtr histogram2d_free h'
return $ H nx ny h''
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_clone" histogram2d_clone :: Hist2DHandle -> IO Hist2DHandle
-----------------------------------------------------------------------------
-- | add 1.0 to the correct bin for each element of the list, fails silently if the value is outside the range
addList :: Histogram2D -> [(Double,Double)] -> Histogram2D
addList h xs = unsafePerformIO $ do
h' <- cloneHistogram2D h
incrementListIO h' xs
return h'
-- | add 1.0 to the correct bin for each element of the vector pair, fails silently if the value is outside the range
addVector :: Histogram2D -> Vector Double -> Vector Double -> Histogram2D
addVector h x y = unsafePerformIO $ do
h' <- cloneHistogram2D h
incrementVectorIO h' x y
return h'
-- add the appropriate weight for each element of the list, fails silently if the value is outside the range
addListWeighted :: Histogram2D -> [(Double,Double,Double)] -> Histogram2D
addListWeighted h xs = unsafePerformIO $ do
h' <- cloneHistogram2D h
accumulateListIO h' xs
return h'
-- add the appropriate weight for each element of the vector pair, fails silently if the value is outside the range
addVectorWeighted :: Histogram2D -> Vector Double -> Vector Double -> Vector Double -> Histogram2D
addVectorWeighted h x y w = unsafePerformIO $ do
h' <- cloneHistogram2D h
accumulateVectorIO h' x y w
return h'
-----------------------------------------------------------------------------
-- | add 1.0 to the correct bin, fails silently if the value is outside the range
incrementIO :: Histogram2D -> Double -> Double -> IO ()
incrementIO (H _ _ h) x y = do
withForeignPtr h (\f -> histogram2d_increment f x y)
return ()
-- | add 1.0 to the correct bin for each element of the vector pair, fails silently if the value is outside the range
incrementVectorIO :: Histogram2D -> Vector Double -> Vector Double -> IO ()
incrementVectorIO (H _ _ h) x y = do
app2 (\xs xp ys yp -> withForeignPtr h (\f -> histogram2d_increment_matrix f xs xp ys yp)) vec x vec y "incrementVector"
return ()
-- | add 1.0 to the correct bin for each element of the list, fails silently if the value is outside the range
incrementListIO :: Histogram2D -> [(Double,Double)] -> IO ()
incrementListIO (H _ _ h) zs = withForeignPtr h (\f -> mapM_ (\(x,y) -> histogram2d_increment f x y) zs)
-- | Adds the weight (third Double) to the bin appropriate for the value (first two Doubles)
accumulateIO :: Histogram2D -> Double -> Double -> Double -> IO ()
accumulateIO (H _ _ h) x y w = do
withForeignPtr h (\f -> histogram2d_accumulate f x y w)
return ()
-- | add the weight (third) to the correct bin for each vector pair element, fails silently if the value is outside the range
accumulateVectorIO :: Histogram2D -> Vector Double -> Vector Double -> Vector Double -> IO ()
accumulateVectorIO (H _ _ h) x y w = do
app3 (\xs xp ys yp ws wp -> withForeignPtr h (\f -> histogram2d_accumulate_matrix f xs xp ys yp ws wp)) vec x vec y vec w "accumulateVector"
return ()
-- | 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 :: Histogram2D -> [(Double,Double,Double)] -> IO ()
accumulateListIO (H _ _ h) zs = withForeignPtr h (\f -> mapM_ (\(x,y,w) -> histogram2d_accumulate f x y w) zs)
-- | returns the contents of the i-th bin
getBin :: Histogram2D -> Int -> Int -> Double
getBin (H _ _ h) bx by = unsafePerformIO $ do
withForeignPtr h (\f -> histogram2d_get f (fromIntegral bx) (fromIntegral by))
-- | returns the upper and lower limits in the first dimension of the i-th bin
getXRange :: Histogram2D -> Int -> (Double,Double)
getXRange (H _ _ h) b = unsafePerformIO $ do
alloca $ \l ->
alloca $ \u -> do
withForeignPtr h (\f -> histogram2d_get_xrange f (fromIntegral b) l u)
l' <- peek l
u' <- peek u
return (l',u')
-- | returns the upper and lower limits in the second dimension of the i-th bin
getYRange :: Histogram2D -> Int -> (Double,Double)
getYRange (H _ _ h) b = unsafePerformIO $ do
alloca $ \l ->
alloca $ \u -> do
withForeignPtr h (\f -> histogram2d_get_yrange f (fromIntegral b) l u)
l' <- peek l
u' <- peek u
return (l',u')
-- | the maximum upper range limit in the first dimension
getXMax :: Histogram2D -> Double
getXMax (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmax
-- | the minimum lower range limit in the first dimension
getXMin :: Histogram2D -> Double
getXMin (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmin
-- | the number of binsin the first dimension
getXBins :: Histogram2D -> Int
getXBins (H bx _ _) = bx
-- | the maximum upper range limit in the first dimension
getYMax :: Histogram2D -> Double
getYMax (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymax
-- | the minimum lower range limit in the first dimension
getYMin :: Histogram2D -> Double
getYMin (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymin
-- | the number of binsin the first dimension
getYBins :: Histogram2D -> Int
getYBins (H _ by _) = by
-- | reset all the bins to zero
reset :: Histogram2D -> IO ()
reset (H _ _ h) = do
withForeignPtr h histogram2d_reset
return ()
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_increment" histogram2d_increment :: Hist2DHandle -> Double -> Double -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_accumulate" histogram2d_accumulate :: Hist2DHandle -> Double -> Double -> Double -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_get" histogram2d_get :: Hist2DHandle -> CInt -> CInt -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_get_xrange" histogram2d_get_xrange :: Hist2DHandle -> CInt -> Ptr Double -> Ptr Double -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_get_yrange" histogram2d_get_yrange :: Hist2DHandle -> CInt -> Ptr Double -> Ptr Double -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xmax" histogram2d_xmax :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xmin" histogram2d_xmin :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_nx" histogram2d_xn :: Hist2DHandle -> IO Int
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ymax" histogram2d_ymax :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ymin" histogram2d_ymin :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ny" histogram2d_yn :: Hist2DHandle -> IO Int
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_reset" histogram2d_reset :: Hist2DHandle -> IO ()
foreign import ccall "histogram-aux.h increment_matrix" histogram2d_increment_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
foreign import ccall "histogram-aux.h accumulate_matrix" histogram2d_accumulate_matrix :: Hist2DHandle -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
-----------------------------------------------------------------------------
-- | find the bin corresponding to the value
find :: Histogram2D -> (Double,Double) -> Maybe (Int,Int)
find (H _ _ h) (x,y) = unsafePerformIO $ do
alloca $ \bx ->
alloca $ \by -> do
err <- withForeignPtr h (\f -> histogram2d_find f x y bx by)
if err == 0
then do
bx' <- peek bx
by' <- peek by
return $ Just $ (fromIntegral bx',fromIntegral by')
else return Nothing
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_find" histogram2d_find :: Hist2DHandle -> Double -> Double -> Ptr CInt -> Ptr CInt -> IO CInt
-----------------------------------------------------------------------------
-- | the maximum value contained in the bins
maxVal :: Histogram2D -> Double
maxVal (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_max_val
-- | the index of the bin containing the maximum value
maxBin :: Histogram2D -> (Int,Int)
maxBin (H _ _ h) = unsafePerformIO $ do
alloca $ \bx ->
alloca $ \by -> do
withForeignPtr h (\f -> histogram2d_max_bin f bx by)
bx' <- peek bx
by' <- peek by
return $ (fromIntegral bx',fromIntegral by')
-- | the minimum value contained in the bins
minVal :: Histogram2D -> Double
minVal (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_min_val
-- | the index of the bin containing the minimum value
minBin :: Histogram2D -> (Int,Int)
minBin (H _ _ h) = unsafePerformIO $ do
alloca $ \bx ->
alloca $ \by -> do
withForeignPtr h (\f -> histogram2d_min_bin f bx by)
bx' <- peek bx
by' <- peek by
return $ (fromIntegral bx',fromIntegral by')
-- | the mean of the values in the first dimension, accuracy limited by bin width
xmean :: Histogram2D -> Double
xmean (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmean
-- | the mean of the values in the second dimension, accuracy limited by bin width
ymean :: Histogram2D -> Double
ymean (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymean
-- | the standard deviation of the values in thee first dimension, accuracy limited by bin width
xstddev :: Histogram2D -> Double
xstddev (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xsigma
-- | the standard deviation of the values in thee first dimension, accuracy limited by bin width
ystddev :: Histogram2D -> Double
ystddev (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ysigma
-- | the covariance of the first and second dimensions
covariance :: Histogram2D -> Double
covariance (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_cov
-- | the sum of the values, accuracy limited by bin width
sum :: Histogram2D -> Double
sum (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_sum
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_max_val" histogram2d_max_val :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_max_bin" histogram2d_max_bin :: Hist2DHandle -> Ptr CInt -> Ptr CInt -> IO ()
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_min_val" histogram2d_min_val :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_min_bin" histogram2d_min_bin :: Hist2DHandle -> Ptr CInt -> Ptr CInt -> IO ()
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xmean" histogram2d_xmean :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ymean" histogram2d_ymean :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_xsigma" histogram2d_xsigma :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_ysigma" histogram2d_ysigma :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_cov" histogram2d_cov :: Hist2DHandle -> IO Double
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_sum" histogram2d_sum :: Hist2DHandle -> IO Double
-----------------------------------------------------------------------------
-- | 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
withForeignPtr h2 $ \p2 -> histogram2d_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 :: Histogram2D -> Histogram2D -> Histogram2D
add d (H _ _ s) = unsafePerformIO $ do
h@(H _ _ d') <- cloneHistogram2D d
check "add" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram2d_add dp sp
return h
-- | subtracts the contents of the bins of the second histogram from the first
subtract :: Histogram2D -> Histogram2D -> Histogram2D
subtract d (H _ _ s) = unsafePerformIO $ do
h@(H _ _ d') <- cloneHistogram2D d
check "subtract" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram2d_sub dp sp
return h
-- | multiplies the contents of the bins of the second histogram by the first
multiply :: Histogram2D -> Histogram2D -> Histogram2D
multiply d (H _ _ s) = unsafePerformIO $ do
h@(H _ _ d') <- cloneHistogram2D d
check "multiply" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram2d_mul dp sp
return h
-- | divides the contents of the bins of the first histogram by the second
divide :: Histogram2D -> Histogram2D -> Histogram2D
divide d (H _ _ s) = unsafePerformIO $ do
h@(H _ _ d') <- cloneHistogram2D d
check "divide" $
withForeignPtr d' $ \dp ->
withForeignPtr s $ \sp -> histogram2d_div dp sp
return h
-- | multiplies the contents of the bins by a constant
scale :: Histogram2D -> Double -> Histogram2D
scale d s = unsafePerformIO $ do
h@(H _ _ d') <- cloneHistogram2D d
check "scale" $
withForeignPtr d' $ (\f -> histogram2d_scale f s)
return h
-- | adds a constant to the contents of the bins
shift :: Histogram2D -> Double -> Histogram2D
shift d s = unsafePerformIO $ do
h@(H _ _ d') <- cloneHistogram2D d
check "shift" $
withForeignPtr d' $ (\f -> histogram2d_shift f s)
return h
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_equal_bins_p" histogram2d_equal_bins :: Hist2DHandle -> Hist2DHandle -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_add" histogram2d_add :: Hist2DHandle -> Hist2DHandle -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_sub" histogram2d_sub :: Hist2DHandle -> Hist2DHandle -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_mul" histogram2d_mul :: Hist2DHandle -> Hist2DHandle -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_div" histogram2d_div :: Hist2DHandle -> Hist2DHandle -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_scale" histogram2d_scale :: Hist2DHandle -> Double -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_shift" histogram2d_shift :: Hist2DHandle -> Double -> IO CInt
-----------------------------------------------------------------------------
-- | write a histogram in the native binary format (may not be portable)
fwriteHistogram2D :: FilePath -> Histogram2D -> IO ()
fwriteHistogram2D fn (H _ _ h) = do
cn <- newCString fn
check "fwriteHistogram2d2D" $
withForeignPtr h $ histogram2d_fwrite cn
free cn
-- | read a histogram in the native binary format, number of bins must be known
freadHistogram2D :: FilePath -> Int -> Int -> IO Histogram2D
freadHistogram2D fn bx by = do
h <- histogram2d_new (fromIntegral bx) (fromIntegral by)
h' <- newForeignPtr histogram2d_free h
cn <- newCString fn
check "freadHistogram2d2D" $
withForeignPtr h' $ histogram2d_fread cn
return $ H bx by h'
-- | saves the histogram with the given formats (%f,%e,%g) for ranges and bins
-- each line comprises: xrange[i] xrange[i+1] xrange[j] xrange[j+1] bin(i,j)
fprintfHistogram2D :: FilePath -> String -> String -> Histogram2D -> IO ()
fprintfHistogram2D fn fr fb (H _ _ h) = do
cn <- newCString fn
cr <- newCString fr
cb <- newCString fb
check "fprintfHistogram2d2D" $
withForeignPtr h $ histogram2d_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
fscanfHistogram2D :: FilePath -> Int -> Int -> IO Histogram2D
fscanfHistogram2D fn bx by = do
h <- histogram2d_new (fromIntegral bx) (fromIntegral by)
h' <- newForeignPtr histogram2d_free h
cn <- newCString fn
check "fscanfHistogram2d2D" $
withForeignPtr h' $ histogram2d_fscanf cn
return $ H bx by h'
foreign import ccall "histogram-aux.h hist2d_fwrite" histogram2d_fwrite :: Ptr CChar -> Hist2DHandle -> IO CInt
foreign import ccall "histogram-aux.h hist2d_fread" histogram2d_fread :: Ptr CChar -> Hist2DHandle -> IO CInt
foreign import ccall "histogram-aux.h hist2d_fprintf" histogram2d_fprintf :: Ptr CChar -> Ptr CChar -> Ptr CChar -> Hist2DHandle -> IO CInt
foreign import ccall "histogram-aux.h hist2d_fscanf" histogram2d_fscanf :: Ptr CChar -> Hist2DHandle -> IO CInt
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_pdf_alloc" histogram2d_pdf_new :: CInt -> CInt -> IO PDFHandle
foreign import ccall "gsl-histogram2d.h &gsl_histogram2d_pdf_free" histogram2d_pdf_free :: FunPtr (PDFHandle -> IO ())
-----------------------------------------------------------------------------
-- | create a histogram PDF from a histogram
fromHistogram2D :: Histogram2D -> Histogram2DPDF
fromHistogram2D (H bx by h) = unsafePerformIO $ do
p <- histogram2d_pdf_new (fromIntegral bx) (fromIntegral by)
p' <- newForeignPtr histogram2d_pdf_free p
withForeignPtr p' $ \p'' ->
withForeignPtr h $ \h' -> do
histogram2d_pdf_init p'' h'
return $ P p'
-- | given a randomm from the uniform distribution [0,1], draw a random sample from the PDF
sample :: Histogram2DPDF -> Double
sample (P p) = unsafePerformIO $ withForeignPtr p $ \p' -> histogram2d_pdf_sample p'
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_pdf_init" histogram2d_pdf_init :: PDFHandle -> Hist2DHandle -> IO CInt
foreign import ccall "gsl-histogram2d.h gsl_histogram2d_pdf_sample" histogram2d_pdf_sample :: PDFHandle -> IO Double
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------