// timeseries analysis utilities

// Data generation

sub cpu(from, to) {
  demo cdn metrics 'cpu' -from from -to to -every :m:
    -nhosts 4 -dos .5 -dos_dur :5m: -ripple .3 -cpu_alpha 0.8
  | filter host ~ 'sjc*'
}

sub daily_cpu(from, to) {
  demo cdn metrics 'cpu' -from from -to to -every :10m:
    -nhosts 4 -dos .3 -dos_dur :15m: -cpu_alpha 0.8 -daily 4
  | filter host ~ 'sjc*'
}

sub bumpy_cpu(from, to) {
  demo cdn metrics 'cpu' -from from -to to -every :5m:
    -nhosts 4 -dos .3 -dos_dur :15s: -daily 3
  | filter host ~ 'sjc*'
}

// exponentially weighted moving avarage smoother.
// alpha [0..1] is the smoothing factor.
reducer ewma(fld, alpha) {
    var yhat;
    var initial = true;

    function update() {
        if (initial) {
            yhat = *fld;
            initial = false;
        } else {
            yhat = alpha * *fld + (1 - alpha) * yhat;
        }
    }
    function result() {
        return yhat;
    }
    function reset() {
        initial = true;
    }
}

// Holt-Winters forecaster, which models an unobserved
// level and trend component in a noisy signal, along with
// optional "seasonal" effects for periodic daily variations and weekly variations
// (scaling factors applied to the current level). the result is a smoothed
// version of the signal which may be used to forecast future values.
//
// Y is the point field to smooth.
// level, trend, daily, weekly are smoothing factors, numbers ranging [0..1] or null.
// They determine how quickly the smoother adjusts to new values as they arrive.
// Setting a factor to 0 freezes the feature at its initial estimate,
// while setting a factor to 1 causes that feature to depend solely
// on the most recent point. Setting a factor to null removes that feature
// from the model (except for level, which must be [0..1])
//
// smooth assumes arriving points are equally spaced in time.
// smooth returns the next-step forecast of the model, an array containing:
// level, trend, daily multiplier, weekly multiplier
// (level is the estimate of the next series value, trend the estimate of the
//  next series slope, etc).
reducer smooth(Y, level, trend, daily, weekly) {
    var lv, tr, lv_1, tr_1;
    var initial = true;

    function update() {
        if (initial) {
            lv = *Y;
            tr = 0;
            lv_1 = *Y;
            tr_1 = 0;
            initial = false;
        } else {
            lv = level * *Y + (1 - level) * (lv_1 + tr_1);
            if (trend != null) {
                tr = trend * (lv - lv_1) + (1 - trend) * tr_1;
            }
            lv_1 = lv;
            tr_1 = tr;
        }
    }
    function result() {
        return [lv +  tr, tr];
    }
    function reset() {
        initial = true;
    }
}
// smoother/estimator for standard deviation of a time series
// alpha [0..1] is the smoothing factor.
reducer smooth_sd(Y, alpha) {
    var sigma2, last;
    var initial = true;

    function update() {
        if (initial) {
            sigma2 = 0;
            last = *Y;
            initial = false;
        } else {
            sigma2 = alpha * (*Y - last) * (*Y - last) + (1 - alpha) * sigma2;
            last = *Y;
        }
    }
    function result() {
        return Math.sqrt(sigma2);
    }
    function reset() {
        initial = true;
    }
}
sub forecastPI(Y, HI, LOW, STDERR, z, alpha, level, trend, daily, weekly) {
    put state = smooth(Y, level, trend, daily, weekly)
    | put yhat = state[0], tr=state[1]
    | put err2 = (*Y - yhat) * (*Y - yhat)
    | put *STDERR = smooth(err2, alpha, null, null, null)[0]
    | put *HI = yhat + z * *STDERR
    | put *LOW = yhat - z * *STDERR
}

sub trend_anomalies() {
}

//cpu -from :2014-01-01: -to :2014-01-01T01:20:00:
//daily_cpu -from :2014-01-01: -to :2014-01-02T01:20:00:
bumpy_cpu -from :2014-01-01T13:00:00:  -to :2014-01-01T21:00:00:
| forecastPI -Y "value" -HI "hi" -LOW "low" -STDERR "stderr"
             -z 2 -alpha .8 -level .8 -trend .8 -daily null -weekly null
             | remove low, hi, state, yhat
| split low, hi, stderr, value, yhat, err2,tr
| @timechart