SHOGUN
v3.2.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 * Written (W) 2011 Alesis Novik 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 00009 */ 00010 #include <shogun/lib/config.h> 00011 00012 #ifdef HAVE_LAPACK 00013 00014 #include <shogun/clustering/GMM.h> 00015 #include <shogun/clustering/KMeans.h> 00016 #include <shogun/distance/EuclideanDistance.h> 00017 #include <shogun/base/Parameter.h> 00018 #include <shogun/mathematics/Math.h> 00019 #include <shogun/mathematics/lapack.h> 00020 #include <shogun/labels/MulticlassLabels.h> 00021 #include <shogun/multiclass/KNN.h> 00022 00023 #include <vector> 00024 00025 using namespace shogun; 00026 using namespace std; 00027 00028 CGMM::CGMM() : CDistribution(), m_components(), m_coefficients() 00029 { 00030 register_params(); 00031 } 00032 00033 CGMM::CGMM(int32_t n, ECovType cov_type) : CDistribution(), m_components(), m_coefficients() 00034 { 00035 m_coefficients.vector=SG_MALLOC(float64_t, n); 00036 m_coefficients.vlen=n; 00037 m_components = vector<CGaussian*>(n); 00038 00039 for (int32_t i=0; i<n; i++) 00040 { 00041 m_components[i]=new CGaussian(); 00042 SG_REF(m_components[i]); 00043 m_components[i]->set_cov_type(cov_type); 00044 } 00045 00046 register_params(); 00047 } 00048 00049 CGMM::CGMM(vector<CGaussian*> components, SGVector<float64_t> coefficients, bool copy) : CDistribution() 00050 { 00051 ASSERT(int32_t(components.size())==coefficients.vlen) 00052 00053 if (!copy) 00054 { 00055 m_components=components; 00056 m_coefficients=coefficients; 00057 for (int32_t i=0; i<int32_t(components.size()); i++) 00058 { 00059 SG_REF(m_components[i]); 00060 } 00061 } 00062 else 00063 { 00064 m_coefficients = coefficients; 00065 m_components = vector<CGaussian*>(components.size()); 00066 00067 for (int32_t i=0; i<int32_t(components.size()); i++) 00068 { 00069 m_components[i]=new CGaussian(); 00070 SG_REF(m_components[i]); 00071 m_components[i]->set_cov_type(components[i]->get_cov_type()); 00072 00073 SGVector<float64_t> old_mean=components[i]->get_mean(); 00074 SGVector<float64_t> new_mean(old_mean.vlen); 00075 memcpy(new_mean.vector, old_mean.vector, old_mean.vlen*sizeof(float64_t)); 00076 m_components[i]->set_mean(new_mean); 00077 00078 SGVector<float64_t> old_d=components[i]->get_d(); 00079 SGVector<float64_t> new_d(old_d.vlen); 00080 memcpy(new_d.vector, old_d.vector, old_d.vlen*sizeof(float64_t)); 00081 m_components[i]->set_d(new_d); 00082 00083 if (components[i]->get_cov_type()==FULL) 00084 { 00085 SGMatrix<float64_t> old_u=components[i]->get_u(); 00086 SGMatrix<float64_t> new_u(old_u.num_rows, old_u.num_cols); 00087 memcpy(new_u.matrix, old_u.matrix, old_u.num_rows*old_u.num_cols*sizeof(float64_t)); 00088 m_components[i]->set_u(new_u); 00089 } 00090 00091 m_coefficients[i]=coefficients[i]; 00092 } 00093 } 00094 00095 register_params(); 00096 } 00097 00098 CGMM::~CGMM() 00099 { 00100 if (!m_components.empty()) 00101 cleanup(); 00102 } 00103 00104 void CGMM::cleanup() 00105 { 00106 for (int32_t i = 0; i < int32_t(m_components.size()); i++) 00107 SG_UNREF(m_components[i]); 00108 00109 m_components = vector<CGaussian*>(); 00110 m_coefficients = SGVector<float64_t>(); 00111 } 00112 00113 bool CGMM::train(CFeatures* data) 00114 { 00115 ASSERT(m_components.size() != 0) 00116 00118 if (data) 00119 { 00120 if (!data->has_property(FP_DOT)) 00121 SG_ERROR("Specified features are not of type CDotFeatures\n") 00122 set_features(data); 00123 } 00124 00125 return true; 00126 } 00127 00128 float64_t CGMM::train_em(float64_t min_cov, int32_t max_iter, float64_t min_change) 00129 { 00130 if (!features) 00131 SG_ERROR("No features to train on.\n") 00132 00133 CDotFeatures* dotdata=(CDotFeatures *) features; 00134 int32_t num_vectors=dotdata->get_num_vectors(); 00135 00136 SGMatrix<float64_t> alpha; 00137 00138 /* compute initialization via kmeans if none is present */ 00139 if (m_components[0]->get_mean().vector==NULL) 00140 { 00141 CKMeans* init_k_means=new CKMeans(int32_t(m_components.size()), new CEuclideanDistance()); 00142 init_k_means->train(dotdata); 00143 SGMatrix<float64_t> init_means=init_k_means->get_cluster_centers(); 00144 00145 alpha=alpha_init(init_means); 00146 00147 SG_UNREF(init_k_means); 00148 00149 max_likelihood(alpha, min_cov); 00150 } 00151 else 00152 alpha=SGMatrix<float64_t>(num_vectors,int32_t(m_components.size())); 00153 00154 int32_t iter=0; 00155 float64_t log_likelihood_prev=0; 00156 float64_t log_likelihood_cur=0; 00157 float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00158 float64_t* logPx=SG_MALLOC(float64_t, num_vectors); 00159 //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen); 00160 00161 while (iter<max_iter) 00162 { 00163 log_likelihood_prev=log_likelihood_cur; 00164 log_likelihood_cur=0; 00165 00166 for (int32_t i=0; i<num_vectors; i++) 00167 { 00168 logPx[i]=0; 00169 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00170 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00171 { 00172 logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]); 00173 logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]); 00174 } 00175 00176 logPx[i]=CMath::log(logPx[i]); 00177 log_likelihood_cur+=logPx[i]; 00178 00179 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00180 { 00181 //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i]; 00182 alpha.matrix[i*m_components.size()+j]=CMath::exp(logPxy[i*m_components.size()+j]-logPx[i]); 00183 } 00184 } 00185 00186 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change) 00187 break; 00188 00189 max_likelihood(alpha, min_cov); 00190 00191 iter++; 00192 } 00193 00194 SG_FREE(logPxy); 00195 SG_FREE(logPx); 00196 00197 return log_likelihood_cur; 00198 } 00199 00200 float64_t CGMM::train_smem(int32_t max_iter, int32_t max_cand, float64_t min_cov, int32_t max_em_iter, float64_t min_change) 00201 { 00202 if (!features) 00203 SG_ERROR("No features to train on.\n") 00204 00205 if (m_components.size()<3) 00206 SG_ERROR("Can't run SMEM with less than 3 component mixture model.\n") 00207 00208 CDotFeatures* dotdata=(CDotFeatures *) features; 00209 int32_t num_vectors=dotdata->get_num_vectors(); 00210 00211 float64_t cur_likelihood=train_em(min_cov, max_em_iter, min_change); 00212 00213 int32_t iter=0; 00214 float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00215 float64_t* logPx=SG_MALLOC(float64_t, num_vectors); 00216 float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00217 float64_t* logPostSum=SG_MALLOC(float64_t, m_components.size()); 00218 float64_t* logPostSum2=SG_MALLOC(float64_t, m_components.size()); 00219 float64_t* logPostSumSum=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2); 00220 float64_t* split_crit=SG_MALLOC(float64_t, m_components.size()); 00221 float64_t* merge_crit=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2); 00222 int32_t* split_ind=SG_MALLOC(int32_t, m_components.size()); 00223 int32_t* merge_ind=SG_MALLOC(int32_t, m_components.size()*(m_components.size()-1)/2); 00224 00225 while (iter<max_iter) 00226 { 00227 memset(logPostSum, 0, m_components.size()*sizeof(float64_t)); 00228 memset(logPostSum2, 0, m_components.size()*sizeof(float64_t)); 00229 memset(logPostSumSum, 0, (m_components.size()*(m_components.size()-1)/2)*sizeof(float64_t)); 00230 for (int32_t i=0; i<num_vectors; i++) 00231 { 00232 logPx[i]=0; 00233 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00234 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00235 { 00236 logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]); 00237 logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]); 00238 } 00239 00240 logPx[i]=CMath::log(logPx[i]); 00241 00242 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00243 { 00244 logPost[i*m_components.size()+j]=logPxy[i*m_components.size()+j]-logPx[i]; 00245 logPostSum[j]+=CMath::exp(logPost[i*m_components.size()+j]); 00246 logPostSum2[j]+=CMath::exp(2*logPost[i*m_components.size()+j]); 00247 } 00248 00249 int32_t counter=0; 00250 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00251 { 00252 for (int32_t k=j+1; k<int32_t(m_components.size()); k++) 00253 { 00254 logPostSumSum[counter]+=CMath::exp(logPost[i*m_components.size()+j]+logPost[i*m_components.size()+k]); 00255 counter++; 00256 } 00257 } 00258 } 00259 00260 int32_t counter=0; 00261 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00262 { 00263 logPostSum[i]=CMath::log(logPostSum[i]); 00264 split_crit[i]=0; 00265 split_ind[i]=i; 00266 for (int32_t j=0; j<num_vectors; j++) 00267 { 00268 split_crit[i]+=(logPost[j*m_components.size()+i]-logPostSum[i]-logPxy[j*m_components.size()+i]+CMath::log(m_coefficients[i]))* 00269 (CMath::exp(logPost[j*m_components.size()+i])/CMath::exp(logPostSum[i])); 00270 } 00271 for (int32_t j=i+1; j<int32_t(m_components.size()); j++) 00272 { 00273 merge_crit[counter]=CMath::log(logPostSumSum[counter])-(0.5*CMath::log(logPostSum2[i]))-(0.5*CMath::log(logPostSum2[j])); 00274 merge_ind[counter]=i*m_components.size()+j; 00275 counter++; 00276 } 00277 } 00278 CMath::qsort_backward_index(split_crit, split_ind, int32_t(m_components.size())); 00279 CMath::qsort_backward_index(merge_crit, merge_ind, int32_t(m_components.size()*(m_components.size()-1)/2)); 00280 00281 bool better_found=false; 00282 int32_t candidates_checked=0; 00283 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00284 { 00285 for (int32_t j=0; j<int32_t(m_components.size()*(m_components.size()-1)/2); j++) 00286 { 00287 if (merge_ind[j]/int32_t(m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%m_components.size()) != split_ind[i]) 00288 { 00289 candidates_checked++; 00290 CGMM* candidate=new CGMM(m_components, m_coefficients, true); 00291 candidate->train(features); 00292 candidate->partial_em(split_ind[i], merge_ind[j]/int32_t(m_components.size()), merge_ind[j]%int32_t(m_components.size()), min_cov, max_em_iter, min_change); 00293 float64_t cand_likelihood=candidate->train_em(min_cov, max_em_iter, min_change); 00294 00295 if (cand_likelihood>cur_likelihood) 00296 { 00297 cur_likelihood=cand_likelihood; 00298 set_comp(candidate->get_comp()); 00299 set_coef(candidate->get_coef()); 00300 00301 for (int32_t k=0; k<int32_t(candidate->get_comp().size()); k++) 00302 { 00303 SG_UNREF(candidate->get_comp()[k]); 00304 } 00305 00306 better_found=true; 00307 break; 00308 } 00309 else 00310 delete candidate; 00311 00312 if (candidates_checked>=max_cand) 00313 break; 00314 } 00315 } 00316 if (better_found || candidates_checked>=max_cand) 00317 break; 00318 } 00319 if (!better_found) 00320 break; 00321 iter++; 00322 } 00323 00324 SG_FREE(logPxy); 00325 SG_FREE(logPx); 00326 SG_FREE(logPost); 00327 SG_FREE(split_crit); 00328 SG_FREE(merge_crit); 00329 SG_FREE(logPostSum); 00330 SG_FREE(logPostSum2); 00331 SG_FREE(logPostSumSum); 00332 SG_FREE(split_ind); 00333 SG_FREE(merge_ind); 00334 00335 return cur_likelihood; 00336 } 00337 00338 void CGMM::partial_em(int32_t comp1, int32_t comp2, int32_t comp3, float64_t min_cov, int32_t max_em_iter, float64_t min_change) 00339 { 00340 CDotFeatures* dotdata=(CDotFeatures *) features; 00341 int32_t num_vectors=dotdata->get_num_vectors(); 00342 00343 float64_t* init_logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00344 float64_t* init_logPx=SG_MALLOC(float64_t, num_vectors); 00345 float64_t* init_logPx_fix=SG_MALLOC(float64_t, num_vectors); 00346 float64_t* post_add=SG_MALLOC(float64_t, num_vectors); 00347 00348 for (int32_t i=0; i<num_vectors; i++) 00349 { 00350 init_logPx[i]=0; 00351 init_logPx_fix[i]=0; 00352 00353 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00354 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00355 { 00356 init_logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]); 00357 init_logPx[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]); 00358 if (j!=comp1 && j!=comp2 && j!=comp3) 00359 { 00360 init_logPx_fix[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]); 00361 } 00362 } 00363 00364 init_logPx[i]=CMath::log(init_logPx[i]); 00365 post_add[i]=CMath::log(CMath::exp(init_logPxy[i*m_components.size()+comp1]-init_logPx[i])+ 00366 CMath::exp(init_logPxy[i*m_components.size()+comp2]-init_logPx[i])+ 00367 CMath::exp(init_logPxy[i*m_components.size()+comp3]-init_logPx[i])); 00368 } 00369 00370 vector<CGaussian*> components(3); 00371 SGVector<float64_t> coefficients(3); 00372 components[0]=m_components[comp1]; 00373 components[1]=m_components[comp2]; 00374 components[2]=m_components[comp3]; 00375 coefficients.vector[0]=m_coefficients.vector[comp1]; 00376 coefficients.vector[1]=m_coefficients.vector[comp2]; 00377 coefficients.vector[2]=m_coefficients.vector[comp3]; 00378 float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2]; 00379 00380 int32_t dim_n=components[0]->get_mean().vlen; 00381 00382 float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]); 00383 float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]); 00384 00385 float64_t noise_mag=SGVector<float64_t>::twonorm(components[0]->get_mean().vector, dim_n)*0.1/ 00386 CMath::sqrt((float64_t)dim_n); 00387 00388 SGVector<float64_t>::add(components[1]->get_mean().vector, alpha1, components[1]->get_mean().vector, alpha2, 00389 components[2]->get_mean().vector, dim_n); 00390 00391 for (int32_t i=0; i<dim_n; i++) 00392 { 00393 components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag; 00394 components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag; 00395 } 00396 00397 coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2]; 00398 coefficients.vector[2]=coefficients.vector[0]*0.5; 00399 coefficients.vector[0]=coefficients.vector[0]*0.5; 00400 00401 if (components[0]->get_cov_type()==FULL) 00402 { 00403 SGMatrix<float64_t> c1=components[1]->get_cov(); 00404 SGMatrix<float64_t> c2=components[2]->get_cov(); 00405 SGVector<float64_t>::add(c1.matrix, alpha1, c1.matrix, alpha2, c2.matrix, dim_n*dim_n); 00406 00407 components[1]->set_d(SGVector<float64_t>(SGMatrix<float64_t>::compute_eigenvectors(c1.matrix, dim_n, dim_n), dim_n)); 00408 components[1]->set_u(c1); 00409 00410 float64_t new_d=0; 00411 for (int32_t i=0; i<dim_n; i++) 00412 { 00413 new_d+=CMath::log(components[0]->get_d().vector[i]); 00414 for (int32_t j=0; j<dim_n; j++) 00415 { 00416 if (i==j) 00417 { 00418 components[0]->get_u().matrix[i*dim_n+j]=1; 00419 components[2]->get_u().matrix[i*dim_n+j]=1; 00420 } 00421 else 00422 { 00423 components[0]->get_u().matrix[i*dim_n+j]=0; 00424 components[2]->get_u().matrix[i*dim_n+j]=0; 00425 } 00426 } 00427 } 00428 new_d=CMath::exp(new_d*(1./dim_n)); 00429 for (int32_t i=0; i<dim_n; i++) 00430 { 00431 components[0]->get_d().vector[i]=new_d; 00432 components[2]->get_d().vector[i]=new_d; 00433 } 00434 } 00435 else if(components[0]->get_cov_type()==DIAG) 00436 { 00437 SGVector<float64_t>::add(components[1]->get_d().vector, alpha1, components[1]->get_d().vector, 00438 alpha2, components[2]->get_d().vector, dim_n); 00439 00440 float64_t new_d=0; 00441 for (int32_t i=0; i<dim_n; i++) 00442 { 00443 new_d+=CMath::log(components[0]->get_d().vector[i]); 00444 } 00445 new_d=CMath::exp(new_d*(1./dim_n)); 00446 for (int32_t i=0; i<dim_n; i++) 00447 { 00448 components[0]->get_d().vector[i]=new_d; 00449 components[2]->get_d().vector[i]=new_d; 00450 } 00451 } 00452 else if(components[0]->get_cov_type()==SPHERICAL) 00453 { 00454 components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+ 00455 alpha2*components[2]->get_d().vector[0]; 00456 00457 components[2]->get_d().vector[0]=components[0]->get_d().vector[0]; 00458 } 00459 00460 CGMM* partial_candidate=new CGMM(components, coefficients); 00461 partial_candidate->train(features); 00462 00463 float64_t log_likelihood_prev=0; 00464 float64_t log_likelihood_cur=0; 00465 int32_t iter=0; 00466 SGMatrix<float64_t> alpha(num_vectors, 3); 00467 float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*3); 00468 float64_t* logPx=SG_MALLOC(float64_t, num_vectors); 00469 //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen); 00470 00471 while (iter<max_em_iter) 00472 { 00473 log_likelihood_prev=log_likelihood_cur; 00474 log_likelihood_cur=0; 00475 00476 for (int32_t i=0; i<num_vectors; i++) 00477 { 00478 logPx[i]=0; 00479 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00480 for (int32_t j=0; j<3; j++) 00481 { 00482 logPxy[i*3+j]=components[j]->compute_log_PDF(v)+CMath::log(coefficients[j]); 00483 logPx[i]+=CMath::exp(logPxy[i*3+j]); 00484 } 00485 00486 logPx[i]=CMath::log(logPx[i]+init_logPx_fix[i]); 00487 log_likelihood_cur+=logPx[i]; 00488 00489 for (int32_t j=0; j<3; j++) 00490 { 00491 //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i]; 00492 alpha.matrix[i*3+j]=CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]); 00493 } 00494 } 00495 00496 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change) 00497 break; 00498 00499 partial_candidate->max_likelihood(alpha, min_cov); 00500 partial_candidate->get_coef().vector[0]=partial_candidate->get_coef().vector[0]*coef_sum; 00501 partial_candidate->get_coef().vector[1]=partial_candidate->get_coef().vector[1]*coef_sum; 00502 partial_candidate->get_coef().vector[2]=partial_candidate->get_coef().vector[2]*coef_sum; 00503 iter++; 00504 } 00505 00506 m_coefficients.vector[comp1]=coefficients.vector[0]; 00507 m_coefficients.vector[comp2]=coefficients.vector[1]; 00508 m_coefficients.vector[comp3]=coefficients.vector[2]; 00509 00510 delete partial_candidate; 00511 SG_FREE(logPxy); 00512 SG_FREE(logPx); 00513 SG_FREE(init_logPxy); 00514 SG_FREE(init_logPx); 00515 SG_FREE(init_logPx_fix); 00516 SG_FREE(post_add); 00517 } 00518 00519 void CGMM::max_likelihood(SGMatrix<float64_t> alpha, float64_t min_cov) 00520 { 00521 CDotFeatures* dotdata=(CDotFeatures *) features; 00522 int32_t num_dim=dotdata->get_dim_feature_space(); 00523 00524 float64_t alpha_sum; 00525 float64_t alpha_sum_sum=0; 00526 float64_t* mean_sum; 00527 float64_t* cov_sum=NULL; 00528 00529 for (int32_t i=0; i<alpha.num_cols; i++) 00530 { 00531 alpha_sum=0; 00532 mean_sum=SG_MALLOC(float64_t, num_dim); 00533 memset(mean_sum, 0, num_dim*sizeof(float64_t)); 00534 00535 for (int32_t j=0; j<alpha.num_rows; j++) 00536 { 00537 alpha_sum+=alpha.matrix[j*alpha.num_cols+i]; 00538 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(j); 00539 SGVector<float64_t>::add(mean_sum, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, mean_sum, v.vlen); 00540 } 00541 00542 for (int32_t j=0; j<num_dim; j++) 00543 mean_sum[j]/=alpha_sum; 00544 00545 m_components[i]->set_mean(SGVector<float64_t>(mean_sum, num_dim)); 00546 00547 ECovType cov_type=m_components[i]->get_cov_type(); 00548 00549 if (cov_type==FULL) 00550 { 00551 cov_sum=SG_MALLOC(float64_t, num_dim*num_dim); 00552 memset(cov_sum, 0, num_dim*num_dim*sizeof(float64_t)); 00553 } 00554 else if(cov_type==DIAG) 00555 { 00556 cov_sum=SG_MALLOC(float64_t, num_dim); 00557 memset(cov_sum, 0, num_dim*sizeof(float64_t)); 00558 } 00559 else if(cov_type==SPHERICAL) 00560 { 00561 cov_sum=SG_MALLOC(float64_t, 1); 00562 cov_sum[0]=0; 00563 } 00564 00565 for (int32_t j=0; j<alpha.num_rows; j++) 00566 { 00567 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(j); 00568 SGVector<float64_t>::add(v.vector, 1, v.vector, -1, mean_sum, v.vlen); 00569 00570 switch (cov_type) 00571 { 00572 case FULL: 00573 cblas_dger(CblasRowMajor, num_dim, num_dim, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, v.vector, 00574 1, (double*) cov_sum, num_dim); 00575 00576 break; 00577 case DIAG: 00578 for (int32_t k=0; k<num_dim; k++) 00579 cov_sum[k]+=v.vector[k]*v.vector[k]*alpha.matrix[j*alpha.num_cols+i]; 00580 00581 break; 00582 case SPHERICAL: 00583 float64_t temp=0; 00584 00585 for (int32_t k=0; k<num_dim; k++) 00586 temp+=v.vector[k]*v.vector[k]; 00587 00588 cov_sum[0]+=temp*alpha.matrix[j*alpha.num_cols+i]; 00589 break; 00590 } 00591 } 00592 00593 switch (cov_type) 00594 { 00595 case FULL: 00596 for (int32_t j=0; j<num_dim*num_dim; j++) 00597 cov_sum[j]/=alpha_sum; 00598 00599 float64_t* d0; 00600 d0=SGMatrix<float64_t>::compute_eigenvectors(cov_sum, num_dim, num_dim); 00601 for (int32_t j=0; j<num_dim; j++) 00602 d0[j]=CMath::max(min_cov, d0[j]); 00603 00604 m_components[i]->set_d(SGVector<float64_t>(d0, num_dim)); 00605 m_components[i]->set_u(SGMatrix<float64_t>(cov_sum, num_dim, num_dim)); 00606 00607 break; 00608 case DIAG: 00609 for (int32_t j=0; j<num_dim; j++) 00610 { 00611 cov_sum[j]/=alpha_sum; 00612 cov_sum[j]=CMath::max(min_cov, cov_sum[j]); 00613 } 00614 00615 m_components[i]->set_d(SGVector<float64_t>(cov_sum, num_dim)); 00616 00617 break; 00618 case SPHERICAL: 00619 cov_sum[0]/=alpha_sum*num_dim; 00620 cov_sum[0]=CMath::max(min_cov, cov_sum[0]); 00621 00622 m_components[i]->set_d(SGVector<float64_t>(cov_sum, 1)); 00623 00624 break; 00625 } 00626 00627 m_coefficients.vector[i]=alpha_sum; 00628 alpha_sum_sum+=alpha_sum; 00629 } 00630 00631 for (int32_t i=0; i<alpha.num_cols; i++) 00632 m_coefficients.vector[i]/=alpha_sum_sum; 00633 } 00634 00635 int32_t CGMM::get_num_model_parameters() 00636 { 00637 return 1; 00638 } 00639 00640 float64_t CGMM::get_log_model_parameter(int32_t num_param) 00641 { 00642 ASSERT(num_param==1) 00643 00644 return CMath::log(m_components.size()); 00645 } 00646 00647 index_t CGMM::get_num_components() const 00648 { 00649 return m_components.size(); 00650 } 00651 00652 CDistribution* CGMM::get_component(index_t index) const 00653 { 00654 return m_components[index]; 00655 } 00656 00657 float64_t CGMM::get_log_derivative(int32_t num_param, int32_t num_example) 00658 { 00659 SG_NOTIMPLEMENTED 00660 return 0; 00661 } 00662 00663 float64_t CGMM::get_log_likelihood_example(int32_t num_example) 00664 { 00665 SG_NOTIMPLEMENTED 00666 return 1; 00667 } 00668 00669 float64_t CGMM::get_likelihood_example(int32_t num_example) 00670 { 00671 float64_t result=0; 00672 00673 ASSERT(features); 00674 ASSERT(features->get_feature_class() == C_DENSE); 00675 ASSERT(features->get_feature_type() == F_DREAL); 00676 00677 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00678 { 00679 SGVector<float64_t> point= ((CDenseFeatures<float64_t>*) features)->get_feature_vector(num_example); 00680 result+=CMath::exp(m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i])); 00681 } 00682 00683 return result; 00684 } 00685 00686 SGVector<float64_t> CGMM::get_nth_mean(int32_t num) 00687 { 00688 ASSERT(num<int32_t(m_components.size())) 00689 return m_components[num]->get_mean(); 00690 } 00691 00692 void CGMM::set_nth_mean(SGVector<float64_t> mean, int32_t num) 00693 { 00694 ASSERT(num<int32_t(m_components.size())) 00695 m_components[num]->set_mean(mean); 00696 } 00697 00698 SGMatrix<float64_t> CGMM::get_nth_cov(int32_t num) 00699 { 00700 ASSERT(num<int32_t(m_components.size())) 00701 return m_components[num]->get_cov(); 00702 } 00703 00704 void CGMM::set_nth_cov(SGMatrix<float64_t> cov, int32_t num) 00705 { 00706 ASSERT(num<int32_t(m_components.size())) 00707 m_components[num]->set_cov(cov); 00708 } 00709 00710 SGVector<float64_t> CGMM::get_coef() 00711 { 00712 return m_coefficients; 00713 } 00714 00715 void CGMM::set_coef(const SGVector<float64_t> coefficients) 00716 { 00717 m_coefficients=coefficients; 00718 } 00719 00720 vector<CGaussian*> CGMM::get_comp() 00721 { 00722 return m_components; 00723 } 00724 00725 void CGMM::set_comp(vector<CGaussian*> components) 00726 { 00727 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00728 { 00729 SG_UNREF(m_components[i]); 00730 } 00731 00732 m_components=components; 00733 00734 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00735 { 00736 SG_REF(m_components[i]); 00737 } 00738 } 00739 00740 SGMatrix<float64_t> CGMM::alpha_init(SGMatrix<float64_t> init_means) 00741 { 00742 CDotFeatures* dotdata=(CDotFeatures *) features; 00743 int32_t num_vectors=dotdata->get_num_vectors(); 00744 00745 SGVector<float64_t> label_num(init_means.num_cols); 00746 00747 for (int32_t i=0; i<init_means.num_cols; i++) 00748 label_num.vector[i]=i; 00749 00750 CKNN* knn=new CKNN(1, new CEuclideanDistance(), new CMulticlassLabels(label_num)); 00751 knn->train(new CDenseFeatures<float64_t>(init_means)); 00752 CMulticlassLabels* init_labels=(CMulticlassLabels*) knn->apply(features); 00753 00754 SGMatrix<float64_t> alpha(num_vectors, int32_t(m_components.size())); 00755 memset(alpha.matrix, 0, num_vectors*m_components.size()*sizeof(float64_t)); 00756 00757 for (int32_t i=0; i<num_vectors; i++) 00758 alpha.matrix[i*m_components.size()+init_labels->get_int_label(i)]=1; 00759 00760 SG_UNREF(init_labels); 00761 00762 return alpha; 00763 } 00764 00765 SGVector<float64_t> CGMM::sample() 00766 { 00767 REQUIRE(m_components.size()>0, "Number of mixture components is %d but " 00768 "must be positive\n", m_components.size()); 00769 float64_t rand_num=CMath::random(float64_t(0), float64_t(1)); 00770 float64_t cum_sum=0; 00771 for (int32_t i=0; i<m_coefficients.vlen; i++) 00772 { 00773 cum_sum+=m_coefficients.vector[i]; 00774 if (cum_sum>=rand_num) 00775 { 00776 SG_DEBUG("Sampling from mixture component %d\n", i); 00777 return m_components[i]->sample(); 00778 } 00779 } 00780 return m_components[m_coefficients.vlen-1]->sample(); 00781 } 00782 00783 SGVector<float64_t> CGMM::cluster(SGVector<float64_t> point) 00784 { 00785 SGVector<float64_t> answer(m_components.size()+1); 00786 answer.vector[m_components.size()]=0; 00787 00788 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00789 { 00790 answer.vector[i]=m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]); 00791 answer.vector[m_components.size()]+=CMath::exp(answer.vector[i]); 00792 } 00793 answer.vector[m_components.size()]=CMath::log(answer.vector[m_components.size()]); 00794 00795 return answer; 00796 } 00797 00798 void CGMM::register_params() 00799 { 00800 //TODO serialization broken 00801 //m_parameters->add((SGVector<CSGObject*>*) &m_components, "m_components", "Mixture components"); 00802 m_parameters->add(&m_coefficients, "m_coefficients", "Mixture coefficients."); 00803 } 00804 00805 #endif