00001
00002
00003
00004
00005
00006
00007
00008
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
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
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
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
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