{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts,
MultiParamTypeClasses, TypeFamilies, CPP #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
module Numeric.Sum (
Summation(..)
, sumVector
, KBNSum(..)
, kbn
, KB2Sum(..)
, kb2
, KahanSum(..)
, kahan
, pairwiseSum
) where
import Control.Arrow ((***))
import Control.DeepSeq (NFData(..))
import Data.Bits (shiftR)
import Data.Data (Typeable, Data)
import Data.Monoid (Monoid(..))
#if MIN_VERSION_base(4,9,0)
import Data.Semigroup (Semigroup(..))
#endif
import Data.Vector.Generic (Vector(..), foldl')
import Control.Monad (liftM)
import Data.Vector.Generic.Mutable (MVector(..))
import qualified Data.Foldable as F
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Generic.Mutable as GM
import qualified Data.Vector.Unboxed as U
class Summation s where
zero :: s
add :: s -> Double -> s
sum :: (F.Foldable f) => (s -> Double) -> f Double -> Double
sum f = f . F.foldl' add zero
{-# INLINE sum #-}
instance Summation Double where
zero = 0
add = (+)
data KahanSum = KahanSum {-# UNPACK #-} !Double {-# UNPACK #-} !Double
deriving (Eq, Show, Typeable, Data)
instance U.Unbox KahanSum
newtype instance U.MVector s KahanSum = MV_KahanSum (U.MVector s (Double, Double))
instance MVector U.MVector KahanSum where
{-# INLINE GM.basicLength #-}
{-# INLINE GM.basicUnsafeSlice #-}
{-# INLINE basicOverlaps #-}
{-# INLINE basicUnsafeNew #-}
{-# INLINE basicInitialize #-}
{-# INLINE basicUnsafeReplicate #-}
{-# INLINE basicUnsafeRead #-}
{-# INLINE basicUnsafeWrite #-}
{-# INLINE basicClear #-}
{-# INLINE basicSet #-}
{-# INLINE GM.basicUnsafeCopy #-}
{-# INLINE basicUnsafeMove #-}
{-# INLINE basicUnsafeGrow #-}
basicLength (MV_KahanSum mvec) = GM.basicLength mvec
basicUnsafeSlice idx len (MV_KahanSum mvec) = MV_KahanSum (GM.basicUnsafeSlice idx len mvec)
basicOverlaps (MV_KahanSum mvec) (MV_KahanSum mvec') = basicOverlaps mvec mvec'
basicUnsafeNew len = MV_KahanSum `liftM` basicUnsafeNew len
basicInitialize (MV_KahanSum mvec) = basicInitialize mvec
basicUnsafeReplicate len val = MV_KahanSum `liftM` basicUnsafeReplicate len ((\ (KahanSum a b) -> (a, b)) val)
basicUnsafeRead (MV_KahanSum mvec) idx = (\ (a, b) -> KahanSum a b) `liftM` basicUnsafeRead mvec idx
basicUnsafeWrite (MV_KahanSum mvec) idx val = basicUnsafeWrite mvec idx ((\ (KahanSum a b) -> (a, b)) val)
basicClear (MV_KahanSum mvec) = basicClear mvec
basicSet (MV_KahanSum mvec) val = basicSet mvec ((\ (KahanSum a b) -> (a, b)) val)
basicUnsafeCopy (MV_KahanSum mvec) (MV_KahanSum mvec') = GM.basicUnsafeCopy mvec mvec'
basicUnsafeMove (MV_KahanSum mvec) (MV_KahanSum mvec') = basicUnsafeMove mvec mvec'
basicUnsafeGrow (MV_KahanSum mvec) len = MV_KahanSum `liftM` basicUnsafeGrow mvec len
newtype instance U.Vector KahanSum = V_KahanSum (U.Vector (Double, Double))
instance Vector U.Vector KahanSum where
{-# INLINE basicUnsafeFreeze #-}
{-# INLINE basicUnsafeThaw #-}
{-# INLINE G.basicLength #-}
{-# INLINE G.basicUnsafeSlice #-}
{-# INLINE basicUnsafeIndexM #-}
{-# INLINE G.basicUnsafeCopy #-}
{-# INLINE elemseq #-}
basicUnsafeFreeze (MV_KahanSum mvec) = V_KahanSum `liftM` basicUnsafeFreeze mvec
basicUnsafeThaw (V_KahanSum vec) = MV_KahanSum `liftM` basicUnsafeThaw vec
basicLength (V_KahanSum vec) = G.basicLength vec
basicUnsafeSlice idx len (V_KahanSum vec) = V_KahanSum (G.basicUnsafeSlice idx len vec)
basicUnsafeIndexM (V_KahanSum vec) idx = (\ (a, b) -> KahanSum a b) `liftM` basicUnsafeIndexM vec idx
basicUnsafeCopy (MV_KahanSum mvec) (V_KahanSum vec) = G.basicUnsafeCopy mvec vec
elemseq (V_KahanSum vec) val = elemseq vec ((\ (KahanSum a b) -> (a, b)) val)
instance Summation KahanSum where
zero = KahanSum 0 0
add = kahanAdd
instance NFData KahanSum where
rnf !_ = ()
instance Monoid KahanSum where
mempty = zero
s `mappend` KahanSum s' _ = add s s'
#if MIN_VERSION_base(4,9,0)
instance Semigroup KahanSum where
(<>) = mappend
#endif
kahanAdd :: KahanSum -> Double -> KahanSum
kahanAdd (KahanSum sum c) x = KahanSum sum' c'
where sum' = sum + y
c' = (sum' - sum) - y
y = x - c
kahan :: KahanSum -> Double
kahan (KahanSum sum _) = sum
data KBNSum = KBNSum {-# UNPACK #-} !Double {-# UNPACK #-} !Double
deriving (Eq, Show, Typeable, Data)
instance U.Unbox KBNSum
newtype instance U.MVector s KBNSum = MV_KBNSum (U.MVector s (Double, Double))
instance MVector U.MVector KBNSum where
{-# INLINE GM.basicLength #-}
{-# INLINE GM.basicUnsafeSlice #-}
{-# INLINE basicOverlaps #-}
{-# INLINE basicUnsafeNew #-}
{-# INLINE basicInitialize #-}
{-# INLINE basicUnsafeReplicate #-}
{-# INLINE basicUnsafeRead #-}
{-# INLINE basicUnsafeWrite #-}
{-# INLINE basicClear #-}
{-# INLINE basicSet #-}
{-# INLINE GM.basicUnsafeCopy #-}
{-# INLINE basicUnsafeMove #-}
{-# INLINE basicUnsafeGrow #-}
basicLength (MV_KBNSum mvec) = GM.basicLength mvec
basicUnsafeSlice idx len (MV_KBNSum mvec) = MV_KBNSum (GM.basicUnsafeSlice idx len mvec)
basicOverlaps (MV_KBNSum mvec) (MV_KBNSum mvec') = basicOverlaps mvec mvec'
basicUnsafeNew len = MV_KBNSum `liftM` basicUnsafeNew len
basicInitialize (MV_KBNSum mvec) = basicInitialize mvec
basicUnsafeReplicate len val = MV_KBNSum `liftM` basicUnsafeReplicate len ((\ (KBNSum a b) -> (a, b)) val)
basicUnsafeRead (MV_KBNSum mvec) idx = (\ (a, b) -> KBNSum a b) `liftM` basicUnsafeRead mvec idx
basicUnsafeWrite (MV_KBNSum mvec) idx val = basicUnsafeWrite mvec idx ((\ (KBNSum a b) -> (a, b)) val)
basicClear (MV_KBNSum mvec) = basicClear mvec
basicSet (MV_KBNSum mvec) val = basicSet mvec ((\ (KBNSum a b) -> (a, b)) val)
basicUnsafeCopy (MV_KBNSum mvec) (MV_KBNSum mvec') = GM.basicUnsafeCopy mvec mvec'
basicUnsafeMove (MV_KBNSum mvec) (MV_KBNSum mvec') = basicUnsafeMove mvec mvec'
basicUnsafeGrow (MV_KBNSum mvec) len = MV_KBNSum `liftM` basicUnsafeGrow mvec len
newtype instance U.Vector KBNSum = V_KBNSum (U.Vector (Double, Double))
instance Vector U.Vector KBNSum where
{-# INLINE basicUnsafeFreeze #-}
{-# INLINE basicUnsafeThaw #-}
{-# INLINE G.basicLength #-}
{-# INLINE G.basicUnsafeSlice #-}
{-# INLINE basicUnsafeIndexM #-}
{-# INLINE G.basicUnsafeCopy #-}
{-# INLINE elemseq #-}
basicUnsafeFreeze (MV_KBNSum mvec) = V_KBNSum `liftM` basicUnsafeFreeze mvec
basicUnsafeThaw (V_KBNSum vec) = MV_KBNSum `liftM` basicUnsafeThaw vec
basicLength (V_KBNSum vec) = G.basicLength vec
basicUnsafeSlice idx len (V_KBNSum vec) = V_KBNSum (G.basicUnsafeSlice idx len vec)
basicUnsafeIndexM (V_KBNSum vec) idx = (\ (a, b) -> KBNSum a b) `liftM` basicUnsafeIndexM vec idx
basicUnsafeCopy (MV_KBNSum mvec) (V_KBNSum vec) = G.basicUnsafeCopy mvec vec
elemseq (V_KBNSum vec) val = elemseq vec ((\ (KBNSum a b) -> (a, b)) val)
instance Summation KBNSum where
zero = KBNSum 0 0
add = kbnAdd
instance NFData KBNSum where
rnf !_ = ()
instance Monoid KBNSum where
mempty = zero
s `mappend` KBNSum s' c' = add (add s s') c'
#if MIN_VERSION_base(4,9,0)
instance Semigroup KBNSum where
(<>) = mappend
#endif
kbnAdd :: KBNSum -> Double -> KBNSum
kbnAdd (KBNSum sum c) x = KBNSum sum' c'
where c' | abs sum >= abs x = c + ((sum - sum') + x)
| otherwise = c + ((x - sum') + sum)
sum' = sum + x
kbn :: KBNSum -> Double
kbn (KBNSum sum c) = sum + c
data KB2Sum = KB2Sum {-# UNPACK #-} !Double
{-# UNPACK #-} !Double
{-# UNPACK #-} !Double
deriving (Eq, Show, Typeable, Data)
instance U.Unbox KB2Sum
newtype instance U.MVector s KB2Sum = MV_KB2Sum (U.MVector s (Double, Double, Double))
instance MVector U.MVector KB2Sum where
{-# INLINE GM.basicLength #-}
{-# INLINE GM.basicUnsafeSlice #-}
{-# INLINE basicOverlaps #-}
{-# INLINE basicUnsafeNew #-}
{-# INLINE basicInitialize #-}
{-# INLINE basicUnsafeReplicate #-}
{-# INLINE basicUnsafeRead #-}
{-# INLINE basicUnsafeWrite #-}
{-# INLINE basicClear #-}
{-# INLINE basicSet #-}
{-# INLINE GM.basicUnsafeCopy #-}
{-# INLINE basicUnsafeMove #-}
{-# INLINE basicUnsafeGrow #-}
basicLength (MV_KB2Sum mvec) = GM.basicLength mvec
basicUnsafeSlice idx len (MV_KB2Sum mvec) = MV_KB2Sum (GM.basicUnsafeSlice idx len mvec)
basicOverlaps (MV_KB2Sum mvec) (MV_KB2Sum mvec') = basicOverlaps mvec mvec'
basicUnsafeNew len = MV_KB2Sum `liftM` basicUnsafeNew len
basicInitialize (MV_KB2Sum mvec) = basicInitialize mvec
basicUnsafeReplicate len val = MV_KB2Sum `liftM` basicUnsafeReplicate len ((\ (KB2Sum a b c) -> (a, b, c)) val)
basicUnsafeRead (MV_KB2Sum mvec) idx = (\ (a, b, c) -> KB2Sum a b c) `liftM` basicUnsafeRead mvec idx
basicUnsafeWrite (MV_KB2Sum mvec) idx val = basicUnsafeWrite mvec idx ((\ (KB2Sum a b c) -> (a, b, c)) val)
basicClear (MV_KB2Sum mvec) = basicClear mvec
basicSet (MV_KB2Sum mvec) val = basicSet mvec ((\ (KB2Sum a b c) -> (a, b, c)) val)
basicUnsafeCopy (MV_KB2Sum mvec) (MV_KB2Sum mvec') = GM.basicUnsafeCopy mvec mvec'
basicUnsafeMove (MV_KB2Sum mvec) (MV_KB2Sum mvec') = basicUnsafeMove mvec mvec'
basicUnsafeGrow (MV_KB2Sum mvec) len = MV_KB2Sum `liftM` basicUnsafeGrow mvec len
newtype instance U.Vector KB2Sum = V_KB2Sum (U.Vector (Double, Double, Double))
instance Vector U.Vector KB2Sum where
{-# INLINE basicUnsafeFreeze #-}
{-# INLINE basicUnsafeThaw #-}
{-# INLINE G.basicLength #-}
{-# INLINE G.basicUnsafeSlice #-}
{-# INLINE basicUnsafeIndexM #-}
{-# INLINE G.basicUnsafeCopy #-}
{-# INLINE elemseq #-}
basicUnsafeFreeze (MV_KB2Sum mvec) = V_KB2Sum `liftM` basicUnsafeFreeze mvec
basicUnsafeThaw (V_KB2Sum vec) = MV_KB2Sum `liftM` basicUnsafeThaw vec
basicLength (V_KB2Sum vec) = G.basicLength vec
basicUnsafeSlice idx len (V_KB2Sum vec) = V_KB2Sum (G.basicUnsafeSlice idx len vec)
basicUnsafeIndexM (V_KB2Sum vec) idx = (\ (a, b, c) -> KB2Sum a b c) `liftM` basicUnsafeIndexM vec idx
basicUnsafeCopy (MV_KB2Sum mvec) (V_KB2Sum vec) = G.basicUnsafeCopy mvec vec
elemseq (V_KB2Sum vec) val = elemseq vec ((\ (KB2Sum a b c) -> (a, b, c)) val)
instance Summation KB2Sum where
zero = KB2Sum 0 0 0
add = kb2Add
instance NFData KB2Sum where
rnf !_ = ()
instance Monoid KB2Sum where
mempty = zero
s `mappend` KB2Sum s' c' cc' = add (add (add s s') c') cc'
#if MIN_VERSION_base(4,9,0)
instance Semigroup KB2Sum where
(<>) = mappend
#endif
kb2Add :: KB2Sum -> Double -> KB2Sum
kb2Add (KB2Sum sum c cc) x = KB2Sum sum' c' cc'
where sum' = sum + x
c' = c + k
cc' | abs c >= abs k = cc + ((c - c') + k)
| otherwise = cc + ((k - c') + c)
k | abs sum >= abs x = (sum - sum') + x
| otherwise = (x - sum') + sum
kb2 :: KB2Sum -> Double
kb2 (KB2Sum sum c cc) = sum + c + cc
sumVector :: (Vector v Double, Summation s) =>
(s -> Double) -> v Double -> Double
sumVector f = f . foldl' add zero
{-# INLINE sumVector #-}
pairwiseSum :: (Vector v Double) => v Double -> Double
pairwiseSum v
| len <= 256 = G.sum v
| otherwise = uncurry (+) . (pairwiseSum *** pairwiseSum) .
G.splitAt (len `shiftR` 1) $ v
where len = G.length v
{-# SPECIALIZE pairwiseSum :: V.Vector Double -> Double #-}
{-# SPECIALIZE pairwiseSum :: U.Vector Double -> Double #-}