/* filtnn.c     Joe Mietus      Apr 1 2011 */

/* filtnn.c     Joe Mietus      Oct 7 2008 */

/* getline renamed to lineget May 18 2010 */

/*

filtnn :

Usage: filtnn [options] filt hwin

  Reads 2 (y,ann) or 3 (x,y,ann) columns from stdin and

  filters NN intervals by marking intervals not within `filt'

  range from the average of intervals

  within a window `hwin' distance on either side of

  current interval.

  NEED TO DO A `sort -n' ON FILTERED SERIES TO RESTORE TIME ORDER

  options are :

  [-n] : print ratio of nnout to nnin to rrin

  [-x min max] : exclude date outside min - max

  [-p] : print excluded data at start/end of hwin buffer

  [-C ann] : change annotation of excluded points to ann (default: X)

*/

#include <stdio.h>

#include <stdlib.h>

#define MAXSTR 256

#define MAXWIN 1024

char *prog, line[MAXSTR], tmp[MAXSTR], ann[MAXWIN+1];

int ncols;

static double t, x[MAXWIN+1], y[MAXWIN+1];

FILE *fp;

main(argc, argv)

int argc;

char * argv[];

{

    char c, lastann;

    int hwin, win, i, j, k, cflag, xflag, pflag, nflag, rrin, nnin, nnout;

    double filt, filtmin, filtmax, min, max, mtmp, sum, av;

    c = 'X';

    cflag = xflag = pflag = nflag = rrin = nnin = nnout = 0;

    sum = 0.0;

    prog = argv[0];

    if (argc < 3) {

        usage();

    exit(1);

    }

    i = 1;

    if ((filt = atof(argv[i])) <= 0) {

            fprintf(stderr, "%s: filt must be greater than 0\n",

                    prog);

            exit(2);

    }

    if ((hwin = atoi(argv[++i])) <= 0) {

            fprintf(stderr, "%s: hwin must be integer greater than 0\n",

                    prog);

            exit(2);

    }

    if ((win = 2*hwin) > MAXWIN) {

            fprintf(stderr, "%s: max hwin = %d\n", prog, MAXWIN/2);

            exit(2);

    }

    for ( ; ++i < argc && *argv[i] == '-'; ) {

        switch(argv[i][1]) {

            case 'n': nflag = 1;

                      break;

            case 'x': if (i+3 > argc) {

                      usage();

              exit(1);

              }

                  min = atof(argv[++i]);

                  max = atof(argv[++i]);

                      if (min > max) {

                  mtmp = min;

              min = max;

              max = mtmp;

                      }

                  xflag = 1;

                      break;

            case 'p': pflag = 1;

                      break;

            case 'C': cflag = 1;

                  c = *argv[++i];

                      break;

            default:  usage();

                      exit(1);

        }

    }

    filtmax = 1.0 + filt;

    filtmin = 1.0 - filt;

    i = 0;

    if (fgets(line, MAXSTR, stdin) == NULL) {

    fprintf(stderr, "%s : no data read\n", prog);

    exit(2);

    }

    if (sscanf(line, "%lf %lf %c %s", &x[i], &y[i], &ann[i], tmp) == 4) {

    fprintf(stderr, "%s : incorrectly formatted data\n", prog);

    exit(2);

    }

    else if (sscanf(line, "%lf %lf %c", &x[i], &y[i], &ann[i]) == 3) {

        ncols = 3;

    }

    else if (sscanf(line, "%lf %c %s", &y[i], &ann[i], tmp) == 3) {

    fprintf(stderr, "%s : incorrectly formatted data\n", prog);

    exit(2);

    }

    else if (sscanf(line, "%lf %c", &y[i], &ann[i]) == 2) {

        x[i] = t = y[i];

        ncols = 2;

    }

    else {

    fprintf(stderr, "%s : incorrectly formatted data\n", prog);

    exit(2);

    }

    rrin++;

    lastann = ann[i];

    if (pflag)

    cprintline(i, c);

    while (i<=win) {

        if (lineget(i) == EOF) {

        fprintf(stderr, "%s : insufficient data\n", prog);

            exit(2);

    }

        rrin++;

    if (ann[i] != 'N' || lastann != 'N') {

        if (pflag)

            printline(i);

            lastann = ann[i];

    }

        else if (xflag && (y[i]<min || y[i]>max)) {

        if (pflag)

            cprintline(i, c);

            lastann = ann[i];

        }

        else {

        if (pflag && i<hwin)

           cprintline(i, c);

        nnin++;

            lastann = ann[i];

        i++;

    }

    }

    j = hwin;

    for (i=0; i<=win; i++)

        sum += y[i];

    sum -= y[j];

    av = sum/win;

    sum += y[j] - y[0];

    if (y[j] <= filtmax*av && y[j] >= filtmin*av)

        nnout++;

    else

        ann[j] = c;

    printline(j);

    i = 0;

    while (lineget(i) != EOF) {

        rrin++;

    if (ann[i] != 'N' || lastann != 'N') {

        printline(i);

            lastann = ann[i];

    }

        else if (xflag && (y[i]<min || y[i]>max)) {

        cprintline(i, c);

            lastann = ann[i];

        }

        else {

        nnin++;

        if (++j > win)

            j = 0;

        sum += y[i];

        sum -= y[j];

        av = sum/win;

        lastann = ann[i];

        if (++i > win)

        i = 0;

        sum += y[j] - y[i];

        if (y[j] <= filtmax*av && y[j] >= filtmin*av) {

            printline(j);

            nnout++;

        }

        else

            cprintline(j, c);

    }

    }

    if (pflag)

        for (i=0; i<hwin; i++) {

            if (++j > win)

                j = 0;

            cprintline(j, c);

        }

    if (nflag)

        fprintf(stderr, "nnout : nnin : rrin = %d : %d : %d = %g : %g = %g\n", 

        nnout, nnin, rrin, (double)nnout/nnin, (double)nnin/rrin,

        (double)nnout/rrin);

    exit(0);

}

lineget(i)

int i;

{

    int j;

    static int n = 2;

    if (fgets(line, MAXSTR, stdin) == NULL)

        return(EOF);

    if (ncols == 2) {

        if ((j = sscanf(line, "%lf %c %s", &y[i], &ann[i], tmp)) != 2) {

            fprintf(stderr, "%s: incorrectly formatted data : ", prog);

        fprintf(stderr, "line = %d\n", n);

        exit(2);

        }

    t += y[i];

    x[i] = t;

    }

    else if (ncols == 3) {

        if ((j = sscanf(line, "%lf %lf %c %s", &x[i], &y[i], &ann[i], tmp)) != 3) {

            fprintf(stderr, "%s: incorrectly formatted data : ", prog);

        fprintf(stderr, "line = %d\n", n);

        exit(2);

        }

    }

    n++;

    return(j);

}

printline(i)

int i;

{

    if (ncols == 3)

        printf("%.9g ", x[i]);

    printf("%.8g %c\n", y[i], ann[i]);

}

cprintline(i, c)

int i;

char c;

{

    if (ncols == 3)

        printf("%.9g ", x[i]);

    printf("%.8g %c\n", y[i], c);

}

usage()

{

    fprintf(stderr, "Usage: %s [options] filt hwin\n", prog);

    fprintf(stderr, " Reads 2 (y,ann) or 3 (x,y,ann) columns from stdin and\n");

    fprintf(stderr, " filters NN intervals by marking intervals not within `filt'\n");

    fprintf(stderr, " range from the average of intervals\n");

    fprintf(stderr, " within a window `hwin' distance on either side of\n");

    fprintf(stderr, " current interval.\n");

    fprintf(stderr, " NEED TO DO A `sort -n' ON FILTERED SERIES");

    fprintf(stderr, " TO RESTORE TIME ORDER\n\n");

    fprintf(stderr, " options are :\n");

    fprintf(stderr, " [-n] : print ratio of nnout to nnin to rrin\n");

    fprintf(stderr, " [-x min max] : exclude date outside min - max\n");

    fprintf(stderr, " [-p] : print excluded data at start/end of hwin buffer\n");

    fprintf(stderr, " [-C ann] : change annotation of excluded points to ann");

    fprintf(stderr, " (default: X)\n");

}