|
e@0
|
1 #!/usr/bin/python2
|
|
e@0
|
2 # -*- coding: utf-8 -*-
|
|
e@0
|
3 """
|
|
e@0
|
4 Created on Thu Apr 23 11:53:17 2015
|
|
e@0
|
5
|
|
e@0
|
6 @author: mmxgn
|
|
e@0
|
7 """
|
|
e@0
|
8
|
|
e@0
|
9 # This file does the cluster estimation and the removal of outliers
|
|
e@0
|
10
|
|
e@0
|
11 from sys import argv, exit
|
|
e@0
|
12 from essentia.standard import YamlInput, YamlOutput
|
|
e@0
|
13 from essentia import Pool
|
|
e@0
|
14 from pca import *
|
|
e@0
|
15
|
|
e@0
|
16 from numpy import *
|
|
e@0
|
17 from sklearn import cluster
|
|
e@0
|
18 from sklearn.metrics import pairwise_distances
|
|
e@0
|
19 from sklearn.cluster import KMeans, MiniBatchKMeans
|
|
e@0
|
20 from matplotlib.pyplot import *
|
|
e@0
|
21 #from sklearn.mixture import GMM
|
|
e@0
|
22 from sklearn.naive_bayes import GaussianNB, MultinomialNB
|
|
e@0
|
23 from scipy.signal import decimate
|
|
e@0
|
24 from sklearn import cross_validation
|
|
e@0
|
25
|
|
e@0
|
26 #from hmmlearn import hmm
|
|
e@0
|
27 from hmmlearn.hmm import GMM
|
|
e@0
|
28 from hmmlearn import hmm
|
|
e@0
|
29 #from adpcm import adm, adm_reconstruct
|
|
e@0
|
30
|
|
e@0
|
31
|
|
e@0
|
32 mse = lambda A,B: ((array(A)-array(B)) ** 2).mean()
|
|
e@0
|
33
|
|
e@0
|
34
|
|
e@0
|
35 def smooth_matrix_1D(X):
|
|
e@0
|
36 window = scipy.signal.gaussian(51,4)
|
|
e@0
|
37 window = window/sum(window)
|
|
e@0
|
38 intermx = zeros((X.shape[0],X.shape[1]+100))
|
|
e@0
|
39 intermx[:, 50:-50] = X
|
|
e@0
|
40
|
|
e@0
|
41 for m in range(0, X.shape[0]):
|
|
e@0
|
42 # print intermx.shape
|
|
e@0
|
43 intermx[m,:] = convolve(ravel(intermx[m,:]), window,'same')
|
|
e@0
|
44
|
|
e@0
|
45 return intermx[:,50:-50]
|
|
e@0
|
46
|
|
e@0
|
47 def adm_reconstruct(codeword, h, dmin=.01, dmax=.28):
|
|
e@0
|
48 x = zeros((1, codeword.shape[1]))
|
|
e@0
|
49
|
|
e@0
|
50 delta1 = dmin
|
|
e@0
|
51 delta2 = dmin
|
|
e@0
|
52 Sum = h
|
|
e@0
|
53
|
|
e@0
|
54 x[0] = h
|
|
e@0
|
55 for i in range(0, codeword.shape[1]):
|
|
e@0
|
56 if codeword[0,i] == 0:
|
|
e@0
|
57 delta1 = dmin
|
|
e@0
|
58 delta2 = dmin
|
|
e@0
|
59
|
|
e@0
|
60 elif codeword[0,i] == 1:
|
|
e@0
|
61 delta2 = dmin
|
|
e@0
|
62 Sum += delta1
|
|
e@0
|
63 delta1 *= 2
|
|
e@0
|
64 if delta1 > dmax:
|
|
e@0
|
65 delta1 = dmax
|
|
e@0
|
66
|
|
e@0
|
67 elif codeword[0,i] == -1:
|
|
e@0
|
68 delta1 = dmin
|
|
e@0
|
69 Sum -= delta2
|
|
e@0
|
70 delta2 *= 2
|
|
e@0
|
71 if delta2 > dmax:
|
|
e@0
|
72 delta2 = dmax
|
|
e@0
|
73 x[0,i] = Sum
|
|
e@0
|
74 return x
|
|
e@0
|
75
|
|
e@0
|
76 def adm(x, dmin=.01, dmax=.28, tol=0.0001):
|
|
e@0
|
77
|
|
e@0
|
78 # Adaptive delta modulation adapted by code:
|
|
e@0
|
79 # (adeltamod.m)
|
|
e@0
|
80 #
|
|
e@0
|
81 # Adaptive Delta Modulator
|
|
e@0
|
82 # by Gandhar Desai (gdesai)
|
|
e@0
|
83 # BITS Pilani Goa Campus
|
|
e@0
|
84 # Date: 28 Sept, 2013
|
|
e@0
|
85
|
|
e@0
|
86 xsig = x
|
|
e@0
|
87
|
|
e@0
|
88 Lx = len(x)
|
|
e@0
|
89
|
|
e@0
|
90 ADMout = zeros((1, Lx))
|
|
e@0
|
91 codevec = zeros((1, Lx))
|
|
e@0
|
92
|
|
e@0
|
93
|
|
e@0
|
94 Sum = x[0]
|
|
e@0
|
95 delta1 = dmin
|
|
e@0
|
96 delta2 = dmin
|
|
e@0
|
97 mult1 = 2
|
|
e@0
|
98 mult2 = 2
|
|
e@0
|
99 for i in range(0, Lx):
|
|
e@0
|
100 #print abs(xsig[i] - Sum)
|
|
e@0
|
101 if (abs(xsig[i] - Sum) < tol):
|
|
e@0
|
102 bit = 0
|
|
e@0
|
103 delta2 = dmin
|
|
e@0
|
104 delta1 = dmin
|
|
e@0
|
105
|
|
e@0
|
106
|
|
e@0
|
107 elif (xsig[i] >= Sum):
|
|
e@0
|
108 bit = 1
|
|
e@0
|
109 delta2 = dmin
|
|
e@0
|
110 Sum += delta1
|
|
e@0
|
111 delta1 *= mult1
|
|
e@0
|
112 if delta1 > dmax:
|
|
e@0
|
113 delta1 = dmax
|
|
e@0
|
114
|
|
e@0
|
115
|
|
e@0
|
116 else:
|
|
e@0
|
117 bit = -1
|
|
e@0
|
118 delta1 = dmin
|
|
e@0
|
119 Sum -= delta2
|
|
e@0
|
120 delta2 *= mult2
|
|
e@0
|
121 if delta2 > dmax:
|
|
e@0
|
122 delta2 = dmax
|
|
e@0
|
123
|
|
e@0
|
124
|
|
e@0
|
125
|
|
e@0
|
126 ADMout[0, i] = Sum
|
|
e@0
|
127 codevec[0, i]= bit
|
|
e@0
|
128
|
|
e@0
|
129 return ADMout,codevec, x[0]
|
|
e@0
|
130
|
|
e@0
|
131 if __name__=="__main__":
|
|
e@0
|
132 if len(argv) != 2:
|
|
e@0
|
133 print "[EE] Wrong number of arguments"
|
|
e@0
|
134 print "[II] Correct syntax is:"
|
|
e@0
|
135 print "[II] \t%s <training_file>"
|
|
e@0
|
136 print "[II] where <training_file> is a .yaml file containing the"
|
|
e@0
|
137 print "[II] features of the dataset (try output2_stage/fulltraining-last.yaml)"
|
|
e@0
|
138 exit(-1)
|
|
e@0
|
139
|
|
e@0
|
140
|
|
e@0
|
141 n_clusters = 25
|
|
e@0
|
142 UpsamplingFactor = 10
|
|
e@0
|
143 dmin = 0.001
|
|
e@0
|
144 dmax = 0.28
|
|
e@0
|
145 tol = 0.001
|
|
e@0
|
146
|
|
e@0
|
147 infile = argv[1]
|
|
e@0
|
148
|
|
e@0
|
149 features_pool = YamlInput(filename = infile)()
|
|
e@0
|
150
|
|
e@0
|
151
|
|
e@0
|
152
|
|
e@0
|
153 feature_captions = features_pool.descriptorNames()
|
|
e@0
|
154 parameter_captions = []
|
|
e@0
|
155
|
|
e@0
|
156
|
|
e@0
|
157 for c in features_pool.descriptorNames():
|
|
e@0
|
158 if c.split('.')[0] == 'parameter':
|
|
e@0
|
159 parameter_captions.append(c)
|
|
e@0
|
160 if c.split('.')[0] == 'metadata' or c.split('.')[0] == 'parameter':
|
|
e@0
|
161 feature_captions.remove(c)
|
|
e@0
|
162
|
|
e@0
|
163
|
|
e@0
|
164
|
|
e@0
|
165 # close('all')
|
|
e@0
|
166
|
|
e@0
|
167 print "[II] Loaded training data from %s (%s) " % (infile, features_pool['metadata.date'][0])
|
|
e@0
|
168 print "[II] %d Features Available: " % len(feature_captions)
|
|
e@0
|
169
|
|
e@0
|
170
|
|
e@0
|
171
|
|
e@0
|
172 print str(feature_captions).replace("', ","\n").replace('[','').replace("'","[II]\t ")[:-7]
|
|
e@0
|
173
|
|
e@0
|
174 nfeatures_in = len(feature_captions)
|
|
e@0
|
175 nparameters_in = len(parameter_captions)
|
|
e@0
|
176 features_vector = zeros((nfeatures_in, len(features_pool[feature_captions[0]])))
|
|
e@0
|
177
|
|
e@0
|
178 parameters_vector = zeros((nparameters_in, len(features_pool[parameter_captions[0]])))
|
|
e@0
|
179
|
|
e@0
|
180
|
|
e@0
|
181 for i in range(0, nfeatures_in):
|
|
e@0
|
182 features_vector[i, :] = features_pool[feature_captions[i]].T
|
|
e@0
|
183 for i in range(0, nparameters_in):
|
|
e@0
|
184 parameters_vector[i, :] = features_pool[parameter_captions[0]].T
|
|
e@0
|
185
|
|
e@0
|
186 print "[II] %d parameters used:" % len(parameter_captions)
|
|
e@0
|
187 print str(parameter_captions).replace("', ","\n").replace('[','').replace("'","[II]\t ")[:-7].replace('parameter.','')
|
|
e@0
|
188
|
|
e@0
|
189 print "[II] Marking silent parts"
|
|
e@0
|
190
|
|
e@0
|
191 silent_parts = zeros((1, len(features_pool[feature_captions[i]].T)))
|
|
e@0
|
192
|
|
e@0
|
193 rms = features_vector[feature_captions.index('rms'), :]
|
|
e@0
|
194
|
|
e@0
|
195 # Implementing Hysteresis Gate -- High threshold is halfway between
|
|
e@0
|
196 # the mean and the max and Low is halfway between the mean dn the min
|
|
e@0
|
197
|
|
e@0
|
198 rms_threshold_mean = mean(rms)
|
|
e@0
|
199
|
|
e@0
|
200 rms_threshold_max = max(rms)
|
|
e@0
|
201 rms_threshold_min = min(rms)
|
|
e@0
|
202
|
|
e@0
|
203 rms_threshold_high = 0.1 * rms_threshold_mean
|
|
e@0
|
204 rms_threshold_low = 0.01 * rms_threshold_mean
|
|
e@0
|
205
|
|
e@0
|
206 for n in range(1, len(rms)):
|
|
e@0
|
207 prev = rms[n-1]
|
|
e@0
|
208 curr = rms[n]
|
|
e@0
|
209
|
|
e@0
|
210 if prev >= rms_threshold_high:
|
|
e@0
|
211 if curr < rms_threshold_low:
|
|
e@0
|
212 silent_parts[0,n] = 1
|
|
e@0
|
213 else:
|
|
e@0
|
214 silent_parts[0,n] = 0
|
|
e@0
|
215 elif prev <= rms_threshold_low:
|
|
e@0
|
216 if curr > rms_threshold_high:
|
|
e@0
|
217 silent_parts[0,n] = 0
|
|
e@0
|
218 else:
|
|
e@0
|
219 silent_parts[0,n] = 1
|
|
e@0
|
220 else:
|
|
e@0
|
221 silent_parts[0,n] = silent_parts[0,n-1]
|
|
e@0
|
222
|
|
e@0
|
223
|
|
e@0
|
224 if silent_parts[0,1] == 1:
|
|
e@0
|
225 silent_parts[0, 0] = 1
|
|
e@0
|
226
|
|
e@0
|
227
|
|
e@0
|
228
|
|
e@0
|
229 active_index = invert(silent_parts.flatten().astype(bool))
|
|
e@0
|
230
|
|
e@0
|
231 # Keep only active parts
|
|
e@0
|
232
|
|
e@0
|
233 # Uncomment this
|
|
e@0
|
234 features_vector = features_vector[:, active_index]
|
|
e@0
|
235
|
|
e@0
|
236 moments_vector = zeros((features_vector.shape[0], 2))
|
|
e@0
|
237
|
|
e@0
|
238 print "[II] Storing moments vector"
|
|
e@0
|
239 for i in range(0, features_vector.shape[0]):
|
|
e@0
|
240 mean_ = mean(features_vector[i,:])
|
|
e@0
|
241 std_ = std(features_vector[i,:], ddof=1)
|
|
e@0
|
242 moments_vector[i,0] = mean_
|
|
e@0
|
243 moments_vector[i,1] = std_
|
|
e@0
|
244
|
|
e@0
|
245 features_vector[i,:] = (features_vector[i,:] - mean_)/std_
|
|
e@0
|
246
|
|
e@0
|
247 features_vector_original = features_vector
|
|
e@0
|
248
|
|
e@0
|
249
|
|
e@0
|
250 print "[II] Extracting PCA configuration "
|
|
e@0
|
251
|
|
e@0
|
252 kernel, q, featurelist = extract_pca_configuration_from_data(features_vector)
|
|
e@0
|
253
|
|
e@0
|
254 print "[II] Optimal number of PCs to keep: %d" % q
|
|
e@0
|
255
|
|
e@0
|
256 feature_captions_array = array(feature_captions)
|
|
e@0
|
257
|
|
e@0
|
258 features_to_keep = list(feature_captions_array[featurelist])
|
|
e@0
|
259 print "[II] Decided to keep %d features:" % len(features_to_keep)
|
|
e@0
|
260 print str(features_to_keep).replace("', ","\n").replace('[','').replace("'","[II]\t ")[:-7]
|
|
e@0
|
261
|
|
e@0
|
262
|
|
e@0
|
263 features_kept_data = features_vector[featurelist,:]
|
|
e@0
|
264
|
|
e@0
|
265 features_vector = (kernel.T*features_kept_data)[0:q,:]
|
|
e@0
|
266
|
|
e@0
|
267 parameters_k_means = KMeans(init='k-means++', n_init=10, max_iter=300, tol=0.0000001, verbose = 0)
|
|
e@0
|
268
|
|
e@0
|
269 print "[II] Trying ADM-coded parameters"
|
|
e@0
|
270 print "[II] Upsampling features and parameters by a factor of %d" % UpsamplingFactor
|
|
e@0
|
271
|
|
e@0
|
272
|
|
e@0
|
273 # Upsampled features and parameters
|
|
e@0
|
274 features_vector_upsampled = smooth_matrix_1D(repeat(features_vector,UpsamplingFactor, axis=1))
|
|
e@0
|
275
|
|
e@0
|
276 # feature_labels_upsampled = repeat(features_clustering_labels,UpsamplingFactor, axis=0)
|
|
e@0
|
277 parameters_vector_upsampled = repeat(parameters_vector,UpsamplingFactor, axis=1)
|
|
e@0
|
278
|
|
e@0
|
279 # parameters_vector_upsampled = smooth_matrix_1D(parameters_vector_upsampled)
|
|
e@0
|
280
|
|
e@0
|
281 parameters_vector_upsampled_adm = matrix(zeros(shape(parameters_vector_upsampled)))
|
|
e@0
|
282 parameters_vector_upsampled_code = matrix(zeros(shape(parameters_vector_upsampled)))
|
|
e@0
|
283 parameters_vector_upsampled_firstval = matrix(zeros((parameters_vector_upsampled.shape[0],1)))
|
|
e@0
|
284
|
|
e@0
|
285 # Reconstructed parameters
|
|
e@0
|
286
|
|
e@0
|
287 parameters_vector_upsampled_reconstructed = matrix(zeros(shape(parameters_vector_upsampled)))
|
|
e@0
|
288
|
|
e@0
|
289
|
|
e@0
|
290
|
|
e@0
|
291
|
|
e@0
|
292 def adm_matrix(X, dmin=0.001,dmax=0.28,tol=0.001):
|
|
e@0
|
293
|
|
e@0
|
294 out = matrix(zeros(shape(X)))
|
|
e@0
|
295 code = matrix(zeros(shape(X)))
|
|
e@0
|
296 firstval = matrix(zeros((X.shape[0], 1)))
|
|
e@0
|
297
|
|
e@0
|
298 for i in range(0, X.shape[0]):
|
|
e@0
|
299 out[i,:], code[i,:], firstval[i,0] = adm(X[i,:],dmin=dmin,dmax=dmax,tol=tol)
|
|
e@0
|
300
|
|
e@0
|
301 return out,code,firstval
|
|
e@0
|
302
|
|
e@0
|
303 # parameters_vector_upsampled_reconstructed[i,:] = adm_reconstruct(parameters_vector_upsampled_code[i,:],parameters_vector_upsampled_firstval[i,0], dmin=dmin,dmax=dmax)
|
|
e@0
|
304
|
|
e@0
|
305 def adm_matrix_reconstruct(code, firstval, dmin=0.001, dmax=0.28):
|
|
e@0
|
306 X = matrix(zeros(shape(code)))
|
|
e@0
|
307 for i in range(0, code.shape[0]):
|
|
e@0
|
308 X[i,:] = adm_reconstruct(code[i,:], firstval[i,0], dmin=dmin, dmax=dmax)
|
|
e@0
|
309
|
|
e@0
|
310 return X
|
|
e@0
|
311
|
|
e@0
|
312
|
|
e@0
|
313 parameters_vector_upsampled_adm, parameters_vector_upsampled_code, parameters_vector_upsampled_firstval = adm_matrix(parameters_vector_upsampled, dmin, dmax, tol)
|
|
e@0
|
314
|
|
e@0
|
315
|
|
e@0
|
316 def diff_and_pad(X):
|
|
e@0
|
317 return concatenate((
|
|
e@0
|
318 zeros((
|
|
e@0
|
319 shape(X)[0],
|
|
e@0
|
320 1
|
|
e@0
|
321 )),
|
|
e@0
|
322 diff(X, axis=1)),
|
|
e@0
|
323 axis=1)
|
|
e@0
|
324
|
|
e@0
|
325
|
|
e@0
|
326 print "[II] Clustering features."
|
|
e@0
|
327 #
|
|
e@0
|
328 features_clustering = GMM(n_components = n_clusters, covariance_type='diag')
|
|
e@0
|
329 #
|
|
e@0
|
330 features_clustering.fit( features_vector_upsampled.T, y=parameters_vector_upsampled_code)
|
|
e@0
|
331 #
|
|
e@0
|
332 features_clustering_means = features_clustering.means_
|
|
e@0
|
333 features_clustering_labels = features_clustering.predict(features_vector_upsampled.T)
|
|
e@0
|
334 features_clustering_sigmas = features_clustering.covars_
|
|
e@0
|
335 #
|
|
e@0
|
336 features_vector_upsampled_estimated = zeros(shape(features_vector_upsampled))
|
|
e@0
|
337 #
|
|
e@0
|
338 #
|
|
e@0
|
339 for n in range(0, len(features_vector_upsampled_estimated[0])):
|
|
e@0
|
340 features_vector_upsampled_estimated[:,n] = features_clustering_means[features_clustering_labels[n]]
|
|
e@0
|
341 #
|
|
e@0
|
342 #
|
|
e@0
|
343 print "[II] Features MSE for %d clusters: %.3f" % (n_clusters, mse(features_vector_upsampled, features_vector_upsampled_estimated))
|
|
e@0
|
344
|
|
e@0
|
345
|
|
e@0
|
346
|
|
e@0
|
347
|
|
e@0
|
348 def cross_validate_classification(data, classes, estimator):
|
|
e@0
|
349 print "[II] Crossvalidating... "
|
|
e@0
|
350 from copy import deepcopy
|
|
e@0
|
351 estimator_fulldata = deepcopy(estimator)
|
|
e@0
|
352 estimator_fulldata.fit(data, classes)
|
|
e@0
|
353
|
|
e@0
|
354 percents = arange(0.1,0.9,0.1)
|
|
e@0
|
355 MSEs = []
|
|
e@0
|
356 labels = estimator.predict(data)
|
|
e@0
|
357
|
|
e@0
|
358 print "[II] for full training-testing: %.2f" % (sum(array(classes==labels).astype(float))/len(labels))
|
|
e@0
|
359
|
|
e@0
|
360 for p in percents:
|
|
e@0
|
361 train,test,trainlabels,testlabels = cross_validation.train_test_split(data,classes,test_size=p,random_state=0)
|
|
e@0
|
362 estimator_ = deepcopy(estimator)
|
|
e@0
|
363 estimator_.fit(train, trainlabels)
|
|
e@0
|
364 labels = estimator.predict(test)
|
|
e@0
|
365 print "[II] for training(%.2f)-testing(%.2f): %.2f" % ((1-p),p,sum(array(testlabels==labels).astype(float))/len(labels))
|
|
e@0
|
366
|
|
e@0
|
367 return MSEs
|
|
e@0
|
368
|
|
e@0
|
369 def cross_validate_clustering(data, estimator):
|
|
e@0
|
370 print "[II] Crossvalidating... "
|
|
e@0
|
371 estimator_fulldata = estimator
|
|
e@0
|
372 estimator_fulldata.fit(data)
|
|
e@0
|
373
|
|
e@0
|
374 # labels = estimator_fulldata.predict(data)
|
|
e@0
|
375 means = estimator_fulldata.means_
|
|
e@0
|
376 # print means
|
|
e@0
|
377
|
|
e@0
|
378 percents = arange(0.1,0.6,0.1)
|
|
e@0
|
379 MSEs = []
|
|
e@0
|
380 reconstructed = zeros(shape(data))
|
|
e@0
|
381 labels = estimator.predict(data)
|
|
e@0
|
382 for n in range(0, len(reconstructed)):
|
|
e@0
|
383 reconstructed[n,:] = means[labels[n]]
|
|
e@0
|
384
|
|
e@0
|
385 MSEs.append(mse(data,reconstructed))
|
|
e@0
|
386 for p in percents:
|
|
e@0
|
387 train,test = cross_validation.train_test_split(data,test_size=p,random_state=0)
|
|
e@0
|
388 train = matrix(train)
|
|
e@0
|
389 test = matrix(test)
|
|
e@0
|
390
|
|
e@0
|
391 estimator.fit(train)
|
|
e@0
|
392 means = estimator.means_
|
|
e@0
|
393 labels = estimator.predict(test)
|
|
e@0
|
394 reconstructed = zeros(shape(test))
|
|
e@0
|
395 for n in range(0, len(reconstructed)):
|
|
e@0
|
396 reconstructed[n,:] = means[labels[n]]
|
|
e@0
|
397
|
|
e@0
|
398 m = mse(test,reconstructed)
|
|
e@0
|
399
|
|
e@0
|
400 print "[II] MSE for clustering crossvalidated data %.2f-%.2f: %.5f" % ((1-p), p, m)
|
|
e@0
|
401 MSEs.append(m)
|
|
e@0
|
402
|
|
e@0
|
403 print "[II] Crossvalidation complete"
|
|
e@0
|
404
|
|
e@0
|
405 return MSEs
|
|
e@0
|
406
|
|
e@0
|
407
|
|
e@0
|
408
|
|
e@0
|
409
|
|
e@0
|
410 # Construct parameters alphabet, each symbol is going to be a different column vector
|
|
e@0
|
411 # in parameter code matrix
|
|
e@0
|
412
|
|
e@0
|
413
|
|
e@0
|
414 def vector_to_states(X):
|
|
e@0
|
415 """
|
|
e@0
|
416 Input: a vector MxN with N samples and M variables
|
|
e@0
|
417 Output: a codeword dictionary `parameters_alphabet',
|
|
e@0
|
418 state_seq, inverse `parameters_alphabet_inv' """
|
|
e@0
|
419
|
|
e@0
|
420
|
|
e@0
|
421 parameters_alphabet = {}
|
|
e@0
|
422 n = 0
|
|
e@0
|
423
|
|
e@0
|
424 for i in range(0, X.shape[1]):
|
|
e@0
|
425 vec = tuple(ravel(X[:,i]))
|
|
e@0
|
426 if vec not in parameters_alphabet:
|
|
e@0
|
427 parameters_alphabet[vec] = n
|
|
e@0
|
428 n += 1
|
|
e@0
|
429
|
|
e@0
|
430 parameters_alphabet_inv = dict([(parameters_alphabet[m],m) for m in parameters_alphabet])
|
|
e@0
|
431
|
|
e@0
|
432 state_seq = array([parameters_alphabet[tuple(ravel(X[:,m]))] for m in range(0, parameters_vector_upsampled_code.shape[1])] )
|
|
e@0
|
433
|
|
e@0
|
434 return state_seq, parameters_alphabet, parameters_alphabet_inv
|
|
e@0
|
435
|
|
e@0
|
436
|
|
e@0
|
437 def states_to_vector(predicted, parameters_alphabet_inv):
|
|
e@0
|
438 estimated = matrix(zeros((len(parameters_alphabet_inv[0]), len(predicted))))
|
|
e@0
|
439 for i in range(0, len(state_seq)):
|
|
e@0
|
440 estimated[:, i] = matrix(parameters_alphabet_inv[predicted[i]]).T
|
|
e@0
|
441
|
|
e@0
|
442 return estimated
|
|
e@0
|
443
|
|
e@0
|
444 state_seq, parameters_alphabet, parameters_alphabet_inv = vector_to_states(parameters_vector_upsampled_code)
|
|
e@0
|
445
|
|
e@0
|
446
|
|
e@0
|
447 parameters_change_variable = matrix(diff_and_pad(parameters_vector_upsampled)!=0).astype(int)
|
|
e@0
|
448
|
|
e@0
|
449 changes_state_seq, changes_parameters_alphabet, changes_parameters_alphabet_inv = vector_to_states(parameters_change_variable)
|
|
e@0
|
450
|
|
e@0
|
451
|
|
e@0
|
452 # This is an hmm that just codes the changes"
|
|
e@0
|
453 # We have only two states, change and stay the same.
|
|
e@0
|
454
|
|
e@0
|
455
|
|
e@0
|
456 parameters_state, parameter_state_alphabet, parameter_state_alphabet_inv = vector_to_states(parameters_vector_upsampled)
|
|
e@0
|
457
|
|
e@0
|
458
|
|
e@0
|
459 print "[II] Testing Gaussian Naive Bayes Classifier"
|
|
e@0
|
460 gnb = GaussianNB()
|
|
e@0
|
461 gnb.fit(features_vector_upsampled.T, parameters_state)
|
|
e@0
|
462
|
|
e@0
|
463 parameters_state_estimated = gnb.predict(features_vector_upsampled.T)
|
|
e@0
|
464
|
|
e@0
|
465 output = states_to_vector(parameters_state_estimated, parameter_state_alphabet_inv)
|
|
e@0
|
466
|
|
e@0
|
467 figure()
|
|
e@0
|
468 subplot(211)
|
|
e@0
|
469 plot(parameters_vector_upsampled.T)
|
|
e@0
|
470 title('Parameter value upsampled by a factor of %d' % UpsamplingFactor)
|
|
e@0
|
471 ylabel('value')
|
|
e@0
|
472 xlabel('frame #')
|
|
e@0
|
473 subplot(212)
|
|
e@0
|
474 #plot(smooth_matrix_1D(output).T)
|
|
e@0
|
475 plot(output.T)
|
|
e@0
|
476 ylabel('value')
|
|
e@0
|
477 xlabel('frame #')
|
|
e@0
|
478 cross_validate_classification(features_vector_upsampled.T, parameters_state, gnb)
|
|
e@0
|
479
|
|
e@0
|
480 print "[II] Trying Multinomial HMM"
|
|
e@0
|
481
|
|
e@0
|
482 # In order to do classification with HMMs, we need to:
|
|
e@0
|
483 # 1. Split the parameters into classes
|
|
e@0
|
484 # 2. Train one model per class
|
|
e@0
|
485 # 3. Feed our data to all the models
|
|
e@0
|
486 # 4. Check which has a better score and assig,n to M
|
|
e@0
|
487
|
|
e@0
|
488
|
|
e@0
|
489 class HmmClassifier:
|
|
e@0
|
490 def __init__(self, N=2, n_components = 1):
|
|
e@0
|
491 self.n_components = n_components
|
|
e@0
|
492 self.chain_size = N
|
|
e@0
|
493 self.hmms_ = []
|
|
e@0
|
494 self.N = N
|
|
e@0
|
495
|
|
e@0
|
496 def fit(self, X, states):
|
|
e@0
|
497 self.n_states = len(unique(states))
|
|
e@0
|
498
|
|
e@0
|
499 for n in range(0, self.n_states):
|
|
e@0
|
500 hmm_ = hmm.GaussianHMM(n_components = self.n_components, covariance_type = 'full')
|
|
e@0
|
501
|
|
e@0
|
502 # Get training data for each class
|
|
e@0
|
503 vector = X[states == n,:]
|
|
e@0
|
504
|
|
e@0
|
505 # Fit the HMM
|
|
e@0
|
506 # print vector
|
|
e@0
|
507 hmm_.fit([vector])
|
|
e@0
|
508
|
|
e@0
|
509 # And append to the list
|
|
e@0
|
510 self.hmms_.append(hmm_)
|
|
e@0
|
511
|
|
e@0
|
512 def predict(self,X):
|
|
e@0
|
513 labels = zeros((X.shape[0],))
|
|
e@0
|
514 N = self.N
|
|
e@0
|
515
|
|
e@0
|
516 m = 0
|
|
e@0
|
517
|
|
e@0
|
518 scores = zeros((1, self.n_states))
|
|
e@0
|
519
|
|
e@0
|
520
|
|
e@0
|
521 while m*N < X.shape[0]:
|
|
e@0
|
522 if (m+1)*N > X.shape[0]:
|
|
e@0
|
523 testdata = X[m*N:,:]
|
|
e@0
|
524 else:
|
|
e@0
|
525 testdata = X[m*N:(m+1)*N,:]
|
|
e@0
|
526
|
|
e@0
|
527 # print testdata
|
|
e@0
|
528
|
|
e@0
|
529 for i in range(0, self.n_states):
|
|
e@0
|
530 scores[0,i] = self.hmms_[i].score(testdata)
|
|
e@0
|
531
|
|
e@0
|
532 if (m+1)*N > X.shape[0]:
|
|
e@0
|
533 labels[m*N:] = argmax(scores)
|
|
e@0
|
534 else:
|
|
e@0
|
535 labels[m*N:(m+1)*N] = argmax(scores)
|
|
e@0
|
536
|
|
e@0
|
537 m+= 1
|
|
e@0
|
538
|
|
e@0
|
539 return labels
|
|
e@0
|
540
|
|
e@0
|
541 N = 150
|
|
e@0
|
542 n_components = 1
|
|
e@0
|
543
|
|
e@0
|
544 hmmc = HmmClassifier(N = N, n_components = n_components)
|
|
e@0
|
545 hmmc.fit(features_vector_upsampled.T, parameters_state)
|
|
e@0
|
546
|
|
e@0
|
547 cross_validate_classification(features_vector_upsampled.T, parameters_state, hmmc)
|
|
e@0
|
548
|
|
e@0
|
549
|
|
e@0
|
550
|
|
e@0
|
551
|
|
e@0
|
552 # hmms_ = []
|
|
e@0
|
553 #
|
|
e@0
|
554 # for n in range(0, len(parameter_state_alphabet)):
|
|
e@0
|
555 # #hmm_ = hmm.GMMHMM(n_components = 1, n_mix = 2)
|
|
e@0
|
556 # hmm_ = hmm.GaussianHMM(n_components = 1,covariance_type = 'full')
|
|
e@0
|
557 #
|
|
e@0
|
558 # # Get training data for each class
|
|
e@0
|
559 # vector = features_vector_upsampled[:,parameters_state == n]
|
|
e@0
|
560 #
|
|
e@0
|
561 # #if vector.shape[1] < n_clusters:
|
|
e@0
|
562 # # hmms_.append(None)
|
|
e@0
|
563 # #else:
|
|
e@0
|
564 #
|
|
e@0
|
565 # hmm_.fit([vector.T])
|
|
e@0
|
566 #
|
|
e@0
|
567 # # Append to the list
|
|
e@0
|
568 #
|
|
e@0
|
569 # hmms_.append(hmm_)
|
|
e@0
|
570 #
|
|
e@0
|
571 # labels = zeros((features_vector_upsampled.shape[1],))
|
|
e@0
|
572 #
|
|
e@0
|
573 # N = 20
|
|
e@0
|
574 # m = 0
|
|
e@0
|
575 #
|
|
e@0
|
576 # while m*N < features_vector_upsampled.shape[1]:
|
|
e@0
|
577 #
|
|
e@0
|
578 # scores = zeros((1, len(parameter_state_alphabet)))
|
|
e@0
|
579 #
|
|
e@0
|
580 # if (m+1)*N > features_vector_upsampled.shape[1]:
|
|
e@0
|
581 # testdata = features_vector_upsampled[:,m*N:]
|
|
e@0
|
582 # else:
|
|
e@0
|
583 # testdata = features_vector_upsampled[:,m*N:(m+1)*N]
|
|
e@0
|
584 #
|
|
e@0
|
585 # for i in range(0, len(parameter_state_alphabet)):
|
|
e@0
|
586 # if hmms_[i] is not None:
|
|
e@0
|
587 # scores[0,i] = hmms_[i].score(testdata.T)
|
|
e@0
|
588 # else:
|
|
e@0
|
589 # scores[0,i] = -100000 # Very large negative score
|
|
e@0
|
590 # if (m+1)*N >= features_vector_upsampled.shape[1]:
|
|
e@0
|
591 # labels[m*N:] = argmax(scores)
|
|
e@0
|
592 # else:
|
|
e@0
|
593 # labels[m*N:(m+1)*N] = argmax(scores)
|
|
e@0
|
594 #
|
|
e@0
|
595 # m += 1
|
|
e@0
|
596
|
|
e@0
|
597
|
|
e@0
|
598 # figure()
|
|
e@0
|
599 #plot(labels.T)
|
|
e@0
|
600
|
|
e@0
|
601
|
|
e@0
|
602 labels = hmmc.predict(features_vector_upsampled.T)
|
|
e@0
|
603 estimated = states_to_vector(labels,parameter_state_alphabet_inv)
|
|
e@0
|
604 plot(estimated.T,'r--')
|
|
e@0
|
605
|
|
e@0
|
606 title('Estimated parameter values')
|
|
e@0
|
607 legend(['Naive Bayes Classifier', 'HMM chain size %d (%.1fms)' % (N, float(N)/UpsamplingFactor*23.0)])
|
|
e@0
|
608
|
|
e@0
|
609 ylabel('value')
|
|
e@0
|
610 xlabel('frame #') |
