/******************************************************************/
/*
 * sampling.c : some useful routines for sampling things
 */
/******************************************************************/

#include <stdio.h>

int SampleDiscrete(float probabilities[], float total, double *x_1, double *pdf)
{
  int i=0;
  float sum, left;
  double sample = *x_1 * total;

  sum = probabilities[0];

  while(sample > sum) 
  {
    i++;
    sum += probabilities[i];
  }

  /* Rescale x_1 */
  left = sum - probabilities[i];
  *x_1 = ((sample - left) / (sum - left));
  *pdf = probabilities[i] / total;

  return i;
}

int DSampleDiscrete(double probabilities[], double total, double *x_1, double *pdf)
{
  int i=0;
  double sum, left;
  double sample = *x_1 * total;

  sum = probabilities[0];

  /*  printf("D Total %g, p 0 %g\n", total, sum); */

  while(sample > sum) 
  {
    i++;
    sum += probabilities[i];
    /* printf("D Sum %g, p %i %g\n", sum, i, probabilities[i]); */
  }

  /* Rescale x_1 */
  left = sum - probabilities[i];
  /* printf("D Left %g, Sample %g\n", left, sample); */

  *x_1 = ((sample - left) / (sum - left));
  *pdf = probabilities[i] / total;
  return i;
}