//
// da -- command line data analysis tool
//
// This tool operates like one-column command-line "Excel" or a "sed" for numerical data.
// It is implemented as a typical unix filter; very convenient for working with large files
// of clock phase & frequency data. I've been using it for years and add options as needed.
//
// The general flow is raw data is read, partitioning or decimation operations performed,
// a series of one or more math or filtering operations are applied, and resulting data is
// then written.
//
// The goal is documented, fully repeatable, and scriptable operations on raw data instead
// of ad-hoc mouse clicks of GUI-based clock analysis tools.
//
// 30-Dec-2010 Tom Van Baak (tvb) www.LeapSecond.com/tools
//

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

char *Tool;
char Ver[] = "2.04";
char Usage[] = "Usage: %s [ops] <input >output\n"
    "\nInput ops:\n"
    "    in f      yes $     no $      begin $   after $   before $  mth k     nth k     col n     fld n     fld n.i\n"
    "    head k    tail k    skip k    take k    trim k\n"
    "    oct       hex       hp        hms       ddmm      comma\n"
    "\nOutput ops:\n"
    "    out f     fmt $     e #.#     f #.#     g #.#     plus\n"
    "\nPrune ops:\n"
    "    lowbar #  highbar # above #   below #   nonzero   wrap #    unwrap #\n"
    "    within #  without # at i      from i    to i      count i   offset #  slope #   drift #   fit k\n"
    "    del i     ins i #   app i #   repl i #  inside # #          outside # #         delete i j\n"
    "\nMath ops:\n"
    "    neg       inv       floor     ceil      int       trunc     frac      abs       sqrt      ln        sin       cos\n"
    "    add #     sub #     rsub #    mul #     div #     rdiv #    mod #     pow #     exp #     log #\n"
    "    +#        -#        *#        /#\n"
    "    rel #     zero      one       first0    first1    min0\n"
    "    sum k     avg k     runavg k  unavg k\n"
    "\nTime (phase), aboslute frequency, relative frequency error:\n"
    "    t2e tau   e2t tau   f2e tau   e2f tau\n"
    "\nOther ops:\n"
    "    sort      uniq      rev       rat2      diff2     delta     rsum\n"
    "    total     mean      rms       sdev      wc        stats     pop\n"
    "    ver       nop       v\n"
    "\nExamples:\n"
    "    da <z.dat fld 2 /2e5 skip 1700 take 20000 fit 1 > clock.dat\n"
    "    da in z.dat fld 2 /200000 from 7480 offset 7.04024e-010 skip 1700 out z3.dat\n"
    ;

// Program global options.
struct {
    char *in, *out, format[100], *yes, *no, *begin, *after, *before;
    long skip, take, mth, nth, col, fld, sub, from, to, count;
    int verbose, comma, hp, oct, hex, hms, ddmm, pop, plus;
} opt;

int command3 (char *cmd, char *arg1, char *arg2);
int command2 (char *cmd, char *arg);
int command1 (char *cmd);

// Display error message and quit.
void fatal_exit (char *format, ...)
{
    char buf[1000];
    va_list vap;
    va_start(vap, format);
    vsprintf(buf, format, vap);
    fprintf(stderr, "%s: ERROR: %s.\n", Tool, buf);
    va_end(vap);
    exit(1);
}

int main (int argc, char *argv[])
{
    Tool = argv[0];
    if (argc == 1 || strstr(argv[1], "?") != NULL) {
        fprintf(stderr, Usage, Tool);
        return 1;
    }

    // Set defaults.
    command2("in", "-");
    command2("fmt", "%.15lg");

    // Process options left to right.
    while (argc > 1) {
        if (argc > 3 && command3(argv[1], argv[2], argv[3])) {
            argc -= 3; argv += 3;
            continue;
        }
        if (argc > 2 && command2(argv[1], argv[2])) {
            argc -= 2; argv += 2;
            continue;
        }
        if (argc > 1 && command1(argv[1])) {
            argc -= 1; argv += 1;
            continue;
        }
        fatal_exit("operation `%s' unknown or missing parameter", argv[1]);
    }

    // Write output if not done already.
    if (opt.out == NULL) {
        command2("out", "-");
    }
    return 0;
}

double *get_data (void);
long get_count (void);
double *new_data (long count);
void set_data (double *newdata, long count);
void write_output (void);
int get_from_to (long *from, long *to);

void fit0 (double *data, long n, double c[1]);
void fit1 (double *data, long n, double c[2]);
void fit2 (double *data, long n, double c[3]);
double mean_sdev (double *data, long n, double *sdev, int pop);

#define CMD(s) (strcmp(cmd, s) == 0)
#define MINUS1(i) ((i) - 1)    // since user row 1 is array index 0
#define PLUS1(i)  ((i) + 1)    // since array index 0 is user row 1
#define FPRINTF opt.verbose && fprintf

// The input file is read at the last possible moment. Operations are performed left
// to right. So some commands are too late if input data has already been read.
void too_late1 (char *cmd)
{
    if (get_count()) {
        fatal_exit("operation `%s' too late (data already read)", cmd);
    }
}
void too_late2 (char *cmd, char *arg)
{
    if (get_count()) {
        fatal_exit("operation `%s %s' too late (data already read)", cmd, arg);
    }
}

// Handle commands with two arguments.

int command3 (char *cmd, char *arg1, char *arg2)
{
    // Keep only points inside or outside bounds.
    if (CMD("inside") || CMD("outside")) {
        double lo_val = atof(arg1);
        double hi_val = atof(arg2);
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        int inside  = CMD("inside");
        int outside = CMD("outside");

        for (i = 0; i < n; i += 1) {
            int within = (d[i] >= lo_val && d[i] <= hi_val);
            if (inside && within || outside && !within) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Delete block of points.
    if (CMD("delete")) {
        long from = MINUS1(atol(arg1));
        long to   = MINUS1(atol(arg2));
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);

        if (from < 0 || from >= n || to < 0 || to >= n || from > to) {
            fatal_exit("operation `%s %s %s' invalid", cmd, arg1, arg2);
        }
        for (i = 0; i < n; i += 1) {
            if (i < from || i > to) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Insert or append one data point.
    if (CMD("ins") || CMD("app")) {
        long here = MINUS1(atol(arg1));
        double datum = atof(arg2);
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n + 1);
        int ins = CMD("ins");
        int app = CMD("app");

        if (here < 0 || here >= n) {
            fatal_exit("operation `%s %s %s' invalid", cmd, arg1, arg2);
        }
        for (i = 0; i < n; i += 1) {
            if (ins && i == here) {
                newdata[k++] = datum;
            }
            newdata[k++] = d[i];
            if (app && i == here) {
                newdata[k++] = datum;
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Replace one data point.
    if (CMD("repl")) {
        long here = MINUS1(atol(arg1));
        double datum = atof(arg2);
        double *d = get_data();
        long n = get_count();

        if (here < 0 || here >= n) {
            fatal_exit("operation `%s %s %s' invalid", cmd, arg1, arg2);
        }
        FPRINTF(stderr, "%s: line %ld: %.15le -> %.15le\n", Tool, PLUS1(here), d[here], datum);
        d[here] = datum;
        return 1;
    }

    return 0;
}

// Handle commands with one argument.

int command2 (char *cmd, char *arg)
{
    double xval = atof(arg);
    long   ival = atol(arg);

    if (CMD("in")) {
        too_late2(cmd, arg);
        opt.in = arg;
        return 1;
    }

    if (CMD("out")) {
        opt.out = arg;
        get_data();
        write_output();
        return 1;
    }

    if (CMD("yes")) {
        too_late2(cmd, arg);
        opt.yes = arg;
        return 1;
    }

    if (CMD("no")) {
        too_late2(cmd, arg);
        opt.no = arg;
        return 1;
    }

    if (CMD("begin")) {
        too_late2(cmd, arg);
        opt.begin = arg;
        return 1;
    }

    if (CMD("after")) {
        too_late2(cmd, arg);
        opt.after = arg;
        return 1;
    }

    if (CMD("before")) {
        too_late2(cmd, arg);
        opt.before = arg;
        return 1;
    }

    if (CMD("fmt")) {
        // Leave room for \n and \0.
        if (strlen(arg) > sizeof opt.format - 2) {
            fatal_exit("operation `%s %s' too long", cmd, arg);
        }
        sprintf(opt.format, "%s\n", arg);
        return 1;
    }

    if (CMD("e") || CMD("f") || CMD("g")) {
        int n, f;
        if (sscanf(arg, "%d.%d\n", &n, &f) == 2) {
            sprintf(opt.format, "%%%d.%dl%s\n", n, f, cmd);
        } else
        if (sscanf(arg, "%d", &n) == 1) {
            sprintf(opt.format, "%%%dl%s\n", n, cmd);
        } else
        if (sscanf(arg, ".%d", &f) == 1) {
            sprintf(opt.format, "%%.%dl%s\n", f, cmd);
        } else
        {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        return 1;
    }

    if (CMD("mth")) {
        too_late2(cmd, arg);
        opt.mth = ival;
        return 1;
    }

    if (CMD("nth")) {
        too_late2(cmd, arg);
        opt.nth = ival;
        return 1;
    }

    if (CMD("col")) {
        if (ival < 1) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        too_late2(cmd, arg);
        opt.col = ival;
        return 1;
    }

    if (CMD("fld")) {
        if (ival < 1) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        too_late2(cmd, arg);
        if (sscanf(arg, "%ld.%ld", &opt.fld, &opt.sub) != 2) {
            opt.fld = ival;
            opt.sub = 0;
        }
        return 1;
    }

    if (CMD("add")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] += xval;
        }
        return 1;
    }

    if (CMD("sub")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] -= xval;
        }
        return 1;
    }

    if (CMD("rsub") || CMD("rub")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = xval - d[i];
        }
        return 1;
    }

    if (CMD("mul")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] *= xval;
        }
        return 1;
    }

    if (CMD("div")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] /= xval;
        }
        return 1;
    }

    if (CMD("rdiv") || CMD("riv")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = xval / d[i];
        }
        return 1;
    }

    if (CMD("mod")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = fmod(d[i], xval);
        }
        return 1;
    }

    if (CMD("pow")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = pow(d[i], xval);
        }
        return 1;
    }

    if (CMD("exp")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = pow(xval, d[i]);
        }
        return 1;
    }

    if (CMD("log")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
            d[i] = (d[i] == 0.0) ? 1e-99 : log(d[i]) / log(xval);
        }
        return 1;
    }

    // Replace points below X with X.
    if (CMD("lowbar")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           if (d[i] < xval) d[i] = xval;
        }
        return 1;
    }

    // Replace points above X with X.
    if (CMD("highbar")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           if (d[i] > xval) d[i] = xval;
        }
        return 1;
    }

    // Keep only points above X.
    if (CMD("above")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (d[i] > xval) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Keep only points below X.
    if (CMD("below")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (d[i] < xval) {
                newdata[k++] = d[i];
            } else {
                FPRINTF(stderr, "%s: line %ld: %.15le\n", Tool, PLUS1(i), d[i]);
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Delete any zero points.
    if (CMD("nozero")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (d[i] != 0.0) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Handle phase wrap.
    if (CMD("wrap")) {
        double *d = get_data();
        long i, n = get_count();
        long cycles = 0;
        for (i = 0; i < n; i += 1) {
            double diff = d[i] - cycles * xval;
            if (diff >=  xval) cycles += 1;
            if (diff <= -xval) cycles -= 1;
            d[i] -= cycles * xval;
        }
        return 1;
    }

    // Handle phase unwrap.
    if (CMD("unwrap")) {
        double *d = get_data();
        long i, n = get_count();
        long cycles = 0;
        double d_prev = d[0];

        // For example, allow unwrap 2^16 or 2p32.
        if (strlen(arg) >= 3 && arg[0] == '2' && (arg[1] == '^' || arg[1] == 'p')) {
            double exponent = atoi(&arg[2]);
            xval = pow(2.0, exponent);
        }

        for (i = 0; i < n; i += 1) {
            double diff = d[i] - d_prev;
            if (diff <= -xval / 3.0) cycles += 1;
            if (diff >=  xval / 3.0) cycles -= 1;
            d_prev = d[i];
            d[i] += cycles * xval;
        }
        return 1;
    }

    if (CMD("rel")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
            d[i] = (d[i] - xval) / xval;
        }
        return 1;
    }

    // If data has not yet been read do skip/take efficiently as file is being read.
    if (get_count() == 0) {
        if (CMD("skip")) {
            opt.skip += ival;
            FPRINTF(stderr, "%s: operation `%s %s' pending.\n", Tool, cmd, arg);
            return 1;
        }
        if (CMD("head") || CMD("take")) {
            opt.take += ival;
            FPRINTF(stderr, "%s: operation `%s %s' pending.\n", Tool, cmd, arg);
            return 1;
        }
    }

    // Relationship among head tail skip trim operations:
    // ---------------------------------------- original data (40 points)
    // ----------                               head 10
    //           ------------------------------ skip 10
    //                               ---------- tail 10
    // ------------------------------           trim 10
    //           ----------                     skip 10 head 10
    //           --------------------           skip 10 trim 10
    //                     ----------           trim 10 tail 10

    if (CMD("head") || CMD("take") || CMD("tail") || CMD("skip") || CMD("trim")) {
        double *newdata, *d = get_data();
        long i, n = get_count();
        long offset = 0, length = n;

        if (CMD("head") || CMD("take")) {
            offset = 0; length = ival;
        }
        if (CMD("skip")) {
            offset = ival; length = n - ival;
        }
        if (CMD("tail")) {
            offset = n - ival; length = ival;
        }
        if (CMD("trim")) {
            offset = 0; length = n - ival;
        }
        if (offset < 0 || offset > n || length < 0 || length > n) {
            fatal_exit("operation `%s %s' invalid count (%ld...%ld)", cmd, arg, 0, n);
        }
        newdata = new_data(length);
        for (i = 0; i < length; i += 1) {
            newdata[i] = d[i + offset];
        }
        set_data(newdata, length);
        return 1;
    }

    // Keep points near to zero.
    if (CMD("within")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (fabs(d[i]) <= xval) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Keep points away from zero.
    if (CMD("without")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (fabs(d[i]) >= xval) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Delete one data point.
    if (CMD("del")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n - 1);
        long here = MINUS1(ival);
        if (here < 0 || here >= n) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        for (i = 0; i < n; i += 1) {
            if (i != here) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // Commands at/from/to/count allow user to specify point or range of points.
    if (CMD("at")) {
        if (ival < 1) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        opt.from = MINUS1(ival);
        opt.count = 1;
        opt.to = 0;
        return 1;
    }

    if (CMD("from")) {
        if (ival < 1) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        opt.from = MINUS1(ival);
        return 1;
    }

    if (CMD("to")) {
        if (ival < 1) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        opt.to = MINUS1(ival);
        opt.count = 0;
        return 1;
    }

    if (CMD("count")) {
        if (ival < 1) {
            fatal_exit("operation `%s %s' invalid", cmd, arg);
        }
        opt.count = ival;
        opt.to = 0;
        return 1;
    }

    // Apply offset, slope, or drift correction to all or partial data set.
    if (CMD("offset") || CMD("slope") || CMD("drift")) {
        double *d = get_data();
        long from, to, i, n = get_count();
        double a0 = 0, a1 = 0, a2 = 0;

        // Get optional partial range of fit.
        if ( !get_from_to(&from, &to) ) {
            fatal_exit("operation `%s %s' invalid from/to/count", cmd, arg);
        }
        if (CMD("offset")) a0 = xval;
        if (CMD("slope"))  a1 = xval;
        if (CMD("drift"))  a2 = xval;

        for (i = from; i < to; i += 1) {
            d[i] += a0;
            a0 += a1;
            a1 += a2;
        }
        return 1;
    }

    if (CMD("fit")) {
        int order = ival;
        double coeff[3], *d = get_data();
        long from, to, x, n = get_count();
        coeff[2] = coeff[1] = coeff[0] = 0.0;

        // Get optional partial range of fit.
        if ( !get_from_to(&from, &to) ) {
            fatal_exit("operation `%s %s' invalid from/to/count", cmd, arg);
        }
        // Fit polynomial based on entire or partial data set.
        switch (order) {
            case  0 : fit0(&d[from], to - from, coeff); break;
            case  1 : fit1(&d[from], to - from, coeff); break;
            case  2 : fit2(&d[from], to - from, coeff); break;
        default :
            fatal_exit("%s `%s' not supported", cmd, arg);
        }

        // Replace data set with residuals.
        for (x = 0; x < n; x += 1) {
            double y = (coeff[2] * x * x) + (coeff[1] * x) + coeff[0];
            d[x] -= y;
        }

        // Echo fit parameters.
        switch (order) {
            case  2 : FPRINTF(stderr, "%s: %s: %14.6lg drift (x2=%lg/d)\n", Tool, cmd,
                              coeff[2], 2.0 * coeff[2] * 86400.0);
            case  1 : FPRINTF(stderr, "%s: %s: %14.6lg slope  (%.6le)\n", Tool, cmd, coeff[1], coeff[1]);
            case  0 : FPRINTF(stderr, "%s: %s: %14.6lg offset (%.6le)\n", Tool, cmd, coeff[0], coeff[0]);
        }
        return 1;
    }

    if (CMD("avg") || CMD("sum")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double sum = 0, *newdata = new_data(n);
        double divisor = CMD("avg") ? ival : 1;

        if (divisor < 1 || ival != xval) {
            fatal_exit("`%s %s' invalid", cmd, arg);
        }
        for (i = 0; i < n; i += 1) {
            sum += d[i];
            if ((i % ival) == ival - 1) {
                newdata[k++] = sum / divisor;
                sum = 0;
            }
        }
        set_data(newdata, k);
        return 1;
    }

    // RUNAVG n -- forward running average.
    if (CMD("runavg")) {
        long weight = ival;
        double *d = get_data();
        long i, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (i == 0) {
                newdata[i] = d[i];
            } else {
                newdata[i] = (newdata[i - 1] * (weight - 1) + d[i]) / weight;
            }
        }
        set_data(newdata, n);
        return 1;
    }

    // UNAVG n -- undo running average.
    if (CMD("unavg")) {
        long weight = ival;
        double *d = get_data();
        long i, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (i == 0) {
                newdata[i] = d[i];
            } else {
                newdata[i] = (weight * d[i]) - ((weight - 1) * d[i - 1]);
            }
        }
        set_data(newdata, n);
        return 1;
    }

    // T2E tau -- convert time (accumulated phase) to frequency error.
    if (CMD("t2e")) {
        double *d = get_data();
        long i, n = get_count();
        double *newdata = new_data(n - 1);
        for (i = 0; (i + 1) < n; i += 1) {
            newdata[i] = (d[i + 1] - d[i]) / xval;
        }
        set_data(newdata, n - 1);
        return 1;
    }

    // E2T tau -- convert frequency error to time (accumulated phase).
    if (CMD("e2t")) {
        double *d = get_data();
        long i, n = get_count();
        double *newdata = new_data(n + 1);
        double ti = 0.0;
        newdata[0] = ti;
        for (i = 0; i < n; i += 1) {
            ti += d[i] * xval;
            newdata[i + 1] = ti;
        }
        set_data(newdata, n + 1);
        return 1;
    }

    // F2E f0 -- convert absolute frequency to relative frequency error.
    if (CMD("f2e")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
            d[i] = (d[i] - xval) / xval;
        }
        return 1;
    }

    // E2F f0 -- convert relative frequency error to absolute frequency.
    if (CMD("e2f")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
            d[i] = xval * (1 + d[i]);
        }
        return 1;
    }

    return 0;
}

// For qsort() function below.
int double_cmp (const void *ap, const void *bp)
{
    double a = *(double *)ap;
    double b = *(double *)bp;
    return (a < b) ? -1 : (a > b) ? 1 : 0;
}

// Handle commands with no argument.

int command1 (char *cmd)
{
    if (CMD("neg")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = -d[i];
        }
        return 1;
    }

    if (CMD("inv")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = 1 / d[i];
        }
        return 1;
    }

    if (CMD("floor")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = floor(d[i]);
        }
        return 1;
    }

    if (CMD("ceil")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = ceil(d[i]);
        }
        return 1;
    }

    // Convert to integer, truncate.
    if (CMD("trunc")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = (long) d[i];
        }
        return 1;
    }

    // Convert to integer, rounding up.
    if (CMD("int")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = (long) (d[i] + 0.5);
        }
        return 1;
    }

    // Retain fraction only.
    if (CMD("frac")) {
        double ignore, *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = modf(d[i], &ignore);
        }
        return 1;
    }

    if (CMD("abs")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = fabs(d[i]);
        }
        return 1;
    }

    if (CMD("sqrt")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = sqrt(d[i]);
        }
        return 1;
    }

    if (CMD("ln")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
            d[i] = (d[i] == 0.0) ? 1e-99 : log(d[i]);
        }
        return 1;
    }

    if (CMD("sin")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = sin(d[i]);
        }
        return 1;
    }

    if (CMD("cos")) {
        double *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
           d[i] = cos(d[i]);
        }
        return 1;
    }

    if (CMD("nop")) {
        get_data();
        return 1;
    }

    if (CMD("ver")) {
        fprintf(stderr, "%s: version %s (%s %s) %s\n", Tool, Ver, __DATE__, __TIME__, "www.LeapSecond.com/tools");
        exit(0);
    }

    if (CMD("pop")) {
        opt.pop = 1;
        return 1;
    }

    if (CMD("plus")) {
        opt.plus = 1;
        return 1;
    }

    if (CMD("v")) {
        opt.verbose ^= 1;
        return 1;
    }

    if (CMD("comma")) {
        too_late1(cmd);
        opt.comma = 1;
        return 1;
    }

    if (CMD("hp")) {
        too_late1(cmd);
        opt.hp = 1;
        return 1;
    }

    if (CMD("oct")) {
        too_late1(cmd);
        opt.oct = 1;
        return 1;
    }

    if (CMD("hex")) {
        too_late1(cmd);
        opt.hex = 1;
        return 1;
    }

    if (CMD("hms")) {
        too_late1(cmd);
        opt.hms = 1;
        return 1;
    }

    if (CMD("ddmm")) {
        too_late1(cmd);
        opt.ddmm = 1;
        return 1;
    }

    if (CMD("wc")) {
        double *newdata = new_data(1);
        get_data();
        newdata[0] = get_count();
        set_data(newdata, 1);
        return 1;
    }

    if (CMD("total") || CMD("mean")) {
        double *d = get_data();
        long n = get_count();
        double mean, *newdata = new_data(1);
        mean = mean_sdev(d, n, NULL, 0);
        if (CMD("mean"))  newdata[0] = mean;
        if (CMD("total")) newdata[0] = mean * n;
        set_data(newdata, 1);
        return 1;
    }

    if (CMD("sdev") || CMD("rms")) {
        double *d = get_data();
        long n = get_count();
        double sdev, *newdata = new_data(1);
        int pop = CMD("rms") ? 0 : opt.pop;
        mean_sdev(d, n, &sdev, pop);
        newdata[0] = sdev;
        set_data(newdata, 1);
        return 1;
    }

    if (CMD("stats")) {
        double *d = get_data();
        long n = get_count();
        double mean, sdev, *newdata;
        mean = mean_sdev(d, n, &sdev, opt.pop);
        newdata = new_data(3);
        newdata[0] = n;
        newdata[1] = mean;
        newdata[2] = sdev;
        set_data(newdata, 3);
        return 1;
    }

    if (CMD("zero") || CMD("norm")) {
        double mean, *d = get_data();
        long i, n = get_count();
        mean = mean_sdev(d, n, NULL, 0);
        for (i = 0; i < n; i += 1) {
            d[i] -= mean;
        }
        return 1;
    }

    if (CMD("one")) {
        double mean, sdev, *d = get_data();
        long i, n = get_count();
        mean = mean_sdev(d, n, &sdev, opt.pop);
        for (i = 0; i < n; i += 1) {
            d[i] = (d[i] - mean) / sdev;
        }
        return 1;
    }

    if (CMD("sort")) {
        double *d = get_data();
        long n = get_count();
        qsort(d, n, sizeof d[0], double_cmp);
        return 1;
    }

    if (CMD("uniq")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n);
        for (i = 0; i < n; i += 1) {
            if (i == 0 || d[i] != d[i - 1]) {
                newdata[k++] = d[i];
            }
        }
        set_data(newdata, k);
        return 1;
    }

    if (CMD("rev")) {
        double temp, *d = get_data();
        long i, j, n = get_count();
        for (i = 0, j = n - 1; i < j; i += 1, j -= 1) {
            temp = d[i];
            d[i] = d[j];
            d[j] = temp;
        }
        return 1;
    }

    // Pair-wise ratio.
    if (CMD("rat2")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n / 2);
        for (i = 0; (i + 1) < n; i += 2) {
            newdata[k++] = d[i] / d[i + 1];
        }
        set_data(newdata, k);
        return 1;
    }

    // Pair-wise difference.
    if (CMD("diff2")) {
        double *d = get_data();
        long i, k = 0, n = get_count();
        double *newdata = new_data(n / 2);
        for (i = 0; (i + 1) < n; i += 2) {
            newdata[k++] = d[i] - d[i + 1];
        }
        set_data(newdata, k);
        return 1;
    }

    // Delta (differentiation).
    if (CMD("delta")) {
        double *d = get_data();
        long i, n = get_count();
        double *newdata = new_data(n - 1);
        for (i = 0; (i + 1) < n; i += 1) {
            newdata[i] = d[i + 1] - d[i];
        }
        set_data(newdata, n - 1);
        return 1;
    }

    // Running sum (integration).
    if (CMD("rsum")) {
        double sum = 0, *d = get_data();
        long i, n = get_count();
        for (i = 0; i < n; i += 1) {
            sum += d[i];
            d[i] = sum;
        }
        return 1;
    }

    // Apply offset so lowest point is 0.
    if (CMD("min0")) {
        double *d = get_data();
        long i, n = get_count();
        double lowest = d[0];
        for (i = 0; i < n; i += 1) {
            if (d[i] < lowest) lowest = d[i];
        }
        for (i = 0; i < n; i += 1) {
            d[i] -= lowest;
        }
        return 1;
    }

    // Scale (subtract) so first is 0.
    if (CMD("first0")) {
        double *d = get_data();
        long i, n = get_count();
        double first = d[0];
        for (i = 0; i < n; i += 1) {
           d[i] -= first;
        }
        return 1;
    }

    // Scale (divide) so first is 1.
    if (CMD("first1")) {
        double *d = get_data();
        long i, n = get_count();
        double first = d[0];
        for (i = 0; i < n; i += 1) {
           d[i] /= first;
        }
        return 1;
    }

    // Allow +N -N *N /N short-cuts for add N sub N mul N div N

    if (cmd[0] == '+' && isdigit(cmd[1])) return command2("add", cmd + 1);
    if (cmd[0] == '-' && isdigit(cmd[1])) return command2("sub", cmd + 1);
    if (cmd[0] == '*' && isdigit(cmd[1])) return command2("mul", cmd + 1);
    if (cmd[0] == '/' && isdigit(cmd[1])) return command2("div", cmd + 1);

    return 0;
}

struct {
    double *array;
    long count;
    long max;
} data;

// String match with regex-like hack.
int match (char *str, char *patt)
{
    // match empty string (blank line)
    if (strcmp(patt, "^$") == 0) {
        return (str[0] == '\0');
    }
    // match exact string
    if (strcmp(patt, "^") == 0) {
        return (strncmp(str, patt + 1, strlen(str) - 1) == 0);
    }
    // match string within
    return (strstr(str, patt) != NULL);
}

// Read input file.
void read_data (FILE *file)
{
    char line[1000];
    long lineno = 0, count = 0;
    size_t size;

    // Allocate small initial array.
    data.max = 128 * 1024;
    size = data.max * sizeof data.array[0];
    if ((data.array = malloc(size)) == NULL) {
        fatal_exit("malloc(%ld) failed", size);
    }
    data.count = 0;

    line[sizeof line - 1] = '\0';
    while (fgets(line, sizeof line, file) != NULL) {
        char *s1, *eol;
        double datum;

        // Sanity check on lines too long.
        if (line[sizeof line - 1] != '\0') {
            fatal_exit("line %ld: line too long", lineno);
        }

        // Avoid partial lines (usually at end of live log file).
        eol = strpbrk(line, "\r\n");
        if (eol == NULL) {
            fprintf(stderr, "%s: line %ld: truncated line (ignored) '%s'\n", Tool, lineno, line);
            break;
        }
        *eol = '\0';

        // Now have valid line, without terminator character(s).
        lineno += 1;
        s1 = line;

        // "yes" -- accept only lines containing this string.
        if (opt.yes && strstr(s1, opt.yes) == NULL) {
            FPRINTF(stderr, "%s: line %ld: skipped (missing `%s')\n", Tool, lineno, opt.yes);
            continue;
        }

        // "no" -- ignore all lines containing this string.
        if (opt.no && strstr(s1, opt.no) != NULL) {
            FPRINTF(stderr, "%s: line %ld: skipped (found `%s')\n", Tool, lineno, opt.no);
            continue;
        }

        // "begin" -- ignore all lines until this string appears.
        if (opt.begin) {
            if (match(s1, opt.begin)) {
                FPRINTF(stderr, "%s: line %ld: begin `%s')\n", Tool, lineno, opt.begin);
                opt.begin = NULL;
            } else {
                continue;
            }
        }

        // "after" -- ignore all lines until after this string appears.
        if (opt.after) {
            if (match(s1, opt.after)) {
                FPRINTF(stderr, "%s: line %ld: after `%s')\n", Tool, lineno, opt.after);
                opt.after = NULL;
            }
            continue;
        }

        // "before" -- stop when this string appears.
        if (opt.before) {
            if (match(s1, opt.before)) {
                FPRINTF(stderr, "%s: line %ld: before `%s')\n", Tool, lineno, opt.before);
                break;
            }
        }

        // Complain and ignore blank lines.
        if (s1[0] == '\0') {
            fprintf(stderr, "%s: line %ld: empty line (ignored)\n", Tool, lineno);
            continue;
        }

        // Ignore HP 53131A statistics (look for colon and space).
        if (opt.hp) {
            if (strstr(s1, ": ") != NULL) {
                continue;
            }
        }

        // Implement data decimation feature.
        count += 1;
        if (opt.mth > 0 && (count % opt.mth) != 1) {
            continue;
        }
        if (opt.nth > 0 && (count % opt.nth) != 0) {
            continue;
        }
        if (opt.skip > 0 && count <= opt.skip) {
            continue;
        }

        // Handle column number (column 1 is first column).
        if (opt.col > 1) {
            int len = strlen(s1);
            if (MINUS1(opt.col) > len) {
                fatal_exit("line %ld: column %d not present", lineno, opt.col);
            }
            s1 += MINUS1(opt.col);
        }

        // Replace commas with spaces.
        if (opt.comma) {
            int i;
            for (i = 0; s1[i] != '\0'; i += 1) {
                if (s1[i] == ',') s1[i] = ' ';
            }
        }

#define SKIP_WHITE(s) while (*s != '\0' &&  isspace(*s)) { s++; }
#define SKIP_BLACK(s) while (*s != '\0' && !isspace(*s)) { s++; }

        // Handle field number (field 1 is first field).
        if (opt.fld > 1) {
            int f;
            SKIP_WHITE(s1)
            for (f = 1; f < opt.fld; f += 1) {
                SKIP_BLACK(s1)
                SKIP_WHITE(s1)
            }
            if (s1[0] == '\0') {
                fatal_exit("line %ld: field %ld not present", lineno, opt.fld);
            }
        }

        // Handle substring of field (field n.1 is first character).
        if (MINUS1(opt.sub) > 0) {
            char *s2 = s1;
            SKIP_BLACK(s2)
            s1 += MINUS1(opt.sub);
            if (s1 >= s2) {
                fatal_exit("line %ld: field %ld.%ld not present", lineno, opt.fld, opt.sub);
            }
        }

        // Handle Hewlett-Packard/hp/Agilent/Keysight 53131A/53132A "user-friendly" data format.
        if (opt.hp) {
            int i = 0, j = 0;
            double scale = 1.0;

            // Remove thousand/thousandth commas (US version HP 53131A).
            do {
                if (s1[i] != ',') s1[j++] = s1[i];
            } while (s1[i++] != '\0');

            // Extract microsecond, second, MHz and Hz unit label.
            if (strstr(s1, " us") != NULL) {
                scale = 1e-6;
            } else if (strstr(s1, " u") != NULL) {
                scale = 1e-6;
            } else if (strstr(s1, " s") != NULL) {
                scale = 1.0;
            } else if (strstr(s1, " MHz") != NULL) {
                scale = 1e6;
            } else if (strstr(s1, " Hz") != NULL) {
                scale = 1.0;
            }
            datum = atof(s1) * scale;
        } else

        if (opt.oct) {
            long i32 = 0;
            sscanf(s1, "%lo", &i32);
            datum = i32;
        } else

        if (opt.hex) {
            unsigned long u32 = 0;
            int offset = (s1[0] == '0' && s1[1] == 'x') ? 2 : 0;
            sscanf(s1 + offset, "%lx", &u32);
            datum = u32;
        } else

        if (opt.hms) {
            long hhmmss, h, m, s;
            hhmmss = atol(s1);
            h = hhmmss / 10000 % 100;
            m = hhmmss /   100 % 100;
            s = hhmmss         % 100;
            datum = (((h * 60) + m) * 60) + s;
        } else

        if (opt.ddmm) {
            double ddmm = atof(s1);
            int deg = (int) (ddmm / 100);
            double min = ddmm - deg * 100;
            datum = deg + min / 60.0;

        } else {
            SKIP_WHITE(s1)
            if ( ! (s1[0] == '+' || s1[0] == '-' || isdigit(s1[0])) ) {
                fatal_exit("line %ld: non-numeric data `%s'", lineno, s1);
            }
            datum = atof(s1);
        }

        if (data.count == data.max) {
            data.max *= 2;
            size = data.max * sizeof data.array[0];
            if ((data.array = realloc(data.array, size)) == NULL) {
                fatal_exit("malloc(%ld) failed", size);
            }
        }
        data.array[data.count++] = datum;
        if (opt.take > 0 && data.count == opt.take) {
            break;
        }
    }
    if (data.count == 0) {
        fatal_exit("no data");
    }
}

// Return address of data array (read input file on first call).
double *get_data (void)
{
    if (data.array == NULL) {
        if (strcmp(opt.in, "-") == 0) {
            read_data(stdin);
        } else {
            FILE *f;
            if ((f = fopen(opt.in, "r")) == NULL) {
                fatal_exit("cannot open input file `%s'", opt.in);
            }
            read_data(f);
            fclose(f);
        }
    }
    return data.array;
}

// Return number of elements in data array.
long get_count (void)
{
    return data.count;
}

// Allocate new data array.
double *new_data (long count)
{
    size_t size = count * sizeof (double);
    double *newdata;

    if (count == 0) {
        fatal_exit("no data");
    }
    if ((newdata = malloc(size)) == NULL) {
        fatal_exit("malloc(%ld) failed", size);
    }
    return newdata;
}

// Commit new data array.
void set_data (double *newdata, long count)
{
    free(data.array);
    data.array = newdata;
    data.count = count;
    data.max = count;
}

// Write output file.
void write_output (void)
{
    FILE *f;
    long i;

    if (data.array == NULL || data.count == 0) {
        fatal_exit("cannot write output (no data)");
    }
    if (strcmp(opt.out, "-") == 0) {
        f = stdout;

    } else {
        if ((f = fopen(opt.out, "w")) == NULL) {
            fatal_exit("cannot open output file `%s'", opt.out);
        }
    }
    for (i =  0; i < data.count; i += 1) {
        if (opt.plus && data.array[i] >= 0) fprintf(f, "+");
        fprintf(f, opt.format, data.array[i]);
    }
    if (f != stdout) {
        fclose(f);
    }
}

// Allow from,to or from,count.
int get_from_to (long *from, long *to)
{
    long count = get_count();
    *from = opt.from;
    *to = count;
    if (opt.to == 0 && opt.count != 0) *to = opt.from + opt.count;
    if (opt.to != 0 && opt.count == 0) *to = opt.to + 1;
    if (*from < 0 || *from > count || *to < 0 || *to > count || *from >= *to) {
        FPRINTF(stderr, "%s: (%ld %ld %ld) from %ld to %ld count %ld\n", Tool, opt.from, opt.to, opt.count, *from, *to, count);
        return 0;
    }
    return 1;
}

double mean_sdev (double *data, long n, double *sdev, int pop)
{
    long i;
    double sum1 = 0, mean;

    for (i = 0; i < n; i += 1) {
        sum1 += data[i];
    }
    mean = sum1 / n;

    if (sdev != NULL) {
        double sum2 = 0;
        for (i = 0; i < n; i += 1) {
            double d = data[i] - mean;
            sum2 += d * d;
        }
        if (pop) {
            *sdev = (n < 1) ? 0 : sqrt(sum2 / n);        // standard deviation (population)
        } else {
            *sdev = (n < 2) ? 0 : sqrt(sum2 / (n - 1));  // standard deviation (sample), default
        }
    }
    return mean;
}

//
// c[0] = mean
//

void fit0 (double *data, long n, double c[1])
{
    long i;
    double sum;

    sum = 0.0;
    for (i = 0; i < n; i += 1) {
        sum += data[i];
    }
    c[0] = sum / n;
}

#define XDATA(i) i
#define YDATA(i) data[i]

//
// Linear least squares fit.
//
// c[0] = offset, c[1] = slope
//

void fit1 (double *data, long n, double c[2])
{
    long i;
    double dx, dy, Sx, Sy, Sdx2, Sdxdy, xmean, ymean;

    Sx = Sy = 0.0;
    for (i = 0; i < n; i += 1) {
        Sx += XDATA(i);
        Sy += YDATA(i);
    }
    xmean = Sx / n;
    ymean = Sy / n;

    Sdx2 = Sdxdy = 0.0;
    for (i = 0; i < n; i += 1) {
        dx = XDATA(i) - xmean;
        dy = YDATA(i) - ymean;
        Sdx2  += dx * dx;
        Sdxdy += dx * dy;
    }
    c[1] = Sdxdy / Sdx2;         // slope
    c[0] = ymean - c[1] * xmean; // offset
}

//
// Quadratic least squares fit.
//
// c[0] = offset, c[1] = slope, c[2] = drift
//

void fit2 (double *data, long n, double c[3])
{
    int i;
    double xmean, ymean, denom;
    double Sx0, Sx1, Sx2, Sx3, Sx4;
    double Sy1, Syx0, Syx1, Syx2;
    double xi, yi;

    // Remove mean from x and y to improve floating point precision.
    Sx1 = Sy1 = 0.0;
    for (i = 0; i < n; i += 1) {
        Sx1 += XDATA(i);
        Sy1 += YDATA(i);
    }
    xmean = Sx1 / n;
    ymean = Sy1 / n;

    // Compute sums.
    Sx0 = Sx1 = Sx2 = Sx3 = Sx4 = 0.0;
    Syx0 = Syx1 = Syx2 = 0.0;

    for (i = 0; i < n; i += 1) {

        xi = XDATA(i) - xmean;
        yi = YDATA(i) - ymean;

        Sx0 += 1.0;
        Sx1 += xi;
        Sx2 += xi * xi;
        Sx3 += xi * xi * xi;
        Sx4 += xi * xi * xi * xi;

        Syx0 += yi;
        Syx1 += yi * xi;
        Syx2 += yi * xi * xi;
    }

    // Compute coefficients.
    denom = Sx0 * (Sx2 * Sx4 - Sx3 * Sx3)
          - Sx1 * (Sx1 * Sx4 - Sx2 * Sx3)
          + Sx2 * (Sx1 * Sx3 - Sx2 * Sx2);

    c[0] = ( Syx0 * (Sx2 * Sx4 - Sx3 * Sx3)
           + Syx1 * (Sx2 * Sx3 - Sx1 * Sx4)
           + Syx2 * (Sx1 * Sx3 - Sx2 * Sx2) ) / denom;

    c[1] = ( Syx0 * (Sx2 * Sx3 - Sx1 * Sx4)
           + Syx1 * (Sx0 * Sx4 - Sx2 * Sx2)
           + Syx2 * (Sx1 * Sx2 - Sx0 * Sx3) ) / denom;

    c[2] = ( Syx0 * (Sx1 * Sx3 - Sx2 * Sx2)
           + Syx1 * (Sx2 * Sx1 - Sx0 * Sx3)
           + Syx2 * (Sx0 * Sx2 - Sx1 * Sx1) ) / denom;

    // Adjust coefficients to back out above mean removal.
    c[0] = c[2] * xmean * xmean - c[1] * xmean + ymean;
    c[1] -= 2 * c[2] * xmean;
}