SHOGUN  v3.0.0
 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  if (m_components[0]->get_mean().vector==NULL)
139  {
140  CKMeans* init_k_means=new CKMeans(int32_t(m_components.size()), new CEuclideanDistance());
141  init_k_means->train(dotdata);
142  SGMatrix<float64_t> init_means=init_k_means->get_cluster_centers();
143 
144  alpha=alpha_init(init_means);
145 
146  SG_UNREF(init_k_means);
147 
148  max_likelihood(alpha, min_cov);
149  }
150  else
151  alpha=SGMatrix<float64_t>(num_vectors,int32_t(m_components.size()));
152 
153  int32_t iter=0;
154  float64_t log_likelihood_prev=0;
155  float64_t log_likelihood_cur=0;
156  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size());
157  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
158  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
159 
160  while (iter<max_iter)
161  {
162  log_likelihood_prev=log_likelihood_cur;
163  log_likelihood_cur=0;
164 
165  for (int32_t i=0; i<num_vectors; i++)
166  {
167  logPx[i]=0;
169  for (int32_t j=0; j<int32_t(m_components.size()); j++)
170  {
171  logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]);
172  logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]);
173  }
174 
175  logPx[i]=CMath::log(logPx[i]);
176  log_likelihood_cur+=logPx[i];
177 
178  for (int32_t j=0; j<int32_t(m_components.size()); j++)
179  {
180  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
181  alpha.matrix[i*m_components.size()+j]=CMath::exp(logPxy[i*m_components.size()+j]-logPx[i]);
182  }
183  }
184 
185  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
186  break;
187 
188  max_likelihood(alpha, min_cov);
189 
190  iter++;
191  }
192 
193  SG_FREE(logPxy);
194  SG_FREE(logPx);
195 
196  return log_likelihood_cur;
197 }
198 
199 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)
200 {
201  if (!features)
202  SG_ERROR("No features to train on.\n")
203 
204  if (m_components.size()<3)
205  SG_ERROR("Can't run SMEM with less than 3 component mixture model.\n")
206 
207  CDotFeatures* dotdata=(CDotFeatures *) features;
208  int32_t num_vectors=dotdata->get_num_vectors();
209 
210  float64_t cur_likelihood=train_em(min_cov, max_em_iter, min_change);
211 
212  int32_t iter=0;
213  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size());
214  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
215  float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.size());
216  float64_t* logPostSum=SG_MALLOC(float64_t, m_components.size());
217  float64_t* logPostSum2=SG_MALLOC(float64_t, m_components.size());
218  float64_t* logPostSumSum=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2);
219  float64_t* split_crit=SG_MALLOC(float64_t, m_components.size());
220  float64_t* merge_crit=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2);
221  int32_t* split_ind=SG_MALLOC(int32_t, m_components.size());
222  int32_t* merge_ind=SG_MALLOC(int32_t, m_components.size()*(m_components.size()-1)/2);
223 
224  while (iter<max_iter)
225  {
226  memset(logPostSum, 0, m_components.size()*sizeof(float64_t));
227  memset(logPostSum2, 0, m_components.size()*sizeof(float64_t));
228  memset(logPostSumSum, 0, (m_components.size()*(m_components.size()-1)/2)*sizeof(float64_t));
229  for (int32_t i=0; i<num_vectors; i++)
230  {
231  logPx[i]=0;
233  for (int32_t j=0; j<int32_t(m_components.size()); j++)
234  {
235  logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]);
236  logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]);
237  }
238 
239  logPx[i]=CMath::log(logPx[i]);
240 
241  for (int32_t j=0; j<int32_t(m_components.size()); j++)
242  {
243  logPost[i*m_components.size()+j]=logPxy[i*m_components.size()+j]-logPx[i];
244  logPostSum[j]+=CMath::exp(logPost[i*m_components.size()+j]);
245  logPostSum2[j]+=CMath::exp(2*logPost[i*m_components.size()+j]);
246  }
247 
248  int32_t counter=0;
249  for (int32_t j=0; j<int32_t(m_components.size()); j++)
250  {
251  for (int32_t k=j+1; k<int32_t(m_components.size()); k++)
252  {
253  logPostSumSum[counter]+=CMath::exp(logPost[i*m_components.size()+j]+logPost[i*m_components.size()+k]);
254  counter++;
255  }
256  }
257  }
258 
259  int32_t counter=0;
260  for (int32_t i=0; i<int32_t(m_components.size()); i++)
261  {
262  logPostSum[i]=CMath::log(logPostSum[i]);
263  split_crit[i]=0;
264  split_ind[i]=i;
265  for (int32_t j=0; j<num_vectors; j++)
266  {
267  split_crit[i]+=(logPost[j*m_components.size()+i]-logPostSum[i]-logPxy[j*m_components.size()+i]+CMath::log(m_coefficients[i]))*
268  (CMath::exp(logPost[j*m_components.size()+i])/CMath::exp(logPostSum[i]));
269  }
270  for (int32_t j=i+1; j<int32_t(m_components.size()); j++)
271  {
272  merge_crit[counter]=CMath::log(logPostSumSum[counter])-(0.5*CMath::log(logPostSum2[i]))-(0.5*CMath::log(logPostSum2[j]));
273  merge_ind[counter]=i*m_components.size()+j;
274  counter++;
275  }
276  }
277  CMath::qsort_backward_index(split_crit, split_ind, int32_t(m_components.size()));
278  CMath::qsort_backward_index(merge_crit, merge_ind, int32_t(m_components.size()*(m_components.size()-1)/2));
279 
280  bool better_found=false;
281  int32_t candidates_checked=0;
282  for (int32_t i=0; i<int32_t(m_components.size()); i++)
283  {
284  for (int32_t j=0; j<int32_t(m_components.size()*(m_components.size()-1)/2); j++)
285  {
286  if (merge_ind[j]/int32_t(m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%m_components.size()) != split_ind[i])
287  {
288  candidates_checked++;
289  CGMM* candidate=new CGMM(m_components, m_coefficients, true);
290  candidate->train(features);
291  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);
292  float64_t cand_likelihood=candidate->train_em(min_cov, max_em_iter, min_change);
293 
294  if (cand_likelihood>cur_likelihood)
295  {
296  cur_likelihood=cand_likelihood;
297  set_comp(candidate->get_comp());
298  set_coef(candidate->get_coef());
299 
300  for (int32_t k=0; k<int32_t(candidate->get_comp().size()); k++)
301  {
302  SG_UNREF(candidate->get_comp()[k]);
303  }
304 
305  better_found=true;
306  break;
307  }
308  else
309  delete candidate;
310 
311  if (candidates_checked>=max_cand)
312  break;
313  }
314  }
315  if (better_found || candidates_checked>=max_cand)
316  break;
317  }
318  if (!better_found)
319  break;
320  iter++;
321  }
322 
323  SG_FREE(logPxy);
324  SG_FREE(logPx);
325  SG_FREE(logPost);
326  SG_FREE(split_crit);
327  SG_FREE(merge_crit);
328  SG_FREE(logPostSum);
329  SG_FREE(logPostSum2);
330  SG_FREE(logPostSumSum);
331  SG_FREE(split_ind);
332  SG_FREE(merge_ind);
333 
334  return cur_likelihood;
335 }
336 
337 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)
338 {
339  CDotFeatures* dotdata=(CDotFeatures *) features;
340  int32_t num_vectors=dotdata->get_num_vectors();
341 
342  float64_t* init_logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size());
343  float64_t* init_logPx=SG_MALLOC(float64_t, num_vectors);
344  float64_t* init_logPx_fix=SG_MALLOC(float64_t, num_vectors);
345  float64_t* post_add=SG_MALLOC(float64_t, num_vectors);
346 
347  for (int32_t i=0; i<num_vectors; i++)
348  {
349  init_logPx[i]=0;
350  init_logPx_fix[i]=0;
351 
353  for (int32_t j=0; j<int32_t(m_components.size()); j++)
354  {
355  init_logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]);
356  init_logPx[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]);
357  if (j!=comp1 && j!=comp2 && j!=comp3)
358  {
359  init_logPx_fix[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]);
360  }
361  }
362 
363  init_logPx[i]=CMath::log(init_logPx[i]);
364  post_add[i]=CMath::log(CMath::exp(init_logPxy[i*m_components.size()+comp1]-init_logPx[i])+
365  CMath::exp(init_logPxy[i*m_components.size()+comp2]-init_logPx[i])+
366  CMath::exp(init_logPxy[i*m_components.size()+comp3]-init_logPx[i]));
367  }
368 
369  vector<CGaussian*> components(3);
370  SGVector<float64_t> coefficients(3);
371  components[0]=m_components[comp1];
372  components[1]=m_components[comp2];
373  components[2]=m_components[comp3];
374  coefficients.vector[0]=m_coefficients.vector[comp1];
375  coefficients.vector[1]=m_coefficients.vector[comp2];
376  coefficients.vector[2]=m_coefficients.vector[comp3];
377  float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2];
378 
379  int32_t dim_n=components[0]->get_mean().vlen;
380 
381  float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]);
382  float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]);
383 
384  float64_t noise_mag=SGVector<float64_t>::twonorm(components[0]->get_mean().vector, dim_n)*0.1/
385  CMath::sqrt((float64_t)dim_n);
386 
387  SGVector<float64_t>::add(components[1]->get_mean().vector, alpha1, components[1]->get_mean().vector, alpha2,
388  components[2]->get_mean().vector, dim_n);
389 
390  for (int32_t i=0; i<dim_n; i++)
391  {
392  components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
393  components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
394  }
395 
396  coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2];
397  coefficients.vector[2]=coefficients.vector[0]*0.5;
398  coefficients.vector[0]=coefficients.vector[0]*0.5;
399 
400  if (components[0]->get_cov_type()==FULL)
401  {
402  SGMatrix<float64_t> c1=components[1]->get_cov();
403  SGMatrix<float64_t> c2=components[2]->get_cov();
404  SGVector<float64_t>::add(c1.matrix, alpha1, c1.matrix, alpha2, c2.matrix, dim_n*dim_n);
405 
406  components[1]->set_d(SGVector<float64_t>(SGMatrix<float64_t>::compute_eigenvectors(c1.matrix, dim_n, dim_n), dim_n));
407  components[1]->set_u(c1);
408 
409  float64_t new_d=0;
410  for (int32_t i=0; i<dim_n; i++)
411  {
412  new_d+=CMath::log(components[0]->get_d().vector[i]);
413  for (int32_t j=0; j<dim_n; j++)
414  {
415  if (i==j)
416  {
417  components[0]->get_u().matrix[i*dim_n+j]=1;
418  components[2]->get_u().matrix[i*dim_n+j]=1;
419  }
420  else
421  {
422  components[0]->get_u().matrix[i*dim_n+j]=0;
423  components[2]->get_u().matrix[i*dim_n+j]=0;
424  }
425  }
426  }
427  new_d=CMath::exp(new_d*(1./dim_n));
428  for (int32_t i=0; i<dim_n; i++)
429  {
430  components[0]->get_d().vector[i]=new_d;
431  components[2]->get_d().vector[i]=new_d;
432  }
433  }
434  else if(components[0]->get_cov_type()==DIAG)
435  {
436  SGVector<float64_t>::add(components[1]->get_d().vector, alpha1, components[1]->get_d().vector,
437  alpha2, components[2]->get_d().vector, dim_n);
438 
439  float64_t new_d=0;
440  for (int32_t i=0; i<dim_n; i++)
441  {
442  new_d+=CMath::log(components[0]->get_d().vector[i]);
443  }
444  new_d=CMath::exp(new_d*(1./dim_n));
445  for (int32_t i=0; i<dim_n; i++)
446  {
447  components[0]->get_d().vector[i]=new_d;
448  components[2]->get_d().vector[i]=new_d;
449  }
450  }
451  else if(components[0]->get_cov_type()==SPHERICAL)
452  {
453  components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+
454  alpha2*components[2]->get_d().vector[0];
455 
456  components[2]->get_d().vector[0]=components[0]->get_d().vector[0];
457  }
458 
459  CGMM* partial_candidate=new CGMM(components, coefficients);
460  partial_candidate->train(features);
461 
462  float64_t log_likelihood_prev=0;
463  float64_t log_likelihood_cur=0;
464  int32_t iter=0;
465  SGMatrix<float64_t> alpha(num_vectors, 3);
466  float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*3);
467  float64_t* logPx=SG_MALLOC(float64_t, num_vectors);
468  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
469 
470  while (iter<max_em_iter)
471  {
472  log_likelihood_prev=log_likelihood_cur;
473  log_likelihood_cur=0;
474 
475  for (int32_t i=0; i<num_vectors; i++)
476  {
477  logPx[i]=0;
479  for (int32_t j=0; j<3; j++)
480  {
481  logPxy[i*3+j]=components[j]->compute_log_PDF(v)+CMath::log(coefficients[j]);
482  logPx[i]+=CMath::exp(logPxy[i*3+j]);
483  }
484 
485  logPx[i]=CMath::log(logPx[i]+init_logPx_fix[i]);
486  log_likelihood_cur+=logPx[i];
487 
488  for (int32_t j=0; j<3; j++)
489  {
490  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
491  alpha.matrix[i*3+j]=CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]);
492  }
493  }
494 
495  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
496  break;
497 
498  partial_candidate->max_likelihood(alpha, min_cov);
499  partial_candidate->get_coef().vector[0]=partial_candidate->get_coef().vector[0]*coef_sum;
500  partial_candidate->get_coef().vector[1]=partial_candidate->get_coef().vector[1]*coef_sum;
501  partial_candidate->get_coef().vector[2]=partial_candidate->get_coef().vector[2]*coef_sum;
502  iter++;
503  }
504 
505  m_coefficients.vector[comp1]=coefficients.vector[0];
506  m_coefficients.vector[comp2]=coefficients.vector[1];
507  m_coefficients.vector[comp3]=coefficients.vector[2];
508 
509  delete partial_candidate;
510  SG_FREE(logPxy);
511  SG_FREE(logPx);
512  SG_FREE(init_logPxy);
513  SG_FREE(init_logPx);
514  SG_FREE(init_logPx_fix);
515  SG_FREE(post_add);
516 }
517 
519 {
520  CDotFeatures* dotdata=(CDotFeatures *) features;
521  int32_t num_dim=dotdata->get_dim_feature_space();
522 
523  float64_t alpha_sum;
524  float64_t alpha_sum_sum=0;
525  float64_t* mean_sum;
526  float64_t* cov_sum=NULL;
527 
528  for (int32_t i=0; i<alpha.num_cols; i++)
529  {
530  alpha_sum=0;
531  mean_sum=SG_MALLOC(float64_t, num_dim);
532  memset(mean_sum, 0, num_dim*sizeof(float64_t));
533 
534  for (int32_t j=0; j<alpha.num_rows; j++)
535  {
536  alpha_sum+=alpha.matrix[j*alpha.num_cols+i];
538  SGVector<float64_t>::add(mean_sum, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, mean_sum, v.vlen);
539  }
540 
541  for (int32_t j=0; j<num_dim; j++)
542  mean_sum[j]/=alpha_sum;
543 
544  m_components[i]->set_mean(SGVector<float64_t>(mean_sum, num_dim));
545 
546  ECovType cov_type=m_components[i]->get_cov_type();
547 
548  if (cov_type==FULL)
549  {
550  cov_sum=SG_MALLOC(float64_t, num_dim*num_dim);
551  memset(cov_sum, 0, num_dim*num_dim*sizeof(float64_t));
552  }
553  else if(cov_type==DIAG)
554  {
555  cov_sum=SG_MALLOC(float64_t, num_dim);
556  memset(cov_sum, 0, num_dim*sizeof(float64_t));
557  }
558  else if(cov_type==SPHERICAL)
559  {
560  cov_sum=SG_MALLOC(float64_t, 1);
561  cov_sum[0]=0;
562  }
563 
564  for (int32_t j=0; j<alpha.num_rows; j++)
565  {
567  SGVector<float64_t>::add(v.vector, 1, v.vector, -1, mean_sum, v.vlen);
568 
569  switch (cov_type)
570  {
571  case FULL:
572  cblas_dger(CblasRowMajor, num_dim, num_dim, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, v.vector,
573  1, (double*) cov_sum, num_dim);
574 
575  break;
576  case DIAG:
577  for (int32_t k=0; k<num_dim; k++)
578  cov_sum[k]+=v.vector[k]*v.vector[k]*alpha.matrix[j*alpha.num_cols+i];
579 
580  break;
581  case SPHERICAL:
582  float64_t temp=0;
583 
584  for (int32_t k=0; k<num_dim; k++)
585  temp+=v.vector[k]*v.vector[k];
586 
587  cov_sum[0]+=temp*alpha.matrix[j*alpha.num_cols+i];
588  break;
589  }
590  }
591 
592  switch (cov_type)
593  {
594  case FULL:
595  for (int32_t j=0; j<num_dim*num_dim; j++)
596  cov_sum[j]/=alpha_sum;
597 
598  float64_t* d0;
599  d0=SGMatrix<float64_t>::compute_eigenvectors(cov_sum, num_dim, num_dim);
600  for (int32_t j=0; j<num_dim; j++)
601  d0[j]=CMath::max(min_cov, d0[j]);
602 
603  m_components[i]->set_d(SGVector<float64_t>(d0, num_dim));
604  m_components[i]->set_u(SGMatrix<float64_t>(cov_sum, num_dim, num_dim));
605 
606  break;
607  case DIAG:
608  for (int32_t j=0; j<num_dim; j++)
609  {
610  cov_sum[j]/=alpha_sum;
611  cov_sum[j]=CMath::max(min_cov, cov_sum[j]);
612  }
613 
614  m_components[i]->set_d(SGVector<float64_t>(cov_sum, num_dim));
615 
616  break;
617  case SPHERICAL:
618  cov_sum[0]/=alpha_sum*num_dim;
619  cov_sum[0]=CMath::max(min_cov, cov_sum[0]);
620 
621  m_components[i]->set_d(SGVector<float64_t>(cov_sum, 1));
622 
623  break;
624  }
625 
626  m_coefficients.vector[i]=alpha_sum;
627  alpha_sum_sum+=alpha_sum;
628  }
629 
630  for (int32_t i=0; i<alpha.num_cols; i++)
631  m_coefficients.vector[i]/=alpha_sum_sum;
632 }
633 
635 {
636  return 1;
637 }
638 
640 {
641  ASSERT(num_param==1)
642 
643  return CMath::log(m_components.size());
644 }
645 
646 float64_t CGMM::get_log_derivative(int32_t num_param, int32_t num_example)
647 {
649  return 0;
650 }
651 
653 {
655  return 1;
656 }
657 
659 {
660  return CMath::exp(get_log_likelihood_example(num_example));
661 }
662 
664 {
665  ASSERT(num<int32_t(m_components.size()))
666  return m_components[num]->get_mean();
667 }
668 
670 {
671  ASSERT(num<int32_t(m_components.size()))
672  m_components[num]->set_mean(mean);
673 }
674 
676 {
677  ASSERT(num<int32_t(m_components.size()))
678  return m_components[num]->get_cov();
679 }
680 
682 {
683  ASSERT(num<int32_t(m_components.size()))
684  m_components[num]->set_cov(cov);
685 }
686 
688 {
689  return m_coefficients;
690 }
691 
692 void CGMM::set_coef(const SGVector<float64_t> coefficients)
693 {
694  m_coefficients=coefficients;
695 }
696 
697 vector<CGaussian*> CGMM::get_comp()
698 {
699  return m_components;
700 }
701 
702 void CGMM::set_comp(vector<CGaussian*> components)
703 {
704  for (int32_t i=0; i<int32_t(m_components.size()); i++)
705  {
707  }
708 
709  m_components=components;
710 
711  for (int32_t i=0; i<int32_t(m_components.size()); i++)
712  {
713  SG_REF(m_components[i]);
714  }
715 }
716 
717 SGMatrix<float64_t> CGMM::alpha_init(SGMatrix<float64_t> init_means)
718 {
719  CDotFeatures* dotdata=(CDotFeatures *) features;
720  int32_t num_vectors=dotdata->get_num_vectors();
721 
722  SGVector<float64_t> label_num(init_means.num_cols);
723 
724  for (int32_t i=0; i<init_means.num_cols; i++)
725  label_num.vector[i]=i;
726 
727  CKNN* knn=new CKNN(1, new CEuclideanDistance(), new CMulticlassLabels(label_num));
728  knn->train(new CDenseFeatures<float64_t>(init_means));
729  CMulticlassLabels* init_labels=(CMulticlassLabels*) knn->apply(features);
730 
731  SGMatrix<float64_t> alpha(num_vectors, int32_t(m_components.size()));
732  memset(alpha.matrix, 0, num_vectors*m_components.size()*sizeof(float64_t));
733 
734  for (int32_t i=0; i<num_vectors; i++)
735  alpha.matrix[i*m_components.size()+init_labels->get_int_label(i)]=1;
736 
737  SG_UNREF(init_labels);
738 
739  return alpha;
740 }
741 
743 {
744  ASSERT(m_components.size())
745  float64_t rand_num=CMath::random(float64_t(0), float64_t(1));
746  float64_t cum_sum=0;
747  for (int32_t i=0; i<m_coefficients.vlen; i++)
748  {
749  cum_sum+=m_coefficients.vector[i];
750  if (cum_sum>=rand_num)
751  return m_components[i]->sample();
752  }
753  return m_components[m_coefficients.vlen-1]->sample();
754 }
755 
757 {
758  SGVector<float64_t> answer(m_components.size()+1);
759  answer.vector[m_components.size()]=0;
760 
761  for (int32_t i=0; i<int32_t(m_components.size()); i++)
762  {
763  answer.vector[i]=m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]);
764  answer.vector[m_components.size()]+=CMath::exp(answer.vector[i]);
765  }
766  answer.vector[m_components.size()]=CMath::log(answer.vector[m_components.size()]);
767 
768  return answer;
769 }
770 
771 void CGMM::register_params()
772 {
773  //TODO serialization broken
774  //m_parameters->add((SGVector<CSGObject*>*) &m_components, "m_components", "Mixture components");
775  m_parameters->add(&m_coefficients, "m_coefficients", "Mixture coefficients.");
776 }
777 
778 #endif

SHOGUN Machine Learning Toolbox - Documentation