SHOGUN  v3.2.0
 全部  命名空间 文件 函数 变量 类型定义 枚举 枚举值 友元 宏定义  
Gaussian.cpp
浏览该文件的文档.
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
SGVector< float64_t > sample()
Definition: Gaussian.cpp:257
float64_t m_constant
Definition: Gaussian.h:225
virtual void set_features(CFeatures *f)
Definition: Distribution.h:157
Gaussian distribution interface.
Definition: Gaussian.h:46
virtual bool train(CFeatures *data=NULL)
Definition: Gaussian.cpp:61
static float64_t randn_double()
Definition: Math.h:607
#define SG_ERROR(...)
Definition: SGIO.h:131
#define SG_NOTIMPLEMENTED
Definition: SGIO.h:141
Parameter * m_parameters
Definition: SGObject.h:482
virtual float64_t compute_log_PDF(SGVector< float64_t > point)
Definition: Gaussian.cpp:113
Base class Distribution from which all methods implementing a distribution are derived.
Definition: Distribution.h:41
ECovType m_cov_type
Definition: Gaussian.h:233
Features that support dot products among other operations.
Definition: DotFeatures.h:41
full covariance
Definition: Gaussian.h:32
spherical covariance
Definition: Gaussian.h:36
virtual SGVector< float64_t > get_mean()
void add(bool *param, const char *name, const char *description="")
Definition: Parameter.cpp:27
SGMatrix< float64_t > m_u
Definition: Gaussian.h:229
#define ASSERT(x)
Definition: SGIO.h:203
virtual SGVector< float64_t > get_mean()
Definition: Gaussian.cpp:152
static SGVector< float64_t > compute_eigenvectors(SGMatrix< float64_t > matrix)
Definition: SGMatrix.cpp:781
virtual float64_t get_log_model_parameter(int32_t num_param)
Definition: Gaussian.cpp:93
static CGaussian * obtain_from_generic(CDistribution *distribution)
Definition: Gaussian.cpp:308
virtual void set_cov(SGMatrix< float64_t > cov)
Definition: Gaussian.cpp:165
double float64_t
Definition: common.h:48
ECovType
Definition: Gaussian.h:29
#define SG_REF(x)
Definition: SGRefObject.h:34
SGVector< float64_t > m_mean
Definition: Gaussian.h:231
index_t num_rows
Definition: SGMatrix.h:301
virtual SGMatrix< float64_t > get_cov()
Definition: Gaussian.cpp:179
#define M_PI
workaround for log2 being a define on cygwin
Definition: Math.h:56
index_t num_cols
Definition: SGMatrix.h:303
virtual ~CGaussian()
Definition: Gaussian.cpp:57
diagonal covariance
Definition: Gaussian.h:34
virtual float64_t get_log_likelihood_example(int32_t num_example)
Definition: Gaussian.cpp:105
#define SG_DEBUG(...)
Definition: SGIO.h:109
int machine_int_t
Definition: common.h:57
The class Features is the base class of all feature objects.
Definition: Features.h:62
static float64_t log(float64_t v)
Definition: Math.h:420
virtual void set_mean(const SGVector< float64_t > mean)
Definition: Gaussian.cpp:157
static float32_t sqrt(float32_t x)
x^0.5
Definition: Math.h:245
bool has_property(EFeatureProperty p) const
Definition: Features.cpp:292
virtual int32_t get_num_model_parameters()
Definition: Gaussian.cpp:79
index_t vlen
Definition: SGVector.h:706
virtual SGMatrix< float64_t > get_cov()
virtual float64_t get_log_derivative(int32_t num_param, int32_t num_example)
Definition: Gaussian.cpp:99
void set_d(const SGVector< float64_t > d)
Definition: Gaussian.cpp:173
SGVector< float64_t > m_d
Definition: Gaussian.h:227

SHOGUN Machine Learning Toolbox - Documentation