SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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  if (cov.do_free)
41  cov.free_matrix();
42 }
43 
45 {
47  switch (m_cov_type)
48  {
49  case FULL:
50  case DIAG:
51  for (int32_t i=0; i<m_d.vlen; i++)
53  break;
54  case SPHERICAL:
56  break;
57  }
58 }
59 
61 {
65 }
66 
68 {
69  // init features with data if necessary and assure type is correct
70  if (data)
71  {
72  if (!data->has_property(FP_DOT))
73  SG_ERROR("Specified features are not of type CDotFeatures\n");
74  set_features(data);
75  }
76  CDotFeatures* dotdata=(CDotFeatures *) data;
77 
79 
80  m_mean=dotdata->get_mean();
81  SGMatrix<float64_t> cov=dotdata->get_cov();
82 
83  decompose_cov(cov);
84  cov.destroy_matrix();
85 
86  init();
87 
88  return true;
89 }
90 
92 {
93  switch (m_cov_type)
94  {
95  case FULL:
97  case DIAG:
98  return m_d.vlen+m_mean.vlen;
99  case SPHERICAL:
100  return 1+m_mean.vlen;
101  }
102  return 0;
103 }
104 
106 {
108  return 0;
109 }
110 
111 float64_t CGaussian::get_log_derivative(int32_t num_param, int32_t num_example)
112 {
114  return 0;
115 }
116 
118 {
120  SGVector<float64_t> v=((CDotFeatures *)features)->get_computed_dot_feature_vector(num_example);
121  float64_t answer=compute_log_PDF(v);
122  v.free_vector();
123  return answer;
124 }
125 
127 {
129  ASSERT(point.vlen == m_mean.vlen);
130  float64_t* difference=SG_MALLOC(float64_t, m_mean.vlen);
131  memcpy(difference, point.vector, sizeof(float64_t)*m_mean.vlen);
132 
133  for (int32_t i = 0; i < m_mean.vlen; i++)
134  difference[i] -= m_mean.vector[i];
135 
136  float64_t answer=m_constant;
137 
138  if (m_cov_type==FULL)
139  {
140  float64_t* temp_holder=SG_MALLOC(float64_t, m_d.vlen);
141  cblas_dgemv(CblasRowMajor, CblasNoTrans, m_d.vlen, m_d.vlen,
142  1, m_u.matrix, m_d.vlen, difference, 1, 0, temp_holder, 1);
143 
144  for (int32_t i=0; i<m_d.vlen; i++)
145  answer+=temp_holder[i]*temp_holder[i]/m_d.vector[i];
146 
147  SG_FREE(temp_holder);
148  }
149  else if (m_cov_type==DIAG)
150  {
151  for (int32_t i=0; i<m_mean.vlen; i++)
152  answer+=difference[i]*difference[i]/m_d.vector[i];
153  }
154  else
155  {
156  for (int32_t i=0; i<m_mean.vlen; i++)
157  answer+=difference[i]*difference[i]/m_d.vector[0];
158  }
159 
160  SG_FREE(difference);
161 
162  return -0.5*answer;
163 }
164 
166 {
168  memset(cov, 0, sizeof(float64_t)*m_mean.vlen*m_mean.vlen);
169 
170  if (m_cov_type==FULL)
171  {
172  if (!m_u.matrix)
173  SG_ERROR("Unitary matrix not set\n");
174 
175  float64_t* temp_holder=SG_MALLOC(float64_t, m_d.vlen*m_d.vlen);
176  float64_t* diag_holder=SG_MALLOC(float64_t, m_d.vlen*m_d.vlen);
177  memset(diag_holder, 0, sizeof(float64_t)*m_d.vlen*m_d.vlen);
178  for(int32_t i=0; i<m_d.vlen; i++)
179  diag_holder[i*m_d.vlen+i]=m_d.vector[i];
180 
181  cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans,
183  diag_holder, m_d.vlen, 0, temp_holder, m_d.vlen);
184  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
185  m_d.vlen, m_d.vlen, m_d.vlen, 1, temp_holder, m_d.vlen,
186  m_u.matrix, m_d.vlen, 0, cov, m_d.vlen);
187 
188  SG_FREE(diag_holder);
189  SG_FREE(temp_holder);
190  }
191  else if (m_cov_type==DIAG)
192  {
193  for (int32_t i=0; i<m_d.vlen; i++)
194  cov[i*m_d.vlen+i]=m_d.vector[i];
195  }
196  else
197  {
198  for (int32_t i=0; i<m_mean.vlen; i++)
199  cov[i*m_mean.vlen+i]=m_d.vector[0];
200  }
201  return SGMatrix<float64_t>(cov, m_mean.vlen, m_mean.vlen, false);//fix needed
202 }
203 
204 void CGaussian::register_params()
205 {
206  m_parameters->add(&m_u, "m_u", "Unitary matrix.");
207  m_parameters->add(&m_d, "m_d", "Diagonal.");
208  m_parameters->add(&m_mean, "m_mean", "Mean.");
209  m_parameters->add(&m_constant, "m_constant", "Constant part.");
210  m_parameters->add((machine_int_t*)&m_cov_type, "m_cov_type", "Covariance type.");
211 }
212 
213 void CGaussian::decompose_cov(SGMatrix<float64_t> cov)
214 {
216  switch (m_cov_type)
217  {
218  case FULL:
221  memcpy(m_u.matrix, cov.matrix, sizeof(float64_t)*cov.num_rows*cov.num_rows);
222 
224  m_d.vlen=cov.num_rows;
225  m_u.num_rows=cov.num_rows;
226  m_u.num_cols=cov.num_rows;
227  break;
228  case DIAG:
230 
231  for (int32_t i=0; i<cov.num_rows; i++)
232  m_d.vector[i]=cov.matrix[i*cov.num_rows+i];
233 
234  m_d.vlen=cov.num_rows;
235  break;
236  case SPHERICAL:
238 
239  m_d.vector[0]=cov.matrix[0];
240  m_d.vlen=1;
241  break;
242  }
243 }
244 
246 {
248  memset(r_matrix, 0, m_mean.vlen*m_mean.vlen*sizeof(float64_t));
249 
250  switch (m_cov_type)
251  {
252  case FULL:
253  case DIAG:
254  for (int32_t i=0; i<m_mean.vlen; i++)
255  r_matrix[i*m_mean.vlen+i]=CMath::sqrt(m_d.vector[i]);
256 
257  break;
258  case SPHERICAL:
259  for (int32_t i=0; i<m_mean.vlen; i++)
260  r_matrix[i*m_mean.vlen+i]=CMath::sqrt(m_d.vector[0]);
261 
262  break;
263  }
264 
265  float64_t* random_vec=SG_MALLOC(float64_t, m_mean.vlen);
266 
267  for (int32_t i=0; i<m_mean.vlen; i++)
268  random_vec[i]=CMath::randn_double();
269 
270  if (m_cov_type==FULL)
271  {
272  float64_t* temp_matrix=SG_MALLOC(float64_t, m_d.vlen*m_d.vlen);
273  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
275  r_matrix, m_d.vlen, 0, temp_matrix, m_d.vlen);
276  SG_FREE(r_matrix);
277  r_matrix=temp_matrix;
278  }
279 
281 
282  cblas_dgemv(CblasRowMajor, CblasNoTrans, m_mean.vlen, m_mean.vlen,
283  1, r_matrix, m_mean.vlen, random_vec, 1, 0, samp, 1);
284 
285  for (int32_t i=0; i<m_mean.vlen; i++)
286  samp[i]+=m_mean.vector[i];
287 
288  SG_FREE(random_vec);
289  SG_FREE(r_matrix);
290 
291  return SGVector<float64_t>(samp, m_mean.vlen, false);//fix needed
292 }
293 
294 #endif

SHOGUN Machine Learning Toolbox - Documentation