#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"


/* 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);
double fonction_objective(const solution s);

static double calculer_inter(const vector *cgs);
static double calculer_intra(const solution s, const vector *cgs);


/* Variables */
solution solution_opt;
solution new_solution;

int n_clusters = 5;
int n_iterations = 1000;


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

	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);

	/* Init benchmark_data */
	init_benchmark_data();

	srand(time(NULL));

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

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

	f0 = fonction_objective(solution_opt);

	print_solution(solution_opt, "initiale");
	printf("F(solution) = %f\n", f0);

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

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

		if ((f1 = fonction_objective(new_solution)) > f0) {
			print_solution(new_solution, "optimale trouvée");
			printf("F(solution) = %f (%f%% improvement)\n", f1, fabs((f1 - f0) / f0 * 100.0));
			free(solution_opt);
			solution_opt = new_solution;
			new_solution = NULL;
			f0 = f1;
		}
	}

	/* Free solution data */
	free(solution_opt);
	free(new_solution);


	/* Free benchmark data */
	for (i = 0; i < (int) N_VECTORS; i++)
		free(benchmark_data[i]);
	free(benchmark_data);

	exit(EXIT_SUCCESS);
}

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

	if (!u || !v) {
		fputs("distance(): u or v is null.\n", stderr);
		abort();
	}

	/* 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;

	/* cgs[n_clusters] est le centre de gravité global */
	vector *cgs = calloc(n_clusters + 1, sizeof(*cgs));

	int *vectors_in_cluster = calloc(n_clusters + 1, sizeof(*vectors_in_cluster));

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

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

		/* Initialiser à 0 */
		while (i < VECTOR_SIZE)
			cgs[cluster][i++] = 0;
		if (cluster < n_clusters)
		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];
			cgs[n_clusters][j] += benchmark_data[i][j];
		}
	}

	vectors_in_cluster[n_clusters] = N_VECTORS;

	for (cluster = 0; cluster < n_clusters + 1; 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;
}

double fonction_objective(const solution s)
{
	int i;
	double **cgs = centres_gravite(s);
	double inter;
	double intra;


	if (!cgs)
		abort();

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

	for (i = 0; i < n_clusters + 1; i++)
		free(cgs[i]);

	free(cgs);

	return inter - intra;
}

static double calculer_inter(const vector *cgs)
{
	double inter = 0.0;
	int cluster;

	for (cluster = 0; cluster < n_clusters; cluster++)
		inter += distance(cgs[n_clusters], cgs[cluster]);

	return inter / (double) n_clusters;
}

static double calculer_intra(const solution s, const 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(cgs[cluster], benchmark_data[i]);
		}

		intra += intra_for_this_cluster / (double) vectors_in_this_cluster;
	}

	return intra / (double) n_clusters;
}