diff --git a/Data/Histogram.hs b/Data/Histogram.hs --- a/Data/Histogram.hs +++ b/Data/Histogram.hs @@ -1,8 +1,6 @@ - {-# LANGUAGE GADTs #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE TypeOperators #-} -{-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Data.Histogram -- Copyright : Copyright (c) 2009, Alexey Khudyakov @@ -30,7 +28,10 @@ -- ** Convert to other data types , asList , asVector - -- * Slicing histograms + -- * Slicing histogram + , sliceByIx + , sliceByVal + -- * Splitting 2D histograms , sliceX , sliceY -- * Modify histogram @@ -136,6 +137,12 @@ histZipSafe :: (Bin bin, Eq bin, Unbox a, Unbox b, Unbox c) => (a -> b -> c) -> Histogram bin a -> Histogram bin b -> Maybe (Histogram bin c) histZipSafe = H.histZipSafe + +sliceByIx :: (Bin1D bin, Unbox a) => Int -> Int -> Histogram bin a -> Histogram bin a +sliceByIx = H.sliceByIx + +sliceByVal :: (Bin1D bin, Unbox a) => BinValue bin -> BinValue bin -> Histogram bin a -> Histogram bin a +sliceByVal = H.sliceByVal -- | Slice 2D histogram along Y axis. This function is fast because it does not require reallocations. sliceY :: (Unbox a, Bin bX, Bin bY) => Histogram (Bin2D bX bY) a -> [(BinValue bY, Histogram bX a)] diff --git a/Data/Histogram/Bin.hs b/Data/Histogram/Bin.hs --- a/Data/Histogram/Bin.hs +++ b/Data/Histogram/Bin.hs @@ -1,10 +1,8 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE DeriveDataTypeable #-} -- Requred for Bin2D conversions +{-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE OverlappingInstances #-} +-- Yes I DO want orphans here +{-# OPTIONS_GHC -fno-warn-orphans #-} -- | -- Module : Data.Histogram.Bin -- Copyright : Copyright (c) 2009, Alexey Khudyakov @@ -16,565 +14,22 @@ -- indices and approximate reverse. module Data.Histogram.Bin ( -- * Type classes - Bin(..) - , Bin1D(..) - , UniformBin1D(..) - , VariableBin1D(..) - , ConvertBin(..) - -- * Bin types - -- ** Integer bins - , BinI(..) - , binI0 - -- ** Integer bins with non-1 size - , BinInt(..) - , binInt - -- ** Enum based bin - , BinEnum(..) - , binEnum - , binEnumFull - -- ** Floating point bins - , BinF - , binF - , binFn - , binFstep - , scaleBinF - -- *** Specialized for Double - , BinD - , binD - , binDn - , binDstep - , scaleBinD - -- ** Log scale point - , LogBinD - , logBinD - -- ** 2D bins - , Bin2D(..) - , (><) - , nBins2D - , toIndex2D - , fmapBinX - , fmapBinY - ) where - -import Control.Monad (liftM, liftM2, liftM3) -import GHC.Float (double2Int) - -import qualified Data.Vector.Generic as G -import Data.Vector.Generic (Vector) -import Data.Typeable (Typeable) -import Text.Read (Read(..)) - -import Data.Histogram.Parse - - - ----------------------------------------------------------------- --- Type classes ----------------------------------------------------------------- - --- | This type represent some abstract data binning algorithms. Such --- algorithm maps sets of values to integer indices. --- --- Following invariant is expected to hold: --- --- > toIndex . fromIndex == id -class Bin b where - -- | Type of value to bin - type BinValue b - -- | Convert from value to index. Function must not fail for any - -- input and should produce out of range indices for invalid input. - toIndex :: b -> BinValue b -> Int - -- | Convert from index to value. Returned value should correspond - -- to center of bin. Definition of center is left for definition - -- of instance. Funtion may fail for invalid indices but - -- encouraged not to do so. - fromIndex :: b -> Int -> BinValue b - -- | Check whether value in range. Values which lay in range must - -- produce valid indices and conversely value which produce - -- valid index must be 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 - -- | Minimal accepted value of histogram - lowerLimit :: b -> BinValue b - -- | Maximal accepted value of histogram - upperLimit :: b -> BinValue b - -- | List of center of bins in ascending order. Default - -- implementation is: - -- - -- > binsList b = G.generate (nBins b) (fromIndex b) - binsList :: Vector v (BinValue b) => b -> v (BinValue b) - binsList b = G.generate (nBins b) (fromIndex b) - -- | List of bins in ascending order. First element of tuple is - -- lower bound second is upper bound of bin - binsListRange :: Vector v (BinValue b, BinValue b) => b -> v (BinValue b, BinValue b) - {-# INLINE binsList #-} - - --- | 1D binning algorithms with variable bin size -class Bin1D b => VariableBin1D b where - -- | Size of n'th bin. - binSizeN :: b -> Int -> BinValue b - - --- | 1D binning algorithms with constant size bins. Constant sized --- bins could be thought as specialization of variable-sized bins --- therefore a superclass constraint. -class VariableBin1D b => UniformBin1D b where - -- | Size of bin. Default implementation just uses 0 bin. - binSize :: b -> BinValue b - binSize b = binSizeN b 0 - - --- | Class for conversion between binning algorithms -class (Bin b, Bin b') => ConvertBin b b' where - -- | Convert bins - convertBin :: b -> b' - ----------------------------------------------------------------- --- Integer bin ----------------------------------------------------------------- --- | Simple binning algorithm which map continous range of bins onto --- indices. Each number correcsponds to different bin --- --- 1. Lower bound (inclusive) --- --- 2. Upper bound (inclusive) -data BinI = BinI - {-# UNPACK #-} !Int -- Lower bound (inclusive) - {-# UNPACK #-} !Int -- Upper bound (inclusive) - deriving (Eq,Typeable) - --- | 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 - fromIndex !(BinI base _) !x = x + base - inRange !(BinI x y) i = i>=x && i<=y - nBins !(BinI x y) = y - x + 1 - {-# INLINE toIndex #-} - {-# INLINE inRange #-} - -instance Bin1D BinI where - lowerLimit (BinI i _) = i - upperLimit (BinI _ i) = i - binsList b@(BinI lo _) = G.enumFromN lo (nBins b) - binsListRange b@(BinI lo _) = G.generate (nBins b) (\i -> let n = lo+i in (n,n)) - {-# INLINE binsList #-} - {-# INLINE binsListRange #-} - -instance VariableBin1D BinI where - binSizeN _ _ = 1 - -instance UniformBin1D BinI where - binSize _ = 1 - -instance Show BinI where - show (BinI lo hi) = unlines [ "# BinI" - , "# Low = " ++ show lo - , "# High = " ++ show hi - ] -instance Read BinI where - readPrec = keyword "BinI" >> liftM2 BinI (value "Low") (value "High") - - - ----------------------------------------------------------------- --- Another form of Integer bin ----------------------------------------------------------------- - --- | Integer bins with size which differ from 1. --- --- 1. Low bound --- --- 2. Bin size --- --- 3. Number of bins -data BinInt = BinInt - {-# UNPACK #-} !Int -- Low bound - {-# UNPACK #-} !Int -- Bin size - {-# UNPACK #-} !Int -- Number of bins - deriving (Eq,Typeable) - --- | 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 - fromIndex !(BinInt base sz _) !x = x * sz + base - inRange !(BinInt base sz n) i = i>=base && i<(base+n*sz) - nBins !(BinInt _ _ n) = n - {-# INLINE toIndex #-} - {-# INLINE inRange #-} - -instance Bin1D BinInt where - lowerLimit (BinInt base _ _) = base - upperLimit (BinInt base sz n) = base + sz * n - 1 - binsListRange b@(BinInt _ sz n) = G.generate n (\i -> let x = fromIndex b i in (x,x + sz - 1)) - -instance VariableBin1D BinInt where - binSizeN (BinInt _ sz _) _ = sz - -instance UniformBin1D BinInt where - binSize (BinInt _ sz _) = sz - -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") - - ----------------------------------------------------------------- --- Enumeration bin ----------------------------------------------------------------- - --- | Bin for types which are instnaces of Enum type class -newtype BinEnum a = BinEnum BinI - deriving (Eq,Typeable) - --- | Create enum based bin -binEnum :: Enum a => a -> a -> BinEnum a -binEnum a b = BinEnum $ BinI (fromEnum a) (fromEnum b) - --- | Use full range of data -binEnumFull :: (Enum a, Bounded a) => BinEnum a -binEnumFull = binEnum minBound maxBound - -instance Enum a => Bin (BinEnum a) where - type BinValue (BinEnum a) = a - toIndex (BinEnum b) = toIndex b . fromEnum - fromIndex (BinEnum b) = toEnum . fromIndex b - inRange (BinEnum b) = inRange b . fromEnum - nBins (BinEnum b) = nBins b - -instance Enum a => Bin1D (BinEnum a) where - lowerLimit (BinEnum b) = toEnum $ lowerLimit b - upperLimit (BinEnum b) = toEnum $ upperLimit b - binsListRange b = G.generate (nBins b) (\n -> let x = fromIndex b n in (x,x)) - {-# INLINE binsListRange #-} - -instance Show (BinEnum a) where - show (BinEnum b) = "# BinEnum\n" ++ show b -instance Read (BinEnum a) where - readPrec = keyword "BinEnum" >> liftM BinEnum readPrec - - - ----------------------------------------------------------------- --- Floating point bin ----------------------------------------------------------------- - --- | Floaintg point bins with equal sizes. --- --- Note that due to GHC bug #2271 this toIndex is really slow (20x --- slowdown with respect to BinD) and use of BinD is recommended --- --- 1. Lower bound --- --- 2. Size of bin --- --- 3. Number of bins -data BinF f = BinF {-# UNPACK #-} !f -- Lower bound - {-# UNPACK #-} !f -- Size of bin - {-# UNPACK #-} !Int -- Number of bins - deriving (Eq,Typeable) - --- | Create bins. -binF :: RealFrac f => - f -- ^ Lower bound of range - -> Int -- ^ Number of bins - -> f -- ^ Upper bound of range - -> BinF f -binF from n to = BinF from ((to - from) / fromIntegral n) n - --- | Create bins. Note that actual upper bound can differ from specified. -binFn :: RealFrac f => - f -- ^ Begin of range - -> f -- ^ Size of step - -> f -- ^ Approximation of end of range - -> BinF f -binFn from step to = BinF from step (round $ (to - from) / step) - --- | Create bins -binFstep :: RealFrac f => - f -- ^ Begin of range - -> f -- ^ Size of step - -> Int -- ^ Number of bins - -> BinF f -binFstep = BinF - --- | '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 RealFrac f => Bin (BinF f) where - type BinValue (BinF f) = f - toIndex !(BinF from step _) !x = floor $ (x-from) / step - fromIndex !(BinF from step _) !i = (step/2) + (fromIntegral i * step) + from - inRange !(BinF from step n) x = x > from && x < from + step*fromIntegral n - nBins !(BinF _ _ n) = n - {-# INLINE toIndex #-} - {-# INLINE inRange #-} - -instance RealFrac f => Bin1D (BinF f) where - lowerLimit (BinF from _ _) = from - upperLimit (BinF from step n) = from + step * fromIntegral n - binsListRange !b@(BinF _ step n) = G.generate n toPair - where - toPair k = (x - step/2, x + step/2) where x = fromIndex b k - {-# INLINE binsListRange #-} - -instance RealFrac f => VariableBin1D (BinF f) where - binSizeN (BinF _ step _) _ = step - -instance RealFrac f => UniformBin1D (BinF f) where - binSize (BinF _ step _) = step - -instance Show f => Show (BinF f) where - show (BinF base step n) = unlines [ "# BinF" - , "# Base = " ++ show base - , "# Step = " ++ show step - , "# N = " ++ show n - ] -instance (Read f, RealFrac f) => Read (BinF f) where - readPrec = keyword "BinF" >> liftM3 BinF (value "Base") (value "Step") (value "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. --- --- 1. Lower bound --- --- 2. Size of bin --- --- 3. Number of bins -data BinD = BinD {-# UNPACK #-} !Double -- Lower bound - {-# UNPACK #-} !Double -- Size of bin - {-# UNPACK #-} !Int -- Number of bins - deriving (Eq,Typeable) - --- | 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) - --- | Create bins -binDstep :: Double -- ^ Begin of range - -> Double -- ^ Size of step - -> Int -- ^ Number of bins - -> BinD -binDstep = BinD - --- | '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 - fromIndex !(BinD from step _) !i = (step/2) + (fromIntegral i * step) + from - inRange !(BinD from step n) x = x > from && x < from + step*fromIntegral n - nBins !(BinD _ _ n) = n - {-# INLINE toIndex #-} - {-# INLINE inRange #-} - -instance Bin1D BinD where - lowerLimit (BinD from _ _) = from - upperLimit (BinD from step n) = from + step * fromIntegral n - binsListRange b@(BinD _ step n) = G.generate n toPair - where - toPair k = (x - step/2, x + step/2) where x = fromIndex b k - {-# INLINE binsListRange #-} - - -instance VariableBin1D BinD where - binSizeN (BinD _ step _) _ = step - -instance UniformBin1D BinD where - binSize (BinD _ step _) = step - -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 = keyword "BinD" >> liftM3 BinD (value "Base") (value "Step") (value "N") - - - ----------------------------------------------------------------- --- Log-scale bin ----------------------------------------------------------------- --- | Logarithmic scale bins. --- --- 1. Lower bound --- --- 2. Upper bound --- --- 2. Increment ratio --- --- 3. Number of bins -data LogBinD = LogBinD - Double -- Low border - Double -- Hi border - Double -- Increment ratio - Int -- Number of bins - deriving (Eq,Typeable) - --- | 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) - fromIndex !(LogBinD base _ step _) !i | i >= 0 = base * step ** (fromIntegral i + 0.5) - | otherwise = -1 / 0 - inRange !(LogBinD lo hi _ _) x = x >= lo && x < hi - nBins !(LogBinD _ _ _ n) = n - {-# INLINE toIndex #-} - {-# INLINE inRange #-} - -instance Bin1D LogBinD where - lowerLimit (LogBinD lo _ _ _) = lo - upperLimit (LogBinD _ hi _ _) = hi - binsListRange (LogBinD base _ step n) = G.unfoldrN n next base - where - next x = let x' = x * step in Just ((x,x'), x') - {-# INLINE binsListRange #-} - -instance VariableBin1D LogBinD where - binSizeN (LogBinD base _ step _) n = let x = base * step ^ n in x*step - x - -instance Show LogBinD where - show (LogBinD lo hi _ n) = - unlines [ "# LogBinD" - , "# Lo = " ++ show lo - , "# N = " ++ show n - , "# Hi = " ++ show hi - ] -instance Read LogBinD where - readPrec = do - keyword "LogBinD" - liftM3 logBinD (value "Lo") (value "N") (value "Hi") - - ----------------------------------------------------------------- --- 2D bin ----------------------------------------------------------------- - --- | 2D bins. binX is binning along X axis and binY is one along Y axis. -data Bin2D binX binY = Bin2D { binX :: !binX -- ^ Binning algorithm for X axis - , binY :: !binY -- ^ Binning algorithm for Y axis - } - deriving (Eq,Typeable) - --- | Alias for 'Bin2D'. -(><) :: binX -> binY -> Bin2D binX binY -(><) = Bin2D - -instance (Bin binX, Bin binY) => Bin (Bin2D binX binY) where - type BinValue (Bin2D binX binY) = (BinValue binX, BinValue binY) - toIndex !(Bin2D bx by) !(x,y) - | inRange bx x = toIndex bx x + toIndex by y * nBins bx - | otherwise = maxBound - 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 - nBins (Bin2D bx by) = nBins bx * nBins by - {-# INLINE toIndex #-} - {-# INLINE inRange #-} - --- | Convert index into pair of indices for X and Y axes -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) -{-# INLINE toIndex2D #-} - --- | 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 + module Data.Histogram.Bin.Classes + , module Data.Histogram.Bin.BinI + , module Data.Histogram.Bin.BinInt + , module Data.Histogram.Bin.BinEnum + , module Data.Histogram.Bin.BinF + , module Data.Histogram.Bin.LogBinD + , module Data.Histogram.Bin.Bin2D + ) where -instance (Show b1, Show b2) => Show (Bin2D b1 b2) where - 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" - keyword "X" - b1 <- readPrec - keyword "Y" - b2 <- readPrec - return $ Bin2D b1 b2 +import Data.Histogram.Bin.Classes +import Data.Histogram.Bin.BinI +import Data.Histogram.Bin.BinInt +import Data.Histogram.Bin.BinEnum +import Data.Histogram.Bin.BinF +import Data.Histogram.Bin.LogBinD +import Data.Histogram.Bin.Bin2D ---------------------------------------------------------------- -- Bin conversion diff --git a/Data/Histogram/Bin/Bin2D.hs b/Data/Histogram/Bin/Bin2D.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/Bin2D.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Data.Histogram.Bin.Bin2D ( + Bin2D(..) + , (><) + , nBins2D + , fmapBinX + , fmapBinY + ) where + +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) + +import Data.Histogram.Bin.Classes +import Data.Histogram.Parse + +-- | 2D bins. binX is binning along X axis and binY is one along Y +-- axis. Data is stored in row-major order +data Bin2D binX binY = Bin2D { binX :: !binX -- ^ Binning algorithm for X axis + , binY :: !binY -- ^ Binning algorithm for Y axis + } + deriving (Eq,Data,Typeable) + +-- | Alias for 'Bin2D'. +(><) :: binX -> binY -> Bin2D binX binY +(><) = Bin2D + +instance (Bin binX, Bin binY) => Bin (Bin2D binX binY) where + type BinValue (Bin2D binX binY) = (BinValue binX, BinValue binY) + toIndex !(Bin2D bx by) !(x,y) + | inRange bx x = toIndex bx x + toIndex by y * nBins bx + | otherwise = maxBound + 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 + nBins (Bin2D bx by) = nBins bx * nBins by + {-# INLINE toIndex #-} + +-- | Convert index into pair of indices for X and Y axes +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) +{-# INLINE toIndex2D #-} + +-- | 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) = 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" + keyword "X" + b1 <- readPrec + keyword "Y" + b2 <- readPrec + return $ Bin2D b1 b2 diff --git a/Data/Histogram/Bin/BinEnum.hs b/Data/Histogram/Bin/BinEnum.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/BinEnum.hs @@ -0,0 +1,50 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +module Data.Histogram.Bin.BinEnum ( + BinEnum(..) + , binEnum + , binEnumFull + ) where + +import Control.Monad (liftM) +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) + +import Data.Histogram.Bin.Classes +import Data.Histogram.Bin.BinI +import Data.Histogram.Parse + +-- | Bin for types which are instnaces of Enum type class +newtype BinEnum a = BinEnum BinI + deriving (Eq,Data,Typeable,GrowBin) + +-- | Create enum based bin +binEnum :: Enum a => a -> a -> BinEnum a +binEnum a b = BinEnum $ binI (fromEnum a) (fromEnum b) + +-- | Use full range of data +binEnumFull :: (Enum a, Bounded a) => BinEnum a +binEnumFull = binEnum minBound maxBound + +instance Enum a => Bin (BinEnum a) where + type BinValue (BinEnum a) = a + toIndex (BinEnum b) = toIndex b . fromEnum + fromIndex (BinEnum b) = toEnum . fromIndex b + inRange (BinEnum b) = inRange b . fromEnum + nBins (BinEnum b) = nBins b + +instance Enum a => IntervalBin (BinEnum a) where + binInterval b x = (n,n) where n = fromIndex b x + +instance Enum a => Bin1D (BinEnum a) where + lowerLimit (BinEnum b) = toEnum $ lowerLimit b + upperLimit (BinEnum b) = toEnum $ upperLimit b + unsafeSliceBin i j (BinEnum b) = BinEnum $ unsafeSliceBin i j b + +instance Show (BinEnum a) where + show (BinEnum b) = "# BinEnum\n" ++ show b +instance Read (BinEnum a) where + readPrec = keyword "BinEnum" >> liftM BinEnum readPrec diff --git a/Data/Histogram/Bin/BinF.hs b/Data/Histogram/Bin/BinF.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/BinF.hs @@ -0,0 +1,193 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Data.Histogram.Bin.BinF ( + -- * Generic and slow + BinF(..) + , binF + , binFn + , binFstep + , scaleBinF + -- * Specialized for Double and fast + , BinD(..) + , binD + , binDn + , binDstep + , scaleBinD + ) where + +import Control.Monad (liftM3) +import GHC.Float (double2Int) +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) + +import Data.Histogram.Bin.Classes +import Data.Histogram.Parse + + +-- | Floaintg point bins with equal sizes. +-- +-- Note that due to GHC bug #2271 this toIndex is really slow (20x +-- slowdown with respect to BinD) and use of BinD is recommended +-- +-- 1. Lower bound +-- +-- 2. Size of bin +-- +-- 3. Number of bins +data BinF f = BinF !f -- Lower bound + !f -- Size of bin + {-# UNPACK #-} !Int -- Number of bins + deriving (Eq,Data,Typeable) + +-- | Create bins. +binF :: RealFrac f => + f -- ^ Lower bound of range + -> Int -- ^ Number of bins + -> f -- ^ Upper bound of range + -> BinF f +binF from n to = BinF from ((to - from) / fromIntegral n) n + +-- | Create bins. Note that actual upper bound can differ from specified. +binFn :: RealFrac f => + f -- ^ Begin of range + -> f -- ^ Size of step + -> f -- ^ Approximation of end of range + -> BinF f +binFn from step to = BinF from step (round $ (to - from) / step) + +-- | Create bins +binFstep :: RealFrac f => + f -- ^ Begin of range + -> f -- ^ Size of step + -> Int -- ^ Number of bins + -> BinF f +binFstep = BinF + +-- | '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 RealFrac f => Bin (BinF f) where + type BinValue (BinF f) = f + toIndex !(BinF from step _) !x = floor $ (x-from) / step + fromIndex !(BinF from step _) !i = (step/2) + (fromIntegral i * step) + from + nBins !(BinF _ _ n) = n + {-# INLINE toIndex #-} + +instance RealFrac f => IntervalBin (BinF f) where + binInterval (BinF from step _) i = (x, x + step) where x = from + step * fromIntegral i + +instance RealFrac f => Bin1D (BinF f) where + lowerLimit (BinF from _ _) = from + upperLimit (BinF from step n) = from + step * fromIntegral n + unsafeSliceBin i j (BinF from step _) = BinF (from + step * fromIntegral i) step (j-i+1) + +instance RealFrac f => GrowBin (BinF f) where + zeroBin (BinF from step _) = BinF from step 0 + appendBin (BinF from step n) = BinF from step (n+1) + prependBin (BinF from step n) = BinF (from-step) step (n+1) + +instance RealFrac f => VariableBin (BinF f) where + binSizeN (BinF _ step _) _ = step + +instance RealFrac f => UniformBin (BinF f) where + binSize (BinF _ step _) = step + +instance Show f => Show (BinF f) where + show (BinF base step n) = unlines [ "# BinF" + , "# Base = " ++ show base + , "# Step = " ++ show step + , "# N = " ++ show n + ] +instance (Read f, RealFrac f) => Read (BinF f) where + readPrec = keyword "BinF" >> liftM3 BinF (value "Base") (value "Step") (value "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. +-- +-- 1. Lower bound +-- +-- 2. Size of bin +-- +-- 3. Number of bins +data BinD = BinD {-# UNPACK #-} !Double -- Lower bound + {-# UNPACK #-} !Double -- Size of bin + {-# UNPACK #-} !Int -- Number of bins + deriving (Eq,Data,Typeable) + +-- | 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) + +-- | Create bins +binDstep :: Double -- ^ Begin of range + -> Double -- ^ Size of step + -> Int -- ^ Number of bins + -> BinD +binDstep = BinD + +-- | '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 + fromIndex !(BinD from step _) !i = (step/2) + (fromIntegral i * step) + from + nBins !(BinD _ _ n) = n + {-# INLINE toIndex #-} + +instance IntervalBin BinD where + binInterval (BinD from step _) i = (x, x + step) where x = from + step * fromIntegral i + +instance Bin1D BinD where + lowerLimit (BinD from _ _) = from + upperLimit (BinD from step n) = from + step * fromIntegral n + unsafeSliceBin i j (BinD from step _) = BinD (from + step * fromIntegral i) step (j-i+1) + +instance GrowBin BinD where + zeroBin (BinD from step _) = BinD from step 0 + appendBin (BinD from step n) = BinD from step (n+1) + prependBin (BinD from step n) = BinD (from-step) step (n+1) + +instance VariableBin BinD where + binSizeN (BinD _ step _) _ = step + +instance UniformBin BinD where + binSize (BinD _ step _) = step + +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 = keyword "BinD" >> liftM3 BinD (value "Base") (value "Step") (value "N") diff --git a/Data/Histogram/Bin/BinI.hs b/Data/Histogram/Bin/BinI.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/BinI.hs @@ -0,0 +1,74 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Data.Histogram.Bin.BinI ( + BinI(..) + , binI + , binI0 + ) where + +import Control.Monad (liftM2) +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) + +import Data.Histogram.Bin.Classes +import Data.Histogram.Parse + + + +-- | Simple binning algorithm which map continous range of bins onto +-- indices. Each number correcsponds to different bin +-- +-- 1. Lower bound (inclusive) +-- +-- 2. Upper bound (inclusive) +data BinI = BinI + {-# UNPACK #-} !Int -- Lower bound (inclusive) + {-# UNPACK #-} !Int -- Upper bound (inclusive) + deriving (Eq,Data,Typeable) + +-- | Safe constructor for BinI. It does checks that upper bound is +-- greater or equal than lower bound +binI :: Int -> Int -> BinI +binI lo hi | lo <= hi = BinI lo hi + | otherwise = error "Data.Histogram.Bin.BinI.binI: invalid paramters" + +-- | Construct BinI with n bins. Indexing starts from 0. n must be positive +binI0 :: Int -> BinI +binI0 n = binI 0 (n - 1) + +instance Bin BinI where + type BinValue BinI = Int + toIndex !(BinI base _) !x = x - base + fromIndex !(BinI base _) !x = x + base + inRange !(BinI x y) i = i>=x && i<=y + nBins !(BinI x y) = y - x + 1 + {-# INLINE toIndex #-} + +instance IntervalBin BinI where + binInterval b i = (n,n) where n = fromIndex b i + +instance Bin1D BinI where + lowerLimit (BinI i _) = i + upperLimit (BinI _ i) = i + unsafeSliceBin i j (BinI l _) = BinI (l+i) (l+j) + +instance VariableBin BinI where + binSizeN _ _ = 1 + +instance UniformBin BinI where + binSize _ = 1 + +instance GrowBin BinI where + zeroBin (BinI l _) = BinI l l + appendBin (BinI l u) = BinI l (u+1) + prependBin (BinI l u) = BinI (l-1) u + +instance Show BinI where + show (BinI lo hi) = unlines [ "# BinI" + , "# Low = " ++ show lo + , "# High = " ++ show hi + ] +instance Read BinI where + readPrec = keyword "BinI" >> liftM2 BinI (value "Low") (value "High") diff --git a/Data/Histogram/Bin/BinInt.hs b/Data/Histogram/Bin/BinInt.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/BinInt.hs @@ -0,0 +1,89 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Data.Histogram.Bin.BinInt ( + BinInt(..) + , binInt + , binIntN + ) where + +import Control.Monad (liftM3) +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) + +import Data.Histogram.Bin.Classes +import Data.Histogram.Parse + + + +-- | Integer bins with size which differ from 1. +-- +-- 1. Low bound +-- +-- 2. Bin size +-- +-- 3. Number of bins +data BinInt = BinInt + {-# UNPACK #-} !Int -- Low bound + {-# UNPACK #-} !Int -- Bin size + {-# UNPACK #-} !Int -- Number of bins + deriving (Eq,Data,Typeable) + +-- FIXME: no sanity checks +-- | 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 + +binIntN :: Int -- ^ Lower bound + -> Int -- ^ Bin size + -> Int -- ^ Upper bound + -> BinInt +binIntN lo n hi + | n > rng = BinInt lo 1 rng + | otherwise = BinInt lo undefined n + where + rng = hi - lo + 1 + + +instance Bin BinInt where + type BinValue BinInt = Int + toIndex !(BinInt base sz _) !x = (x - base) `div` sz + fromIndex !(BinInt base sz _) !x = x * sz + base + nBins !(BinInt _ _ n) = n + {-# INLINE toIndex #-} + +instance IntervalBin BinInt where + binInterval b i = (n, n + binSize b - 1) where n = fromIndex b i + +instance Bin1D BinInt where + lowerLimit (BinInt base _ _) = base + upperLimit (BinInt base sz n) = base + sz * n - 1 + unsafeSliceBin i j (BinInt base sz _) = BinInt (base + i*sz) sz (j-i+1) + +instance GrowBin BinInt where + zeroBin (BinInt l sz _) = BinInt l sz 0 + appendBin (BinInt l sz n) = BinInt l sz (n+1) + prependBin (BinInt l sz n) = BinInt (l-sz) sz (n+1) + +instance VariableBin BinInt where + binSizeN (BinInt _ sz _) _ = sz + +instance UniformBin BinInt where + binSize (BinInt _ sz _) = sz + +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") diff --git a/Data/Histogram/Bin/Classes.hs b/Data/Histogram/Bin/Classes.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/Classes.hs @@ -0,0 +1,127 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +-- | +-- Module : Data.Histogram.Bin +-- Copyright : Copyright (c) 2011, Alexey Khudyakov +-- License : BSD3 +-- Maintainer : Alexey Khudyakov +-- Stability : experimental +-- +-- Type classes for binning algorithms. This is mapping from set of +-- interest to integer indices and approximate reverse. +module Data.Histogram.Bin.Classes ( + -- * Bin type class + Bin(..) + , binsCenters + -- * 1D bins + , IntervalBin(..) + , Bin1D(..) + , sliceBin + , VariableBin(..) + , UniformBin(..) + , GrowBin(..) + -- * Conversion + , ConvertBin(..) + ) where + +import qualified Data.Vector.Generic as G +import Data.Vector.Generic (Vector) + + +-- | This type represent some abstract data binning algorithms. It +-- maps sets/intervals of values of type 'BinValue b' to integer +-- indices. +-- +-- Following invariant is expected to hold: +-- +-- > toIndex . fromIndex == id +class Bin b where + -- | Type of value to bin + type BinValue b + -- | Convert from value to index. Function must not fail for any + -- input and should produce out of range indices for invalid input. + toIndex :: b -> BinValue b -> Int + -- | Convert from index to value. Returned value should correspond + -- to center of bin. Definition of center is left for definition + -- of instance. Funtion may fail for invalid indices but + -- encouraged not to do so. + fromIndex :: b -> Int -> BinValue b + -- | Total number of bins. + nBins :: b -> Int + -- | Check whether value in range. Have default + -- implementation. Should satisfy: + -- inRange b x ⇔ toIndex b x ∈ [0,nBins b) + inRange :: b -> BinValue b -> Bool + inRange b x = i >= 0 && i < nBins b where i = toIndex b x + +-- | Return vector of bin centers +binsCenters :: (Bin b, Vector v (BinValue b)) => b -> v (BinValue b) +binsCenters b = G.generate (nBins b) (fromIndex b) +{-# INLINE binsCenters #-} + +---------------------------------------------------------------- +-- 1D bins +---------------------------------------------------------------- + +-- | For binning algorithms which work with bin values which have some +-- natural ordering and every bin is continous interval. +class Bin b => IntervalBin b where + -- | Interval for n'th bin + binInterval :: b -> Int -> (BinValue b, BinValue b) + -- | List of all bins. Could be overridden for efficiency. + binsList :: Vector v (BinValue b, BinValue b) => b -> v (BinValue b, BinValue b) + binsList b = G.generate (nBins b) (binInterval b) + {-# INLINE binsList #-} + + +-- | IntervalBin for which domain is single finite interval +class IntervalBin b => Bin1D b where + -- | Minimal accepted value of histogram + lowerLimit :: b -> BinValue b + -- | Maximal accepted value of histogram + upperLimit :: b -> BinValue b + -- | Slice bin by indices. This function doesn't perform any checks + -- and may produce invalid bin + unsafeSliceBin :: Int -> Int -> b -> b + +-- | Slice bin using indices +sliceBin :: Bin1D b => Int -> Int -> b -> b +sliceBin i j b + | i < 0 || j < 0 || i > j || i >= n || j >= n = error "sliceBin: bad slice" + | otherwise = unsafeSliceBin i j b + where + n = nBins b + +-- | Binning algorithm which individual +class Bin1D b => GrowBin b where + -- | Set numbers to zero. By convention bins are shrinked to lower bound + zeroBin :: b -> b + -- | Append one bin at upper bound + appendBin :: b -> b + -- | Prepend one bin at lower bin + prependBin :: b -> b + +---- Bin sizes ------------------------------------------------ + +-- | 1D binning algorithms with variable bin size +class Bin b => VariableBin b where + -- | Size of n'th bin. + binSizeN :: b -> Int -> BinValue b + + +-- | 1D binning algorithms with constant size bins. Constant sized +-- bins could be thought as specialization of variable-sized bins +-- therefore a superclass constraint. +class VariableBin b => UniformBin b where + -- | Size of bin. Default implementation just uses 0th bin. + binSize :: b -> BinValue b + binSize b = binSizeN b 0 + + +---- Conversion ------------------------------------------------ + +-- | Class for conversion between binning algorithms. +class (Bin b, Bin b') => ConvertBin b b' where + -- | Convert bins + convertBin :: b -> b' diff --git a/Data/Histogram/Bin/Extra.hs b/Data/Histogram/Bin/Extra.hs --- a/Data/Histogram/Bin/Extra.hs +++ b/Data/Histogram/Bin/Extra.hs @@ -3,6 +3,7 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Data.Histogram.Bin -- Copyright : Copyright (c) 2010, Alexey Khudyakov @@ -20,14 +21,16 @@ , permuteBin ) where -import Control.Applicative -import Control.Monad -- (forM_,liftM2) +import Control.Applicative ((<$>), Applicative(..)) +import Control.Monad (forM_,liftM2, guard) +import Control.Monad.ST (ST) -import qualified Data.Vector.Generic as G import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as M import Data.Vector.Generic ((!)) -import Text.Read +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) import Data.Histogram.Bin import Data.Histogram.Parse @@ -54,6 +57,7 @@ -- | Binning for 2D enumerations newtype BinEnum2D i = BinEnum2D (Bin2D BinI BinI) + deriving (Eq,Data,Typeable) -- | Construct indexed bin binEnum2D :: Enum2D i => i -> i -> BinEnum2D i @@ -88,6 +92,7 @@ , permuteTo :: U.Vector Int -- ^ Maps original bin's indices to new indices , permuteFrom :: U.Vector Int -- ^ Inverse of pervious table } + deriving (Eq,Data,Typeable) instance Bin b => Bin (BinPermute b) where type BinValue (BinPermute b) = BinValue b @@ -96,26 +101,13 @@ inRange (BinPermute b _ _) x = inRange b x nBins = nBins . permutedBin -instance (Bin1D b) => Bin1D (BinPermute b) where - lowerLimit = lowerLimit . permutedBin - upperLimit = upperLimit . permutedBin - binsList (BinPermute b _ a) = res - where - res = G.generate (nBins b) fun - arr = binsList b `asTypeOf` res - fun i = arr ! (a ! i) - binsListRange (BinPermute b _ a) = res - where - res = G.generate (nBins b) fun - arr = binsListRange b `asTypeOf` res - fun i = arr ! (a ! i) - {-# INLINE binsList #-} - {-# INLINE binsListRange #-} +instance IntervalBin b => IntervalBin (BinPermute b) where + binInterval b i = binInterval (permutedBin b) (permuteFrom b ! i) -instance VariableBin1D b => VariableBin1D (BinPermute b) where +instance VariableBin b => VariableBin (BinPermute b) where binSizeN b i = binSizeN (permutedBin b) (permuteFrom b ! i) -instance UniformBin1D b => UniformBin1D (BinPermute b) where +instance UniformBin b => UniformBin (BinPermute b) where binSize = binSize . permutedBin @@ -151,6 +143,7 @@ return a where n = U.length v + writeInvert :: M.MVector s Int -> Int -> ST s () writeInvert a i | j >= 0 && j < n = M.write a j i | otherwise = return () where j = v ! i diff --git a/Data/Histogram/Bin/LogBinD.hs b/Data/Histogram/Bin/LogBinD.hs new file mode 100644 --- /dev/null +++ b/Data/Histogram/Bin/LogBinD.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Data.Histogram.Bin.LogBinD ( + -- * Generic and slow + LogBinD(..) + , logBinD + ) where + +import Control.Monad (liftM3) +import GHC.Float (double2Int) +import Data.Typeable (Typeable) +import Data.Data (Data) +import Text.Read (Read(..)) + +import Data.Histogram.Bin.Classes +import Data.Histogram.Parse +-- | Logarithmic scale bins. +-- +-- 1. Lower bound +-- +-- 2. Increment ratio +-- +-- 3. Number of bins +data LogBinD = LogBinD + Double -- Low border + Double -- Increment ratio + Int -- Number of bins + deriving (Eq,Data,Typeable) + +-- | Create log-scale bins. +logBinD :: Double -> Int -> Double -> LogBinD +logBinD lo n hi = LogBinD lo ((hi/lo) ** (1 / fromIntegral n)) n + +-- Fast variant of flooor +floorD :: Double -> Int +floorD x | x < 0 = double2Int x - 1 + | otherwise = double2Int x +{-# INLINE floorD #-} + +instance Bin LogBinD where + type BinValue LogBinD = Double + toIndex !(LogBinD base step _) !x = floorD $ logBase step (x / base) + fromIndex !(LogBinD base step _) !i | i >= 0 = base * step ** (fromIntegral i + 0.5) + | otherwise = -1 / 0 + nBins !(LogBinD _ _ n) = n + {-# INLINE toIndex #-} + +instance IntervalBin LogBinD where + binInterval (LogBinD base step _) i = (x, x*step) where x = base * step ** (fromIntegral i) + +instance Bin1D LogBinD where + lowerLimit (LogBinD lo _ _) = lo + upperLimit (LogBinD lo r n) = lo * r ^ n + unsafeSliceBin i j (LogBinD from step _) = LogBinD (from * step ^ i) step (j-i+1) + +instance VariableBin LogBinD where + binSizeN (LogBinD base step _) n = let x = base * step ^ n in x*step - x + +instance Show LogBinD where + show b = + unlines [ "# LogBinD" + , "# Lo = " ++ show (lowerLimit b) + , "# N = " ++ show (nBins b) + , "# Hi = " ++ show (upperLimit b) + ] +instance Read LogBinD where + readPrec = do + keyword "LogBinD" + liftM3 logBinD (value "Lo") (value "N") (value "Hi") diff --git a/Data/Histogram/Fill.hs b/Data/Histogram/Fill.hs --- a/Data/Histogram/Fill.hs +++ b/Data/Histogram/Fill.hs @@ -13,6 +13,7 @@ , (<<-) , (<<-|) , (< h a b -> (h a b -> h a' b) -> h a' b +h <<-$ f = f h +{-# INLINE (<<-$) #-} + -- | Modify output of histogram. In fact it's same as '<$>' but have opposite fixity (-<<) :: HistBuilder h => (b -> b') -> h a b -> h a b' (-<<) = modifyOut @@ -113,6 +113,7 @@ infixl 5 <<- infixl 5 <<-| infixl 5 < Histogram bin Float forceFloat = id - ----------------------------------------------------------------- --- | Join list of builders into one builders -joinHBuilderMonoidM :: (PrimMonad m, Monoid b) => [HBuilderM m a b] -> HBuilderM m a b -joinHBuilderMonoidM = mconcat -{-# INLINE joinHBuilderMonoidM #-} -{-# DEPRECATED joinHBuilderMonoidM "Use mconcat instead. Will be removed in 0.5" #-} - --- | Join list of builders -joinHBuilderMonoid :: Monoid b => [HBuilder a b] -> HBuilder a b -joinHBuilderMonoid = mconcat -{-# INLINE joinHBuilderMonoid #-} -{-# DEPRECATED joinHBuilderMonoid "Use mconcat instead. Will be removed in 0.5" #-} - -treeHBuilderMonoidM :: (PrimMonad m, Monoid b') => - [HBuilderM m a b -> HBuilderM m a' b'] -> HBuilderM m a b -> HBuilderM m a' b' -treeHBuilderMonoidM fs h = joinHBuilderMonoidM $ map ($ h) fs -{-# INLINE treeHBuilderMonoidM #-} -{-# DEPRECATED treeHBuilderMonoidM "Will be removed in 0.5" #-} - -treeHBuilderMonoid :: Monoid b' => [HBuilder a b -> HBuilder a' b'] -> HBuilder a b -> HBuilder a' b' -treeHBuilderMonoid fs h = joinHBuilderMonoid $ map ($ h) fs -{-# INLINE treeHBuilderMonoid #-} -{-# DEPRECATED treeHBuilderMonoid "Will be removed in 0.5" #-} diff --git a/Data/Histogram/Generic.hs b/Data/Histogram/Generic.hs --- a/Data/Histogram/Generic.hs +++ b/Data/Histogram/Generic.hs @@ -25,7 +25,10 @@ -- ** Convert to other data types , asList , asVector - -- * Slicing histograms + -- * Slicing histogram + , sliceByIx + , sliceByVal + -- * Splitting 2D histograms , sliceX , sliceY -- * Modify histogram @@ -37,13 +40,11 @@ import Control.Applicative ((<$>),(<*>)) import Control.Arrow ((***)) -import Control.Monad (ap, forM_) -import Control.Monad.ST (runST) +import Control.Monad (ap) -import qualified Data.Vector.Generic.Mutable as M import qualified Data.Vector.Generic as G import Data.Typeable (Typeable1(..), Typeable2(..), mkTyConApp, mkTyCon) -import Data.Vector.Generic (Vector) +import Data.Vector.Generic (Vector,(!)) import Text.Read import Data.Histogram.Bin @@ -186,6 +187,17 @@ f2 (x,x') (y,y') = (f x y, f x' y') +sliceByIx :: (Bin1D bin, Vector v a) => Int -> Int -> Histogram v bin a -> Histogram v bin a +sliceByIx i j (Histogram b _ v) = + Histogram (sliceBin i j b) Nothing (G.slice i (j - i + 1) v) + +sliceByVal :: (Bin1D bin, Vector v a) => BinValue bin -> BinValue bin -> Histogram v bin a -> Histogram v bin a +sliceByVal x y h + | inRange b x && inRange b y = sliceByIx (toIndex b x) (toIndex b y) h + | otherwise = error "sliceByVal: Values are out of range" + where + b = bins h + -- | 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] @@ -201,6 +213,4 @@ (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 + mkArray x = G.generate ny (\y -> a ! (y*nx + x)) diff --git a/histogram-fill.cabal b/histogram-fill.cabal --- a/histogram-fill.cabal +++ b/histogram-fill.cabal @@ -1,5 +1,5 @@ Name: histogram-fill -Version: 0.4 +Version: 0.5 Cabal-Version: >= 1.6 License: BSD3 License-File: LICENSE @@ -26,7 +26,15 @@ Data.Histogram.Generic Data.Histogram.Fill Data.Histogram.Bin + Data.Histogram.Bin.Classes + Data.Histogram.Bin.BinI + Data.Histogram.Bin.BinInt + Data.Histogram.Bin.BinEnum + Data.Histogram.Bin.BinF + Data.Histogram.Bin.LogBinD + Data.Histogram.Bin.Bin2D Data.Histogram.Bin.Extra Data.Histogram.ST Other-modules: Data.Histogram.Parse Ghc-options: -O2 -Wall + Ghc-prof-options: -auto-all