SHOGUN
v2.0.0
|
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