Programming Praxis 10 May 2013

/*
Programming Praxis 2013-05-10
see
http://programmingpraxis.com/2013/05/10/mindcipher/#comments
*/

#include <iostream>
#include <vector>
#include <stdexcept>

//some doodads I'll need.
#include <boost/numeric/ublas/matrix.hpp>


/*
    I'll start with #1.

    1: Coin Flip: On Monday, you flip a coin all day. You start flipping it until you see the pattern Head, Tail, Head. You record the number of flips required to reach this pattern, and start flipping again (and counting up from 1 again) until you see that pattern again, you record the second number, and start again. At the end of the day you average all of the numbers you’ve recorded. On Tuesday you do the EXACT same thing except you flip until you see the pattern Head, Tail, Tail.

    Will Monday’s number be higher than Tuesday’s, equal to Tuesday’s, or lower than Tuesday’s?

    Obviously the answer is 'I can't possibly know in advance'. But this is probably supposed to be about long-term behaviour, so
    I'll start by getting the means I want.

    The Markov state "the last three flips" can be represented by a 3-bit binary number, which is how I'll proceed.
*/

typedef boost::numeric::ublas::matrix<double,boost::numeric::ublas::row_major,std::vector<double>> matrix_t;
typedef boost::numeric::ublas::vector<double,std::vector<double>> vector_t;

//define some row ops
matrix_t& swap_rows(matrix_t::size_type i,matrix_t::size_type j, matrix_t& m){
    if(i == j) return m;
    for(matrix_t::size_type k = 0; k < m.size2() ; ++k){
        double temp = m(i,k);
        m(i,k) = m(j,k);
        m(j,k) = temp;
    }
    return m;
}

matrix_t& scale_row(matrix_t::size_type i,double scale,matrix_t& m){
    for(matrix_t::size_type k = 0; k < m.size2() ; ++k){
        m(i,k) *= scale;
    }
    return m;
}

matrix_t& addmultiple_rows(matrix_t::size_type i,matrix_t::size_type j,double scale,matrix_t& m){
    for(matrix_t::size_type k = 0 ; k < m.size2() ; ++k){
        m(i,k) += scale*m(j,k);
    }
    return m;
}

//Perform row reduction.
void row_reduce(double e,matrix_t& m,vector_t& v){
    typedef matrix_t::size_type st;

    st row = 0;
    st pivot = 0;
    double temp;
    while(pivot < m.size2()){
        //find the biggest element in the current column.
        st max = row;
        for(st i = row+1 ; i < m.size1() ; ++i)
            if(abs(m(i,pivot)) > abs(m(max,pivot))) max = i;
        swap_rows(row,max,m);
        //also swap in the vector
        temp = v(row);
        v(row) = v(max);
        v(max) = temp;
        //check it's worth operating on.
        if(abs(m(row,pivot)) > e){
            //scale things.
            v(row) *= 1.0/m(row,pivot);
            scale_row(row,1.0/m(row,pivot),m);
            //clean out
            for(st i = 0 ; i < m.size1() ; ++i){
                if(i != row){
                    double scale = -m(i,pivot);
                    v(i) += scale * v(row);
                    addmultiple_rows(i,row,scale,m);
                }
            }
            ++row;
        }
        ++pivot;
    }
}

double mean_time(unsigned char target_state){
    if(target_state>8) throw std::invalid_argument("mean_time parameter not in range.");

    //Build the bits I need for production of mean hitting times.
    //Hitting {target_state}.
    matrix_t mean_matrix(9,9,0.0);
    for(int i = 0 ; i < 8 ; ++i){
        //set up the ith row...
        if(i == target_state) mean_matrix(i,i) = 1.0;
        else{
            int u = (i << 1) & 7;
            int v = u | 1;
            mean_matrix(i,u) += 0.5;
            mean_matrix(i,v) += 0.5;
            mean_matrix(i,i) += -1.0;
        }
        //set up the start row while I'm here.
        mean_matrix(8,i) = 0.125;
    }
    //finish the start row.
    mean_matrix(8,8) = -1.0;

    vector_t rhs(9,-1.0);
    rhs(target_state) = 0.0;
    rhs(8) = -3.0;  //Essentially the first step is three.

    //triangularise mean_matrix.
    row_reduce(1e-10,mean_matrix,rhs);

    //The matrices are all invertible here, so there's no need to perform the next step.
    return rhs(8);
}

int main(){
    for(int i = 0; i < 8 ; ++i)
        std::cout << "Mean time to hit " << ((i&4)?"H":"T") << ((i&2)?"H":"T") << ((i&1)?"H":"T") << " = "<< mean_time(i) << std::endl;
    std::cin.get();
    return 0;
}