//
// tistat - Time Interval Statistics
//
// 18-Dec-2004 /tvb
//

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define VERSION "1.3"

char Usage[] =
    "\n"
    "Usage: tistat [options] [output] < input\n"
    "\n"
    "Options:\n"
    "    /avg:#   average every N data points\n"
    "    /add:#   input offset\n"
    "    /mul:#   input scale factor\n"
    "    /nth:#   take only every Nth data point\n"
    "    /pct:#   outlier percent\n"
    "    /q       quiet (do not report outliers)\n"
    "    /quad    quadratic fit\n"
    "    /tau:#   set ADEV tau value\n"
    "    /units:s display units (e.g., us ns ps)\n"
    "    /v       verbose (show all statistics)\n"
    "    /?       /help\n"
    "\n"
    "Description:\n"
    "    Calculate basic statistics for phase data including mean,\n"
    "    standard deviation, slope, and correlation coefficient.\n"
    "    \n"
    "    Note that phase slope is frequency offset.\n"
    "    \n"
    "    The correlation coefficient for a good fit is close to 1.0.\n"
    "    When it is significantly less than that the fit is poor and\n"
    "    the frequency offset is not certain and should be ignored\n"
    "    until more data can gathered to establish a clear trend.\n"
    "    \n"
    "    Allan deviation calculations are based on 1 Hz data.\n"
    "    \n"
    "    Can be used directly with raw data from HP 53131A, 53132A,\n"
    "    SR 620, TSC 5110A, and other TI (time interval) counters.\n"
    "    \n"
    "    Phase residuals are written to the optional output file.\n"
    "    \n"
    "Examples:\n"
    "    tistat < log16626.txt\n"
    "    tistat /mul:1-3 /units:ns < log16626.txt\n"
    "    tistat /mul:1-3 /units:ns resid.txt < log16626.txt\n"
    "    \n"
    "Version: %s\n"
    ;

#include "./stat_lib.h"

//
// Define program options.
//

struct {
    double Add;
    long Avg;
    double Mul;
    long Nth;
    double Pct;
    int Quiet;
    int Quad;
    long Tau;
    char *Units;
    int Verbose;
} Options;

int GetOptions (int argc, char *argv[]);
void read_data (PSTAT s);
int outlier (double x, double err);
int valid_data (char *s);
void display_stats(PSTAT s);

double Tau;

int main (int argc, char *argv[])
{
    int args;
    char *resid_path;
    STAT stats, resid;

    // Set option defaults

    Options.Units = "unit?";
    Options.Add = 0.0;
    Options.Avg = 1;
    Options.Mul = 1.0;
    Options.Nth = 1;
    Options.Pct = 5.0;
    Options.Tau = 1;

    // Get options

    args = GetOptions(argc, argv);
    argc -= args;
    argv += args;

    // Get optional residuals output file

    resid_path = NULL;
    if (argc > 1) {
        resid_path = argv[1];
        argc -= 1;
        argv += 1;
    }
    if (argc > 1) {
        fprintf(stderr, "Unexpected argument: %s\n", argv[1]);
        exit(1);
    }

    // Read raw data

    stat_init(&stats);

    stats.offset = Options.Add;
    stats.scale = Options.Mul;
    Tau = Options.Nth * Options.Avg * Options.Tau;
    stats.units = Options.Units;

    read_data(&stats);

    // Compute and display statistics

    stat_calc(&stats);
    if (Options.Verbose) {
        stat_dump(&stats);
    }

    display_stats(&stats);

    //
    // Create residual data set from calculated slope (frequency offset).
    // - standard deviation of residuals is of interest
    // - mean, slope, and r^2 of residuals should be close to zero
    // - adev will be the same (it is invariant to frequency offset)
    //

    if (resid_path) {
        FILE *f;
        long i;

        f = fopen(resid_path, "w");
        if (f == NULL) {
            fprintf(stderr, "Open failed: %s\n", argv[1]);
            exit(1);
        }

        stat_line_fit(&stats, &resid);
        stat_calc(&resid);
        if (Options.Verbose > 1) {
            stat_dump(&resid);
        }
        for (i = 0; i < (long)resid.count; i += 1) {
            fprintf(f, "%.8le\n", resid.data[i]);
        }
        stat_free(&resid);
        fclose(f);
    }

    //
    // Perform quadratic fit and display frequency offset (phase slope)
    // and frequency drift rate.
    //

    if (Options.Quad) {
        stat_quad_fit(&stats, &resid);
        printf("%16.4le frequency offset\n", stats.slope2 / Tau);
        printf("%16.4le frequency drift / day\n", stats.drift2 * 86400 / Tau / Tau);
        printf("%16.4le sdev\n", stats.sdev2);

        if (Options.Verbose) {
            stat_dump(&resid);
        }
        stat_free(&resid);
    }

    stat_free(&stats);

    return 0;
}

//
// Parse options.
//

int GetOptions (int argc, char *argv[])
{
    int args = 0;

    while (argc > 1 && argv[1][0] == '/') {

        if (strcmp(argv[1], "/v") == 0) {
            Options.Verbose += 1;
        } else

        if (strcmp(argv[1], "/q") == 0) {
            Options.Quiet = 1;
        } else

        if (strcmp(argv[1], "/quad") == 0) {
            Options.Quad = 1;
        } else

        if (strncmp(argv[1], "/add:", 5) == 0) {
            sscanf(&argv[1][5], "%lf", &Options.Add);
        } else

        if (strncmp(argv[1], "/avg:", 5) == 0) {
            sscanf(&argv[1][5], "%ld", &Options.Avg);
        } else

        if (strncmp(argv[1], "/nth:", 5) == 0) {
            sscanf(&argv[1][5], "%ld", &Options.Nth);
        } else

        if (strncmp(argv[1], "/mul:", 5) == 0) {
            sscanf(&argv[1][5], "%lf", &Options.Mul);
        } else

        if (strncmp(argv[1], "/pct:", 5) == 0) {
            sscanf(&argv[1][5], "%lf", &Options.Pct);
        } else

        if (strncmp(argv[1], "/tau:", 5) == 0) {
            sscanf(&argv[1][5], "%ld", &Options.Tau);
        } else

        if (strncmp(argv[1], "/units:", 7) == 0) {
            Options.Units = &argv[1][7];
        } else

        if (strcmp(argv[1], "/?") == 0 || strcmp(argv[1], "/help") == 0) {
            fprintf(stderr, Usage, VERSION);
            exit(0);

        } else {
            fprintf(stderr, "Unknown option: %s\n", argv[1]);
            exit(1);
        }
        args += 1;
        argc -= 1;
        argv += 1;
    }

    return args;
}

//
// Read input file.
// - check for obviously bad data
// - optionally take only one of every N data points
// - optionally average M points together into one mean
// - or do both (N * M)
//

void read_data (PSTAT s)
{
    long avg, nth;
    char line[1000];
    double phase, sum;

    avg = nth = 0;
    sum = 0.0;
    while (fgets(line, sizeof line, stdin) != NULL) {
        if (valid_data(line)) {
            if (sscanf(line, "%lf", &phase) == 1) {
                if ( ! outlier(phase, Options.Pct) ) {
                    nth += 1;
                    if ((nth % Options.Nth) == 0) {
                        sum += phase;
                        avg += 1;
                        if ((avg % Options.Avg) == 0) {
                            stat_data(s, sum / Options.Avg);
                            sum = 0.0;
                        }
                    }
                }
            }
        }
    }
}

//
// Process HP 53131A data format.
//

int valid_data (char *s)
{
    char *p;

    // Remove user comments
    if (s[0] == '#') {
        return 0;
    }

    // Remove HP 53131A statistics
    if (strchr(s, ':') != 0) {
        // printf("Comment: %s", s);
        return 0;
    }

    // Remove embedded commas (HP 53131A)
    while ((p = strchr(s, ',')) != NULL) {
        strcpy(p, p + 1);
    }

    // Remove HP 53131A units
    if ((p = strstr(s, "us")) != NULL) {
        *p = '\0';
    }

    return 1;
}

//
// Simple outlier dectection.
//

#define AVG_WIN 10

int outlier (double v, double pct)
{
    static long count = 0;
    static double hist[AVG_WIN];
    static long line = 0;

    double avg, sum;
    int i;

    if (pct == 0.0) {
        return 0;
    }
    line += 1;
    if (count >= AVG_WIN) {
        sum = 0;
        for (i = 0; i < AVG_WIN; i += 1) {
            sum += hist[i];
        }
        avg = sum / AVG_WIN;
        if (fabs(1.0 - v / avg) >= (pct / 100.0)) {
            if (!Options.Quiet) {
                fprintf(stderr, "** Line %ld: ", line);
                fprintf(stderr, "outlier (%g) exceeds %g%% of running mean (%g)\n",
                        v, pct, avg);
            }
            return 1;
        }
    }

    hist[count % AVG_WIN] = v;
    count += 1;

    return 0;
}

//
// Print relevant statistics
//

#define WARNING1 "  <-- WARNING: correlation too low"
#define WARNING2 "  <-- WARNING: S/N too low"
#define SUGGEST "\n** Check for glitches or gather a lot more data and retry **"

void display_stats(PSTAT s)
{
    char *warn1, *warn2;

    printf("%16.0lf count\n", s->count);

    if (s->count > 1) {
        printf("%16.8le %s mean\n", s->mean, s->units);
        printf("%16.8le %s range\n", s->range, s->units);
    }

    if (s->count > 2) {
        printf("%16.8le %s sdev\n", s->sdev, s->units);
        printf("%16.8le %s sdev_resid\n", s->steyx, s->units);
    }

    if (s->count > 3) {
        printf("%16.8le %s adev\n", s->adev / Tau, s->units);
        printf("%16.4le %s/s slope\n", s->slope / Tau, s->units);
        printf("%16.4le %s/s slope_sdev\n", s->slope_sdev, s->units);

        warn1 = (s->r < 0.5) ? WARNING1 : "";
        printf("%16.4lf r^2%s\n", s->r * s->r, warn1);

        warn2 = (s->sn < 2.0) ? WARNING2 : "";
        printf("%16.1lf slope:noise ratio%s\n", s->sn, warn2);

        if (*warn1 || *warn2) {
            printf("%s\n", SUGGEST);
        }
    }
}

// poor man's makefile...

#include "./stat_lib.c"