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/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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation