--- /dev/null
+
+# clueless-kmean is a variant of k-mean which enforces balanced
+# distribution of classes in every cluster
+#
+# Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+# Written by Francois Fleuret <francois.fleuret@idiap.ch>
+#
+# This file is part of clueless-kmean.
+#
+# clueless-kmean is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# version 3 as published by the Free Software Foundation.
+#
+# clueless-kmean is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with selector. If not, see <http://www.gnu.org/licenses/>.
+
+ifeq ($(STATIC),yes)
+ LDFLAGS=-static -lm -ljpeg -lpng -lz -lglpk
+else
+ LDFLAGS= -lm -ljpeg -lpng -lz -lglpk
+endif
+
+ifeq ($(DEBUG),yes)
+ OPTIMIZE_FLAG = -ggdb3 -DDEBUG -fno-omit-frame-pointer
+else
+ OPTIMIZE_FLAG = -ggdb3 -O3
+endif
+
+ifeq ($(PROFILE),yes)
+ PROFILE_FLAG = -pg
+endif
+
+CXXFLAGS = -Wall $(OPTIMIZE_FLAG) $(PROFILE_FLAG) $(CXXGLPK)
+
+# LDFLAGS=-lglpk
+
+all: clueless-kmean
+
+clueless-kmean: \
+ misc.o \
+ sample_set.o \
+ clusterer.o \
+ clueless-kmean.o
+ $(CXX) $(CXXFLAGS) -o $@ $^ $(LDFLAGS)
+
+clean:
+ rm -f *.o clueless-kmean
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef ARRAYS_H
+#define ARRAYS_H
+
+template<class T>
+T inline sqdist(int dim, T *x, T *y) {
+ T result = 0;
+ for(int k = 0; k < dim; k++) result += sq(x[k] - y[k]);
+ return result;
+}
+
+template<class T>
+T inline sum(int dim, T *x) {
+ T result = 0;
+ for(int k = 0; k < dim; k++) result += x[k];
+ return result;
+}
+
+template<class T>
+T inline dotprod(int dim, T *x, T *y) {
+ T result = 0;
+ for(int k = 0; k < dim; k++) result += x[k] * y[k];
+ return result;
+}
+
+template<class T>
+void fill_vector(int dim, T *vector, T v) {
+ for(int k = 0; k < dim; k++) {
+ vector[k] = v;
+ }
+}
+
+template<class T>
+void copy_vector(int dim, T *vector_dst, T *vector_src) {
+ for(int k = 0; k < dim; k++) {
+ vector_dst[k] = vector_src[k];
+ }
+}
+
+template<class T>
+T **allocate_array(int a, int b) {
+ T *whole = new T[a * b];
+ T **array = new T *[a];
+ for(int k = 0; k < a; k++) {
+ array[k] = whole;
+ whole += b;
+ }
+ return array;
+}
+
+template<class T>
+void deallocate_array(T **array) {
+ if(array) {
+ delete[] array[0];
+ delete[] array;
+ }
+}
+
+template<class T>
+void fill_array(int a, int b, T **array, T v) {
+ for(int k = 0; k < a * b; k++) {
+ array[0][k] = v;
+ }
+}
+
+template<class T>
+T ***allocate_volume(int a, int b, int c) {
+ T *whole = new T[a * b * c];
+ T **column = new T *[a * b];
+ T ***volume = new T **[a];
+
+ for(int k = 0; k < a; k++) {
+ volume[k] = column;
+ for(int l = 0; l < b; l++) {
+ column[l] = whole;
+ whole += c;
+ }
+ column += b;
+ }
+
+ return volume;
+}
+
+template<class T>
+void deallocate_volume(T ***volume) {
+ if(volume) {
+ delete[] volume[0][0];
+ delete[] volume[0];
+ delete[] volume;
+ }
+}
+
+template<class T>
+void fill_volume(int a, int b, int c, T ***volume, T v) {
+ for(int k = 0; k < a * b * c; k++) {
+ volume[0][0][k] = v;
+ }
+}
+
+#endif
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <iostream>
+#include <fstream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+#include <glpk.h>
+
+using namespace std;
+
+#include "misc.h"
+#include "arrays.h"
+#include "sample_set.h"
+#include "clusterer.h"
+
+void generate_toy_problem(SampleSet *sample_set) {
+ int dim = 2;
+ int nb_points = 1000;
+
+ sample_set->resize(dim, nb_points);
+ sample_set->nb_classes = 2;
+
+ for(int n = 0; n < nb_points; n++) {
+ sample_set->labels[n] = int(drand48() * 2);
+ if(sample_set->labels[n] == 0) {
+ sample_set->points[n][0] = (2 * drand48() - 1) * 0.8;
+ sample_set->points[n][1] = - 0.6 + (2 * drand48() - 1) * 0.4;
+ } else {
+ sample_set->points[n][0] = (2 * drand48() - 1) * 0.4;
+ sample_set->points[n][1] = 0.6 + (2 * drand48() - 1) * 0.4;
+ }
+ }
+}
+
+int main(int argc, char **argv) {
+ SampleSet sample_set;
+ Clusterer clusterer;
+ int nb_clusters = 3;
+
+ generate_toy_problem(&sample_set);
+
+ {
+ ofstream out("points.dat");
+ for(int n = 0; n < sample_set.nb_points; n++) {
+ out << sample_set.labels[n];
+ for(int d = 0; d < sample_set.dim; d++) {
+ out << " " << sample_set.points[n][d];
+ }
+ out << endl;
+ }
+ }
+
+ int *associated_clusters = new int[sample_set.nb_points];
+
+ glp_term_out(0);
+
+ clusterer.train(nb_clusters, sample_set.dim,
+ sample_set.nb_points, sample_set.points,
+ sample_set.nb_classes, sample_set.labels,
+ associated_clusters);
+
+ {
+ ofstream out("associated_clusters.dat");
+ for(int n = 0; n < sample_set.nb_points; n++) {
+ out << associated_clusters[n];
+ for(int d = 0; d < sample_set.dim; d++) {
+ out << " " << sample_set.points[n][d];
+ }
+ out << endl;
+ }
+ }
+
+ {
+ ofstream out("clusters.dat");
+ for(int k = 0 ; k < clusterer._nb_clusters; k++) {
+ out << k;
+ for(int d = 0; d < sample_set.dim; d++) {
+ out << " " << clusterer._cluster_means[k][d];
+ }
+ for(int d = 0; d < sample_set.dim; d++) {
+ out << " " << 2 * sqrt(clusterer._cluster_var[k][d]);
+ }
+ out << endl;
+ }
+ }
+
+ delete[] associated_clusters;
+}
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "clusterer.h"
+#include <glpk.h>
+
+Clusterer::Clusterer() {
+ _cluster_means = 0;
+ _cluster_var = 0;
+}
+
+Clusterer::~Clusterer() {
+ deallocate_array<scalar_t>(_cluster_means);
+ deallocate_array<scalar_t>(_cluster_var);
+}
+
+scalar_t Clusterer::baseline_cluster_association(int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ scalar_t **gamma) {
+ int *associated_clusters = new int[nb_points];
+ scalar_t total_dist = 0;
+
+ 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));
+ }
+
+ if(k == 0 || dist <= lowest_dist) {
+ lowest_dist = dist;
+ associated_clusters[n] = k;
+ }
+ }
+
+ total_dist += lowest_dist;
+ }
+
+ for(int n = 0; n < nb_points; n++) {
+ for(int k = 0; k < _nb_clusters; k++) {
+ gamma[n][k] = 0.0;
+ }
+ gamma[n][associated_clusters[n]] = 1.0;
+ }
+
+ delete[] associated_clusters;
+
+ return total_dist;
+}
+
+scalar_t Clusterer::baseline_lp_cluster_association(int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ scalar_t **gamma) {
+ glp_prob *lp;
+
+ int *ia = new int[nb_points * _nb_clusters + 1];
+ int *ja = new int[nb_points * _nb_clusters + 1];
+ scalar_t *ar = new scalar_t[nb_points * _nb_clusters + 1];
+
+ lp = glp_create_prob();
+
+ glp_set_prob_name(lp, "baseline_lp_cluster_association");
+ glp_set_obj_dir(lp, GLP_MIN);
+
+ glp_add_rows(lp, nb_points);
+
+ for(int n = 1; n <= nb_points; n++) {
+ glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
+ }
+
+ 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 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_col_bnds(lp, i, GLP_DB, 0.0, 1.0);
+ }
+ }
+
+ int l = 1;
+
+ for(int n = 1; n <= nb_points; n++) {
+ for(int k = 1; k <= _nb_clusters; k++) {
+ ia[l] = n;
+ ja[l] = n + nb_points * (k - 1);
+ ar[l] = 1.0;
+ l++;
+ }
+ }
+
+ glp_load_matrix(lp, nb_points * _nb_clusters, ia, ja, ar);
+
+ glp_simplex(lp, NULL);
+
+ scalar_t total_dist = glp_get_obj_val(lp);
+
+ for(int k = 1; k <= _nb_clusters; k++) {
+ for(int n = 1; n <= nb_points; n++) {
+ int i = n + nb_points * (k - 1);
+ gamma[n-1][k-1] = glp_get_col_prim(lp, i);
+ }
+ }
+
+ delete[] ia;
+ delete[] ja;
+ delete[] ar;
+
+ glp_delete_prob(lp);
+
+ return total_dist;
+}
+
+scalar_t Clusterer::uninformative_lp_cluster_association(int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ scalar_t **gamma) {
+ glp_prob *lp;
+
+ int nb_coeffs = nb_points * _nb_clusters + nb_points * _nb_clusters;
+
+ int *ia = new int[nb_coeffs + 1];
+ int *ja = new int[nb_coeffs + 1];
+ scalar_t *ar = new scalar_t[nb_coeffs + 1];
+
+ scalar_t *nb_samples_per_class = new scalar_t[nb_classes];
+ for(int c = 0; c < nb_classes; c++) {
+ nb_samples_per_class[c] = 0.0;
+ }
+
+ for(int n = 0; n < nb_points; n++) {
+ nb_samples_per_class[labels[n]] += 1.0;
+ }
+
+ lp = glp_create_prob();
+
+ glp_set_prob_name(lp, "uninformative_lp_cluster_association");
+ glp_set_obj_dir(lp, GLP_MIN);
+
+ glp_add_rows(lp, nb_points + _nb_clusters * nb_classes);
+
+ for(int n = 1; n <= nb_points; n++) {
+ glp_set_row_bnds(lp, n, GLP_FX, 1.0, 1.0);
+ }
+
+ 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);
+ }
+ }
+
+ 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);
+
+ 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, r, dist);
+ glp_set_col_bnds(lp, r, GLP_DB, 0.0, 1.0);
+ }
+ }
+
+ int l = 1;
+
+ for(int n = 1; n <= nb_points; n++) {
+ for(int k = 1; k <= _nb_clusters; k++) {
+ int row = n;
+ ia[l] = row;
+ ja[l] = nb_points * (k - 1) + n;
+ ar[l] = 1.0;
+ l++;
+ }
+ }
+
+ 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;
+ for(int n = 1; n <= nb_points; n++) {
+ if(labels[n-1] == c - 1) {
+ ia[l] = row;
+ ja[l] = (k-1) * nb_points + n;
+ ar[l] = 1.0;
+ l++;
+ }
+ }
+ }
+ }
+
+ ASSERT(l == nb_coeffs + 1);
+
+ glp_load_matrix(lp, nb_coeffs, ia, ja, ar);
+
+ glp_simplex(lp, NULL);
+
+ scalar_t total_dist = glp_get_obj_val(lp);
+
+ for(int k = 1; k <= _nb_clusters; k++) {
+ for(int n = 1; n <= nb_points; n++) {
+ int i = n + nb_points * (k - 1);
+ gamma[n-1][k-1] = glp_get_col_prim(lp, i);
+ }
+ }
+
+ // { // ******************************* 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;
+ delete[] ar;
+ glp_delete_prob(lp);
+
+ return total_dist;
+}
+
+void Clusterer::baseline_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++) {
+ _cluster_means[k][d] = 0.0;
+ _cluster_var[k][d] = 0.0;
+ }
+
+ scalar_t sum_gamma = 0;
+ for(int n = 0; n < nb_points; n++) {
+ sum_gamma += gamma[n][k];
+ for(int d = 0; d < _dim; d++) {
+ _cluster_means[k][d] += gamma[n][k] * points[n][d];
+ _cluster_var[k][d] += gamma[n][k] * sq(points[n][d]);
+ }
+ }
+
+ ASSERT(sum_gamma >= 1);
+
+ 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);
+ } 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;
+ }
+ used[l] = 1;
+ }
+ delete[] used;
+}
+
+void Clusterer::train(int nb_clusters, int dim,
+ int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ int *cluster_associations) {
+ deallocate_array<scalar_t>(_cluster_means);
+ deallocate_array<scalar_t>(_cluster_var);
+ _nb_clusters = nb_clusters;
+ _dim = dim;
+ _cluster_means = allocate_array<scalar_t>(_nb_clusters, _dim);
+ _cluster_var = allocate_array<scalar_t>(_nb_clusters, _dim);
+
+ scalar_t **gammas = allocate_array<scalar_t>(nb_points, _nb_clusters);
+
+ if(nb_clusters > nb_points) abort();
+
+ initialize_clusters(nb_points, points);
+
+ scalar_t pred_total_distance, total_distance = FLT_MAX;
+ int nb_rounds = 0;
+
+ do {
+ pred_total_distance = total_distance;
+ 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);
+ cout << "TRAIN " << nb_rounds << " " << total_distance << endl;
+ baseline_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];
+ }
+ }
+ 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;
+ }
+ }
+ }
+
+ deallocate_array<scalar_t>(gammas);
+}
+
+int Clusterer::cluster(scalar_t *point) {
+ scalar_t lowest_dist = 0;
+ int associated_cluster = -1;
+
+ 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] - point[d]) / (2 * _cluster_var[k][d]);
+ dist += 0.5 * log(_cluster_var[k][d]);
+ ASSERT(!isnan(dist) && !isinf(dist));
+ }
+
+ if(k == 0 || dist <= lowest_dist) {
+ lowest_dist = dist;
+ associated_cluster = k;
+ }
+ }
+
+ ASSERT(associated_cluster >= 0);
+
+ return associated_cluster;
+}
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef CLUSTERER_H
+#define CLUSTERER_H
+
+#include "misc.h"
+#include "arrays.h"
+
+class Clusterer {
+public:
+ const static int max_nb_iterations = 10;
+ const static scalar_t min_iteration_improvement = 0.999;
+
+ int _nb_clusters;
+ int _dim;
+ scalar_t **_cluster_means, **_cluster_var;
+
+ void initialize_clusters(int nb_points, scalar_t **points);
+
+ scalar_t baseline_cluster_association(int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ scalar_t **gamma);
+
+ scalar_t baseline_lp_cluster_association(int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ scalar_t **gamma);
+
+ scalar_t uninformative_lp_cluster_association(int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ scalar_t **gamma);
+
+ void baseline_update_clusters(int nb_points, scalar_t **points, scalar_t **gamma);
+
+public:
+ Clusterer();
+ ~Clusterer();
+ void train(int nb_clusters, int dim,
+ int nb_points, scalar_t **points,
+ int nb_classes, int *labels,
+ int *cluster_associations);
+
+ int cluster(scalar_t *point);
+};
+
+#endif
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include <fstream>
+
+using namespace std;
+
+#include "misc.h"
+
+char *next_word(char *buffer, char *r, int buffer_size) {
+ char *s;
+ s = buffer;
+
+ if(r != 0) {
+ while((*r == ' ') || (*r == '\t') || (*r == ',')) r++;
+ if(*r == '"') {
+ r++;
+ while((*r != '"') && (*r != '\0') &&
+ (s<buffer+buffer_size-1))
+ *s++ = *r++;
+ if(*r == '"') r++;
+ } else {
+ while((*r != '\r') && (*r != '\n') && (*r != '\0') &&
+ (*r != '\t') && (*r != ' ') && (*r != ',')) {
+ if(s == buffer + buffer_size) {
+ cerr << "Buffer overflow in next_word." << endl;
+ exit(1);
+ }
+ *s++ = *r++;
+ }
+ }
+
+ while((*r == ' ') || (*r == '\t') || (*r == ',')) r++;
+ if((*r == '\0') || (*r=='\r') || (*r=='\n')) r = 0;
+ }
+ *s = '\0';
+
+ return r;
+}
+
+scalar_t discrete_entropy(int *n, int nb) {
+ scalar_t s = 0, t = 0;
+ for(int k = 0; k < nb; k++) if(n[k] > 0) {
+ s += n[k] * log(scalar_t(n[k]));
+ t += n[k];
+ }
+ return (log(t) - s/scalar_t(t))/log(2.0);
+}
+
+void random_permutation(int *val, int nb) {
+ for(int k = 0; k < nb; k++) val[k] = k;
+ int i, t;
+ for(int k = 0; k < nb - 1; k++) {
+ i = int(drand48() * (nb - k)) + k;
+ t = val[i];
+ val[i] = val[k];
+ val[k] = t;
+ }
+}
+
+void tag_subset(bool *val, int nb_total, int nb_to_tag) {
+ ASSERT(nb_to_tag <= nb_total);
+ int index[nb_total];
+ random_permutation(index, nb_total);
+ for(int n = 0; n < nb_total; n++) val[n] = false;
+ for(int n = 0; n < nb_to_tag; n++) val[index[n]] = true;
+}
+
+int compare_couple(const void *a, const void *b) {
+ if(((Couple *) a)->value < ((Couple *) b)->value) return -1;
+ else if(((Couple *) a)->value > ((Couple *) b)->value) return 1;
+ else return 0;
+}
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef MISC_H
+#define MISC_H
+
+#include <iostream>
+#include <cmath>
+#include <fstream>
+#include <cfloat>
+#include <stdlib.h>
+#include <string.h>
+
+using namespace std;
+
+typedef double scalar_t;
+// typedef float scalar_t;
+
+const int buffer_size = 1024;
+
+using namespace std;
+
+#ifdef DEBUG
+#define ASSERT(x) if(!(x)) { \
+ std::cerr << "ASSERT FAILED IN " << __FILE__ << ":" << __LINE__ << endl; \
+ abort(); \
+}
+#else
+#define ASSERT(x)
+#endif
+
+template<class T>
+T smooth_min(T x, T y) {
+ T z = exp(x - y);
+ return 0.5 * (x + y - (x - y)/(1 + 1/z) - (y - x)/(1 + z));
+}
+
+template <class T>
+void write_var(ostream *os, const T *x) { os->write((char *) x, sizeof(T)); }
+
+template <class T>
+void read_var(istream *is, T *x) { is->read((char *) x, sizeof(T)); }
+
+template <class T>
+void grow(int *nb_max, int nb, T** current, int factor) {
+ ASSERT(*nb_max > 0);
+ if(nb == *nb_max) {
+ T *tmp = new T[*nb_max * factor];
+ memcpy(tmp, *current, *nb_max * sizeof(T));
+ delete[] *current;
+ *current = tmp;
+ *nb_max *= factor;
+ }
+}
+
+template <class T>
+inline T sq(T x) {
+ return x * x;
+}
+
+inline scalar_t log2(scalar_t x) {
+ return log(x)/log(2.0);
+}
+
+inline scalar_t xi(scalar_t x) {
+ if(x <= 0.0) return 0.0;
+ else return - x * log(x)/log(2.0);
+}
+
+scalar_t discrete_entropy(int *n, int nb);
+
+char *basename(char *name);
+
+char *next_word(char *buffer, char *r, int buffer_size);
+
+void random_permutation(int *val, int nb);
+void tag_subset(bool *val, int nb_total, int nb_to_tag);
+
+struct Couple {
+ int index;
+ scalar_t value;
+};
+
+int compare_couple(const void *a, const void *b);
+
+#endif
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#include "sample_set.h"
+#include "arrays.h"
+
+SampleSet::SampleSet() {
+ labels = 0;
+ points = 0;
+}
+
+SampleSet::~SampleSet() {
+ delete[] labels;
+ deallocate_array<scalar_t>(points);
+}
+
+void SampleSet::resize(int d, int np) {
+ delete[] labels;
+ deallocate_array<scalar_t>(points);
+ nb_points = np;
+ dim = d;
+ points = allocate_array<scalar_t>(nb_points, dim);
+ labels = new int[nb_points];
+}
+
--- /dev/null
+/*
+ * clueless-kmean is a variant of k-mean which enforces balanced
+ * distribution of classes in every cluster
+ *
+ * Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+ * Written by Francois Fleuret <francois.fleuret@idiap.ch>
+ *
+ * This file is part of clueless-kmean.
+ *
+ * clueless-kmean is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * version 3 as published by the Free Software Foundation.
+ *
+ * clueless-kmean is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with selector. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef SAMPLE_SET_H
+#define SAMPLE_SET_H
+
+#include "misc.h"
+
+class SampleSet {
+public:
+ int nb_points, dim;
+ int nb_classes, *labels;
+ scalar_t **points;
+
+ SampleSet();
+ ~SampleSet();
+ void resize(int d, int np);
+};
+
+#endif
--- /dev/null
+#!/bin/bash
+
+
+# clueless-kmean is a variant of k-mean which enforces balanced
+# distribution of classes in every cluster
+#
+# Copyright (c) 2013 Idiap Research Institute, http://www.idiap.ch/
+# Written by Francois Fleuret <francois.fleuret@idiap.ch>
+#
+# This file is part of clueless-kmean.
+#
+# clueless-kmean is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# version 3 as published by the Free Software Foundation.
+#
+# clueless-kmean is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with selector. If not, see <http://www.gnu.org/licenses/>.
+
+set -e
+
+make -j -k
+
+./clueless-kmean
+
+CLUSTER1=($(grep ^0 clusters.dat))
+CLUSTER2=($(grep ^1 clusters.dat))
+CLUSTER3=($(grep ^2 clusters.dat))
+
+gnuplot <<EOF
+set terminal pngcairo truecolor size 1024,768
+set output "result.png"
+set size ratio 1
+set key out vert
+set key left top
+set object 1 ellipse center ${CLUSTER1[1]}, ${CLUSTER1[2]} size ${CLUSTER1[3]}, ${CLUSTER1[4]} angle 0 front fs empty bo 0 lw 1
+set object 2 ellipse center ${CLUSTER2[1]}, ${CLUSTER2[2]} size ${CLUSTER2[3]}, ${CLUSTER2[4]} angle 0 front fs empty bo 0 lw 1
+set object 3 ellipse center ${CLUSTER3[1]}, ${CLUSTER3[2]} size ${CLUSTER3[3]}, ${CLUSTER3[4]} angle 0 front fs empty bo 0 lw 1
+plot [-1.2:1.2][-1.2:1.2] "< grep ^0 associated_clusters.dat" using 2:3 w p lc rgb "#e00000" pt 6 ps 2.0 title "Cluster 1", \
+ "< grep ^1 associated_clusters.dat" using 2:3 w p lc rgb "#00c000" pt 6 ps 2.0 title "Cluster 2", \
+ "< grep ^2 associated_clusters.dat" using 2:3 w p lc rgb "#0000c0" pt 6 ps 2.0 title "Cluster 3", \
+ "< grep ^0 points.dat" using 2:3 w p lc rgb "#e00000" pt 7 ps 1.0 title "Class 1", \
+ "< grep ^1 points.dat" using 2:3 w p lc rgb "#00c000" pt 7 ps 1.0 title "Class 2", \
+ "< grep ^0 clusters.dat" using 2:3 w p lc rgb "#ffffff" pt 2 lw 9 ps 4.0 notitle, \
+ "< grep ^0 clusters.dat" using 2:3 w p lc rgb "#e00000" pt 2 lw 4 ps 4.0 title "Centroid 1", \
+ "< grep ^1 clusters.dat" using 2:3 w p lc rgb "#ffffff" pt 2 lw 9 ps 4.0 notitle, \
+ "< grep ^1 clusters.dat" using 2:3 w p lc rgb "#00c000" pt 2 lw 4 ps 4.0 title "Centroid 2", \
+ "< grep ^2 clusters.dat" using 2:3 w p lc rgb "#ffffff" pt 2 lw 9 ps 4.0 notitle, \
+ "< grep ^2 clusters.dat" using 2:3 w p lc rgb "#0000c0" pt 2 lw 4 ps 4.0 title "Centroid 3"
+EOF
projects / clueless-kmeans.git / commitdiff