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