/*
   pa2.c

   Second of two implementation files for the computation of
   piece-wise analytic fits based on sample data.  Implements those
   definitions in the header file pa.h that require a random number
   generator.  The random number generator should be available as a
   global function "double rnd(void*)".

   Ramses van Zon, December 4, 2009

   Reference:
   [] Ramses van Zon and Jeremy Schofield, "Constructing smooth
      potentials of mean force, radial distribution functions and
      probability densities from sampled data", 
      arxiv.0912.0465 [cond-mat.stat-mech].
*/

#include <stdlib.h>
#include "pa.h"

extern double rnd(void* g); /* will pick up any c function called 
                               rnd that takes a pointer (of any kind) 
                               and returns a double */

/* perform a resampling of data for use with a bootstrap estimator. */ 
void paboot( int n, const double *r, void *s, double *x )
{
  int i;
  for (i = 0; i < n; ++i) 
    x[i] = r[(int)(n*rnd(s))];
}

/* perform a resampling of data for use with a block-jackknife. */
void pajackboot( int j, int m, int n, const double *r, void *s, double *x )
{
  int i, k, nb, jnb;
  nb = n/m;
  jnb = j*nb; 
  for (i = 0; i< n-nb; ++i) {
    k = (int)((n-nb)*rnd(s));
    if (k > jnb)
      k += nb;
    x[i] = r[k];
  }
}

/* perform a weighted resampling of data. */
double pabias( int n, const double *r, double (*w)(double),
               void *s, double *x ) 
{
  double nrm;
  if (w) {
    int     i, k;
    double  pmax, *p;
    p = (double*)malloc(n*sizeof(double));
    pmax = 0.0;
    nrm = 0.0;
    for (i = 0; i < n; ++i) {
      p[i] = w(r[i]);
      nrm += p[i];
      if (p[i] > pmax)
        pmax = p[i];
    }
    nrm /= n;
    for (i = 0; i < n; ++i)  
      p[i] /= pmax;
    for (i = 0; i < n; ++i) {
      do {
        k    = (int)(n*rnd(s));
        x[i] = r[k];
      } while (rnd(s) >= p[k]);
    }
    free(p);
    return nrm;
  } else {
    paboot(n, r, s, x);  
    return 1.0;
  }
}

/* perform a weighted resampling of data for use with a jackknife. */
double pajbias( int j, int m, int n, const double *r,
                double (*w)(double), void *s, double *x ) 
{
  double nrm;
  if (w) {
    int     i, k, nb,jnb;
    double  pmax, *p;
    nb   = n/m;
    jnb  = j*nb;
    p    = (double*)malloc(n*sizeof(double));
    pmax = 0.0;
    nrm  = 0.0;
    for ( i = 0; i < n; ++i) {
      p[i] = w(r[i]);
      nrm += p[i];
      if (p[i] > pmax)
        pmax = p[i];
    }
    nrm /= n;
    for ( i = 0; i < n; ++i)
      p[i] /= pmax;
    for (i = 0; i < n-nb; ++i) {
      do {
        k = (int)((n-nb)*rnd(s));
        if (k > jnb)
          k += nb;
        x[i] = r[k];
      } while (rnd(s) >= p[k]);
    }
    free(p);
    return nrm;
  } else {
    pajackboot(j,m,n,r,s,x);
    return 1.0;
  }
}