/******************************************************************************
* PDPCTOOLS.C                                                                 *
*                                                                             *
* Various helpful (and a few necessary) funtions for earthquake analysis.     *
*                                                                             *
* James Holliday                                                              *
* University of California - Davis                                            *
*                                                                             *
******************************************************************************/


/* Standard include files */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

/* Custom include files */
#include "pdpctools.h"
#include "fftw3.h"


/******************************************************************************
* EQERROR                                                                     *
*                                                                             *
*      This function acts as a standard error reporting mechanism.            *
*      (Standard in that it will be identical in all these routines.)         *
*                                                                             *
******************************************************************************/
void inline eqerror(char* message, char* routine)
{
  fprintf(stderr,"*** %s in %s().   Exiting!\n",message,routine);
  exit(-1);
}


/******************************************************************************
* DAYS_FROM_DATE                                                              *
*                                                                             *
*      This function takes in a two DATE objects and calculates how many      *
*      days have passed between them.                                         *
*                                                                             *
******************************************************************************/
int days_from_date(DATE start, DATE end)
{
  int sYear = 0, sMonth = 0, sDay = 0;
  int eYear = 0, eMonth = 0, eDay = 0;
  int days = 0;
  int daysinmonth[13] = {0,31,28,31,30,31,30,31,31,30,31,30,31};
  int i = 0;

  /* Break the DATE structure into year/month/day objects */
  sYear  = (int)floor(start / 10000.0);
  sMonth = (int)floor(start /   100.0) - (100.0 * floor(start / 10000.0) );
  sDay   = (int)floor(start /     1.0) - (100.0 * floor(start /   100.0) );

  eYear  = (int)floor(end / 10000.0);
  eMonth = (int)floor(end /   100.0) - (100.0 * floor(end / 10000.0) );
  eDay   = (int)floor(end /     1.0) - (100.0 * floor(end /   100.0) );

  /* Check that start comes before end */
  if ( end < start )  eqerror("End date predates start date","DAYS_FROM_DATE");

  /* Calculate the number of days from Jan 1, SYEAR to Dec 31, EYEAR */
  /* Treat the years as normal year--we'll handle leap years next. */
  days = (eYear - sYear + 1)*365;

  /* Tack on an extra day for every leap year in the interval. */
  for (i=sYear; i<=eYear; i++)
  {
    if ( (i%4 == 0) && !( (i%100 == 0) && (i%400 != 0) ) )  days++;
  }

  /* Subtract off days from Jan 1 to SMONTH SDAY.  Check for leap year. */
  for (i=1; i<sMonth; i++)  days -= daysinmonth[i];
  for (i=1; i<sDay  ; i++)  days--;

  if ( ( (sMonth < 2) || ( (sMonth == 2) && (sDay == 29) ) ) &&
       (sYear%4 == 0) && !((sYear%100 == 0) && (sYear%400 != 0)) ) days--;

  /* Subtract off days from EMONTH EDAY to Dec 31.  Check for leap year. */
  for (i=eMonth+1; i<=12;                  i++)  days -= daysinmonth[i];
  for (i=eDay  +1; i<=daysinmonth[eMonth]; i++)  days--;

  /* Don't double count Feb 29 */
  if ( (eMonth <= 2) && !( (eMonth == 2) && (eDay == 29) ) &&
       (eYear%4 == 0) && !((eYear%100 == 0) && (eYear%400 != 0)) ) days--;

  /* Return our result. */

  return days;
}


/******************************************************************************
* BIN_DATA                                                                    *
*                                                                             *
*      This function takes in a catalog, bin matrix, and parameter struct     *
*      and fills the matrix--thus creating n(X_i;T).  The total number        *
*      of binned quakes is returned.                                          *
*                                                                             *
*      The catalog is expected to be in the following format (6 items/line):  *
*      "year"  "month"  "day"  "lon"  "lat"  "mag".  Also, events are assumed *
*      to be time sorted.  Wrong formats may cause unpredictable behavior.    *
*                                                                             *
******************************************************************************/
int bin_data(FILE *cat, int **bins, eqparams p)
{
  /* Catalog variables */
  int    year  = 0;
  int    month = 0;
  int    day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Number of quakes within our window */
  int c = 0;

  /* Bin location */
  int t = 0;
  int x = 0;
  int y = 0;

  DATE date = 0;

  /* Start reading from the catalog file and parsing the data. */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
                 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    /* Create DATE object */
    date = (year*10000) + (month*100) + (day*1);

    /* If the quake date is after our stop date, return early  */
    if ( date > p.D2 ) return c;

    /* If the quake falls in our parameter window, bin it! */
    if ( (lon >= p.lon_min) && (lon  < p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  < p.lat_max ) &&
	 (mag >= p.min_mag) && (date >= p.D0      ) )
    {
      /* Increment the bin number */
      t = (int)floor(days_from_date(p.D0,date)/p.timestep);
      x = (int)floor((lon-p.lon_min)/p.binsize);
      y = (int)floor((lat-p.lat_min)/p.binsize);

      bins[t][y*p.numLon + x] += 1.0;

      /* Increment the count */
      c++;
    }
  }

  /* Return the number of valid quakes */
  return c;
}


int bin_data_f(FILE *cat, float **bins, eqparams p)
{
  /* Catalog variables */
  int    year  = 0;
  int    month = 0;
  int    day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Number of quakes within our window */
  int c = 0;

  /* Bin location */
  int t = 0;
  int x = 0;
  int y = 0;

  DATE date = 0;

  /* Start reading from the catalog file and parsing the data. */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
                 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    /* Create DATE object */
    date = (year*10000) + (month*100) + (day*1);
    /* If the quake date is after our stop date, return early  */
    if ( date > p.D2 ) return c;

    /* If the quake falls in our parameter window, bin it! */
    if ( (lon >= p.lon_min) && (lon  < p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  < p.lat_max ) &&
	 (mag >= p.min_mag) && (date >= p.D0      ) )
    {
      /* Increment the bin number */
      t = (int)floor(days_from_date(p.D0,date)/p.timestep);
      x = (int)floor((lon-p.lon_min)/p.binsize);
      y = (int)floor((lat-p.lat_min)/p.binsize);

      bins[y*p.numLon + x][t] += 1.0 / (p.timestep * p.binsize * p.binsize);

      /* Increment the count */
      c++;
    }
  }

  /* Return the number of valid quakes */
  return c;
}

int bin_data_c(FILE *cat, cdouble **bins, eqparams p)
{
  /* Catalog variables */
  int    year  = 0;
  int    month = 0;
  int    day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Number of quakes within our window */
  int c = 0;

  /* Bin location */
  int t = 0;
  int x = 0;
  int y = 0;

  DATE date = 0;

  /* Start reading from the catalog file and parsing the data. */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
                 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    /* Create DATE object */
    date = (year*10000) + (month*100) + (day*1);

    /* If the quake date is after our stop date, return early  */
    if ( date > p.D2 ) return c;

    /* If the quake falls in our parameter window, bin it! */
    if ( (lon >= p.lon_min) && (lon  < p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  < p.lat_max ) &&
	 (mag >= p.min_mag) && (date >= p.D0      ) )
    {
      /* Increment the bin number */
      t = (int)floor(days_from_date(p.D0,date)/p.timestep);
      x = (int)floor((lon-p.lon_min)/p.binsize);
      y = (int)floor((lat-p.lat_min)/p.binsize);

      bins[y*p.numLon + x][t] += 1.0 / (p.timestep * p.binsize * p.binsize);

      /* Increment the count */
      c++;
    }
  }

  /* Return the number of valid quakes */
  return c;
}

/******************************************************************************
* HILBERT_TRANSFORM                                                           *
*                                                                             *
*      This function takes in an array (N sets) of vectors (length L) and     *
*      performs a hilbert transform.  The vectors are assumed to be purely    *
*      real (stored as complex doubles) and the hilbert transorm component    *
*      is returned in the imaginary component.                                *
*                                                                             *
*      This function requires FFTW be installed.                              *
*                                                                             *
******************************************************************************/
void hilbert_transform(cdouble **vector, int N, int L)
{
  /* Calculation variables */
  cdouble old  = 0.0;
  cdouble prev = 0.0;
  cdouble next = 0.0;

  /* FFTW variables */
  fftw_plan forward;
  fftw_plan reverse;
  fftw_complex* tem = NULL;
  fftw_complex* out = NULL;

  /* Iteration variables */
  int i,j;

  /* Allocate space for FFT calculations and create plan files */
  tem = fftw_malloc(sizeof(fftw_complex)*N*L);
  out = fftw_malloc(sizeof(fftw_complex)*N*L);

  forward = fftw_plan_dft_2d(N,L,vector[0],tem,FFTW_FORWARD,FFTW_ESTIMATE);
  reverse = fftw_plan_dft_2d(N,L,tem,out,FFTW_BACKWARD,FFTW_ESTIMATE);

  /* Smooth the input vector before applying FFT */
  for (i=0; i<N; i++)
  {
    for (j=0; j<L; j++)
    {
      if (j == 0)  prev = vector[i][j];
      else         prev = old;
  
      old  = vector[i][j];
  
      if (j == L-1)  next = vector[i][j];
      else           next = vector[i][j+1];
  
      vector[i][j] = (0.25*prev) + (0.50*old) + (0.25*next);
    }
  }

  /* Execute FFT transform */
  fftw_execute( forward );
  for (j=0; j<L*N; j++)
  {
    if ((j%L) <  L/2) tem[j] *= I/(L*N);
    if ((j%L) >= L/2) tem[j] /= I*(L*N);
  }
  fftw_execute( reverse );

  /* Add the transform to the original vector */
  for (i=0; i<N; i++)
    for (j=0; j<L; j++)
      vector[i][j] += I * creal(out[i*L + j]);

  /* Delete plan files and free up memory */
  fftw_destroy_plan(reverse);
  fftw_destroy_plan(forward);
  fftw_free(out);
  fftw_free(tem);

  /* That's all, folks */
  return;
}


/******************************************************************************
* SAVE_BIG_EVENTS                                                             *
*                                                                             *
*      This function takes in a catalog, date integers, and minimum           *
*      magnitude and saves to file events with magnitude greater than M.      *
*                                                                             *
******************************************************************************/
int save_big_events(FILE *catalog, char* filename, eqparams p, double M)
{
  /* Number of post T2 events */
  int count = 0;

  /* Catalog variables */
  int   year  = 0;
  int   month = 0;
  int   day   = 0;
  DATE  date  = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Output file */
  FILE* file = NULL;

  /* Open the file for output */
  file = fopen(filename,"w");
  if ( file == NULL )  eqerror("Could not open output file","SAVE_BIG_EVENTS");

  /* Start reading from the catalog file and parsing the data. */
  rewind(catalog);
  while ( fscanf(catalog,"%i %i %i %lf %lf %lf",
                 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    /* Create date integer */
    date = (year * 10000) + (month * 100) + (day * 1);

    /* If the quake falls in our parameter window, save it! */
    if ( (lon >= p.lon_min) && (lon  <= p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  <= p.lat_max ) &&
	 (mag >= M        ) && (date >= p.D1      ) &&
	 (date < p.D3     ) )
    {

      if(1)
      {
	/* D1 <= date <= D2 */
	if ( date <= p.D2 )  fprintf(file,"%f %f %f i\n",lon,lat,mag);

	/* D2 < date */
	if ( date >  p.D2 )
	{
	  fprintf(file,"%f %f %f c\n",lon,lat,mag);
	  count++;
	}
      }
      else
      {
	if (date <= p.D2-10000)
	  fprintf(file,"%f %f %f i\n",lon,lat,mag);
	if (date>p.D2-10000 && date<p.D2)
	  fprintf(file,"%f %f %f x\n",lon,lat,mag);
	if (date>p.D2 && date<p.D2+10000)
	  fprintf(file,"%f %f %f c\n",lon,lat,mag);
      }
    }
  }

  /* Get outta Dodge */
  return count;
}


/******************************************************************************
* NORMALIZE                                                                   *
*                                                                             *
*      This function takes in a vector of doubles and normalizes it           *
*      (subtracts the mean and divides by total length).                      *
*                                                                             *
******************************************************************************/
double normalize(double *vector, int num)
{
  /* Storage and iterating variables */
  double sum  = 0.0;
  double ave  = 0.0;
  double norm = 0.0;
  int i = 0;

  /* Loop thru once to get the mean */
  for (i=0; i<num; i++)  sum += vector[i];
  ave = sum / num;

  /* Loop thru again to remove the mean and get the normalization */
  sum = 0.0;
  for (i=0; i<num; i++)
  {
    vector[i] = vector[i] - ave;
    sum += vector[i] * vector[i];
  }

  /* Check for ZERO-vector */
  if ( sum )  norm = sqrt(sum);
  else        norm = 1.0;

  /* Loop thru one last time and normalize the vector */
  for (i=0; i<num; i++)  vector[i] = vector[i] / norm;

  /* That's all, folks */
  return norm;
}


double normalize_f(float *vector, int num)
{
  /* Storage and iterating variables */
  double sum  = 0.0;
  double ave  = 0.0;
  double norm = 0.0;
  int i = 0;

  /* Loop thru once to get the mean */
  for (i=0; i<num; i++)  sum += vector[i];
  ave = sum / num;

  /* Loop thru again to remove the mean and get the normalization */
  sum = 0.0;
  for (i=0; i<num; i++)
  {
    vector[i] = vector[i] - ave;
    sum += vector[i] * vector[i];
  }

  /* Check for ZERO-vector */
  if ( sum )  norm = sqrt(sum);
  else        norm = 1.0;

  /* Loop thru one last time and normalize the vector */
  for (i=0; i<num; i++)  vector[i] = vector[i] / norm;

  /* That's all, folks */
  return norm;
}

double normalize_c(cdouble *vector, int num)
{
  /* Storage and iterating variables */
  cdouble sum = 0.0;
  cdouble ave = 0.0;
  double norm = 0.0;
  int i = 0;

  /* Loop thru once to get the mean */
  for (i=0; i<num; i++)  sum += vector[i];
  ave = sum / num;

  /* Loop thru again to remove the mean and get the normalization */
  sum = 0.0;
  for (i=0; i<num; i++)
  {
    vector[i] = vector[i] - ave;
    sum += creal( vector[i] * conj(vector[i]) );
  }

  /* Check for ZERO-vector */
  if ( sum )  norm = csqrt(sum);
  else        norm = 1.0;

  /* Loop thru one last time and normalize the vector */
  for (i=0; i<num; i++)  vector[i] = vector[i] / norm;

  /* That's all, folks */
  return norm;
}


/******************************************************************************
* CENTER                                                                      *
*                                                                             *
*      This function takes in a vector of doubles and centers it              *
*      by removing the mean.                                                  *
*                                                                             *
******************************************************************************/
void center(double *vector, int num)
{
  /* Storage and iterating variables */
  double sum = 0.0;
  int i = 0;

  /* Loop thru once to get the mean */
  for (i=0; i<num; i++)  sum += vector[i];

  /* Loop thru again to remove the mean */
  for (i=0; i<num; i++)  vector[i] = vector[i] - (sum / num);

  /* That's all, folks */
  return;
}

void center_c(cdouble *vector, int num)
{
  /* Storage and iterating variables */
  cdouble sum = 0.0;
  int i = 0;

  /* Loop thru once to get the mean */
  for (i=0; i<num; i++)  sum += vector[i];

  /* Loop thru again to remove the mean */
  for (i=0; i<num; i++)  vector[i] = vector[i] - (sum / num);

  /* That's all, folks */
  return;
}

/******************************************************************************
* SAVE_GMT                                                                    *
*                                                                             *
*      This function takes in a spatial vector of doubles and saves a three-  *
*      column data file ready for GMT.  The data column is scaled so that     *
*      the values all lie between -1 and 1                                    *
*                                                                             *
******************************************************************************/
void save_gmt(char* filename, double *vector, eqparams params)
{
  /* Storage and iterating variables */
  FILE* file = NULL;
  double x = 0.0;
  double y = 0.0;
  int i = 0;
  int j = 0;
  int N = 0;

  /* Open the file for output */
  file = fopen(filename,"w");
  if ( file == NULL )  eqerror("Could not open output file","SAVE_GMT");

  if(0)
  {
    int ii;
    double ave = 0.0;
    for (ii=0; ii<params.numBin; ii++)
      if (vector[ii] != 0.0) { ave += vector[ii]; N++; }
    ave /= N;
    for (j=0; j<params.numLat; j++)
      for (i=0; i<params.numLon; i++)
      {
	x = params.lon_min + i*params.binsize + 0.5*params.binsize;
	y = params.lat_min + j*params.binsize + 0.5*params.binsize;

	if (vector[j*params.numLon + i] != 0.0)
	  fprintf(file,"%f %f %.15f\n",
		  x,y,(vector[j*params.numLon + i])-ave)/(2*ave*params.numBin);
	else
	  fprintf(file,"%f %f %.15f\n",x,y,(vector[j*params.numLon + i]));
      }

    fclose(file);
    return;
  }

  /* Loop thru the data and save the output */
  for (j=0; j<params.numLat; j++)
    for (i=0; i<params.numLon; i++)
    {
      x = params.lon_min + i*params.binsize + 0.5*params.binsize;
      y = params.lat_min + j*params.binsize + 0.5*params.binsize;
      if(1)
      {
	fprintf(file,"%f %f %.15f\n",x,y,vector[j*params.numLon + i]);
      }
      else
      {
	double value = 0.0;
	int ii,jj;
	for (ii = MAX(0,i-1); ii < MIN(params.numLon,i+2); ii++) {
	  for (jj = MAX(0,j-1); jj < MIN(params.numLat,j+2); jj++) {
	    //value = MAX(value,vector[jj*params.numLon + ii]);
	    value += MAX(0,vector[jj*params.numLon + ii])/9.0;
	  }
	}
	fprintf(file,"%f %f %.15f\n",x,y,value);
      }
    }

  /* Close the file and get outta Dodge */
  fclose(file);
  return;
}

void ___save_gmt(char* filename, double *vector, eqparams params)
{
  /* Storage and iterating variables */
  FILE* file = NULL;
  double x = 0.0;
  double y = 0.0;
  int i = 0;
  int j = 0;

  /* Open the file for output */
  file = fopen(filename,"w");
  if ( file == NULL )  eqerror("Could not open output file","SAVE_GMT");

  if(0)
  {
    int ii;
    double ave = 0.0;
    for (ii=0; ii<params.numBin; ii++)  ave += vector[ii];
    ave /= params.numBin;
    for (j=0; j<params.numLat; j++)
      for (i=0; i<params.numLon; i++)
      {
	x = params.lon_min + i*params.binsize + 0.5*params.binsize;
	y = params.lat_min + j*params.binsize + 0.5*params.binsize;

	fprintf(file,"%f %f %.15f\n",
		x,y,(vector[j*params.numLon + i])+ave)/(2*ave*params.numBin);
      }

    fclose(file);
    return;
  }

  /* Loop thru the data and save the output */
  for (j=0; j<params.numLat; j++)
    for (i=0; i<params.numLon; i++)
    {
      x = params.lon_min + i*params.binsize + 0.5*params.binsize;
      y = params.lat_min + j*params.binsize + 0.5*params.binsize;
      if(1)
      {
	fprintf(file,"%f %f %.15f\n",x,y,vector[j*params.numLon + i]);
      }
      else
      {
	double value = 0.0;
	int ii,jj;
	for (ii = MAX(0,i-1); ii < MIN(params.numLon,i+2); ii++) {
	  for (jj = MAX(0,j-1); jj < MIN(params.numLat,j+2); jj++) {
	    //value = MAX(value,vector[jj*params.numLon + ii]);
	    value += MAX(0,vector[jj*params.numLon + ii])/9.0;
	  }
	}
	fprintf(file,"%f %f %.15f\n",x,y,value);
      }
    }

  /* Close the file and get outta Dodge */
  fclose(file);
  return;
}


/******************************************************************************
* CALC_LIKELIHOOD                                                             *
*                                                                             *
*      These functions take in a forcast vector and calculates the (log)      *
*      likelihood value.                                                      *
*                                                                             *
******************************************************************************/
int fact(int N) { return (N <= 1) ? 1 : N*fact(N-1); }
double LOG(double x) { return (x>0) ? log(x) : -9999999; }

double calc_likelihood_relm(double *S, eqparams p, FILE *cat)
{
  /* Calculation Variables */
  double log_L  = 0.0;
  double log_L0 = 0.0;
  int   N      = 0;
  double pTot   = 0.0;

  /* Catalog variables */
  DATE  date  = 0;
  int   year  = 0;
  int   month = 0;
  int   day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Storage Variables */
  double* probs  = NULL;
  int*   actual = NULL;

  /* Iterating Variables */
  int i = 0;
  int j = 0;

  /* Allocate and initialize memory */
  probs  = (double*)calloc(p.numBin , sizeof(double));
  actual = (   int*)calloc(p.numBin , sizeof(   int));

  if ( !probs  )  eqerror("Not enough memory for probs","CALC_LIKELIHOOD");
  if ( !actual )  eqerror("Not enough memory for actuals","CALC_LIKELIHOOD");

  for (i=0; i<p.numBin; i++)
  {
    probs[i]  = MIN( S[i] , 1.0 );
    actual[i] = 0;

    pTot += probs[i];
  }

  /* Scan the catalog and count the number of post-T1/T2 events */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
		 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    if ( (lon >= p.lon_min) && (lon  <= p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  <= p.lat_max ) &&
	 (mag >= 5.0      )                       )
    {
      i = (int)floor((lon-p.lon_min)/p.binsize);
      j = (int)floor((lat-p.lat_min)/p.binsize);

      date = (year*10000) + (month*100) + (day*1);

      /* Find "actuals" */
      if ( (date > p.D2) && (date <= p.D3) )
      {
	actual[i + j*p.numLon] += 1;
	N++;
      }
    }
  }

  /* Scale the calculated probabilities so that they integrate to N_events. */
  /* Calculate the log likelihood using the RELMTEST method.                */
  for (i=0; i<p.numBin; i++)
  {
    probs[i] *= N / pTot;

    log_L += actual[i]*LOG(probs[i]) - log(fact(actual[i])) - probs[i];

    if ( actual[i] != 0 )
      log_L0 += actual[i]*log(actual[i]) - log(fact(actual[i])) - actual[i];
  }

  /* Free up memory and return our result */
  free(probs);
  free(actual);

  return (log_L - log_L0);
}

double calc_likelihood_relm_trunc(double *S, eqparams p, FILE *cat)
{
  /* Calculation Variables */
  double log_L  = 0.0;
  double log_L0 = 0.0;
  int   N      = 0;
  double pTot   = 0.0;

  /* Catalog variables */
  DATE  date  = 0;
  int   year  = 0;
  int   month = 0;
  int   day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Storage Variables */
  double* probs  = NULL;
  int*   actual = NULL;

  /* Iterating Variables */
  int i = 0;
  int j = 0;

  /* Allocate and initialize memory */
  probs  = (double*)calloc(p.numBin , sizeof(double));
  actual = (   int*)calloc(p.numBin , sizeof(   int));

  if ( !probs  )  eqerror("Not enough memory for probs","CALC_LIKELIHOOD");
  if ( !actual )  eqerror("Not enough memory for actuals","CALC_LIKELIHOOD");

  for (i=0; i<p.numBin; i++)
  {
    /* Originally we normalized to the standard deviation.           */
    /* The "historgram cut" plays a big role in this case.           */
    /* We can scale up by multiplying by N and put it back in play.  */
    /* (N vs sqrt(N) because the probability has been squared.)      */
    probs[i] = MIN( S[i]*p.numBin , 1.0      ); /* Probs = [0 , 1.0] */

    pTot += probs[i];
  }

  /* Scan the catalog and count the number of post-T1/T2 events */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
		 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    if ( (lon >= p.lon_min) && (lon  <= p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  <= p.lat_max ) &&
	 (mag >= 5.0      )                       )
    {
      i = (int)floor((lon-p.lon_min)/p.binsize);
      j = (int)floor((lat-p.lat_min)/p.binsize);

      date = (year*10000) + (month*100) + (day*1);

      /* Find "actuals" */
      if ( (date > p.D2) && (date <= p.D3) )
      {
	actual[i + j*p.numLon] += 1;
	N++;
      }
    }
  }

  /* Scale the calculated probabilities so that they integrate to N_events. */
  /* Calculate the log likelihood using the RELMTEST method.                */
  for (i=0; i<p.numBin; i++)
  {
    probs[i] *= N / pTot;

    log_L += actual[i]*LOG(probs[i]) - log(fact(actual[i])) - probs[i];

    if ( actual[i] != 0 )
      log_L0 += actual[i]*log(actual[i]) - log(fact(actual[i])) - actual[i];
  }

  /* Free up memory and return our result */
  free(probs);
  free(actual);

  return (log_L - log_L0);
}

double calc_likelihood_kristy(double *S, eqparams p, FILE *cat)
{
  /* Calculation Variables */
  double L = 0.0, log_L = 0.0;
  double Ptot = 0.0;
  double R = 0.0;
  double sigma = p.binsize * p.binsize;
  double pLat = 0.0, pLon = 0.0;
  int flagged = 0;

  /* Catalog variables */
  DATE  date  = 0;
  int   year  = 0;
  int   month = 0;
  int   day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Bin location */
  int x = 0;
  int y = 0;
  int z = 0;

  /* Storage and iterating variables */
  double *probs = NULL;

  double alon[100] = {0.0};
  double alat[100] = {0.0};
  int   nactual   = 0;

  double elon[100] = {0.0};
  double elat[100] = {0.0};
  int   nevents   = 0;

  int i = 0;
  int j = 0;

  /* Allocate and initialize memory */
  probs = (double*)calloc(p.numBin , sizeof(double));
  if ( !probs )  eqerror("Not enough memory for probs","CALC_LIKELIHOOD");

  for (i=0; i<p.numBin; i++)
  {
    probs[i] = MIN( S[i] , 1.0 );  /* Probs = [0 , 1.0] */
    Ptot += probs[i];
  }

  /* Scan the catalog and count the number of post-T1/T2 events */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
		 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    if ( (lon >= p.lon_min) && (lon  <= p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  <= p.lat_max ) &&
	 (mag >= 5.0      )                       )
    {
      x = (int)floor((lon-p.lon_min)/p.binsize);
      y = (int)floor((lat-p.lat_min)/p.binsize);

      date = (year*10000) + (month*100) + (day*1);

      /* Find "actuals" */
      if ( (date > p.D1) && (date < p.D2) )
      {
	alon[nactual] = lon;
	alat[nactual] = lat;
	nactual++;
      }

      /* Find "actuals" */
      if ( (date > p.D2) && (date <= p.D3) )
      {
	elon[nevents] = lon;
	elat[nevents] = lat;
	nevents++;
      }
    }
  }

  /* Remove actuals from data */
  for (i=0; i<nactual; i++)
  {
    for (x=0; x<p.numLon; x++)
    for (y=0; y<p.numLat; y++)
    {
      z = y*p.numLon + x;

      pLon = p.lon_min + x*p.binsize + 0.5*p.binsize;
      pLat = p.lat_min + y*p.binsize + 0.5*p.binsize;

      R = pow(alat[i]-pLat,2.0) + pow(alon[i]-pLon,2.0)*pow(cos(alat[i]),2.0);

      if (R < (4.0 * sigma) )
      {
	Ptot -= probs[z];
	probs[z] = 0.0;
      }
    }
  }

  /* Calculate LOG(L) */
  for (i=0; i<nevents; i++)
  {
    L = 0.0;
    flagged = 0;

    for (j=0; j<nactual; j++)
    {
      R = pow(alat[j]-elat[i],2.)+pow(alon[j]-elon[i],2.)*pow(cos(alat[j]),2.);

      if (R < (4.0 * sigma) )  flagged = 1;
    }

    if ( !flagged )
    {
      for (x=0; x<p.numLon; x++)
      for (y=0; y<p.numLat; y++)
      {
	z = y*p.numLon + x;

	pLon = p.lon_min + x*p.binsize + 0.5*p.binsize;
	pLat = p.lat_min + y*p.binsize + 0.5*p.binsize;

	R = pow(pLat-elat[i],2.0) + pow(pLon-elon[i],2.0)*pow(cos(pLat),2.0);
	L += probs[z]/Ptot * (1.0/sigma) * exp(-R/sigma);
      }
    }

    if (L > 0.0)  log_L += log(L);
  }

  /* Free up memory and return our result */
  free(probs);

  return log_L;
}

double calc_likelihood_kristy_trunc(double *S, eqparams p, FILE *cat)
{
  /* Calculation Variables */
  double L = 0.0, log_L = 0.0;
  double Ptot = 0.0;
  double R = 0.0;
  double sigma = p.binsize * p.binsize;
  double pLat = 0.0, pLon = 0.0;
  int flagged = 0;

  /* Catalog variables */
  DATE  date  = 0;
  int   year  = 0;
  int   month = 0;
  int   day   = 0;
  double lon   = 0.0;
  double lat   = 0.0;
  double mag   = 0.0;

  /* Bin location */
  int x = 0;
  int y = 0;
  int z = 0;

  /* Storage and iterating variables */
  double *probs = NULL;

  double alon[100] = {0.0};
  double alat[100] = {0.0};
  int   nactual   = 0;

  double elon[100] = {0.0};
  double elat[100] = {0.0};
  int   nevents   = 0;

  int i = 0;
  int j = 0;

  /* Allocate and initialize memory */
  probs = (double*)calloc(p.numBin , sizeof(double));
  if ( !probs )  eqerror("Not enough memory for probs","CALC_LIKELIHOOD");

  for (i=0; i<p.numBin; i++)
  {
    /* Originally we normalized to the standard deviation.           */
    /* The "historgram cut" plays a big role in this case.           */
    /* We can scale up by multiplying by N and put it back in play.  */
    /* (N vs sqrt(N) because the probability has been squared.)      */
    probs[i] = MIN( S[i]*p.numBin , 1.0 ); /* Probs = [0 , 1.0]      */
    Ptot += probs[i];
  }

  /* Scan the catalog and count the number of post-T1/T2 events */
  rewind(cat);
  while ( fscanf(cat,"%i %i %i %lf %lf %lf",
		 &year,&month,&day,&lon,&lat,&mag) != EOF )
  {
    if ( (lon >= p.lon_min) && (lon  <= p.lon_max ) &&
	 (lat >= p.lat_min) && (lat  <= p.lat_max ) &&
	 (mag >= 5.0      )                       )
    {
      x = (int)floor((lon-p.lon_min)/p.binsize);
      y = (int)floor((lat-p.lat_min)/p.binsize);

      date = (year*10000) + (month*100) + (day*1);

      /* Find "actuals" */
      if ( (date > p.D1) && (date < p.D2) )
      {
	alon[nactual] = lon;
	alat[nactual] = lat;
	nactual++;
      }

      /* Find "actuals" */
      if ( (date > p.D2) && (date <= p.D3) )
      {
	elon[nevents] = lon;
	elat[nevents] = lat;
	nevents++;
      }
    }
  }

  /* Remove actuals from data */
  for (i=0; i<nactual; i++)
  {
    for (x=0; x<p.numLon; x++)
    for (y=0; y<p.numLat; y++)
    {
      z = y*p.numLon + x;

      pLon = p.lon_min + x*p.binsize + 0.5*p.binsize;
      pLat = p.lat_min + y*p.binsize + 0.5*p.binsize;

      R = pow(alat[i]-pLat,2.0) + pow(alon[i]-pLon,2.0)*pow(cos(alat[i]),2.0);

      if (R < (4.0 * sigma) )
      {
	Ptot -= probs[z];
	probs[z] = 0.0;
      }
    }
  }

  /* Calculate LOG(L) */
  for (i=0; i<nevents; i++)
  {
    L = 0.0;
    flagged = 0;

    for (j=0; j<nactual; j++)
    {
      R = pow(alat[j]-elat[i],2.)+pow(alon[j]-elon[i],2.)*pow(cos(alat[j]),2.);

      if (R < (4.0 * sigma) )  flagged = 1;
    }

    if ( !flagged )
    {
      for (x=0; x<p.numLon; x++)
      for (y=0; y<p.numLat; y++)
      {
	z = y*p.numLon + x;

	pLon = p.lon_min + x*p.binsize + 0.5*p.binsize;
	pLat = p.lat_min + y*p.binsize + 0.5*p.binsize;

	R = pow(pLat-elat[i],2.0) + pow(pLon-elon[i],2.0)*pow(cos(pLat),2.0);
	L += probs[z]/Ptot * (1.0/sigma) * exp(-R/sigma);
      }
    }

    if (L > 0.0)  log_L += log(L);
  }

  /* Free up memory and return our result */
  free(probs);

  return log_L;
}