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

SHOGUN Machine Learning Toolbox - Documentation