SHOGUN  6.1.3
Gaussian.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  * Written (W) 2014 Parijat Mazumdar
9  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
10  */
11 #include <shogun/lib/config.h>
12 
13 #include <shogun/base/Parameter.h>
19 
20 using namespace shogun;
21 using namespace linalg;
22 
23 CGaussian::CGaussian() : CDistribution(), m_constant(0), m_d(), m_u(), m_mean(), m_cov_type(FULL)
24 {
25  register_params();
26 }
27 
30  : CDistribution()
31 {
32  ASSERT(mean.vlen==cov.num_rows)
33  ASSERT(cov.num_rows==cov.num_cols)
36  m_cov_type=cov_type;
37 
38  m_mean=mean;
39 
40  if (cov.num_rows==1)
41  m_cov_type=SPHERICAL;
42 
43  decompose_cov(cov);
44  init();
45  register_params();
46 }
47 
49 {
51  switch (m_cov_type)
52  {
53  case FULL:
54  case DIAG:
55  for (int32_t i=0; i<m_d.vlen; i++)
57  break;
58  case SPHERICAL:
60  break;
61  }
62 }
63 
65 {
66 }
67 
69 {
70  // init features with data if necessary and assure type is correct
71  if (data)
72  {
73  if (!data->has_property(FP_DOT))
74  SG_ERROR("Specified features are not of type CDotFeatures\n")
75  set_features(data);
76  }
77 
78  CDotFeatures* dotdata=(CDotFeatures *) data;
79  m_mean=dotdata->get_mean();
80  SGMatrix<float64_t> cov=dotdata->get_cov();
81  decompose_cov(cov);
82  init();
83  return true;
84 }
85 
87 {
88  switch (m_cov_type)
89  {
90  case FULL:
92  case DIAG:
93  return m_d.vlen+m_mean.vlen;
94  case SPHERICAL:
95  return 1+m_mean.vlen;
96  }
97  return 0;
98 }
99 
101 {
103  return 0;
104 }
105 
106 float64_t CGaussian::get_log_derivative(int32_t num_param, int32_t num_example)
107 {
109  return 0;
110 }
111 
113 {
115  SGVector<float64_t> v=((CDotFeatures *)features)->get_computed_dot_feature_vector(num_example);
116  float64_t answer=compute_log_PDF(v);
117  return answer;
118 }
119 
121 {
122  CDotFeatures* dotdata=dynamic_cast<CDotFeatures *>(features);
123  REQUIRE(
124  dotdata, "dynamic cast from CFeatures to CDotFeatures returned NULL\n");
125  int32_t num_dim=dotdata->get_dim_feature_space();
126 
127  // compute mean
128 
129  float64_t alpha_k_sum=0;
130  SGVector<float64_t> mean(num_dim);
131  mean.fill_vector(mean.vector, mean.vlen, 0);
132  for (int32_t i = 0; i < len; i++)
133  {
134  alpha_k_sum+=alpha_k[i];
136  linalg::add(v, mean, mean, alpha_k[i], 1.0);
137  }
138 
139  linalg::scale(mean, mean, 1.0 / alpha_k_sum);
140 
141  set_mean(mean);
142 
143  // compute covariance matrix
144 
145  SGMatrix<float64_t> cov_sum;
146  ECovType cov_type=get_cov_type();
147  if (cov_type==FULL)
148  {
149  cov_sum = SGMatrix<float64_t>(num_dim, num_dim);
150  cov_sum.zero();
151  }
152  else if(cov_type==DIAG)
153  {
154  cov_sum = SGMatrix<float64_t>(1, num_dim);
155  cov_sum.zero();
156  }
157  else if(cov_type==SPHERICAL)
158  {
159  cov_sum = SGMatrix<float64_t>(1, 1);
160  cov_sum.zero();
161  }
162 
163  for (int32_t j=0; j<len; j++)
164  {
166  linalg::add(v, mean, v, -1.0, 1.0);
167 
168  switch (cov_type)
169  {
170  case FULL:
171 #ifdef HAVE_LAPACK
172  cblas_dger(
173  CblasRowMajor, num_dim, num_dim, alpha_k[j], v.vector, 1,
174  v.vector, 1, (double*)cov_sum.matrix, num_dim);
175 #else
176  linalg::dger<float64_t>(alpha_k[j], v, v, cov_sum);
177 #endif
178  break;
179  case DIAG:
180  for (int32_t k = 0; k < num_dim; k++)
181  cov_sum(1, k) += v.vector[k] * v.vector[k] * alpha_k[j];
182 
183  break;
184  case SPHERICAL:
185  float64_t temp = 0;
186 
187  temp = linalg::dot(v, v);
188 
189  cov_sum(0, 0) += temp * alpha_k[j];
190  break;
191  }
192  }
193 
194  switch (cov_type)
195  {
196  case FULL:
197  {
198  linalg::scale(cov_sum, cov_sum, 1 / alpha_k_sum);
199 
200  SGVector<float64_t> d0(num_dim);
201 #ifdef HAVE_LAPACK
203  cov_sum.matrix, num_dim, num_dim);
204 #else
205  // FIXME use eigenvectors computeation warpper by micmn
206  typename SGMatrix<float64_t>::EigenMatrixXtMap eig = cov_sum;
207  typename SGVector<float64_t>::EigenVectorXtMap eigenvalues_eig = d0;
208 
209  Eigen::EigenSolver<typename SGMatrix<float64_t>::EigenMatrixXt> solver(
210  eig);
211  eigenvalues_eig = solver.eigenvalues().real();
212 #endif
213 
214  set_d(d0);
215  set_u(cov_sum);
216 
217  break;
218  }
219  case DIAG:
220  linalg::scale(cov_sum, cov_sum, 1 / alpha_k_sum);
221 
222  set_d(cov_sum.get_row_vector(0));
223 
224  break;
225 
226  case SPHERICAL:
227  cov_sum[0] /= alpha_k_sum * num_dim;
228 
229  set_d(cov_sum.get_row_vector(0));
230 
231  break;
232  }
233 
234  return alpha_k_sum;
235 }
236 
238 {
240  ASSERT(point.vlen == m_mean.vlen)
241  SGVector<float64_t> difference = point.clone();
242 
243  linalg::add(difference, m_mean, difference, -1.0, 1.0);
244 
245  float64_t answer=m_constant;
246 
247  if (m_cov_type==FULL)
248  {
249  SGVector<float64_t> temp_holder(m_d.vlen);
250  temp_holder.zero();
251 #ifdef HAVE_LAPACK
252  cblas_dgemv(
253  CblasRowMajor, CblasNoTrans, m_d.vlen, m_d.vlen, 1, m_u.matrix,
254  m_d.vlen, difference, 1, 0, temp_holder, 1);
255 #else
256  linalg::dgemv<float64_t>(1, m_u, false, difference, 0, temp_holder);
257 #endif
258 
259  for (int32_t i=0; i<m_d.vlen; i++)
260  answer+=temp_holder[i]*temp_holder[i]/m_d.vector[i];
261  }
262  else if (m_cov_type==DIAG)
263  {
264  for (int32_t i=0; i<m_mean.vlen; i++)
265  answer+=difference[i]*difference[i]/m_d.vector[i];
266  }
267  else
268  {
269  for (int32_t i=0; i<m_mean.vlen; i++)
270  answer += difference[i] * difference[i] / m_d.vector[0];
271  }
272 
273  return -0.5 * answer;
274 }
275 
277 {
278  return m_mean;
279 }
280 
282 {
283  if (mean.vlen==1)
285 
286  m_mean=mean;
287 }
288 
290 {
291  ASSERT(cov.num_rows==cov.num_cols)
292  ASSERT(cov.num_rows==m_mean.vlen)
293  decompose_cov(cov);
294  init();
295 }
296 
298 {
299  m_d = d;
300  init();
301 }
302 
304 {
306  cov.zero();
307 
308  if (m_cov_type==FULL)
309  {
310  if (!m_u.matrix)
311  SG_ERROR("Unitary matrix not set\n")
312 
313  SGMatrix<float64_t> temp_holder(m_mean.vlen, m_mean.vlen);
314  SGMatrix<float64_t> diag_holder(m_mean.vlen, m_mean.vlen);
315  diag_holder.zero();
316  for (int32_t i = 0; i < m_d.vlen; i++)
317  diag_holder(i, i) = m_d.vector[i];
318 #ifdef HAVE_LAPACK
319  cblas_dgemm(
320  CblasRowMajor, CblasTrans, CblasNoTrans, m_d.vlen, m_d.vlen,
321  m_d.vlen, 1, m_u.matrix, m_d.vlen, diag_holder.matrix, m_d.vlen, 0,
322  temp_holder.matrix, m_d.vlen);
323  cblas_dgemm(
324  CblasRowMajor, CblasNoTrans, CblasNoTrans, m_d.vlen, m_d.vlen,
325  m_d.vlen, 1, temp_holder.matrix, m_d.vlen, m_u.matrix, m_d.vlen, 0,
326  cov.matrix, m_d.vlen);
327 #else
328  linalg::dgemm<float64_t>(
329  1, m_u, diag_holder, true, false, 0, temp_holder);
330  linalg::dgemm<float64_t>(1, temp_holder, m_u, false, false, 0, cov);
331 #endif
332  }
333  else if (m_cov_type == DIAG)
334  {
335  for (int32_t i = 0; i < m_d.vlen; i++)
336  cov(i, i) = m_d.vector[i];
337  }
338  else
339  {
340  for (int32_t i = 0; i < m_mean.vlen; i++)
341  cov(i, i) = m_d.vector[0];
342  }
343  return cov;
344 }
345 
346 void CGaussian::register_params()
347 {
348  SG_ADD(&m_u, "m_u", "Unitary matrix.",MS_NOT_AVAILABLE);
349  SG_ADD(&m_d, "m_d", "Diagonal.",MS_NOT_AVAILABLE);
350  SG_ADD(&m_mean, "m_mean", "Mean.",MS_NOT_AVAILABLE);
351  SG_ADD(&m_constant, "m_constant", "Constant part.",MS_NOT_AVAILABLE);
352  SG_ADD((machine_int_t*)&m_cov_type, "m_cov_type", "Covariance type.",MS_NOT_AVAILABLE);
353 }
354 
355 void CGaussian::decompose_cov(SGMatrix<float64_t> cov)
356 {
357  switch (m_cov_type)
358  {
359  case FULL:
360  {
362  m_u = cov.clone();
364 #ifdef HAVE_LAPACK
366  m_u.matrix, cov.num_rows, cov.num_rows);
367 #else
368  // FIXME use eigenvectors computeation warpper by micmn
370  typename SGVector<float64_t>::EigenVectorXtMap eigenvalues_eig = m_d;
371 
372  Eigen::EigenSolver<typename SGMatrix<float64_t>::EigenMatrixXt> solver(
373  eig);
374  eigenvalues_eig = solver.eigenvalues().real();
375 #endif
376  break;
377  }
378  case DIAG:
380  for (int32_t i = 0; i < cov.num_rows; i++)
381  m_d[i] = cov.matrix[i * cov.num_rows + i];
382 
383  break;
384  case SPHERICAL:
386  m_d.vector[0] = cov.matrix[0];
387  break;
388  }
389 }
390 
392 {
393  SG_DEBUG("Entering\n");
395  r_matrix.zero();
396 
397  switch (m_cov_type)
398  {
399  case FULL:
400  case DIAG:
401  for (int32_t i = 0; i < m_mean.vlen; i++)
402  r_matrix(i, i) = CMath::sqrt(m_d.vector[i]);
403 
404  break;
405  case SPHERICAL:
406  for (int32_t i = 0; i < m_mean.vlen; i++)
407  r_matrix(i, i) = CMath::sqrt(m_d.vector[0]);
408 
409  break;
410  }
411 
412  SGVector<float64_t> random_vec(m_mean.vlen);
413 
414  for (int32_t i = 0; i < m_mean.vlen; i++)
415  random_vec.vector[i] = CMath::randn_double();
416 
417  if (m_cov_type == FULL)
418  {
419  SGMatrix<float64_t> temp_matrix(m_d.vlen, m_d.vlen);
420  temp_matrix.zero();
421 #ifdef HAVE_LAPACK
422  cblas_dgemm(
423  CblasRowMajor, CblasNoTrans, CblasNoTrans, m_d.vlen, m_d.vlen,
424  m_d.vlen, 1, m_u.matrix, m_d.vlen, r_matrix.matrix, m_d.vlen, 0,
425  temp_matrix.matrix, m_d.vlen);
426 #else
427  linalg::dgemm<float64_t>(
428  1, m_u, r_matrix, false, false, 0, temp_matrix);
429 #endif
430  r_matrix = temp_matrix;
431  }
432 
434 
435 #ifdef HAVE_LAPACK
436  cblas_dgemv(
437  CblasRowMajor, CblasNoTrans, m_mean.vlen, m_mean.vlen, 1,
438  r_matrix.matrix, m_mean.vlen, random_vec.vector, 1, 0, samp.vector, 1);
439 #else
440  linalg::dgemv<float64_t>(1.0, r_matrix, false, random_vec, 0.0, samp);
441 #endif
442  for (int32_t i = 0; i < m_mean.vlen; i++)
443  samp.vector[i] += m_mean.vector[i];
444 
445  SG_DEBUG("Leaving\n");
446  return samp;
447 }
448 
450 {
451  if (!distribution)
452  return NULL;
453 
454  CGaussian* casted=dynamic_cast<CGaussian*>(distribution);
455  if (!casted)
456  return NULL;
457 
458  /* since an additional reference is returned */
459  SG_REF(casted);
460  return casted;
461 }
SGVector< float64_t > sample()
Definition: Gaussian.cpp:391
float64_t m_constant
Definition: Gaussian.h:235
void set_u(SGMatrix< float64_t > u)
Definition: Gaussian.h:203
SGVector< T > get_row_vector(index_t row) const
Definition: SGMatrix.cpp:1211
virtual void set_features(CFeatures *f)
Definition: Distribution.h:160
virtual SGMatrix< float64_t > get_cov(bool copy_data_for_speed=true)
static T sqrt(T x)
Definition: Math.h:428
void scale(SGVector< T > &a, SGVector< T > &result, T alpha=1)
Gaussian distribution interface.
Definition: Gaussian.h:47
ECovType get_cov_type()
Definition: Gaussian.h:159
virtual bool train(CFeatures *data=NULL)
Definition: Gaussian.cpp:68
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)
virtual float64_t compute_log_PDF(SGVector< float64_t > point)
Definition: Gaussian.cpp:237
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)
ECovType m_cov_type
Definition: Gaussian.h:243
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
virtual int32_t get_dim_feature_space() const =0
spherical covariance
Definition: Gaussian.h:37
SGMatrix< T > clone() const
Definition: SGMatrix.cpp:330
virtual float64_t update_params_em(float64_t *alpha_k, int32_t len)
Definition: Gaussian.cpp:120
virtual SGVector< float64_t > get_mean()
SGMatrix< float64_t > m_u
Definition: Gaussian.h:239
#define ASSERT(x)
Definition: SGIO.h:176
virtual SGVector< float64_t > get_mean()
Definition: Gaussian.cpp:276
static SGVector< float64_t > compute_eigenvectors(SGMatrix< float64_t > matrix)
Definition: SGMatrix.cpp:1006
virtual float64_t get_log_model_parameter(int32_t num_param)
Definition: Gaussian.cpp:100
static CGaussian * obtain_from_generic(CDistribution *distribution)
Definition: Gaussian.cpp:449
virtual void set_cov(SGMatrix< float64_t > cov)
Definition: Gaussian.cpp:289
double float64_t
Definition: common.h:60
ECovType
Definition: Gaussian.h:30
SGVector< float64_t > m_mean
Definition: Gaussian.h:241
index_t num_rows
Definition: SGMatrix.h:495
virtual SGMatrix< float64_t > get_cov()
Definition: Gaussian.cpp:303
#define M_PI
Definition: Math.h:50
SGVector< T > clone() const
Definition: SGVector.cpp:262
static void fill_vector(T *vec, int32_t len, T value)
Definition: SGVector.cpp:279
index_t num_cols
Definition: SGMatrix.h:497
virtual ~CGaussian()
Definition: Gaussian.cpp:64
diagonal covariance
Definition: Gaussian.h:35
virtual float64_t get_log_likelihood_example(int32_t num_example)
Definition: Gaussian.cpp:112
#define SG_DEBUG(...)
Definition: SGIO.h:106
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
int machine_int_t
Definition: common.h:69
The class Features is the base class of all feature objects.
Definition: Features.h:69
static float64_t log(float64_t v)
Definition: Math.h:714
SGVector< float64_t > get_computed_dot_feature_vector(int32_t num)
virtual void set_mean(const SGVector< float64_t > mean)
Definition: Gaussian.cpp:281
#define SG_ADD(...)
Definition: SGObject.h:93
bool has_property(EFeatureProperty p) const
Definition: Features.cpp:295
virtual int32_t get_num_model_parameters()
Definition: Gaussian.cpp:86
index_t vlen
Definition: SGVector.h:571
virtual float64_t get_log_derivative(int32_t num_param, int32_t num_example)
Definition: Gaussian.cpp:106
void set_d(const SGVector< float64_t > d)
Definition: Gaussian.cpp:297
SGVector< float64_t > m_d
Definition: Gaussian.h:237

SHOGUN Machine Learning Toolbox - Documentation