//========================================================= file = rs.c =====
//=  Program to compute the R/S statistic for a series X of size N          =
//=    > Used to estimate the self-similarity Hurst parameter (H)           =
//===========================================================================
//=  Notes:                                                                 =
//=    1) Input from input file "in.dat" to stdin (see example below)       =
//=        * Comments are bounded by "&" characters at the beginning and    =
//=          end of the comment block                                       =
//=    2) Output is to stdout                                               =
//=    3) X should have a "large number" of values for a "good" R/S value   =
//=       to be computed.                                                   =
//=-------------------------------------------------------------------------=
//=  Example "in.dat" file:                                                 =
//=                                                                         =
//=    & Sample series of data which can be integers or reals.              =
//=      There are 6 values in this file. &                                 =
//=    12                                                                   =
//=    56                                                                   =
//=    99                                                                   =
//=   111                                                                   =
//=    87                                                                   =
//=    99                                                                   =
//=-------------------------------------------------------------------------=
//=  Example output (for above "in.dat"):                                   =
//=                                                                         =
//=   ---------------------------------------------------- rs.c -----       =
//=     R/S = 2.557638 for series X of 6 values                             =
//=   ---------------------------------------------------------------       =
//=-------------------------------------------------------------------------=
//=  Build: bcc32 rs.c, cl rs.c, gcc rs.c -lm                               =
//=-------------------------------------------------------------------------=
//=  Execute: rs < in.dat                                                   =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (09/16/98) - Genesis                                      =
//=           KJC (09/05/00) - Removed inside double loop per advice from   =
//=                            Sunwoo Lee.  Program runs much faster.       =
//=           KJC (10/31/00) - Fixed a one-off error in main loop.  See     =
//=                            "Fix #1" tag.                                =
//=           KJC (05/20/02) - Renamed to rs.c                              =
//===========================================================================

//----- Include files -------------------------------------------------------
#include <stdio.h>                 // Needed for printf() and feof()
#include <math.h>                  // Needed for pow() 
#include <stdlib.h>                // Needed for exit() and atof()
#include <string.h>                // Needed for strcmp()

//----- Defines -------------------------------------------------------------
#define MAX_SIZE 1000000L          // Maximum size of time series data array

//----- Globals -------------------------------------------------------------
double     X[MAX_SIZE];            // Time series read from "in.dat"
long int   N;                      // Number of values in in.dat

//----- Function prototypes -------------------------------------------------
void   load_X_array(void);         // Load X array from "in.dat"
double compute_rs(void);           // Compute R/S for X of length N

//===========================================================================
//=  Main program                                                           =
//===========================================================================
void main(void)
{
  double  rs_value;                // Computed R/S value

  // Load the series X
  printf("---------------------------------------------------- rs.c -----\n");
  load_X_array();

  // Compute R/S value for series X of length N
  rs_value = compute_rs();

  // Output R/S value
  printf("  R/S = %f for series X of %ld values \n", rs_value, N);
  printf("---------------------------------------------------------------\n");
}

//===========================================================================
//=  Function to load X array from stdin and determine N                    =
//===========================================================================
void load_X_array(void)
{
  char      temp_string[1024];     // Temporary string variable

  // Read all values into X
  N = 0;
  while(1)
  {
    scanf("%s", temp_string);
    if (feof(stdin)) goto end;

    // This handles a comment bounded by "&" symbols
    while (strcmp(temp_string, "&") == 0)
    {
      do
      {
        scanf("%s", temp_string);
        if (feof(stdin)) goto end;
      } while (strcmp(temp_string, "&") != 0);
      scanf("%s", temp_string);
      if (feof(stdin)) goto end;
    }

    // Enter value in array and increment array index
    X[N] = atof(temp_string);
    N++;

    // Check if MAX_SIZE data values exceeded
    if (N >= MAX_SIZE)
    {
      printf("*** ERROR - greater than %ld data values \n", MAX_SIZE);
      exit(1);
    }
  }

  // End-of-file escape
  end:

  return;
}

//===========================================================================
//=  Function to compute R/S value for series X of length N                 =
//===========================================================================
double compute_rs()
{
  double     mom1;      // First moment
  double     mom2;      // Second moment
  double     x_bar;     // Mean (X_bar value)
  double     s;         // Standard deviation (S value)
  double     w;         // W value
  double     r;         // R value
  double     min_w;     // Minimum W value
  double     max_w;     // Maximum W value
  double     rs_value;  // R/S value to be returned
  double     sum;       // Temporary sum value
  long int   i, j;      // Loop counters

  // Loop to compute mean and standard deviation of X
  mom1 = mom2 = 0.0;
  for (i=0; i<N; i++)
  {
    mom1 = mom1 + (X[i] / N);
    mom2 = mom2 + (pow(X[i], 2.0) / N);
  }
  x_bar = mom1;
  s = sqrt(mom2 - pow(mom1, 2.0));

  // Double loop to find minimum and maximum W values
  min_w = max_w = 0.0;
  sum = 0.0;
  for (i=0; i<N; i++)
  {
    sum = sum + X[i];

    w = sum - ((i+1) * x_bar);         // Fix #1
    if (w > max_w) max_w = w;
    if (w < min_w) min_w = w;
  }

  // Compute R value as maximum W minus minimum W
  r = max_w - min_w;

  // Compute R/S value
  rs_value = r / s;

  // Return R/S value
  return(rs_value);
}