Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

main.c 6.6KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322
  1. #include <assert.h>
  2. #include <math.h>
  3. #include <stdbool.h>
  4. #include <stdio.h>
  5. #include <stdlib.h>
  6. #include <string.h>
  7. #include <time.h>
  8. #include <unistd.h>
  9. #include "benchmark-data.h"
  10. /* Types */
  11. typedef int *solution;
  12. /* Functions */
  13. double distance(const vector u, const vector v);
  14. bool trouver_solution_initiale(void);
  15. bool solution_valide(const solution s);
  16. void print_solution(const solution s, const char *name);
  17. vector *centres_gravite(const solution s);
  18. int solution_voisine(solution *dest, const solution src);
  19. double fonction_objective(const solution s, double *inter, double *intra);
  20. static double calculer_inter(const vector *cgs);
  21. static double calculer_intra(const solution s, const vector *cgs);
  22. /* Variables */
  23. solution solution_opt;
  24. solution new_solution;
  25. int n_clusters = 5;
  26. int n_iterations = 1000;
  27. int main(int argc, char *argv[])
  28. {
  29. int c;
  30. int i;
  31. double f0, inter, intra;;
  32. while ((c = getopt(argc, argv, ":c:i:h")) != -1) {
  33. int k;
  34. switch (c) {
  35. case 'c':
  36. k = atoi(optarg);
  37. if (k > 1)
  38. n_clusters = k;
  39. break;
  40. case 'i':
  41. k = atoi(optarg);
  42. if (k > 0)
  43. n_iterations = k;
  44. break;
  45. case 'h':
  46. printf("Usage: %s [-c N_CLUSTERS] [-i N_ITERATIONS]", argv[0]);
  47. exit(EXIT_SUCCESS);
  48. break;
  49. default:
  50. exit(EXIT_FAILURE);
  51. }
  52. }
  53. printf("Using %d clusters\n", n_clusters);
  54. /* Init benchmark_data */
  55. init_benchmark_data();
  56. srand(time(NULL));
  57. if (!(solution_opt = calloc(N_VECTORS, sizeof(solution_opt[0]))))
  58. abort();
  59. /* Solution initiale */
  60. while (!trouver_solution_initiale())
  61. continue; /* do it again */
  62. f0 = fonction_objective(solution_opt, &inter, &intra);
  63. print_solution(solution_opt, "Solution initiale");
  64. printf("F = %f - %f = %f\n", inter, intra, f0);
  65. for (i = 0; i < n_iterations; i++) {
  66. double f1;
  67. while (!solution_voisine(&new_solution, solution_opt))
  68. continue; /* do it again */
  69. if ((f1 = fonction_objective(new_solution, &inter, &intra)) > f0) {
  70. print_solution(new_solution, "Nouvelle solution optimale trouvée");
  71. printf("F = %f - %f = %f (amélioration de %f%% par rapport à la dernière solution %f)\n",
  72. inter, intra, f1, fabs((f1 - f0) / f0 * 100.0), f0);
  73. free(solution_opt);
  74. solution_opt = new_solution;
  75. new_solution = NULL;
  76. f0 = f1;
  77. }
  78. }
  79. /* Free solution data */
  80. free(solution_opt);
  81. free(new_solution);
  82. /* Free benchmark data */
  83. for (i = 0; i < (int) N_VECTORS; i++)
  84. free(benchmark_data[i]);
  85. free(benchmark_data);
  86. exit(EXIT_SUCCESS);
  87. }
  88. double distance(const vector u, const vector v)
  89. {
  90. size_t i = 0;
  91. double sum = 0.0;
  92. if (!u || !v) {
  93. fputs("distance(): u or v is null.\n", stderr);
  94. abort();
  95. }
  96. /* d(u, v) = sqrt( (u1 - v1)^2 + (u2 - v2)^2 + ... + (un - vn)^2 ) */
  97. for (; i < VECTOR_SIZE; i++)
  98. sum += pow(u[i] - v[i], 2);
  99. return sqrt(sum);
  100. }
  101. bool trouver_solution_initiale(void)
  102. {
  103. size_t i;
  104. for (i = 0; i < N_VECTORS; i++)
  105. solution_opt[i] = rand() % n_clusters;
  106. return solution_valide(solution_opt);
  107. }
  108. bool solution_valide(const solution s)
  109. {
  110. /* This function returns true if the solution that it found is "valid",
  111. * i.e., it has at least one element in each cluster. */
  112. bool *cluster_empty = calloc(n_clusters, sizeof *cluster_empty);
  113. bool valid = true;
  114. int i;
  115. for (i = 0; i < n_clusters; i++)
  116. cluster_empty[i] = true;
  117. for (i = 0; (size_t) i < N_VECTORS; i++)
  118. cluster_empty[s[i]] = false;
  119. for (i = 0; valid && i < n_clusters; i++)
  120. valid = (valid && !cluster_empty[i]);
  121. free(cluster_empty);
  122. return valid;
  123. }
  124. void print_solution(const solution s, const char *name)
  125. {
  126. size_t i = 0;
  127. printf("%s: [\n", name ? name : "Solution");
  128. while (i < N_VECTORS) {
  129. const char *end;
  130. if (i + 1 == N_VECTORS)
  131. end = "]\n";
  132. else if ((i + 1) % 20 == 0)
  133. end = ",\n";
  134. else
  135. end = ", ";
  136. printf("%d%s", s[i++], end);
  137. }
  138. }
  139. vector *centres_gravite(const solution s)
  140. {
  141. /* Calculer le centre de gravité pour les vecteurs du cluster cluster selon la solution s */
  142. int cluster;
  143. size_t i = 0;
  144. /* cgs[n_clusters] est le centre de gravité global */
  145. vector *cgs = calloc(n_clusters + 1, sizeof(*cgs));
  146. int *vectors_in_cluster = calloc(n_clusters + 1, sizeof(*vectors_in_cluster));
  147. if (!cgs || !vectors_in_cluster)
  148. abort();
  149. for (cluster = 0; cluster < n_clusters + 1; cluster++) {
  150. if (!(cgs[cluster] = calloc(VECTOR_SIZE, sizeof(*cgs[cluster]))))
  151. abort();
  152. /* Initialiser à 0 */
  153. while (i < VECTOR_SIZE)
  154. cgs[cluster][i++] = 0;
  155. if (cluster < n_clusters)
  156. vectors_in_cluster[cluster] = 0;
  157. }
  158. for (i = 0; i < N_VECTORS; i++) {
  159. size_t j;
  160. cluster = s[i];
  161. ++vectors_in_cluster[cluster];
  162. for (j = 0; j < VECTOR_SIZE; j++) {
  163. cgs[cluster][j] += benchmark_data[i][j];
  164. cgs[n_clusters][j] += benchmark_data[i][j];
  165. }
  166. }
  167. vectors_in_cluster[n_clusters] = N_VECTORS;
  168. for (cluster = 0; cluster < n_clusters + 1; cluster++)
  169. for (i = 0; i < VECTOR_SIZE; i++)
  170. cgs[cluster][i] /= (double) vectors_in_cluster[cluster];
  171. free(vectors_in_cluster);
  172. return cgs;
  173. }
  174. int solution_voisine(solution *dest, const solution src)
  175. {
  176. /* This function changes about 10% of the solution randomly */
  177. size_t i;
  178. const int a = 10000; /* b = 10% of a */
  179. const int b = 1000;
  180. bool valid;
  181. if (!(*dest = calloc(N_VECTORS, sizeof(**dest))))
  182. abort();
  183. memcpy(*dest, src, N_VECTORS * sizeof(*src));
  184. assert(n_clusters > 1);
  185. for (i = 0; i < N_VECTORS; i++) {
  186. if (rand() % a < b) {
  187. int new_cluster = rand() % (n_clusters - 1);
  188. /* To avoid the case where the cluster isn't actually changed. */
  189. new_cluster += (new_cluster >= src[i]);
  190. (*dest)[i] = new_cluster;
  191. }
  192. }
  193. if (!(valid = solution_valide(*dest))) {
  194. free(*dest);
  195. *dest = NULL;
  196. }
  197. return valid;
  198. }
  199. double fonction_objective(const solution s, double *inter, double *intra)
  200. {
  201. int i;
  202. double **cgs = centres_gravite(s);
  203. if (!cgs)
  204. abort();
  205. *inter = calculer_inter(cgs);
  206. *intra = calculer_intra(s, cgs);
  207. for (i = 0; i < n_clusters + 1; i++)
  208. free(cgs[i]);
  209. free(cgs);
  210. return *inter - *intra;
  211. }
  212. static double calculer_inter(const vector *cgs)
  213. {
  214. double inter = 0.0;
  215. int cluster;
  216. for (cluster = 0; cluster < n_clusters; cluster++)
  217. inter += distance(cgs[n_clusters], cgs[cluster]);
  218. return inter / (double) n_clusters;
  219. }
  220. static double calculer_intra(const solution s, const vector *cgs)
  221. {
  222. double intra = 0.0;
  223. int cluster;
  224. for (cluster = 0; cluster < n_clusters; cluster++) {
  225. double intra_for_this_cluster = 0.0;
  226. int vectors_in_this_cluster = 0;
  227. size_t i;
  228. for (i = 0; i < N_VECTORS; i++) {
  229. if (s[i] != cluster)
  230. continue;
  231. ++vectors_in_this_cluster;
  232. intra_for_this_cluster += distance(cgs[cluster], benchmark_data[i]);
  233. }
  234. intra += intra_for_this_cluster / (double) vectors_in_this_cluster;
  235. }
  236. return intra / (double) n_clusters;
  237. }