-
Notifications
You must be signed in to change notification settings - Fork 3
/
hmm_bt.py
635 lines (587 loc) · 22.6 KB
/
hmm_bt.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
#!/usr/bin/python
#
# Utilities for BT-HMM_
#
import numpy as np
from random import *
import matplotlib.pyplot as plt
import editdistance as ed #sudo pip install editdistance
from tqdm import tqdm
import os
import sys
#sudo pip install scikit-learn # dep for hmmlearn
#pip install -U --user hmmlearn
from hmmlearn import hmm
import random as random
import abt_constants as ac
import abtclass as abtc
def outputAmat(A,title,names,of=sys.stdout):
print >> of, title # eg, "Original A matrix:"
for i in range(A.shape[0]):
print >> of, '{0: <7}'.format(names[i]),
for j in range(A.shape[1]):
print >> of, '{:.5f} '.format(A[i,j]),
print >> of, '\n'
def A_row_check(A,of):
# this is like a test but will not stop the program (see A_row_test() below)
print >> of, "A-matrix row check"
eps = 1.0E-6 # accuracy
result_ok = True
for i in range(A.shape[0]):
r = 0
for j in range(A.shape[1]):
if A[i,j] < 0.0:
r += 10000000
result_OK = False
r += A[i,j]
print >> of, i,r
if abs(r-1.0) > eps:
result_ok = False
print >> of, 'Problem: row ',i,' of A-matrix sum is != 1.0 -or- row contains a P<0 sum = ', r
return result_ok
def A_row_test(A,of):
# an actual test using asserts
eps = 1.0E-6 # accuracy
#print 'A-matrix row test'
for i in range(A.shape[0]):
r = 0
for j in range(A.shape[1]):
assert A[i,j] >= 0.0, ' A matrix Prob value < 0.0!'
assert A[i,j] <= 1.0, ' A matrix Prob value > 1!'
r += A[i,j]
#print 'assertion:', i,r
assert abs(r-1.0) < eps, 'Assert Problem: a row sum of A-matrix is != 1.0, sum = '+str(r)
#def A_non_zero(A):
#'''
#make sure trans matrix has no 0.0 elements
#'''
#eps = 1.0E-5 # min value
##print 'A-matrix row test'
#for i in range(A.shape[0]):
#for j in range(A.shape[1]):
#if(A[i,j] < eps):
#A[i,j] = eps
#sum = 0.0
#for j in range(A.shape[1]):
#sum += A[i,j]
#for j in range(A.shape[1]):
#A[i,j] /= sum
def HMM_setup(model, toler=0.01, maxiter=20): # New: setup model.B: discrete emission probs.
#print 'Size: A: ', A.shape
# components = states features = observation symbols
l = model.A.shape[0]
#print 'len(Pi): ', len(Pi), l
#M = hmm.GaussianHMM(n_components=l, covariance_type='diag', n_iter=maxiter, tol=toler, init_params='')
#M.typestring = 'GaussianHMM'
# fit all params: params='ste'
# fit only A matrix: params='t'
M = hmm.MultinomialHMM(n_components=l, n_iter=maxiter, params='t', tol=toler, init_params='',verbose=True)
M.typestring = 'MultinomialHMM'
M.n_features = ac.NSYMBOLS
M.n_components = l
M.startprob_ = model.Pi
M.transmat_ = model.A
############################# Gaussian emissions
# set emissionprob below
# .means and .covars are for GaussianHMM()
#tmpmeans = []
#for i in range(len(names)):
#tmpmeans.append( [ outputs[names[i]] ] )
#M.means_ = np.array(tmpmeans)
#M.means_ = 0.5*np.ones(l).reshape([l,1]) # is this a bug??? what about \delta\mu * i???
#m = np.zeros(model.n)
#for i,n in enumerate(model.names):
#m[i] = model.outputs[n]
#m.shape = [l,1]
#M.means_ = m
#print 'means shape: ', M.means_.shape
#tmpcovars = model.sigma * np.ones((l))
#tmpcovars.shape = [l,1]
#M.covars_ = np.array(tmpcovars)
#########################
#sig = 2.0 # HACK!!!
############################# Multinomial emissions
# setup discrete model.B for MultinomialHMM()
# set up a obs density with mean=model.outputs[n].
for i,n in enumerate(model.names): # names of "components" or "states"
tmp_leaf = abtc.aug_leaf(0.500) # dummy leaf to use SetObsDensity() method
tmp_leaf.set_Obs_Density(model.outputs[n], ac.sig)
#print 'mean: ', model.outputs[n]
for j in range(ac.NSYMBOLS): # "features" or "output symbols"
model.B[i,j] = tmp_leaf.Obs[j] # guarantees same P's as ABT(!)
M.emissionprob_ = np.array(model.B.copy()) # docs unclear on this name!!!!
return M
# Kullback-Leibler divergence
def KL_diverge(M,i,j):
return KL_basic(M.emissionprob_[i,:], M.emissionprob_[j,:])
def KL_basic(A,B): # just two vector probs
N = len(A)
assert N==len(B), 'mismatched prob distribs'
pr = 0.0
for k in range(N):
pr += A[k]*np.log2(B[k]/A[k])
return -1.0*pr
def Shannon_Entropy(A): # A is a discrete prob density
se = 0.0
for i in range(len(A)):
if A[i] > 1.0e-8: # avoid numerical undeflow of log2
se -= A[i]*np.log2(A[i])
return se
## JS divergence between all pairs of states
def JS_ALL(M):
ws = 0.0
B1 = M.emissionprob_[0,:]
N = len(B1) # of symbols
n = M.n_components # of distribs/states
#print 'N DISTRIBS', n
#print 'test:', Shannon_Entropy(M.emissionprob_[2,:])
#muh = 0.0
#for j in range(N):
#muh += float(j)*M.emissionprob_[2,j]
#print 'mu: ', muh
wpi = 1.0/n
wP = np.zeros(N)
for i in range(N): # iterate over the symbols
for j in range(n): # sum over the distribs/states
#if M.emissionprob_[j,i] > 1.0e-6:
#print 'adding: ', i, M.emissionprob_[j,i]
wP[i] += M.emissionprob_[j,i]
wP[i] *= wpi
s = 0.0
for j in range(n):
s += wpi*Shannon_Entropy(M.emissionprob_[j,:])
jsa = Shannon_Entropy(wP) - s
return jsa
## Jensen-Shannon divergence
def JS_diverge(HMM,i,j):
B1 = HMM.emissionprob_[i,:]
B2 = HMM.emissionprob_[j,:]
#print 'shape: ', np.shape(B1)[0]
M = np.zeros(np.shape(B1)[0])
#print 'shape(M): ', np.shape(M)
for k in range(len(B2)):
M[k] = 0.500*(B1[k]+B2[k])
return 0.500*(KL_basic(B1,M)+KL_basic(B2,M))
#def JS_divergeALL(HMMs): # an array of HMMs
def HMM_Project(M,i,j):
'''
Compute the projection of observation densities between states i,j for model M
'''
N = M.emissionprob_.shape[1] # number of output symbols
pr = 0.0
for k in range(N): # emission symbols
pr += M.emissionprob_[i,k]*M.emissionprob_[j,k]
return pr
def HMM_ProjectAll(M):
'''
Compute the projection of observation densities between ALL states for model M
'''
N = M.emissionprob_.shape[1]
l = M.transmat_.shape[0]
pr = 0.0
for i in range(l): # state i
for j in range(l): # state j
for k in range(N): # emission symbols
if i != j: # not the diagonal
pr += M.emissionprob_[i,k]*M.emissionprob_[j,k]
pr *= (2/float(l*(l-1))) # number of off diagonal elements
return pr
def HMM_model_sizes_check(M):
#print 'Model size check:'
#print 'Transmatrix: ',M.transmat_.shape
#print 'Outputs: ',M.means_.shape
#print 'done...'
#quit()
l = M.transmat_.shape[0]
fs = 'Your HMM has inconsistent model sizes and will not run: quitting'
assert M.transmat_.shape == (l,l), fs
print 'Model components/states: ', M.n_components, ' features/symbols: ', M.n_features
assert M.emissionprob_.shape == (l, ac.NSYMBOLS), fs+': emissionprob_'
#assert M.means_.shape == (l,1), fs
# Replace ABT transition probabilities with
# random values (only the non-zero elements tho).
def HMM_ABT_to_random(M):
# A matrix
A = M.transmat_
[r1, c1] = A.shape
r1 -= 2 # don't perturb for output states: Os and Of
for r in range(r1):
flag = -1
for c in range(c1):
# second non-zero element of row
#print 'looking at element: ',r,c
#print 'flag = ', flag
if flag > 0 and A[r][c] > 0:
A[r][c] = 1.0 - flag
#print 'setting second element to', 1.0 - flag
# first non-zero element of row
elif A[r][c] > 0:
if abs(A[r][c] - 1.0) < 0.000001: # don't mess with 1.0 transitions
continue
A[r][c] = random.random() # Uniform(0.0-1.0)
flag = A[r][c] # store value (for use above)
M.transmat_ = A # maybe unnecessary??
#
# initialize A-matrix to all NxN elements random
# (subject to Sum(row) == 1)
#
def HMM_fully_random(model):
A_rand = model.A.copy()
[rn,cn] = A_rand.shape
for r in range(rn):
rsum = 0.0
for c in range(cn):
A_rand[r][c] = random.random()
rsum += A_rand[r][c]
for c in range(cn): # normalize the rows
A_rand[r][c] /= rsum
# randomize means of the output observations
B = np.zeros(model.n)
for i,n in enumerate(model.names):
B[i] = int( 0.5 + ac.FIRSTSYMBOL + (ac.NSYMBOLS-ac.FIRSTSYMBOL)*random.random() )
B.shape = [rn,1] # match req. of hmmlearn
return A_rand, B
# apply a delta (random +-) to the elements of A
# subject to sum of row = 1.]
#
# NEW: if d > 5 initialize A matrix to RANDOM values
#
def HMM_perturb(M, d, model=abtc.model(1)):
''' M = an hmmlearn HMM object
d = perturbation (0 < d < 1.0)
model = model parameters
'''
assert len(model.names) > 1, 'HMM_perturb() [in hmm_bt.py] must be called with a model (2nd argument)'
A = M.transmat_
#np.save("M_trans",M.transmat_)
#np.save("Means",M.means_)
[r1, c1] = A.shape
r1 -= 2 # don't perturb for Os and Of states
for r in range(r1): # go through the rows
flag = -1
rowcnt = 0 # how many non-zero elements in this row?
for c in range(c1):
if A[r][c] > 0.0000001:
rowcnt += 1
##assert rowcnt > 1, 'only 1 non-zero element should not occur ('+str(rowcnt)+') - quitting'
if rowcnt == 2: # ABT type models and SLR models
for c in range(c1):
# second non-zero element of row
#print 'looking at element: ',r,c
#print 'flag = ', flag
if flag > 0 and A[r][c] > 0:
A[r][c] = 1.0 - flag
#print 'setting second element to', 1.0 - flag
# first non-zero element of row
elif A[r][c] > 0.0000001:
if abs(A[r][c] - 1.0) < 0.000001: # don't mess with 1.0 transitions
continue
change = randsign() * d
#print 'Applying change 1.0 + ',change
pbef = A[r][c]
A[r][c] *= (1.0 + change)
paft = A[r][c]
#print 'Actual Change: ', (paft-pbef)/pbef
if A[r][c] > 0.9999:
A[r][c] = 0.9999 # don't allow going to 1.0 or above
if A[r][c] < 0.0000001: # don't allow negative
A[r][c] = 0.0000001
flag = A[r][c] # store value (for use above)
elif rowcnt == 3: #ABT + duration type models
flag = 0
for c in range(c1):
if c > r: # only above diagonal (don't change A[r,r]
if flag > 0 and A[r][c] > 0:
# we've found the second transition to one of two next states
A[r][c] = (1.0-A[r][r]) - flag # keep sum == 1.0
elif A[r][c] > 0.0000001:
if abs(A[r][c] - 1.0) < 0.000001: # don't mess with 1.0 transitions
continue
# we've found the first transition to one of two next states
change = randsign() * d
#print 'Applying change 1.0 + ',change
pbef = A[r][c]
A[r][c] *= (1.0 + change)
paft = A[r][c]
#print 'Actual Change: ', (paft-pbef)/pbef
if A[r][c] > 0.9999:
A[r][c] = 0.9999 # don't allow going to 1.0 or above
if A[r][c] < 0.0000001: # don't allow negative
A[r][c] = 0.0000001
flag = A[r][c] # store value (for use above)
else:
#print 'I dont know how to perturb ', rowcnt, ' non-zero values in a row'
#quit()
pass # if there is a 1.0 transition (100% success or fail) just leave it alone
# Perturb B matrix means. Each mean must be perturbed by same amount, not by a 1+delta as above
# because before, some states had bigger probability errors than others.
sigma = 2.0 # HACK
bdelta = 2 * d * sigma # factor of 2 just feels better(!)
#
# New coding for Multinomial - explicit probs over the symbol integers
# (some kind of "shift"??) or just regenerate.
#
if M.typestring == 'GaussianHMM':
B = M.means_
for i,b in enumerate(B):
#B[i] = int(0.5 + b * (1.0 + randsign() * d))
B[i] += randsign() * bdelta
M.means_ = B
elif M.typestring == 'MultinomialHMM':
############################# Multinomial emissions
# setup discrete model.B for MultinomialHMM()
for i,n in enumerate(model.names):
tmp_leaf = abtc.aug_leaf(0.500) # dummy leaf to use SetObsDensity() method
#perturb the mean by bdelta before generating the emission probs.
newmean = model.outputs[n] + randsign() * bdelta
#print 'Setting mean for state ',i, 'from ' , model.outputs[n], ' to ', newmean
tmp_leaf.set_Obs_Density(newmean, ac.sig)
for j in range(ac.NSYMBOLS):
model.B[i,j] = tmp_leaf.Obs[j] # guarantees same P's as ABT(!)
M.emissionprob_ = np.array(model.B.copy()) # docs unclear on this name!!!!
else:
print 'Unknown model typestring'
quit()
return
def randsign():
a = random.random()
if a > 0.500:
return 1
else:
return -1
# read in observation sequences data file
def read_obs_seqs(logf):
#logf = open(fn,'r')
X = [] # state names
Y = [] # observations
Ls =[] # lengths
seq = [] # current state seq
os = [] # current obs seq
for line in logf:
#print '>>>',line
line = line.strip()
if line == '---':
Ls.append(len(os))
os = []
else:
[state, obs ] = line.split(',')
X.append(state)
Y.append([int(obs)])
os.append([int(obs)])
Y=np.array(Y).reshape(-1,1) # make 2D
Ls = np.array(Ls)
logf.close()
return [X,Y,Ls]
######################################################
#
# Print an A matrix comparison/diff report
#
def Adiff_Report(A1,A2,names,of=sys.stdout):
[e,e2,em,N2,im,jm,anoms,erasures] = Adiff(A1, A2, names)
N = len(names)
print >> of, 'RMS A-matrix error: {:.3f}'.format(e)
print >> of, 'RMS A-matrix error: {:.8f} (only the {:d}/{:d} non zero elements)'.format(e2,N2,N*N)
print >> of, 'Max abs A-matrix error: {:.3f} (at {:d} to {:d})'.format(em,im,jm)
if len(anoms) == 0:
anoms = 'None'
print >> of, 'Anomalies: ', anoms
if len(erasures) == 0:
erasures = 'None'
print >> of, 'Erasures : ', erasures
def Adiff(A1,A2,names): # from 8/28
e = 0
em = -9.9999E100
imax = np.nan
jmax = np.nan
e2 = 0 # avg error of NON ZERO elements
N = A1.shape[0]
assert A1.shape == A2.shape, 'Adiff: A-matrix size mismatch!'
assert not np.isnan(A1).all(), 'A-matrix contains nan'
assert not np.isnan(A2).all(), 'A-matrix contains nan'
#print 'Adiff: A shape: ', A1.shape
#print "A1: "
#print A1
#print "A2: "
#print A2
N2 = 0 # count the non-zero Aij entries
# should be 2(l+2) for ABT-based matrices
anoms = [] #identification
erasures = []
for i in range(N):
for j in range(N):
e1 = (A1[i,j]-A2[i,j])**2
#print 'error: ', e1,i,j
#print 'A1[ij] ',A1[i,j], ' A2[ij] ',A2[i,j], (A1[i,j]-A2[i,j])
if(e1 > em):
em = e1
imax = i
jmax = j
#print "storing emax: ", em, i,j
if(A1[i,j] > 0.000001):
e2 += e1
N2 += 1
e += e1
if(A1[i,j]==0.0 and A2[i,j]>0.0): # for regular but not random perturbations
anoms.append([i,j])
if(A1[i,j]>0.0 and A2[i,j] < 0.0000001):
erasures.append([names[i],names[j]])
e = np.sqrt(e/(N*N)) # div total number of Aij elements
e2 = np.sqrt(e2/N2) # RMS error of NON zero A[1]ij
em = np.sqrt(em) # Max error
#print 'imax, jmax; ', imax, jmax
return [e,e2,em,N2,imax,jmax,anoms,erasures]
###############################################
# compare two A-matrices (NEW)
#
#def Adiff(A1,A2,names):
#e = 0
#e_abs_total = 0.0
#em = -99999.9
#e2 = 0 # avg error of NON ZERO elementstoler=0.01, maxiter=20
#N = A1.shape[0]
##print 'Adiff: A shape: ', A1.shape
#N2 = 0 # count the non-zero Aij entries
## should be 2(l+2) of course
#anoms = [] #identification
#erasures = []
#for i in range(N-2): # skip last two rows which are 1.000
#for j in range(N):
#e1 = (A1[i,j]-A2[i,j])**2
#ea = abs(A1[i,j]-A2[i,j])
##print 'error: ', e1,i,j
##print 'A1[ij] ',A1[i,j], ' A2[ij] ',A2[i,j], (A1[i,j]-A2[i,j])
#if(ea > em): # should be absolute error not e^2
#em = ea
#imax = i+1 # change from array index to state numbers
#jmax = j+1
##print "storing emax: ", em, i,j
#if(A1[i,j] > 0.000001):
#e2 += ea # accumulate error for non-zero Aij
#N2 += 1
#e_abs_total += ea
#if(A1[i,j]==0.0 and A2[i,j]>0.0): # NOTE: implies direction btwn A1 and A2
#anoms.append([i,j])
#if(A1[i,j]>0.0 and A2[i,j] < 0.0000001):
#erasures.append([names[i],names[j]])
#e = (e_abs_total/(N*N)) # div total number of Aij elements
#e2 = (e2/N2) # RMS error of NON zero Aij
##print 'imax, jmax; ', imax, jmax
#return [e,e2,em,N2,imax,jmax,anoms,erasures]
###############################################################
# Evaluation of Viterbi
#
# Blake's new start
#
def Veterbi_Eval(p,x,names,l,statenos):
'''
p = state sequence estimates (concatenated state seqs, np.array() of numbers)
x = true state sequences (list of numbers)
names = list of state names (Nx1)
l = lengths of each state sequence in 'x'
statenos =
'''
nseqs = len(l)
maxd = -999999
count = 0
totald = 0
states_visited = set(x)
assert len(names) == len(states_visited), 'Viterbi Evaluation: wrong states/state-names'
p1 = list(p)
offset = 0
for i in range(len(l)): # iterate over sequences
st1 = ''
st2 = ''
for j in range(l[i]): # iterate over symbols in seq i
st1 += str(p[offset+j])
st2 += str(x[offset+j])
d1 = ed.eval(st1,st2)
d = d1 / float(len(st2)) # Str edit distance / length of true seq.
if d > maxd:
maxd = d
totald += d
count += d1
offset+=j+1
avgd = totald / float(nseqs)
return [avgd, maxd, count]
# avgd: average sed per symbol
# maxd: max sed for any seq
# count: total sed for all seqs
##############################################
#Forward Pass
##############################################
def Foward_eval(obs,l,M):
# counter = 0
# os.system('clear')
# print "######################################",obs.shape
# obs_sequence = 0
# while counter < len(l):
# obs_slice = obs[counter:(l[obs_sequence]-1)]
# foward = np.zeros((M.transmat_.shape[0],l[obs_sequence]))
# for s in range(M.transmat_.shape[0]):
# forward[s][0] = M.startprob_[s] * obs_slice[0,s]
# for t in range (1,len(obs.slice)):
# for s in range (M.transmat_.shape[0]):
# for last_step in range(M.transmat_.shape[0]):
# foward[s,t] += foward[t-1,last_step]*M.transmat_[last_step,s]*obs_slice[t,s]
counter = 0
logprob = 0
log_avg = 0
for i in range(len(l)):
sample = obs[counter:counter+l[i]]
logprob += M.score(sample,[l[i]])
counter += l[i]
log_avg = logprob/len(l)
return log_avg
######################################################
#Plotter
######################################################
def Plotter(master,y):
#Foward_eval
ffig, fax = plt.subplots(1,1)
for run in range(master.shape[1]):
labeler = "Run Number "+ str(run)
fax.plot(y,master[0,run,:,0],label = labeler)
fax.set(ylabel='Average Log Probability',xlabel='HMM Delta',title = 'Forward')
fax.grid()
vfig, vax = plt.subplots(1,1)
vfig2, vax2 = plt.subplots(1,1)
for run in range(master.shape[1]):
labeler = "Run Number "+ str(run)
vax.plot(y,master[1,run,:,0], label = labeler)
vax.set(ylabel='Total Edit Distance',xlabel='HMM Delta',title = 'Viterbi')
vax.grid()
for run in range(master.shape[1]):
labeler = "Run Number "+ str(run)
vax2.plot(y,master[1,run,:,2], label = labeler)
vax2.set(ylabel='Number of Exact Matches',xlabel='HMM Delta',title = 'Viterbi')
vax2.grid()
bfig, bax = plt.subplots(1,1)
bfig2,bax2 = plt.subplots(1,1)
bfig3,bax3 = plt.subplots(1,1)
for run in range(master.shape[1]):
labeler = "Run Number "+ str(run)
bax.plot(y,master[2,run,:,0], label = labeler)
bax.set(ylabel='Eaverage Distance',xlabel='HMM Delta',title = 'BaumWelch')
bax.grid()
for run in range(master.shape[1]):
labeler = "Run Number "+ str(run)
bax2.plot(y,master[2,run,:,1], label = labeler)
bax2.set(ylabel='EAinfty',xlabel='HMM Delta',title = 'BaumWelch')
bax2.grid()
for run in range(master.shape[1]):
labeler = "Run Number "+ str(run)
bax3.plot(y,master[2,run,:,2], label = labeler)
bax3.set(ylabel='EMax',xlabel='HMM Delta',title = 'BaumWelch')
bax3.grid()
plt.legend(handles=[fax,vax,bax ])
ffig.savefig("forward.png")
vfig.savefig("Veterbi.png")
vfig2.savefig("Veterbi2.png")
bfig.savefig("BaumWelch.png")
bfig2.savefig("BaumWelch2.png")
bfig3.savefig("BaumWelch3.png")
#print "shapes:"
#print "outputs", len(outputs)
#print "means_", (M.means_.shape)
#print "covars", (tmpcovars.shape)
#print "trans", (M.transmat_.shape)