00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <shogun/lib/config.h>
00015 #ifdef HAVE_EIGEN3
00016 #ifdef HAVE_LAPACK
00017
00018 #include <shogun/io/SGIO.h>
00019 #include <shogun/regression/GaussianProcessRegression.h>
00020 #include <shogun/mathematics/lapack.h>
00021 #include <shogun/mathematics/Math.h>
00022 #include <shogun/kernel/Kernel.h>
00023 #include <shogun/labels/RegressionLabels.h>
00024 #include <shogun/features/CombinedFeatures.h>
00025
00026 using namespace shogun;
00027
00028 CGaussianProcessRegression::CGaussianProcessRegression()
00029 : CMachine()
00030 {
00031 init();
00032 }
00033
00034 CGaussianProcessRegression::CGaussianProcessRegression(CInferenceMethod* inf, CFeatures* data, CLabels* lab)
00035 : CMachine()
00036 {
00037 init();
00038
00039 set_labels(lab);
00040 set_features(data);
00041 set_method(inf);
00042 }
00043
00044 void CGaussianProcessRegression::init()
00045 {
00046
00047 m_features = NULL;
00048 m_method = NULL;
00049 m_data = NULL;
00050 m_return = GP_RETURN_MEANS;
00051
00052 SG_ADD((CSGObject**) &m_features, "features", "Feature object.",
00053 MS_NOT_AVAILABLE);
00054 SG_ADD((CSGObject**) &m_method, "inference_method", "Inference Method.",
00055 MS_AVAILABLE);
00056 }
00057
00058 void CGaussianProcessRegression::update_kernel_matrices()
00059 {
00060 CKernel* kernel = NULL;
00061
00062 if (m_method)
00063 kernel = m_method->get_kernel();
00064
00065 if (kernel)
00066 {
00067 float64_t m_scale = m_method->get_scale();
00068
00069 CFeatures* latent_features = m_method->get_latent_features();
00070
00071 if (latent_features)
00072 kernel->init(latent_features, m_data);
00073 else
00074 kernel->init(m_data, m_data);
00075
00076
00077 m_k_trts = kernel->get_kernel_matrix();
00078
00079 for (index_t i = 0; i < m_k_trts.num_rows; i++)
00080 {
00081 for (index_t j = 0; j < m_k_trts.num_cols; j++)
00082 m_k_trts(i,j) *= (m_scale*m_scale);
00083 }
00084
00085 kernel->init(m_data, m_data);
00086
00087 m_k_tsts = kernel->get_kernel_matrix();
00088
00089 for (index_t i = 0; i < m_k_tsts.num_rows; i++)
00090 {
00091 for (index_t j = 0; j < m_k_tsts.num_cols; j++)
00092 m_k_tsts(i,j) *= (m_scale*m_scale);
00093 }
00094
00095 kernel->remove_lhs_and_rhs();
00096
00097 SG_UNREF(kernel);
00098 SG_UNREF(latent_features);
00099 }
00100 }
00101
00102 CRegressionLabels* CGaussianProcessRegression::apply_regression(CFeatures* data)
00103 {
00104
00105 if (data)
00106 {
00107 if(data->get_feature_class() == C_COMBINED)
00108 {
00109 CDotFeatures* feat =
00110 (CDotFeatures*)((CCombinedFeatures*)data)->
00111 get_first_feature_obj();
00112
00113 if (!feat->has_property(FP_DOT))
00114 SG_ERROR("Specified features are not of type CFeatures\n");
00115
00116 if (feat->get_feature_class() != C_DENSE)
00117 SG_ERROR("Expected Simple Features\n");
00118
00119 if (feat->get_feature_type() != F_DREAL)
00120 SG_ERROR("Expected Real Features\n");
00121
00122 SG_UNREF(feat);
00123 }
00124
00125 else
00126 {
00127 if (!data->has_property(FP_DOT))
00128 SG_ERROR("Specified features are not of type CFeatures\n");
00129
00130 if (data->get_feature_class() != C_DENSE)
00131 SG_ERROR("Expected Simple Features\n");
00132
00133 if (data->get_feature_type() != F_DREAL)
00134 SG_ERROR("Expected Real Features\n");
00135 }
00136
00137 SG_UNREF(m_data);
00138 SG_REF(data);
00139 m_data = (CFeatures*)data;
00140 update_kernel_matrices();
00141 }
00142
00143 else if (!m_data)
00144 SG_ERROR("No testing features!\n");
00145
00146 else if (update_parameter_hash())
00147 update_kernel_matrices();
00148
00149 if (m_return == GP_RETURN_COV)
00150 {
00151 CRegressionLabels* result =
00152 new CRegressionLabels(get_covariance_vector());
00153
00154 return result;
00155 }
00156
00157 if (m_return == GP_RETURN_MEANS)
00158 {
00159 CRegressionLabels* result =
00160 new CRegressionLabels(get_mean_vector());
00161
00162 return result;
00163 }
00164
00165 else
00166 {
00167
00168 SGVector<float64_t> mean_vector = get_mean_vector();
00169 SGVector<float64_t> cov_vector = get_covariance_vector();
00170
00171 index_t size = mean_vector.vlen+cov_vector.vlen;
00172
00173 SGVector<float64_t> result_vector(size);
00174
00175 for (index_t i = 0; i < size; i++)
00176 {
00177 if (i < mean_vector.vlen)
00178 result_vector[i] = mean_vector[i];
00179 else
00180 result_vector[i] = cov_vector[i-mean_vector.vlen];
00181 }
00182
00183 CRegressionLabels* result =
00184 new CRegressionLabels(result_vector);
00185
00186 return result;
00187 }
00188
00189 }
00190
00191 bool CGaussianProcessRegression::train_machine(CFeatures* data)
00192 {
00193 return false;
00194 }
00195
00196
00197 SGVector<float64_t> CGaussianProcessRegression::get_mean_vector()
00198 {
00199
00200 SGVector<float64_t> m_alpha = m_method->get_alpha();
00201
00202 CMeanFunction* mean_function = m_method->get_mean();
00203
00204 SGMatrix<float64_t> features;
00205 if(m_data->get_feature_class() == C_COMBINED)
00206 {
00207 features = ((CDotFeatures*)((CCombinedFeatures*)m_data)->
00208 get_first_feature_obj())->
00209 get_computed_dot_feature_matrix();
00210 }
00211
00212 else
00213 {
00214 features = ((CDotFeatures*)m_data)->
00215 get_computed_dot_feature_matrix();
00216 }
00217
00218 if (!mean_function)
00219 SG_ERROR("Mean function is NULL!\n");
00220
00221
00222 SGVector<float64_t> means = mean_function->get_mean_vector(features);
00223
00224 SGVector< float64_t > result_vector(features.num_cols);
00225
00226
00227 cblas_dgemv(CblasColMajor, CblasTrans, m_k_trts.num_rows,
00228 m_alpha.vlen, 1.0, m_k_trts.matrix,
00229 m_k_trts.num_cols, m_alpha.vector, 1, 0.0,
00230 result_vector.vector, 1);
00231
00232 for (index_t i = 0; i < result_vector.vlen; i++)
00233 result_vector[i] += means[i];
00234
00235 CLikelihoodModel* lik = m_method->get_model();
00236
00237 result_vector = lik->evaluate_means(result_vector);
00238
00239 SG_UNREF(lik);
00240 SG_UNREF(mean_function);
00241
00242 return result_vector;
00243 }
00244
00245
00246 SGVector<float64_t> CGaussianProcessRegression::get_covariance_vector()
00247 {
00248
00249 if (!m_data)
00250 SG_ERROR("No testing features!\n");
00251
00252 SGVector<float64_t> diagonal = m_method->get_diagonal_vector();
00253
00254 if (diagonal.vlen > 0)
00255 {
00256 SGVector<float64_t> diagonal2(m_data->get_num_vectors());
00257
00258 SGMatrix<float64_t> temp1(m_data->get_num_vectors(), diagonal.vlen);
00259
00260 SGMatrix<float64_t> m_L = m_method->get_cholesky();
00261
00262 SGMatrix<float64_t> temp2(m_L.num_rows, m_L.num_cols);
00263
00264 for (index_t i = 0; i < diagonal.vlen; i++)
00265 {
00266 for (index_t j = 0; j < m_data->get_num_vectors(); j++)
00267 temp1(j,i) = diagonal[i]*m_k_trts(j,i);
00268 }
00269
00270 for (index_t i = 0; i < diagonal2.vlen; i++)
00271 diagonal2[i] = 0;
00272
00273 memcpy(temp2.matrix, m_L.matrix,
00274 m_L.num_cols*m_L.num_rows*sizeof(float64_t));
00275
00276
00277 temp2.transpose_matrix(temp2.matrix, temp2.num_rows, temp2.num_cols);
00278
00279 SGVector<int32_t> ipiv(temp2.num_cols);
00280
00281
00282 clapack_dgetrf(CblasColMajor, temp2.num_rows, temp2.num_cols,
00283 temp2.matrix, temp2.num_cols, ipiv.vector);
00284
00285 clapack_dgetrs(CblasColMajor, CblasNoTrans,
00286 temp2.num_rows, temp1.num_cols, temp2.matrix,
00287 temp2.num_cols, ipiv.vector, temp1.matrix,
00288 temp1.num_cols);
00289
00290 for (index_t i = 0; i < temp1.num_rows; i++)
00291 {
00292 for (index_t j = 0; j < temp1.num_cols; j++)
00293 temp1(i,j) = temp1(i,j)*temp1(i,j);
00294 }
00295
00296 for (index_t i = 0; i < temp1.num_cols; i++)
00297 {
00298 diagonal2[i] = 0;
00299
00300 for (index_t j = 0; j < temp1.num_rows; j++)
00301 diagonal2[i] += temp1(j,i);
00302 }
00303
00304
00305 SGVector<float64_t> result(m_k_tsts.num_cols);
00306
00307
00308 for (index_t i = 0; i < m_k_tsts.num_cols; i++)
00309 result[i] = m_k_tsts(i,i) - diagonal2[i];
00310
00311 CLikelihoodModel* lik = m_method->get_model();
00312
00313 SG_UNREF(lik);
00314
00315 return lik->evaluate_variances(result);
00316 }
00317
00318 else
00319 {
00320 SGMatrix<float64_t> m_L = m_method->get_cholesky();
00321
00322 SGMatrix<float64_t> temp1(m_k_trts.num_rows, m_k_trts.num_cols);
00323 SGMatrix<float64_t> temp2(m_L.num_rows, m_L.num_cols);
00324
00325 memcpy(temp1.matrix, m_k_trts.matrix,
00326 m_k_trts.num_cols*m_k_trts.num_rows*sizeof(float64_t));
00327
00328 memcpy(temp2.matrix, m_L.matrix,
00329 m_L.num_cols*m_L.num_rows*sizeof(float64_t));
00330
00331 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m_L.num_rows,
00332 m_k_trts.num_cols, m_L.num_rows, 1.0, m_L.matrix, m_L.num_cols,
00333 m_k_trts.matrix, m_k_trts.num_cols, 0.0, temp1.matrix,
00334 temp1.num_cols);
00335
00336 for (index_t i = 0; i < temp1.num_rows; i++)
00337 {
00338 for (index_t j = 0; j < temp1.num_cols; j++)
00339 temp1(i,j) *= m_k_trts(i,j);
00340 }
00341
00342 SGVector<float64_t> temp3(temp2.num_cols);
00343
00344 for (index_t i = 0; i < temp1.num_cols; i++)
00345 {
00346 float64_t sum = 0;
00347 for (index_t j = 0; j < temp1.num_rows; j++)
00348 sum += temp1(j,i);
00349 temp3[i] = sum;
00350 }
00351
00352 SGVector<float64_t> result(m_k_tsts.num_cols);
00353
00354 for (index_t i = 0; i < m_k_tsts.num_cols; i++)
00355 result[i] = m_k_tsts(i,i) + temp3[i];
00356
00357 CLikelihoodModel* lik = m_method->get_model();
00358
00359 SG_UNREF(lik);
00360
00361 return lik->evaluate_variances(result);
00362 }
00363 }
00364
00365
00366 CGaussianProcessRegression::~CGaussianProcessRegression()
00367 {
00368 SG_UNREF(m_features);
00369 SG_UNREF(m_method);
00370 SG_UNREF(m_data);
00371 }
00372
00373 void CGaussianProcessRegression::set_kernel(CKernel* k)
00374 {
00375 m_method->set_kernel(k);
00376 }
00377
00378 bool CGaussianProcessRegression::load(FILE* srcfile)
00379 {
00380 SG_SET_LOCALE_C;
00381 SG_RESET_LOCALE;
00382 return false;
00383 }
00384
00385 CKernel* CGaussianProcessRegression::get_kernel()
00386 {
00387 return m_method->get_kernel();
00388 }
00389
00390 bool CGaussianProcessRegression::save(FILE* dstfile)
00391 {
00392 SG_SET_LOCALE_C;
00393 SG_RESET_LOCALE;
00394 return false;
00395 }
00396 #endif
00397 #endif