Cosmetics.
[clueless-kmeans.git] / 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       glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
181     }
182   }
183
184   glp_add_cols(lp, nb_points * _nb_clusters);
185
186   for(int k = 1; k <= _nb_clusters; k++) {
187     for(int n = 1; n <= nb_points; n++) {
188       int r = n + nb_points * (k - 1);
189
190       scalar_t dist = 0;
191
192       for(int d = 0; d < _dim; d++) {
193         dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
194         dist += 0.5 * log(_cluster_var[k-1][d]);
195       }
196
197       glp_set_obj_coef(lp, r, dist);
198       glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
199     }
200   }
201
202   int l = 1;
203
204   for(int n = 1; n <= nb_points; n++) {
205     for(int k = 1; k <= _nb_clusters; k++) {
206       int row = n;
207       ia[l] = row;
208       ja[l] = nb_points * (k - 1) + n;
209       ar[l] = 1.0;
210       l++;
211     }
212   }
213
214   for(int k = 1; k <= _nb_clusters; k++) {
215     for(int c = 1; c <= nb_classes; c++) {
216       int row = nb_points + (k - 1) * nb_classes + c;
217       for(int n = 1; n <= nb_points; n++) {
218         if(labels[n-1] == c - 1) {
219           ia[l] = row;
220           ja[l] = (k-1)  * nb_points + n;
221           ar[l] = 1.0;
222           l++;
223         }
224       }
225     }
226   }
227
228   ASSERT(l == nb_coeffs + 1);
229
230   glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
231
232   glp_simplex(lp, NULL);
233
234   scalar_t total_dist = glp_get_obj_val(lp);
235
236   for(int k = 1; k <= _nb_clusters; k++) {
237     for(int n = 1; n <= nb_points; n++) {
238       int i = n + nb_points * (k - 1);
239       gamma[n-1][k-1] = glp_get_col_prim(lp, i);
240     }
241   }
242
243   delete[] nb_samples_per_class;
244   delete[] ia;
245   delete[] ja;
246   delete[] ar;
247   glp_delete_prob(lp);
248
249   return total_dist;
250 }
251
252 void Clusterer::update_clusters(int nb_points, scalar_t **points, scalar_t **gamma)  {
253   for(int k = 0; k < _nb_clusters; k++) {
254
255     for(int d = 0; d < _dim; d++) {
256       _cluster_means[k][d] = 0.0;
257       _cluster_var[k][d] = 0.0;
258     }
259
260     scalar_t sum_gamma = 0;
261     for(int n = 0; n < nb_points; n++) {
262       sum_gamma += gamma[n][k];
263       for(int d = 0; d < _dim; d++) {
264         _cluster_means[k][d] += gamma[n][k] * points[n][d];
265         _cluster_var[k][d] += gamma[n][k] * sq(points[n][d]);
266       }
267     }
268
269     ASSERT(sum_gamma >= 1);
270
271     for(int d = 0; d < _dim; d++) {
272       if(sum_gamma >= 2) {
273         _cluster_var[k][d] = (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
274       } else {
275         _cluster_var[k][d] = 1;
276       }
277       _cluster_var[k][d] = max(0.01, _cluster_var[k][d]);
278       _cluster_means[k][d] /= sum_gamma;
279     }
280   }
281 }
282
283 void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
284   int *used = new int[nb_points];
285   for(int k = 0; k < nb_points; k++) { used[k] = 0; }
286   for(int k = 0; k < _nb_clusters; k++) {
287     int l;
288     do { l = int(drand48() * nb_points); } while(used[l]);
289     for(int d = 0; d < _dim; d++) {
290       _cluster_means[k][d] = points[l][d];
291       _cluster_var[k][d] = 1.0;
292     }
293     used[l] = 1;
294   }
295   delete[] used;
296 }
297
298 void Clusterer::train(int mode,
299                       int nb_clusters, int dim,
300                       int nb_points, scalar_t **points,
301                       int nb_classes, int *labels,
302                       int *cluster_associations) {
303   deallocate_array<scalar_t>(_cluster_means);
304   deallocate_array<scalar_t>(_cluster_var);
305   _nb_clusters = nb_clusters;
306   _dim = dim;
307   _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
308   _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
309
310   scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
311
312   if(nb_clusters > nb_points) abort();
313
314   initialize_clusters(nb_points, points);
315
316   scalar_t pred_total_distance, total_distance = FLT_MAX;
317   int nb_rounds = 0;
318
319   do {
320     pred_total_distance = total_distance;
321
322     switch(mode) {
323
324     case STANDARD_ASSOCIATION:
325     total_distance =
326       baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
327       break;
328
329     case STANDARD_LP_ASSOCIATION:
330     total_distance =
331       baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
332       break;
333
334     case UNINFORMATIVE_LP_ASSOCIATION:
335     total_distance =
336       uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
337       break;
338
339     default:
340       cerr << "Unknown sample-cluster association mode." << endl;
341       abort();
342     }
343
344     cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
345     update_clusters(nb_points, points, gammas);
346     nb_rounds++;
347   } while(total_distance < min_iteration_improvement * pred_total_distance &&
348           nb_rounds < max_nb_iterations);
349
350   if(cluster_associations) {
351     for(int n = 0; n < nb_points; n++) {
352       for(int k = 0; k < _nb_clusters; k++) {
353         if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
354           cluster_associations[n] = k;
355         }
356       }
357     }
358   }
359
360   deallocate_array<scalar_t>(gammas);
361 }
362
363 int Clusterer::cluster(scalar_t *point) {
364   scalar_t lowest_dist = 0;
365   int associated_cluster = -1;
366
367   for(int k = 0; k < _nb_clusters; k++) {
368     scalar_t dist = 0;
369
370     for(int d = 0; d < _dim; d++) {
371       dist += sq(_cluster_means[k][d] - point[d]) / (2 * _cluster_var[k][d]);
372       dist += 0.5 * log(_cluster_var[k][d]);
373       ASSERT(!isnan(dist) && !isinf(dist));
374     }
375
376     if(k == 0 || dist <= lowest_dist) {
377       lowest_dist = dist;
378       associated_cluster = k;
379     }
380   }
381
382   ASSERT(associated_cluster >= 0);
383
384   return associated_cluster;
385 }