//-------------------------------------------------------------------
//	hemoxfer.cpp
//
//	This program utilizes research on the relative frequency of
//	ABO blood types in the United States to compute an estimate 
//	for the liklihood of a donor-recipient mismatch between two
//	randomly selected individuals, using a computer simulation. 
//
//		to compile:  $ g++ hemoxfer.cpp -o hemoxfer
//		to execute:  $ ./hemoxfer
//
//	Reference: 
//	Online research statistics: "Blood Types in the U.S.," 
//	Stanford University School of Medicine (2011).
//	
//	programmer: ALLAN CRUSE
//	written on: 03 FEB 2011
//-------------------------------------------------------------------

#include <stdio.h>	// for printf() 
#include <stdlib.h>	// for drand48(), srand48()
#include <time.h>	// for time()

#define MAX_REPS	1000000		// one million repetitions

#define TYPE_O		(0.374 + 0.066)
#define TYPE_A		(0.357 + 0.063)
#define TYPE_B		(0.085 + 0.015)
#define TYPE_AB		(0.034 + 0.006)

int mismatch[4][4] = {
			{ 0, 0, 0, 0 },		// O ---> O, A, B, AB
			{ 1, 0, 1, 0 }, 	// A ---> A, AB
			{ 1, 1, 0, 0 },		// B ---> B, AB
			{ 1, 1, 1, 0 } 		// AB --> AB
		     };		


int arrayindex( double x )
{
	if ( x < TYPE_O ) return 0;
	else if ( x < TYPE_O + TYPE_A ) return 1;
	else if ( x < TYPE_O + TYPE_A + TYPE_B ) return 2;
	else return 3;
}

char legend[] = "Computer simulation of random blood transfusions in the U.S."; 

int main( int argc, char **argv )
{
	// display the program legend
	printf( "\n\n %8s %s\n\n", " ", legend );

	// seed the random number generator
	srand48( time( NULL ) );

	// perform the demonstration three times
	for (int demo = 0; demo < 3; demo++)
		{
		int	n_samples = 0;
		int	n_rejects = 0;
		for (int rep = 0; rep < MAX_REPS; rep++)
			{
			double	x = drand48();	// choose a donor
			double	y = drand48();	// choose a recipient
			int	i = arrayindex( x );
			int	j = arrayindex( y );
			if ( mismatch[i][j] > 0 ) n_rejects += 1;
			n_samples += 1;
			
			switch ( n_samples )
				{
				case 100:
				case 1000:
				case 10000:
				case 100000:
				case 1000000:
				float	prob = (n_rejects*1.0)/(n_samples*1.0);
				printf( "\t\ttrials=%-8d ", n_samples );
				printf( "rejection-probability " );
				printf( "= %1.3f \n", prob );
				}
			}
		printf( "\n" );
		}
	printf( "\n" );
}