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 Foreign hiding(shift)
import Foreign.C.Types(CInt,CChar)
import Foreign.C.String(newCString)
import Prelude hiding(subtract,sum)
data Hist2D
type Hist2DHandle = Ptr Hist2D
data Histogram2D = H { hxdim :: !Int
, hydim :: !Int
, hist :: !(ForeignPtr Hist2D) }
data PDF
type PDFHandle = Ptr PDF
data Histogram2DPDF = P { pdf :: !(ForeignPtr PDF)}
instance Eq Histogram2D where
(==) = equalBins
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 ())
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'
fromLimitsIO :: Int -> Int
-> (Double,Double)
-> (Double,Double)
-> IO Histogram2D
fromLimitsIO nx ny (lx,ux) (uy,ly) = do
h <- histogram2d_new (fromIntegral nx) (fromIntegral ny)
h' <- newForeignPtr histogram2d_free h
check "set_ranges_uniform" $ 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
fromRanges :: Vector Double
-> Vector Double
-> [(Double,Double)]
-> Histogram2D
fromRanges rx ry d = unsafePerformIO $ do
h <- fromRangesIO rx ry
incrementListIO h d
return h
fromLimits :: Int -> Int
-> (Double,Double)
-> (Double,Double)
-> [(Double,Double)]
-> Histogram2D
fromLimits nx ny rx ry d = unsafePerformIO $ do
h <- fromLimitsIO nx ny rx ry
incrementListIO h d
return h
toMatrix :: Histogram2D -> (Vector Double,Vector Double,Matrix Double)
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
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
addList :: Histogram2D -> [(Double,Double)] -> Histogram2D
addList h xs = unsafePerformIO $ do
h' <- cloneHistogram2D h
incrementListIO h' xs
return h'
addVector :: Histogram2D -> Vector Double -> Vector Double -> Histogram2D
addVector h x y = unsafePerformIO $ do
h' <- cloneHistogram2D h
incrementVectorIO h' x y
return h'
addListWeighted :: Histogram2D -> [(Double,Double,Double)] -> Histogram2D
addListWeighted h xs = unsafePerformIO $ do
h' <- cloneHistogram2D h
accumulateListIO h' xs
return h'
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'
incrementIO :: Histogram2D -> Double -> Double -> IO ()
incrementIO (H _ _ h) x y = do
check "increment" $ withForeignPtr h (\f -> histogram2d_increment f x y)
return ()
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 ()
incrementListIO :: Histogram2D -> [(Double,Double)] -> IO ()
incrementListIO (H _ _ h) zs = withForeignPtr h (\f -> mapM_ (\(x,y) -> histogram2d_increment f x y) zs)
accumulateIO :: Histogram2D -> Double -> Double -> Double -> IO ()
accumulateIO (H _ _ h) x y w = do
check "accumulate" $ withForeignPtr h (\f -> histogram2d_accumulate f x y w)
return ()
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 ()
accumulateListIO :: Histogram2D -> [(Double,Double,Double)] -> IO ()
accumulateListIO (H _ _ h) zs = withForeignPtr h (\f -> mapM_ (\(x,y,w) -> histogram2d_accumulate f x y w) zs)
getBin :: Histogram2D -> Int -> Int -> Double
getBin (H _ _ h) bx by = unsafePerformIO $ do
withForeignPtr h (\f -> histogram2d_get f (fromIntegral bx) (fromIntegral by))
getXRange :: Histogram2D -> Int -> (Double,Double)
getXRange (H _ _ h) b = unsafePerformIO $ do
alloca $ \l ->
alloca $ \u -> do
check "get_xrange" $ withForeignPtr h (\f -> histogram2d_get_xrange f (fromIntegral b) l u)
l' <- peek l
u' <- peek u
return (l',u')
getYRange :: Histogram2D -> Int -> (Double,Double)
getYRange (H _ _ h) b = unsafePerformIO $ do
alloca $ \l ->
alloca $ \u -> do
check "get_yrange" $ withForeignPtr h (\f -> histogram2d_get_yrange f (fromIntegral b) l u)
l' <- peek l
u' <- peek u
return (l',u')
getXMax :: Histogram2D -> Double
getXMax (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmax
getXMin :: Histogram2D -> Double
getXMin (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmin
getXBins :: Histogram2D -> Int
getXBins (H bx _ _) = bx
getYMax :: Histogram2D -> Double
getYMax (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymax
getYMin :: Histogram2D -> Double
getYMin (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymin
getYBins :: Histogram2D -> Int
getYBins (H _ by _) = by
reset :: Histogram2D -> IO ()
reset (H _ _ h) = withForeignPtr h histogram2d_reset
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 :: 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
maxVal :: Histogram2D -> Double
maxVal (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_max_val
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')
minVal :: Histogram2D -> Double
minVal (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_min_val
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')
xmean :: Histogram2D -> Double
xmean (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xmean
ymean :: Histogram2D -> Double
ymean (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ymean
xstddev :: Histogram2D -> Double
xstddev (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_xsigma
ystddev :: Histogram2D -> Double
ystddev (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_ysigma
covariance :: Histogram2D -> Double
covariance (H _ _ h) = unsafePerformIO $ withForeignPtr h histogram2d_cov
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
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
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
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
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
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
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
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
fwriteHistogram2D :: FilePath -> Histogram2D -> IO ()
fwriteHistogram2D fn (H _ _ h) = do
cn <- newCString fn
check "fwriteHistogram2d2D" $
withForeignPtr h $ histogram2d_fwrite cn
free cn
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'
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 ()
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 ())
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
check "pdf_init" $ histogram2d_pdf_init p'' h'
return $ P p'
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