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::distance_to_centroid(scalar_t *x, int k) {
39 for(int d = 0; d < _dim; d++) {
40 dist += sq(_cluster_means[k][d] - x[d]) / (2 * _cluster_var[k][d]);
41 dist += 0.5 * log(_cluster_var[k][d]);
42 ASSERT(!isnan(dist) && !isinf(dist));
47 void Clusterer::initialize_clusters(int nb_points, scalar_t **points) {
48 int *used = new int[nb_points];
49 for(int k = 0; k < nb_points; k++) { used[k] = 0; }
50 for(int k = 0; k < _nb_clusters; k++) {
52 do { l = int(drand48() * nb_points); } while(used[l]);
53 for(int d = 0; d < _dim; d++) {
54 _cluster_means[k][d] = points[l][d];
55 _cluster_var[k][d] = 1.0;
62 scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
63 int nb_classes, int *labels,
65 int *associated_clusters = new int[nb_points];
66 scalar_t total_dist = 0;
68 for(int n = 0; n < nb_points; n++) {
69 scalar_t lowest_dist = 0;
70 for(int k = 0; k < _nb_clusters; k++) {
71 scalar_t dist = distance_to_centroid(points[n], k);
72 if(k == 0 || dist <= lowest_dist) {
74 associated_clusters[n] = k;
78 total_dist += lowest_dist;
81 for(int n = 0; n < nb_points; n++) {
82 for(int k = 0; k < _nb_clusters; k++) {
85 gamma[n][associated_clusters[n]] = 1.0;
88 delete[] associated_clusters;
93 scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **points,
94 int nb_classes, int *labels,
98 int *ia = new int[nb_points * _nb_clusters + 1];
99 int *ja = new int[nb_points * _nb_clusters + 1];
100 scalar_t *ar = new scalar_t[nb_points * _nb_clusters + 1];
102 lp = glp_create_prob();
104 glp_set_prob_name(lp, "baseline_lp_cluster_association");
105 glp_set_obj_dir(lp, GLP_MIN);
107 glp_add_rows(lp, nb_points);
109 for(int n = 1; n <= nb_points; n++) {
110 glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
113 glp_add_cols(lp, nb_points * _nb_clusters);
114 for(int k = 1; k <= _nb_clusters; k++) {
115 for(int n = 1; n <= nb_points; n++) {
116 int i = n + nb_points * (k - 1);
117 glp_set_obj_coef(lp, i, distance_to_centroid(points[n-1], k-1));
118 glp_set_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
124 for(int n = 1; n <= nb_points; n++) {
125 for(int k = 1; k <= _nb_clusters; k++) {
127 ja[l] = n + nb_points * (k - 1);
133 glp_load_matrix(lp, nb_points * _nb_clusters, ia, ja, ar);
135 glp_simplex(lp, NULL);
137 scalar_t total_dist = glp_get_obj_val(lp);
139 for(int k = 1; k <= _nb_clusters; k++) {
140 for(int n = 1; n <= nb_points; n++) {
141 int i = n + nb_points * (k - 1);
142 gamma[n-1][k-1] = glp_get_col_prim(lp, i);
155 scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
156 int nb_classes, int *labels,
160 int nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
162 int *ia = new int[nb_coeffs + 1];
163 int *ja = new int[nb_coeffs + 1];
164 scalar_t *ar = new scalar_t[nb_coeffs + 1];
166 scalar_t *nb_samples_per_class = new scalar_t[nb_classes];
167 for(int c = 0; c < nb_classes; c++) {
168 nb_samples_per_class[c] = 0.0;
171 for(int n = 0; n < nb_points; n++) {
172 nb_samples_per_class[labels[n]] += 1.0;
175 lp = glp_create_prob();
177 glp_set_prob_name(lp, "uninformative_lp_cluster_association");
178 glp_set_obj_dir(lp, GLP_MIN);
180 // We have one constraint per point and one per cluster/class
182 glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
184 // (A) For each point, the constraint is that the sum of its
185 // association coefficients is equal to 1.0
187 for(int n = 1; n <= nb_points; n++) {
188 glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
191 // (B) For each pair cluster/class, the sum of the association
192 // coefficient to this cluster for this class is equal to the number
193 // of sample of that class, divided by the number of clusters
195 for(int k = 1; k <= _nb_clusters; k++) {
196 for(int c = 1; c <= nb_classes; c++) {
197 int row = nb_points + (k - 1) * nb_classes + c;
198 scalar_t tau = nb_samples_per_class[c-1] / scalar_t(_nb_clusters);
199 glp_set_row_bnds(lp, row, GLP_FX, tau, tau);
203 // Each one of the constraints above involve a linear combination of
204 // all the association coefficients
206 glp_add_cols(lp, nb_points * _nb_clusters);
208 for(int k = 1; k <= _nb_clusters; k++) {
209 for(int n = 1; n <= nb_points; n++) {
210 int r = n + nb_points * (k - 1);
212 // scalar_t dist = 0;
214 // for(int d = 0; d < _dim; d++) {
215 // dist += sq(_cluster_means[k-1][d] - points[n-1][d]) / (2 * _cluster_var[k-1][d]);
216 // dist += 0.5 * log(_cluster_var[k-1][d]);
219 // The LP weight on this association coefficient for the global
220 // loss is the normalized distance of that sample to the
221 // centroid of that cluster
223 glp_set_obj_coef(lp, r, distance_to_centroid(points[n-1], k-1));
225 // And this association coefficient is in [0,1]
227 glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
233 // We build the matrix of the LP coefficients
235 // The sums of the association coefficients per points for the
236 // constraints (A) above.
238 for(int n = 1; n <= nb_points; n++) {
239 for(int k = 1; k <= _nb_clusters; k++) {
242 ja[l] = nb_points * (k - 1) + n;
248 // And the sums of coefficients for each pair class/cluster for
249 // constraint (B) above.
251 for(int k = 1; k <= _nb_clusters; k++) {
252 for(int c = 1; c <= nb_classes; c++) {
253 int row = nb_points + (k - 1) * nb_classes + c;
254 for(int n = 1; n <= nb_points; n++) {
255 if(labels[n-1] == c - 1) {
257 ja[l] = (k-1) * nb_points + n;
265 ASSERT(l == nb_coeffs + 1);
267 glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
269 // Now a miracle occurs
271 glp_simplex(lp, NULL);
273 // We retrieve the result
275 scalar_t total_dist = glp_get_obj_val(lp);
277 for(int k = 1; k <= _nb_clusters; k++) {
278 for(int n = 1; n <= nb_points; n++) {
279 int i = n + nb_points * (k - 1);
280 gamma[n-1][k-1] = glp_get_col_prim(lp, i);
284 delete[] nb_samples_per_class;
293 void Clusterer::update_clusters(int nb_points, scalar_t **points, scalar_t **gamma) {
294 for(int k = 0; k < _nb_clusters; k++) {
296 for(int d = 0; d < _dim; d++) {
297 _cluster_means[k][d] = 0.0;
298 _cluster_var[k][d] = 0.0;
301 scalar_t sum_gamma = 0;
302 for(int n = 0; n < nb_points; n++) {
303 sum_gamma += gamma[n][k];
304 for(int d = 0; d < _dim; d++) {
305 _cluster_means[k][d] += gamma[n][k] * points[n][d];
306 _cluster_var[k][d] += gamma[n][k] * sq(points[n][d]);
310 ASSERT(sum_gamma >= 1);
312 for(int d = 0; d < _dim; d++) {
314 _cluster_var[k][d] = (_cluster_var[k][d] - sq(_cluster_means[k][d]) / sum_gamma) / (sum_gamma - 1);
316 _cluster_var[k][d] = 1;
318 _cluster_var[k][d] = max(0.01, _cluster_var[k][d]);
319 _cluster_means[k][d] /= sum_gamma;
324 void Clusterer::train(int mode,
325 int nb_clusters, int dim,
326 int nb_points, scalar_t **points,
327 int nb_classes, int *labels,
328 int *cluster_associations) {
329 deallocate_array<scalar_t>(_cluster_means);
330 deallocate_array<scalar_t>(_cluster_var);
331 _nb_clusters = nb_clusters;
333 _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
334 _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
336 scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
338 if(nb_clusters > nb_points) abort();
340 initialize_clusters(nb_points, points);
342 scalar_t pred_total_distance, total_distance = FLT_MAX;
346 pred_total_distance = total_distance;
350 case STANDARD_ASSOCIATION:
352 baseline_cluster_association(nb_points, points, nb_classes, labels, gammas);
355 case STANDARD_LP_ASSOCIATION:
357 baseline_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
360 case UNINFORMATIVE_LP_ASSOCIATION:
362 uninformative_lp_cluster_association(nb_points, points, nb_classes, labels, gammas);
366 cerr << "Unknown sample-cluster association mode." << endl;
370 cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
371 update_clusters(nb_points, points, gammas);
373 } while(total_distance < min_iteration_improvement * pred_total_distance &&
374 nb_rounds < max_nb_iterations);
376 if(cluster_associations) {
377 for(int n = 0; n < nb_points; n++) {
378 for(int k = 0; k < _nb_clusters; k++) {
379 if(k == 0 || gammas[n][k] > gammas[n][cluster_associations[n]]) {
380 cluster_associations[n] = k;
386 deallocate_array<scalar_t>(gammas);
389 int Clusterer::cluster(scalar_t *point) {
390 scalar_t lowest_dist = 0;
391 int associated_cluster = -1;
393 for(int k = 0; k < _nb_clusters; k++) {
396 for(int d = 0; d < _dim; d++) {
397 dist += sq(_cluster_means[k][d] - point[d]) / (2 * _cluster_var[k][d]);
398 dist += 0.5 * log(_cluster_var[k][d]);
399 ASSERT(!isnan(dist) && !isinf(dist));
402 if(k == 0 || dist <= lowest_dist) {
404 associated_cluster = k;
408 ASSERT(associated_cluster >= 0);
410 return associated_cluster;