SHOGUN  6.1.3
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/base/Parameter.h>
15 #include <shogun/clustering/GMM.h>
22 #include <shogun/multiclass/KNN.h>
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 
34 {
36  m_components = vector<CGaussian*>(n);
37 
38  for (int32_t i=0; i<n; i++)
39  {
40  m_components[i]=new CGaussian();
41  SG_REF(m_components[i]);
42  m_components[i]->set_cov_type(cov_type);
43  }
44 
45  register_params();
46 }
47 
48 CGMM::CGMM(vector<CGaussian*> components, SGVector<float64_t> coefficients, bool copy) : CDistribution()
49 {
50  ASSERT(int32_t(components.size())==coefficients.vlen)
51 
52  if (!copy)
53  {
54  m_components=components;
55  m_coefficients=coefficients;
56  for (int32_t i=0; i<int32_t(components.size()); i++)
57  {
58  SG_REF(m_components[i]);
59  }
60  }
61  else
62  {
63  m_coefficients = coefficients;
64  m_components = vector<CGaussian*>(components.size());
65 
66  for (int32_t i=0; i<int32_t(components.size()); i++)
67  {
68  m_components[i]=new CGaussian();
69  SG_REF(m_components[i]);
70  m_components[i]->set_cov_type(components[i]->get_cov_type());
71 
72  SGVector<float64_t> old_mean=components[i]->get_mean();
73  SGVector<float64_t> new_mean = old_mean.clone();
74  m_components[i]->set_mean(new_mean);
75 
76  SGVector<float64_t> old_d=components[i]->get_d();
77  SGVector<float64_t> new_d = old_d.clone();
78  m_components[i]->set_d(new_d);
79 
80  if (components[i]->get_cov_type()==FULL)
81  {
82  SGMatrix<float64_t> old_u=components[i]->get_u();
83  SGMatrix<float64_t> new_u = old_u.clone();
84  m_components[i]->set_u(new_u);
85  }
86 
87  m_coefficients[i]=coefficients[i];
88  }
89  }
90 
91  register_params();
92 }
93 
95 {
96  if (!m_components.empty())
97  cleanup();
98 }
99 
101 {
102  for (int32_t i = 0; i < int32_t(m_components.size()); i++)
104 
105  m_components = vector<CGaussian*>();
107 }
108 
110 {
111  ASSERT(m_components.size() != 0)
112 
114  if (data)
115  {
116  if (!data->has_property(FP_DOT))
117  SG_ERROR("Specified features are not of type CDotFeatures\n")
118  set_features(data);
119  }
120 
121  return true;
122 }
123 
124 float64_t CGMM::train_em(float64_t min_cov, int32_t max_iter, float64_t min_change)
125 {
126  if (!features)
127  SG_ERROR("No features to train on.\n")
128 
129  CDotFeatures* dotdata=(CDotFeatures *) features;
130  int32_t num_vectors=dotdata->get_num_vectors();
131 
132  SGMatrix<float64_t> alpha;
133 
134  /* compute initialization via kmeans if none is present */
135  if (m_components[0]->get_mean().vector==NULL)
136  {
137  CKMeans* init_k_means=new CKMeans(int32_t(m_components.size()), new CEuclideanDistance());
138  init_k_means->train(dotdata);
139  SGMatrix<float64_t> init_means=init_k_means->get_cluster_centers();
140 
141  alpha=alpha_init(init_means);
142 
143  SG_UNREF(init_k_means);
144 
145  max_likelihood(alpha, min_cov);
146  }
147  else
148  alpha=SGMatrix<float64_t>(num_vectors,int32_t(m_components.size()));
149 
150  int32_t iter=0;
151  float64_t log_likelihood_prev=0;
152  float64_t log_likelihood_cur=0;
153  SGVector<float64_t> logPxy(num_vectors * m_components.size());
154  SGVector<float64_t> logPx(num_vectors);
155  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
156 
157  while (iter<max_iter)
158  {
159  log_likelihood_prev=log_likelihood_cur;
160  log_likelihood_cur=0;
161 
162  for (int32_t i=0; i<num_vectors; i++)
163  {
164  logPx[i]=0;
166  for (int32_t j=0; j<int32_t(m_components.size()); j++)
167  {
168  logPxy[index_t(i * m_components.size() + j)] =
169  m_components[j]->compute_log_PDF(v) +
171  logPx[i] +=
172  CMath::exp(logPxy[index_t(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(
182  logPxy[index_t(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  return log_likelihood_cur;
195 }
196 
197 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)
198 {
199  if (!features)
200  SG_ERROR("No features to train on.\n")
201 
202  if (m_components.size()<3)
203  SG_ERROR("Can't run SMEM with less than 3 component mixture model.\n")
204 
205  CDotFeatures* dotdata=(CDotFeatures *) features;
206  int32_t num_vectors=dotdata->get_num_vectors();
207 
208  float64_t cur_likelihood=train_em(min_cov, max_em_iter, min_change);
209 
210  int32_t iter=0;
211  SGVector<float64_t> logPxy(num_vectors * m_components.size());
212  SGVector<float64_t> logPx(num_vectors);
213  SGVector<float64_t> logPost(num_vectors * m_components.size());
214  SGVector<float64_t> logPostSum(m_components.size());
215  SGVector<float64_t> logPostSum2(m_components.size());
216  SGVector<float64_t> logPostSumSum(
217  m_components.size() * (m_components.size() - 1) / 2);
218  SGVector<float64_t> split_crit(m_components.size());
219  SGVector<float64_t> merge_crit(
220  m_components.size() * (m_components.size() - 1) / 2);
221  SGVector<int32_t> split_ind(m_components.size());
222  SGVector<int32_t> merge_ind(
223  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[index_t(i * m_components.size() + j)] =
237  m_components[j]->compute_log_PDF(v) +
239  logPx[i] +=
240  CMath::exp(logPxy[index_t(i * m_components.size() + j)]);
241  }
242 
243  logPx[i]=CMath::log(logPx[i]);
244 
245  for (int32_t j=0; j<int32_t(m_components.size()); j++)
246  {
247  logPost[index_t(i * m_components.size() + j)] =
248  logPxy[index_t(i * m_components.size() + j)] - logPx[i];
249  logPostSum[j] +=
250  CMath::exp(logPost[index_t(i * m_components.size() + j)]);
251  logPostSum2[j] += CMath::exp(
252  2 * logPost[index_t(i * m_components.size() + j)]);
253  }
254 
255  int32_t counter=0;
256  for (int32_t j=0; j<int32_t(m_components.size()); j++)
257  {
258  for (int32_t k=j+1; k<int32_t(m_components.size()); k++)
259  {
260  logPostSumSum[counter] += CMath::exp(
261  logPost[index_t(i * m_components.size() + j)] +
262  logPost[index_t(i * m_components.size() + k)]);
263  counter++;
264  }
265  }
266  }
267 
268  int32_t counter=0;
269  for (int32_t i=0; i<int32_t(m_components.size()); i++)
270  {
271  logPostSum[i]=CMath::log(logPostSum[i]);
272  split_crit[i]=0;
273  split_ind[i]=i;
274  for (int32_t j=0; j<num_vectors; j++)
275  {
276  split_crit[i] +=
277  (logPost[index_t(j * m_components.size() + i)] -
278  logPostSum[i] -
279  logPxy[index_t(j * m_components.size() + i)] +
281  (CMath::exp(logPost[index_t(j * m_components.size() + i)]) /
282  CMath::exp(logPostSum[i]));
283  }
284  for (int32_t j=i+1; j<int32_t(m_components.size()); j++)
285  {
286  merge_crit[counter]=CMath::log(logPostSumSum[counter])-(0.5*CMath::log(logPostSum2[i]))-(0.5*CMath::log(logPostSum2[j]));
287  merge_ind[counter]=i*m_components.size()+j;
288  counter++;
289  }
290  }
292  split_crit.vector, split_ind.vector, int32_t(m_components.size()));
294  merge_crit.vector, merge_ind.vector,
295  int32_t(m_components.size() * (m_components.size() - 1) / 2));
296 
297  bool better_found=false;
298  int32_t candidates_checked=0;
299  for (int32_t i=0; i<int32_t(m_components.size()); i++)
300  {
301  for (int32_t j=0; j<int32_t(m_components.size()*(m_components.size()-1)/2); j++)
302  {
303  if (merge_ind[j]/int32_t(m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%m_components.size()) != split_ind[i])
304  {
305  candidates_checked++;
306  CGMM* candidate=new CGMM(m_components, m_coefficients, true);
307  candidate->train(features);
308  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);
309  float64_t cand_likelihood=candidate->train_em(min_cov, max_em_iter, min_change);
310 
311  if (cand_likelihood>cur_likelihood)
312  {
313  cur_likelihood=cand_likelihood;
314  set_comp(candidate->get_comp());
315  set_coef(candidate->get_coef());
316 
317  for (int32_t k=0; k<int32_t(candidate->get_comp().size()); k++)
318  {
319  SG_UNREF(candidate->get_comp()[k]);
320  }
321 
322  better_found=true;
323  break;
324  }
325  else
326  delete candidate;
327 
328  if (candidates_checked>=max_cand)
329  break;
330  }
331  }
332  if (better_found || candidates_checked>=max_cand)
333  break;
334  }
335  if (!better_found)
336  break;
337  iter++;
338  }
339 
340  return cur_likelihood;
341 }
342 
343 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)
344 {
345  CDotFeatures* dotdata=(CDotFeatures *) features;
346  int32_t num_vectors=dotdata->get_num_vectors();
347 
348  SGVector<float64_t> init_logPxy(num_vectors * m_components.size());
349  SGVector<float64_t> init_logPx(num_vectors);
350  SGVector<float64_t> init_logPx_fix(num_vectors);
351  SGVector<float64_t> post_add(num_vectors);
352 
353  for (int32_t i=0; i<num_vectors; i++)
354  {
355  init_logPx[i]=0;
356  init_logPx_fix[i]=0;
357 
359  for (int32_t j=0; j<int32_t(m_components.size()); j++)
360  {
361  init_logPxy[index_t(i * m_components.size() + j)] =
362  m_components[j]->compute_log_PDF(v) +
364  init_logPx[i] +=
365  CMath::exp(init_logPxy[index_t(i * m_components.size() + j)]);
366  if (j!=comp1 && j!=comp2 && j!=comp3)
367  {
368  init_logPx_fix[i] += CMath::exp(
369  init_logPxy[index_t(i * m_components.size() + j)]);
370  }
371  }
372 
373  init_logPx[i]=CMath::log(init_logPx[i]);
374  post_add[i] = CMath::log(
375  CMath::exp(
376  init_logPxy[index_t(i * m_components.size() + comp1)] -
377  init_logPx[i]) +
378  CMath::exp(
379  init_logPxy[index_t(i * m_components.size() + comp2)] -
380  init_logPx[i]) +
381  CMath::exp(
382  init_logPxy[index_t(i * m_components.size() + comp3)] -
383  init_logPx[i]));
384  }
385 
386  vector<CGaussian*> components(3);
387  SGVector<float64_t> coefficients(3);
388  components[0]=m_components[comp1];
389  components[1]=m_components[comp2];
390  components[2]=m_components[comp3];
391  coefficients.vector[0]=m_coefficients.vector[comp1];
392  coefficients.vector[1]=m_coefficients.vector[comp2];
393  coefficients.vector[2]=m_coefficients.vector[comp3];
394  float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2];
395 
396  int32_t dim_n=components[0]->get_mean().vlen;
397 
398  float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]);
399  float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]);
400 
401  float64_t noise_mag=SGVector<float64_t>::twonorm(components[0]->get_mean().vector, dim_n)*0.1/
402  CMath::sqrt((float64_t)dim_n);
403 
404  auto temp_mean = components[2]->get_mean();
405  auto temp_mean_result = components[1]->get_mean();
406  linalg::add(temp_mean_result, temp_mean, temp_mean_result, alpha1, alpha2);
407  components[1]->set_mean(temp_mean_result);
408 
409  for (int32_t i=0; i<dim_n; i++)
410  {
411  components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
412  components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag;
413  }
414 
415  coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2];
416  coefficients.vector[2]=coefficients.vector[0]*0.5;
417  coefficients.vector[0]=coefficients.vector[0]*0.5;
418 
419  if (components[0]->get_cov_type()==FULL)
420  {
421  SGMatrix<float64_t> c1=components[1]->get_cov();
422  SGMatrix<float64_t> c2=components[2]->get_cov();
423  linalg::add(c1, c2, c1, alpha1, alpha2);
424 
425  components[1]->set_d(SGVector<float64_t>(SGMatrix<float64_t>::compute_eigenvectors(c1.matrix, dim_n, dim_n), dim_n));
426  components[1]->set_u(c1);
427 
428  float64_t new_d=0;
429  for (int32_t i=0; i<dim_n; i++)
430  {
431  new_d+=CMath::log(components[0]->get_d().vector[i]);
432  for (int32_t j=0; j<dim_n; j++)
433  {
434  if (i==j)
435  {
436  components[0]->get_u().matrix[i*dim_n+j]=1;
437  components[2]->get_u().matrix[i*dim_n+j]=1;
438  }
439  else
440  {
441  components[0]->get_u().matrix[i*dim_n+j]=0;
442  components[2]->get_u().matrix[i*dim_n+j]=0;
443  }
444  }
445  }
446  new_d=CMath::exp(new_d*(1./dim_n));
447  for (int32_t i=0; i<dim_n; i++)
448  {
449  components[0]->get_d().vector[i]=new_d;
450  components[2]->get_d().vector[i]=new_d;
451  }
452  }
453  else if(components[0]->get_cov_type()==DIAG)
454  {
455  auto result_d = components[1]->get_d();
456  auto temp_d = components[2]->get_d();
457  linalg::add(result_d, temp_d, result_d, alpha1, alpha2);
458  components[1]->set_d(result_d);
459 
460  float64_t new_d=0;
461  for (int32_t i=0; i<dim_n; i++)
462  {
463  new_d+=CMath::log(components[0]->get_d().vector[i]);
464  }
465  new_d=CMath::exp(new_d*(1./dim_n));
466  for (int32_t i=0; i<dim_n; i++)
467  {
468  components[0]->get_d().vector[i]=new_d;
469  components[2]->get_d().vector[i]=new_d;
470  }
471  }
472  else if(components[0]->get_cov_type()==SPHERICAL)
473  {
474  components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+
475  alpha2*components[2]->get_d().vector[0];
476 
477  components[2]->get_d().vector[0]=components[0]->get_d().vector[0];
478  }
479 
480  CGMM* partial_candidate=new CGMM(components, coefficients);
481  partial_candidate->train(features);
482 
483  float64_t log_likelihood_prev=0;
484  float64_t log_likelihood_cur=0;
485  int32_t iter=0;
486  SGMatrix<float64_t> alpha(num_vectors, 3);
487  SGVector<float64_t> logPxy(num_vectors * 3);
488  SGVector<float64_t> logPx(num_vectors);
489  //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen);
490 
491  while (iter<max_em_iter)
492  {
493  log_likelihood_prev=log_likelihood_cur;
494  log_likelihood_cur=0;
495 
496  for (int32_t i=0; i<num_vectors; i++)
497  {
498  logPx[i]=0;
500  for (int32_t j=0; j<3; j++)
501  {
502  logPxy[i*3+j]=components[j]->compute_log_PDF(v)+CMath::log(coefficients[j]);
503  logPx[i]+=CMath::exp(logPxy[i*3+j]);
504  }
505 
506  logPx[i]=CMath::log(logPx[i]+init_logPx_fix[i]);
507  log_likelihood_cur+=logPx[i];
508 
509  for (int32_t j=0; j<3; j++)
510  {
511  //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i];
512  alpha.matrix[i*3+j]=CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]);
513  }
514  }
515 
516  if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change)
517  break;
518 
519  partial_candidate->max_likelihood(alpha, min_cov);
520  partial_candidate->get_coef().vector[0]=partial_candidate->get_coef().vector[0]*coef_sum;
521  partial_candidate->get_coef().vector[1]=partial_candidate->get_coef().vector[1]*coef_sum;
522  partial_candidate->get_coef().vector[2]=partial_candidate->get_coef().vector[2]*coef_sum;
523  iter++;
524  }
525 
526  m_coefficients.vector[comp1]=coefficients.vector[0];
527  m_coefficients.vector[comp2]=coefficients.vector[1];
528  m_coefficients.vector[comp3]=coefficients.vector[2];
529 
530  delete partial_candidate;
531 }
532 
534 {
535  CDotFeatures* dotdata=(CDotFeatures *) features;
536  int32_t num_dim=dotdata->get_dim_feature_space();
537 
538  float64_t alpha_sum;
539  float64_t alpha_sum_sum=0;
540 
541  for (int32_t i=0; i<alpha.num_cols; i++)
542  {
543  alpha_sum=0;
544  SGVector<float64_t> mean_sum(num_dim);
545  linalg::zero(mean_sum);
546 
547  for (int32_t j=0; j<alpha.num_rows; j++)
548  {
549  alpha_sum+=alpha.matrix[j*alpha.num_cols+i];
551  linalg::add(
552  v, mean_sum, mean_sum, alpha.matrix[j * alpha.num_cols + i],
553  1.0);
554  }
555 
556  linalg::scale(mean_sum, mean_sum, 1.0 / alpha_sum);
557 
558  m_components[i]->set_mean(mean_sum);
559 
560  SGMatrix<float64_t> cov_sum;
561 
562  ECovType cov_type = m_components[i]->get_cov_type();
563  if (cov_type==FULL)
564  {
565  cov_sum = SGMatrix<float64_t>(num_dim, num_dim);
566  linalg::zero(cov_sum);
567  }
568  else if(cov_type==DIAG)
569  {
570  cov_sum = SGMatrix<float64_t>(1, num_dim);
571  linalg::zero(cov_sum);
572  }
573  else if(cov_type==SPHERICAL)
574  {
575  cov_sum = SGMatrix<float64_t>(1, 1);
576  linalg::zero(cov_sum);
577  }
578 
579  for (int32_t j=0; j<alpha.num_rows; j++)
580  {
582 
583  linalg::add(v, mean_sum, v, 1.0, -1.0);
584  switch (cov_type)
585  {
586  case FULL:
587  cblas_dger(
588  CblasRowMajor, num_dim, num_dim,
589  alpha.matrix[j * alpha.num_cols + i], v.vector, 1,
590  v.vector, 1, (double*)cov_sum.matrix, num_dim);
591 
592  break;
593  case DIAG:
594  {
595  auto temp_matrix = SGMatrix<float64_t>(v.vector, 1, v.vlen);
596  auto temp_result = linalg::matrix_prod(
597  temp_matrix, temp_matrix, true, false);
598  cov_sum = temp_result.get_diagonal_vector().clone();
600  cov_sum, cov_sum, alpha.matrix[j * alpha.num_cols + i]);
601  }
602 
603  break;
604  case SPHERICAL:
605  float64_t temp = 0;
606 
607  temp = linalg::dot(v, v);
608 
609  cov_sum(0, 0) +=
610  temp * alpha.matrix[j * alpha.num_cols + i];
611  break;
612  }
613  }
614 
615  switch (cov_type)
616  {
617  case FULL:
618  {
619  linalg::scale(cov_sum, cov_sum, 1.0 / alpha_sum);
620 
623  for (int32_t j = 0; j < num_dim; j++)
624  d0[j] = CMath::max(min_cov, d0[j]);
625 
626  m_components[i]->set_d(d0);
627  m_components[i]->set_u(cov_sum);
628 
629  break;
630  }
631  case DIAG:
632  for (int32_t j = 0; j < num_dim; j++)
633  {
634  cov_sum(0, j) /= alpha_sum;
635  cov_sum(0, j) = CMath::max(min_cov, cov_sum(0, j));
636  }
637 
638  m_components[i]->set_d(cov_sum.get_row_vector(0));
639 
640  break;
641  case SPHERICAL:
642  cov_sum[0] /= alpha_sum * num_dim;
643  cov_sum[0] = CMath::max(min_cov, cov_sum[0]);
644 
645  m_components[i]->set_d(cov_sum.get_row_vector(0));
646 
647  break;
648  }
649 
650  m_coefficients.vector[i]=alpha_sum;
651  alpha_sum_sum+=alpha_sum;
652  }
653 
654  linalg::scale(m_coefficients, m_coefficients, 1.0 / alpha_sum_sum);
655 }
656 
658 {
659  return 1;
660 }
661 
663 {
664  ASSERT(num_param==1)
665 
666  return CMath::log(m_components.size());
667 }
668 
670 {
671  return m_components.size();
672 }
673 
675 {
676  return m_components[index];
677 }
678 
679 float64_t CGMM::get_log_derivative(int32_t num_param, int32_t num_example)
680 {
682  return 0;
683 }
684 
686 {
688  return 1;
689 }
690 
692 {
693  float64_t result=0;
694 
695  ASSERT(features);
698 
699  for (int32_t i=0; i<int32_t(m_components.size()); i++)
700  {
701  SGVector<float64_t> point= ((CDenseFeatures<float64_t>*) features)->get_feature_vector(num_example);
702  result+=CMath::exp(m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]));
703  }
704 
705  return result;
706 }
707 
709 {
710  ASSERT(num<int32_t(m_components.size()))
711  return m_components[num]->get_mean();
712 }
713 
715 {
716  ASSERT(num<int32_t(m_components.size()))
717  m_components[num]->set_mean(mean);
718 }
719 
721 {
722  ASSERT(num<int32_t(m_components.size()))
723  return m_components[num]->get_cov();
724 }
725 
727 {
728  ASSERT(num<int32_t(m_components.size()))
729  m_components[num]->set_cov(cov);
730 }
731 
733 {
734  return m_coefficients;
735 }
736 
737 void CGMM::set_coef(const SGVector<float64_t> coefficients)
738 {
739  m_coefficients=coefficients;
740 }
741 
742 vector<CGaussian*> CGMM::get_comp()
743 {
744  return m_components;
745 }
746 
747 void CGMM::set_comp(vector<CGaussian*> components)
748 {
749  for (int32_t i=0; i<int32_t(m_components.size()); i++)
750  {
752  }
753 
754  m_components=components;
755 
756  for (int32_t i=0; i<int32_t(m_components.size()); i++)
757  {
758  SG_REF(m_components[i]);
759  }
760 }
761 
762 SGMatrix<float64_t> CGMM::alpha_init(SGMatrix<float64_t> init_means)
763 {
764  CDotFeatures* dotdata=(CDotFeatures *) features;
765  int32_t num_vectors=dotdata->get_num_vectors();
766 
767  SGVector<float64_t> label_num(init_means.num_cols);
768 
769  for (int32_t i=0; i<init_means.num_cols; i++)
770  label_num[i] = i;
771 
772  CKNN* knn=new CKNN(1, new CEuclideanDistance(), new CMulticlassLabels(label_num));
773  knn->train(new CDenseFeatures<float64_t>(init_means));
774  CMulticlassLabels* init_labels=(CMulticlassLabels*) knn->apply(features);
775 
776  SGMatrix<float64_t> alpha(num_vectors, int32_t(m_components.size()));
777  linalg::zero(alpha);
778 
779  for (int32_t i=0; i<num_vectors; i++)
780  alpha[i * m_components.size() + init_labels->get_int_label(i)] = 1;
781 
782  SG_UNREF(init_labels);
783 
784  return alpha;
785 }
786 
788 {
789  REQUIRE(m_components.size()>0, "Number of mixture components is %d but "
790  "must be positive\n", m_components.size());
791  float64_t rand_num = CMath::random(0.0, 1.0);
792  float64_t cum_sum=0;
793  for (int32_t i=0; i<m_coefficients.vlen; i++)
794  {
795  cum_sum+=m_coefficients.vector[i];
796  if (cum_sum>=rand_num)
797  {
798  SG_DEBUG("Sampling from mixture component %d\n", i);
799  return m_components[i]->sample();
800  }
801  }
802  return m_components[m_coefficients.vlen-1]->sample();
803 }
804 
806 {
807  SGVector<float64_t> answer(m_components.size()+1);
808  answer.vector[m_components.size()]=0;
809 
810  for (int32_t i=0; i<int32_t(m_components.size()); i++)
811  {
812  answer.vector[i]=m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]);
813  answer.vector[m_components.size()]+=CMath::exp(answer.vector[i]);
814  }
815  answer.vector[m_components.size()]=CMath::log(answer.vector[m_components.size()]);
816 
817  return answer;
818 }
819 
820 void CGMM::register_params()
821 {
822  //TODO serialization broken
823  //m_parameters->add((SGVector<CSGObject*>*) &m_components, "m_components", "Mixture components");
824  m_parameters->add(&m_coefficients, "m_coefficients", "Mixture coefficients.");
825 }
826 
827 #endif
virtual float64_t get_likelihood_example(int32_t num_example)
Definition: GMM.cpp:691
SGVector< T > get_row_vector(index_t row) const
Definition: SGMatrix.cpp:1211
virtual int32_t get_num_model_parameters()
Definition: GMM.cpp:657
SGVector< float64_t > m_coefficients
Definition: GMM.h:249
float64_t train_smem(int32_t max_iter=100, int32_t max_cand=5, float64_t min_cov=1e-9, int32_t max_em_iter=1000, float64_t min_change=1e-9)
Definition: GMM.cpp:197
static T max(T a, T b)
Definition: Math.h:149
std::vector< CGaussian * > m_components
Definition: GMM.h:247
int32_t index_t
Definition: common.h:72
virtual void set_features(CFeatures *f)
Definition: Distribution.h:160
virtual float64_t get_log_derivative(int32_t num_param, int32_t num_example)
Definition: GMM.cpp:679
static T sqrt(T x)
Definition: Math.h:428
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
void max_likelihood(SGMatrix< float64_t > alpha, float64_t min_cov)
Definition: GMM.cpp:533
Gaussian distribution interface.
Definition: Gaussian.h:47
virtual int32_t get_num_vectors() const =0
virtual void set_nth_mean(SGVector< float64_t > mean, int32_t num)
Definition: GMM.cpp:714
Definition: basetag.h:132
void add(SGVector< T > &a, SGVector< T > &b, SGVector< T > &result, T alpha=1, T beta=1)
static float64_t randn_double()
Definition: Math.h:924
#define SG_ERROR(...)
Definition: SGIO.h:128
#define REQUIRE(x,...)
Definition: SGIO.h:181
#define SG_NOTIMPLEMENTED
Definition: SGIO.h:138
T dot(const SGVector< T > &a, const SGVector< T > &b)
Parameter * m_parameters
Definition: SGObject.h:609
virtual SGVector< float64_t > get_coef()
Definition: GMM.cpp:732
Base class Distribution from which all methods implementing a distribution are derived.
Definition: Distribution.h:44
std::enable_if<!std::is_same< T, complex128_t >::value, float64_t >::type mean(const Container< T > &a)
Features that support dot products among other operations.
Definition: DotFeatures.h:44
#define SG_REF(x)
Definition: SGObject.h:52
full covariance
Definition: Gaussian.h:33
static uint64_t random()
Definition: Math.h:811
virtual int32_t get_dim_feature_space() const =0
spherical covariance
Definition: Gaussian.h:37
Multiclass Labels for multi-class classification.
SGMatrix< T > clone() const
Definition: SGMatrix.cpp:330
void add(bool *param, const char *name, const char *description="")
Definition: Parameter.cpp:38
virtual ~CGMM()
Definition: GMM.cpp:94
#define ASSERT(x)
Definition: SGIO.h:176
SGVector< float64_t > sample()
Definition: GMM.cpp:787
static SGVector< float64_t > compute_eigenvectors(SGMatrix< float64_t > matrix)
Definition: SGMatrix.cpp:1006
KMeans clustering, partitions the data into k (a-priori specified) clusters.
Definition: KMeans.h:45
void zero(Container< T > &a)
float64_t train_em(float64_t min_cov=1e-9, int32_t max_iter=1000, float64_t min_change=1e-9)
Definition: GMM.cpp:124
virtual SGVector< float64_t > get_nth_mean(int32_t num)
Definition: GMM.cpp:708
CDistribution * get_component(index_t index) const
Definition: GMM.cpp:674
int32_t get_int_label(int32_t idx)
virtual std::vector< CGaussian * > get_comp()
Definition: GMM.cpp:742
double float64_t
Definition: common.h:60
ECovType
Definition: Gaussian.h:30
index_t num_rows
Definition: SGMatrix.h:495
virtual EFeatureClass get_feature_class() const =0
void cleanup()
Definition: GMM.cpp:100
Class KNN, an implementation of the standard k-nearest neigbor classifier.
Definition: KNN.h:74
SGVector< T > clone() const
Definition: SGVector.cpp:262
index_t num_cols
Definition: SGMatrix.h:497
diagonal covariance
Definition: Gaussian.h:35
void matrix_prod(SGMatrix< T > &A, SGVector< T > &b, SGVector< T > &result, bool transpose=false)
virtual void set_coef(const SGVector< float64_t > coefficients)
Definition: GMM.cpp:737
virtual void set_nth_cov(SGMatrix< float64_t > cov, int32_t num)
Definition: GMM.cpp:726
SGVector< float64_t > cluster(SGVector< float64_t > point)
Definition: GMM.cpp:805
#define SG_UNREF(x)
Definition: SGObject.h:53
#define SG_DEBUG(...)
Definition: SGIO.h:106
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
virtual float64_t get_log_likelihood_example(int32_t num_example)
Definition: GMM.cpp:685
virtual float64_t get_log_model_parameter(int32_t num_param)
Definition: GMM.cpp:662
The class Features is the base class of all feature objects.
Definition: Features.h:69
virtual SGMatrix< float64_t > get_nth_cov(int32_t num)
Definition: GMM.cpp:720
static float64_t exp(float64_t x)
Definition: Math.h:551
virtual bool train(CFeatures *data=NULL)
Definition: Machine.cpp:43
static float64_t log(float64_t v)
Definition: Math.h:714
SGMatrix< float64_t > get_cluster_centers()
Definition: KMeansBase.cpp:227
SGVector< float64_t > get_computed_dot_feature_vector(int32_t num)
SGVector< T > get_diagonal_vector() const
Definition: SGMatrix.cpp:1223
virtual bool train(CFeatures *data=NULL)
Definition: GMM.cpp:109
bool has_property(EFeatureProperty p) const
Definition: Features.cpp:295
virtual void set_comp(std::vector< CGaussian * > components)
Definition: GMM.cpp:747
Gaussian Mixture Model interface.
Definition: GMM.h:38
class EuclideanDistance
static void qsort_backward_index(T1 *output, T2 *index, int32_t size)
Definition: Math.h:2066
virtual EFeatureType get_feature_type() const =0
index_t get_num_components() const
Definition: GMM.cpp:669
index_t vlen
Definition: SGVector.h:571
static T twonorm(const T *x, int32_t len)
|| x ||_2
virtual CLabels * apply(CFeatures *data=NULL)
Definition: Machine.cpp:159

SHOGUN Machine Learning Toolbox - Documentation