module Gauge.Analysis
(
Outliers(..)
, OutlierEffect(..)
, OutlierVariance(..)
, SampleAnalysis(..)
, analyseSample
, scale
, analyseMean
, countOutliers
, classifyOutliers
, noteOutliers
, outlierVariance
, resolveAccessors
, validateAccessors
, regress
) where
import Data.Monoid
import Control.Arrow (second)
import Control.Monad (unless, when)
import Gauge.IO.Printf (note, prolix)
import Gauge.Measurement (secs, threshold)
import Gauge.Monad (Gauge, getGen, getOverhead, askConfig, gaugeIO)
import Gauge.Types
import Data.Int (Int64)
import Data.Maybe (fromJust)
import Statistics.Function (sort)
import Statistics.Quantile (weightedAvg, Sorted(..))
import Statistics.Regression (bootstrapRegress, olsRegress)
import Statistics.Resampling (Estimator(..),resample)
import Statistics.Sample (mean)
import Statistics.Sample.KernelDensity (kde)
import Statistics.Types (Sample)
import System.Random.MWC (GenIO)
import qualified Data.List as List
import qualified Data.Map as Map
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Statistics.Resampling.Bootstrap as B
import qualified Statistics.Types as B
import Prelude
classifyOutliers :: Sample -> Outliers
classifyOutliers sa = U.foldl' ((. outlier) . mappend) mempty ssa
where outlier e = Outliers {
samplesSeen = 1
, lowSevere = if e <= loS && e < hiM then 1 else 0
, lowMild = if e > loS && e <= loM then 1 else 0
, highMild = if e >= hiM && e < hiS then 1 else 0
, highSevere = if e >= hiS && e > loM then 1 else 0
}
!loS = q1 (iqr * 3)
!loM = q1 (iqr * 1.5)
!hiM = q3 + (iqr * 1.5)
!hiS = q3 + (iqr * 3)
q1 = weightedAvg 1 4 (Sorted ssa)
q3 = weightedAvg 3 4 (Sorted ssa)
ssa = sort sa
iqr = q3 q1
outlierVariance
:: B.Estimate B.ConfInt Double
-> B.Estimate B.ConfInt Double
-> Double
-> OutlierVariance
outlierVariance µ σ a = OutlierVariance effect desc varOutMin
where
( effect, desc ) | varOutMin < 0.01 = (Unaffected, "no")
| varOutMin < 0.1 = (Slight, "slight")
| varOutMin < 0.5 = (Moderate, "moderate")
| otherwise = (Severe, "severe")
varOutMin = (minBy varOut 1 (minBy cMax 0 µgMin)) / σb2
varOut c = (ac / a) * (σb2 ac * σg2) where ac = a c
σb = B.estPoint σ
µa = B.estPoint µ / a
µgMin = µa / 2
σg = min (µgMin / 4) (σb / sqrt a)
σg2 = σg * σg
σb2 = σb * σb
minBy f q r = min (f q) (f r)
cMax x = fromIntegral (floor (2 * k0 / (k1 + sqrt det)) :: Int)
where
k1 = σb2 a * σg2 + ad
k0 = a * ad
ad = a * d
d = k * k where k = µa x
det = k1 * k1 4 * σg2 * k0
countOutliers :: Outliers -> Int64
countOutliers (Outliers _ a b c d) = a + b + c + d
analyseMean :: Sample
-> Int
-> Gauge Double
analyseMean a iters = do
let µ = mean a
_ <- note "mean is %s (%d iterations)\n" (secs µ) iters
noteOutliers . classifyOutliers $ a
return µ
scale :: Double
-> SampleAnalysis -> SampleAnalysis
scale f s@SampleAnalysis{..} = s {
anMean = B.scale f anMean
, anStdDev = B.scale f anStdDev
}
analyseSample :: Int
-> String
-> V.Vector Measured
-> Gauge (Either String Report)
analyseSample i name meas = do
Config{..} <- askConfig
overhead <- getOverhead
let ests = [Mean,StdDev]
stime = measure (measTime . rescale) .
G.filter ((>= threshold) . measTime) . G.map fixTime .
G.tail $ meas
fixTime m = m { measTime = measTime m overhead / 2 }
n = G.length meas
s = G.length stime
_ <- prolix "bootstrapping with %d of %d samples (%d%%)\n" s n ((s * 100) `quot` n)
gen <- getGen
ers <- (sequence <$>) . mapM (\(ps,r) -> regress gen ps r meas) $ ((["iters"],"time"):regressions)
case ers of
Left err -> pure $ Left err
Right rs -> do
resamps <- gaugeIO $ resample gen ests resamples stime
let [estMean,estStdDev] = B.bootstrapBCA confInterval stime resamps
ov = outlierVariance estMean estStdDev (fromIntegral n)
an = SampleAnalysis
{ anRegress = rs
, anOverhead = overhead
, anMean = estMean
, anStdDev = estStdDev
, anOutlierVar = ov
}
return $ Right $ Report
{ reportNumber = i
, reportName = name
, reportKeys = measureKeys
, reportMeasured = meas
, reportAnalysis = an
, reportOutliers = classifyOutliers stime
, reportKDEs = [uncurry (KDE "time") (kde 128 stime)]
}
regress :: GenIO
-> [String]
-> String
-> V.Vector Measured
-> Gauge (Either String Regression)
regress gen predNames respName meas
| G.null meas = pure $ Left "no measurements"
| otherwise = case validateAccessors predNames respName of
Left err -> pure $ Left err
Right accs -> do
let unmeasured = [n | (n, Nothing) <- map (second ($ G.head meas)) accs]
if not (null unmeasured)
then pure $ Left $ "no data available for " ++ renderNames unmeasured
else do
let (r:ps) = map ((`measure` meas) . (fromJust .) . snd) accs
Config{..} <- askConfig
(coeffs,r2) <- gaugeIO $ bootstrapRegress gen resamples confInterval olsRegress ps r
pure $ Right $ Regression
{ regResponder = respName
, regCoeffs = Map.fromList (zip (predNames ++ ["y"]) (G.toList coeffs))
, regRSquare = r2
}
singleton :: [a] -> Bool
singleton [_] = True
singleton _ = False
resolveAccessors :: [String]
-> Either String [(String, Measured -> Maybe Double)]
resolveAccessors names =
case unresolved of
[] -> Right [(n, a) | (n, Just (a,_)) <- accessors]
_ -> Left $ "unknown metric " ++ renderNames unresolved
where
unresolved = [n | (n, Nothing) <- accessors]
accessors = flip map names $ \n -> (n, Map.lookup n measureAccessors)
validateAccessors :: [String]
-> String
-> Either String [(String, Measured -> Maybe Double)]
validateAccessors predNames respName = do
when (null predNames) $
Left "no predictors specified"
let names = respName:predNames
dups = map head . filter (not . singleton) .
List.group . List.sort $ names
unless (null dups) $
Left $ "duplicated metric " ++ renderNames dups
resolveAccessors names
renderNames :: [String] -> String
renderNames = List.intercalate ", " . map show
noteOutliers :: Outliers -> Gauge ()
noteOutliers o = do
let frac n = (100::Double) * fromIntegral n / fromIntegral (samplesSeen o)
check :: Int64 -> Double -> String -> Gauge ()
check k t d = when (frac k > t) $
note " %d (%.1g%%) %s\n" k (frac k) d
outCount = countOutliers o
when (outCount > 0) $ do
_ <- note "found %d outliers among %d samples (%.1g%%)\n"
outCount (samplesSeen o) (frac outCount)
check (lowSevere o) 0 "low severe"
check (lowMild o) 1 "low mild"
check (highMild o) 1 "high mild"
check (highSevere o) 0 "high severe"