/*
    Copyright 2023 Harm Jetten

    This file is part of Moby Dam.

    Moby Dam is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Moby Dam is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Moby Dam.  If not, see <http://www.gnu.org/licenses/>.
*/

/* weightstats.c: show statistics of generated weights */

#include "train.h"

bool debug_info;               /* print extra debug info */
double limit_val = 1e9;

/* calculate weight statistics */
/* file = weights file path/name */
/* file contains one floating point weight value per line */
static void weight_stats(char *file)
{
    FILE *fp;
    char line[80];
    char *str;
    int i, n, nz, neg, lim;
    double w, lo, hi;
    double delta, mean, delta2, m2;

    fp = fopen(file, "r");
    if (fp == NULL)
    {
        fprintf(stderr, "error opening %s\n", file);
        return;
    }
    nz = neg = lim = 0;
    lo = hi = 0.0;
    mean = m2 = 0.0;
    for (i = 0; ; i++)
    {
        str = fgets(line, sizeof line, fp);
        if (str == NULL)
        {
            /* reached end of input file */
            fclose(fp);
            break;
        }
        while (*str != '\0')      /* get rid of any brackets present */
        {
            if (*str == '[' || *str == ']')
            {
                *str = ' ';
            }
            str++;
        }
        n = sscanf(line, "%lf", &w);
        if (n != 1)
        {
            /* reached end of input file */
            fclose(fp);
            break;
        }

        /* find extremes */
        lo = min(lo, w);
        hi = max(hi, w);

        /* count non-zero weights */
        if (w != 0.0)
        {
            nz++;

            if (w < 0.0)
            {
                neg++;
            }

            if (w < -limit_val || w > limit_val)
            {
                lim++;
            }

            /* one-pass variance/stdev calculation */
            if (nz == 1)
            {
                mean = w;
                m2 = 0.0;
            }
            else
            {
                delta = w - mean;
                mean += delta/nz;
                delta2 = w - mean;
                m2 += delta*delta2;
            }
        }
    }

    /* print interesting statistics */
    printf("%s: count=%d nonzero=%d neg=%d pos=%d lim=%d low=%lf high=%lf "
           "nz_mean=%lf nz_var=%lf nz_stdev=%lf\n",
           file, i, nz, neg, nz - neg, lim, lo, hi, mean, m2/i, sqrt(m2/i));
    return;
}

/* the program entry point */
int main(int argc, char *argv[])
{
    int opt;

    while (TRUE)
    {
        opt = getopt(argc, argv, "dx:");
        if (opt == -1) 
        {
            break;
        }
        switch (opt) 
        {
        case 'd':
            debug_info = TRUE;
            break;
        case 'x':
            limit_val = atof(optarg);
            break;
        default:
            printf("Usage: %s [-d] [-x limit] file...\n", argv[0]);
            printf("  -d = print lots of extra debug info\n"
                   "  -x limit = count weights whose absolute value exceeds given limit\n"
                   "  file... = one or more weights files\n");
            exit(EXIT_FAILURE);
        }
    }

    while (optind < argc)
    {
        weight_stats(argv[optind++]);
    }
    return EXIT_SUCCESS;
}
