foldl-incremental 0.1.1.0 → 0.2.0.0

raw patch · 8 files changed

+664/−55 lines, 8 filesdep +QuickCheckdep +containersdep +deepseqdep −hastachedep ~foldldep ~foldl-incrementalPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, containers, deepseq, histogram-fill, mwc-random, pipes, vector

Dependencies removed: hastache

Dependency ranges changed: foldl, foldl-incremental

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +4/−0
• foldl-incremental.cabal +26/−20
• src/Control/Foldl/Incremental.hs +121/−31
• src/Control/Foldl/Incremental/Histogram.hs +198/−0
• src/Control/Foldl/Incremental/Simple.hs +62/−0
• src/Data/Histogram/Adaptable.hs +126/−0
• src/Data/Histogram/Bin/BinDU.hs +126/−0
• test/test.hs +1/−4

@@ -1,3 +1,7 @@ 0.1.1 --- * Removed Increment from API++0.2.0+---+* inclusion of Histogram folds

@@ -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 +  <https://hackage.haskell.org/package/foldl foldl> 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

@@ -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:++<http://queue.acm.org/detail.cfm?id=2534976 Online Algorithms in High-frequency Trading>+<http://www.google.com.au/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=0CCwQFjAA&url=http%3A%2F%2Fwww-personal.umich.edu%2F~annastef%2Fpapers_Long_ctrl%2FJournalPaperMassGrade_Final.pdf&ei=kq1ZU76IFIPGkAXtqoHACQ&usg=AFQjCNHG7gVfpVoh1g9gjWUcK3Hb22JIlg&sig2=_YxAjKZFo_UEXoFWJafUsw&bvm=bv.65397613,d.dGI  recursive least squares with exponential forgetting>+ -}  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

@@ -0,0 +1,198 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++{-| incremental folds of 'Histogram's from the <https://hackage.haskell.org/package/histogram-fill histogram-fill> 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 #-}

@@ -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+

@@ -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)

@@ -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)

@@ -5,9 +5,6 @@ main = defaultMain tests  tests :: TestTree-tests = testGroup "Tests" [tests']--tests' :: TestTree-tests' = testGroup "tests"+tests = testGroup "incremental tests"   [ TestIncremental.quickTests   ]