/* Discrete-event simulation of Markov model teenager    */
/*                                                       */
/* Usage: cc -o teensim teensim.c -lm                    */
/*        ./teensim                                      */
/*                                                       */
/* 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 ASLEEP 0
#define GRUMPY 1
#define HUNGRY 2
#define JOY 3
#define TIRED 4

/* Transition rates for Markov chain model */
#define LAMBDA01 3.0   /* transition rate from Asleep to Grumpy */
#define LAMBDA02 5.0   /* transition rate from Asleep to Hungry */
#define LAMBDA04 2.0   /* transition rate from Asleep to Tired */
#define LAMBDA14 4.0   /* transition rate from Grumpy to Tired */
#define LAMBDA21 4.0   /* transition rate from Hungry to Grumpy */
#define LAMBDA23 1.0   /* transition rate from Hungry to Joy */
#define LAMBDA31 10.0  /* transition rate from Joy to Grumpy */
#define LAMBDA34 10.0  /* transition rate from Joy to Tired */
#define LAMBDA40 3.0   /* transition rate from Tired to Asleep */

#define END_OF_TIME 10000.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;
    float p0, p1, p2, p3, p4; 
    float option1, option2, option3;
    int current_state, next_state;
    int i, j, counts[5][5];

    /* 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;
    p0 = 0.0;
    p1 = 0.0;
    p2 = 0.0;
    p3 = 0.0;
    p4 = 0.0;
    current_state = ASLEEP;
    next_state = ASLEEP;
    nextevent = now;
    for( i = 0; i < 5; i++ )
      for( j = 0; j < 5; j++ )
	counts[i][j] = 0;
    
    while( now < END_OF_TIME )
      {
	now = nextevent;
#ifdef DEBUG
	printf("Top of loop: crediting state %d with %f from %f to %f\n",
	       current_state, now - last, last, now);
#endif	
	/* update statistics */
	if( current_state == ASLEEP )
	  p0 += now - last;
	else if( current_state == GRUMPY )
	  p1 += now - last;
	else if( current_state == HUNGRY )
	  p2 += now - last;
	else if( current_state == JOY )
	  p3 += now - last;
	else if( current_state == TIRED )
	  p4 += now - last;
	else printf("Unknown teenager state %d at end!!!\n", current_state);
	
	/* make the transition */
#ifdef DEBUG
	printf("Middle of loop: switching from state %d to state %d\n",
	       current_state, next_state);
#endif
	counts[current_state][next_state]++;
	current_state = next_state;
	last = now;

	/* determine next transition */
	if( current_state == ASLEEP )
	  {
	    /* determine what happens sooner */
	    option1 = Exponential(1.0/LAMBDA01);
	    option2 = Exponential(1.0/LAMBDA02);
	    option3 = Exponential(1.0/LAMBDA04);
#ifdef DEBUG
	    printf("ASLEEP options are %f to G or %f to H or %f to T\n",
		   option1, option2, option3);
#endif
	    if( (option1 < option2) && (option1 < option3) )
	      {
		nextevent = now + option1;
		next_state = GRUMPY;
	      }
	    else if( (option2 < option1) && (option2 < option3) )
	      {
		nextevent = now + option2;
		next_state = HUNGRY;
	      }
	    else
	      {
		nextevent = now + option3;
		next_state = TIRED;
	      }
	  }
	else if( current_state == GRUMPY )
	  {
	    nextevent = now + Exponential(1.0/LAMBDA14);
	    next_state = TIRED;
	  }
	else if( current_state == HUNGRY )
	  {
	    /* determine what happens sooner */
	    option1 = Exponential(1.0/LAMBDA21);
	    option2 = Exponential(1.0/LAMBDA23);
#ifdef DEBUG
	    printf("HUNGRY options are %f to G or %f to J\n", option1, option2);
#endif
	    if( option1 < option2 )
	      {
		nextevent = now + option1;
		next_state = GRUMPY;
	      }
	    else
	      {
		nextevent = now + option2;
		next_state = JOY;
	      }
	  }
	else if( current_state == JOY )
	  {
	    /* determine what happens sooner */
	    option1 = Exponential(1.0/LAMBDA31);
	    option2 = Exponential(1.0/LAMBDA34);
#ifdef DEBUG
	    printf("JOY options are %f to G or %f to T\n", option1, option2);
#endif
	    if( option1 < option2 )
	      {
		nextevent = now + option1;
		next_state = GRUMPY;
	      }
	    else
	      {
		nextevent = now + option2;
		next_state = TIRED;
	      }
	  }
	else if( current_state == TIRED )
	  {
	    nextevent = now + Exponential(1.0/LAMBDA40);
	    next_state = ASLEEP;
	  }
	else printf("Unknown teenager 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("Done big loop: crediting state %d with %f from %f to %f\n",
	       current_state, now - last, last, now);
#endif	
    /* calculate and report statistics */
    if( current_state == ASLEEP )
      p0 += now - last;
    else if( current_state == GRUMPY )
      p1 += now - last;
    else if( current_state == HUNGRY )
      p2 += now - last;
    else if( current_state == JOY )
      p3 += now - last;
    else if( current_state == TIRED )
      p4 += now - last;
    else printf("Unknown teenager state %d at end!!!\n", current_state);

    /* convert state time into a normalized probability */
    p0 /= now;
    p1 /= now;
    p2 /= now;
    p3 /= now;
    p4 /= now;

    printf("\n");
    printf("Simulated Markov Chain Model of Teenager\n");
    printf("\n");
    printf("  Prob(Asleep): p0 = %f\n", p0);
    printf("  Prob(Grumpy): p1 = %f\n", p1);
    printf("  Prob(Hungry): p2 = %f\n", p2);
    printf("  Prob(Joy):    p3 = %f\n", p3);
    printf("  Prob(Tired):  p4 = %f\n", p4);
    printf("\n");
    printf("Sanity check: p0 + p1 + p2 + p3 + p4 = %f\n",
	   p0 + p1 + p2 + p3 + p4);
    printf("\n");

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