meop/main.c

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include <unistd.h>

#include "benchmark-data.h"

/* TODO: Use typedefs for solution and vector instead of using int (*s)[230]
 * and double (*v)[42]. */

/* Types */
typedef int *solution;


/* Functions */
double distance(const vector *u, const vector *v);
bool trouver_solution_initiale(void);
bool solution_valide(const solution s);
void print_solution(const solution s, const char *name);
vector *centres_gravite(const solution s);
int solution_voisine(solution *dest, const solution src);
#if 0
double fonction_objective(const solution s);

static double calculer_inter(const solution s, vector *cgs);
#endif /* 0 */
static double calculer_intra(const solution s, vector *cgs);

solution solution_opt;
solution new_solution;

int n_clusters = 5;
int n_iterations = 1000;

int main(int argc, char *argv[])
{
	int c;
	int i;

	while ((c = getopt(argc, argv, ":c:i:h")) != -1) {
		int k;
		switch (c) {
		case 'c':
			k = atoi(optarg);
			if (k > 1)
				n_clusters = k;
			break;
		case 'i':
			k = atoi(optarg);
			if (k > 0)
				n_iterations = k;
			break;
		case 'h':
			printf("Usage: %s [-c N_CLUSTERS] [-i N_ITERATIONS]", argv[0]);
			exit(EXIT_SUCCESS);
			break;
		default:
			exit(EXIT_FAILURE);
		}
	}

	printf("Using %d clusters\n", n_clusters);

	srand(time(NULL));

	if (!(solution_opt = calloc(N_VECTORS, sizeof(solution_opt[0]))))
		abort();

	/* Solution initiale */
	while (!trouver_solution_initiale())
		continue; /* do it again */

	print_solution(solution_opt, "initiale");

	for (i = 0; i < n_iterations; i++) {

		while (!solution_voisine(&new_solution, solution_opt))
			continue; /* do it again */

		/* TODO: Calculer les centres de gravités, les distances moyennes inter-/intra-cluster,
		 * et faire le reste de l'algorithme */
		print_solution(new_solution, "new solution");
	}

	exit(EXIT_SUCCESS);
}

double distance(const vector *u, const vector *v)
{
	size_t i = 0;
	double sum = 0.0;

	/* d(u, v) = sqrt( (u1 - v1)^2 + (u2 - v2)^2 + ... + (un - vn)^2 ) */

	for (; i < VECTOR_SIZE; i++)
		sum += pow(u[i] - v[i], 2);

	return sqrt(sum);
}

bool trouver_solution_initiale(void)
{
	size_t i;

	for (i = 0; i < N_VECTORS; i++)
		solution_opt[i] = rand() % n_clusters;

	return solution_valide(solution_opt);
}

bool solution_valide(const solution s)
{
	/* This function returns true if the solution that it found is "valid",
	 * i.e., it has at least one element in each cluster. */

	bool *cluster_empty = calloc(n_clusters, sizeof *cluster_empty);
	bool valid = true;
	int i;

	for (i = 0; i < n_clusters; i++)
		cluster_empty[i] = true;

	for (i = 0; (size_t) i < N_VECTORS; i++)
		cluster_empty[s[i]] = false;

	for (i = 0; valid && i < n_clusters; i++)
		valid = (valid && !cluster_empty[i]);

	free(cluster_empty);

	return valid;
}

void print_solution(const solution s, const char *name)
{
	size_t i = 0;

	printf("Solution%s%s: [\n", name ? " " : "", name ? name : "");

	while (i < N_VECTORS) {
		const char *end;

		if (i + 1 == N_VECTORS)
			end = "]\n";
		else if ((i + 1) % 20 == 0)
			end = ",\n";
		else
			end = ", ";

		printf("%d%s", s[i++], end);
	}
}

vector *centres_gravite(const solution s)
{
	/* Calculer le centre de gravité pour les vecteurs du cluster cluster selon la solution s */

	int cluster;
	size_t i = 0;
	vector *cgs = calloc(n_clusters, sizeof(*cgs));
	int *vectors_in_cluster = calloc(n_clusters, sizeof(*vectors_in_cluster));

	if (!cgs || !vectors_in_cluster)
		abort();

	for (cluster = 0; cluster < n_clusters; cluster++) {
		if (!(cgs[cluster] = calloc(VECTOR_SIZE, sizeof(*cgs[cluster]))))
			abort();

		/* Initialiser à 0 */
		while (i < VECTOR_SIZE)
			cgs[cluster][i++] = 0;
		vectors_in_cluster[cluster] = 0;
	}

	for (i = 0; i < N_VECTORS; i++) {
		size_t j;

		cluster = s[i];
		++vectors_in_cluster[cluster];

		for (j = 0; j < VECTOR_SIZE; j++)
			cgs[cluster][j] += benchmark_data[i][j];
	}

	for (cluster = 0; cluster < n_clusters; cluster++)
		for (i = 0; i < VECTOR_SIZE; i++)
			cgs[cluster][i] /= (double) vectors_in_cluster[cluster];

	free(vectors_in_cluster);

	return cgs;
}

int solution_voisine(solution *dest, const solution src)
{
	/* This function changes about 10% of the solution randomly */

	size_t i;
	const int a = 10000; /* b = 10% of a */
	const int b =  1000;
	bool valid;

	if (!(*dest = calloc(N_VECTORS, sizeof(**dest))))
		abort();

	memcpy(*dest, src, N_VECTORS * sizeof(*src));

	assert(n_clusters > 1);

	for (i = 0; i < N_VECTORS; i++) {
		if (rand() % a < b) {
			int new_cluster = rand() % (n_clusters - 1);
			/* To avoid the case where the cluster isn't actually changed. */
			new_cluster += (new_cluster >= src[i]);
			(*dest)[i] = new_cluster;
		}
	}

	if (!(valid = solution_valide(*dest))) {
		free(*dest);
		*dest = NULL;
	}

	return valid;
}

#if 0
double fonction_objective(const solution s)
{
	/* NOTE: this is very ineffecient on processor time, because
	 * each of these functions re-calculates the same centre_gravite
	 * twice. */

	double **cgs = calloc(n_clusters, sizeof *cgs);
	double inter;
	double intra;


	if (!cgs)
		abort();

	inter = calculer_inter(s, cgs);
	inter = calculer_inter(s, cgs);

	return inter - intra;
}

static double calculer_inter(const solution s, vector *cgs)
{
	const size_t n = sizeof(*s) / sizeof((*s)[0]);
	double inter = 0.0;
	int cluster;

	for (cluster = 0; cluster < n_clusters; cluster++) {
		double intra_for_this_cluster = 0.0;
		int vectors_in_this_cluster = 0;
		size_t i;

		centre_gravite(cluster, s, &cg);
	}

}

static double calculer_intra(const solution s, vector *cgs)
{
	double intra = 0.0;
	int cluster;

	for (cluster = 0; cluster < n_clusters; cluster++) {
		double intra_for_this_cluster = 0.0;
		int vectors_in_this_cluster = 0;
		size_t i;

		for (i = 0; i < N_VECTORS; i++) {
			if (s[i] != cluster)
				continue;

			++vectors_in_this_cluster;
			intra_for_this_cluster += distance((const vector *) cgs[cluster],
					(const vector *) benchmark_data[i]);
		}

		intra += intra_for_this_cluster / (double) vectors_in_this_cluster;
	}

	return intra / (double) n_clusters;
}
#endif /* 0 */