/* Discrete-event simulation of Baseball with DoublePlays   */
/*                                                       */
/* Usage: cc -o bbsim2 bbsim2.c -lm                      */
/*        ./bbsim2                                       */
/*                                                       */
/* Written by Carey Williamson, University of Calgary    */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h> /* use man 3 log for the math functions */

/* States in Markov chain model */
#define NOBODY_OUT 0
#define ONE_OUT 1
#define TWO_OUT 2
#define THREE_OUT 3

/* what types of double plays are allowed */
#define FIRST 1
#define SECOND 1

/* Transition rates for Markov chain model */
#define LAMBDA1 0.20   /* transition rate to get first out */
#define LAMBDA2 0.25   /* transition rate to get second out */
#define LAMBDA3 0.25   /* transition rate to get third out */
#define LAMBDA4 0.50   /* transition rate for TV commercial */

#define END_OF_TIME 100000.0 /* Time far in the future */

/* #define DEBUG 1   /* verbose debugging */

/***********************************************************************/
/*                 RANDOM NUMBER GENERATION STUFF                      */
/***********************************************************************/

/* Parameters for random number generation. */
#define MAX_INT 2147483647       /* Maximum positive integer 2^31 - 1 */

/* Generate a random floating point value uniformly distributed in [0,1] */
float Uniform01()
  {
    float randnum;

    /* get a random positive integer from random() */
    randnum = (float) 1.0 * random();

    /* divide by max int to get something in 0..1  */
    randnum = (float) randnum / (1.0 * MAX_INT); 

    return( randnum );
  }

/* Generate a random floating point number from an exponential    */
/* distribution with mean mu.                                     */

float Exponential(mu)
    float mu;
  {
    float randnum, ans;

    randnum = Uniform01();
    ans = -(mu) * log(randnum);

    return( ans );
  }

/***********************************************************************/
/*                 MAIN PROGRAM                                        */
/***********************************************************************/

int main()
  {
    float now, last, nextevent, option1, option2;
    float p0, p1, p2, p3; 
    int current_state, next_state;
    int i, j, counts[4][4];
    float lambdaD;
    int doubleplays;

    printf("Enter double-play rate (lambdaD): ");
    scanf("%f", &lambdaD);

    /* Seed the PRNG */
    srandom(1234567);       /* for a fixed seed every time you run it */
    srandom(time(NULL));    /* for different seed every time you run it */

    /* Initialization */
    now = 0.0;
    last = 0.0;
    nextevent = 0.0;
    p0 = 0.0;
    p1 = 0.0;
    p2 = 0.0;
    p3 = 0.0;
    doubleplays = 0;
    current_state = NOBODY_OUT;
    for( i = 0; i < 4; i++ )
      for( j = 0; j < 4; j++ )
	counts[i][j] = 0;
    
    while( now < END_OF_TIME )
      {
	now = nextevent;
#ifdef DEBUG
	printf("Top of loop: game was in state %d for time %f from %f to %f\n",
	       current_state, now - last, last, now);
#endif	
	/* update statistics */
	if( current_state == NOBODY_OUT )
	  p0 += now - last;
	else if( current_state == ONE_OUT )
	  p1 += now - last;
	else if( current_state == TWO_OUT )
	  p2 += now - last;
	else if( current_state == THREE_OUT )
	  p3 += now - last;
	else printf("Unknown baseball state %d here!!!\n", current_state);
	
	/* make the transition */
#ifdef DEBUG
	printf("Middle of loop: changing from state %d to state %d at time %f\n",
	       current_state, next_state, now);
#endif	
	counts[current_state][next_state]++;
	current_state = next_state;
	last = now;

	/* determine next transition */
	if( current_state == NOBODY_OUT )
	  {
#ifdef FIRST
	    /* determine what happens sooner */
	    option1 = Exponential(1.0/lambdaD);
	    option2 = Exponential(1.0/LAMBDA1);
#ifdef DEBUG
	    printf("Options are %f or %f\n", option1, option2);
#endif
	    if( option1 < option2 )
	      {
#ifdef DEBUG
		printf("DP1!!!\n");
#endif
		doubleplays++;
		nextevent = now + option1;
		next_state = TWO_OUT;
	      }
	    else
	      {
		nextevent = now + option2;
		next_state = ONE_OUT;
	      }
#else
	    nextevent = now + Exponential(1.0/LAMBDA1);
	    next_state = ONE_OUT;
#endif
	  }
	else if( current_state == ONE_OUT )
	  {
#ifdef SECOND
	    /* determine what happens sooner */
	    option1 = Exponential(1.0/lambdaD);
	    option2 = Exponential(1.0/LAMBDA2);
#ifdef DEBUG
	    printf("Options are %f or %f\n", option1, option2);
#endif
	    if( option1 < option2 )
	      {
#ifdef DEBUG
		printf("DP2!!!\n");
#endif
		doubleplays++;
		nextevent = now + option1;
		next_state = THREE_OUT;
	      }
	    else
	      {
		nextevent = now + option2;
		next_state = TWO_OUT;
	      }
#else
	    nextevent = now + Exponential(1.0/LAMBDA2);
	    next_state = TWO_OUT;
#endif
	  }
	else if( current_state == TWO_OUT )
	  {
	    nextevent = now + Exponential(1.0/LAMBDA3);
	    next_state = THREE_OUT;
	  }
	else if( current_state == THREE_OUT )
	  {
	    nextevent = now + Exponential(1.0/LAMBDA4);
	    next_state = NOBODY_OUT;
	  }
	else printf("Unknown baseball state %d at time %f!!!\n",
		    current_state, now);
	
#ifdef DEBUG
	printf("Bottom of loop: next state will be %d at time %f\n", next_state, nextevent);
#endif
      }

#ifdef DEBUG
    printf("\n");
    printf("End time: %f   p0: %f   p1: %f   p2: %f   p3: %f\n",
	   now, p0, p1, p2, p3);
#endif
    /* convert state time into a normalized probability */
    p0 /= now;
    p1 /= now;
    p2 /= now;
    p3 /= now;

    printf("\n");
    printf("Simulated Markov Chain Model of Baseball with Double Plays\n");
    printf("\n");
    printf("  Prob(nobody out): p0 = %f\n", p0);
    printf("  Prob(one out):    p1 = %f\n", p1);
    printf("  Prob(two out):    p2 = %f\n", p2);
    printf("  Prob(three out):  p3 = %f\n", p3);
    printf("\n");
    printf("Sanity check: p0 + p1 + p2 + p3 = %f\n",
	   p0 + p1 + p2 + p3);
    printf("\n");

    printf("Double play rate: %f\n", lambdaD);
    printf("Double plays observed: %d\n", doubleplays);

    /* print out matric with counts of state transitions */
    printf("\n");
    printf("Matrix with State Transition Counts\n");
    printf("Outs:   0     1     2     3 \n");
    printf("----------------------------\n");
    for( i = 0; i < 4; i++ )
      {
	printf(" %d: ", i);
	for( j = 0; j < 4; j++ )
	  {
	    printf(" %4d ", counts[i][j]);
	  }
	printf("\n");
      }
  }