clusterer.cc

   1 /*
   2  *  clueless-kmean is a variant of k-mean which enforces balanced
   3  *  distribution of classes in every cluster
   4  *
   5  *  Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
   6  *  Written by Francois Fleuret <francois.fleuret@idiap.ch>
   7  *
   8  *  This file is part of clueless-kmean.
   9  *
  10  *  clueless-kmean is free software: you can redistribute it and/or
  11  *  modify it under the terms of the GNU General Public License
  12  *  version 3 as published by the Free Software Foundation.
  13  *
  14  *  clueless-kmean is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  *  General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with selector.  If not, see <http://www.gnu.org/licenses/>.
  21  *
  22  */
  23
  24 #include "clusterer.h"
  25 #include <glpk.h>
  26
  27 Clusterer::Clusterer() {
  28   _cluster_means = 0;
  29   _cluster_var = 0;
  30 }
  31
  32 Clusterer::~Clusterer() {
  33   deallocate_array<scalar_t>(_cluster_means);
  34   deallocate_array<scalar_t>(_cluster_var);
  35 }
  36
  37 scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
  38                                                  int nb_classes, int *labels,
  39                                                  scalar_t **gamma)  {
  40   int *associated_clusters = new int[nb_points];
  41   scalar_t total_dist = 0;
  42
  43   for(int n = 0; n < nb_points; n++) {
  44     scalar_t lowest_dist = 0;
  45     for(int k = 0; k < _nb_clusters; k++) {
  46       scalar_t dist = 0;
  47
  48       for(int d = 0; d < _dim; d++) {
  49         dist += sq(_cluster_means[k][d] - points[n][d]) / (2 * _cluster_var[k][d]);
  50         dist += 0.5 * log(_cluster_var[k][d]);
  51         ASSERT(!isnan(dist) && !isinf(dist));
  52       }
  53
  54       if(k == 0 || dist <= lowest_dist) {
  55         lowest_dist = dist;
  56         associated_clusters[n] = k;
  57       }
  58     }
  59
  60     total_dist += lowest_dist;
  61   }
  62
  63   for(int n = 0; n < nb_points; n++) {
  64     for(int k = 0; k < _nb_clusters; k++) {
  65       gamma[n][k] = 0.0;
  66     }
  67     gamma[n][associated_clusters[n]] = 1.0;
  68   }
  69
  70   delete[] associated_clusters;
  71
  72   return total_dist;
  73 }
  74
  75 scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **points,
  76                                                     int nb_classes, int *labels,
  77                                                     scalar_t **gamma)  {
  78   glp_prob *lp;
  79
  80   int *ia = new int[nb_points * _nb_clusters + 1];
  81   int *ja = new int[nb_points * _nb_clusters + 1];
  82   scalar_t *ar = new scalar_t[nb_points * _nb_clusters + 1];
  83
  84   lp = glp_create_prob();
  85
  86   glp_set_prob_name(lp, "baseline_lp_cluster_association");
  87   glp_set_obj_dir(lp, GLP_MIN);
  88
  89   glp_add_rows(lp, nb_points);
  90
  91   for(int n = 1; n <= nb_points; n++) {
  92     glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
  93   }
  94
  95   glp_add_cols(lp, nb_points * _nb_clusters);
  96   for(int k = 1; k <= _nb_clusters; k++) {
  97     for(int n = 1; n <= nb_points; n++) {
  98       int i = n + nb_points * (k - 1);
  99
 100       scalar_t dist = 0;
 101
 102       for(int d = 0; d < _dim; d++) {
 103         dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
 104         dist += 0.5 * log(_cluster_var[k-1][d]);
 105       }
 106
 107       glp_set_obj_coef(lp, i, dist);
 108       glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
 109     }
 110   }
 111
 112   int l = 1;
 113
 114   for(int n = 1; n <= nb_points; n++) {
 115     for(int k = 1; k <= _nb_clusters; k++) {
 116       ia[l] = n;
 117       ja[l] = n + nb_points * (k - 1);
 118       ar[l] = 1.0;
 119       l++;
 120     }
 121   }
 122
 123   glp_load_matrix(lp, nb_points * _nb_clusters, ia, ja, ar);
 124
 125   glp_simplex(lp, NULL);
 126
 127   scalar_t total_dist = glp_get_obj_val(lp);
 128
 129   for(int k = 1; k <= _nb_clusters; k++) {
 130     for(int n = 1; n <= nb_points; n++) {
 131       int i = n + nb_points * (k - 1);
 132       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
 133     }
 134   }
 135
 136   delete[] ia;
 137   delete[] ja;
 138   delete[] ar;
 139
 140   glp_delete_prob(lp);
 141
 142   return total_dist;
 143 }
 144
 145 scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
 146                                                          int nb_classes, int *labels,
 147                                                          scalar_t **gamma)  {
 148   glp_prob *lp;
 149
 150   int nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
 151
 152   int *ia = new int[nb_coeffs + 1];
 153   int *ja = new int[nb_coeffs + 1];
 154   scalar_t *ar = new scalar_t[nb_coeffs + 1];
 155
 156   scalar_t *nb_samples_per_class = new scalar_t[nb_classes];
 157   for(int c = 0; c < nb_classes; c++) {
 158     nb_samples_per_class[c] = 0.0;
 159   }
 160
 161   for(int n = 0; n < nb_points; n++) {
 162     nb_samples_per_class[labels[n]] += 1.0;
 163   }
 164
 165   lp = glp_create_prob();
 166
 167   glp_set_prob_name(lp, "uninformative_lp_cluster_association");
 168   glp_set_obj_dir(lp, GLP_MIN);
 169
 170   glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
 171
 172   for(int n = 1; n <= nb_points; n++) {
 173     glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
 174   }
 175
 176   for(int k = 1; k <= _nb_clusters; k++) {
 177     for(int c = 1; c <= nb_classes; c++) {
 178       int row = nb_points + (k - 1) * nb_classes + c;
 179       scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
 180       // cout << "tau " << k << " " << c << " " << tau << endl;
 181       glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
 182     }
 183   }
 184
 185   glp_add_cols(lp, nb_points * _nb_clusters);
 186
 187   for(int k = 1; k <= _nb_clusters; k++) {
 188     for(int n = 1; n <= nb_points; n++) {
 189       int r = n + nb_points * (k - 1);
 190
 191       scalar_t dist = 0;
 192
 193       for(int d = 0; d < _dim; d++) {
 194         dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
 195         dist += 0.5 * log(_cluster_var[k-1][d]);
 196       }
 197
 198       glp_set_obj_coef(lp, r, dist);
 199       glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
 200     }
 201   }
 202
 203   int l = 1;
 204
 205   for(int n = 1; n <= nb_points; n++) {
 206     for(int k = 1; k <= _nb_clusters; k++) {
 207       int row = n;
 208       ia[l] = row;
 209       ja[l] = nb_points * (k - 1) + n;
 210       ar[l] = 1.0;
 211       l++;
 212     }
 213   }
 214
 215   for(int k = 1; k <= _nb_clusters; k++) {
 216     for(int c = 1; c <= nb_classes; c++) {
 217       int row = nb_points + (k - 1) * nb_classes + c;
 218       for(int n = 1; n <= nb_points; n++) {
 219         if(labels[n-1] == c - 1) {
 220           ia[l] = row;
 221           ja[l] = (k-1)  * nb_points + n;
 222           ar[l] = 1.0;
 223           l++;
 224         }
 225       }
 226     }
 227   }
 228
 229   ASSERT(l == nb_coeffs + 1);
 230
 231   glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
 232
 233   glp_simplex(lp, NULL);
 234
 235   scalar_t total_dist = glp_get_obj_val(lp);
 236
 237   for(int k = 1; k <= _nb_clusters; k++) {
 238     for(int n = 1; n <= nb_points; n++) {
 239       int i = n + nb_points * (k - 1);
 240       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
 241     }
 242   }
 243
 244   // { // ******************************* START ***************************
 245 // #warning Test code added on 2013 Feb 07 20:32:05
 246     // // for(int n = 0; n < nb_points; n++) {
 247     // // scalar_t sum = 0;
 248     // // for(int k = 0; k < _nb_clusters; k++) {
 249     // // ASSERT(gamma[n][k] >= 0 && gamma[n][k] <= 1);
 250     // // sum += gamma[n][k];
 251     // // }
 252     // // cout << sum << endl;
 253     // // }
 254
 255     // scalar_t *sum_gamma = new scalar_t[nb_classes];
 256
 257     // for(int k = 0; k < _nb_clusters; k++) {
 258       // for(int c = 0; c < nb_classes; c++) { sum_gamma[c] = 0.0; }
 259       // for(int n = 0; n < nb_points; n++) {
 260         // sum_gamma[labels[n]] += gamma[n][k];
 261       // }
 262       // cout << "CLUSTER" << k;
 263       // for(int c = 0; c < nb_classes; c++) {
 264         // cout << " " << sum_gamma[c];
 265       // }
 266       // cout << endl;
 267     // }
 268
 269     // delete sum_gamma;
 270
 271   // } // ******************************** END ****************************
 272
 273   delete[] nb_samples_per_class;
 274   delete[] ia;
 275   delete[] ja;
 276   delete[] ar;
 277   glp_delete_prob(lp);
 278
 279   return total_dist;
 280 }
 281
 282 void Clusterer::baseline_update_clusters(int nb_points, scalar_t **points, scalar_t **gamma)  {
 283   for(int k = 0; k < _nb_clusters; k++) {
 284
 285     for(int d = 0; d < _dim; d++) {
 286       _cluster_means[k][d] = 0.0;
 287       _cluster_var[k][d] = 0.0;
 288     }
 289
 290     scalar_t sum_gamma = 0;
 291     for(int n = 0; n < nb_points; n++) {
 292       sum_gamma += gamma[n][k];
 293       for(int d = 0; d < _dim; d++) {
 294         _cluster_means[k][d] += gamma[n][k] * points[n][d];
 295         _cluster_var[k][d] += gamma[n][k] * sq(points[n][d]);
 296       }
 297     }
 298
 299     ASSERT(sum_gamma >= 1);
 300
 301     for(int d = 0; d < _dim; d++) {
 302       if(sum_gamma >= 2) {
 303         _cluster_var[k][d] = (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
 304       } else {
 305         _cluster_var[k][d] = 1;
 306       }
 307       _cluster_var[k][d] = max(0.01, _cluster_var[k][d]);
 308       _cluster_means[k][d] /= sum_gamma;
 309     }
 310   }
 311 }
 312
 313 void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
 314   int *used = new int[nb_points];
 315   for(int k = 0; k < nb_points; k++) { used[k] = 0; }
 316   for(int k = 0; k < _nb_clusters; k++) {
 317     int l;
 318     do { l = int(drand48() * nb_points); } while(used[l]);
 319     for(int d = 0; d < _dim; d++) {
 320       _cluster_means[k][d] = points[l][d];
 321       _cluster_var[k][d] = 1.0;
 322     }
 323     used[l] = 1;
 324   }
 325   delete[] used;
 326 }
 327
 328 void Clusterer::train(int nb_clusters, int dim,
 329                       int nb_points, scalar_t **points,
 330                       int nb_classes, int *labels,
 331                       int *cluster_associations) {
 332   deallocate_array<scalar_t>(_cluster_means);
 333   deallocate_array<scalar_t>(_cluster_var);
 334   _nb_clusters = nb_clusters;
 335   _dim = dim;
 336   _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
 337   _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
 338
 339   scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
 340
 341   if(nb_clusters > nb_points) abort();
 342
 343   initialize_clusters(nb_points, points);
 344
 345   scalar_t pred_total_distance, total_distance = FLT_MAX;
 346   int nb_rounds = 0;
 347
 348   do {
 349     pred_total_distance = total_distance;
 350     total_distance =
 351       // baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
 352       // baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
 353       uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
 354     cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
 355     baseline_update_clusters(nb_points, points, gammas);
 356     nb_rounds++;
 357   } while(total_distance < min_iteration_improvement * pred_total_distance &&
 358           nb_rounds < max_nb_iterations);
 359
 360   {
 361     cout << "TOTAL_NB_SAMPLES";
 362     for(int c = 0; c < nb_classes; c++) {
 363       int nb_samples = 0;
 364       for(int n = 0; n < nb_points; n++) {
 365         if(labels[n] == c) {
 366           nb_samples++;
 367         }
 368       }
 369       cout << " " << nb_samples;
 370     }
 371     cout << endl;
 372
 373     for(int k = 0; k < _nb_clusters; k++) {
 374       cout << "CLUSTER_GAMMA_SUM " << k << " :";
 375       for(int c = 0; c < nb_classes; c++) {
 376         scalar_t sum = 0.0;
 377         for(int n = 0; n < nb_points; n++) {
 378           if(labels[n] == c) {
 379             sum += gammas[n][k];
 380           }
 381         }
 382         cout << " " << sum;
 383       }
 384       cout << endl;
 385     }
 386   }
 387
 388   for(int n = 0; n < nb_points; n++) {
 389     for(int k = 0; k < _nb_clusters; k++) {
 390       if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
 391         cluster_associations[n] = k;
 392       }
 393     }
 394   }
 395
 396   deallocate_array<scalar_t>(gammas);
 397 }
 398
 399 int Clusterer::cluster(scalar_t *point) {
 400   scalar_t lowest_dist = 0;
 401   int associated_cluster = -1;
 402
 403   for(int k = 0; k < _nb_clusters; k++) {
 404     scalar_t dist = 0;
 405
 406     for(int d = 0; d < _dim; d++) {
 407       dist += sq(_cluster_means[k][d] - point[d]) / (2 * _cluster_var[k][d]);
 408       dist += 0.5 * log(_cluster_var[k][d]);
 409       ASSERT(!isnan(dist) && !isinf(dist));
 410     }
 411
 412     if(k == 0 || dist <= lowest_dist) {
 413       lowest_dist = dist;
 414       associated_cluster = k;
 415     }
 416   }
 417
 418   ASSERT(associated_cluster >= 0);
 419
 420   return associated_cluster;
 421 }