/* Numerical solution to Markov chain model of a teenager */

/* Usage: cc -o teen teen.c -lm                           */
/*        ./teen                                          */

/* Written by Carey Williamson, January 2019 */

#include <stdio.h>
#include <math.h>

/* Choose which of the following parameter settings to use */
#define REAL 1

/* Transition rates for Markov chain model */
#ifdef REAL
#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 */
#endif

#ifdef TEST0
#define LAMBDA01 1.0   /* transition rate from Asleep to Grumpy */
#define LAMBDA02 1.0   /* transition rate from Asleep to Hungry */
#define LAMBDA04 1.0   /* transition rate from Asleep to Tired */
#define LAMBDA14 1.0   /* transition rate from Grumpy to Tired */
#define LAMBDA21 1.0   /* transition rate from Hungry to Grumpy */
#define LAMBDA23 1.0   /* transition rate from Hungry to Joy */
#define LAMBDA31 1.0   /* transition rate from Joy to Grumpy */
#define LAMBDA34 1.0   /* transition rate from Joy to Tired */
#define LAMBDA40 1.0   /* transition rate from Tired to Asleep */
#endif

#ifdef TEST1
#define LAMBDA01 0.333 /* transition rate from Asleep to Grumpy */
#define LAMBDA02 0.333 /* transition rate from Asleep to Hungry */
#define LAMBDA04 0.333 /* transition rate from Asleep to Tired */
#define LAMBDA14 1.0   /* transition rate from Grumpy to Tired */
#define LAMBDA21 0.5   /* transition rate from Hungry to Grumpy */
#define LAMBDA23 0.5   /* transition rate from Hungry to Joy */
#define LAMBDA31 0.5   /* transition rate from Joy to Grumpy */
#define LAMBDA34 0.5   /* transition rate from Joy to Tired */
#define LAMBDA40 1.0   /* transition rate from Tired to Asleep */
#endif

#ifdef TEST2
#define LAMBDA01 1.0   /* transition rate from Asleep to Grumpy */
#define LAMBDA02 1.667 /* transition rate from Asleep to Hungry */
#define LAMBDA04 0.667 /* transition rate from Asleep to Tired */
#define LAMBDA14 4.0   /* transition rate from Grumpy to Tired */
#define LAMBDA21 2.0   /* transition rate from Hungry to Grumpy */
#define LAMBDA23 0.5   /* transition rate from Hungry to Joy */
#define LAMBDA31 5.0   /* transition rate from Joy to Grumpy */
#define LAMBDA34 5.0   /* transition rate from Joy to Tired */
#define LAMBDA40 3.0   /* transition rate from Tired to Asleep */
#endif

#ifdef TEST3
#define LAMBDA01 2.0   /* transition rate from Asleep to Grumpy */
#define LAMBDA02 3.333 /* transition rate from Asleep to Hungry */
#define LAMBDA04 1.333 /* transition rate from Asleep to Tired */
#define LAMBDA14 8.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 6.0   /* transition rate from Tired to Asleep */
#endif

#define DEBUG 1   /* verbose debugging */

int main ()
  {
    float p0, p1, p2, p3, p4; 
    float ratio, sum;

    /* Assume p0 = 1 initially, and do relative probs. Renormalize later. */

    p0 = 1.0;
#ifdef DEBUG
    printf("p0 is %f\n", p0);
#endif

    ratio = LAMBDA02;
    ratio /= (LAMBDA21 + LAMBDA23);
    p2 = p0 * ratio;
#ifdef DEBUG
    printf("p2 is %f times p0\n", ratio);
#endif

    ratio = LAMBDA23;
    ratio /= (LAMBDA31 + LAMBDA34);
    p3 = p2 * ratio;
#ifdef DEBUG
    printf("p3 is %f times p2\n", ratio);
#endif

    p1 = LAMBDA01 * p0 + LAMBDA21 * p2 + LAMBDA31 * p3;
    p1 /= LAMBDA14;
#ifdef DEBUG
    printf("p1 is %f\n", p1);
#endif

    p4 = LAMBDA04 * p0 + LAMBDA14 * p1 + LAMBDA34 * p3;
    p4 /= LAMBDA40;
#ifdef DEBUG
    printf("p4 is %f\n", p4);
#endif

    sum = p0 + p1 + p2 + p3 + p4;
#ifdef DEBUG
    printf("sum is %f\n", sum);
#endif

    /* Renormalize the probabilities now */
    p0 /= sum;
    p1 /= sum;
    p2 /= sum;
    p3 /= sum;
    p4 /= sum;
    
    printf("\n");
    printf("Analytical 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");
  }