import numpy as np

#to do with real data; I would have to 
#1. import the observation sequence where each obs represents a tract with a vector = price frequency distribution
#2. calculate emission probabilities for each tract (run rand.py for the size of tract dataset)

obs = ('G100200P', 'G250400P', 'G450600P', 'G650800P', 'G9002000P')

states = ('Residential', 'Commercial', 'Manufacturing')
 
initprob = {'Residential': 0.3, 'Commercial': 0.3, 'Manufacturing': 0.4}
 
 #emisprob was defined by Gaussian Probability Density Function (see rand.py)
emisprob = {
   'Residential' : {'G100200P': 0.37, 'G250400P': 0.37, 'G450600P': 0.59, 'G650800P': 0.83, 'G9002000P': 0.05},
   'Commercial' : {'G100200P': 0.58, 'G250400P': 0.31, 'G450600P': 0.83, 'G650800P': 0.43, 'G9002000P': 0.49},
   'Manufacturing': {'G100200P': 0.15, 'G250400P': 0.97, 'G450600P': 0.85, 'G650800P': 0.97, 'G9002000P': 0.31}
   }
 
transprob = {
   'Residential' : {'Residential': 0.5, 'Commercial': 0.4, 'Manufacturing': 0.1},
   'Commercial' : {'Residential': 0.3, 'Commercial': 0.5, 'Manufacturing': 0.2},
   'Manufacturing' : {'Residential': 0.1, 'Commercial': 0.4, 'Manufacturing': 0.5}
   }

def viterbi(obs, states, initprob, emisprob, transprob):
	U = [{}] #initialize an empty set
	path = {} #initialize an empty path

#initialize 
	for z in states:
		U[0][z] = initprob[z] * emisprob[z][obs[0]]
		path[z] = [z]

#When t > 0, perform Viterbi
	#for time through length of observation
	for	t in range(1, len(obs)):
		#U=MIU, this is the saved MAX
		U.append({})
		newpath = {}
	
		for z in states:
			#Recursion, gives max stored in "prob".
			(prob, state) = max((emisprob[z][obs[t]] * transprob[z0][z] * U[t-1][z0], z0) for z0 in states)
			#update U as prob
			U[t][z] = prob
			#update newpath. Keep track of max seq path at each step
			newpath[z] = path[state] + [z] 
		#throw away old path
		path = newpath
	n = 0
	if len(obs) !=1:
		n=t
	print_dptable(U)
	(prob, state) = max((U[n][z], z) for z in states)
	return (prob, path[state])

def print_dptable(U):
    s = "    " + " ".join(("%5d" % i) for i in range(len(U))) + "\n"
    for z in U[0]:
        s += "%.5s: " % z
        s += " ".join("%.5s" % ("%f" % u[z]) for u in U)
        s += "\n"
    print(s)

def run():
    return viterbi(obs, states, initprob, emisprob, transprob)
print(run())