SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GMM.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011 Alesis Novik
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 #include <shogun/lib/config.h>
11 
12 #ifdef HAVE_LAPACK
13 
14 #include <shogun/clustering/GMM.h>
17 #include <shogun/base/Parameter.h>
21 #include <shogun/multiclass/KNN.h>
22 
23 #include <vector>
24 
25 using namespace shogun;
26 using namespace std;
27 
28 CGMM::CGMM() : CDistribution(), m_components(), m_coefficients()
29 {
30  register_params();
31 }
32 
33 CGMM::CGMM(int32_t n, ECovType cov_type) : CDistribution(), m_components(), m_coefficients()
34 {
35  m_coefficients.vector=SG_MALLOC(float64_t, n);
37  m_components = vector<CGaussian*>(n);
38 
39  for (int32_t i=0; i<n; i++)
40  {
41  m_components[i]=new CGaussian();
42  SG_REF(m_components[i]);
43  m_components[i]->set_cov_type(cov_type);
44  }
45 
46  register_params();
47 }
48 
49 CGMM::CGMM(vector<CGaussian*> components, SGVector<float64_t> coefficients, bool copy) : CDistribution()
50 {
51  ASSERT(int32_t(components.size())==coefficients.vlen)
52 
53  if (!copy)
54  {
55  m_components=components;
56  m_coefficients=coefficients;
57  for (int32_t i=0; i<int32_t(components.size()); i++)
58  {
59  SG_REF(m_components[i]);
60  }
61  }
62  else
63  {
64  m_coefficients = coefficients;
65  m_components = vector<CGaussian*>(components.size());
66 
67  for (int32_t i=0; i<int32_t(components.size()); i++)
68  {
69  m_components[i]=new CGaussian();
70  SG_REF(m_components[i]);
71  m_components[i]->set_cov_type(components[i]->get_cov_type());
72 
73  SGVector<float64_t> old_mean=components[i]->get_mean();
74  SGVector<float64_t> new_mean(old_mean.vlen);
75  memcpy(new_mean.vector, old_mean.vector, old_mean.vlen*sizeof(float64_t));
76  m_components[i]->set_mean(new_mean);
77 
78  SGVector<float64_t> old_d=components[i]->get_d();
79  SGVector<float64_t> new_d(old_d.vlen);
80  memcpy(new_d.vector, old_d.vector, old_d.vlen*sizeof(float64_t));
81  m_components[i]->set_d(new_d);
82 
83  if (components[i]->get_cov_type()==FULL)
84  {
85  SGMatrix<float64_t> old_u=components[i]->get_u();
86  SGMatrix<float64_t> new_u(old_u.num_rows, old_u.num_cols);
87  memcpy(new_u.matrix, old_u.matrix, old_u.num_rows*old_u.num_cols*sizeof(float64_t));
88  m_components[i]->set_u(new_u);
89  }
90 
91  m_coefficients[i]=coefficients[i];
92  }
93  }
94 
95  register_params();
96 }
97 
99 {
100  if (!m_components.empty())
101  cleanup();
102 }
103 
105 {
106  for (int32_t i = 0; i < int32_t(m_components.size()); i++)
108 
109  m_components = vector<CGaussian*>();
111 }
112 
114 {
115  ASSERT(m_components.size() != 0)
116 
118  if (data)
119  {
120  if (!data->has_property(FP_DOT))
121  SG_ERROR("Specified features are not of type CDotFeatures\n")
122  set_features(data);
123  }
124 
125  return true;
126 }
127 
128 float64_t CGMM::train_em(float64_t min_cov, int32_t max_iter, float64_t min_change)
129 {
130  if (!features)
131  SG_ERROR("No features to train on.\n")
132 
133  CDotFeatures* dotdata=(CDotFeatures *) features;
134  int32_t num_vectors=dotdata->get_num_vectors();
135 
136  SGMatrix<float64_t> alpha;
137 
138  /* compute initialization via kmeans if none is present */
139  if (m_components[0]->get_mean().vector==NULL)
140  {
141  CKMeans* init_k_means=new CKMeans(int32_t(m_components.size()), new CEuclideanDistance());
142  init_k_means->train(dotdata);
143  SGMatrix<float64_t> init_means=init_k_means->get_cluster_centers();
144 
145  alpha=alpha_init(init_means);
146 
147  SG_UNREF(init_k_means);
148 
149  max_likelihood(alpha, min_cov);
150  }
151  else
152  alpha=SGMatrix<float64_t>(num_vectors,int32_t(m_components.size()));
153 
154  int32_t iter=0;
155  float64_t log_likelihood_prev=0;
156  float64_t log_likelihood_cur=0;
157  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size());
158  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
159  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
160 
161  while (iter<max_iter)
162  {
163  log_likelihood_prev=log_likelihood_cur;
164  log_likelihood_cur=0;
165 
166  for (int32_t i=0; i<num_vectors; i++)
167  {
168  logPx[i]=0;
170  for (int32_t j=0; j<int32_t(m_components.size()); j++)
171  {
172  logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]);
173  logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]);
174  }
175 
176  logPx[i]=CMath::log(logPx[i]);
177  log_likelihood_cur+=logPx[i];
178 
179  for (int32_t j=0; j<int32_t(m_components.size()); j++)
180  {
181  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
182  alpha.matrix[i*m_components.size()+j]=CMath::exp(logPxy[i*m_components.size()+j]-logPx[i]);
183  }
184  }
185 
186  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
187  break;
188 
189  max_likelihood(alpha, min_cov);
190 
191  iter++;
192  }
193 
194  SG_FREE(logPxy);
195  SG_FREE(logPx);
196 
197  return log_likelihood_cur;
198 }
199 
200 float64_t CGMM::train_smem(int32_t max_iter, int32_t max_cand, float64_t min_cov, int32_t max_em_iter, float64_t min_change)
201 {
202  if (!features)
203  SG_ERROR("No features to train on.\n")
204 
205  if (m_components.size()<3)
206  SG_ERROR("Can't run SMEM with less than 3 component mixture model.\n")
207 
208  CDotFeatures* dotdata=(CDotFeatures *) features;
209  int32_t num_vectors=dotdata->get_num_vectors();
210 
211  float64_t cur_likelihood=train_em(min_cov, max_em_iter, min_change);
212 
213  int32_t iter=0;
214  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size());
215  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
216  float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.size());
217  float64_t* logPostSum=SG_MALLOC(float64_t, m_components.size());
218  float64_t* logPostSum2=SG_MALLOC(float64_t, m_components.size());
219  float64_t* logPostSumSum=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2);
220  float64_t* split_crit=SG_MALLOC(float64_t, m_components.size());
221  float64_t* merge_crit=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2);
222  int32_t* split_ind=SG_MALLOC(int32_t, m_components.size());
223  int32_t* merge_ind=SG_MALLOC(int32_t, m_components.size()*(m_components.size()-1)/2);
224 
225  while (iter<max_iter)
226  {
227  memset(logPostSum, 0, m_components.size()*sizeof(float64_t));
228  memset(logPostSum2, 0, m_components.size()*sizeof(float64_t));
229  memset(logPostSumSum, 0, (m_components.size()*(m_components.size()-1)/2)*sizeof(float64_t));
230  for (int32_t i=0; i<num_vectors; i++)
231  {
232  logPx[i]=0;
234  for (int32_t j=0; j<int32_t(m_components.size()); j++)
235  {
236  logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]);
237  logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]);
238  }
239 
240  logPx[i]=CMath::log(logPx[i]);
241 
242  for (int32_t j=0; j<int32_t(m_components.size()); j++)
243  {
244  logPost[i*m_components.size()+j]=logPxy[i*m_components.size()+j]-logPx[i];
245  logPostSum[j]+=CMath::exp(logPost[i*m_components.size()+j]);
246  logPostSum2[j]+=CMath::exp(2*logPost[i*m_components.size()+j]);
247  }
248 
249  int32_t counter=0;
250  for (int32_t j=0; j<int32_t(m_components.size()); j++)
251  {
252  for (int32_t k=j+1; k<int32_t(m_components.size()); k++)
253  {
254  logPostSumSum[counter]+=CMath::exp(logPost[i*m_components.size()+j]+logPost[i*m_components.size()+k]);
255  counter++;
256  }
257  }
258  }
259 
260  int32_t counter=0;
261  for (int32_t i=0; i<int32_t(m_components.size()); i++)
262  {
263  logPostSum[i]=CMath::log(logPostSum[i]);
264  split_crit[i]=0;
265  split_ind[i]=i;
266  for (int32_t j=0; j<num_vectors; j++)
267  {
268  split_crit[i]+=(logPost[j*m_components.size()+i]-logPostSum[i]-logPxy[j*m_components.size()+i]+CMath::log(m_coefficients[i]))*
269  (CMath::exp(logPost[j*m_components.size()+i])/CMath::exp(logPostSum[i]));
270  }
271  for (int32_t j=i+1; j<int32_t(m_components.size()); j++)
272  {
273  merge_crit[counter]=CMath::log(logPostSumSum[counter])-(0.5*CMath::log(logPostSum2[i]))-(0.5*CMath::log(logPostSum2[j]));
274  merge_ind[counter]=i*m_components.size()+j;
275  counter++;
276  }
277  }
278  CMath::qsort_backward_index(split_crit, split_ind, int32_t(m_components.size()));
279  CMath::qsort_backward_index(merge_crit, merge_ind, int32_t(m_components.size()*(m_components.size()-1)/2));
280 
281  bool better_found=false;
282  int32_t candidates_checked=0;
283  for (int32_t i=0; i<int32_t(m_components.size()); i++)
284  {
285  for (int32_t j=0; j<int32_t(m_components.size()*(m_components.size()-1)/2); j++)
286  {
287  if (merge_ind[j]/int32_t(m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%m_components.size()) != split_ind[i])
288  {
289  candidates_checked++;
290  CGMM* candidate=new CGMM(m_components, m_coefficients, true);
291  candidate->train(features);
292  candidate->partial_em(split_ind[i], merge_ind[j]/int32_t(m_components.size()), merge_ind[j]%int32_t(m_components.size()), min_cov, max_em_iter, min_change);
293  float64_t cand_likelihood=candidate->train_em(min_cov, max_em_iter, min_change);
294 
295  if (cand_likelihood>cur_likelihood)
296  {
297  cur_likelihood=cand_likelihood;
298  set_comp(candidate->get_comp());
299  set_coef(candidate->get_coef());
300 
301  for (int32_t k=0; k<int32_t(candidate->get_comp().size()); k++)
302  {
303  SG_UNREF(candidate->get_comp()[k]);
304  }
305 
306  better_found=true;
307  break;
308  }
309  else
310  delete candidate;
311 
312  if (candidates_checked>=max_cand)
313  break;
314  }
315  }
316  if (better_found || candidates_checked>=max_cand)
317  break;
318  }
319  if (!better_found)
320  break;
321  iter++;
322  }
323 
324  SG_FREE(logPxy);
325  SG_FREE(logPx);
326  SG_FREE(logPost);
327  SG_FREE(split_crit);
328  SG_FREE(merge_crit);
329  SG_FREE(logPostSum);
330  SG_FREE(logPostSum2);
331  SG_FREE(logPostSumSum);
332  SG_FREE(split_ind);
333  SG_FREE(merge_ind);
334 
335  return cur_likelihood;
336 }
337 
338 void CGMM::partial_em(int32_t comp1, int32_t comp2, int32_t comp3, float64_t min_cov, int32_t max_em_iter, float64_t min_change)
339 {
340  CDotFeatures* dotdata=(CDotFeatures *) features;
341  int32_t num_vectors=dotdata->get_num_vectors();
342 
343  float64_t* init_logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size());
344  float64_t* init_logPx=SG_MALLOC(float64_t, num_vectors);
345  float64_t* init_logPx_fix=SG_MALLOC(float64_t, num_vectors);
346  float64_t* post_add=SG_MALLOC(float64_t, num_vectors);
347 
348  for (int32_t i=0; i<num_vectors; i++)
349  {
350  init_logPx[i]=0;
351  init_logPx_fix[i]=0;
352 
354  for (int32_t j=0; j<int32_t(m_components.size()); j++)
355  {
356  init_logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]);
357  init_logPx[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]);
358  if (j!=comp1 && j!=comp2 && j!=comp3)
359  {
360  init_logPx_fix[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]);
361  }
362  }
363 
364  init_logPx[i]=CMath::log(init_logPx[i]);
365  post_add[i]=CMath::log(CMath::exp(init_logPxy[i*m_components.size()+comp1]-init_logPx[i])+
366  CMath::exp(init_logPxy[i*m_components.size()+comp2]-init_logPx[i])+
367  CMath::exp(init_logPxy[i*m_components.size()+comp3]-init_logPx[i]));
368  }
369 
370  vector<CGaussian*> components(3);
371  SGVector<float64_t> coefficients(3);
372  components[0]=m_components[comp1];
373  components[1]=m_components[comp2];
374  components[2]=m_components[comp3];
375  coefficients.vector[0]=m_coefficients.vector[comp1];
376  coefficients.vector[1]=m_coefficients.vector[comp2];
377  coefficients.vector[2]=m_coefficients.vector[comp3];
378  float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2];
379 
380  int32_t dim_n=components[0]->get_mean().vlen;
381 
382  float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]);
383  float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]);
384 
385  float64_t noise_mag=SGVector<float64_t>::twonorm(components[0]->get_mean().vector, dim_n)*0.1/
386  CMath::sqrt((float64_t)dim_n);
387 
388  SGVector<float64_t>::add(components[1]->get_mean().vector, alpha1, components[1]->get_mean().vector, alpha2,
389  components[2]->get_mean().vector, dim_n);
390 
391  for (int32_t i=0; i<dim_n; i++)
392  {
393  components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
394  components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
395  }
396 
397  coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2];
398  coefficients.vector[2]=coefficients.vector[0]*0.5;
399  coefficients.vector[0]=coefficients.vector[0]*0.5;
400 
401  if (components[0]->get_cov_type()==FULL)
402  {
403  SGMatrix<float64_t> c1=components[1]->get_cov();
404  SGMatrix<float64_t> c2=components[2]->get_cov();
405  SGVector<float64_t>::add(c1.matrix, alpha1, c1.matrix, alpha2, c2.matrix, dim_n*dim_n);
406 
407  components[1]->set_d(SGVector<float64_t>(SGMatrix<float64_t>::compute_eigenvectors(c1.matrix, dim_n, dim_n), dim_n));
408  components[1]->set_u(c1);
409 
410  float64_t new_d=0;
411  for (int32_t i=0; i<dim_n; i++)
412  {
413  new_d+=CMath::log(components[0]->get_d().vector[i]);
414  for (int32_t j=0; j<dim_n; j++)
415  {
416  if (i==j)
417  {
418  components[0]->get_u().matrix[i*dim_n+j]=1;
419  components[2]->get_u().matrix[i*dim_n+j]=1;
420  }
421  else
422  {
423  components[0]->get_u().matrix[i*dim_n+j]=0;
424  components[2]->get_u().matrix[i*dim_n+j]=0;
425  }
426  }
427  }
428  new_d=CMath::exp(new_d*(1./dim_n));
429  for (int32_t i=0; i<dim_n; i++)
430  {
431  components[0]->get_d().vector[i]=new_d;
432  components[2]->get_d().vector[i]=new_d;
433  }
434  }
435  else if(components[0]->get_cov_type()==DIAG)
436  {
437  SGVector<float64_t>::add(components[1]->get_d().vector, alpha1, components[1]->get_d().vector,
438  alpha2, components[2]->get_d().vector, dim_n);
439 
440  float64_t new_d=0;
441  for (int32_t i=0; i<dim_n; i++)
442  {
443  new_d+=CMath::log(components[0]->get_d().vector[i]);
444  }
445  new_d=CMath::exp(new_d*(1./dim_n));
446  for (int32_t i=0; i<dim_n; i++)
447  {
448  components[0]->get_d().vector[i]=new_d;
449  components[2]->get_d().vector[i]=new_d;
450  }
451  }
452  else if(components[0]->get_cov_type()==SPHERICAL)
453  {
454  components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+
455  alpha2*components[2]->get_d().vector[0];
456 
457  components[2]->get_d().vector[0]=components[0]->get_d().vector[0];
458  }
459 
460  CGMM* partial_candidate=new CGMM(components, coefficients);
461  partial_candidate->train(features);
462 
463  float64_t log_likelihood_prev=0;
464  float64_t log_likelihood_cur=0;
465  int32_t iter=0;
466  SGMatrix<float64_t> alpha(num_vectors, 3);
467  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*3);
468  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
469  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
470 
471  while (iter<max_em_iter)
472  {
473  log_likelihood_prev=log_likelihood_cur;
474  log_likelihood_cur=0;
475 
476  for (int32_t i=0; i<num_vectors; i++)
477  {
478  logPx[i]=0;
480  for (int32_t j=0; j<3; j++)
481  {
482  logPxy[i*3+j]=components[j]->compute_log_PDF(v)+CMath::log(coefficients[j]);
483  logPx[i]+=CMath::exp(logPxy[i*3+j]);
484  }
485 
486  logPx[i]=CMath::log(logPx[i]+init_logPx_fix[i]);
487  log_likelihood_cur+=logPx[i];
488 
489  for (int32_t j=0; j<3; j++)
490  {
491  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
492  alpha.matrix[i*3+j]=CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]);
493  }
494  }
495 
496  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
497  break;
498 
499  partial_candidate->max_likelihood(alpha, min_cov);
500  partial_candidate->get_coef().vector[0]=partial_candidate->get_coef().vector[0]*coef_sum;
501  partial_candidate->get_coef().vector[1]=partial_candidate->get_coef().vector[1]*coef_sum;
502  partial_candidate->get_coef().vector[2]=partial_candidate->get_coef().vector[2]*coef_sum;
503  iter++;
504  }
505 
506  m_coefficients.vector[comp1]=coefficients.vector[0];
507  m_coefficients.vector[comp2]=coefficients.vector[1];
508  m_coefficients.vector[comp3]=coefficients.vector[2];
509 
510  delete partial_candidate;
511  SG_FREE(logPxy);
512  SG_FREE(logPx);
513  SG_FREE(init_logPxy);
514  SG_FREE(init_logPx);
515  SG_FREE(init_logPx_fix);
516  SG_FREE(post_add);
517 }
518 
520 {
521  CDotFeatures* dotdata=(CDotFeatures *) features;
522  int32_t num_dim=dotdata->get_dim_feature_space();
523 
524  float64_t alpha_sum;
525  float64_t alpha_sum_sum=0;
526  float64_t* mean_sum;
527  float64_t* cov_sum=NULL;
528 
529  for (int32_t i=0; i<alpha.num_cols; i++)
530  {
531  alpha_sum=0;
532  mean_sum=SG_MALLOC(float64_t, num_dim);
533  memset(mean_sum, 0, num_dim*sizeof(float64_t));
534 
535  for (int32_t j=0; j<alpha.num_rows; j++)
536  {
537  alpha_sum+=alpha.matrix[j*alpha.num_cols+i];
539  SGVector<float64_t>::add(mean_sum, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, mean_sum, v.vlen);
540  }
541 
542  for (int32_t j=0; j<num_dim; j++)
543  mean_sum[j]/=alpha_sum;
544 
545  m_components[i]->set_mean(SGVector<float64_t>(mean_sum, num_dim));
546 
547  ECovType cov_type=m_components[i]->get_cov_type();
548 
549  if (cov_type==FULL)
550  {
551  cov_sum=SG_MALLOC(float64_t, num_dim*num_dim);
552  memset(cov_sum, 0, num_dim*num_dim*sizeof(float64_t));
553  }
554  else if(cov_type==DIAG)
555  {
556  cov_sum=SG_MALLOC(float64_t, num_dim);
557  memset(cov_sum, 0, num_dim*sizeof(float64_t));
558  }
559  else if(cov_type==SPHERICAL)
560  {
561  cov_sum=SG_MALLOC(float64_t, 1);
562  cov_sum[0]=0;
563  }
564 
565  for (int32_t j=0; j<alpha.num_rows; j++)
566  {
568  SGVector<float64_t>::add(v.vector, 1, v.vector, -1, mean_sum, v.vlen);
569 
570  switch (cov_type)
571  {
572  case FULL:
573  cblas_dger(CblasRowMajor, num_dim, num_dim, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, v.vector,
574  1, (double*) cov_sum, num_dim);
575 
576  break;
577  case DIAG:
578  for (int32_t k=0; k<num_dim; k++)
579  cov_sum[k]+=v.vector[k]*v.vector[k]*alpha.matrix[j*alpha.num_cols+i];
580 
581  break;
582  case SPHERICAL:
583  float64_t temp=0;
584 
585  for (int32_t k=0; k<num_dim; k++)
586  temp+=v.vector[k]*v.vector[k];
587 
588  cov_sum[0]+=temp*alpha.matrix[j*alpha.num_cols+i];
589  break;
590  }
591  }
592 
593  switch (cov_type)
594  {
595  case FULL:
596  for (int32_t j=0; j<num_dim*num_dim; j++)
597  cov_sum[j]/=alpha_sum;
598 
599  float64_t* d0;
600  d0=SGMatrix<float64_t>::compute_eigenvectors(cov_sum, num_dim, num_dim);
601  for (int32_t j=0; j<num_dim; j++)
602  d0[j]=CMath::max(min_cov, d0[j]);
603 
604  m_components[i]->set_d(SGVector<float64_t>(d0, num_dim));
605  m_components[i]->set_u(SGMatrix<float64_t>(cov_sum, num_dim, num_dim));
606 
607  break;
608  case DIAG:
609  for (int32_t j=0; j<num_dim; j++)
610  {
611  cov_sum[j]/=alpha_sum;
612  cov_sum[j]=CMath::max(min_cov, cov_sum[j]);
613  }
614 
615  m_components[i]->set_d(SGVector<float64_t>(cov_sum, num_dim));
616 
617  break;
618  case SPHERICAL:
619  cov_sum[0]/=alpha_sum*num_dim;
620  cov_sum[0]=CMath::max(min_cov, cov_sum[0]);
621 
622  m_components[i]->set_d(SGVector<float64_t>(cov_sum, 1));
623 
624  break;
625  }
626 
627  m_coefficients.vector[i]=alpha_sum;
628  alpha_sum_sum+=alpha_sum;
629  }
630 
631  for (int32_t i=0; i<alpha.num_cols; i++)
632  m_coefficients.vector[i]/=alpha_sum_sum;
633 }
634 
636 {
637  return 1;
638 }
639 
641 {
642  ASSERT(num_param==1)
643 
644  return CMath::log(m_components.size());
645 }
646 
648 {
649  return m_components.size();
650 }
651 
653 {
654  return m_components[index];
655 }
656 
657 float64_t CGMM::get_log_derivative(int32_t num_param, int32_t num_example)
658 {
660  return 0;
661 }
662 
664 {
666  return 1;
667 }
668 
670 {
671  float64_t result=0;
672 
673  ASSERT(features);
676 
677  for (int32_t i=0; i<int32_t(m_components.size()); i++)
678  {
679  SGVector<float64_t> point= ((CDenseFeatures<float64_t>*) features)->get_feature_vector(num_example);
680  result+=CMath::exp(m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]));
681  }
682 
683  return result;
684 }
685 
687 {
688  ASSERT(num<int32_t(m_components.size()))
689  return m_components[num]->get_mean();
690 }
691 
693 {
694  ASSERT(num<int32_t(m_components.size()))
695  m_components[num]->set_mean(mean);
696 }
697 
699 {
700  ASSERT(num<int32_t(m_components.size()))
701  return m_components[num]->get_cov();
702 }
703 
705 {
706  ASSERT(num<int32_t(m_components.size()))
707  m_components[num]->set_cov(cov);
708 }
709 
711 {
712  return m_coefficients;
713 }
714 
715 void CGMM::set_coef(const SGVector<float64_t> coefficients)
716 {
717  m_coefficients=coefficients;
718 }
719 
720 vector<CGaussian*> CGMM::get_comp()
721 {
722  return m_components;
723 }
724 
725 void CGMM::set_comp(vector<CGaussian*> components)
726 {
727  for (int32_t i=0; i<int32_t(m_components.size()); i++)
728  {
730  }
731 
732  m_components=components;
733 
734  for (int32_t i=0; i<int32_t(m_components.size()); i++)
735  {
736  SG_REF(m_components[i]);
737  }
738 }
739 
740 SGMatrix<float64_t> CGMM::alpha_init(SGMatrix<float64_t> init_means)
741 {
742  CDotFeatures* dotdata=(CDotFeatures *) features;
743  int32_t num_vectors=dotdata->get_num_vectors();
744 
745  SGVector<float64_t> label_num(init_means.num_cols);
746 
747  for (int32_t i=0; i<init_means.num_cols; i++)
748  label_num.vector[i]=i;
749 
750  CKNN* knn=new CKNN(1, new CEuclideanDistance(), new CMulticlassLabels(label_num));
751  knn->train(new CDenseFeatures<float64_t>(init_means));
752  CMulticlassLabels* init_labels=(CMulticlassLabels*) knn->apply(features);
753 
754  SGMatrix<float64_t> alpha(num_vectors, int32_t(m_components.size()));
755  memset(alpha.matrix, 0, num_vectors*m_components.size()*sizeof(float64_t));
756 
757  for (int32_t i=0; i<num_vectors; i++)
758  alpha.matrix[i*m_components.size()+init_labels->get_int_label(i)]=1;
759 
760  SG_UNREF(init_labels);
761 
762  return alpha;
763 }
764 
766 {
767  REQUIRE(m_components.size()>0, "Number of mixture components is %d but "
768  "must be positive\n", m_components.size());
769  float64_t rand_num=CMath::random(float64_t(0), float64_t(1));
770  float64_t cum_sum=0;
771  for (int32_t i=0; i<m_coefficients.vlen; i++)
772  {
773  cum_sum+=m_coefficients.vector[i];
774  if (cum_sum>=rand_num)
775  {
776  SG_DEBUG("Sampling from mixture component %d\n", i);
777  return m_components[i]->sample();
778  }
779  }
780  return m_components[m_coefficients.vlen-1]->sample();
781 }
782 
784 {
785  SGVector<float64_t> answer(m_components.size()+1);
786  answer.vector[m_components.size()]=0;
787 
788  for (int32_t i=0; i<int32_t(m_components.size()); i++)
789  {
790  answer.vector[i]=m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]);
791  answer.vector[m_components.size()]+=CMath::exp(answer.vector[i]);
792  }
793  answer.vector[m_components.size()]=CMath::log(answer.vector[m_components.size()]);
794 
795  return answer;
796 }
797 
798 void CGMM::register_params()
799 {
800  //TODO serialization broken
801  //m_parameters->add((SGVector<CSGObject*>*) &m_components, "m_components", "Mixture components");
802  m_parameters->add(&m_coefficients, "m_coefficients", "Mixture coefficients.");
803 }
804 
805 #endif

SHOGUN Machine Learning Toolbox - Documentation