/* 
 * Selecting an outcome according to a given probability distribution 
 * is a common task faced by simulation programs. The following program 
 * demonstrates one way of accomplishing this task. In the main program, 
 * two random experiments are simulated by calling the function 
 * "distro_select()", defined at the bottom, with appropriate arguments.
 */


#include <stdio.h>
#include <stdlib.h>

int distro_select(int numOfOutcomes, double *probDistro);

int main(){
	
	int i, n;

	// Tossing a coin results in 2 random outcomes: heads or tail
	// Probability distribution of a fair coin
	double faircoin_probdistro[] = {
		1.0 / 2.0, 
		1.0 / 2.0
	};

	// Tossing a die results in 6 different outcomes
	// Probability distribution of a standard die
	double dice_probdistro[] = {
		1.0 / 6.0, 
		1.0 / 6.0, 
		1.0 / 6.0, 
		1.0 / 6.0, 
		1.0 / 6.0, 
		1.0 / 6.0 
	};

	//
	// Provide a seed to the random generator
	// Type "man rand()" on the command prompt to see more details
	//
	srand(time(NULL));

	//
	// Now do some simulations
	//
	printf("\n");
	printf("Tossing a fair coin 100 times: \n");
	for(i = 0; i < 100; i++){
		n = distro_select(2, faircoin_probdistro);
		printf("%d ", n);
	}
	printf("\n");

	printf("\n");
	printf("Throwing a fair die 100 times: \n");
	for(i = 0; i < 100; i++){
		n = distro_select(6, dice_probdistro);
		printf("%d ", n);
	}
	printf("\n\n");

	return 0;
}

/*
  Select an outcome according to a given probability distribution
  of a discrete random variable X with a finite set of outcomes
  x0, x1, ..., xn. The probability of each outcome is represented
  by the array probDistro[0], probDistro[1], ..., probDistro[n].

   numOfOutcomes: number of possible outcomes
   probDistro: pointer to the first element in the probability distribution

   Return a possitive index i between 0 and numOfOutcomes representing the
   i-th outcome if succeeded. If this function fails it returns -1. 
*/  
int distro_select(int numOfOutcomes, double *probDistro){
	
	int i;
	double range[numOfOutcomes];
	double x;

	//
	// Construct intervals representing range of each outcome.
        // Because the sum over i of probDistro[i] equals 1, each interval  
	// uniquely represents an outcome. For example, the outcome
	// x0 is represented by (0,probDistro[0]), the next outcome is
	// represented by (probDistro[0],probDistro[0]+probDistro[1]), etc.
	//
	range[0] = probDistro[0];
	for(i = 1; i < numOfOutcomes; i++){
		range[i] = range[i-1] + probDistro[i];
	}

	//
	// Select an outcome using the random generator
	//
	x = (double)rand() / (double)RAND_MAX;
	for(i = 0; i < numOfOutcomes; i++){
		if (x < range[i]){
			return i;
		}
	}
	
	return -1;		

}