//
// VCO EFC linearity analysis tool.
//
// - Input file is volt,freq data from SR620.
// - Plot output file(s) with Excel.
// - Hz(V) dF(V) Hz/V(V) dF/V(V)
// - Usage:
//   efc1 log12509.txt
//   efc1 log12509.txt /f0:15728264
//   efc1 log12509.txt /f0:15728264 /box:20 /col:2 NUL rate.txt
//
// 17-Apr-2003 tvb
//

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

//
// Define statistics structures.
//

typedef struct {
    double *data;
    double count;
    double min;
    double max;
    double range;
    double first;
    double last;
    double span;
    double mean;
    double std_dev;
    double normal;
    char *units;
} STAT1, *PSTAT1;

typedef struct {
    PSTAT1 x;
    PSTAT1 y;
    double slope;
    double offset;
    double chi_square;
    char *units;
} STAT2, *PSTAT2;

//
// Global data (zero initialized).
//

STAT1 stEFC, stFREQ, stRESID, stRATE;
STAT2 stSLOPE;

FILE *Out1;
FILE *Out2;

struct {
    double f0;
    long box;
    int column;
} Options;

//
// Define forward-referenced function prototypes.
//

FILE *OpenFile (char *Path, char *Comment);
long ReadData2 (PSTAT1 x, PSTAT1 y);
void Stat1 (PSTAT1 s, char *title);
void Stat2 (PSTAT2 s, char *title);
void Stat3 (PSTAT2 s, PSTAT1 resid);
void Stat4 (PSTAT1 x, PSTAT1 y, PSTAT1 rate);

//
// Read data and generate statistics.
//

void main (int argc, char *argv[])
{
    long count;

    //
    // Parse options.
    //

    Options.f0 = 0;
    Options.box = 50;
    Options.column = -1;
    while (argc > 1 && argv[1][0] == '/') {
        if (strncmp(argv[1], "/f0:", 4) == 0) {
            sscanf(&argv[1][4], "%lf", &Options.f0);
            fprintf(stderr, "Normalizing data with: %lg\n", Options.f0);
        } else

        if (strncmp(argv[1], "/box:", 5) == 0) {
            sscanf(&argv[1][5], "%ld", &Options.box);
            fprintf(stderr, "Calculate slope using box size: %ld\n", Options.box);
        } else

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

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

    //
    // Open optional output files.
    //

    if (argc > 1) {
        Out1 = OpenFile(argv[1], "residuals");
    }

    if (argc > 2) {
        Out2 = OpenFile(argv[2], "EFC rate");
    }

    //
    // Read data into EFC and FREQ arrays.
    //

    stEFC.units = "V";
    stFREQ.units = "Hz";
    stFREQ.normal = Options.f0;
    count = ReadData2(&stEFC, &stFREQ);

    //
    // Compute basic statistics.
    //

    Stat1(&stEFC, "EFC data set:");
    Stat1(&stFREQ, "Frequency data set:");

    //
    // Compute V to F transfer function (freq/volts slope)
    //

    stSLOPE.x = &stEFC;
    stSLOPE.y = &stFREQ;
    stSLOPE.units = "Hz/V";
    Stat2(&stSLOPE, "EFC slope:");

    //
    // Compute residuals after linear fit.
    //

    stRESID.count = count;
    stRESID.units = "Hz (resid)";
    stRESID.data = (double *) malloc(sizeof(double) * count);
    Stat3(&stSLOPE, &stRESID);
    Stat1(&stRESID, "Frequency residuals output:");

    //
    // Differentiate FREQ, examine variations in VtoF slope (EFC gain).
    //

    if (Options.box != 0) {
        stRATE.count = count / Options.box;
        stRATE.units = "Hz/V";
        stRATE.data = (double *) malloc(sizeof(double) * count / Options.box);
        Stat4(&stEFC, &stFREQ, &stRATE);
        Stat1(&stRATE, "EFC rate output:");
    }

    return;
}

//
// Common function to open output file.
//

FILE *OpenFile (char *Path, char *Comment)
{
    FILE *f;

    if (strcmp(Path, "-") == 0) {
        f = stdout;
        fprintf(stderr, "Writing %s to stdout\n", Comment);

    } else {
        f = fopen(Path, "w");
        if (f == NULL) {
            fprintf(stderr, "Open failed: %s\n", Path);
            exit(1);
        }
        fprintf(stderr, "Writing %s to file: %s\n", Comment, Path);
    }

    return f;
}

//
// Read raw data into two auto-sizing arrays.
//

long ReadData2 (PSTAT1 x, PSTAT1 y)
{
    long index;
    long limit;
    char Line[100];
    int n;

    x->data = NULL;
    y->data = NULL;
    limit = 0;

    index = 0;
    while (fgets(Line, sizeof(Line), stdin) != NULL) {
        if (index == limit) {
            limit = 1000 + limit * 2;
            x->data = (double *) realloc(x->data, sizeof(double) * limit);
            y->data = (double *) realloc(y->data, sizeof(double) * limit);
        }
        n = sscanf(Line, "%lf %lf", &x->data[index], &y->data[index]);
        if (n != 2) {
            fprintf(stderr, "SKIPPING: input format error: %s\n", Line);
        }
        if (Options.f0 != 0.0) {
            y->data[index] = (y->data[index] - Options.f0) / Options.f0;
        }
        index += 1;
    }

    x->count = index;
    y->count = index;

    return index;
}

char *fmt1 = "%12s: %15.6lf %s\n";
char *fmt1e = "%12s: %15.6le %s\n";

void DISPLAY1(char *name, double value, char *units)
{
    if (value >= 0.1) {
        printf(fmt1, name, value, units);

    } else {
        printf(fmt1e, name, value, units);
    }
}

char *fmt2 = "%.3lf ";
char *fmt2e = "%.3le ";

void DISPLAY2(FILE *f, double value)
{
    if (value >= 0.1) {
        fprintf(f, fmt2, value);

    } else {
        fprintf(f, fmt2e, value);
    }
}

#define SQUARE(x) ((x) * (x))

//
// Calculate basic statistics.
//

void Stat1 (PSTAT1 s, char *title)
{
    long i, n;
    double x, dx, sum;

    n = (long) s->count;
    s->min = s->max = sum = s->data[0];
    for (i = 1; i < n; i += 1) {
        x = s->data[i];
        if (x < s->min) s->min = x;
        if (x > s->max) s->max = x;
        sum += x;
    }
    s->mean = sum / n;
    s->range = s->max - s->min;
    s->first = s->data[0];
    s->last = s->data[n - 1];
    s->span = s->last - s->first;

    sum = 0.0;
    for (i = 0; i < n; i += 1) {
        dx = s->data[i] - s->mean;
        sum += SQUARE(dx);
    }

    if (n > 1) {
        s->std_dev = sqrt(sum / (n - 1));
    } else {
        s->std_dev = 0;
    }

    printf("\n%s\n", title);
    DISPLAY1("count", s->count, "");
    // DISPLAY1("first", s->first, s->units);
    // DISPLAY1("last", s->last, s->units);
    // DISPLAY1("span", s->span, s->units);
    DISPLAY1("min", s->min, s->units);
    DISPLAY1("max", s->max, s->units);
    DISPLAY1("range", s->range, s->units);
    DISPLAY1("mean", s->mean, s->units);
    DISPLAY1("std_dev", s->std_dev, s->units);
    if (s->normal) {
        DISPLAY1("normalized", s->normal, s->units);
    }

    return;
}

//
// Calculate slope and chi square (linear least squares fit).
//

void Stat2 (PSTAT2 s, char *title)
{
    long i, n;
    double dx, dy;
    double Sdx2, Sdy2, Sdxdy;

    n = (long) s->x->count;
    Sdx2 = Sdy2 = Sdxdy = 0.0;
    for (i = 0; i < n; i += 1) {
        dx = s->x->data[i] - s->x->mean;
        dy = s->y->data[i] - s->y->mean;
        Sdx2 += SQUARE(dx);
        Sdy2 += SQUARE(dy);
        Sdxdy += dx * dy;
    }

    if (Sdx2 != 0.0 && Sdy2 != 0.0) {
        s->slope = Sdxdy / Sdx2;
        s->offset = s->y->mean - (s->slope * s->x->mean);
        s->chi_square = SQUARE(Sdxdy / sqrt(Sdx2 * Sdy2));
    } else {
        s->slope = 0.0;
        s->offset = 0.0;
        s->chi_square = 0.0;
    }

    printf("\n%s\n", title);
    DISPLAY1("slope", s->slope, s->units);
    DISPLAY1("offset", s->offset, s->units);
    DISPLAY1("chi_square", s->chi_square, "");

    return;
}

//
// Remove slope, compute residuals, and write output file.
//

void Stat3 (PSTAT2 s, PSTAT1 resid)
{
    long i, n;
    double predict;

    n = (long) s->x->count;
    for (i = 0; i < n; i += 1) {
        predict = (s->slope * s->x->data[i]) + s->offset;
        resid->data[i] = s->y->data[i] - predict;
        if (Out1) {
            if (Options.column & 1) DISPLAY2(Out1, s->x->data[i]);
            if (Options.column & 2) DISPLAY2(Out1, s->y->data[i]);
            if (Options.column & 4) DISPLAY2(Out1, predict);
            if (Options.column & 8) DISPLAY2(Out1, resid->data[i]);
            if (Options.column & 15) fprintf(Out1, "\n");
        }
    }

    return;
}

//
// Compute rate and write output file.
//

void Stat4 (PSTAT1 x, PSTAT1 y, PSTAT1 rate)
{
    long i, n, box;
    double Lsq (double *x, double *y, long n);

    n = (long) x->count;
    box = Options.box;
    for (i = 0; i < n / box; i += 1) {
        rate->data[i] = Lsq(&x->data[i * box], &y->data[i * box], box);
        if (Out2) {
            if (Options.column & 1) DISPLAY2(Out2, x->data[i * box]);
            if (Options.column & 2) DISPLAY2(Out2, rate->data[i]);
            if (Options.column & 3) fprintf(Out2, "\n");
        }
    }

    return;
}

double Lsq (double *x, double *y, long n)
{
    long i;
    double dx, dy;
    double Sdx2, Sdy2, Sdxdy;
    double xmean, ymean;
    double Mean (double *x, long n);

    xmean = Mean(x, n);
    ymean = Mean(y, n);
    Sdx2 = Sdy2 = Sdxdy = 0.0;
    for (i = 0; i < n; i += 1) {
        dx = x[i] - xmean;
        dy = y[i] - ymean;
        Sdx2 += SQUARE(dx);
        Sdy2 += SQUARE(dy);
        Sdxdy += dx * dy;
    }

    if (Sdx2 == 0.0 || Sdy2 == 0.0) {
        return 0.0;
    }

    return (Sdxdy / Sdx2);
}

double Mean (double *x, long n)
{
    long i;
    double sum;

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

    return (sum / n);
}