/* tsad4.cpp - Time Series Anomaly Detector ver. 4.

(C) 2004, Matt Mahoney.  This is free software distributed under terms
of the GNU General Public License, http://www.gnu.org/licenses/gpl.txt

tsad4 assigns anomaly scores to time series data.  It reads N samples
training of a time series and records the path formed by the low-pass
filtered data and its first and second derivatives in M dimensional space, 
scaled to fit a unit M-cube.  The anomaly score of the remaining test data is
the square of the distance to the closest point on the path.  The
path is represented by a piecewise linear approximation with K segments.
The data may be undersampled by reading only every S samples.

Usage: dir/b | tsad4 N C T K M S P > result

For example: dir/b tek*.csv | tsad4 1000 2 5 20 3 12 1 1

The input is a list of training files and test files containing time
series data formatted as columns of numbers separated by other characters
such as spaces, tabs, or commas.  Parameters are integers with the
following meanings:

N = number of training samples.  All the rest are test samples.

C = Input column number.  The leftmost column is 1.  All other
    columns are ignored.  Default is 1.

T = Filtering time constant (in samples).  The input data x[i] is smoothed
    to y[i] by low pass filtering: y[i] = ((T-1) * y[i-1] + x[i]) / T.
    Default is 5.

K = Number of vertices to approximate the training path.  Default is N.

M = Number of dimensions (data plus M-1 derivatives), range 1-4.

R = State transitions to test.  The state is the closest path segment.
    Normally the state advances when the next segment is closer.
    Setting R=1 allows skipping one or moving back one.  R=2 tests random
    segments and jumps if closer.  R=3 does both tests.  R=4 tests
    every segment.  R+8 is as above but approximates segments by moving
    them 1/4 the distance to the removed vertices.

S = Read every S samples.  Default is 1 (don't skip).  N and T refer
    to the number of samples after skipping S-1 of every S.

P = number of paths in the training data.  Each path is N samples for
    a total of N*P training points.
    Default is 1.  If P > 1 then the anomaly score is the square of the
    distance from the smallest M-dimensional box bounding the P nearest
    points to the test point in each of the paths.  Each path is
    approximated to K vertices.

tsad4 computes a model of N training points in M dimensions by piecewise
linear approximation of the path formed by the input and its first M-1
derivatives (scaled to a unit M-cube) to K vertices.  The anomaly score
of a test point is the square of the distance to this approximated path.

The file tsad4.lst is created containing a table of M+2 columns as follows:

  state a x dx ddx dddx

For the first K rows, each (x,dx,ddx,dddx) point is a vertex in
the segmented approximation, state is the vertex number (0 to K-1),
a is the sample number (0 to N-1) corresponding to this vertex,
x is the input (filtered), and dx, ddx, and dddx are the first 3
derivatives of x.

For the test data (after N samples), (x,dx,ddx,dddx) is the test point,
state is the current state, or closest path segment to the test point,
and a is the anonaly score.

Each point (x,dx,ddx,dddx) is low pass filtered twice with time
constant T samples.  The anomaly score is filtered once.  This
is shown schematically in Fig. 1.

                    x                   dx                  ddx
input -> LP -> LP --+--> D -> LP -> LP --+--> D -> LP -> LP -+
                    |                    |                   |
                    V (col. 1)           V (col. 3)          V (col. 5)
                +-----------------------------------------------+       a
                |               Anomaly Detector                |-> LP -->
                +-----------------------------------------------+ (col. 2)

        Fig. 1.  Anomaly detection configuration and output for M = 3.

In Fig. 1, D is a differentiator, such that D(t) = x[t] - x[t-1], and
LP is a low pass filter, such that LP(t) = ((T-1)LP(t-1) + x[t]) / T,
where x[t] is the t'th input.

Standard output prints a table of maximum and total anomaly scores for
each test file.  The maximum and total are after filtering once as shown
in Fig. 1.  This filter state is reset to 0 after each file.  No data
is output for the training portion of the file.



To compile: g++ tsad4.cpp

*/

#include <cstdio>
#include <cstdlib>
#include <cctype>
#include <cstring>
#include <cmath>
#include <string>
#include <vector>
#include <ctime>
#include <algorithm>
#include <cassert>
using namespace std;

// functoid elapsedTime() returns time in seconds since the previous call
class ElapsedTime {
  clock_t t;  // time of last call
public:
  ElapsedTime() {
    t=clock();
  }
  double operator()() {  // seconds since previous call
    clock_t c=clock();
    double r=double(c-t)/CLOCKS_PER_SEC;
    t=c;
    return r;
  }
} elapsedTime;

// Read line from f into s until EOF
bool getline(FILE* f, string& s) {
  s="";
  if (!f)
    return false;
  int c;
  while((c=getc(f))!=EOF && c!='\n')
    s+=c;
  return c!=EOF || s!="";
}

// Parse numbers from s, return number in column col >= 1
double split(const string& s, int col) {
  int b=0, e=0, i=1;  // begin, end of number, column
  while (true) {
    while (b<int(s.size()) && !strchr("0123456789-.", s[b]))
      ++b;
    e=b;
    while (e<int(s.size()) && strchr("0123456789-.+eE", s[e]))
      ++e;
    if (e>b) {
      if (i==col)
        return atof(s.substr(b, e-b).c_str());
    }
    else
      return 0;
    b=e;
    ++i;
  }
}

// Square
inline double sqr(double x) {return x*x;}

// Globals
int M = 3;  // dimensions
int R = 4;  // flags: 1=test 4 segments, 2=random, 4=all, 8=adjust
const int MMAX = 4;  // max M
double scale[MMAX] = {1,1,1,1};  // Model scale

// Square of the distance from point A to point B in M dimensions
double dist2(const double a[], const double b[]) {
  double d=0;
  for (int i=0; i<M; ++i)
    d+=sqr((a[i]-b[i])*scale[i]);
  return d;
}

// Square of the distance from line segment AB to point C in M dimensions
// Put the closest point on AB in p[] (unless p == 0)
double dist2(const double a[], const double b[], const double c[],
    double *p = 0) {
  double t=0;  // Closest point on line AB to C where 0=A, 1=B
  double d=dist2(a, b);  // Length of AB squared
  if (d>0) {
    double n=0;
    for (int i=0; i<M; ++i)
      n+=(a[i]-c[i])*(a[i]-b[i])*sqr(scale[i]);
    t=n/d;
  }
  if (t<0) t=0;  // Clip to endpoints of AB
  if (t>1) t=1;
  d=0;  // Result
  for (int i=0; i<M; ++i) {
    const double pi=(1-t)*a[i]+t*b[i];
    if (p)
      p[i]=pi;
    d+=sqr((pi-c[i])*scale[i]);
  }
  return d;
}

/* Shrinker.shrink(x, N, K) replaces the path of N M-dimensional vertices
in x[0..N-1][0..M-1] with a piecewise linear approximation of K vertices,
placing the result back in x[0..K-1][0..M-1].  The vertices are assumed
to be scaled by scale[0..M-1] to fit a unit M-dimensional cube.
2 <= K <= N.  0 < M <= MMAX.

The approximation is performed by iteratively removing the vertex with
the lowest error in the range 1..N-2 (keeping the endpoints) until there
are K vertices remaining.  If A,B,C are 3 successive vertices in a path,
then error(B) = |AC|d^2 where |AC| is the length of the line segment
between A and C, and d^2 is the square of the Euclidean distance from B
to the closest point on the line segment AC (possibly an endpoint).

Approximation is performed in O(N log N) time by storing the vertices
in a doubly linked heap.  Each vertex (except the first and last, which
are outside the heap) has an index into x, an error, and pointers to
the previous and next vertex in the path.  The heap is a binary tree
initially with N-2 elements such that a parent node has lower error
than its children.  The heap is an array of elements 1..N-2 (with the
fixed endpoints at 0 and N-1), such that the root node is at 1 and the
children of node i are 2i and 2i+1.  When vertices in the heap are moved,
the linked list pointers to them are updated.  The algorithm is as follows:

  construct doubly linked list of N vertices in heap[0..N-1]
  set heap size = N-2 (elements 1 to N-1)
  compute errors for each vertex in heap
  for (i = N/2..1) heapify(i)
  repeat N-K times
    remove root node from heap and linked list
    for (i = removed.prev, removed.next)
      compute_error(i)
      if error has increased then heapify(i)
      else decrease_key(i)
  traverse K nodes of linked list and store in x

heapify(i) ensures that node i has lower error than its 2 children.
If not, then node i is swapped with the smaller of its children and the
updated child is heapified.

decrease_key(i) ensures that node i has higher error than its parent
(except for the root, i = 1).  If not, then it is swapped with the parent
and decrease_key(i/2) is called on the parent.

compute_error(x) sets the error on a vertex as described above.

exch(i,j) swaps nodes i and j, updating the linked list pointers of
their neighbors.

print(f) prints the model to FILE* f.

*/

class Shrinker {
  struct Vertex {
    double err;  // heap key, smallest at root
    int i;  // doubly linked list key and index into x
    Vertex *prev, *next; // linked list pointers, prev->i < i < next->i
    Vertex(): err(0), i(0), prev(0), next(0) {}
    void compute_error(double x[][MMAX]) {
      if (prev && next)
        err = dist2(x[prev->i], x[next->i], x[i])
            * sqrt(dist2(x[prev->i], x[next->i]));
    }
  };
  Vertex* heap;  // [N] with list endpoints at 0 and N-1
  int size;  // Heap size, N-2 shrinking down to K-2

  // swap heap[i] and heap[j] maintaining links
  void exch(int i, int j) {
    if (i==j)
      return;
    if (heap[j].next == &heap[i])
      swap(j, i);
    if (heap[i].next == &heap[j]) {
      int h=heap[i].prev - heap;  // h - i - j - k becomes h - j - i - k
      int k=heap[j].next - heap;
      swap(heap[i], heap[j]);
      heap[h].next = &heap[j];
      heap[j].next = &heap[i];
      heap[i].next = &heap[k];
      heap[k].prev = &heap[i];
      heap[i].prev = &heap[j];
      heap[j].prev = &heap[h];
    }
    else {  // i and j are not adjacent
      swap(heap[i], heap[j]);
      heap[i].prev->next = heap[i].next->prev = &heap[i];
      heap[j].prev->next = heap[j].next->prev = &heap[j];
    }
  }

  static inline int left(int i) {return i+i;}
  static inline int right(int i) {return i+i+1;}
  static inline int parent(int i) {return (unsigned int)i>>1;}

  // Ensure heap[i].err is smaller than children
  void heapify(int i) {
    assert(i>0 && i<=size);
    while (true) {
      int smallest=i;
      if (left(i)<=size && heap[left(i)].err < heap[smallest].err)
        smallest = left(i);
      if (right(i)<=size && heap[right(i)].err < heap[smallest].err)
        smallest = right(i);
      if (i==smallest)
        break;
      else {
        exch(i, smallest);
        i = smallest;
      }
    }
  }

  // Ensure heap[i].err is larger than parent
  void decrease_key(int i) {
    assert(i>0 && i<=size);
    while (i>1) {
      if (heap[i].err < heap[parent(i)].err) {
        exch(i, parent(i));
        i = parent(i);
      }
      else
        break;
    }
  }

  // Endure heap[i].err is at the right level
  void fix(int i) {
    heapify(i);
    decrease_key(i);
  }
public:
  void shrink(double x[][MMAX], int N, int K);
  void print(FILE* f, double x[][MMAX], int K) const;
};

// Print the model
void Shrinker::print(FILE* f, double x[][MMAX], int K) const {
  Vertex *p=heap;
  for (int i=0; i<K; ++i) {
    fprintf(f, "%d %d", i, p->i);
    for (int j=0; j<M; ++j)
      fprintf(f, " %g", x[i][j]);
    fprintf(f, "\n");
    p=p->next;
  }
}

void Shrinker::shrink(double x[][MMAX], int N, int K) {

  // Initialize doubly linked heap
  heap = new Vertex[N];
  size = N-2;
  for (int i=0; i<N; ++i)
    heap[i].i = i;
  for (int i=0; i<N-1; ++i) {  // make linked list
    heap[i].next = &heap[i+1];
    heap[i+1].prev = &heap[i];
  }
  for (int i=1; i<N-1; ++i)  // compute errors except at endpoints
    heap[i].compute_error(x);
  for (int i=size/2; i>=1; --i)  // make heap at 1..N-2
    heapify(i);

  // Remove vertices
  while (size>0 && size>K-2) {
    exch(1, size--);  // pop smallest err from heap
    Vertex& v = heap[size+1];

    // Move neighbors toward removed vertex
    if (R&8) {
      double p[MMAX];  // closest point on path to removed vertex
      dist2(x[v.prev->i], x[v.next->i], x[v.i], p);
      for (int i=0; i<M; ++i) {
        const double d=(x[v.i][i]-p[i])*0.25;
        x[v.prev->i][i]+=d;
        x[v.next->i][i]+=d;
      }
    }

    // Unlink vertex from list
    v.prev->next = v.next;
    v.next->prev = v.prev;
    heapify(1);

    // Update err of neighbors of removed or adjusted vertices
    if (v.prev > heap) {
      v.prev->compute_error(x);
      fix(v.prev - heap);
      if ((R&8) && v.prev->prev > heap) {
        v.prev->prev->compute_error(x);
        fix(v.prev->prev - heap);
      }
    }
    if (v.next <= heap+size) {
      v.next->compute_error(x);
      fix(v.next - heap);
      if ((R&8) && v.next->next <= heap+size) {
        v.next->next->compute_error(x);
        fix(v.next->next - heap);
      }
    }
  }

  // Save back to x
  Vertex *p=heap;
  for (int i=0; i<K; ++i) {
    assert(p);
    int j=p->i;
    assert(j>=i && j<N);
    if (j>i)
      memmove(x[i], x[j], M*sizeof(double));
    p=p->next;
  }
  assert(!p);
}

int main(int argc, char **argv) {

  // Print help message
  if (argc<2) {
    printf("tsad4 v5 time series anomaly detector\n"
      "(C) 2004, Matt Mahoney.  Free under GPL: http://www.gnu.org/licenses/gpl.txt\n"
      "\n"
      "Usage: dir/b | tsad4 N C T K M R S P\n"
      "\n"
      "Input: training and test files containing columns of numbers\n"
      "N = number of training samples per path; the rest is test\n"
      "C = column number to use (starting at 1); default is 1\n"
      "T = filter time constant; default is 5 samples\n"
      "K = number of vertices in model approximation; default is N\n"
      "M = number of dimensions (M-1 derivatives); range 1-4; default is 3\n"
      "R = state transitions to test: 1=prev and 2nd, 2=random, 4=all (default)\n"
      "    R+8 as above but approximates segments 1/4 distance to removed points\n"
      "S = subsample rate; default is 1\n"
      "P = number of training paths, default is 1\n"
      "\n"
      "Output is maximum and total anomaly score per file\n"
      "tsad4.out is in M+2 columns: state, anomaly score, data, derivatives\n"
      "tsad4.mod lists the model vertices: time, M coordinates\n");
    return 1;
  }
  elapsedTime();  // start timer

  // Input parameters
  const int N = atoi(argv[1]);                    // training set size
  const int C = argc>2 ? atoi(argv[2]) : 1;       // column
  const double T = argc>3 ? atof(argv[3]) : 5.0;  // filter time constant
  const int K = argc>4 ? atoi(argv[4]) : N;       // model segments
  M = argc>5 ? atoi(argv[5]) : 3;                 // dimensions
  if (M>MMAX) M=MMAX;
  if (M<1) M=1;
  R = argc>6 ? atoi(argv[6]) : 4;                 // transitions to test
  const int S = argc>7 ? atoi(argv[7]) : 1;       // sample rate
  const int P = argc>8 ? atoi(argv[8]) : 1;
  printf("N=%d C=%d T=%1.2f K=%d M=%d R=%d S=%d P=%d\n",
    N, C, T, K, M, R, S, P);

  // Open output
  FILE* out=fopen("tsad4.out", "w");
  if (!out) {
    perror("tsad4.out");
    exit(1);
  }

  // Read files
  string filename;
  double v[MMAX]={0};  // current value and M-1 derivatives
  double fv[MMAX]={0}; // filtered v
  double ffv[MMAX]={0}; // filtered fv
  double (*x)[MMAX] = new double[N*P][MMAX];  // training data, derivatives
  memset(x, 0, N*MMAX*P*sizeof(double));
  double amax=0, asum=0, fa=0;  // max, sum, filtered anomaly score
  const double DT=1-1/T;
  int i=0;  // sample number
  vector<int> state(P);  // closest path segment to test data
  vector<double> lo(M), hi(M);  // bounding box of nearest paths to ffv
  double* p=new double[M];  // closest point on j'th path
  while (getline(stdin, filename)) {
    FILE* f=fopen(filename.c_str(), "r");
    if (!f)
      perror(filename.c_str());
    amax=asum=fa=0;
    fill(state.begin(), state.end(), 0);
    string s;
    while (getline(f, s)) {
      ++i;

      // Filter and take derivatives
      v[0]=split(s, C);
      for (int j=0; j<M; ++j) {
        fv[j]=fv[j]*DT+v[j]*(1-DT);
        double ffv1=ffv[j];
        ffv[j]=ffv[j]*DT+fv[j]*(1-DT);
        if (j<M-1)
          v[j+1]=ffv[j]-ffv1;
      }

      // Store training data
      if (i<=N*P) {
        for (int j=0; j<M; ++j)
          x[i-1][j]=ffv[j];
      }

      // At end of training scale path to unit cube and approximate
      if (i==N*P) {
        printf("%1.2f sec. to read %d training samples\n",
            elapsedTime(), N*P);
        for (int j=0; j<M; ++j) {
          double xmin=x[0][j];
          double xmax=x[0][j];
          for (int k=1; k<N; ++k) {
            if (x[k][j]<xmin)
              xmin=x[k][j];
            if (x[k][j]>xmax)
              xmax=x[k][j];
          }
          if (xmax>xmin)
            scale[j]=1/(xmax-xmin);
          else
            scale[j]=0;
          printf("scale[%d] = %f (%f to %f)\n",
            j, scale[j], xmin, xmax);
        }
        printf("%1.2f sec. to scale %d dimensions to unit cube\n",
          elapsedTime(), M);

        // Approximate N*P vertices with K*P vertices
        FILE* mod=fopen("tsad4.mod", "w");
        Shrinker sh;
        for (int j=0; j<P; ++j) {
          sh.shrink(x+j*N, N, K);
          if (mod) sh.print(mod, x+j*N, K);
        }

        // Print the model
        printf("%1.2f sec. to approximate with %d vertices\n",
          elapsedTime(), K);
        if (mod) {
          fclose(mod);
          printf("%d vertex model saved to tsad4.mod\n", K);
        }
        else
          perror("tsad4.mod");
      }
      if (i==N*P+1) {
        printf(
          "\nTest file                             "
          "Anomaly Scores:  Maximum      Total\n");
      }

      // Compute anomaly score
      double dist=0;
      if (i>N*P) {

        // Find closest point on each of P paths, update state[P]
        for (int j=0; j<P; ++j) {
          dist=dist2((x+j*N)[state[j]], (x+j*N)[state[j]+1], ffv);
          double d=dist;
          int newstate=state[j];
          if (state[j]<K-2  // test if next state is closer
              && (d=dist2((x+j*N)[state[j]+1], (x+j*N)[state[j]+2], ffv))
              < dist) {
            newstate=state[j]+1;
            dist=d;
          }
          if ((R&1) && state[j]<K-3  // test state after next
              && (d=dist2((x+j*N)[state[j]+2], (x+j*N)[state[j]+3], ffv))
              < dist) {
            newstate=state[j]+2;
            dist=d;
          }
          if ((R&1) && state[j]>0  // test previous state
              && (d=dist2((x+j*N)[state[j]-1], (x+j*N)[state[j]], ffv))
                  < dist) {
            newstate=state[j]-1;
            dist=d;
          }
          if ((R&2) && K>2) {  // test random state
            int r=rand()%(K-1);
            if ((d=dist2((x+j*N)[r], (x+j*N)[r+1], ffv)) < dist) {
              dist=d;
              newstate=r;
            }
          }
          if (R&4) {  // test all states
            newstate=state[j];
            for (int k=0; k<K-1; ++k) {
              if ((d=dist2((x+j*N)[k], (x+j*N)[k+1], ffv)) < dist) {
                dist=d;
                newstate=k;
              }
            }
          }
          state[j]=newstate;
        }

        // Find box bounding the nearest points on all paths
        for (int j=0; j<P; ++j) {
          dist2((x+j*N)[state[j]], (x+j*N)[state[j]+1], ffv, p);
          for (int k=0; k<M; ++k) {
            if (j==0 || p[k]<lo[k]) lo[k]=p[k];
            if (j==0 || p[k]>hi[k]) hi[k]=p[k];
          }
        }

        // Set p to closest point to ffv in box
        copy(ffv, ffv+M, p);
        for (int j=0; j<M; ++j) {
          if (p[j]<lo[j]) p[j]=lo[j];
          if (p[j]>hi[j]) p[j]=hi[j];
        }

        // Anomaly score is distance squared from ffv to p
        dist=dist2(ffv, p);
        fa=fa*DT+dist*(1-DT);

        // Print test results
        fprintf(out, "%g", fa);
        for (int j=0; j<P; ++j)
          fprintf(out, " %d", state[j]);
        for (int j=0; j<M; ++j)
          fprintf(out, " %g", ffv[j]);
        fprintf(out, "\n");
        asum+=fa;
        if (fa>amax) amax=fa;
      }

      // Subsample
      for (int j=1; j<S; ++j)
        getline(f, s);
    }
    if (f) {
      fclose(f);
      printf("%-50s %12.6f %12.6f\n", filename.c_str(), amax, asum);
    }
  }
  printf("%1.2f sec. to test %d samples\n",
    elapsedTime(), i-N*P);
  fclose(out);
  printf("Anomaly scores saved to tsad4.out\n");
  return 0;
}