SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 #include <shogun/lib/config.h>
11 
12 #ifdef HAVE_LAPACK
13 
16 #include <shogun/base/Parameter.h>
17 
18 using namespace shogun;
19 
20 CGaussian::CGaussian() : CDistribution(), m_constant(0), m_d(), m_u(), m_mean(), m_cov_type(FULL)
21 {
22  register_params();
23 }
24 
26  ECovType cov_type) : CDistribution(), m_d(), m_u(), m_cov_type(cov_type)
27 {
28  ASSERT(mean.vlen==cov.num_rows)
29  ASSERT(cov.num_rows==cov.num_cols)
30 
31  m_mean=mean;
32 
33  if (cov.num_rows==1)
35 
36  decompose_cov(cov);
37  init();
38  register_params();
39 }
40 
42 {
44  switch (m_cov_type)
45  {
46  case FULL:
47  case DIAG:
48  for (int32_t i=0; i<m_d.vlen; i++)
50  break;
51  case SPHERICAL:
53  break;
54  }
55 }
56 
58 {
59 }
60 
62 {
63  // init features with data if necessary and assure type is correct
64  if (data)
65  {
66  if (!data->has_property(FP_DOT))
67  SG_ERROR("Specified features are not of type CDotFeatures\n")
68  set_features(data);
69  }
70 
71  CDotFeatures* dotdata=(CDotFeatures *) data;
72  m_mean=dotdata->get_mean();
73  SGMatrix<float64_t> cov=dotdata->get_cov();
74  decompose_cov(cov);
75  init();
76  return true;
77 }
78 
80 {
81  switch (m_cov_type)
82  {
83  case FULL:
85  case DIAG:
86  return m_d.vlen+m_mean.vlen;
87  case SPHERICAL:
88  return 1+m_mean.vlen;
89  }
90  return 0;
91 }
92 
94 {
96  return 0;
97 }
98 
99 float64_t CGaussian::get_log_derivative(int32_t num_param, int32_t num_example)
100 {
102  return 0;
103 }
104 
106 {
108  SGVector<float64_t> v=((CDotFeatures *)features)->get_computed_dot_feature_vector(num_example);
109  float64_t answer=compute_log_PDF(v);
110  return answer;
111 }
112 
114 {
116  ASSERT(point.vlen == m_mean.vlen)
117  float64_t* difference=SG_MALLOC(float64_t, m_mean.vlen);
118  memcpy(difference, point.vector, sizeof(float64_t)*m_mean.vlen);
119 
120  for (int32_t i = 0; i < m_mean.vlen; i++)
121  difference[i] -= m_mean.vector[i];
122 
123  float64_t answer=m_constant;
124 
125  if (m_cov_type==FULL)
126  {
127  float64_t* temp_holder=SG_MALLOC(float64_t, m_d.vlen);
128  cblas_dgemv(CblasRowMajor, CblasNoTrans, m_d.vlen, m_d.vlen,
129  1, m_u.matrix, m_d.vlen, difference, 1, 0, temp_holder, 1);
130 
131  for (int32_t i=0; i<m_d.vlen; i++)
132  answer+=temp_holder[i]*temp_holder[i]/m_d.vector[i];
133 
134  SG_FREE(temp_holder);
135  }
136  else if (m_cov_type==DIAG)
137  {
138  for (int32_t i=0; i<m_mean.vlen; i++)
139  answer+=difference[i]*difference[i]/m_d.vector[i];
140  }
141  else
142  {
143  for (int32_t i=0; i<m_mean.vlen; i++)
144  answer+=difference[i]*difference[i]/m_d.vector[0];
145  }
146 
147  SG_FREE(difference);
148 
149  return -0.5*answer;
150 }
151 
153 {
154  return m_mean;
155 }
156 
158 {
159  if (mean.vlen==1)
161 
162  m_mean=mean;
163 }
164 
166 {
167  ASSERT(cov.num_rows==cov.num_cols)
168  ASSERT(cov.num_rows==m_mean.vlen)
169  decompose_cov(cov);
170  init();
171 }
172 
174 {
175  m_d = d;
176  init();
177 }
178 
180 {
181  float64_t* cov=SG_MALLOC(float64_t, m_mean.vlen*m_mean.vlen);
182  memset(cov, 0, sizeof(float64_t)*m_mean.vlen*m_mean.vlen);
183 
184  if (m_cov_type==FULL)
185  {
186  if (!m_u.matrix)
187  SG_ERROR("Unitary matrix not set\n")
188 
189  float64_t* temp_holder=SG_MALLOC(float64_t, m_d.vlen*m_d.vlen);
190  float64_t* diag_holder=SG_MALLOC(float64_t, m_d.vlen*m_d.vlen);
191  memset(diag_holder, 0, sizeof(float64_t)*m_d.vlen*m_d.vlen);
192  for(int32_t i=0; i<m_d.vlen; i++)
193  diag_holder[i*m_d.vlen+i]=m_d.vector[i];
194 
195  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
197  diag_holder, m_d.vlen, 0, temp_holder, m_d.vlen);
198  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
199  m_d.vlen, m_d.vlen, m_d.vlen, 1, temp_holder, m_d.vlen,
200  m_u.matrix, m_d.vlen, 0, cov, m_d.vlen);
201 
202  SG_FREE(diag_holder);
203  SG_FREE(temp_holder);
204  }
205  else if (m_cov_type==DIAG)
206  {
207  for (int32_t i=0; i<m_d.vlen; i++)
208  cov[i*m_d.vlen+i]=m_d.vector[i];
209  }
210  else
211  {
212  for (int32_t i=0; i<m_mean.vlen; i++)
213  cov[i*m_mean.vlen+i]=m_d.vector[0];
214  }
215  return SGMatrix<float64_t>(cov, m_mean.vlen, m_mean.vlen, false);//fix needed
216 }
217 
218 void CGaussian::register_params()
219 {
220  m_parameters->add(&m_u, "m_u", "Unitary matrix.");
221  m_parameters->add(&m_d, "m_d", "Diagonal.");
222  m_parameters->add(&m_mean, "m_mean", "Mean.");
223  m_parameters->add(&m_constant, "m_constant", "Constant part.");
224  m_parameters->add((machine_int_t*)&m_cov_type, "m_cov_type", "Covariance type.");
225 }
226 
227 void CGaussian::decompose_cov(SGMatrix<float64_t> cov)
228 {
229  switch (m_cov_type)
230  {
231  case FULL:
233  memcpy(m_u.matrix, cov.matrix, sizeof(float64_t)*cov.num_rows*cov.num_rows);
234 
236  m_d.vlen=cov.num_rows;
237  m_u.num_rows=cov.num_rows;
238  m_u.num_cols=cov.num_rows;
239  break;
240  case DIAG:
241  m_d.vector=SG_MALLOC(float64_t, cov.num_rows);
242 
243  for (int32_t i=0; i<cov.num_rows; i++)
244  m_d.vector[i]=cov.matrix[i*cov.num_rows+i];
245 
246  m_d.vlen=cov.num_rows;
247  break;
248  case SPHERICAL:
249  m_d.vector=SG_MALLOC(float64_t, 1);
250 
251  m_d.vector[0]=cov.matrix[0];
252  m_d.vlen=1;
253  break;
254  }
255 }
256 
258 {
259  SG_DEBUG("Entering\n");
260  float64_t* r_matrix=SG_MALLOC(float64_t, m_mean.vlen*m_mean.vlen);
261  memset(r_matrix, 0, m_mean.vlen*m_mean.vlen*sizeof(float64_t));
262 
263  switch (m_cov_type)
264  {
265  case FULL:
266  case DIAG:
267  for (int32_t i=0; i<m_mean.vlen; i++)
268  r_matrix[i*m_mean.vlen+i]=CMath::sqrt(m_d.vector[i]);
269 
270  break;
271  case SPHERICAL:
272  for (int32_t i=0; i<m_mean.vlen; i++)
273  r_matrix[i*m_mean.vlen+i]=CMath::sqrt(m_d.vector[0]);
274 
275  break;
276  }
277 
278  float64_t* random_vec=SG_MALLOC(float64_t, m_mean.vlen);
279 
280  for (int32_t i=0; i<m_mean.vlen; i++)
281  random_vec[i]=CMath::randn_double();
282 
283  if (m_cov_type==FULL)
284  {
285  float64_t* temp_matrix=SG_MALLOC(float64_t, m_d.vlen*m_d.vlen);
286  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
288  r_matrix, m_d.vlen, 0, temp_matrix, m_d.vlen);
289  SG_FREE(r_matrix);
290  r_matrix=temp_matrix;
291  }
292 
293  float64_t* samp=SG_MALLOC(float64_t, m_mean.vlen);
294 
295  cblas_dgemv(CblasRowMajor, CblasNoTrans, m_mean.vlen, m_mean.vlen,
296  1, r_matrix, m_mean.vlen, random_vec, 1, 0, samp, 1);
297 
298  for (int32_t i=0; i<m_mean.vlen; i++)
299  samp[i]+=m_mean.vector[i];
300 
301  SG_FREE(random_vec);
302  SG_FREE(r_matrix);
303 
304  SG_DEBUG("Leaving\n");
305  return SGVector<float64_t>(samp, m_mean.vlen, false);//fix needed
306 }
307 
309 {
310  if (!distribution)
311  return NULL;
312 
313  CGaussian* casted=dynamic_cast<CGaussian*>(distribution);
314  if (!casted)
315  return NULL;
316 
317  /* since an additional reference is returned */
318  SG_REF(casted);
319  return casted;
320 }
321 
322 #endif // HAVE_LAPACK

SHOGUN Machine Learning Toolbox - Documentation