histogram-fill 0.1.0 → 0.2.0

raw patch · 8 files changed

+1062/−353 lines, 8 filesdep +vectordep −uvectorPVP ok

version bump matches the API change (PVP)

Dependencies added: vector

Dependencies removed: uvector

API changes (from Hackage documentation)

Files

• Data/Histogram.hs +92/−81
• Data/Histogram/Bin.hs +388/−46
• Data/Histogram/Bin/Extra.hs +76/−0
• Data/Histogram/Fill.hs +259/−97
• Data/Histogram/Generic.hs +191/−0
• Data/Histogram/Parse.hs +1/−3
• Data/Histogram/ST.hs +45/−121
• histogram-fill.cabal +10/−5

@@ -1,6 +1,8 @@+ {-# LANGUAGE GADTs #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-} -- | -- Module     : Data.Histogram -- Copyright  : Copyright (c) 2009, Alexey Khudyakov <alexey.skladnoy@gmail.com>@@ -11,119 +13,128 @@ -- Immutable histograms.   module Data.Histogram ( -- * Immutable histogram-                        Histogram(..)-                      , module Data.Histogram.Bin-                      , mapHist-                      , histBin-                      , histData-                      , underflows-                      , overflows-                      , outOfRange-                      , readHistogram-                      -- * Conversion-                      , asList-                      , asPairVector-                      , asVectorPairs-                      -- * Slicing-                      , sliceY-                      , sliceX-                      ) where+    -- * Data type+    Histogram+  , module Data.Histogram.Bin+  , histogram+  , histogramUO+    -- * Read histograms from string+  , readHistogram+  , readFileHistogram+    -- * Accessors+  , bins+  , histData+  , underflows+  , overflows+  , outOfRange+    -- ** Convert to other data types+  , asList+  , asVector+    -- * Slicing histograms+  , sliceX+  , sliceY+    -- * Modify histogram+  , histMap+  , histMapBin+  , histZip+  ) where -import Control.Arrow ((***))-import Control.Monad (ap)-import Data.Array.Vector-import Text.Read-import Text.ParserCombinators.ReadPrec (readPrec_to_S)+import qualified Data.Vector.Unboxed    as U+import Data.Vector.Unboxed (Unbox,Vector) +import qualified Data.Histogram.Generic as H import Data.Histogram.Bin-import Data.Histogram.Parse   -- | Immutable histogram. Histogram consists of binning algorithm, --   optional number of under and overflows, and data. -data Histogram bin a where-    Histogram :: (Bin bin, UA a) => -                 bin-              -> Maybe (a,a)-              -> UArr a-              -> Histogram bin a+type Histogram bin a = H.Histogram U.Vector bin a +-- | Create histogram from binning algorithm and vector with+-- data. Overflows are set to Nothing. +--+-- Number of bins and vector size must match.+histogram :: (Unbox a, Bin bin) => bin -> Vector a -> Histogram bin a+histogram = H.histogram -instance (Show a, Show (BinValue bin), Show bin) => Show (Histogram bin a) where-    show h@(Histogram bin uo _) = "# Histogram\n" ++ showUO uo ++ show bin ++-                                  (unlines $ map showT $ asList h)-        where-          showT (x,y) = show x ++ "\t" ++ show y-          showUO (Just (u,o)) = "# Underflows = " ++ show u ++ "\n" ++-                                "# Overflows  = " ++ show o ++ "\n"-          showUO Nothing      = "# Underflows = \n" ++-                                "# Overflows  = \n"+-- | Create histogram from binning algorithm and vector with data. +--+-- Number of bins and vector size must match.+histogramUO :: (Unbox a, Bin bin) => bin -> Maybe (a,a) -> Vector a -> Histogram bin a+histogramUO = H.histogramUO -histHeader :: (Read bin, Read a, Bin bin, UA a) => ReadPrec (UArr a -> Histogram bin a)-histHeader = do-  keyword "Histogram"-  u   <- maybeValue "Underflows"-  o   <- maybeValue "Overflows"-  bin <- readPrec-  return $ Histogram bin ((,) `fmap` u `ap` o) +----------------------------------------------------------------+-- Instances & reading histograms from strings +----------------------------------------------------------------++ -- | Convert String to histogram. Histogram do not have Read instance --   because of slowness of ReadP-readHistogram :: (Read bin, Read a, Bin bin, UA a) => String -> Histogram bin a-readHistogram str = -    let [(h,rest)] = readPrec_to_S histHeader 0 str -        xs = map last . filter (not . null) . map words . lines $ rest-    in h (toU $ map read xs)+readHistogram :: (Read bin, Read a, Bin bin, Unbox a) => String -> Histogram bin a+readHistogram = H.readHistogram --- | fmap lookalike. It's not possible to create Functor instance---   because of UA restriction.-mapHist :: UA b => (a -> b) -> Histogram bin a -> Histogram bin b-mapHist f (Histogram bin uo a) = Histogram bin (fmap (f *** f) uo) (mapU f a)+-- | Read histogram from file.+readFileHistogram :: (Read bin, Read a, Bin bin, Unbox a) => FilePath -> IO (Histogram bin a)+readFileHistogram = H.readFileHistogram +----------------------------------------------------------------+-- Accessors & conversion+----------------------------------------------------------------+ -- | Histogram bins-histBin :: Histogram bin a -> bin-histBin (Histogram bin _ _) = bin+bins :: Histogram bin a -> bin+bins = H.bins  -- | Histogram data as vector-histData :: Histogram bin a -> UArr a-histData (Histogram _ _ a) = a+histData :: Histogram bin a -> Vector a+histData = H.histData  -- | Number of underflows underflows :: Histogram bin a -> Maybe a-underflows (Histogram _ uo _) = fmap fst uo+underflows = H.underflows  -- | Number of overflows overflows :: Histogram bin a -> Maybe a-overflows (Histogram _ uo _) = fmap snd uo+overflows = H.overflows  -- | Underflows and overflows outOfRange :: Histogram bin a -> Maybe (a,a)-outOfRange (Histogram _ uo _) = uo+outOfRange = H.outOfRange  -- | Convert histogram to list.-asList :: Histogram bin a -> [(BinValue bin, a)]-asList (Histogram bin _ arr) = map (fromIndex bin) [0..] `zip` fromU arr+asList :: (Unbox a, Bin bin) => Histogram bin a -> [(BinValue bin, a)]+asList = H.asList --- | Convert to pair of vectors-asPairVector :: UA (BinValue bin) => Histogram bin a -> (UArr (BinValue bin), UArr a)-asPairVector (Histogram bin _ a) = (toU $ map (fromIndex bin) [0 .. nBins bin], a)+-- | Convert histogram to vector+asVector :: (Bin bin, Unbox a, Unbox (BinValue bin), Unbox (BinValue bin,a)) +         => Histogram bin a -> Vector (BinValue bin, a) +asVector = H.asVector --- | Convert to vector of pairs-asVectorPairs :: UA (BinValue bin) => Histogram bin a -> UArr ((BinValue bin) :*: a)-asVectorPairs h@(Histogram _ _ _) = uncurry zipU . asPairVector $ h+----------------------------------------------------------------+-- Modify histograms+---------------------------------------------------------------- +-- | fmap lookalike. It's not possible to create Functor instance+--   because of class restrictions+histMap :: (Unbox a, Unbox b) => (a -> b) -> Histogram bin a -> Histogram bin b+histMap = H.histMap++-- | Apply function to histogram bins. Function must not change number of bins.+--   If it does error is thrown.+histMapBin :: (Bin bin, Bin bin') => (bin -> bin') -> Histogram bin a -> Histogram bin' a+histMapBin = H.histMapBin++-- | Zip two histograms together. Bins of histograms must be equal+--   otherwise error will be called.+histZip :: (Bin bin, Eq bin, Unbox a, Unbox b, Unbox c) =>+           (a -> b -> c) -> Histogram bin a -> Histogram bin b -> Histogram bin c+histZip = H.histZip+            -- | Slice 2D histogram along Y axis. This function is fast because it does not require reallocations.-sliceY :: (Bin bX, Bin bY) => Histogram (Bin2D bX bY) a -> [(BinValue bY, Histogram bX a)]-sliceY (Histogram b@(Bin2D bX _) _ a) = map mkHist $ init [0, nBins bX .. nBins b]-    where-      mkHist i = ( snd $ fromIndex b i-                 , Histogram bX Nothing (sliceU a i (nBins bX)) )+sliceY :: (Unbox a, Bin bX, Bin bY) => Histogram (Bin2D bX bY) a -> [(BinValue bY, Histogram bX a)]+sliceY = H.sliceY  -- | Slice 2D histogram along X axis.-sliceX :: (Bin bX, Bin bY) => Histogram (Bin2D bX bY) a -> [(BinValue bX, Histogram bY a)]-sliceX (Histogram b@(Bin2D bX bY) _ a) = map mkHist $ init [0 .. nx]-    where-      nx = nBins bX-      n  = nBins b-      mkHist i = ( fst $ fromIndex b i-                 , Histogram bY Nothing (toU $ map (indexU a) [i,i+nx .. n-1]) )+sliceX :: (Unbox a, Bin bX, Bin bY) => Histogram (Bin2D bX bY) a -> [(BinValue bX, Histogram bY a)]+sliceX = H.sliceX

@@ -1,3 +1,4 @@+{-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE GADTs        #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE BangPatterns #-}@@ -11,26 +12,61 @@ -- Binning algorithms. This is mapping from set of interest to integer -- indices and approximate reverse.  -module Data.Histogram.Bin ( -- * Type class+module Data.Histogram.Bin ( -- * Type classes                             Bin(..)-                          -- * Integer bins+                          , Bin1D(..)+                          , Indexable(..)+                          , Indexable2D(..)+                          -- * Bin types+                          -- ** Integer bins                           , BinI(..)-                          -- * Floating point bins+                          , binI0+                          -- ** Integer bins with non-1 size+                          , BinInt+                          , binInt+                          -- ** Indexed bins +                          , BinIx(BinIx,unBinIx)+                          , binIx+                          -- ** Floating point bins                           , BinF                           , binF                           , binFn-                          -- * 2D bins+                          , binI2binF+                          , scaleBinF+                          -- *** Specialized for Double +                          , BinD+                          , binD+                          , binDn+                          , binI2binD+                          , scaleBinD+                          -- ** Log scale point+                          , LogBinD +                          , logBinD+                          -- ** 2D bins                           , Bin2D(..)                           , (><)+                          , nBins2D+                          , toIndex2D+                          , binX+                          , binY+                          , fmapBinX+                          , fmapBinY+                          -- ** 2D indexed bins+                          , BinIx2D (unBinIx2D)+                          , binIx2D                           ) where +import Control.Monad import Data.Histogram.Parse import Text.Read (Read(..)) ----- | Abstract binning algorithm. Following invariant is expected to hold: +import GHC.Float (double2Int)+----------------------------------------------------------------+-- | Abstract binning algorithm. It provides way to map some values+-- onto continous range of integer values starting from zero.  -- +-- Following invariant is expected to hold: +--  -- > toIndex . fromIndex == id --  -- Reverse is not nessearily true. @@ -39,47 +75,172 @@     type BinValue b     -- | Convert from value to index. No bound checking performed     toIndex :: b -> BinValue b -> Int-    {-# INLINE toIndex #-}     -- | Convert from index to value.      fromIndex :: b -> Int -> BinValue b +    -- | Check whether value in range.+    inRange :: b -> BinValue b -> Bool     -- | Total number of bins     nBins :: b -> Int +----------------------------------------------------------------+-- | One dimensional binning algorithm. It means that bin values have+-- some inherent ordering. For example all binning algorithms for real+-- numbers could be members or this type class whereas binning+-- algorithms for R^2 could not. +class Bin b => Bin1D b where+    -- | List of center of bins in ascending order.+    binsList :: b -> [BinValue b]+    -- | List of bins in ascending order.+    binsListRange :: b -> [(BinValue b, BinValue b)]  ----------------------------------------------------------------+-- | Indexable is value which could be converted to and from Int+-- without information loss.+--+-- Always true+--+-- > deindex . index = id+--+-- Only if Int is in range+--+-- > index . deindex = id+class Indexable a where+    -- | Convert value to index+    index :: a -> Int +    -- | Convert index to value+    deindex :: Int -> a++instance Indexable Int where+    index   = id+    deindex = id++----------------------------------------------------------------+-- | This type class is same as Indexable but for 2D values.+class Indexable2D a where+    -- | Convert value to index+    index2D :: a -> (Int,Int)+    -- | Convert index to value+    deindex2D :: (Int,Int) -> a++instance (Indexable a, Indexable b) => Indexable2D (a,b) where+    index2D   (x,y) = (index x,   index y)+    deindex2D (i,j) = (deindex i, deindex j)++---------------------------------------------------------------- -- Integer bin+----------------------------------------------------------------+-- | Simple binning algorithm which map continous range of bins onto+-- indices. Each number correcsponds to different bin+data BinI = BinI {-# UNPACK #-} !Int {-# UNPACK #-} !Int+            deriving Eq --- | Integer bins. This is inclusive interval [from,to]-data BinI = BinI !Int !Int+-- | Construct BinI with n bins. Indexing starts from 0+binI0 :: Int -> BinI+binI0 n = BinI 0 (n-1)  instance Bin BinI where     type BinValue BinI = Int     toIndex   !(BinI base _) !x = x - base+    {-# INLINE toIndex #-}     fromIndex !(BinI base _) !x = x + base+    inRange   !(BinI x y) i     = i>=x && i<=y+    {-# INLINE inRange #-}     nBins     !(BinI x y) = y - x + 1 +instance Bin1D BinI where+    binsList (BinI lo hi) = [lo .. hi]+    binsListRange b = zip (binsList b) (binsList b)+ instance Show BinI where     show (BinI lo hi) = unlines [ "# BinI"                                 , "# Low  = " ++ show lo                                 , "# High = " ++ show hi                                 ]- instance Read BinI where-    readPrec = do-      keyword "BinI"-      l <- value "Low"-      h <- value "High"-      return $ BinI l h+    readPrec = keyword "BinI" >> liftM2 BinI (value "Low") (value "High") +----------------------------------------------------------------+-- Another form of Integer bin+---------------------------------------------------------------- +-- | Integer bins with size which differ from 1.+data BinInt = BinInt +              {-# UNPACK #-} !Int -- Low bound+              {-# UNPACK #-} !Int -- Bin size+              {-# UNPACK #-} !Int -- Number of bins+              deriving Eq++-- | Construct BinInt.+binInt :: Int                   -- ^ Lower bound+       -> Int                   -- ^ Bin size+       -> Int                   -- ^ Upper bound+       -> BinInt+binInt lo n hi = BinInt lo n nb+  where+    nb = (hi-lo) `div` n ++instance Bin BinInt where+    type BinValue BinInt = Int+    toIndex   !(BinInt base sz _) !x = (x - base) `div` sz+    {-# INLINE toIndex #-}+    fromIndex !(BinInt base sz _) !x = x * sz + base+    inRange   !(BinInt base sz n) i  = i>=base && i<(base+n*sz)+    {-# INLINE inRange #-}+    nBins     !(BinInt _ _ n) = n++instance Show BinInt where+    show (BinInt base sz n) = +      unlines [ "# BinInt"+              , "# Base = " ++ show base+              , "# Step = " ++ show sz+              , "# Bins = " ++ show n+              ]++instance Read BinInt where+    readPrec = keyword "BinInt" >> liftM3 BinInt (value "Base") (value "Step") (value "Bins")+ ------------------------------------------------------------------- Floating point bin+-- Bins for indexables+---------------------------------------------------------------- +-- | Binning for indexable values+newtype BinIx i = BinIx { unBinIx :: BinI }+                  deriving Eq++-- | Construct indexed bin+binIx :: Indexable i => i -> i -> BinIx i+binIx lo hi = BinIx $ BinI (index lo) (index hi)++instance Indexable i => Bin (BinIx i) where+    type BinValue (BinIx i) = i+    toIndex   (BinIx b) x = toIndex b (index x)+    fromIndex (BinIx b) i = deindex (fromIndex b i)+    inRange   (BinIx b) x = inRange b (index x)+    nBins (BinIx b) = nBins b++instance Indexable i => Bin1D (BinIx i) where+    binsList (BinIx b) = map deindex (binsList b)+    binsListRange b    = let bins = binsList b in zip bins bins++instance (Show i, Indexable i) => Show (BinIx i) where+    show (BinIx (BinI lo hi)) = unlines [ "# BinIx"+                                        , "# Low  = " ++ show (deindex lo :: i)+                                        , "# High = " ++ show (deindex hi :: i)+                                        ]+instance (Read i, Indexable i) => Read (BinIx i) where+    readPrec = keyword "BinIx" >> liftM2 binIx (value "Low") (value "High")++----------------------------------------------------------------+-- Floating point bin+---------------------------------------------------------------- -- | Floaintg point bins with equal sizes. data BinF f where     BinF :: RealFrac f => !f -> !f -> !Int -> BinF f  --- | Create bins +instance Eq f => Eq (BinF f) where+    (BinF lo hi n) == (BinF lo' hi' n') = lo == lo'  && hi == hi' && n == n'+                                          +-- | Create bins. binF :: RealFrac f =>          f   -- ^ Lower bound of range      -> Int -- ^ Number of bins@@ -95,21 +256,37 @@       -> BinF f  binFn from step to = BinF from step (round $ (to - from) / step) +-- | Convert BinI to BinF+binI2binF :: RealFrac f => BinI -> BinF f+binI2binF b@(BinI i _) = BinF (fromIntegral i) 1 (nBins b)++-- | 'scaleBinF a b' scales BinF using linear transform 'a+b*x'+scaleBinF :: RealFrac f => f -> f -> BinF f -> BinF f+scaleBinF a b (BinF base step n) +    | b > 0     = BinF (a + b*base) (b*step) n+    | otherwise = error $ "scaleBinF: b must be positive (b = "++show b++")"+ instance Bin (BinF f) where     type BinValue (BinF f) = f      toIndex   !(BinF from step _) !x = floor $ (x-from) / step+    {-# INLINE toIndex #-}     fromIndex !(BinF from step _) !i = (step/2) + (fromIntegral i * step) + from +    inRange   !(BinF from step n) x  = x > from && x < from + step*fromIntegral n+    {-# INLINE inRange #-}     nBins     !(BinF _ _ n) = n-    {-# SPECIALIZE instance Bin (BinF Double) #-}-    {-# SPECIALIZE instance Bin (BinF Float) #-} +instance Bin1D (BinF f) where+    binsList b@(BinF _ _ n) = map (fromIndex b) [0..n-1]+    binsListRange b@(BinF _ step _) = map toPair (binsList b)+        where+          toPair x = (x - step/2, x + step/2)+ instance Show f => Show (BinF f) where     show (BinF base step n) = unlines [ "# BinF"-                                  , "# Base = " ++ show base-                                  , "# Step = " ++ show step-                                  , "# N    = " ++ show n-                                  ]-+                                      , "# Base = " ++ show base+                                      , "# Step = " ++ show step+                                      , "# N    = " ++ show n+                                      ] instance (Read f, RealFrac f) => Read (BinF f) where     readPrec = do       keyword "BinF"@@ -118,40 +295,173 @@       n    <- value "N"       return $ BinF base step n +----------------------------------------------------------------+-- Floating point bin /Specialized for Double+----------------------------------------------------------------+-- | Floaintg point bins with equal sizes. If you work with Doubles+-- this data type should be used instead of BinF.+data BinD = BinD {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Int +instance Eq BinD where+    (BinD lo hi n) == (BinD lo' hi' n') = lo == lo'  && hi == hi' && n == n'+                                          +-- | Create bins.+binD :: Double -- ^ Lower bound of range+     -> Int    -- ^ Number of bins+     -> Double -- ^ Upper bound of range+     -> BinD+binD from n to = BinD from ((to - from) / fromIntegral n) n++-- | Create bins. Note that actual upper bound can differ from specified.+binDn :: Double -- ^ Begin of range+      -> Double -- ^ Size of step+      -> Double -- ^ Approximation of end of range+      -> BinD+binDn from step to = BinD from step (round $ (to - from) / step)++-- | Convert BinI to BinF+binI2binD :: BinI -> BinD+binI2binD b@(BinI i _) = BinD (fromIntegral i) 1 (nBins b)++-- | 'scaleBinF a b' scales BinF using linear transform 'a+b*x'+scaleBinD :: Double -> Double -> BinD -> BinD+scaleBinD a b (BinD base step n) +    | b > 0     = BinD (a + b*base) (b*step) n+    | otherwise = error $ "scaleBinF: b must be positive (b = "++show b++")"++-- Fast variant of flooor+floorD :: Double -> Int+floorD x | x < 0     = double2Int x - 1+         | otherwise = double2Int x+{-# INLINE floorD #-}++instance Bin BinD where+    type BinValue BinD = Double+    toIndex   !(BinD from step _) !x = floorD $ (x-from) / step+    {-# INLINE toIndex #-}+    fromIndex !(BinD from step _) !i = (step/2) + (fromIntegral i * step) + from +    inRange   !(BinD from step n) x  = x > from && x < from + step*fromIntegral n+    {-# INLINE inRange #-}+    nBins     !(BinD _ _ n) = n++instance Bin1D BinD where+    binsList b@(BinD _ _ n) = map (fromIndex b) [0..n-1]+    binsListRange b@(BinD _ step _) = map toPair (binsList b)+        where+          toPair x = (x - step/2, x + step/2)++instance Show BinD where+    show (BinD base step n) = unlines [ "# BinD"+                                      , "# Base = " ++ show base+                                      , "# Step = " ++ show step+                                      , "# N    = " ++ show n+                                      ]+instance Read BinD where+    readPrec = do+      keyword "BinD"+      base <- value "Base"+      step <- value "Step"+      n    <- value "N"+      return $ BinD base step n++ ----------------------------------------------------------------+-- Log-scale bin+----------------------------------------------------------------+-- | Logarithmic scale bins.+data LogBinD = LogBinD+               Double -- Low border+               Double -- Hi border+               Double -- Increment ratio+               Int    -- Number of bins+               deriving Eq++-- | Create log-scale bins. +logBinD :: Double -> Int -> Double -> LogBinD+logBinD lo n hi = LogBinD lo hi ((hi/lo) ** (1 / fromIntegral n)) n++instance Bin LogBinD where+    type BinValue LogBinD = Double+    toIndex   !(LogBinD base _ step _) !x = floorD $ logBase step (x / base)+    {-# INLINE toIndex #-}+    fromIndex !(LogBinD base _ step _) !i = base * step ^ i+    inRange   !(LogBinD lo hi _ _) x  = x >= lo && x < hi+    {-# INLINE inRange #-}+    nBins     !(LogBinD _ _ _ n) = n++instance Show LogBinD where+    show (LogBinD lo hi step n) = +        unlines [ "# LogBinD"+                , "# Lo   = " ++ show lo+                , "# Hi   = " ++ show hi+                , "# Step = " ++ show step+                , "# N    = " ++ show n+                ]++---------------------------------------------------------------- -- 2D bin+---------------------------------------------------------------- --- | 2D bins. bin1 is binning along X axis and bin2 is one along Y axis. -data Bin2D bin1 bin2 = Bin2D bin1 bin2+-- | 2D bins. binX is binning along X axis and binY is one along Y axis. +data Bin2D binX binY = Bin2D !binX !binY+                       deriving Eq  -- | Alias for 'Bin2D'.-(><) :: bin1 -> bin2 -> Bin2D bin1 bin2+(><) :: binX -> binY -> Bin2D binX binY (><) = Bin2D -instance (Bin bin1, Bin bin2) => Bin (Bin2D bin1 bin2) where-    type BinValue (Bin2D bin1 bin2) = (BinValue bin1, BinValue bin2)+-- | Get binning algorithm along X axis+binX :: Bin2D bx by -> bx+binX !(Bin2D bx _) = bx -    toIndex   (Bin2D bx by) (x,y) -        | ix < 0 || ix >= rx || iy < 0 || iy >= ry = maxBound-        | otherwise                                = ix + iy*rx-        where-          ix = toIndex bx x-          iy = toIndex by y-          rx = nBins bx-          ry = nBins by+-- | Get binning algorithm along Y axis+binY :: Bin2D bx by -> by+binY !(Bin2D _ by) = by -    fromIndex (Bin2D bx by) i = let (iy,ix) = divMod i (nBins bx)-                                in  (fromIndex bx ix, fromIndex by iy)+instance (Bin binX, Bin binY) => Bin (Bin2D binX binY) where+    type BinValue (Bin2D binX binY) = (BinValue binX, BinValue binY)+    toIndex b@(Bin2D bx by) (x,y) +        | inRange b (x,y) = toIndex bx x + (toIndex by y)*(fromIntegral $ nBins bx)+        | otherwise       = maxBound+    {-# INLINE toIndex #-}+    fromIndex b@(Bin2D bx by) i = let (ix,iy) = toIndex2D b i+                                  in  (fromIndex bx ix, fromIndex by iy)+    inRange (Bin2D bx by) (x,y) = inRange bx x && inRange by y+    {-# INLINE inRange #-}+    nBins (Bin2D bx by) = (nBins bx) * (nBins by) -    nBins (Bin2D b1 b2) = (nBins b1) * (nBins b2)+toIndex2D :: (Bin binX, Bin binY) => Bin2D binX binY -> Int -> (Int,Int)+toIndex2D b i = let (iy,ix) = divMod i (nBins $ binX b) in (ix,iy) +-- | 2-dimensional size of binning algorithm+nBins2D :: (Bin bx, Bin by) => Bin2D bx by -> (Int,Int)+nBins2D (Bin2D bx by) = (nBins bx, nBins by)++-- | Apply function to X binning algorithm. If new binning algorithm+--   have different number of bins will fail.+fmapBinX :: (Bin bx, Bin bx') => (bx -> bx') -> Bin2D bx by -> Bin2D bx' by+fmapBinX f (Bin2D bx by) +    | nBins bx' /= nBins bx = error "fmapBinX: new binnig algorithm has different number of bins"+    | otherwise             = Bin2D bx' by+    where +      bx' = f bx++-- | Apply function to Y binning algorithm. If new binning algorithm+--   have different number of bins will fail.+fmapBinY ::(Bin by, Bin by') => (by -> by') -> Bin2D bx by -> Bin2D bx by'+fmapBinY f (Bin2D bx by)+    | nBins by' /= nBins by = error "fmapBinY: new binnig algorithm has different number of bins"+    | otherwise             = Bin2D bx by'+    where +      by' = f by+ instance (Show b1, Show b2) => Show (Bin2D b1 b2) where-    show (Bin2D b1 b2) = "# Bin2D\n" ++-                         "# X\n" ++ -                         show b1 ++-                         "# Y\n" ++-                         show b2+    show (Bin2D b1 b2) = concat [ "# Bin2D\n"+                                , "# X\n"+                                , show b1+                                , "# Y\n"+                                , show b2+                                ] instance (Read b1, Read b2) => Read (Bin2D b1 b2) where     readPrec = do       keyword "Bin2D"@@ -160,3 +470,35 @@       keyword "Y"       b2 <- readPrec       return $ Bin2D b1 b2+++----------------------------------------------------------------+-- Indexed 2D bins+----------------------------------------------------------------+-- | Binning for 2D indexable value+newtype BinIx2D i = BinIx2D {unBinIx2D :: (Bin2D BinI BinI) }++-- | Construct indexed bin+binIx2D :: Indexable2D i => i -> i -> BinIx2D i+binIx2D lo hi = let (ix,iy) = index2D lo+                    (jx,jy) = index2D hi+                in BinIx2D $ BinI ix jx >< BinI iy jy++instance Indexable2D i => Bin (BinIx2D i) where+    type BinValue (BinIx2D i) = i+    toIndex   (BinIx2D b) x = toIndex b (index2D x)+    fromIndex (BinIx2D b) i = deindex2D $ fromIndex b i+    inRange   (BinIx2D b) x = inRange b (index2D x)+    nBins     (BinIx2D b)   = nBins b++instance (Show i, Indexable2D i) => Show (BinIx2D i) where+    show (BinIx2D b) = unlines [ "# BinIx2D"+                               , "# Low  = " ++ show (deindex2D (fromIndex b 0            ) :: i)+                               , "# High = " ++ show (deindex2D (fromIndex b (nBins b - 1)) :: i)+                               ]+instance (Read i, Indexable2D i) => Read (BinIx2D i) where+    readPrec = do+      keyword "BinIx2D"+      l <- value "Low"+      h <- value "High"+      return $ binIx2D l h

@@ -0,0 +1,76 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+-- |+-- Module     : Data.Histogram.Bin+-- Copyright  : Copyright (c) 2010, Alexey Khudyakov <alexey.skladnoy@gmail.com>+-- License    : BSD3+-- Maintainer : Alexey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability  : experimental+--+-- Extra binning algorithms++module Data.Histogram.Bin.Extra ( BinPermute(permutedBin, permuteTo, permuteFrom)+                                , permuteBin+                                ) where++import Control.Applicative+import Control.Monad (forM_)+import qualified Data.Vector.Unboxed         as U+import qualified Data.Vector.Unboxed.Mutable as M+import Data.Vector.Unboxed ((!))+import Text.Read++import Data.Histogram.Bin+import Data.Histogram.Parse++-- | Direct permutation of indices. +data BinPermute b = BinPermute { permutedBin :: b+                               , permuteTo   :: U.Vector Int+                               , permuteFrom :: U.Vector Int+                               }+instance Bin b => Bin (BinPermute b) where+    type BinValue (BinPermute b) = BinValue b+    toIndex   (BinPermute b to _)   x = to ! toIndex b x+    fromIndex (BinPermute b _ from) i = fromIndex b (from ! i)+    inRange   (BinPermute b _ _) x = inRange b x+    nBins     (BinPermute b _ _) = nBins b++instance Show b => Show (BinPermute b) where+    show (BinPermute b to _) = unlines [ "# BinPermute"+                                       , "# Permutation = " ++ show (U.toList to)+                                       ] ++ show b++instance Read BinI => Read (BinPermute BinI) where+    readPrec = do keyword "BinPermute"+                  to   <- U.fromList <$> value "Permutation"+                  from <- case checkPermutation (invertPermutation to) of+                            Just v  -> return v+                            Nothing -> fail "Invalid permutation"+                  b  <- readPrec +                  return $ BinPermute b to from++-- Check whether this viable permutation+checkPermutation :: U.Vector Int -> Maybe (U.Vector Int)+checkPermutation v | U.any bad v = Nothing+                   | otherwise   = Just v+                   where+                     n     = U.length v+                     bad i = i < 0 || i >= n++-- Calculate inverse permutation                     +invertPermutation :: U.Vector Int -> U.Vector Int+invertPermutation v = U.create $ do a <- M.newWith n (-1)+                                    forM_ [0..n-1] (writeInvert a)+                                    return a+  where+    n = U.length v+    writeInvert a i | j >= 0 && j < n = M.write a j i+                    | otherwise       = return ()+                    where j = v ! i++-- | Constuct bin permutation from function.+permuteBin :: Bin b => (Int -> Int) -> b -> Maybe (BinPermute b)+permuteBin f b = BinPermute b <$> checkPermutation to <*> checkPermutation (invertPermutation to)+    where+      to   = U.map f $ U.enumFromN 0 (nBins b)

@@ -1,188 +1,350 @@ {-# LANGUAGE GADTs        #-} {-# LANGUAGE Rank2Types   #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-} -- | -- Module     : Data.Histogram.Fill -- Copyright  : Copyright (c) 2009, Alexey Khudyakov <alexey.skladnoy@gmail.com> -- License    : BSD3 -- Maintainer : Alexey Khudyakov <alexey.skladnoy@gmail.com> -- Stability  : experimental--- +-- -- Module with algorithms for histogram filling. This is pure wrapper -- around stateful histograms.--- -module Data.Histogram.Fill ( -- * Type classes & wrappers-                             HBuilderCl(..)+--+module Data.Histogram.Fill ( -- * Type classes+                             HistBuilder(..)+                             -- * Histogram builders+                             -- ** Stateful+                           , HBuilderST+                           , feedOne+                           , freezeHBuilderST+                           , joinHBuilderST+                           , joinHBuilderSTList+                           , treeHBuilderST+                             -- ** IO based+                           , HBuilderIO+                           , feedOneIO+                           , freezeHBuilderIO+                           , joinHBuilderIO+                           , joinHBuilderIOList+                           , treeHBuilderIO+                             -- ** Stateless                            , HBuilder-                           , builderList-                           , builderListWrap--                           -- * Fill routines-                           , createHistograms--                           -- * Histogram constructors +                           , joinHBuilder+                           , joinHBuilderList+                           , treeHBuilder+                             -- ** Conversion between builders+                           , toBuilderST+                           , toBuilderIO+                           , builderSTtoIO+                           -- * Fill histograms+                           , fillBuilder+                           -- * Histogram constructors                            , module Data.Histogram.Bin-                           , mkHist+                           -- ** Fixed weigth histograms                            , mkHist1-                           , mkHistWgh+                           , mkHist+                           , mkHistMaybe+                           -- ** Weighted histograms                            , mkHistWgh1-                           , mkHistMonoid+                           , mkHistWgh+                           , mkHistWghMaybe+                           -- ** Histograms with monoidal bins                            , mkHistMonoid1+                           , mkHistMonoid+                           , mkHistMonoidMaybe+                           -- * Auxillary functions                            , forceInt                            , forceDouble                            , forceFloat-                           -- * Internals-                           , HistBuilder                            ) where -import Control.Monad.ST (ST)-import Data.Monoid      (Monoid, mempty)+import Control.Applicative ((<$>))+import Control.Monad       (when)+import Control.Monad.ST  -import Data.Array.Vector+import Data.Monoid         (Monoid, mempty)+import Data.Vector.Unboxed (Unbox)+ import Data.Histogram import Data.Histogram.Bin import Data.Histogram.ST  -------------------------------------------------------------------- | Create and fill histogram(s).-createHistograms :: Monoid b =>-                    HBuilder a b -- ^ Instructions how to fill histograms-                 -> [a]          -- ^ List of data to fill histogram with-                 -> b            -- ^ Result-createHistograms h xs = fillHistograms (runBuilder h) xs-+-- Type class ----------------------------------------------------------------  -- | Histogram builder typeclass. Instance of this class contain --   instructions how to build histograms.-class HBuilderCl h where +class HistBuilder h where     -- | Convert input type of histogram from a to a'-    modifyIn  :: (a' -> a) -> h a b -> h a' b -    -- | Convert output of histogram +    modifyIn  :: (a' -> a) -> h a b -> h a' b+    -- | Make input function accept value only if it's Just a.+    modifyMaybe :: h a b -> h (Maybe a) b+    -- | Add cut to histogram. Only put value histogram if condition is true.+    addCut    :: (a -> Bool) -> h a b -> h a b+    -- | Convert output of histogram     modifyOut :: (b -> b') -> h a b -> h a  b'-    -- | Create stateful histogram from instructions. Histograms could-    --   be filled either in the ST monad or with createHistograms-    runBuilder :: h a b -> ST s (Accum s a b) - ----------------------------------------------------------------+-- ST based builder+---------------------------------------------------------------- --- | Abstract histogram builder. All real builders should be wrapper---   in this type-data HBuilder a b where-    MkHBuilder :: HBuilderCl h => h a b -> HBuilder a b+-- | Stateful histogram builder.+data HBuilderST s a b = HBuilderST { hbInput  :: a -> ST s ()+                                   , hbOutput :: ST s b+                                   } -instance HBuilderCl HBuilder where -    modifyIn  f (MkHBuilder h) = MkHBuilder $ modifyIn f h-    modifyOut g (MkHBuilder h) = MkHBuilder $ modifyOut g h -    runBuilder  (MkHBuilder h) = runBuilder h+instance HistBuilder (HBuilderST s) where+    modifyIn  f h = h { hbInput  = hbInput h . f }+    addCut    f h = h { hbInput  = \x -> when (f x) (hbInput h x) }+    modifyMaybe h = h { hbInput  = modified } +        where modified (Just x) = hbInput h x+              modified Nothing  = return ()+    modifyOut f h = h { hbOutput = f `fmap` hbOutput h } +instance Functor (HBuilderST s a) where+    fmap = modifyOut +-- | Put one value into histogram+feedOne :: HBuilderST s a b -> a -> ST s ()+feedOne = hbInput++-- | Create stateful histogram from instructions. Histograms could+--   be filled either in the ST monad or with createHistograms+freezeHBuilderST :: HBuilderST s a b -> ST s b+freezeHBuilderST = hbOutput+++-- | Join list of builders into one builder+joinHBuilderST :: [HBuilderST s a b] -> HBuilderST s a [b]+joinHBuilderST hs = HBuilderST { hbInput  = \x -> mapM_ (flip hbInput x) hs+                               , hbOutput = mapM hbOutput hs+                               }++-- | Join list of builders into one builders+joinHBuilderSTList :: [HBuilderST s a [b]] -> HBuilderST s a [b]+joinHBuilderSTList = fmap concat . joinHBuilderST++treeHBuilderST :: [HBuilderST s a b -> HBuilderST s a' b'] -> HBuilderST s a b -> HBuilderST s a' [b']+treeHBuilderST fs h = joinHBuilderST $ map ($ h) fs+ ----------------------------------------------------------------+-- IO based+---------------------------------------------------------------- --- List of histograms. -newtype HBuilderList a b = HBuilderList [HBuilder a b]+-- | Stateful histogram builder.+data HBuilderIO a b = HBuilderIO { hbInputIO  :: a -> IO ()+                                 , hbOutputIO :: IO b+                                 } --- | Wrap list of histogram builders into HBuilder.-builderList :: [HBuilder a b] -> HBuilder a [b]-builderList = MkHBuilder . modifyOut (:[]) . HBuilderList+instance HistBuilder (HBuilderIO) where+    modifyIn  f h = h { hbInputIO  = hbInputIO h . f }+    addCut    f h = h { hbInputIO  = \x -> when (f x) (hbInputIO h x) }+    modifyMaybe h = h { hbInputIO  = modified } +        where modified (Just x) = hbInputIO h x+              modified Nothing  = return ()+    modifyOut f h = h { hbOutputIO = f `fmap` hbOutputIO h } --- | Wrap list of histogram builders into HBuilder and do not change return type.-builderListWrap :: [HBuilder a b] -> HBuilder a b-builderListWrap = MkHBuilder . HBuilderList+instance Functor (HBuilderIO a) where+    fmap = modifyOut -instance HBuilderCl HBuilderList where-    modifyIn  f (HBuilderList l) = HBuilderList $ map (modifyIn f) l-    modifyOut g (HBuilderList l) = HBuilderList $ map (modifyOut g) l-    runBuilder (HBuilderList l)  = accumList $ map runBuilder l+-- | Put one value into histogram+feedOneIO :: HBuilderIO a b -> a -> IO ()+feedOneIO = hbInputIO +-- | Create stateful histogram from instructions. Histograms could+--   be filled either in the ST monad or with createHistograms+freezeHBuilderIO :: HBuilderIO a b -> IO b+freezeHBuilderIO = hbOutputIO++-- | Join list of builders into one builder+joinHBuilderIO :: [HBuilderIO a b] -> HBuilderIO a [b]+joinHBuilderIO hs = HBuilderIO { hbInputIO  = \x -> mapM_ (flip hbInputIO x) hs+                               , hbOutputIO = mapM hbOutputIO hs+                               }++-- | Join list of builders into one builders+joinHBuilderIOList :: [HBuilderIO a [b]] -> HBuilderIO a [b]+joinHBuilderIOList = fmap concat . joinHBuilderIO++treeHBuilderIO :: [HBuilderIO a b -> HBuilderIO a' b'] -> HBuilderIO a b -> HBuilderIO a' [b']+treeHBuilderIO fs h = joinHBuilderIO $ map ($ h) fs+ ----------------------------------------------------------------+-- Stateless +---------------------------------------------------------------- --- | Generic histogram builder. -data HistBuilder a b where-    HistBuilder :: (Bin bin, UA val) =>-                   bin                                                -- Bin type-                -> val                                                -- Zero element-                -> (forall s . a -> HistogramST s bin val -> ST s ()) -- Input function-                -> (Histogram bin val -> b)                           -- Output function-                -> HistBuilder a b+-- | Stateless histogram builder+newtype HBuilder a b = HBuilder { toBuilderST :: (forall s . ST s (HBuilderST s a b)) } -instance HBuilderCl HistBuilder where-    modifyIn  f (HistBuilder bin z inp out) = HistBuilder bin z (inp . f) out-    modifyOut g (HistBuilder bin z inp out) = HistBuilder bin z  inp (g . out)-    runBuilder  (HistBuilder bin z inp out) = accumHist inp out =<< newHistogramST z bin+instance HistBuilder (HBuilder) where+    modifyIn  f (HBuilder h) = HBuilder (modifyIn  f <$> h)+    addCut    f (HBuilder h) = HBuilder (addCut    f <$> h)+    modifyMaybe (HBuilder h) = HBuilder (modifyMaybe <$> h)+    modifyOut f (HBuilder h) = HBuilder (modifyOut f <$> h) +instance Functor (HBuilder a) where+    fmap = modifyOut +-- | Join list of builders+joinHBuilder :: [HBuilder a b] -> HBuilder a [b]+joinHBuilder hs = HBuilder (joinHBuilderST <$> mapM toBuilderST hs) +-- | Join list of builders+joinHBuilderList :: [HBuilder a [b]] -> HBuilder a [b]+joinHBuilderList = modifyOut concat . joinHBuilder++treeHBuilder :: [HBuilder a b -> HBuilder a' b'] -> HBuilder a b -> HBuilder a' [b']+treeHBuilder fs h = joinHBuilder $ map ($ h) fs+ ------------------------------------------------------------------- Histogram constructors +-- Conversions ---------------------------------------------------------------- --- | Function used to restrict type of histrogram.-forceInt :: Histogram bin Int -> Histogram bin Int-forceInt = id+-- | Convert ST base builder to IO based one+builderSTtoIO :: HBuilderST RealWorld a b -> HBuilderIO a b+builderSTtoIO (HBuilderST i o) = HBuilderIO (stToIO . i) (stToIO o) --- | Function used to restrict type of histrogram.-forceDouble :: Histogram bin Double -> Histogram bin Double-forceDouble = id+-- | Convert stateless builder to IO based one+toBuilderIO :: HBuilder a b -> IO (HBuilderIO a b)+toBuilderIO h = builderSTtoIO `fmap` stToIO (toBuilderST h) --- | Function used to restrict type of histrogram.-forceFloat :: Histogram bin Float -> Histogram bin Float-forceFloat = id+----------------------------------------------------------------+-- Actual filling of histograms+---------------------------------------------------------------- +fillBuilder :: HBuilder a b -> [a] -> b+fillBuilder hb xs = +    runST $ do h <- toBuilderST hb+               mapM_ (feedOne h) xs+               freezeHBuilderST h+  +----------------------------------------------------------------+-- Histogram constructors+----------------------------------------------------------------+ -- | Create histogram builder which take single item as input. Each---   item has weight 1. To set type of bin 'force*' function could be used.-mkHist1 :: (Bin bin, UA val, Num val) =>+--   item has weight 1.+mkHist1 :: (Bin bin, Unbox val, Num val) =>            bin                      -- ^ Bin information         -> (Histogram bin val -> b) -- ^ Output function          -> (a -> BinValue bin)      -- ^ Input function         -> HBuilder a b-mkHist1 bin out inp = MkHBuilder $ HistBuilder bin 0 (flip fillOne . inp) out+mkHist1 bin out inp = HBuilder $ do+  acc <- newMHistogram 0 bin+  return $ HBuilderST { hbInput  = fillOne acc . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }  -- | Create histogram builder which take many items as input. Each---   item has weight 1. To set type of bin 'force*' function could be---   used.-mkHist :: (Bin bin, UA val, Num val) =>+--   item has weight 1.+mkHist :: (Bin bin, Unbox val, Num val) =>           bin                      -- ^ Bin information        -> (Histogram bin val -> b) -- ^ Output function        -> (a -> [BinValue bin])    -- ^ Input function         -> HBuilder a b-mkHist bin out inp = MkHBuilder $ HistBuilder bin 0 fill out-    where-      fill a h = mapM_ (fillOne h) $ inp a+mkHist bin out inp = HBuilder $ do+  acc <- newMHistogram 0 bin+  return $ HBuilderST { hbInput  = mapM_ (fillOne acc) . inp+                      , hbOutput = fmap out (freezeHist acc)+                      } +-- | Create histogram builder which at most one item as input. Each+--   item has weight 1. +mkHistMaybe :: (Bin bin, Unbox val, Num val) =>+          bin                         -- ^ Bin information+       -> (Histogram bin val -> b)    -- ^ Output function+       -> (a -> Maybe (BinValue bin)) -- ^ Input function +       -> HBuilder a b+mkHistMaybe bin out inp = HBuilder $ do+  acc <- newMHistogram 0 bin+  return $ HBuilderST { hbInput  = maybe (return ()) (fillOne acc) . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }+ -- | Create histogram with weighted bin. Takes one item at time. -mkHistWgh1 :: (Bin bin, UA val, Num val) =>+mkHistWgh1 :: (Bin bin, Unbox val, Num val) =>               bin                        -- ^ Bin information           -> (Histogram bin val -> b)    -- ^ Output function           -> (a -> (BinValue bin, val))  -- ^ Input function           -> HBuilder a b-mkHistWgh1 bin out inp = MkHBuilder $ HistBuilder bin 0 (flip fillOneW . inp) out+mkHistWgh1 bin out inp = HBuilder $ do+  acc <- newMHistogram 0 bin+  return $ HBuilderST { hbInput  = fillOneW acc . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }  -- | Create histogram with weighted bin. Takes many items at time.-mkHistWgh :: (Bin bin, UA val, Num val) => +mkHistWgh :: (Bin bin, Unbox val, Num val) =>               bin                          -- ^ Bin information           -> (Histogram bin val  -> b)    -- ^ Output function           -> (a -> [(BinValue bin, val)]) -- ^ Input function           -> HBuilder a b-mkHistWgh bin out inp = MkHBuilder $ HistBuilder bin 0 fill out-    where-      fill a h = mapM_ (fillOneW h) $ inp a+mkHistWgh bin out inp = HBuilder $ do+  acc <- newMHistogram 0 bin+  return $ HBuilderST { hbInput  = mapM_ (fillOneW acc) . inp+                      , hbOutput = fmap out (freezeHist acc)+                      } +-- | Create histogram with weighted bin. Takes many items at time.+mkHistWghMaybe :: (Bin bin, Unbox val, Num val) => +                  bin                              -- ^ Bin information+               -> (Histogram bin val  -> b)        -- ^ Output function+               -> (a -> Maybe (BinValue bin, val)) -- ^ Input function+               -> HBuilder a b+mkHistWghMaybe bin out inp = HBuilder $ do+  acc <- newMHistogram 0 bin+  return $ HBuilderST { hbInput  = maybe (return ()) (fillOneW acc) . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }+ -- | Create histogram with monoidal bins-mkHistMonoid1 :: (Bin bin, UA val, Monoid val) =>+mkHistMonoid1 :: (Bin bin, Unbox val, Monoid val) =>               bin                         -- ^ Bin information           -> (Histogram bin val -> b)     -- ^ Output function           -> (a -> (BinValue bin, val))   -- ^ Input function           -> HBuilder a b-mkHistMonoid1 bin out inp = MkHBuilder $ HistBuilder bin mempty (flip fillMonoid . inp) out+mkHistMonoid1 bin out inp = HBuilder $ do+  acc <- newMHistogram mempty bin+  return $ HBuilderST { hbInput  = fillMonoid acc . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }  -- | Create histogram with monoidal bins. Takes many items at time.-mkHistMonoid :: (Bin bin, UA val, Monoid val) =>+mkHistMonoid :: (Bin bin, Unbox val, Monoid val) =>               bin                         -- ^ Bin information           -> (Histogram bin val -> b)     -- ^ Output function           -> (a -> [(BinValue bin, val)]) -- ^ Input function           -> HBuilder a b-mkHistMonoid bin out inp = MkHBuilder $ HistBuilder bin mempty fill out-    where-      fill a h = mapM_ (fillMonoid h) $ inp a+mkHistMonoid bin out inp = HBuilder $ do+  acc <- newMHistogram mempty bin+  return $ HBuilderST { hbInput  = mapM_ (fillMonoid acc) . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }++-- | Create histogram with monoidal bins+mkHistMonoidMaybe :: (Bin bin, Unbox val, Monoid val) =>+                     bin                              -- ^ Bin information+                  -> (Histogram bin val -> b)         -- ^ Output function+                  -> (a -> Maybe (BinValue bin, val)) -- ^ Input function+                  -> HBuilder a b+mkHistMonoidMaybe bin out inp = HBuilder $ do+  acc <- newMHistogram mempty bin+  return $ HBuilderST { hbInput  = maybe (return ()) (fillMonoid acc) . inp+                      , hbOutput = fmap out (freezeHist acc)+                      }++----------------------------------------------------------------++-- | Function used to restrict type of histrogram.+forceInt :: Histogram bin Int -> Histogram bin Int+forceInt = id++-- | Function used to restrict type of histrogram.+forceDouble :: Histogram bin Double -> Histogram bin Double+forceDouble = id++-- | Function used to restrict type of histrogram.+forceFloat :: Histogram bin Float -> Histogram bin Float+forceFloat = id

@@ -0,0 +1,191 @@+{-# LANGUAGE FlexibleContexts #-}+-- |+-- Module     : Data.Histogram+-- Copyright  : Copyright (c) 2009, Alexey Khudyakov <alexey.skladnoy@gmail.com>+-- License    : BSD3+-- Maintainer : Alexey Khudyakov <alexey.skladnoy@gmail.com>+-- Stability  : experimental+-- +-- Generic immutable histograms. +module Data.Histogram.Generic ( +    -- * Data type+    Histogram+  , module Data.Histogram.Bin+  , histogram+  , histogramUO+    -- * Read histograms from string+  , readHistogram+  , readFileHistogram+    -- * Accessors+  , bins+  , histData+  , underflows+  , overflows+  , outOfRange+    -- ** Convert to other data types+  , asList+  , asVector+    -- * Slicing histograms+  , sliceX+  , sliceY+    -- * Modify histogram+  , histMap+  , histMapBin+  , histZip+  ) where++import Control.Applicative ((<$>),(<*>))+import Control.Arrow       ((***))+import Control.Monad       (ap, forM_)+import Control.Monad.ST    (runST)++import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Generic as G+import Data.Vector.Generic  (Vector)+import Text.Read++import Data.Histogram.Bin+import Data.Histogram.Parse++----------------------------------------------------------------+-- Data type and smart constructors+----------------------------------------------------------------++-- | Immutable histogram. Histogram consists of binning algorithm,+--   optional number of under and overflows, and data. +data Histogram v bin a = Histogram bin (Maybe (a,a)) (v a)+                         deriving Eq++-- | Create histogram from binning algorithm and vector with+-- data. Overflows are set to Nothing. +--+-- Number of bins and vector size must match.+histogram :: (Vector v a, Bin bin) => bin -> v a -> Histogram v bin a+histogram b v | nBins b == G.length v = Histogram b Nothing v+              | otherwise             = error "histogram: number of bins and vector size doesn't match"+++-- | Create histogram from binning algorithm and vector with data. +--+-- Number of bins and vector size must match.+histogramUO :: (Vector v a, Bin bin) => bin -> Maybe (a,a) -> v a -> Histogram v bin a+histogramUO b uo v | nBins b == G.length v = Histogram b uo v+                   | otherwise             = error "histogram: number of bins and vector size doesn't match"+++----------------------------------------------------------------+-- Instances & reading histograms from strings +----------------------------------------------------------------++instance (Show a, Show (BinValue bin), Show bin, Bin bin, Vector v a) => Show (Histogram v bin a) where+    show h@(Histogram bin uo _) = "# Histogram\n" ++ showUO uo ++ show bin +++                                  (unlines $ map showT $ asList h)+        where+          showT (x,y) = show x ++ "\t" ++ show y+          showUO (Just (u,o)) = "# Underflows = " ++ show u ++ "\n" +++                                "# Overflows  = " ++ show o ++ "\n"+          showUO Nothing      = "# Underflows = \n" +++                                "# Overflows  = \n"++-- Parse histogram header+histHeader :: (Read bin, Read a, Bin bin, Vector v a) => ReadPrec (v a -> Histogram v bin a)+histHeader = do+  keyword "Histogram"+  u   <- maybeValue "Underflows"+  o   <- maybeValue "Overflows"+  bin <- readPrec+  return $ Histogram bin ((,) `fmap` u `ap` o)++-- | Convert String to histogram. Histogram do not have Read instance+--   because of slowness of ReadP+readHistogram :: (Read bin, Read a, Bin bin, Vector v a) => String -> Histogram v bin a+readHistogram str = +    let (h,rest) = case readPrec_to_S histHeader 0 str of+                     [x] -> x+                     _   -> error "Cannot parse histogram header"+        xs = map (unwords . tail) . filter (not . null) . map words . lines $ rest+    in h (G.fromList $ map read xs)++-- | Read histogram from file.+readFileHistogram :: (Read bin, Read a, Bin bin, Vector v a) => FilePath -> IO (Histogram v bin a)+readFileHistogram fname = readHistogram `fmap` readFile fname++----------------------------------------------------------------+-- Accessors & conversion+----------------------------------------------------------------++-- | Histogram bins+bins :: Histogram v bin a -> bin+bins (Histogram bin _ _) = bin++-- | Histogram data as vector+histData :: Histogram v bin a -> v a+histData (Histogram _ _ a) = a++-- | Number of underflows+underflows :: Histogram v bin a -> Maybe a+underflows (Histogram _ uo _) = fst <$> uo++-- | Number of overflows+overflows :: Histogram v bin a -> Maybe a+overflows (Histogram _ uo _) = snd <$> uo++-- | Underflows and overflows+outOfRange :: Histogram v bin a -> Maybe (a,a)+outOfRange (Histogram _ uo _) = uo++-- | Convert histogram to list.+asList :: (Vector v a, Bin bin) => Histogram v bin a -> [(BinValue bin, a)]+asList (Histogram bin _ arr) = map (fromIndex bin) [0..] `zip` G.toList arr++-- | Convert histogram to vector+asVector :: (Bin bin, Vector v a, Vector v (BinValue bin), Vector v (BinValue bin,a)) +         => Histogram v bin a -> v (BinValue bin, a) +asVector (Histogram bin _ arr) = G.zip (G.generate (nBins bin) (fromIndex bin) ) arr++----------------------------------------------------------------+-- Modify histograms+----------------------------------------------------------------++-- | fmap lookalike. It's not possible to create Functor instance+--   because of class restrictions+histMap :: (Vector v a, Vector v b) => (a -> b) -> Histogram v bin a -> Histogram v bin b+histMap f (Histogram bin uo a) = Histogram bin (fmap (f *** f) uo) (G.map f a)++-- | Apply function to histogram bins. Function must not change number of bins.+--   If it does error is thrown.+histMapBin :: (Bin bin, Bin bin') => (bin -> bin') -> Histogram v bin a -> Histogram v bin' a+histMapBin f (Histogram bin uo a)+    | nBins bin == nBins bin' = Histogram (f bin) uo a+    | otherwise               = error "Number of bins doesn't match"+    where+      bin' = bin++-- | Zip two histograms together. Bins of histograms must be equal+--   otherwise error will be called.+histZip :: (Bin bin, Eq bin, Vector v a, Vector v b, Vector v c) =>+           (a -> b -> c) -> Histogram v bin a -> Histogram v bin b -> Histogram v bin c+histZip f (Histogram bin uo v) (Histogram bin' uo' v')+    | bin /= bin' = error "histZip: bins are different"+    | otherwise   = Histogram bin (f2 <$> uo <*> uo') (G.zipWith f v v')+      where+        f2 (x,x') (y,y') = (f x y, f x' y')+           +-- | Slice 2D histogram along Y axis. This function is fast because it does not require reallocations.+sliceY :: (Vector v a, Bin bX, Bin bY) => Histogram v (Bin2D bX bY) a -> [(BinValue bY, Histogram v bX a)]+sliceY (Histogram b _ a) = map mkSlice [0 .. ny-1]+    where+      (nx, ny) = nBins2D b+      mkSlice i = ( fromIndex (binY b) i+                  , Histogram (binX b) Nothing (G.slice nx (nx*i) a) )++-- | Slice 2D histogram along X axis.+sliceX :: (Vector v a, Bin bX, Bin bY) => Histogram v (Bin2D bX bY) a -> [(BinValue bX, Histogram v bY a)]+sliceX (Histogram b _ a) = map mkSlice [0 .. nx-1]+    where+      (nx, ny)  = nBins2D b+      mkSlice i = ( fromIndex (binX b) i+                  , Histogram (binY b) Nothing (mkArray i))+      mkArray x = runST $ do arr <- M.new ny+                             forM_ [0 .. ny-1] $ \y -> M.write arr y (a G.! (y*nx + x))+                             G.unsafeFreeze arr

@@ -7,7 +7,6 @@  import Text.Read import Text.ParserCombinators.ReadP    (ReadP, many, satisfy, char, string)-import Text.ParserCombinators.ReadPrec  -- Whitespaces ws :: ReadP String@@ -28,8 +27,7 @@  getVal :: Read a => ReadPrec a getVal = do x <- readPrec-            lift eol -            return x+            lift eol >> return x  -- Key value pair value :: Read a => String -> ReadPrec a

@@ -11,160 +11,84 @@ -- Mutable histograms.  module Data.Histogram.ST ( -- * Mutable histograms-                           HistogramST(..)-                         , newHistogramST+                           MHistogram+                         , newMHistogram                          , fillOne                          , fillOneW                          , fillMonoid                          , freezeHist--                         -- * Accumulators-                         , Accumulator(..)-                         , Accum(Accum)--                         , accumList-                         , accumHist--                         , fillHistograms                          ) where   import Control.Monad.ST -import Data.Array.Vector import Data.Monoid+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as MU+import qualified Data.Vector.Generic as G  import Data.Histogram-import Data.Histogram.Bin - ---------------------------------------------------------------- -- Mutable histograms ----------------------------------------------------------------  -- | Mutable histogram.-data HistogramST s bin a where-    HistogramST :: (Bin bin, UA a) => -                   bin-                -> MUArr a s -- Over/underflows-                -> MUArr a s -- Data-                -> HistogramST s bin a+data MHistogram s bin a where+    MHistogram :: (Bin bin, MU.Unbox a) => +                  bin            -- Binning+               -> MU.MVector s a -- Over/underflows+               -> MU.MVector s a -- Data+               -> MHistogram s bin a  -- | Create new mutable histogram. All bins are set to zero element as --   passed to function.-newHistogramST :: (Bin bin, UA a) => a -> bin -> ST s (HistogramST s bin a)-newHistogramST zero bin = do-  uo <- newMU 2-  writeMU uo 0 zero >> writeMU uo 1 zero-  a <- newMU (nBins bin)-  mapM_ (\i -> writeMU a i zero) [0 .. (lengthMU a) - 1]-  return $ HistogramST bin uo a+newMHistogram :: (Bin bin, U.Unbox a) => a -> bin -> ST s (MHistogram s bin a)+newMHistogram zero bin = do+  uo <- MU.newWith 2 zero+  a  <- MU.newWith (nBins bin) zero+  return $ MHistogram bin uo a+{-# INLINE newMHistogram #-}  -- | Put one value into histogram-fillOne :: Num a => HistogramST s bin a -> BinValue bin -> ST s ()-fillOne (HistogramST bin uo arr) x-    | i < 0             = writeMU uo  0 . (+1)  =<< readMU uo 0-    | i >= lengthMU arr = writeMU uo  1 . (+1)  =<< readMU uo 1-    | otherwise         = writeMU arr i . (+1)  =<< readMU arr i+fillOne :: Num a => MHistogram s bin a -> BinValue bin -> ST s ()+fillOne (MHistogram bin uo arr) x+    | i < 0              = MU.unsafeWrite uo  0 . (+1)  =<< MU.unsafeRead uo 0+    | i >= MU.length arr = MU.unsafeWrite uo  1 . (+1)  =<< MU.unsafeRead uo 1+    | otherwise          = MU.unsafeWrite arr i . (+1)  =<< MU.unsafeRead arr i     where       i = toIndex bin x+{-# INLINE fillOne #-}  -- | Put one value into histogram with weight-fillOneW :: Num a => HistogramST s bin a -> (BinValue bin, a) -> ST s ()-fillOneW (HistogramST bin uo arr) (x,w)-    | i < 0             = writeMU uo  0 . (+w)  =<< readMU uo 0-    | i >= lengthMU arr = writeMU uo  1 . (+w)  =<< readMU uo 1-    | otherwise         = writeMU arr i . (+w)  =<< readMU arr i+fillOneW :: Num a => MHistogram s bin a -> (BinValue bin, a) -> ST s ()+fillOneW (MHistogram bin uo arr) (x,w)+    | i < 0              = MU.unsafeWrite uo  0 . (+w)  =<< MU.unsafeRead uo 0+    | i >= MU.length arr = MU.unsafeWrite uo  1 . (+w)  =<< MU.unsafeRead uo 1+    | otherwise          = MU.unsafeWrite arr i . (+w)  =<< MU.unsafeRead arr i     where       i = toIndex bin x+{-# INLINE fillOneW #-}   -- | Put one monoidal element-fillMonoid :: Monoid a => HistogramST s bin a -> (BinValue bin, a) -> ST s ()-fillMonoid (HistogramST bin uo arr) (x,m)-    | i < 0             = writeMU uo  1 . (flip mappend m)  =<< readMU uo  0-    | i >= lengthMU arr = writeMU uo  1 . (flip mappend m)  =<< readMU uo  1-    | otherwise         = writeMU arr i . (flip mappend m)  =<< readMU arr i-    where+fillMonoid :: Monoid a => MHistogram s bin a -> (BinValue bin, a) -> ST s ()+fillMonoid (MHistogram bin uo arr) (x,m)+    | i < 0              = MU.unsafeWrite uo  1 . (flip mappend m)  =<< MU.unsafeRead uo  0+    | i >= MU.length arr = MU.unsafeWrite uo  1 . (flip mappend m)  =<< MU.unsafeRead uo  1+    | otherwise          = MU.unsafeWrite arr i . (flip mappend m)  =<< MU.unsafeRead arr i+    where        i = toIndex bin x+{-# fillMonoid #-}  -- | Create immutable histogram from mutable one. This operation involve copying.-freezeHist :: HistogramST s bin a -> ST s (Histogram bin a)-freezeHist (HistogramST bin uo arr) = do-  [u,o] <- fromU `fmap` unsafeFreezeAllMU uo -- Is it safe???+freezeHist :: MHistogram s bin a -> ST s (Histogram bin a)+freezeHist (MHistogram bin uo arr) = do+  u <- MU.unsafeRead uo 0+  o <- MU.unsafeRead uo 1   -- Copy array-  let len = lengthMU arr-  tmp  <- newMU len-  memcpyOffMU arr tmp 0 0 len-  a    <- unsafeFreezeAllMU tmp-  return $ Histogram bin (Just (u,o)) a----------------------------------------------------------------------- Accumulator typeclass-------------------------------------------------------------------- | This is class with accumulation semantics. It's used to fill many---   histogram at once. It accept values of type a and return data of type b.-class Accumulator h where-    -- | Put one element into accumulator-    putOne  :: h s a b -> a   -> ST s () -    -- | Extract data from historam-    extract :: Monoid b => (h s a b) -> ST s b---- | Put many elements in histogram(s) at once -putMany :: Accumulator h => h s a b -> [a] -> ST s () -putMany !h = mapM_ (putOne h) ---- | Put all values into histogram and return result-fillHistograms :: Monoid b => (forall s . ST s (Accum s a b)) -> [a] -> b-fillHistograms h xs = runST $ do h' <- h-                                 putMany h' xs-                                 extract h'--------------------------------------------------------------------- GADT wrapper -------------------------------------------------------------------- | Abstract wrapper for histograms. -data Accum s a b where-    Accum :: Accumulator h => h s a b -> Accum s a b--instance Accumulator Accum where-    putOne  !(Accum h) !x = putOne h x -    extract !(Accum h)    = extract h---------------------------------------------------------------------- List of histograms------------------------------------------------------------------newtype AccumList s a b = AccumList [Accum s a b]- --- | Wrap list of histograms into one 'Accum'-accumList :: [ST s (Accum s a b)] -> ST s (Accum s a b)-accumList l = (Accum . AccumList) `fmap` sequence l--instance Accumulator AccumList where-    putOne  !(AccumList l) !x = mapM_ (flip putOne $ x) l -    extract !(AccumList l)    = mconcat `fmap` mapM extract l ---------------------------------------------------------------------- Generic histogram ------------------------------------------------------------------data AccumHist s a b where-    AccumHist :: (Bin bin) =>-                 (a -> HistogramST s bin val -> ST s ())-              -> (Histogram bin val -> b)-              -> HistogramST s bin val-              -> AccumHist s a b---- | Accumulator for arbitrary 'HistogramST' based histogram-accumHist :: (Bin bin) =>-             (a -> HistogramST s bin val -> ST s ())-          -> (Histogram bin val -> b)-          -> HistogramST s bin val-          -> ST s (Accum s a b)-accumHist inp out h = return . Accum $ AccumHist inp out h+  tmp  <- MU.new (MU.length arr)+  MU.copy tmp arr+  a    <- G.unsafeFreeze tmp+  return $ histogramUO bin (Just (u,o)) a+{-# INLINE freezeHist #-} -instance Accumulator AccumHist where-    putOne  !(AccumHist inp _ st) !x = inp x st-    extract !(AccumHist _ out st)    = out `fmap` freezeHist st

@@ -1,6 +1,6 @@ Name:           histogram-fill-Version:        0.1.0-Cabal-Version:  >= 1.2+Version:        0.2.0+Cabal-Version:  >= 1.6 License:        BSD3 License-File:   LICENSE Author:         Alexey Khudyakov@@ -15,12 +15,17 @@   .   This is very much work in progress so expect API breakage in future relesases. +source-repository head+  type:     hg+  location: http://bitbucket.org/Shimuuar/histogram-fill  Library-  Build-Depends:        base >=3 && <5, uvector >=0.1 && <0.2+  Build-Depends:        base >=3 && <5, vector   Exposed-modules:      Data.Histogram-                        Data.Histogram.Fill +                        Data.Histogram.Generic+                        Data.Histogram.Fill                         Data.Histogram.Bin+                        Data.Histogram.Bin.Extra                         Data.Histogram.ST   Other-modules:        Data.Histogram.Parse-  Ghc-options:          -O2 -Wall -auto-all+  Ghc-options:          -O2 -Wall