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/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
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
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
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
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
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