/* ---------------------------------------------------------------------- 
**  Figure 16.8:  Computing better statistics (main()).
**  Figure 16.6:  Using Chauvenet's ratio.
**  Figure 16.5:  Searching the Chauvenet table  (chauv_ratio()).  
**  Figure 16.9:  From file to masked array.
**  Figure 16.10: The print_stats() function.
**  Figure 16.11: The chauv_ar() function: Creating the mask array.
**  Figure 16.12: The stats() function.
**  Figure 16.13: The print_data() function.
** ---------------------------------------------------------------------- 
*/
#include "tools.h"

#define N  500          /* maximum number of data values  */
typedef struct VALUE { double data; bool mask; } value_t;
typedef struct STATS { double mean, var, stdev; } stat_t;

int    read_file   ( value_t x[], int nmax );
void   print_data  ( value_t x[], int n );
stat_t stats       ( int n, int n_good, value_t x[] );
void   print_stats ( string when, int n, stat_t st );
int    chauv_ar    ( int n, value_t x[], stat_t st );
bool   chauv_test  ( double x, stat_t st, double ratio );
double chauv_ratio ( int n );

/* ---------------------------------------------------------------------- */
void main( void )
{   
    value_t x[N];            /* array for the data values */
    int num;                 /* initial number of data values */
    int n_good;              /* # of values after eliminating "bad" ones. */
    stat_t answer;           /* Statistical results for array x. */ 

    banner();
    num = read_file( x, N );                /* read data from file */
    if (num < 3) fatal( " Too few data items for meaningful analysis." ); 

    answer = stats( num, num, x );          /* compute initial stats */
    print_stats( "before", num, answer );   /* and print them */

    n_good = chauv_ar( num, x, answer );    /* test values */
    answer = stats( num, n_good, x );       /* recompute stats */
    print_stats( "after", n_good, answer ); /* print new stats */

    print_data( x, num );                   /* print results */
    bye();
}

/* ---------------------------------------------------------------------- */
/* Assumes n, the actual number of data points, is 3 or greater.          */
double                          /* Return Chauvenet discrimination ratio. */
chauv_ratio( int n )
{
    int k;                              /* index for search loop */
    typedef struct {int n_pts; double ratio; } table_type;
    #define L  9                       /* Length of Chauvenet's table */
    static const table_type chauv[L] = { {  3, 1.38}, {  5, 1.65}, { 10, 1.96},
                                         { 15, 2.13}, { 25, 2.33}, { 50, 2.57},
                                         {100, 2.81}, {300, 3.14}, {500, 3.29}
                                        };

    /* Compare actual number of data points to numbers in table. */
    /* Stop search when proper place in table is located.  */
    for( k=0; k<L; ++k ) {
        if (n < chauv[k].n_pts) break;
    }
    return chauv[k-1].ratio;    /* We've gone too far; back up one slot. */
}

/* ---------------------------------------------------------------------- */
bool
chauv_test( double x, stat_t st, double ratio )
{
    double dr = fabs( (x - st.mean)/st.stdev ); /* deviation ratio */
    /* Compare ratio to ratio from Chauvenet's table.           */
    /* Return true if dr is small enough, false otherwise.      */
    return ratio > dr;
}

/* --------------------------------------------------------------------- */
int
read_file( value_t x[], int nmax )
{
    int k;              /* subscript variable for loop */
    char filename[80];  /* name of input file */
    stream fin;         /* stream for input file */
    int flag;           /* for fscanf() return value */

    printf( " Name of data file: " );
    scanf( "%s", filename );
    fin = fopen( filename, "r" );
    if (fin == NULL) fatal( "\n Fatal error: cannot open %s", filename );
    printf( " Reading data from %s.\n", fin );  

    for (k = 0; k < nmax; ) {   /* Read data, initialize "good" indicators. */
        flag = fscanf( fin, "%lg", &x[k].data );
        if (flag == 1) {
            x[k].mask = true;   /* Point is "good" until proven bad. */
            ++k;                /* All is well -- count the item. */
        }
        else if (feof( fin )) break;        
        else {
            say( " - -  Bad data while reading slot %i: ", k );
            clean_and_log( fin, stderr );
        }
    }
    fclose( fin );
    return k;
}

/* --------------------------------------------------------------------- */
void
print_stats( string when, int n, stat_t st )
{
    printf( " Stats %s eliminating bad values:\n"
            "\t For %i data values \n\t mean = %.2f \n"
            "\t variance = %.2f \n\t standard deviation = %.2f \n\n",
            when, n, st.mean, st.var, st.stdev );
}

/* ---------------------------------------------------------------------- */
int                                 /* Return number of values to keep */
chauv_ar( int n, value_t x[ ], stat_t st )
{   
    double ratio;                   /* Chauvenet ratio for n values */
    int k;                          /* subscript for data table */
    int n_good = 0;                 /* number of good values */
    
    ratio = chauv_ratio(n);         /* Get proper ratio for n values */

    for (k = 0; k < n; ++k) {       /* Analyze values one at a time */  
        x[k].mask = chauv_test( x[k].data, st, ratio );
        if (x[k].mask) ++n_good;    /* count number of keepers */
    }
    return n_good;
}

/* --------------------------------------------------------------------- */
stat_t 
stats( int n, int n_good, value_t x[] )
{
    int k;                /* subscript for data array */
    double sum1;          /* sum of values */
    double sum2;          /* sum of squares */
    double divisor;
    stat_t local;         /* local storage for answers */

   for (sum1 = k = 0; k < n; ++k)           /* Sum "good" data values. */
        if ( x[k].mask ) sum1 += x[k].data;  

    local.mean = sum1 / n_good;             /* Mean of "good" values. */

    for (sum2 = k = 0; k < n; ++k)          /* Sum difference^2 of "goods". */
        if (x[k].mask) sum2 += pow( (x[k].data - local.mean), 2 );

    if (n_good < 20) divisor = n_good-1;  
    else             divisor = n_good;
    local.var = sum2 / divisor;     
    local.stdev = sqrt( local.var );
    return local;                   /* Return struct containing the stats. */
}

/* --------------------------------------------------------------------- */
void
print_data( value_t x[], int n )
{
    int k;              /* loop counter */
    string answer;      /* set to "good" or "bad" for output */

    printf( " Evaluation of data:\n" );     
    for (k = 0; k < n; ++k) {                    
        if ( x[k].mask ) answer = "good";
        else answer = "bad";
        printf( "\t x[%2i] = %.2f is %s\n", k, x[k].data, answer );
    }
}