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 * 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 if (m_components[0]->get_mean().vector==NULL) 00139 { 00140 CKMeans* init_k_means=new CKMeans(int32_t(m_components.size()), new CEuclideanDistance()); 00141 init_k_means->train(dotdata); 00142 SGMatrix<float64_t> init_means=init_k_means->get_cluster_centers(); 00143 00144 alpha=alpha_init(init_means); 00145 00146 SG_UNREF(init_k_means); 00147 00148 max_likelihood(alpha, min_cov); 00149 } 00150 else 00151 alpha=SGMatrix<float64_t>(num_vectors,int32_t(m_components.size())); 00152 00153 int32_t iter=0; 00154 float64_t log_likelihood_prev=0; 00155 float64_t log_likelihood_cur=0; 00156 float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00157 float64_t* logPx=SG_MALLOC(float64_t, num_vectors); 00158 //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen); 00159 00160 while (iter<max_iter) 00161 { 00162 log_likelihood_prev=log_likelihood_cur; 00163 log_likelihood_cur=0; 00164 00165 for (int32_t i=0; i<num_vectors; i++) 00166 { 00167 logPx[i]=0; 00168 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00169 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00170 { 00171 logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]); 00172 logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]); 00173 } 00174 00175 logPx[i]=CMath::log(logPx[i]); 00176 log_likelihood_cur+=logPx[i]; 00177 00178 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00179 { 00180 //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i]; 00181 alpha.matrix[i*m_components.size()+j]=CMath::exp(logPxy[i*m_components.size()+j]-logPx[i]); 00182 } 00183 } 00184 00185 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change) 00186 break; 00187 00188 max_likelihood(alpha, min_cov); 00189 00190 iter++; 00191 } 00192 00193 SG_FREE(logPxy); 00194 SG_FREE(logPx); 00195 00196 return log_likelihood_cur; 00197 } 00198 00199 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) 00200 { 00201 if (!features) 00202 SG_ERROR("No features to train on.\n"); 00203 00204 if (m_components.size()<3) 00205 SG_ERROR("Can't run SMEM with less than 3 component mixture model.\n"); 00206 00207 CDotFeatures* dotdata=(CDotFeatures *) features; 00208 int32_t num_vectors=dotdata->get_num_vectors(); 00209 00210 float64_t cur_likelihood=train_em(min_cov, max_em_iter, min_change); 00211 00212 int32_t iter=0; 00213 float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00214 float64_t* logPx=SG_MALLOC(float64_t, num_vectors); 00215 float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00216 float64_t* logPostSum=SG_MALLOC(float64_t, m_components.size()); 00217 float64_t* logPostSum2=SG_MALLOC(float64_t, m_components.size()); 00218 float64_t* logPostSumSum=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2); 00219 float64_t* split_crit=SG_MALLOC(float64_t, m_components.size()); 00220 float64_t* merge_crit=SG_MALLOC(float64_t, m_components.size()*(m_components.size()-1)/2); 00221 int32_t* split_ind=SG_MALLOC(int32_t, m_components.size()); 00222 int32_t* merge_ind=SG_MALLOC(int32_t, m_components.size()*(m_components.size()-1)/2); 00223 00224 while (iter<max_iter) 00225 { 00226 memset(logPostSum, 0, m_components.size()*sizeof(float64_t)); 00227 memset(logPostSum2, 0, m_components.size()*sizeof(float64_t)); 00228 memset(logPostSumSum, 0, (m_components.size()*(m_components.size()-1)/2)*sizeof(float64_t)); 00229 for (int32_t i=0; i<num_vectors; i++) 00230 { 00231 logPx[i]=0; 00232 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00233 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00234 { 00235 logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]); 00236 logPx[i]+=CMath::exp(logPxy[i*m_components.size()+j]); 00237 } 00238 00239 logPx[i]=CMath::log(logPx[i]); 00240 00241 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00242 { 00243 logPost[i*m_components.size()+j]=logPxy[i*m_components.size()+j]-logPx[i]; 00244 logPostSum[j]+=CMath::exp(logPost[i*m_components.size()+j]); 00245 logPostSum2[j]+=CMath::exp(2*logPost[i*m_components.size()+j]); 00246 } 00247 00248 int32_t counter=0; 00249 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00250 { 00251 for (int32_t k=j+1; k<int32_t(m_components.size()); k++) 00252 { 00253 logPostSumSum[counter]+=CMath::exp(logPost[i*m_components.size()+j]+logPost[i*m_components.size()+k]); 00254 counter++; 00255 } 00256 } 00257 } 00258 00259 int32_t counter=0; 00260 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00261 { 00262 logPostSum[i]=CMath::log(logPostSum[i]); 00263 split_crit[i]=0; 00264 split_ind[i]=i; 00265 for (int32_t j=0; j<num_vectors; j++) 00266 { 00267 split_crit[i]+=(logPost[j*m_components.size()+i]-logPostSum[i]-logPxy[j*m_components.size()+i]+CMath::log(m_coefficients[i]))* 00268 (CMath::exp(logPost[j*m_components.size()+i])/CMath::exp(logPostSum[i])); 00269 } 00270 for (int32_t j=i+1; j<int32_t(m_components.size()); j++) 00271 { 00272 merge_crit[counter]=CMath::log(logPostSumSum[counter])-(0.5*CMath::log(logPostSum2[i]))-(0.5*CMath::log(logPostSum2[j])); 00273 merge_ind[counter]=i*m_components.size()+j; 00274 counter++; 00275 } 00276 } 00277 CMath::qsort_backward_index(split_crit, split_ind, int32_t(m_components.size())); 00278 CMath::qsort_backward_index(merge_crit, merge_ind, int32_t(m_components.size()*(m_components.size()-1)/2)); 00279 00280 bool better_found=false; 00281 int32_t candidates_checked=0; 00282 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00283 { 00284 for (int32_t j=0; j<int32_t(m_components.size()*(m_components.size()-1)/2); j++) 00285 { 00286 if (merge_ind[j]/int32_t(m_components.size()) != split_ind[i] && int32_t(merge_ind[j]%m_components.size()) != split_ind[i]) 00287 { 00288 candidates_checked++; 00289 CGMM* candidate=new CGMM(m_components, m_coefficients, true); 00290 candidate->train(features); 00291 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); 00292 float64_t cand_likelihood=candidate->train_em(min_cov, max_em_iter, min_change); 00293 00294 if (cand_likelihood>cur_likelihood) 00295 { 00296 cur_likelihood=cand_likelihood; 00297 set_comp(candidate->get_comp()); 00298 set_coef(candidate->get_coef()); 00299 00300 for (int32_t k=0; k<int32_t(candidate->get_comp().size()); k++) 00301 { 00302 SG_UNREF(candidate->get_comp()[k]); 00303 } 00304 00305 better_found=true; 00306 break; 00307 } 00308 else 00309 delete candidate; 00310 00311 if (candidates_checked>=max_cand) 00312 break; 00313 } 00314 } 00315 if (better_found || candidates_checked>=max_cand) 00316 break; 00317 } 00318 if (!better_found) 00319 break; 00320 iter++; 00321 } 00322 00323 SG_FREE(logPxy); 00324 SG_FREE(logPx); 00325 SG_FREE(logPost); 00326 SG_FREE(split_crit); 00327 SG_FREE(merge_crit); 00328 SG_FREE(logPostSum); 00329 SG_FREE(logPostSum2); 00330 SG_FREE(logPostSumSum); 00331 SG_FREE(split_ind); 00332 SG_FREE(merge_ind); 00333 00334 return cur_likelihood; 00335 } 00336 00337 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) 00338 { 00339 CDotFeatures* dotdata=(CDotFeatures *) features; 00340 int32_t num_vectors=dotdata->get_num_vectors(); 00341 00342 float64_t* init_logPxy=SG_MALLOC(float64_t, num_vectors*m_components.size()); 00343 float64_t* init_logPx=SG_MALLOC(float64_t, num_vectors); 00344 float64_t* init_logPx_fix=SG_MALLOC(float64_t, num_vectors); 00345 float64_t* post_add=SG_MALLOC(float64_t, num_vectors); 00346 00347 for (int32_t i=0; i<num_vectors; i++) 00348 { 00349 init_logPx[i]=0; 00350 init_logPx_fix[i]=0; 00351 00352 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00353 for (int32_t j=0; j<int32_t(m_components.size()); j++) 00354 { 00355 init_logPxy[i*m_components.size()+j]=m_components[j]->compute_log_PDF(v)+CMath::log(m_coefficients[j]); 00356 init_logPx[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]); 00357 if (j!=comp1 && j!=comp2 && j!=comp3) 00358 { 00359 init_logPx_fix[i]+=CMath::exp(init_logPxy[i*m_components.size()+j]); 00360 } 00361 } 00362 00363 init_logPx[i]=CMath::log(init_logPx[i]); 00364 post_add[i]=CMath::log(CMath::exp(init_logPxy[i*m_components.size()+comp1]-init_logPx[i])+ 00365 CMath::exp(init_logPxy[i*m_components.size()+comp2]-init_logPx[i])+ 00366 CMath::exp(init_logPxy[i*m_components.size()+comp3]-init_logPx[i])); 00367 } 00368 00369 vector<CGaussian*> components(3); 00370 SGVector<float64_t> coefficients(3); 00371 components[0]=m_components[comp1]; 00372 components[1]=m_components[comp2]; 00373 components[2]=m_components[comp3]; 00374 coefficients.vector[0]=m_coefficients.vector[comp1]; 00375 coefficients.vector[1]=m_coefficients.vector[comp2]; 00376 coefficients.vector[2]=m_coefficients.vector[comp3]; 00377 float64_t coef_sum=coefficients.vector[0]+coefficients.vector[1]+coefficients.vector[2]; 00378 00379 int32_t dim_n=components[0]->get_mean().vlen; 00380 00381 float64_t alpha1=coefficients.vector[1]/(coefficients.vector[1]+coefficients.vector[2]); 00382 float64_t alpha2=coefficients.vector[2]/(coefficients.vector[1]+coefficients.vector[2]); 00383 00384 float64_t noise_mag=SGVector<float64_t>::twonorm(components[0]->get_mean().vector, dim_n)*0.1/ 00385 CMath::sqrt((float64_t)dim_n); 00386 00387 SGVector<float64_t>::add(components[1]->get_mean().vector, alpha1, components[1]->get_mean().vector, alpha2, 00388 components[2]->get_mean().vector, dim_n); 00389 00390 for (int32_t i=0; i<dim_n; i++) 00391 { 00392 components[2]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag; 00393 components[0]->get_mean().vector[i]=components[0]->get_mean().vector[i]+CMath::randn_double()*noise_mag; 00394 } 00395 00396 coefficients.vector[1]=coefficients.vector[1]+coefficients.vector[2]; 00397 coefficients.vector[2]=coefficients.vector[0]*0.5; 00398 coefficients.vector[0]=coefficients.vector[0]*0.5; 00399 00400 if (components[0]->get_cov_type()==FULL) 00401 { 00402 SGMatrix<float64_t> c1=components[1]->get_cov(); 00403 SGMatrix<float64_t> c2=components[2]->get_cov(); 00404 SGVector<float64_t>::add(c1.matrix, alpha1, c1.matrix, alpha2, c2.matrix, dim_n*dim_n); 00405 00406 components[1]->set_d(SGVector<float64_t>(SGMatrix<float64_t>::compute_eigenvectors(c1.matrix, dim_n, dim_n), dim_n)); 00407 components[1]->set_u(c1); 00408 00409 float64_t new_d=0; 00410 for (int32_t i=0; i<dim_n; i++) 00411 { 00412 new_d+=CMath::log(components[0]->get_d().vector[i]); 00413 for (int32_t j=0; j<dim_n; j++) 00414 { 00415 if (i==j) 00416 { 00417 components[0]->get_u().matrix[i*dim_n+j]=1; 00418 components[2]->get_u().matrix[i*dim_n+j]=1; 00419 } 00420 else 00421 { 00422 components[0]->get_u().matrix[i*dim_n+j]=0; 00423 components[2]->get_u().matrix[i*dim_n+j]=0; 00424 } 00425 } 00426 } 00427 new_d=CMath::exp(new_d*(1./dim_n)); 00428 for (int32_t i=0; i<dim_n; i++) 00429 { 00430 components[0]->get_d().vector[i]=new_d; 00431 components[2]->get_d().vector[i]=new_d; 00432 } 00433 } 00434 else if(components[0]->get_cov_type()==DIAG) 00435 { 00436 SGVector<float64_t>::add(components[1]->get_d().vector, alpha1, components[1]->get_d().vector, 00437 alpha2, components[2]->get_d().vector, dim_n); 00438 00439 float64_t new_d=0; 00440 for (int32_t i=0; i<dim_n; i++) 00441 { 00442 new_d+=CMath::log(components[0]->get_d().vector[i]); 00443 } 00444 new_d=CMath::exp(new_d*(1./dim_n)); 00445 for (int32_t i=0; i<dim_n; i++) 00446 { 00447 components[0]->get_d().vector[i]=new_d; 00448 components[2]->get_d().vector[i]=new_d; 00449 } 00450 } 00451 else if(components[0]->get_cov_type()==SPHERICAL) 00452 { 00453 components[1]->get_d().vector[0]=alpha1*components[1]->get_d().vector[0]+ 00454 alpha2*components[2]->get_d().vector[0]; 00455 00456 components[2]->get_d().vector[0]=components[0]->get_d().vector[0]; 00457 } 00458 00459 CGMM* partial_candidate=new CGMM(components, coefficients); 00460 partial_candidate->train(features); 00461 00462 float64_t log_likelihood_prev=0; 00463 float64_t log_likelihood_cur=0; 00464 int32_t iter=0; 00465 SGMatrix<float64_t> alpha(num_vectors, 3); 00466 float64_t* logPxy=SG_MALLOC(float64_t, num_vectors*3); 00467 float64_t* logPx=SG_MALLOC(float64_t, num_vectors); 00468 //float64_t* logPost=SG_MALLOC(float64_t, num_vectors*m_components.vlen); 00469 00470 while (iter<max_em_iter) 00471 { 00472 log_likelihood_prev=log_likelihood_cur; 00473 log_likelihood_cur=0; 00474 00475 for (int32_t i=0; i<num_vectors; i++) 00476 { 00477 logPx[i]=0; 00478 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(i); 00479 for (int32_t j=0; j<3; j++) 00480 { 00481 logPxy[i*3+j]=components[j]->compute_log_PDF(v)+CMath::log(coefficients[j]); 00482 logPx[i]+=CMath::exp(logPxy[i*3+j]); 00483 } 00484 00485 logPx[i]=CMath::log(logPx[i]+init_logPx_fix[i]); 00486 log_likelihood_cur+=logPx[i]; 00487 00488 for (int32_t j=0; j<3; j++) 00489 { 00490 //logPost[i*m_components.vlen+j]=logPxy[i*m_components.vlen+j]-logPx[i]; 00491 alpha.matrix[i*3+j]=CMath::exp(logPxy[i*3+j]-logPx[i]+post_add[i]); 00492 } 00493 } 00494 00495 if (iter>0 && log_likelihood_cur-log_likelihood_prev<min_change) 00496 break; 00497 00498 partial_candidate->max_likelihood(alpha, min_cov); 00499 partial_candidate->get_coef().vector[0]=partial_candidate->get_coef().vector[0]*coef_sum; 00500 partial_candidate->get_coef().vector[1]=partial_candidate->get_coef().vector[1]*coef_sum; 00501 partial_candidate->get_coef().vector[2]=partial_candidate->get_coef().vector[2]*coef_sum; 00502 iter++; 00503 } 00504 00505 m_coefficients.vector[comp1]=coefficients.vector[0]; 00506 m_coefficients.vector[comp2]=coefficients.vector[1]; 00507 m_coefficients.vector[comp3]=coefficients.vector[2]; 00508 00509 delete partial_candidate; 00510 SG_FREE(logPxy); 00511 SG_FREE(logPx); 00512 SG_FREE(init_logPxy); 00513 SG_FREE(init_logPx); 00514 SG_FREE(init_logPx_fix); 00515 SG_FREE(post_add); 00516 } 00517 00518 void CGMM::max_likelihood(SGMatrix<float64_t> alpha, float64_t min_cov) 00519 { 00520 CDotFeatures* dotdata=(CDotFeatures *) features; 00521 int32_t num_dim=dotdata->get_dim_feature_space(); 00522 00523 float64_t alpha_sum; 00524 float64_t alpha_sum_sum=0; 00525 float64_t* mean_sum; 00526 float64_t* cov_sum=NULL; 00527 00528 for (int32_t i=0; i<alpha.num_cols; i++) 00529 { 00530 alpha_sum=0; 00531 mean_sum=SG_MALLOC(float64_t, num_dim); 00532 memset(mean_sum, 0, num_dim*sizeof(float64_t)); 00533 00534 for (int32_t j=0; j<alpha.num_rows; j++) 00535 { 00536 alpha_sum+=alpha.matrix[j*alpha.num_cols+i]; 00537 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(j); 00538 SGVector<float64_t>::add(mean_sum, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, mean_sum, v.vlen); 00539 } 00540 00541 for (int32_t j=0; j<num_dim; j++) 00542 mean_sum[j]/=alpha_sum; 00543 00544 m_components[i]->set_mean(SGVector<float64_t>(mean_sum, num_dim)); 00545 00546 ECovType cov_type=m_components[i]->get_cov_type(); 00547 00548 if (cov_type==FULL) 00549 { 00550 cov_sum=SG_MALLOC(float64_t, num_dim*num_dim); 00551 memset(cov_sum, 0, num_dim*num_dim*sizeof(float64_t)); 00552 } 00553 else if(cov_type==DIAG) 00554 { 00555 cov_sum=SG_MALLOC(float64_t, num_dim); 00556 memset(cov_sum, 0, num_dim*sizeof(float64_t)); 00557 } 00558 else if(cov_type==SPHERICAL) 00559 { 00560 cov_sum=SG_MALLOC(float64_t, 1); 00561 cov_sum[0]=0; 00562 } 00563 00564 for (int32_t j=0; j<alpha.num_rows; j++) 00565 { 00566 SGVector<float64_t> v=dotdata->get_computed_dot_feature_vector(j); 00567 SGVector<float64_t>::add(v.vector, 1, v.vector, -1, mean_sum, v.vlen); 00568 00569 switch (cov_type) 00570 { 00571 case FULL: 00572 cblas_dger(CblasRowMajor, num_dim, num_dim, alpha.matrix[j*alpha.num_cols+i], v.vector, 1, v.vector, 00573 1, (double*) cov_sum, num_dim); 00574 00575 break; 00576 case DIAG: 00577 for (int32_t k=0; k<num_dim; k++) 00578 cov_sum[k]+=v.vector[k]*v.vector[k]*alpha.matrix[j*alpha.num_cols+i]; 00579 00580 break; 00581 case SPHERICAL: 00582 float64_t temp=0; 00583 00584 for (int32_t k=0; k<num_dim; k++) 00585 temp+=v.vector[k]*v.vector[k]; 00586 00587 cov_sum[0]+=temp*alpha.matrix[j*alpha.num_cols+i]; 00588 break; 00589 } 00590 } 00591 00592 switch (cov_type) 00593 { 00594 case FULL: 00595 for (int32_t j=0; j<num_dim*num_dim; j++) 00596 cov_sum[j]/=alpha_sum; 00597 00598 float64_t* d0; 00599 d0=SGMatrix<float64_t>::compute_eigenvectors(cov_sum, num_dim, num_dim); 00600 for (int32_t j=0; j<num_dim; j++) 00601 d0[j]=CMath::max(min_cov, d0[j]); 00602 00603 m_components[i]->set_d(SGVector<float64_t>(d0, num_dim)); 00604 m_components[i]->set_u(SGMatrix<float64_t>(cov_sum, num_dim, num_dim)); 00605 00606 break; 00607 case DIAG: 00608 for (int32_t j=0; j<num_dim; j++) 00609 { 00610 cov_sum[j]/=alpha_sum; 00611 cov_sum[j]=CMath::max(min_cov, cov_sum[j]); 00612 } 00613 00614 m_components[i]->set_d(SGVector<float64_t>(cov_sum, num_dim)); 00615 00616 break; 00617 case SPHERICAL: 00618 cov_sum[0]/=alpha_sum*num_dim; 00619 cov_sum[0]=CMath::max(min_cov, cov_sum[0]); 00620 00621 m_components[i]->set_d(SGVector<float64_t>(cov_sum, 1)); 00622 00623 break; 00624 } 00625 00626 m_coefficients.vector[i]=alpha_sum; 00627 alpha_sum_sum+=alpha_sum; 00628 } 00629 00630 for (int32_t i=0; i<alpha.num_cols; i++) 00631 m_coefficients.vector[i]/=alpha_sum_sum; 00632 } 00633 00634 int32_t CGMM::get_num_model_parameters() 00635 { 00636 return 1; 00637 } 00638 00639 float64_t CGMM::get_log_model_parameter(int32_t num_param) 00640 { 00641 ASSERT(num_param==1); 00642 00643 return CMath::log(m_components.size()); 00644 } 00645 00646 float64_t CGMM::get_log_derivative(int32_t num_param, int32_t num_example) 00647 { 00648 SG_NOTIMPLEMENTED; 00649 return 0; 00650 } 00651 00652 float64_t CGMM::get_log_likelihood_example(int32_t num_example) 00653 { 00654 SG_NOTIMPLEMENTED; 00655 return 1; 00656 } 00657 00658 float64_t CGMM::get_likelihood_example(int32_t num_example) 00659 { 00660 return CMath::exp(get_log_likelihood_example(num_example)); 00661 } 00662 00663 SGVector<float64_t> CGMM::get_nth_mean(int32_t num) 00664 { 00665 ASSERT(num<int32_t(m_components.size())); 00666 return m_components[num]->get_mean(); 00667 } 00668 00669 void CGMM::set_nth_mean(SGVector<float64_t> mean, int32_t num) 00670 { 00671 ASSERT(num<int32_t(m_components.size())); 00672 m_components[num]->set_mean(mean); 00673 } 00674 00675 SGMatrix<float64_t> CGMM::get_nth_cov(int32_t num) 00676 { 00677 ASSERT(num<int32_t(m_components.size())); 00678 return m_components[num]->get_cov(); 00679 } 00680 00681 void CGMM::set_nth_cov(SGMatrix<float64_t> cov, int32_t num) 00682 { 00683 ASSERT(num<int32_t(m_components.size())); 00684 m_components[num]->set_cov(cov); 00685 } 00686 00687 SGVector<float64_t> CGMM::get_coef() 00688 { 00689 return m_coefficients; 00690 } 00691 00692 void CGMM::set_coef(const SGVector<float64_t> coefficients) 00693 { 00694 m_coefficients=coefficients; 00695 } 00696 00697 vector<CGaussian*> CGMM::get_comp() 00698 { 00699 return m_components; 00700 } 00701 00702 void CGMM::set_comp(vector<CGaussian*> components) 00703 { 00704 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00705 { 00706 SG_UNREF(m_components[i]); 00707 } 00708 00709 m_components=components; 00710 00711 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00712 { 00713 SG_REF(m_components[i]); 00714 } 00715 } 00716 00717 SGMatrix<float64_t> CGMM::alpha_init(SGMatrix<float64_t> init_means) 00718 { 00719 CDotFeatures* dotdata=(CDotFeatures *) features; 00720 int32_t num_vectors=dotdata->get_num_vectors(); 00721 00722 SGVector<float64_t> label_num(init_means.num_cols); 00723 00724 for (int32_t i=0; i<init_means.num_cols; i++) 00725 label_num.vector[i]=i; 00726 00727 CKNN* knn=new CKNN(1, new CEuclideanDistance(), new CMulticlassLabels(label_num)); 00728 knn->train(new CDenseFeatures<float64_t>(init_means)); 00729 CMulticlassLabels* init_labels=(CMulticlassLabels*) knn->apply(features); 00730 00731 SGMatrix<float64_t> alpha(num_vectors, int32_t(m_components.size())); 00732 memset(alpha.matrix, 0, num_vectors*m_components.size()*sizeof(float64_t)); 00733 00734 for (int32_t i=0; i<num_vectors; i++) 00735 alpha.matrix[i*m_components.size()+init_labels->get_int_label(i)]=1; 00736 00737 SG_UNREF(init_labels); 00738 00739 return alpha; 00740 } 00741 00742 SGVector<float64_t> CGMM::sample() 00743 { 00744 ASSERT(m_components.size()); 00745 float64_t rand_num=CMath::random(float64_t(0), float64_t(1)); 00746 float64_t cum_sum=0; 00747 for (int32_t i=0; i<m_coefficients.vlen; i++) 00748 { 00749 cum_sum+=m_coefficients.vector[i]; 00750 if (cum_sum>=rand_num) 00751 return m_components[i]->sample(); 00752 } 00753 return m_components[m_coefficients.vlen-1]->sample(); 00754 } 00755 00756 SGVector<float64_t> CGMM::cluster(SGVector<float64_t> point) 00757 { 00758 SGVector<float64_t> answer(m_components.size()+1); 00759 answer.vector[m_components.size()]=0; 00760 00761 for (int32_t i=0; i<int32_t(m_components.size()); i++) 00762 { 00763 answer.vector[i]=m_components[i]->compute_log_PDF(point)+CMath::log(m_coefficients[i]); 00764 answer.vector[m_components.size()]+=CMath::exp(answer.vector[i]); 00765 } 00766 answer.vector[m_components.size()]=CMath::log(answer.vector[m_components.size()]); 00767 00768 return answer; 00769 } 00770 00771 void CGMM::register_params() 00772 { 00773 m_parameters->add((SGVector<CSGObject*>*) &m_components, "m_components", "Mixture components"); 00774 m_parameters->add(&m_coefficients, "m_coefficients", "Mixture coefficients."); 00775 } 00776 00777 #endif