diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown --- a/CHANGELOG.markdown +++ b/CHANGELOG.markdown @@ -1,3 +1,7 @@ 0.1.1 --- * Removed Increment from API + +0.2.0 +--- +* inclusion of Histogram folds diff --git a/foldl-incremental.cabal b/foldl-incremental.cabal --- a/foldl-incremental.cabal +++ b/foldl-incremental.cabal @@ -1,5 +1,5 @@ Name: foldl-incremental -Version: 0.1.1.0 +Version: 0.2.0.0 Author: Tony Day Maintainer: tonyday567@gmail.com License: MIT @@ -15,13 +15,13 @@ README.markdown CHANGELOG.markdown Synopsis: incremental folds -Description: This library provides incremental statistical folds based upon the - foldl libray. An incremental statistical fold can be thought of as +Description: Incremental statistical folds based upon the + libray. An incremental statistical fold can be thought of as exponentially-weighting statistics designed to be efficient computations over a Foldable. It supplies "incrementalize" which turns any unary function into a - "Fold". As a reference, \"incrementalize id\" is an exponentially-weighted moving average. + "Fold". As a reference, `incrementalize id` is an exponentially-weighted moving average. Homepage: https://github.com/tonyday567/foldl-incremental Bug-Reports: https://github.com/tonyday567/foldl-incremental/issues Tested-With: GHC==7.6.3 @@ -31,10 +31,18 @@ Library Exposed-Modules: Control.Foldl.Incremental + , Control.Foldl.Incremental.Histogram + , Control.Foldl.Incremental.Simple + , Data.Histogram.Adaptable + , Data.Histogram.Bin.BinDU - Build-Depends: base >= 4 && < 5, - foldl >= 1.0.3 && < 2 - + Build-Depends: base >= 4 && < 5 + , containers + , deepseq + , foldl + , histogram-fill + , vector + Default-Language: Haskell2010 HS-Source-Dirs: src @@ -43,22 +51,20 @@ Main-Is: test.hs Default-Language: Haskell2010 HS-Source-Dirs: src, test - Build-Depends: base >= 4 && < 5, - bytestring >= 0.10.0.2, - foldl >= 1.0.3 && < 2, - tasty >= 0.7 && < 1, - tasty-golden >= 2.2.0.2 && < 3, - tasty-quickcheck >= 0.8, - tasty-hunit >= 0.4.1 && < 5, - foldl-incremental >= 0.1.0.1 + Build-Depends: base >= 4 && < 5, tasty >= 0.7 && < 1, tasty-golden >= 2.2.0.2 && < 3, tasty-quickcheck >= 0.8, tasty-hunit >= 0.4.1 && < 5, foldl-incremental >= 0.2.0, QuickCheck >= 2.7.5 + , bytestring >= 0.10.0.2 + , containers >= 0.5.0.0 + , foldl >= 1.0.3 && < 2 + , histogram-fill >= 0.8.4.1 + , mwc-random >= 0.13.1.1 + , pipes >= 4.1.1 + , vector >= 0.10.0.1 Benchmark bench Type: exitcode-stdio-1.0 Main-Is: bench.hs Default-Language: Haskell2010 HS-Source-Dirs: test - Build-Depends: base >= 4 && < 5, - foldl >= 1.0.3 && < 2, - hastache == 0.5.1, - criterion >= 0.8.0.1, - foldl-incremental >= 0.1.0.1 + Build-Depends: base >= 4 && < 5, criterion >= 0.8.0.1, foldl-incremental >= 0.2.0 + , containers >= 0.5.0.0 + , foldl >= 1.0.3 && < 2 diff --git a/src/Control/Foldl/Incremental.hs b/src/Control/Foldl/Incremental.hs --- a/src/Control/Foldl/Incremental.hs +++ b/src/Control/Foldl/Incremental.hs @@ -23,28 +23,48 @@ The folds represent incremental statistics such as moving averages`. -The stream of moving averages with a `rate` of 0.1 is: +The stream of moving averages with a forgetting `rate` of 0.9 is: ->>> L.scan (incMa 0.1) [1.. 10] +>>> L.scan (ma 0.9) [1..10] +[NaN,1.0,1.5263157894736843,2.070110701107011,2.6312881651642916,3.2097140484969837,3.805217699371904,4.4175932632947745,5.046601250122929,5.691970329383086,6.3533993278762955] or if you just want the moving average at the end. ->>> L.fold (incMa 0.1) [1..10] +>>> L.fold (ma 0.9) [1..10] +6.3533993278762955 +The simple average is obtained via a decay rate of 1.0 (ie no decay) + +>>> L.fold (ma 1.0) [1..10] +5.5 + + +further reading: + + + + -} module Control.Foldl.Incremental ( -- * incrementalize incrementalize -- * common incremental folds - , incMa - , incAbs - , incSq - , incStd + , ma + , absma + , sqma + , std + , cov + , corr + , length + , beta + , alpha + , autocorr ) where -import Control.Applicative ((<$>), (<*>)) -import Control.Foldl (Fold(..)) +import Control.Applicative ((<$>), (<*>)) +import Control.Foldl (Fold(..),premap) +import Prelude hiding (length) -- | An Increment is the incremental state within an exponential moving average fold. data Increment = Increment @@ -53,7 +73,7 @@ , _rate :: {-# UNPACK #-} !Double } deriving (Show) -{-| Incrementalize takes a function and turns it into a `Control.Foldl.Fold` where the step is an Increment iso to the typical step in an exponential moving average calculation. +{-| Incrementalize takes a function and turns it into a `Control.Foldl.Fold` where the step is an Increment similar to the typical step in an exponential moving average calculation. >>> incrementalize id @@ -67,20 +87,36 @@ >>> std r = (\s ss -> sqrt (ss - s**2)) <$> incrementalize id r <*> incrementalize (*2) r -The rate is the parameter regulating the discount of current state and the introduction of the current value. +incrementalize works with any function that produces a double. A correlation fold of a tuple is quite intuitive: +>>> cov r = (\xy xbar ybar -> xy - xbar * ybar) <$> incrementalize (uncurry (*)) r <*> incrementalize fst r <*> incrementalize snd r +>>> corr r = (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov r <*> L.premap fst (std r) <*> L.premap snd (std r) + +The rate is the parameter regulating the discount (or forgetting) of current state and the introduction of the current value. + >>> incrementalize id 1 -tracks the sum/average of an entire Foldable. +tracks the sum/average of an entire Foldable. In other words, prior values are never forgotten. >>> incrementalize id 0 -produces the latest value (ie current state is discounted (or decays) to zero) +produces the latest value (ie current state is discounted (or decays) to zero). In other words, prior values are immediately forgotten. -A exponential moving average with a duration of 10 (the average lag of the values effecting the calculation) is +A exponential moving average with an exponetially-weighted length (duration if its a time series) of 10 (the average lag of the values effecting the calculation) is ->>> incrementalize id (1/10) +>>> incrementalize id (1 - 1/10) +>>> L.fold (length 0.9) [1..100] +9.999734386011127 + +There is no particular reason for different parts to have the same rate. A standard deviation where mean is expected to be static (eg equal to the unconditional sample average) would be: + +>>> std' r = (\s ss -> sqrt (ss - s**2)) <$> incrementalize id 1 <*> incrementalize (*2) r + +and a standard deviation with a prior for the mean (eg ignoring sample averges) would be: + +>>> std'' mean r = incrementalize (\x -> x*2 - mean**2) r + -} incrementalize :: (a -> Double) -> Double -> Fold a Double incrementalize f r = Fold step (Increment 0 0 r) (\(Increment a c _) -> a / c) @@ -88,22 +124,76 @@ step (Increment n d r') n' = Increment (r' * n + f n') (r' * d + 1) r' {-# INLINABLE incrementalize #-} --- | moving average fold -incMa :: Double -> Fold Double Double -incMa = incrementalize id -{-# INLINABLE incMa #-} +-- | incremental average +ma :: Double -> Fold Double Double +ma = incrementalize id +{-# INLINABLE ma #-} --- | moving absolute average -incAbs :: Double -> Fold Double Double -incAbs = incrementalize abs -{-# INLINABLE incAbs #-} +-- | incremental absolute average +absma :: Double -> Fold Double Double +absma = incrementalize abs +{-# INLINABLE absma #-} --- | moving average square -incSq :: Double -> Fold Double Double -incSq = incrementalize (\x -> x*x) -{-# INLINABLE incSq #-} +-- | incremental average square +sqma :: Double -> Fold Double Double +sqma = incrementalize (\x -> x*x) +{-# INLINABLE sqma #-} --- | moving standard deviation -incStd :: Double -> Fold Double Double -incStd rate = (\s ss -> sqrt (ss - s**2)) <$> incMa rate <*> incSq rate -{-# INLINABLE incStd #-} +-- | incremental standard deviation +std :: Double -> Fold Double Double +std rate = (\s ss -> sqrt (ss - s**2)) <$> ma rate <*> sqma rate +{-# INLINABLE std #-} + +-- | incremental covariance +cov :: Double -> Fold (Double, Double) Double +cov r = (\xy xbar ybar -> xy - xbar * ybar) <$> incrementalize (uncurry (*)) r <*> incrementalize fst r <*> incrementalize snd r +{-# INLINABLE cov #-} + +-- | incremental corelation +corr :: Double -> Fold (Double, Double) Double +corr r = (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov r <*> premap fst (std r) <*> premap snd (std r) +{-# INLINABLE corr #-} + +-- | the exponentially weighted length of a rate, which is 1/(1-rate) at infinity +length :: Double -> Fold a Double +length r = Fold step (Increment 0 0 r) (\(Increment _ d _) -> d) + where + step (Increment _ d r') _ = Increment 1 (r' * d + 1) r' +{-# INLINABLE length #-} + +-- | the beta in a simple linear regression of `snd` on `fst` +beta :: Double -> Fold (Double, Double) Double +beta r = (/) <$> cov r <*> premap snd (std r) +{-# INLINABLE beta #-} + +-- | the alpha in a simple linear regression of `snd` on `fst` +alpha :: Double -> Fold (Double, Double) Double +alpha r = (\y b x -> y - b * x) <$> premap fst (ma r) <*> beta r <*> premap snd (ma r) +{-# INLINABLE alpha #-} + +{-| autocorrelation is a slippery concept. This method starts with the concept that there is an underlying random error process (e), and autocorrelation is a process on top of that ie for a one-step correlation relationship. + +value@t = e@t + k * e@t-1 + +where k is the autocorrelation. + +There are thus two decay rates needed: one for the average being considered to be the dependent variable, and one for the decay of the correlation calculation between the most recent value and the moving average. + +>>> L.fold (autoCorr 0 1) + +Would estimate the one-step autocorrelation relationship of the previous value and the current value over the entire sample set. + +-} +autocorr :: Double -> Double -> Fold Double Double +autocorr maR corrR = + case ma maR of + (Fold maStep maBegin maDone) -> + case corr corrR of + (Fold corrStep corrBegin corrDone) -> + let begin = (maBegin, corrBegin) + step (maAcc,corrAcc) a = (maStep maAcc a, + if isNaN (maDone maAcc) + then corrAcc + else corrStep corrAcc (maDone maAcc, a)) + done = corrDone . snd in + Fold step begin done diff --git a/src/Control/Foldl/Incremental/Histogram.hs b/src/Control/Foldl/Incremental/Histogram.hs new file mode 100644 --- /dev/null +++ b/src/Control/Foldl/Incremental/Histogram.hs @@ -0,0 +1,198 @@ +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +{-| incremental folds of 'Histogram's from the library. + +-} + +module Control.Foldl.Incremental.Histogram ( + -- * Incrementalize + incrementalizeHist + , incrementalizeHist2D + -- * Common Histogram Folds + , incHist + , incHist2D + , incAdaptiveHist + ) where + +import Control.Foldl (Fold(..)) +import qualified Control.Foldl as L +import Data.Histogram +import Data.Histogram.Adaptable +import Data.Histogram.Bin.BinDU +import Data.Histogram.Fill +import Data.Vector.Unboxed ((//)) +import qualified Data.Vector.Unboxed as V +import GHC.Float (double2Int) + +-- | FIXME: make Increment a class and IncrementHist an instance. +data IncrementHist a = IncrementHist + { _adder :: Histogram a Double + , _counter :: {-# UNPACK #-} !Double + , _rate :: {-# UNPACK #-} !Double + } + +{-| incrementalizeHist takes a function governing the input to the histogram. + +>>> incrementalizeHist (const 1) + +is the usual boiler-plate meaning of histogram. + +-} +incrementalizeHist :: BinD -> (Double -> Double) -> Double -> L.Fold Double (Histogram BinD Double) +incrementalizeHist b f r = L.Fold step begin done + where + step (IncrementHist x' d r') a = IncrementHist + (Data.Histogram.zip (\x0 y0 -> x0 * r' + y0) x' (mkOne b f a)) + (d * r' + 1) + r' + begin = IncrementHist (fillBuilderVec (mkSimple b) V.empty) 0 r + done (IncrementHist a c _) = + if c /= 0 + then Data.Histogram.map (/c) a + else a + mkOne :: BinD -> (Double -> Double) -> Double -> Histogram BinD Double + mkOne b' f' a = histogram b' d + where + i = toIndex b a + d = V.replicate (nBins b) 0 // [(i,f' a)] + +-- 2D histogram folding +data IncrementHist2D = IncrementHist2D + { _adder2D :: Histogram (Bin2D BinD BinD) Double + , _counter2D :: {-# UNPACK #-} !Double + , _rate2D :: {-# UNPACK #-} !Double + } deriving (Show) + +-- | 2D version +incrementalizeHist2D :: Bin2D BinD BinD -> ((Double,Double) -> Double) -> Double -> L.Fold (Double,Double) (Histogram (Bin2D BinD BinD) Double) +incrementalizeHist2D b f r = L.Fold step begin done + where + step (IncrementHist2D x' d r') a = IncrementHist2D + (Data.Histogram.zip (\x0 y0 -> x0 * r' + y0) x' (mkOne2D b f a)) + (d * r' + 1) + r' + begin = IncrementHist2D (fillBuilderVec (mkSimple b) V.empty) 0 r + done (IncrementHist2D a c _) = + if c /= 0 + then Data.Histogram.map (/c) a + else a + mkOne2D :: Bin2D BinD BinD -> ((Double,Double) -> Double) -> (Double,Double) -> Histogram (Bin2D BinD BinD) Double + mkOne2D b' f' a = histogram b' d + where + i = toIndex b a + d = V.replicate (nBins b) 0 // [(i,f' a)] + +{-| incremental histogram with pre-defined bins + +>>> import Control.Foldl.Incremental +>>> import qualified Control.Foldl as L + +>>> let b = binDn 0 2 12 +>>> L.fold (incHist b 0.9) [1..10] +-} +incHist :: BinD -> Double -> Fold Double (Histogram BinD Double) +incHist b = incrementalizeHist b (const 1) +{-# INLINABLE incHist #-} + +-- | incremental 2D histogram +incHist2D :: Bin2D BinD BinD -> Double -> Fold (Double,Double) (Histogram (Bin2D BinD BinD) Double) +incHist2D b = incrementalizeHist2D b (const 1) + +{-| adaptable histogram fold + +TODO: integrate Histogram.Adaptable upstream + +incHist requires a pre-specified bin, which in turn requires an initial pass over the stream to determine the data ranges. + +For a one pass histogram fold, we require an incremental approach to bin creation, which, in turn, requires some way of creating a histogram from scratch. + +'Data.Histogram.Adaptable' and 'Data.Histogram.Bin.BinDU' is a draft solution to enable a one-pass at histogram creation. + + +This function takes + +- a maximum frequency (thresh) for a bin, which, when triggered causes a bin to be split (at the moment using a uniform distribution assumption which is pretty bad). +- a minimum bin size (grain). bins are further constrained to be multiples of this. + +- a maximum number of bins, which, when triggered, causes bins to be merged. + +>>> L.fold (incAdaptiveHist 0.2 1.0 10 1.0) [1..1000] + +provides a histogram with no bin more than 20% frequency size, with a minimum bin size of 1, with at most 10 bins, and a decay rate of 1.0 + +-} +incAdaptiveHist :: Double -> Double -> Int -> Double -> Fold Double (Histogram BinDU Double) +incAdaptiveHist thresh grain maxBins rate = Fold step begin done + where + step (IncrementHist h d _) a = + let h' = Data.Histogram.map (*rate) (step' h a) in + IncrementHist h' ((d+1)*rate) rate + step' x a = + checkMaxBins maxBins + (maybeSlice thresh grain (checkMaxBins maxBins (addOne grain x a)) a) + where + checkMaxBins max' x' = + if nBins (bins x') > max' + then checkMaxBins max' (mergeSmallest x') + else x' + begin = IncrementHist (histogram (unsafeBinDU V.empty) V.empty) 0 rate + done (IncrementHist h c _) = + if c /= 0 + then Data.Histogram.map (/c) h + else h + +-- helpers +addOne :: Double -> Histogram BinDU Double -> Double -> Histogram BinDU Double +addOne grain x a + | V.length (cuts (bins x)) == 0 = addFirst grain a + | a < lowerLimit (bins x) = addLower grain x a + | a >= upperLimit (bins x) = addUpper grain x a + | otherwise = addMiddle x a + +addMiddle :: (Bin bin, V.Unbox a, Num a) => Histogram bin a -> BinValue bin -> Histogram bin a +addMiddle h a = histogram (bins h) (histData h // [(i, 1 + h `atI` i)]) + where + i = toIndex (bins h) a + +addFirst :: Double -> Double -> Histogram BinDU Double +addFirst grain a = histogram bin (V.fromList [1.0]) + where + a' = roundD grain a + bin = binDU (V.fromList [a' - grain/2, a' + grain/2]) + +addLower :: Double -> Histogram BinDU Double-> Double -> Histogram BinDU Double +addLower grain x a = x1 + where + a' = roundD grain a + x0 = sliceAt x (a' - grain/2) + x1 = addMiddle x0 a + +addUpper :: Double -> Histogram BinDU Double-> Double -> Histogram BinDU Double +addUpper grain x a = x1 + where + a' = roundD grain a + x0 = sliceAt x (a' + grain/2) + x1 = addMiddle x0 a + +maybeSlice :: Double-> Double-> Histogram BinDU Double-> BinValue BinDU-> Histogram BinDU Double +maybeSlice thresh' grain' x' a' = + let freq = x' `atV` a' + i = toIndex (bins x') a' + (l,u) = binInterval (bins x') i + numGrains = round ((u-l) / grain') :: Int + center = l + grain' * fromIntegral (round (fromIntegral numGrains / 2)) + in + if binSizeN (bins x') i > grain' + 3e-11 && + freq / V.sum (histData x') > thresh' + then sliceAt x' center + else x' + +roundD :: Double -> Double -> Double +roundD grain x = if frac > 0.5 then whole + grain else whole + where + whole = floorDD + frac = (x - floorDD) / grain + floorDD = fromIntegral (floorD (x / grain)) * grain + floorD x' | x' < 0 = double2Int x' - 1 + | otherwise = double2Int x' +{-# INLINE roundD #-} diff --git a/src/Control/Foldl/Incremental/Simple.hs b/src/Control/Foldl/Incremental/Simple.hs new file mode 100644 --- /dev/null +++ b/src/Control/Foldl/Incremental/Simple.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE NoMonomorphismRestriction #-} + +{-| Simple moving average style folds -} + +module Control.Foldl.Incremental.Simple ( + -- * incrementalize + incrementalizeSimple + -- * common simple folds + , sma + , sabsma + , ssqma + , sstd + ) where + +import Control.Applicative +import Control.Foldl as L +import Data.Sequence (ViewR(EmptyR, (:>)), (<|)) +import qualified Data.Sequence as Seq + +() :: (Fractional c, Applicative f) => f c -> f c -> f c +() = liftA2 (/) + +(<+>) :: (Fractional c, Applicative f) => f c -> f c -> f c +(<+>) = liftA2 (+) + +maybeSum :: Fold (Maybe Double) (Maybe Double) +maybeSum = Fold (<+>) (Just 0) id + +{-| +Incrementalize takes a function and turns it into a `Control.Foldl.Fold` where the step is an Increment similar to the typical step in a simple moving average calculation. +-} +incrementalizeSimple :: (a -> Double) -> Int -> Fold a Double +incrementalizeSimple f n = L.Fold step begin done + where + begin = Seq.replicate n Nothing + av x = L.fold maybeSum x pure (fromIntegral n) + done x = case av x of + Nothing -> 0/0 + Just x' -> x' + step x a = Just (f a) <| pop x + pop x = case Seq.viewr x of + EmptyR -> x + x'' :> _ -> x'' + +-- | a simple moving average +sma :: Int -> L.Fold Double Double +sma = incrementalizeSimple id + +-- | simple squared moving average +ssqma :: Int -> L.Fold Double Double +ssqma = incrementalizeSimple (**2) + +-- | simple absolute moving average +sabsma :: Int -> L.Fold Double Double +sabsma = incrementalizeSimple abs + +-- | simple standard deviation +sstd :: Int -> Fold Double Double +sstd n = (\s ss -> sqrt (ss - s**2)) <$> sma n <*> ssqma n + diff --git a/src/Data/Histogram/Adaptable.hs b/src/Data/Histogram/Adaptable.hs new file mode 100644 --- /dev/null +++ b/src/Data/Histogram/Adaptable.hs @@ -0,0 +1,126 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE MultiWayIf #-} + +-- Immutable histogram specialised to Unboxed Double, that can add and delete bins +module Data.Histogram.Adaptable ( + -- * Immutable adaptable histograms + HistogramDU + , sliceAt + , insertAt + , mergeAtCut + , smallestCutContiguous + , smallestCutSingle + , mergeSmallest + , mergeSmallestSingle + ) where + +import qualified Data.Vector.Unboxed as V +import Data.Vector.Unboxed ((!),(++)) + +import Data.Histogram +import Data.Histogram.Bin.BinDU + +import Prelude hiding ((++)) + +-- | Immutable Adaptable histogram. +type HistogramDU = Histogram BinDU Double + +sliceAt :: Histogram BinDU Double -> Double -> Histogram BinDU Double +sliceAt h x + | V.length (cuts (bins h)) == 0 = addFirst + | x < lowerLimit (bins h) = addLower + | x > upperLimit (bins h) = addUpper + | otherwise = addMiddle + where + b = bins h + n = nBins b + b' = addCut b x + v = histData h + i = toIndex b x + freq = h `atV` x + r = binInterval b i + size = binSizeN b i + slice0 = freq * (x - fst r)/size + slice1 = freq * (snd r - x)/size + addFirst = histogram (unsafeBinDU (V.fromList [x])) V.empty + addLower = histogram b' (V.singleton 0 ++ v) + addUpper = histogram b' (v ++ V.singleton 0) + addMiddle = histogram b' (V.concat + [ V.take i v + , V.fromList [slice0,slice1] + , V.drop (min (i+1) (n+1)) v + ]) + + +-- | dont interpolate the bin values +insertAt :: Histogram BinDU Double -> Double -> Histogram BinDU Double +insertAt h x + | V.length (cuts (bins h)) == 0 = addFirst + | x < lowerLimit (bins h) = addLower + | x > upperLimit (bins h) = addUpper + | otherwise = addMiddle + where + b = bins h + n = nBins b + b' = addCut b x + v = histData h + i = toIndex b x + addFirst = histogram (unsafeBinDU (V.fromList [x])) V.empty + addLower = histogram b' (V.singleton 0 ++ v) + addUpper = histogram b' (v ++ V.singleton 0) + addMiddle = histogram b' (V.concat + [ V.take (i+1) v + , V.singleton 0 + , V.drop (min (i+1) (n+1)) v + ]) + +mergeAtCut :: HistogramDU -> Int -> HistogramDU +mergeAtCut h i + | i <0 || i>n = error "Data.Histogram.HistogramA': outside index range" + | i == 0 = case h `atI` 0 of + 0 -> histogram (deleteCut (bins h) i) (V.drop 1 v) + _ -> error "Data.Histogram.HistogramA': can't delete outer bin with non-zero frequency" + | i == n = case h `atI` (i-1) of + 0 -> histogram (deleteCut (bins h) i) (V.init v) + _ -> error "Data.Histogram.HistogramA': can't delete outer bin with non-zero frequency" + | otherwise = histogram b' v' + where + n = nBins (bins h) + b' = deleteCut (bins h) i + v = histData h + v' = V.concat + [ V.take (max 0 (i-1)) v + , V.singleton (v ! (i - 1) + v ! i) + , V.drop (min (i+1) n) v + ] + +smallestCutContiguous :: (Bin b, V.Unbox a, Ord a, Num a) => Histogram b a -> Int +smallestCutContiguous h + | v ! 0 == 0 = 0 + | v ! n == 0 = n+1 + | otherwise = 1 + V.minIndex (V.zipWith (+) (V.init v) (V.tail v)) + where + v = histData h + n = nBins (bins h) - 1 + +smallestCutSingle :: (Bin b, V.Unbox a, Ord a, Num a) => Histogram b a -> Int +smallestCutSingle h + | v ! 0 == 0 = 0 + | v ! n == 0 = n+1 + | mini == 0 = 1 + | mini == n = n + | v!(mini-1) < v!(mini+1) = mini + | otherwise = mini+1 + where + v = histData h + n = nBins (bins h) - 1 + mini = V.minIndex v + +mergeSmallest :: HistogramDU -> HistogramDU +mergeSmallest h = mergeAtCut h (smallestCutContiguous h) + +mergeSmallestSingle :: HistogramDU -> HistogramDU +mergeSmallestSingle h = mergeAtCut h (smallestCutSingle h) diff --git a/src/Data/Histogram/Bin/BinDU.hs b/src/Data/Histogram/Bin/BinDU.hs new file mode 100644 --- /dev/null +++ b/src/Data/Histogram/Bin/BinDU.hs @@ -0,0 +1,126 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TypeFamilies #-} + +{- +λ: let bin = binDU (V.fromList [0,1,2]) + +λ: bin +# BinDU cuts +0.0 1.0 2.0 + +λ: addCut bin 3 +# BinDU cuts +0.0 1.0 2.0 3.0 + +λ: deleteCut bin 0 +# BinDU cuts +1.0 2.0 + +-} + +module Data.Histogram.Bin.BinDU ( + -- * Specialized to Double, Unboxed Vectors + BinDU(..) + , binDU + , cuts + , unsafeBinDU + , AdaptableBin(..) + ) where + +import Control.DeepSeq (NFData(..)) +import Data.Data (Data,Typeable) +import Data.Vector.Unboxed (Vector,(!)) +import qualified Data.Vector.Unboxed as VU +import Data.Maybe + +import Data.Histogram.Bin.Classes + +-- | Double bins of unequal sizes. +-- Bins are defined by a vector of cuts marking bounadries between bins (The entire range is continuous. There are n+1 cuts for n bins +-- Cuts are assumed to be in ascending order +-- Specialized on Data.Vector.Unboxed +-- TODO: Generic Vector type. +-- Type paramter: +-- +-- [@v@] type of vector used to define bin cuts + +data BinDU = BinDU !(Vector Double) -- vector of cuts + deriving (Data,Typeable,Eq) + +-- | Create bins unsafely +unsafeBinDU :: Vector Double -- ^ cuts + -> BinDU +unsafeBinDU = BinDU + +binDU :: Vector Double -- ^ cuts + -> BinDU +binDU c + | VU.length c < 2 = error "Data.Histogram.Bin.BinDU.binDU': nonpositive number of bins" + | VU.any (uncurry (>)) (VU.zip (VU.init c) (VU.drop 1 c)) = error "Data.Histogram.Bin.BinDU.binDU': cuts not in ascending order" + | otherwise = BinDU c + +cuts :: BinDU -> Vector Double +cuts (BinDU c) = c + +instance Bin BinDU where + type BinValue BinDU = Double + toIndex (BinDU c) !x = case VU.findIndex (>x) c of + Nothing -> error "Data.Histogram.Bin.BinDU.toIndex: above range" + Just i -> case i of + 0 -> error "Data.Histogram.Bin.BinDU.toIndex: below range" + _ -> i-1 + + fromIndex (BinDU c) !i + | i >= VU.length c - 1 = + error "Data.Histogram.Bin.BinDU.fromIndex: above range" + | otherwise = ((c ! i) + (c ! (i+1)))/2 + + nBins (BinDU c) = if VU.length c < 2 then 0 else VU.length c - 1 + {-# INLINE toIndex #-} + +instance IntervalBin BinDU where + binInterval (BinDU c) i = (c ! i, c ! (i+1)) + +instance Bin1D BinDU where + lowerLimit (BinDU c) = VU.head c + upperLimit (BinDU c) = VU.last c + +instance SliceableBin BinDU where + unsafeSliceBin i j (BinDU c) = BinDU (VU.drop i $ VU.take (j-i) c) + +instance VariableBin BinDU where + binSizeN (BinDU c) !i = c ! (i+1) - c ! i + +-- | Equality is up to 3e-11 (2/3th of digits) +instance BinEq BinDU where + binEq (BinDU c) (BinDU c') + = isNothing (VU.find (\(d,d') -> d - d' > eps * abs d) $ VU.zip c c') + where + eps = 3e-11 + +instance Show BinDU where + show (BinDU c) = "# BinDU cuts\n" ++ concat (fmap showCut $ VU.toList c) ++ "\n\n" + where + showCut x = show x ++ "\t" + +instance NFData BinDU + +-- | Binning algorithms which support adaption. +class Bin b => AdaptableBin b where + -- | delete a bin + deleteCut :: b -> Int -> b + -- | add a new bin + addCut :: b -> Double -> b + +instance AdaptableBin BinDU where + deleteCut (BinDU c) !i + | VU.length c <= 2 = + error "Data.Histogram.Bin.BinDU.deletBin: deleting single bin" + | otherwise = BinDU (VU.take i c VU.++ VU.drop (i+1) c) + + addCut (BinDU c) !x = BinDU (VU.concat [VU.take i c, VU.singleton x, VU.drop i c]) + where + i = fromMaybe (VU.length c) (VU.findIndex (> x) c) diff --git a/test/test.hs b/test/test.hs --- a/test/test.hs +++ b/test/test.hs @@ -5,9 +5,6 @@ main = defaultMain tests tests :: TestTree -tests = testGroup "Tests" [tests'] - -tests' :: TestTree -tests' = testGroup "tests" +tests = testGroup "incremental tests" [ TestIncremental.quickTests ]