GaussianProcessRegression.cpp

Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Copyright (C) 2012 Jacob Walker
00008  * 
00009  * Code adapted from Gaussian Process Machine Learning Toolbox
00010  * http://www.gaussianprocess.org/gpml/code/matlab/doc/
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         //K(X_test, X_train)
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     //Here we multiply K*^t by alpha to receive the mean predictions.
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         //Solve L^T V = K(Train,Test)*Diagonal Vector Entries for V
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         //Subtract V from K(Test,Test) to get covariances.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation