# -*- coding: utf-8 -*-
"""
Created on Sun Nov 06 08:49:00 2016

@author: wderousse
"""

import sys
import xlrd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

np.set_printoptions(threshold=np.nan)
print sys.version
print __name__

filename='linear_6.xlsx'
print
print 'Opening',filename
ss=xlrd.open_workbook(filename,on_demand=True)
sheet_index=0
print ' Worksheet',ss.sheet_names()[sheet_index]
sh=ss.sheet_by_index(sheet_index)

print ' Reading data values'
stoch=[]
for row_index in range(0,sh.nrows):
    stoch.append([sh.cell_value(row_index,0),sh.cell_value(row_index,1),sh.cell_value(row_index,2)])
stoch=np.matrix(stoch)
print ' Stochastic model size:',len(stoch)
print
state0=[]
for row_index in range(0,sh.nrows):
    state0.append([sh.cell_value(row_index,4)])
state0=np.array(state0)
print ' State model size:',len(state0)

print
print 'Stoch:'
print stoch
print 'State x0:'
print state0
print 'State x1:'
state1=stoch*state0
print state1
print 'State x2:'
state2=stoch*state1
print state2

vector_sub1=[]
vector_sub2=[]
vector_sub3=[]
vector_error_sub1=[]
vector_error_sub2=[]
vector_error_sub3=[]
print
vector_sub1.append(state0[0])
vector_sub2.append(state0[1])
vector_sub3.append(state0[2])
vector_error_sub1.append([0])
vector_error_sub2.append([0])
vector_error_sub3.append([0])
max_error=1e-3
ss_state_prev=state0
ss_state_now=stoch*ss_state_prev
i=0
while abs(ss_state_prev[0]-ss_state_now[0]) > max_error:
    vector_sub1.append(ss_state_now[0])
    vector_sub2.append(ss_state_now[1])
    vector_sub3.append(ss_state_now[2])
    vector_error_sub1.append(ss_state_prev[0]-ss_state_now[0])
    vector_error_sub2.append(ss_state_prev[1]-ss_state_now[1])
    vector_error_sub3.append(ss_state_prev[2]-ss_state_now[2])
    i=i+1
    ss_state_prev=ss_state_now
    ss_state_now=stoch*ss_state_prev
print
print 'Max Error: ', max_error
print 'Final Error: ', ss_state_now[0]-ss_state_prev[0]
print 'Number of Iterations: ', i
print
print 'Steady State Vector:'
print ss_state_now

t=range(0,i+1)
plt.plot(t,vector_sub1,t,vector_sub2,t,vector_sub3)
plt.ylabel('# people')
plt.xlabel('time interval')
plt.show()
plt.clf()  
plt.plot(t,vector_error_sub1,t,vector_error_sub2,t,vector_error_sub3)
plt.ylabel('error')
plt.xlabel('time interval')
plt.show()