diff --git a/Data/Histogram.hs b/Data/Histogram.hs --- a/Data/Histogram.hs +++ b/Data/Histogram.hs @@ -1,6 +1,8 @@ + {-# LANGUAGE GADTs #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Data.Histogram -- Copyright : Copyright (c) 2009, Alexey Khudyakov @@ -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 diff --git a/Data/Histogram/Bin.hs b/Data/Histogram/Bin.hs --- a/Data/Histogram/Bin.hs +++ b/Data/Histogram/Bin.hs @@ -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 diff --git a/Data/Histogram/Bin/Extra.hs b/Data/Histogram/Bin/Extra.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/Extra.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE FlexibleContexts #-} +-- | +-- Module : Data.Histogram.Bin +-- Copyright : Copyright (c) 2010, Alexey Khudyakov +-- License : BSD3 +-- Maintainer : Alexey Khudyakov +-- 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) diff --git a/Data/Histogram/Fill.hs b/Data/Histogram/Fill.hs --- a/Data/Histogram/Fill.hs +++ b/Data/Histogram/Fill.hs @@ -1,188 +1,350 @@ {-# LANGUAGE GADTs #-} {-# LANGUAGE Rank2Types #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} -- | -- Module : Data.Histogram.Fill -- Copyright : Copyright (c) 2009, Alexey Khudyakov -- License : BSD3 -- Maintainer : Alexey Khudyakov -- 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 diff --git a/Data/Histogram/Generic.hs b/Data/Histogram/Generic.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Generic.hs @@ -0,0 +1,191 @@ +{-# LANGUAGE FlexibleContexts #-} +-- | +-- Module : Data.Histogram +-- Copyright : Copyright (c) 2009, Alexey Khudyakov +-- License : BSD3 +-- Maintainer : Alexey Khudyakov +-- 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 diff --git a/Data/Histogram/Parse.hs b/Data/Histogram/Parse.hs --- a/Data/Histogram/Parse.hs +++ b/Data/Histogram/Parse.hs @@ -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 diff --git a/Data/Histogram/ST.hs b/Data/Histogram/ST.hs --- a/Data/Histogram/ST.hs +++ b/Data/Histogram/ST.hs @@ -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 diff --git a/histogram-fill.cabal b/histogram-fill.cabal --- a/histogram-fill.cabal +++ b/histogram-fill.cabal @@ -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