// Arup Guha
// 9/3/09
// Example of simulating a Markov process
// This program simulates the Markov chain shown in class on 9/1/09.

#include <stdio.h>
#include <time.h>

#define NUM_STEPS 1000000

int main() {

    srand(time(0));
        
    // I am hard-coding this example. It would be better style to read this in
    // from a file and have it be more general.
    double transMatrix[3][3] = {{0, 1.0/3, 2.0/3},{0,.5, .5},{1.0/3,1.0/3,1.0/3}};
    
    // Will keep track of how many steps we are in each state.
    int counter[3];    
    int i;
    for (i=0; i<3; i++)
        counter[i] = 0;
    
    // Making state 0 our start state.    
    int curState = 0;
    
    for (i=0; i<NUM_STEPS; i++) {
    
        // Add 1 here stating that at this time step we were at curState.
        counter[curState]++;
        
        // Gives me a random number in [0, 1).
        double prob = rand()/32768.0;
        
        int newState = 0;
        double totalProb = 0;
        
        // I am adding up the probabilities, starting from the left, on
        // a single row of the transition matrix. Essentially, my goal is 
        // this: to find the state such that the probability of being in that
        // state or less exceeds prob.
        while (totalProb < prob) {
              
            // Add in this next probability.  
            totalProb += transMatrix[curState][newState];
            
            // Go to the next state.
            newState++;
        }
        
        // We overshot our state by 1 since the loop only exists after going
        // past prob, so we go back to the previous state.
        newState--;
        
        // Set up our new state as the current one for the next iteration.
        curState = newState;
    }
    
    // Print out the probability distribution of ending up in each state.
    for (i=0; i<3; i++) 
        printf("Probability of being in state %d = %.3lf.\n", 
                                i, (double)counter[i]/NUM_STEPS);
                                
    system("PAUSE");
    return 0;
       
}