X-Git-Url: https://fleuret.org/cgi-bin/gitweb/gitweb.cgi?a=blobdiff_plain;f=clusterer.cc;h=04a9af4eb35f243effc4192897ce5e9a5373eedb;hb=04d2b44ba34a811e1fab0b90d38ebd06cd918c52;hp=42af81b7ced3cc5f299c59911ff63965ad1a4920;hpb=2455f83ba251602d5e04640067094f09f03aaa3d;p=clueless-kmeans.git

diff --git a/clusterer.cc b/clusterer.cc
index 42af81b..04a9af4 100644
--- a/clusterer.cc
+++ b/clusterer.cc
@@ -34,6 +34,33 @@ Clusterer::~Clusterer() {
   deallocate_array<scalar_t>(_cluster_var);
 }
 
+scalar_t Clusterer::distance_to_centroid(scalar_t *x, int k) {
+  // We take the variance into account + the normalization term. This
+  // is between k-mean and EM with a diagonal covariance
+  scalar_t dist = 0;
+  for(int d = 0; d < _dim; d++) {
+    dist += sq(_cluster_means[k][d] - x[d]) / (2 * _cluster_var[k][d]);
+    dist += 0.5 * log(_cluster_var[k][d]);
+    ASSERT(!isnan(dist) && !isinf(dist));
+  }
+  return dist;
+}
+
+void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
+  int *used = new int[nb_points];
+  for(int k = 0; k < nb_points; k++) { used[k] = 0; }
+  for(int k = 0; k < _nb_clusters; k++) {
+    int l;
+    do { l = int(drand48() * nb_points); } while(used[l]);
+    for(int d = 0; d < _dim; d++) {
+      _cluster_means[k][d] = points[l][d];
+      _cluster_var[k][d] = 1.0;
+    }
+    used[l] = 1;
+  }
+  delete[] used;
+}
+
 scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
                                                  int nb_classes, int *labels,
                                                  scalar_t **gamma)  {
@@ -43,14 +70,7 @@ scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **point
   for(int n = 0; n < nb_points; n++) {
     scalar_t lowest_dist = 0;
     for(int k = 0; k < _nb_clusters; k++) {
-      scalar_t dist = 0;
-
-      for(int d = 0; d < _dim; d++) {
-        dist += sq(_cluster_means[k][d] - points[n][d]) / (2 * _cluster_var[k][d]);
-        dist += 0.5 * log(_cluster_var[k][d]);
-        ASSERT(!isnan(dist) && !isinf(dist));
-      }
-
+      scalar_t dist = distance_to_centroid(points[n], k);
       if(k == 0 || dist <= lowest_dist) {
         lowest_dist = dist;
         associated_clusters[n] = k;
@@ -96,15 +116,7 @@ scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **po
   for(int k = 1; k <= _nb_clusters; k++) {
     for(int n = 1; n <= nb_points; n++) {
       int i = n + nb_points * (k - 1);
-
-      scalar_t dist = 0;
-
-      for(int d = 0; d < _dim; d++) {
-        dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
-        dist += 0.5 * log(_cluster_var[k-1][d]);
-      }
-
-      glp_set_obj_coef(lp, i, dist);
+      glp_set_obj_coef(lp, i, distance_to_centroid(points[n-1], k-1));
       glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
     }
   }
@@ -167,41 +179,58 @@ scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t
   glp_set_prob_name(lp, "uninformative_lp_cluster_association");
   glp_set_obj_dir(lp, GLP_MIN);
 
+  // We have one constraint per point and one per cluster/class
+
   glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
 
+  // (A) For each point, the constraint is that the sum of its
+  // association coefficients is equal to 1.0
+
   for(int n = 1; n <= nb_points; n++) {
-    glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
+    int row = n;
+    glp_set_row_bnds(lp, row, GLP_FX, 1.0, 1.0);
   }
 
+  // (B) For each pair cluster/class, the sum of the association
+  // coefficient to this cluster for this class is equal to the number
+  // of sample of that class, divided by the number of clusters
+
   for(int k = 1; k <= _nb_clusters; k++) {
     for(int c = 1; c <= nb_classes; c++) {
       int row = nb_points + (k - 1) * nb_classes + c;
       scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
-      // cout << "tau " << k << " " << c << " " << tau << endl;
       glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
     }
   }
 
+  // Each one of the constraints above involve a linear combination of
+  // all the association coefficients
+
   glp_add_cols(lp, nb_points * _nb_clusters);
 
   for(int k = 1; k <= _nb_clusters; k++) {
     for(int n = 1; n <= nb_points; n++) {
-      int r = n + nb_points * (k - 1);
+      int col = n + nb_points * (k - 1);
 
-      scalar_t dist = 0;
+      // The LP weight on this association coefficient for the global
+      // loss is the normalized distance of that sample to the
+      // centroid of that cluster
 
-      for(int d = 0; d < _dim; d++) {
-        dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
-        dist += 0.5 * log(_cluster_var[k-1][d]);
-      }
+      glp_set_obj_coef(lp, col, distance_to_centroid(points[n-1], k-1));
+
+      // And all the association coefficient is in [0,1]
 
-      glp_set_obj_coef(lp, r, dist);
-      glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
+      glp_set_col_bnds(lp, col, GLP_DB, 0.0, 1.0);
     }
   }
 
   int l = 1;
 
+  // We build the matrix of the LP coefficients
+
+  // The sums of the association coefficients per points for the
+  // constraints (A) above.
+
   for(int n = 1; n <= nb_points; n++) {
     for(int k = 1; k <= _nb_clusters; k++) {
       int row = n;
@@ -212,6 +241,9 @@ scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t
     }
   }
 
+  // And the sums of coefficients for each pair class/cluster for
+  // constraint (B) above.
+
   for(int k = 1; k <= _nb_clusters; k++) {
     for(int c = 1; c <= nb_classes; c++) {
       int row = nb_points + (k - 1) * nb_classes + c;
@@ -230,8 +262,12 @@ scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t
 
   glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
 
+  // Now a miracle occurs
+
   glp_simplex(lp, NULL);
 
+  // We retrieve the result
+
   scalar_t total_dist = glp_get_obj_val(lp);
 
   for(int k = 1; k <= _nb_clusters; k++) {
@@ -241,35 +277,6 @@ scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t
     }
   }
 
-  // { // ******************************* START ***************************
-// #warning Test code added on 2013 Feb 07 20:32:05
-    // // for(int n = 0; n < nb_points; n++) {
-    // // scalar_t sum = 0;
-    // // for(int k = 0; k < _nb_clusters; k++) {
-    // // ASSERT(gamma[n][k] >= 0 && gamma[n][k] <= 1);
-    // // sum += gamma[n][k];
-    // // }
-    // // cout << sum << endl;
-    // // }
-
-    // scalar_t *sum_gamma = new scalar_t[nb_classes];
-
-    // for(int k = 0; k < _nb_clusters; k++) {
-      // for(int c = 0; c < nb_classes; c++) { sum_gamma[c] = 0.0; }
-      // for(int n = 0; n < nb_points; n++) {
-        // sum_gamma[labels[n]] += gamma[n][k];
-      // }
-      // cout << "CLUSTER" << k;
-      // for(int c = 0; c < nb_classes; c++) {
-        // cout << " " << sum_gamma[c];
-      // }
-      // cout << endl;
-    // }
-
-    // delete sum_gamma;
-
-  // } // ******************************** END ****************************
-
   delete[] nb_samples_per_class;
   delete[] ia;
   delete[] ja;
@@ -279,7 +286,7 @@ scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t
   return total_dist;
 }
 
-void Clusterer::baseline_update_clusters(int nb_points, scalar_t **points, scalar_t **gamma)  {
+void Clusterer::update_clusters(int nb_points, scalar_t **points, scalar_t **gamma)  {
   for(int k = 0; k < _nb_clusters; k++) {
 
     for(int d = 0; d < _dim; d++) {
@@ -300,32 +307,20 @@ void Clusterer::baseline_update_clusters(int nb_points, scalar_t **points, scala
 
     for(int d = 0; d < _dim; d++) {
       if(sum_gamma >= 2) {
-        _cluster_var[k][d] = (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
+        _cluster_var[k][d] =
+          (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
+        _cluster_var[k][d] = max(scalar_t(min_cluster_variance), _cluster_var[k][d]);
       } else {
         _cluster_var[k][d] = 1;
       }
-      _cluster_var[k][d] = max(0.01, _cluster_var[k][d]);
-      _cluster_means[k][d] /= sum_gamma;
-    }
-  }
-}
 
-void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
-  int *used = new int[nb_points];
-  for(int k = 0; k < nb_points; k++) { used[k] = 0; }
-  for(int k = 0; k < _nb_clusters; k++) {
-    int l;
-    do { l = int(drand48() * nb_points); } while(used[l]);
-    for(int d = 0; d < _dim; d++) {
-      _cluster_means[k][d] = points[l][d];
-      _cluster_var[k][d] = 1.0;
+      _cluster_means[k][d] /= sum_gamma;
     }
-    used[l] = 1;
   }
-  delete[] used;
 }
 
-void Clusterer::train(int nb_clusters, int dim,
+void Clusterer::train(int mode,
+                      int nb_clusters, int dim,
                       int nb_points, scalar_t **points,
                       int nb_classes, int *labels,
                       int *cluster_associations) {
@@ -347,48 +342,41 @@ void Clusterer::train(int nb_clusters, int dim,
 
   do {
     pred_total_distance = total_distance;
+
+    switch(mode) {
+
+    case STANDARD_ASSOCIATION:
+    total_distance =
+      baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
+      break;
+
+    case STANDARD_LP_ASSOCIATION:
+    total_distance =
+      baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
+      break;
+
+    case UNINFORMATIVE_LP_ASSOCIATION:
     total_distance =
-      // baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
-      // baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
       uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
+      break;
+
+    default:
+      cerr << "Unknown sample-cluster association mode." << endl;
+      abort();
+    }
+
     cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
-    baseline_update_clusters(nb_points, points, gammas);
+    update_clusters(nb_points, points, gammas);
     nb_rounds++;
   } while(total_distance < min_iteration_improvement * pred_total_distance &&
           nb_rounds < max_nb_iterations);
 
-  {
-    cout << "TOTAL_NB_SAMPLES";
-    for(int c = 0; c < nb_classes; c++) {
-      int nb_samples = 0;
-      for(int n = 0; n < nb_points; n++) {
-        if(labels[n] == c) {
-          nb_samples++;
-        }
-      }
-      cout << " " << nb_samples;
-    }
-    cout << endl;
-
-    for(int k = 0; k < _nb_clusters; k++) {
-      cout << "CLUSTER_GAMMA_SUM " << k << " :";
-      for(int c = 0; c < nb_classes; c++) {
-        scalar_t sum = 0.0;
-        for(int n = 0; n < nb_points; n++) {
-          if(labels[n] == c) {
-            sum += gammas[n][k];
-          }
+  if(cluster_associations) {
+    for(int n = 0; n < nb_points; n++) {
+      for(int k = 0; k < _nb_clusters; k++) {
+        if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
+          cluster_associations[n] = k;
         }
-        cout << " " << sum;
-      }
-      cout << endl;
-    }
-  }
-
-  for(int n = 0; n < nb_points; n++) {
-    for(int k = 0; k < _nb_clusters; k++) {
-      if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
-        cluster_associations[n] = k;
       }
     }
   }