2 * clueless-kmean is a variant of k-mean which enforces balanced
3 * distribution of classes in every cluster
5 * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
6 * Written by Francois Fleuret <francois.fleuret@idiap.ch>
8 * This file is part of clueless-kmean.
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.
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.
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/>.
24 #include "clusterer.h"
27 Clusterer::Clusterer() {
32 Clusterer::~Clusterer() {
33 deallocate_array<scalar_t>(_cluster_means);
34 deallocate_array<scalar_t>(_cluster_var);
37 scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
38 int nb_classes, int *labels,
40 int *associated_clusters = new int[nb_points];
41 scalar_t total_dist = 0;
43 for(int n = 0; n < nb_points; n++) {
44 scalar_t lowest_dist = 0;
45 for(int k = 0; k < _nb_clusters; k++) {
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));
54 if(k == 0 || dist <= lowest_dist) {
56 associated_clusters[n] = k;
60 total_dist += lowest_dist;
63 for(int n = 0; n < nb_points; n++) {
64 for(int k = 0; k < _nb_clusters; k++) {
67 gamma[n][associated_clusters[n]] = 1.0;
70 delete[] associated_clusters;
75 scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **points,
76 int nb_classes, int *labels,
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];
84 lp = glp_create_prob();
86 glp_set_prob_name(lp, "baseline_lp_cluster_association");
87 glp_set_obj_dir(lp, GLP_MIN);
89 glp_add_rows(lp, nb_points);
91 for(int n = 1; n <= nb_points; n++) {
92 glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
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);
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]);
107 glp_set_obj_coef(lp, i, dist);
108 glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
114 for(int n = 1; n <= nb_points; n++) {
115 for(int k = 1; k <= _nb_clusters; k++) {
117 ja[l] = n + nb_points * (k - 1);
123 glp_load_matrix(lp, nb_points * _nb_clusters, ia, ja, ar);
125 glp_simplex(lp, NULL);
127 scalar_t total_dist = glp_get_obj_val(lp);
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);
145 scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
146 int nb_classes, int *labels,
150 int nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
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];
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;
161 for(int n = 0; n < nb_points; n++) {
162 nb_samples_per_class[labels[n]] += 1.0;
165 lp = glp_create_prob();
167 glp_set_prob_name(lp, "uninformative_lp_cluster_association");
168 glp_set_obj_dir(lp, GLP_MIN);
170 glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
172 for(int n = 1; n <= nb_points; n++) {
173 glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
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);
185 glp_add_cols(lp, nb_points * _nb_clusters);
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);
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]);
198 glp_set_obj_coef(lp, r, dist);
199 glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
205 for(int n = 1; n <= nb_points; n++) {
206 for(int k = 1; k <= _nb_clusters; k++) {
209 ja[l] = nb_points * (k - 1) + n;
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) {
221 ja[l] = (k-1) * nb_points + n;
229 ASSERT(l == nb_coeffs + 1);
231 glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
233 glp_simplex(lp, NULL);
235 scalar_t total_dist = glp_get_obj_val(lp);
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);
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];
252 // // cout << sum << endl;
255 // scalar_t *sum_gamma = new scalar_t[nb_classes];
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];
262 // cout << "CLUSTER" << k;
263 // for(int c = 0; c < nb_classes; c++) {
264 // cout << " " << sum_gamma[c];
271 // } // ******************************** END ****************************
273 delete[] nb_samples_per_class;
282 void Clusterer::baseline_update_clusters(int nb_points, scalar_t **points, scalar_t **gamma) {
283 for(int k = 0; k < _nb_clusters; k++) {
285 for(int d = 0; d < _dim; d++) {
286 _cluster_means[k][d] = 0.0;
287 _cluster_var[k][d] = 0.0;
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]);
299 ASSERT(sum_gamma >= 1);
301 for(int d = 0; d < _dim; d++) {
303 _cluster_var[k][d] = (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
305 _cluster_var[k][d] = 1;
307 _cluster_var[k][d] = max(0.01, _cluster_var[k][d]);
308 _cluster_means[k][d] /= sum_gamma;
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++) {
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;
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;
336 _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
337 _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
339 scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
341 if(nb_clusters > nb_points) abort();
343 initialize_clusters(nb_points, points);
345 scalar_t pred_total_distance, total_distance = FLT_MAX;
349 pred_total_distance = 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);
357 } while(total_distance < min_iteration_improvement * pred_total_distance &&
358 nb_rounds < max_nb_iterations);
361 cout << "TOTAL_NB_SAMPLES";
362 for(int c = 0; c < nb_classes; c++) {
364 for(int n = 0; n < nb_points; n++) {
369 cout << " " << nb_samples;
373 for(int k = 0; k < _nb_clusters; k++) {
374 cout << "CLUSTER_GAMMA_SUM " << k << " :";
375 for(int c = 0; c < nb_classes; c++) {
377 for(int n = 0; n < nb_points; n++) {
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;
396 deallocate_array<scalar_t>(gammas);
399 int Clusterer::cluster(scalar_t *point) {
400 scalar_t lowest_dist = 0;
401 int associated_cluster = -1;
403 for(int k = 0; k < _nb_clusters; k++) {
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));
412 if(k == 0 || dist <= lowest_dist) {
414 associated_cluster = k;
418 ASSERT(associated_cluster >= 0);
420 return associated_cluster;