/* ------------------------------------------------------------------
 * Generate a string of letters S letters long, with a pressure P,
 * 1 - 50, of adding randomness to "Hello world".
 * 0 is no change from  "Hello world"
 * 1  should be close to "Hello world"
 * 50 should be closer to random
 * The mangled version of "Hello world" will be inserted into another
 * string of length S. 
 * --------------------------------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <malloc.h>
#include <time.h>

#define REFARRAY_LEN 11

void Usage  (void) ;
void Info   (void) ;

void  ScrambleInPlace      (int len, int* thestring, int pressure) ; 

/* ---------------------------------------------------------------
 * I want a reusable pearson's r program, and that should take
 * arrays of digits. 
 * But what I want for this test program is to work with a character
 * array. It turns out to be really very quite tediously awkwardly
 * gauche to keep inserting \0 and checking for it just for the
 * same of printing. So the arrays will be converted to printable
 * strings at the last moment.
 */
float pearsons_r            (int N, int* X, int* Y)  ;


float HelloWorldCorrelator (int* l_t_o, int* thestring, int* refstring) ; 
                                               // l_t_o means length, then offset. It's a 
                                               // pass by reference. Details in HelloWorldCorrelator


/* *****************************************************************************
 *                                  MAIN
 * ************************************************************************** */
int main (int argc, char *argv[])
{
  // reference "Hello world" 
  int HWR [REFARRAY_LEN] = {72, 101, 108, 108, 111, 32, 119, 111, 114, 108, 100 } ;
  // array of  "Hello world" to mangle
  int HW  [REFARRAY_LEN] = {72, 101, 108, 108, 111, 32, 119, 111, 114, 108, 100 } ;
  int STR[128] ;  // re-determined later by S
  int i, j, l, S, P ; 
  int l_t_o ;
  int printcsv = 0 ;

  for (i=1;i<argc;i++) {
    if (argv[i][0] == '-') {
      switch (argv[i][1]) {
        case 'S': S = atoi(argv[++i]) ; break ;
        case 'P': P = atoi(argv[++i]) ; break ;
        case 'c': printcsv = 1 ; break ;
  } } }

  /* ------------- input checking --------------------------*/ 
  if (S > 128)            { S = 128 ;}           // S can't be > 128
  if (S < REFARRAY_LEN)   { S = REFARRAY_LEN ;}  // S can't be shorter than HW
  if (P > 50)  { P = 50 ;}  // P can't be > 50
  if (P < 0)   { P =  0 ;}  // P can't be < 0
  /* ------------- input checking --------------------------*/ 


  l_t_o = S ;
  srand (time (NULL));
  ScrambleInPlace(REFARRAY_LEN, HW, P) ;
  //       0 to 10
  for (j = 0; j < S ; j++ ) { STR[j] = 33 + (rand() % (126-33)) ; }
  i = rand() % (S-10) ;  // Find a starting point between 0 and end of S to insert HW
  j = 0 ;
  int k = i ;
  while (i < k+11) {
   STR[i++] = HW[j++] ;
  }

  // pearsons_r for crosscheck w/spreadsheet
  if (printcsv) 
  {
    printf ("\nSTR [\n") ; for (i=0;i<S-1;i++) { printf("%d,\n",STR[i]) ;} printf ("%d\n]\n",STR[i]) ;
    printf ("\nHW  [\n") ; for (i=0;i<10;i++)  { printf("%d,\n",HW[i]) ;}  printf ("%d\n]\n",HW[i]) ;
  }

  printf ("\nCC max with HWR is %1.3f", HelloWorldCorrelator (&l_t_o, STR, HWR) ) ;
  printf (" at the fragment starting at offset %d. Offsets start at 0.\n",l_t_o) ;

  printf ("\n") ; 
  return (0) ;
}



/* ---------------------------------------------------------------------- */
void ScrambleInPlace (int len, int* s, int pressure) 
{
  /* This is a multipass scramble. It does a per character random walk, 
   * P times. The range is the ASCII character set 32 - 126, so that's 
   * the space character, plus the printable characters. 
   * The walk saturates top and bottom. So random is 3-way: increment, 
   * decrement, and no-change; unless it's at the wall of 32 or 126, and
   * in that case it's 2-way: stay or decrement (at 32), stay or 
FIXME comment is wrong
   * increment (at 126). 
   * The string is scrambled in place.
   */
  int i, j, r;

  if (pressure > 0) {
    for (j=0;j<len;j++) {
      r = rand()%(2*pressure +1 ) - pressure ; 
      s[j] = s[j] + rand()%(2*pressure +1) - pressure ;
      if      (s[j] < 32)  { s[j] =  32; }
      else if (s[j] > 126) { s[j] = 126; }
    }
  }
}



/* ---------------------------------------------------------------------- */
float HelloWorldCorrelator (int* l_t_o, int* s, int* ref) 
{
  /* This routine calculates the correlation coefficient, CC, of the reference
   * string ref and all possible same-length, in-order, no-wrap strings in s.
   * I say "strings", but these are character arrays limited to values 
   * 32 to 126.
   * It keeps the highest |CC| and returns it. It also returns, via reference,
   * the index of the first element in s where the reference ref showed the
   * highest CC. It does that by overwriting l_t_o.
   */
  float CC = 0.0 ;
  float CCtemp = 0.0 ;

  int i, j, indx ;
  indx = 0 ;

  int thisY[REFARRAY_LEN] ;

  // i scoots the array r from 0 to l_t_o - length(HW) along s
  for (i=0; i<(*l_t_o - REFARRAY_LEN); i++)
  {
    for (j=0; j<REFARRAY_LEN; j++)
    { 
      thisY[j] = s[i+j] ; 
    }

    CCtemp = pearsons_r(REFARRAY_LEN, ref, thisY) ;

    if (fabs(CCtemp) > fabs(CC)) 
    {
      CC = fabs(CCtemp) ;
      indx = i ;
    }

  }
  *l_t_o = indx ;
  return CC ;
}



/* ----------------------------------------------------------------------
 * Straight-up Pearson’s r correlation coefficient.
 * Assumptions:
 * 1. Expects that X and Y actually have N entries. Not responsible for segfault.
 * 2. Expects you to ensure that N*X[i]max*Y[i]max < double. Overflow not checked.
 */
float pearsons_r (int N, int* X, int* Y) 
{
  float CC ;
  int i ;

  double XY[REFARRAY_LEN], Xsquared[REFARRAY_LEN], Ysquared[REFARRAY_LEN];
  double sumX, sumY, XYsum, sumXsqrd, sumYsqrd;
  float  numer, denom;

  sumX = sumY = XYsum = sumXsqrd = sumYsqrd = 0;
  CC = 0.0 ;

  for (i = 0; i<N ; i++)
  {
    XY[i]       = X[i] * Y[i];
    Xsquared[i] = X[i] * X[i];
    Ysquared[i] = Y[i] * Y[i];
    sumX        = sumX + X[i];
    sumY        = sumY + Y[i];
    XYsum       = XYsum + XY[i];
    sumXsqrd    = sumXsqrd + Xsquared[i];
    sumYsqrd    = sumYsqrd + Ysquared[i];
  }

  numer = 1.0 * ((N * XYsum) - (sumX * sumY));
  denom = 1.0 * ((N * sumXsqrd - sumX * sumX) * (N * sumYsqrd - sumY * sumY));
  CC = numer / sqrt(denom);

  return CC ;
}