module Network.Mux.DeltaQ.TraceStats
  ( step
  , OneWayDeltaQSample (..)
  , constructSample
  , StatsA
  , initialStatsA
  ) where

import           Data.IntMap.Strict (IntMap)
import qualified Data.IntMap.Strict as IM
import           Data.Maybe
import           Data.Word (Word32)

import           Control.Monad.Class.MonadTime.SI
import           Network.Mux.DeltaQ.TraceStatsSupport
import           Network.Mux.DeltaQ.TraceTypes
import           Network.Mux.Types

-- the per observation processing step
step :: RemoteClockModel        -- ^ Remote clock timestamp
     -> Time                    -- ^ Local clock timestamp
     -> Int                     -- ^ the number of octets in the
                                --   observed outcome
     -> StatsA                  -- ^ accumulation state
     -> (Maybe OneWayDeltaQSample, StatsA)
step :: RemoteClockModel
-> Time -> Key -> StatsA -> (Maybe OneWayDeltaQSample, StatsA)
step RemoteClockModel
remoteTS Time
localTS Key
obsSize StatsA
s
 | forall a. Maybe a -> Bool
isNothing (StatsA -> Maybe (Word32, Time)
referenceTimePoint StatsA
s) -- first observation this sample period
   = RemoteClockModel
-> Time -> Key -> StatsA -> (Maybe OneWayDeltaQSample, StatsA)
step RemoteClockModel
remoteTS Time
localTS Key
obsSize
          (StatsA
s { referenceTimePoint :: Maybe (Word32, Time)
referenceTimePoint = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$! (RemoteClockModel -> Word32
unRemoteClockModel RemoteClockModel
remoteTS, Time
localTS)
             , nextSampleAt :: Time
nextSampleAt       = DiffTime
sampleInterval DiffTime -> Time -> Time
`addTime` Time
localTS
             , timeLastObs :: Time
timeLastObs        = Time
localTS -- for single observation in sample case
             })

 | Time
localTS forall a. Ord a => a -> a -> Bool
<= StatsA -> Time
nextSampleAt StatsA
s    -- still in a sample period
  = let refTimePoint :: (Word32, Time)
refTimePoint = case StatsA -> Maybe (Word32, Time)
referenceTimePoint StatsA
s of
          Just (Word32, Time)
a  -> (Word32, Time)
a
          Maybe (Word32, Time)
Nothing -> forall a. HasCallStack => [Char] -> a
error [Char]
"step: missing referenceTimePoint"
        transitTime :: SISec
transitTime       = (Word32, Time) -> RemoteClockModel -> Time -> SISec
calcTransitTime (Word32, Time)
refTimePoint RemoteClockModel
remoteTS Time
localTS
    in (forall a. Maybe a
Nothing, StatsA -> Time -> Key -> SISec -> StatsA
recordObservation StatsA
s Time
localTS Key
obsSize SISec
transitTime)
 | Bool
otherwise                    -- need to start next sample period
  = let sample :: OneWayDeltaQSample
sample  = StatsA -> OneWayDeltaQSample
constructSample StatsA
s
        (Maybe OneWayDeltaQSample
_, StatsA
s') = RemoteClockModel
-> Time -> Key -> StatsA -> (Maybe OneWayDeltaQSample, StatsA)
step RemoteClockModel
remoteTS Time
localTS Key
obsSize StatsA
initialStatsA
    in (forall a. a -> Maybe a
Just OneWayDeltaQSample
sample, StatsA
s')

-- Calculate the transit time by transforming the remotely reported
-- emit time into local clock domain then calculating differences.
calcTransitTime :: (Word32, Time)
                -> RemoteClockModel
                -> Time
                -> SISec
calcTransitTime :: (Word32, Time) -> RemoteClockModel -> Time -> SISec
calcTransitTime (Word32
remoteRefTS, Time
localRefTS) RemoteClockModel
remoteTS' Time
localTS
  = let remoteTS :: Word32
remoteTS
          = RemoteClockModel -> Word32
unRemoteClockModel RemoteClockModel
remoteTS'
        remoteClockDiffAsTimeDiff :: Word32 -> DiffTime
remoteClockDiffAsTimeDiff
          = (DiffTime
remoteClockPrecision forall a. Num a => a -> a -> a
*) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Fractional a => Rational -> a
fromRational forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (Integral a, Num b) => a -> b
fromIntegral
        correctedEmitTime :: Time
correctedEmitTime
          | Word32
remoteTS forall a. Ord a => a -> a -> Bool
>= Word32
remoteRefTS
          = (Word32 -> DiffTime
remoteClockDiffAsTimeDiff forall a b. (a -> b) -> a -> b
$ Word32
remoteTS forall a. Num a => a -> a -> a
- Word32
remoteRefTS)
            DiffTime -> Time -> Time
`addTime` Time
localRefTS
          | Bool
otherwise -- wrap has occurred
          = (Word32 -> DiffTime
remoteClockDiffAsTimeDiff forall a b. (a -> b) -> a -> b
$ forall a. Bounded a => a
maxBound forall a. Num a => a -> a -> a
- (Word32
remoteRefTS forall a. Num a => a -> a -> a
- Word32
remoteTS))
            DiffTime -> Time -> Time
`addTime` Time
localRefTS
    in Float -> SISec
S forall a b. (a -> b) -> a -> b
$! forall a. Fractional a => Rational -> a
fromRational (forall a. Real a => a -> Rational
toRational (Time
localTS Time -> Time -> DiffTime
`diffTime` Time
correctedEmitTime))

recordObservation :: StatsA -> Time -> Int -> SISec -> StatsA
recordObservation :: StatsA -> Time -> Key -> SISec -> StatsA
recordObservation StatsA
s Time
obsTime Key
obsSize SISec
transitTime
  = let f :: Maybe PerSizeRecord -> Maybe PerSizeRecord
f Maybe PerSizeRecord
Nothing  = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$! SISec -> PerSizeRecord
makePerSizeRecord SISec
transitTime
        f (Just PerSizeRecord
a) = forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$! SISec -> PerSizeRecord
makePerSizeRecord SISec
transitTime forall a. Semigroup a => a -> a -> a
<> PerSizeRecord
a
    in StatsA
s { timeLastObs :: Time
timeLastObs     = Time
obsTime
         , numObservations :: Key
numObservations = forall a. Enum a => a -> a
succ (StatsA -> Key
numObservations StatsA
s)
         , observables :: IntMap PerSizeRecord
observables     = forall a. (Maybe a -> Maybe a) -> Key -> IntMap a -> IntMap a
IM.alter Maybe PerSizeRecord -> Maybe PerSizeRecord
f Key
obsSize (StatsA -> IntMap PerSizeRecord
observables StatsA
s)
         }

-- This might benefit from some strictness analysis, what are the
-- likely usage patterns?, do we want a single use collapse the
-- collective set of thunks or not?
--
-- There is the issue of "bias" (in its statistical meaning) here. The
-- approach here is pragmatic, we are going to use the uncorrected
-- sample standard deviation here as it has a defined solution for a
-- single sample.
--
-- Given that the consumer of this statistic also has access to the
-- population size, they could reconstruct the underlying measures and
-- take it from there.
--
-- We return `NaN` for the appropriate statistics when the population
-- is empty.
constructSample :: StatsA -> OneWayDeltaQSample
constructSample :: StatsA -> OneWayDeltaQSample
constructSample StatsA
sa = OneWaySample
  { duration :: Double
duration       = forall a. Fractional a => Rational -> a
fromRational forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Real a => a -> Rational
toRational forall a b. (a -> b) -> a -> b
$
                     forall b a. b -> (a -> b) -> Maybe a -> b
maybe DiffTime
0 (\(Word32
_,Time
a) -> StatsA -> Time
timeLastObs StatsA
sa Time -> Time -> DiffTime
`diffTime` Time
a)
                             (StatsA -> Maybe (Word32, Time)
referenceTimePoint StatsA
sa)
  , sumPackets :: Key
sumPackets     = Key
population
  , sumTotalSDU :: Key
sumTotalSDU    = Key
totalSDUOctets
  , estDeltaQS :: Double
estDeltaQS     = Double -> Double
normCheck Double
dQSEst
  , estDeltaQVMean :: Double
estDeltaQVMean = Double -> Double
normCheck forall a b. (a -> b) -> a -> b
$ Double
vSum forall a. Fractional a => a -> a -> a
/ Double
pop
  , estDeltaQVVar :: Double
estDeltaQVVar  = Double -> Double
normCheck forall a b. (a -> b) -> a -> b
$ (Double
vSum2 forall a. Num a => a -> a -> a
- Double
vSum forall a. Num a => a -> a -> a
* Double
vSum forall a. Fractional a => a -> a -> a
/ Double
pop) forall a. Fractional a => a -> a -> a
/ Double
pop
  , estR :: Double
estR           = Double -> Double
normCheck Double
rEst
  , sizeDist :: [Char]
sizeDist       = forall a. Show a => a -> [Char]
show [ (Key
a,PerSizeRecord -> Key
count PerSizeRecord
b, let S Float
mtt = PerSizeRecord -> SISec
minTransitTime PerSizeRecord
b in Float
mtt)
                          | (Key
a, PerSizeRecord
b) <- forall a. IntMap a -> [(Key, a)]
IM.toAscList (StatsA -> IntMap PerSizeRecord
observables StatsA
sa)
                          , PerSizeRecord -> Key
count PerSizeRecord
b forall a. Ord a => a -> a -> Bool
> Key
0]
  }
  where
    -- the sample population size
    population :: Key
population = StatsA -> Key
numObservations StatsA
sa
    pop :: Double
pop        = forall a b. (Integral a, Num b) => a -> b
fromIntegral Key
population
    -- Handle the empty population condition
    normCheck :: Double -> Double
normCheck Double
x
      | forall a. IntMap a -> Key
IM.size (StatsA -> IntMap PerSizeRecord
observables StatsA
sa) forall a. Ord a => a -> a -> Bool
> Key
1 = Double
x
      | Bool
otherwise                    = Double
nan
    -- gather the data for the G,S estimations
    (Key
totalSDUOctets, [(Key, SISec)]
minSRev)
      = forall a b. (Key -> a -> b -> b) -> b -> IntMap a -> b
IM.foldrWithKey Key
-> PerSizeRecord -> (Key, [(Key, SISec)]) -> (Key, [(Key, SISec)])
accum (Key
0, []) forall a b. (a -> b) -> a -> b
$ StatsA -> IntMap PerSizeRecord
observables StatsA
sa
    accum :: Key
-> PerSizeRecord -> (Key, [(Key, SISec)]) -> (Key, [(Key, SISec)])
accum Key
nOctets PerSizeRecord
psr (Key
sumSize, [(Key, SISec)]
minS)
      = ( Key
sumSize forall a. Num a => a -> a -> a
+ (PerSizeRecord -> Key
count PerSizeRecord
psr) forall a. Num a => a -> a -> a
* Key
nOctets
        , (Key
nOctets, PerSizeRecord -> SISec
minTransitTime PerSizeRecord
psr) forall a. a -> [a] -> [a]
: [(Key, SISec)]
minS)
    -- fit a line to get the G,S estimation
    (Double
dQGEst, Double
dQSEst, Double
rEst) = [(Key, SISec)] -> (Double, Double, Double)
estimateGS [(Key, SISec)]
minSRev
    -- normalise all the observations
    normalisedObservations :: IntMap PerSizeRecord
normalisedObservations
      = let norm :: a -> SISec
norm a
n = Float -> SISec
S forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Fractional a => Rational -> a
fromRational forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Real a => a -> Rational
toRational
                     forall a b. (a -> b) -> a -> b
$ Double
dQGEst forall a. Num a => a -> a -> a
+ (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n) forall a. Num a => a -> a -> a
* Double
dQSEst
        in forall a b. (Key -> a -> b) -> IntMap a -> IntMap b
IM.mapWithKey (\Key
k -> SISec -> PerSizeRecord -> PerSizeRecord
normalisePSR (forall {a}. Integral a => a -> SISec
norm Key
k)) (StatsA -> IntMap PerSizeRecord
observables StatsA
sa)
    -- calculate the total population V stats
    (Double
vSum, Double
vSum2)
      = let S  Float
v  = SISec
vSum'
            S2 Float
v2 = SISec2
vSum2'
        in (forall a. Fractional a => Rational -> a
fromRational forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Real a => a -> Rational
toRational forall a b. (a -> b) -> a -> b
$ Float
v, forall a. Fractional a => Rational -> a
fromRational forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Real a => a -> Rational
toRational forall a b. (a -> b) -> a -> b
$ Float
v2)
    (SISec
vSum', SISec2
vSum2')
      = forall a b. (a -> b -> b) -> b -> IntMap a -> b
IM.foldr PerSizeRecord -> (SISec, SISec2) -> (SISec, SISec2)
vCalc (SISec
0,SISec2
0) IntMap PerSizeRecord
normalisedObservations
    vCalc :: PerSizeRecord -> (SISec, SISec2) -> (SISec, SISec2)
vCalc PerSizeRecord
psr (SISec
x, SISec2
x2)
      = (SISec
x forall a. Num a => a -> a -> a
+ PerSizeRecord -> SISec
sumTransitTime PerSizeRecord
psr, SISec2
x2 forall a. Num a => a -> a -> a
+ PerSizeRecord -> SISec2
sumTransitTimeSq PerSizeRecord
psr)

-- | One way measurement for interval. Note that the fields are lazy
--   here so that only calculation necessary to satisfy strictness of
--   use occurs.
data OneWayDeltaQSample = OneWaySample
  { OneWayDeltaQSample -> Double
duration       :: Double -- SI Seconds of activity captured
  , OneWayDeltaQSample -> Key
sumPackets     :: Int
  , OneWayDeltaQSample -> Key
sumTotalSDU    :: Int
  , OneWayDeltaQSample -> Double
estDeltaQS     :: Double -- octets per second
  , OneWayDeltaQSample -> Double
estDeltaQVMean :: Double -- SI Seconds
  , OneWayDeltaQSample -> Double
estDeltaQVVar  :: Double
  , OneWayDeltaQSample -> Double
estR           :: Double -- R estimate
  , OneWayDeltaQSample -> [Char]
sizeDist       :: String -- temporary to show size distribution
  }

-- | Statistics accumulator. Strict evaluation used to keep the memory
--   footprint strictly bounded.
data StatsA = StatsA
  { -- per sample
    StatsA -> Maybe (Word32, Time)
referenceTimePoint :: !(Maybe (Word32, Time))
  , StatsA -> Time
nextSampleAt       :: !Time
    -- per observation
  , StatsA -> Key
numObservations    :: !Int
  , StatsA -> Time
timeLastObs        :: !Time
  , StatsA -> IntMap PerSizeRecord
observables        :: !(IntMap PerSizeRecord)
  }

-- This _may_ not be the best representation, but it does appear to be
-- an adequate one. There are known issues with numerical stability for
-- this representation approach in floating point arithmetic where the
-- values being measured are "large" and the variability in the sampled
-- population is "small" (i.e the inherent rounding effect of floating
-- point arithmetic has an effect).
--
-- This is very unlikely to cause an issue here as:
--
--   a) the modulo model of the RemoteClockModel (and hence the
--   establishment of a clock reference offset for a given sample)
--   means that we are only ever recording differences - thus any
--   absolute clock differences get factored out.
--
--   b) the transit delay for a measurement will be small, (probably
--   not even credible) worst case ~10^3 / 10^4 seconds, the finite
--   mantissa of IEEE754
--   (https://en.wikipedia.org/wiki/IEEE_754#Representation_and_encoding_in_memory)
--   even for 32 bit representation (24bits / 7.22 decimal digits)
--   represents an bound on the inherent measured population
--   variability.
--
--   c) longer term clock drift is covered here by the re-establishing
--   of the clock reference offsets every sampling period. Given a
--   reasonable sampling period (say 10 seconds) clock drift ( <
--   100ppm) can't amount to a significant error over such a period.
--
-- To conclude, reasonable model of delay is < 1second, the precision
-- of delay measurement is 10^-6 - this all fits nicely within a
-- IEEE754 32bit representation with its 7.22 decimal digit
-- mantissa. Haskell `Float`s are adequate for this purpose.

data PerSizeRecord = PSR
  { PerSizeRecord -> SISec
minTransitTime   :: !SISec
  , PerSizeRecord -> Key
count            :: !Int
  , PerSizeRecord -> SISec
sumTransitTime   :: !SISec
  , PerSizeRecord -> SISec2
sumTransitTimeSq :: !SISec2
  }

instance Semigroup PerSizeRecord where
  PerSizeRecord
a <> :: PerSizeRecord -> PerSizeRecord -> PerSizeRecord
<> PerSizeRecord
b = PSR { minTransitTime :: SISec
minTransitTime   = (PerSizeRecord -> SISec
minTransitTime PerSizeRecord
a) forall a. Ord a => a -> a -> a
`min` (PerSizeRecord -> SISec
minTransitTime PerSizeRecord
b)
               , count :: Key
count            = PerSizeRecord -> Key
count PerSizeRecord
a forall a. Num a => a -> a -> a
+ PerSizeRecord -> Key
count PerSizeRecord
b
               , sumTransitTime :: SISec
sumTransitTime   = PerSizeRecord -> SISec
sumTransitTime PerSizeRecord
a forall a. Num a => a -> a -> a
+ PerSizeRecord -> SISec
sumTransitTime PerSizeRecord
b
               , sumTransitTimeSq :: SISec2
sumTransitTimeSq = PerSizeRecord -> SISec2
sumTransitTimeSq PerSizeRecord
a forall a. Num a => a -> a -> a
+ PerSizeRecord -> SISec2
sumTransitTimeSq PerSizeRecord
b
               }

-- | Normalise given the calculated G,S for the size
normalisePSR :: SISec -> PerSizeRecord -> PerSizeRecord
normalisePSR :: SISec -> PerSizeRecord -> PerSizeRecord
normalisePSR SISec
norm PerSizeRecord
psr
  = let adj :: SISec
adj  = (forall a b. (Integral a, Num b) => a -> b
fromIntegral (PerSizeRecord -> Key
count PerSizeRecord
psr) forall a. Num a => a -> a -> a
* SISec
norm)
        stt' :: SISec
stt' = (PerSizeRecord -> SISec
sumTransitTime PerSizeRecord
psr) forall a. Num a => a -> a -> a
- SISec
adj
        ttt :: SISec -> SISec -> SISec2
ttt (S Float
a) (S Float
b)
             = Float -> SISec2
S2 forall a b. (a -> b) -> a -> b
$ Float
a forall a. Num a => a -> a -> a
* Float
b
    in PerSizeRecord
psr { minTransitTime :: SISec
minTransitTime   = PerSizeRecord -> SISec
minTransitTime   PerSizeRecord
psr forall a. Num a => a -> a -> a
- SISec
norm
           , sumTransitTime :: SISec
sumTransitTime   = SISec
stt'
           , sumTransitTimeSq :: SISec2
sumTransitTimeSq = PerSizeRecord -> SISec2
sumTransitTimeSq PerSizeRecord
psr
                                forall a. Num a => a -> a -> a
- SISec
norm SISec -> SISec -> SISec2
`ttt` (SISec
2 forall a. Num a => a -> a -> a
* SISec
stt' forall a. Num a => a -> a -> a
+ SISec
norm)
           }


-- | Initial StatsA
initialStatsA :: StatsA
initialStatsA :: StatsA
initialStatsA = StatsA
  { referenceTimePoint :: Maybe (Word32, Time)
referenceTimePoint = forall a. Maybe a
Nothing
  , nextSampleAt :: Time
nextSampleAt       = Time
noTime
  , numObservations :: Key
numObservations    = Key
0
  , timeLastObs :: Time
timeLastObs        = Time
noTime
  , observables :: IntMap PerSizeRecord
observables        = forall a. IntMap a
IM.empty
  }
  where
    noTime :: Time
noTime = DiffTime -> Time
Time DiffTime
0

makePerSizeRecord :: SISec -> PerSizeRecord
makePerSizeRecord :: SISec -> PerSizeRecord
makePerSizeRecord SISec
tt = PSR
  { minTransitTime :: SISec
minTransitTime   = SISec
tt
  , count :: Key
count            = Key
1
  , sumTransitTime :: SISec
sumTransitTime   = SISec
tt
  , sumTransitTimeSq :: SISec2
sumTransitTimeSq = SISec -> SISec2
squareSISec SISec
tt
  }

-- May want to make this a configuration variable

-- NOTE this interval must be less than the wrap around time of the
-- `RemoteClockModel`. The remote clock model has a precision of
-- `remoteClockPrecision`.
sampleInterval :: DiffTime
sampleInterval :: DiffTime
sampleInterval = DiffTime -> DiffTime
check DiffTime
10
  where
    check :: DiffTime -> DiffTime
check DiffTime
n
     | DiffTime
n forall a. Ord a => a -> a -> Bool
> DiffTime
0 Bool -> Bool -> Bool
&& DiffTime
n forall a. Ord a => a -> a -> Bool
< DiffTime
wrapInterval
       = DiffTime
n
     | Bool
otherwise
       = forall a. HasCallStack => [Char] -> a
error [Char]
"Infeasible sampleInterval"
    wrapInterval :: DiffTime
wrapInterval
      = DiffTime
remoteClockPrecision forall a. Num a => a -> a -> a
* (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ RemoteClockModel -> Word32
unRemoteClockModel forall a. Bounded a => a
maxBound)

nan :: Double
nan :: Double
nan = Double
0forall a. Fractional a => a -> a -> a
/Double
0