#include <stdio.h>
#include <math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>
/* #include "block.h"
 */


size_t mod (size_t a, size_t b) {
  return a  - (a/b) * b;
}




/* This function takes a data array as input, and returns a new array where
 * the original data has been avereged in blocks. If for instance the original
 * data contains 1000 elements, and blocks=10, the function will return:
 *    
 * {sum(data[i],i=0..99) / 100, sum(data[i],i=100,199)/100, ... ,sum(data[i],i=900,999)/100}
 *  
 * When dealing with data with temporal correlations this is a simple way to 
 * "make" independent data. The routine can be called by itself, but is also
 * used by 'gsl_stats_mean_std_error', 'gsl_stats_jackknife_estimate', and
 * 'gsl_stats_bootstrap_estimate'.
 */

double * gsl_block_data(const double data[], const size_t stride, const size_t size, const size_t blocks) {
  double * block_data   = NULL;
  size_t blocksize      = size / blocks;
  size_t largeblocks    = mod(size,blocks);
  size_t startindex     = 0;   
  size_t this_blocksize;
  int blocknr;
  

  if (blocks <= size) {
    block_data = (double *) calloc(blocks , sizeof(double));
    for (blocknr = 0; blocknr < blocks; blocknr++) {
      this_blocksize  = (blocknr < largeblocks) ? 1 : 0;
      this_blocksize += blocksize;
      block_data[blocknr] = gsl_stats_mean(&data[startindex] , stride, this_blocksize);
      startindex += this_blocksize;
    }
  } else GSL_ERROR_VAL("blocks > size in 'gsl_block_data'", GSL_EINVAL,0);
  return block_data;
}




/* Calculates the standard error in addition to the mean value of the
 * data in the array 'data'. If the variable 'blocks' is greater than
 * zero the data are firsted collected in blocks with the function
 * 'gsl_block_data'.  */

void gsl_stats_mean_std_error(double data[], const size_t stride, const size_t size, const size_t blocks, double *mean, double *std_error) {
  double * local_data;
  size_t   local_size   = size;
  size_t   local_stride = stride;

  if (blocks > 0) {
    local_data   = gsl_block_data(data,stride,size,blocks);
    local_size   = blocks;
    local_stride = 1;
  } else local_data = data;

  *mean      = gsl_stats_mean(local_data, local_stride, local_size);
  *std_error = gsl_stats_sd_m(local_data, local_stride, local_size, *mean) / sqrt(local_size);

  if (blocks > 0) free(local_data);
}




/* Calculates a jack-knife estimate of the estimator 'estimator' on the 
 * data in variable 'data'. The variable 'blocks' works the same way as 
 * in 'gsl_stats_mean_std_error'. The function 'estimator' must be written
 * in the calling program. 
 * 
 * A small example:
 * ----------------
 *
 * #include <stdio.h>
 * #include <math.h>
 * #include "block.h"
 * 
 * 
 * double estimator(double data[], size_t size) {
 *   int i;
 *   double estimate;
 *   estimate = 0.00;
 *   for (i=0; i<size; i++) {
 *     estimate += data[i] * log(data[i]);
 *   }
 *   return estimate / size;
 * }
 * 
 * 
 * 
 * 
 * int main (void) {
 *   double mean,error;
 *   double data[8] = {0.25, 1.75, 1.00, 1.00, 2.50, 17.50, 45.00, 10.00};
 *   gsl_stats_jackknife_estimate(data,1,8,0,estimator,&mean,&error);
 *   printf("Mean = %10.7f +- %10.7f \n",mean,error);}
 *
 *
 * This little example will calculate <x*ln(x)> with a standard error.
 */

void gsl_stats_jackknife_estimate(double data[], const size_t stride, const size_t size, const size_t blocks, double (estimator)(double *,size_t), double *exp_value, double *std_error) {
  double * local_data;
  size_t   local_size   = size;
  size_t   local_stride = stride;

  double * jack_estimate;
  double * tmp_data;
  int jacknr,index;

  if (blocks > 0) {
    local_data   = gsl_block_data(data,stride,size,blocks);
    local_size   = blocks;
    local_stride = 1;
  } else local_data = data;

  tmp_data      = (double *) calloc(local_size*2, sizeof(double));
  jack_estimate = (double *) calloc(local_size  , sizeof(double));
  
  for (index=0; index <local_size; index++) {
    tmp_data[index]              = local_data[index*local_stride];
    tmp_data[index + local_size] = tmp_data[index];
  }

  for (jacknr=0; jacknr < local_size; jacknr++) 
    jack_estimate[jacknr] = estimator(&tmp_data[jacknr],local_size-1);
					    
  *exp_value  = gsl_stats_mean(jack_estimate,1,local_size);
  *std_error  = gsl_stats_sd_with_fixed_mean(jack_estimate,1,local_size,*exp_value) * sqrt(local_size - 1);

  free(jack_estimate); free(tmp_data); 
  if (blocks > 0) free(local_data);
}




void gsl_stats_bootstrap_estimate(double data[], const size_t stride, const size_t size, const size_t blocks, const int bootstrap_replica, const unsigned long int random_seed, double (estimator)(double *,size_t), double *exp_value, double *std_error) {
  double * local_data;
  size_t   local_size   = size;
  size_t   local_stride = stride;

  double * bootstrap_estimate;
  double * tmp_data;
  int index,bootstrap_nr;

  static gsl_rng * bootstrap_rng = NULL;

  if (blocks > 0) {
    local_data   = gsl_block_data(data,stride,size,blocks);
    local_size   = blocks;
    local_stride = 1;
  } else local_data = data;

  bootstrap_estimate = (double *) calloc(bootstrap_replica, sizeof(double));
  tmp_data           = (double *) calloc(local_size       , sizeof(double));  

  if (bootstrap_rng == NULL) 
    bootstrap_rng = gsl_rng_alloc(gsl_rng_mrg);
  if (random_seed > 0)
    gsl_rng_set(bootstrap_rng, random_seed);
  
  for (bootstrap_nr=0; bootstrap_nr < bootstrap_replica; bootstrap_nr++) {
    for (index=0; index < local_size; index ++) 
      tmp_data[index] = local_data[local_stride * gsl_rng_uniform_int(bootstrap_rng, local_size)];
    bootstrap_estimate[bootstrap_nr] = estimator(tmp_data,local_size);
  }
  
  *exp_value = gsl_stats_mean(bootstrap_estimate,1,bootstrap_replica);
  *std_error = gsl_stats_sd_with_fixed_mean(bootstrap_estimate,1,bootstrap_replica,*exp_value);
  free(bootstrap_estimate); free(tmp_data);
  if (blocks > 0) free(local_data);
}