gts/partition.c

   1 /* GTS - Library for the manipulation of triangulated surfaces
   2  * Copyright (C) 1999 Stéphane Popinet
   3  *
   4  * This library is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Library General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  * License along with this library; if not, write to the
  16  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  * Boston, MA 02111-1307, USA.
  18  */
  19
  20 #include <math.h>
  21
  22 #include "gts.h"
  23
  24 /* #define DEBUG */
  25
  26 /* Graph partition */
  27
  28 /**
  29  * gts_graph_partition_edges_cut:
  30  * @partition: a list of @GtsGraph representing a partition.
  31  *
  32  * Returns: the number of edges cut by the partition.
  33  */
  34 guint gts_graph_partition_edges_cut (GSList * partition)
  35 {
  36   guint cuts = 0;
  37
  38   while (partition) {
  39     cuts += gts_graph_edges_cut (partition->data);
  40     partition = partition->next;
  41   }
  42
  43   return cuts/2;
  44 }
  45
  46 /**
  47  * gts_graph_partition_edges_cut_weight:
  48  * @partition: a list of @GtsGraph representing a partition.
  49  *
  50  * Returns: the total weight of the edges cut by the partition.
  51  */
  52 gfloat gts_graph_partition_edges_cut_weight (GSList * partition)
  53 {
  54   gfloat weight = 0.;
  55
  56   while (partition) {
  57     weight += gts_graph_edges_cut_weight (partition->data);
  58     partition = partition->next;
  59   }
  60
  61   return weight/2.;
  62 }
  63
  64 /**
  65  * gts_graph_partition_print_stats:
  66  * @partition: a list of @GtsGraph representing a partition.
  67  * @fp: a file pointer.
  68  *
  69  * Writes to @fp a summary of the properties of @partition.
  70  */
  71 void gts_graph_partition_print_stats (GSList * partition,
  72                                       FILE * fp)
  73 {
  74   GtsRange weight;
  75   GSList * i;
  76
  77   g_return_if_fail (partition != NULL);
  78   g_return_if_fail (fp != NULL);
  79
  80   gts_range_init (&weight);
  81   i = partition;
  82   while (i) {
  83     gts_range_add_value (&weight, gts_graph_weight (i->data));
  84     i = i->next;
  85   }
  86   gts_range_update (&weight);
  87
  88   fprintf (fp,
  89            "# parts: %d\n"
  90            "#   edge cuts: %5d edge cuts weight: %5g\n"
  91            "#   weight: ",
  92            g_slist_length (partition),
  93            gts_graph_partition_edges_cut (partition),
  94            gts_graph_partition_edges_cut_weight (partition));
  95   gts_range_print (&weight, fp);
  96   fputc ('\n', fp);
  97 }
  98
  99 /**
 100  * gts_graph_partition_balance:
 101  * @partition: a list of @GtsGraph representing a partition.
 102  *
 103  * Returns: the difference between the maximum and the minimum weight
 104  * of the graphs in @partition.
 105  */
 106 gfloat gts_graph_partition_balance (GSList * partition)
 107 {
 108   gfloat wmin = G_MAXFLOAT;
 109   gfloat wmax = - G_MAXFLOAT;
 110
 111   g_return_val_if_fail (partition != NULL, 0.);
 112
 113   while (partition) {
 114     gfloat weight = gts_graph_weight (partition->data);
 115     if (weight < wmin)
 116       wmin = weight;
 117     if (weight > wmax)
 118       wmax = weight;
 119     partition = partition->next;
 120   }
 121   return wmax - wmin;
 122 }
 123
 124 /**
 125  * gts_graph_partition_clone:
 126  * @partition: a list of @GtsGraph representing a partition.
 127  *
 128  * Returns: a new partition clone of @partition (i.e. a list of new
 129  * graphs clones of the graphs in @partition).
 130  */
 131 GSList * gts_graph_partition_clone (GSList * partition)
 132 {
 133   GSList * cparts = NULL;
 134
 135   while (partition) {
 136     cparts =
 137       g_slist_prepend (cparts,
 138                        gts_object_clone (GTS_OBJECT (partition->data)));
 139     partition = partition->next;
 140   }
 141   return cparts;
 142 }
 143
 144 /**
 145  * gts_graph_partition_destroy:
 146  * @partition: a list of @GtsGraph representing a partition.
 147  *
 148  * Destroys all the graphs in @partition and frees @partition.
 149  */
 150 void gts_graph_partition_destroy (GSList * partition)
 151 {
 152   GSList * i = partition;
 153
 154   while (i) {
 155     gts_object_destroy (GTS_OBJECT (i->data));
 156     i = i->next;
 157   }
 158   g_slist_free (partition);
 159 }
 160
 161 static void find_smallest_degree (GtsGNode * n, gpointer * data)
 162 {
 163   GtsGNode ** nmin = data[0];
 164   GtsGraph * g = data[1];
 165   guint * min = data[2];
 166   guint degree = gts_gnode_degree (n, g);
 167
 168   if (degree < *min) {
 169     *min = degree;
 170     *nmin = n;
 171   }
 172 }
 173
 174 static gint graph_comp_weight (GtsGraph * g1, GtsGraph * g2)
 175 {
 176   if (gts_graph_weight (g1) > gts_graph_weight (g2))
 177     return 1;
 178   return -1;
 179 }
 180
 181 static void partition_update (GSList * list, GtsGraph * g)
 182 {
 183   GSList * i;
 184   GtsGraph * g1;
 185   GtsHeap * size_heap;
 186   gboolean reinit = TRUE;
 187
 188   /* initialize traversals */
 189   i = list;
 190   while (i) {
 191     GtsGNode * seed = GTS_OBJECT (i->data)->reserved;
 192     GTS_OBJECT (seed)->reserved =
 193       gts_graph_traverse_new (g, seed, GTS_BREADTH_FIRST, reinit);
 194     reinit = FALSE;
 195     i = i->next;
 196   }
 197
 198   size_heap = gts_heap_new ((GCompareFunc) graph_comp_weight);
 199   i = list;
 200   while (i) {
 201     gts_heap_insert (size_heap, i->data);
 202     i = i->next;
 203   }
 204   while ((g1 = gts_heap_remove_top (size_heap))) {
 205     GtsGraphTraverse * t = GTS_OBJECT (GTS_OBJECT (g1)->reserved)->reserved;
 206     GtsGNode * n = gts_graph_traverse_next (t);
 207     if (n) {
 208       gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 209       gts_heap_insert (size_heap, g1);
 210     }
 211   }
 212   gts_heap_destroy (size_heap);
 213
 214   /* destroy traversals */
 215   i = list;
 216   while (i) {
 217     GtsGNode * seed = GTS_OBJECT (i->data)->reserved;
 218     gts_graph_traverse_destroy (GTS_OBJECT (seed)->reserved);
 219     GTS_OBJECT (seed)->reserved = NULL;
 220     i = i->next;
 221   }
 222 }
 223
 224 static void better_seed (GtsGNode * n, gpointer * data)
 225 {
 226   guint * sum = data[0];
 227   GtsGNode ** seed = data[1];
 228   GtsGraph * g = data[2];
 229   guint sum1 = gts_graph_distance_sum (g, n);
 230
 231   if (sum1 < *sum) {
 232     *sum = sum1;
 233     *seed = n;
 234   }
 235 }
 236
 237 static GtsGNode * graph_new_seed (GtsGraph * g, GtsGNode * seed)
 238 {
 239   guint sum = gts_graph_distance_sum (g, seed);
 240   gpointer data[3];
 241   GtsGNode * new_seed = seed;
 242
 243   data[0] = &sum;
 244   data[1] = &new_seed;
 245   data[2] = g;
 246   gts_gnode_foreach_neighbor (seed, g, (GtsFunc) better_seed, data);
 247
 248   return new_seed;
 249 }
 250
 251 /**
 252  * gts_graph_bubble_partition:
 253  * @g: a #GtsGraph.
 254  * @np: number of partitions.
 255  * @niter: the maximum number of iterations.
 256  * @step_info: a #GtsFunc or %NULL.
 257  * @data: user data to pass to @step_info.
 258  *
 259  * An implementation of the "bubble partitioning algorithm" of
 260  * Diekmann, Preis, Schlimbach and Walshaw (2000). The maximum number
 261  * of iteration on the positions of the graph growing seeds is
 262  * controlled by @niter.
 263  *
 264  * If not %NULL @step_info is called after each iteration on the seeds
 265  * positions passing the partition (a GSList) as argument.
 266  *
 267  * Returns: a list of @np new #GtsGraph representing the partition.
 268  */
 269 GSList * gts_graph_bubble_partition (GtsGraph * g,
 270                                      guint np,
 271                                      guint niter,
 272                                      GtsFunc step_info,
 273                                      gpointer data)
 274 {
 275   GSList * list = NULL, * seeds = NULL;
 276   GtsGNode * seed = NULL;
 277   guint min = G_MAXINT/2 - 1;
 278   gpointer info[3];
 279   GtsGraph * g1;
 280   gboolean changed = TRUE;
 281
 282   g_return_val_if_fail (g != NULL, NULL);
 283   g_return_val_if_fail (np > 0, NULL);
 284
 285   info[0] = &seed;
 286   info[1] = g;
 287   info[2] = &min;
 288   gts_container_foreach (GTS_CONTAINER (g),
 289                          (GtsFunc) find_smallest_degree,
 290                          info);
 291   if (seed == NULL)
 292     return NULL;
 293
 294   g1 = GTS_GRAPH (gts_object_new (GTS_OBJECT (g)->klass));
 295   gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
 296   list = g_slist_prepend (list, g1);
 297   GTS_OBJECT (g1)->reserved = seed;
 298   seeds = g_slist_prepend (seeds, seed);
 299
 300   while (--np && seed)
 301     if ((seed = gts_graph_farthest (g, seeds))) {
 302       g1 = GTS_GRAPH (gts_object_new (GTS_OBJECT (g)->klass));
 303       gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
 304       list = g_slist_prepend (list, g1);
 305       GTS_OBJECT (g1)->reserved = seed;
 306       seeds = g_slist_prepend (seeds, seed);
 307     }
 308   g_slist_free (seeds);
 309
 310   partition_update (list, g);
 311
 312   while (changed && niter--) {
 313     GSList * i;
 314
 315     changed = FALSE;
 316     i = list;
 317     while (i) {
 318       GtsGraph * g1 = i->data;
 319       GtsGNode * seed = GTS_OBJECT (g1)->reserved;
 320       GtsGNode * new_seed = graph_new_seed (g1, seed);
 321       if (new_seed != seed) {
 322         changed = TRUE;
 323         GTS_OBJECT (g1)->reserved = new_seed;
 324       }
 325       i = i->next;
 326     }
 327
 328     if (changed) {
 329       i = list;
 330       while (i) {
 331         GtsGraph * g1 = i->data;
 332         GtsGNode * seed = GTS_OBJECT (g1)->reserved;
 333
 334         gts_object_destroy (GTS_OBJECT (g1));
 335         i->data = g1 = GTS_GRAPH (gts_object_new (GTS_OBJECT (g)->klass));
 336         gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
 337         GTS_OBJECT (g1)->reserved = seed;
 338         i = i->next;
 339       }
 340       partition_update (list, g);
 341       if (step_info)
 342         (* step_info) (list, data);
 343     }
 344   }
 345   g_slist_foreach (list, (GFunc) gts_object_reset_reserved, NULL);
 346   return list;
 347 }
 348
 349 /* Graph bisection */
 350
 351 static gdouble node_cost (GtsGNode * n, gpointer * data)
 352 {
 353   GtsGraph * g = data[0];
 354   GtsGraph * g1 = data[1];
 355   GSList * i = GTS_SLIST_CONTAINER (n)->items;
 356   gdouble cost = 0.;
 357
 358   while (i) {
 359     GtsGEdge * e = i->data;
 360     GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, e);
 361
 362     if (gts_containee_is_contained (GTS_CONTAINEE (n1), GTS_CONTAINER (g))) {
 363       if (gts_containee_is_contained (GTS_CONTAINEE (n1), GTS_CONTAINER (g1)))
 364         cost -= gts_gedge_weight (e);
 365       else
 366         cost += gts_gedge_weight (e);
 367     }
 368     i = i->next;
 369   }
 370
 371   return cost;
 372 }
 373
 374 static void add_neighbor (GtsGNode * n, GtsEHeap * heap)
 375 {
 376   if (GTS_OBJECT (n)->reserved == n)
 377     return;
 378   if (GTS_OBJECT (n)->reserved)
 379     gts_eheap_remove (heap, GTS_OBJECT (n)->reserved);
 380   GTS_OBJECT (n)->reserved = gts_eheap_insert (heap, n);
 381 }
 382
 383 static void add_unused (GtsGNode * n, GtsGraph * g2)
 384 {
 385   if (GTS_OBJECT (n)->reserved)
 386     GTS_OBJECT (n)->reserved = NULL;
 387   else
 388     gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
 389 }
 390
 391 static gdouble degree_cost (GtsGNode * n, GtsGraph * g)
 392 {
 393   return gts_gnode_degree (n, g);
 394 }
 395
 396 static void add_seed (GtsGNode * n, GtsEHeap * heap)
 397 {
 398   gts_eheap_insert (heap, n);
 399 }
 400
 401 static void boundary_node1 (GtsGNode * n, GtsGraphBisection * bg)
 402 {
 403   GSList * i = GTS_SLIST_CONTAINER (n)->items;
 404
 405   while (i) {
 406     GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
 407     if (gts_containee_is_contained (GTS_CONTAINEE (n1),
 408                                     GTS_CONTAINER (bg->g2))) {
 409       g_hash_table_insert (bg->bg1, n, n1);
 410       return;
 411     }
 412     i = i->next;
 413   }
 414 }
 415
 416 static void boundary_node2 (GtsGNode * n, GtsGraphBisection * bg)
 417 {
 418   GSList * i = GTS_SLIST_CONTAINER (n)->items;
 419
 420   while (i) {
 421     GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
 422     if (gts_containee_is_contained (GTS_CONTAINEE (n1),
 423                                     GTS_CONTAINER (bg->g1))) {
 424       g_hash_table_insert (bg->bg2, n, n1);
 425       return;
 426     }
 427     i = i->next;
 428   }
 429 }
 430
 431 static void check_bg (GtsGNode * n, gpointer * data)
 432 {
 433   GHashTable * bg = data[0];
 434   GtsGraph * g = data[1];
 435   gboolean * ok = data[2];
 436   guint * nb = data[3];
 437   guint nn = gts_gnode_degree (n, g);
 438
 439   if (nn > 0)
 440     (*nb)++;
 441   if ((nn > 0 && !g_hash_table_lookup (bg, n)) ||
 442       (nn == 0 && g_hash_table_lookup (bg, n))) {
 443     g_warning ("nn: %d lookup: %p\n",
 444                nn, g_hash_table_lookup (bg, n));
 445     *ok = FALSE;
 446   }
 447 }
 448
 449 /**
 450  * gts_graph_bisection_check:
 451  * @bg: a #GtsGraphBisection.
 452  *
 453  * Checks that the boundary of @bg is correctly defined (used for
 454  * debugging purposes).
 455  *
 456  * Returns: %TRUE if @bg is ok, %FALSE otherwise.
 457  */
 458 gboolean gts_graph_bisection_check (GtsGraphBisection * bg)
 459 {
 460   gboolean ok = TRUE;
 461   guint nb;
 462   gpointer data[4];
 463
 464   g_return_val_if_fail (bg != NULL, FALSE);
 465
 466   nb = 0;
 467   data[0] = bg->bg1;
 468   data[1] = bg->g2;
 469   data[2] = &ok;
 470   data[3] = &nb;
 471   gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) check_bg, data);
 472   g_return_val_if_fail (g_hash_table_size (bg->bg1) == nb, FALSE);
 473
 474   nb = 0;
 475   data[0] = bg->bg2;
 476   data[1] = bg->g1;
 477   gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) check_bg, data);
 478   g_return_val_if_fail (g_hash_table_size (bg->bg2) == nb, FALSE);
 479
 480   return ok;
 481 }
 482
 483 /**
 484  * gts_graph_ggg_bisection:
 485  * @g: a #GtsGraph.
 486  * @ntry: the number of randomly selected initial seeds.
 487  *
 488  * An implementation of the "Greedy Graph Growing" algorithm of
 489  * Karypis and Kumar (1997).
 490  *
 491  * @ntry randomly chosen seeds are used and the best partition is retained.
 492  *
 493  * Returns: a new #GtsGraphBisection of @g.
 494  */
 495 GtsGraphBisection * gts_graph_ggg_bisection (GtsGraph * g, guint ntry)
 496 {
 497   gfloat size, bestcost = G_MAXFLOAT, smin;
 498   GtsGraph * bestg1 = NULL, * bestg2 = NULL;
 499   gboolean balanced = FALSE;
 500   GtsEHeap * degree_heap;
 501   GtsGNode * seed;
 502   GtsGraphBisection * bg;
 503
 504   g_return_val_if_fail (g != NULL, NULL);
 505
 506   bg = g_malloc (sizeof (GtsGraphBisection));
 507   bg->g = g;
 508
 509   size = gts_graph_weight (g)/2.;
 510   smin = 0.9*size;
 511
 512   degree_heap = gts_eheap_new ((GtsKeyFunc) degree_cost, g);
 513   gts_eheap_freeze (degree_heap);
 514   gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_seed, degree_heap);
 515   gts_eheap_thaw (degree_heap);
 516
 517   while (ntry && ((seed = gts_eheap_remove_top (degree_heap, NULL)))) {
 518     GtsGraph * g1, * g2;
 519     GtsGNode * n;
 520     gdouble cost;
 521     gpointer data[2];
 522     GtsEHeap * heap;
 523
 524     g1 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
 525                         g->node_class, g->edge_class);
 526     g2 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
 527                         g->node_class, g->edge_class);
 528
 529     data[0] = g;
 530     data[1] = g1;
 531     heap = gts_eheap_new ((GtsKeyFunc) node_cost, data);
 532
 533     gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
 534     GTS_OBJECT (seed)->reserved = seed;
 535     gts_gnode_foreach_neighbor (seed, g, (GtsFunc) add_neighbor, heap);
 536
 537     while ((n = gts_eheap_remove_top (heap, &cost)))
 538       if (gts_graph_weight (g1) + gts_gnode_weight (n) <= size) {
 539         gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 540         GTS_OBJECT (n)->reserved = n;
 541         gts_gnode_foreach_neighbor (n, g, (GtsFunc) add_neighbor, heap);
 542       }
 543       else
 544         GTS_OBJECT (n)->reserved = NULL;
 545     gts_eheap_destroy (heap);
 546
 547     gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_unused, g2);
 548
 549     cost = gts_graph_edges_cut_weight (g1);
 550     if (!bestg1 ||
 551         (!balanced && gts_graph_weight (g1) >= smin) ||
 552         (cost < bestcost && gts_graph_weight (g1) >= smin)) {
 553       if (bestg1)
 554         bestcost = cost;
 555       if (bestg1)
 556         gts_object_destroy (GTS_OBJECT (bestg1));
 557       if (bestg2)
 558         gts_object_destroy (GTS_OBJECT (bestg2));
 559       bestg1 = g1;
 560       bestg2 = g2;
 561       if (gts_graph_weight (g1) >= smin)
 562         balanced = TRUE;
 563     }
 564     else {
 565       gts_object_destroy (GTS_OBJECT (g1));
 566       gts_object_destroy (GTS_OBJECT (g2));
 567     }
 568
 569     ntry--;
 570   }
 571   gts_eheap_destroy (degree_heap);
 572
 573 #ifdef DEBUG
 574   fprintf (stderr, "bestcost: %5g g1: %5g|%5d g2: %5g|%5d\n",
 575            bestcost,
 576            gts_graph_weight (bestg1),
 577            gts_container_size (GTS_CONTAINER (bestg1)),
 578            gts_graph_weight (bestg2),
 579            gts_container_size (GTS_CONTAINER (bestg2)));
 580 #endif
 581
 582   g_assert (bestg1 != NULL);
 583   bg->g1 = bestg1;
 584   g_assert (bestg2 != NULL);
 585   bg->g2 = bestg2;
 586
 587   /* boundary nodes */
 588   bg->bg1 = g_hash_table_new (NULL, NULL);
 589   gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) boundary_node1, bg);
 590   bg->bg2 = g_hash_table_new (NULL, NULL);
 591   gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) boundary_node2, bg);
 592
 593   return bg;
 594 }
 595
 596 /**
 597  * gts_graph_bfgg_bisection:
 598  * @g: a #GtsGraph.
 599  * @ntry: the number of randomly selected initial seeds.
 600  *
 601  * An implementation of a "Breadth-First Graph Growing" algorithm.
 602  *
 603  * @ntry randomly chosen seeds are used and the best partition is retained.
 604  *
 605  * Returns: a new #GtsGraphBisection of @g.
 606  */
 607 GtsGraphBisection * gts_graph_bfgg_bisection (GtsGraph * g, guint ntry)
 608 {
 609   gfloat size, bestcost = G_MAXFLOAT, smin;
 610   GtsGraph * bestg1 = NULL, * bestg2 = NULL;
 611   GtsEHeap * degree_heap;
 612   GtsGNode * seed;
 613   GtsGraphBisection * bg;
 614
 615   g_return_val_if_fail (g != NULL, NULL);
 616
 617   bg = g_malloc (sizeof (GtsGraphBisection));
 618   bg->g = g;
 619
 620   size = gts_graph_weight (g)/2.;
 621   smin = 0.9*size;
 622
 623   degree_heap = gts_eheap_new ((GtsKeyFunc) degree_cost, g);
 624   gts_eheap_freeze (degree_heap);
 625   gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_seed, degree_heap);
 626   gts_eheap_thaw (degree_heap);
 627
 628   while (ntry && ((seed = gts_eheap_remove_top (degree_heap, NULL)))) {
 629     GtsGraph * g1, * g2;
 630     GtsGNode * n;
 631     gdouble cost;
 632     GtsGraphTraverse * t = gts_graph_traverse_new (g, seed,
 633                                                    GTS_BREADTH_FIRST, TRUE);
 634
 635     g1 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
 636                         g->node_class, g->edge_class);
 637     g2 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
 638                         g->node_class, g->edge_class);
 639
 640     while ((n = gts_graph_traverse_next (t)))
 641       if (gts_graph_weight (g1) + gts_gnode_weight (n) <= size) {
 642         gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 643         GTS_OBJECT (n)->reserved = n;
 644       }
 645     gts_graph_traverse_destroy (t);
 646
 647     gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_unused, g2);
 648
 649     cost = gts_graph_edges_cut_weight (g1);
 650     if (!bestg1 || (cost < bestcost && gts_graph_weight (g1) >= smin)) {
 651       if (bestg1)
 652         bestcost = cost;
 653       if (bestg1)
 654         gts_object_destroy (GTS_OBJECT (bestg1));
 655       if (bestg2)
 656         gts_object_destroy (GTS_OBJECT (bestg2));
 657       bestg1 = g1;
 658       bestg2 = g2;
 659     }
 660     else {
 661       gts_object_destroy (GTS_OBJECT (g1));
 662       gts_object_destroy (GTS_OBJECT (g2));
 663     }
 664
 665     ntry--;
 666   }
 667   gts_eheap_destroy (degree_heap);
 668
 669 #ifdef DEBUG
 670   fprintf (stderr, "bestcost: %5g g1: %5g|%5d g2: %5g|%5d\n",
 671            bestcost,
 672            gts_graph_weight (bestg1),
 673            gts_container_size (GTS_CONTAINER (bestg1)),
 674            gts_graph_weight (bestg2),
 675            gts_container_size (GTS_CONTAINER (bestg2)));
 676 #endif
 677
 678   bg->g1 = bestg1;
 679   bg->g2 = bestg2;
 680
 681   /* boundary nodes */
 682   bg->bg1 = g_hash_table_new (NULL, NULL);
 683   gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) boundary_node1, bg);
 684   bg->bg2 = g_hash_table_new (NULL, NULL);
 685   gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) boundary_node2, bg);
 686
 687   return bg;
 688 }
 689
 690 static gdouble node_move_cost1 (GtsGNode * n, GtsGraphBisection * bg)
 691 {
 692   return gts_gnode_move_cost (n, bg->g1, bg->g2);
 693 }
 694
 695 static gdouble node_move_cost2 (GtsGNode * n, GtsGraphBisection * bg)
 696 {
 697   return gts_gnode_move_cost (n, bg->g2, bg->g1);
 698 }
 699
 700 static void build_heap (GtsGNode * n, GtsEHeap * heap)
 701 {
 702   GTS_OBJECT (n)->reserved = gts_eheap_insert (heap, n);
 703 }
 704
 705 /**
 706  * gts_graph_bisection_kl_refine:
 707  * @bg: a #GtsGraphBisection.
 708  * @mmax: the maximum number of unsuccessful successive moves.
 709  *
 710  * An implementation of the simplified Kernighan-Lin algorithm for
 711  * graph bisection refinement as described in Karypis and Kumar
 712  * (1997).
 713  *
 714  * The algorithm stops if @mmax consecutive modes do not lead to a
 715  * decrease in the number of edges cut. This last @mmax moves are
 716  * undone.
 717  *
 718  * Returns: the decrease in the weight of the edges cut by the bisection.
 719  */
 720 gdouble gts_graph_bisection_kl_refine (GtsGraphBisection * bg,
 721                                        guint mmax)
 722 {
 723   GtsEHeap * h1, * h2;
 724   GtsGNode * n;
 725   guint nm = 0, i;
 726   GtsGNode ** moves;
 727   gdouble bestcost = 0., totalcost = 0., best_balance;
 728
 729   g_return_val_if_fail (bg != NULL, 0.);
 730   g_return_val_if_fail (mmax > 0, 0.);
 731
 732   h1 = gts_eheap_new ((GtsKeyFunc) node_move_cost1, bg);
 733   gts_eheap_freeze (h1);
 734   gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) build_heap, h1);
 735   gts_eheap_thaw (h1);
 736
 737   h2 = gts_eheap_new ((GtsKeyFunc) node_move_cost2, bg);
 738   gts_eheap_freeze (h2);
 739   gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) build_heap, h2);
 740   gts_eheap_thaw (h2);
 741
 742   moves = g_malloc (sizeof (GtsGNode *)*mmax);
 743   best_balance = fabs (gts_graph_weight (bg->g1) - gts_graph_weight (bg->g2));
 744
 745   do {
 746     GtsGraph * g1, * g2;
 747     gdouble cost;
 748
 749     if (gts_graph_weight (bg->g1) > gts_graph_weight (bg->g2)) {
 750       n = gts_eheap_remove_top (h1, &cost);
 751       g1 = bg->g1;
 752       g2 = bg->g2;
 753     }
 754     else {
 755       n = gts_eheap_remove_top (h2, &cost);
 756       g1 = bg->g2;
 757       g2 = bg->g1;
 758     }
 759     if (n) {
 760       GSList * i;
 761
 762       GTS_OBJECT (n)->reserved = NULL;
 763       gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
 764       gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 765
 766       totalcost += cost;
 767       if (totalcost < bestcost) {
 768         bestcost = totalcost;
 769         nm = 0;
 770       }
 771       else if (totalcost == bestcost) {
 772         gdouble balance = fabs (gts_graph_weight (g1) - gts_graph_weight (g2));
 773
 774         if (balance < best_balance) {
 775           best_balance = balance;
 776           nm = 0;
 777         }
 778       }
 779       else
 780         moves[nm++] = n;
 781
 782       i = GTS_SLIST_CONTAINER (n)->items;
 783       while (i) {
 784         GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
 785         if (GTS_OBJECT (n1)->reserved &&
 786             gts_containee_is_contained (GTS_CONTAINEE (n1),
 787                                         GTS_CONTAINER (bg->g))) {
 788           GtsEHeap * h =
 789             gts_containee_is_contained (GTS_CONTAINEE (n1),
 790                                         GTS_CONTAINER (bg->g1)) ? h1 : h2;
 791           gts_eheap_remove (h, GTS_OBJECT (n1)->reserved);
 792           GTS_OBJECT (n1)->reserved = gts_eheap_insert (h, n1);
 793         }
 794         i = i->next;
 795       }
 796     }
 797   } while (n && nm < mmax);
 798
 799   gts_eheap_foreach (h1, (GFunc) gts_object_reset_reserved, NULL);
 800   gts_eheap_foreach (h2, (GFunc) gts_object_reset_reserved, NULL);
 801   gts_eheap_destroy (h1);
 802   gts_eheap_destroy (h2);
 803
 804   /* undo last nm moves */
 805   for (i = 0; i < nm; i++) {
 806     GtsGNode * n = moves[i];
 807     GtsGraph * g1 =
 808       gts_containee_is_contained (GTS_CONTAINEE (n),
 809                                   GTS_CONTAINER (bg->g1)) ? bg->g1 : bg->g2;
 810     GtsGraph * g2 = g1 == bg->g1 ? bg->g2 : bg->g1;
 811
 812     gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
 813     gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 814   }
 815   g_free (moves);
 816
 817   return bestcost;
 818 }
 819
 820 static void build_bheap (GtsGNode * n, GtsGNode * n1, GtsEHeap * heap)
 821 {
 822   GTS_OBJECT (n)->reserved = gts_eheap_insert (heap, n);
 823 }
 824
 825 static void update_neighbors (GtsGNode * n, GtsGraphBisection * bg,
 826                               GtsEHeap * h1, GtsEHeap * h2)
 827 {
 828   GSList * i;
 829
 830   i = GTS_SLIST_CONTAINER (n)->items;
 831   while (i) {
 832     GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
 833     if (gts_containee_is_contained (GTS_CONTAINEE (n1),
 834                                     GTS_CONTAINER (bg->g))) {
 835       GtsEHeap * h;
 836       GtsGraph /* * g1,*/ * g2;
 837       GHashTable * bg1;
 838
 839       if (gts_containee_is_contained (GTS_CONTAINEE (n1),
 840                                       GTS_CONTAINER (bg->g1))) {
 841         h = h1;
 842         //g1 = bg->g1;
 843         g2 = bg->g2;
 844         bg1 = bg->bg1;
 845       }
 846       else {
 847         h = h2;
 848         //g1 = bg->g2;
 849         g2 = bg->g1;
 850         bg1 = bg->bg2;
 851       }
 852       g_hash_table_remove (bg1, n1);
 853       if (h && GTS_OBJECT (n1)->reserved && GTS_OBJECT (n1)->reserved != n1) {
 854         gts_eheap_remove (h, GTS_OBJECT (n1)->reserved);
 855         GTS_OBJECT (n1)->reserved = NULL;
 856       }
 857       if (gts_gnode_degree (n1, g2)) {
 858         g_hash_table_insert (bg1, n1, n1);
 859         if (h && GTS_OBJECT (n1)->reserved != n1)
 860           GTS_OBJECT (n1)->reserved = gts_eheap_insert (h, n1);
 861       }
 862     }
 863     i = i->next;
 864   }
 865 }
 866
 867 /**
 868  * gts_graph_bisection_bkl_refine:
 869  * @bg: a #GtsGraphBisection.
 870  * @mmax: the maximum number of unsuccessful successive moves.
 871  * @imbalance: the maximum relative imbalance allowed between the
 872  * weights of both halves of the partition.
 873  *
 874  * An implementation of the simplified boundary Kernighan-Lin
 875  * algorithm for graph bisection refinement as described in Karypis
 876  * and Kumar (1997).
 877  *
 878  * The algorithm stops if @mmax consecutive modes do not lead to a
 879  * decrease in the number of edges cut. This last @mmax moves are
 880  * undone.
 881  *
 882  * Returns: the decrease in the weight of the edges cut by the bisection.
 883  */
 884 gdouble gts_graph_bisection_bkl_refine (GtsGraphBisection * bg,
 885                                         guint mmax,
 886                                         gfloat imbalance)
 887 {
 888   GtsEHeap * h1, * h2;
 889   GtsGNode * n;
 890   guint nm = 0, i;
 891   GtsGNode ** moves;
 892   gdouble bestcost = 0., totalcost = 0., best_balance;
 893   gboolean balanced = FALSE;
 894
 895   g_return_val_if_fail (bg != NULL, 0.);
 896   g_return_val_if_fail (mmax > 0, 0.);
 897   g_return_val_if_fail (imbalance >= 0. && imbalance <= 1., 0.);
 898
 899   h1 = gts_eheap_new ((GtsKeyFunc) node_move_cost1, bg);
 900   gts_eheap_freeze (h1);
 901   g_hash_table_foreach (bg->bg1, (GHFunc) build_bheap, h1);
 902   gts_eheap_thaw (h1);
 903
 904   h2 = gts_eheap_new ((GtsKeyFunc) node_move_cost2, bg);
 905   gts_eheap_freeze (h2);
 906   g_hash_table_foreach (bg->bg2, (GHFunc) build_bheap, h2);
 907   gts_eheap_thaw (h2);
 908
 909   moves = g_malloc (sizeof (GtsGNode *)*mmax);
 910   imbalance *= gts_graph_weight (bg->g);
 911   best_balance = fabs (gts_graph_weight (bg->g1) - gts_graph_weight (bg->g2));
 912   if (best_balance <= imbalance)
 913     balanced = TRUE;
 914
 915   do {
 916     GtsGraph * g1, * g2;
 917     GHashTable * bg1, * bg2;
 918     gdouble cost;
 919
 920     if (gts_graph_weight (bg->g1) > gts_graph_weight (bg->g2)) {
 921       n = gts_eheap_remove_top (h1, &cost);
 922       g1 = bg->g1;
 923       g2 = bg->g2;
 924       bg1 = bg->bg1;
 925       bg2 = bg->bg2;
 926     }
 927     else {
 928       n = gts_eheap_remove_top (h2, &cost);
 929       g1 = bg->g2;
 930       g2 = bg->g1;
 931       bg1 = bg->bg2;
 932       bg2 = bg->bg1;
 933     }
 934     if (n) {
 935       gdouble balance;
 936
 937       GTS_OBJECT (n)->reserved = n;
 938       gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
 939       gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 940       g_hash_table_remove (bg1, n);
 941       if (gts_gnode_degree (n, g1))
 942         g_hash_table_insert (bg2, n, n);
 943
 944       update_neighbors (n, bg, h1, h2);
 945
 946       totalcost += cost;
 947       balance = fabs (gts_graph_weight (g1) - gts_graph_weight (g2));
 948
 949       if (!balanced && balance <= imbalance) {
 950         bestcost = totalcost;
 951         best_balance = balance;
 952         balanced = TRUE;
 953         nm = 0;
 954       }
 955       else if (totalcost < bestcost &&
 956                (balance < best_balance || balance <= imbalance)) {
 957         bestcost = totalcost;
 958         best_balance = balance;
 959         nm = 0;
 960       }
 961       else if (totalcost == bestcost && balance < best_balance) {
 962         best_balance = balance;
 963         nm = 0;
 964       }
 965       else
 966         moves[nm++] = n;
 967     }
 968   } while (n && nm < mmax);
 969
 970   gts_container_foreach (GTS_CONTAINER (bg->g),
 971                          (GtsFunc) gts_object_reset_reserved, NULL);
 972   gts_eheap_destroy (h1);
 973   gts_eheap_destroy (h2);
 974
 975   /* undo last nm moves */
 976   for (i = 0; i < nm; i++) {
 977     GtsGNode * n = moves[i];
 978     GtsGraph * g1, * g2;
 979     GHashTable * bg1, * bg2;
 980
 981     if (gts_containee_is_contained (GTS_CONTAINEE (n),
 982                                     GTS_CONTAINER (bg->g1))) {
 983       g1 = bg->g1;
 984       g2 = bg->g2;
 985       bg1 = bg->bg1;
 986       bg2 = bg->bg2;
 987     }
 988     else {
 989       g1 = bg->g2;
 990       g2 = bg->g1;
 991       bg1 = bg->bg2;
 992       bg2 = bg->bg1;
 993     }
 994
 995     gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
 996     gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
 997     g_hash_table_remove (bg1, n);
 998     if (gts_gnode_degree (n, g1))
 999       g_hash_table_insert (bg2, n, n);
1000
1001     update_neighbors (n, bg, NULL, NULL);
1002   }
1003   g_free (moves);
1004
1005   return bestcost;
1006 }
1007
1008 /* Multilevel partitioning */
1009
1010 static void bisection_children (GtsGNodeSplit * ns, GtsGraphBisection * bg)
1011 {
1012   GtsGraph * g, * g1;
1013   GHashTable * bbg;
1014   GtsGNode * n1 = GTS_GNODE_SPLIT_N1 (ns);
1015   GtsGNode * n2 = GTS_GNODE_SPLIT_N2 (ns);
1016
1017   if (gts_containee_is_contained (GTS_CONTAINEE (ns->n),
1018                                   GTS_CONTAINER (bg->g1))) {
1019     g = bg->g1;
1020     g1 = bg->g2;
1021     bbg = bg->bg1;
1022   }
1023   else {
1024     g = bg->g2;
1025     g1 = bg->g1;
1026     bbg = bg->bg2;
1027   }
1028
1029   gts_allow_floating_gnodes = TRUE;
1030   gts_container_remove (GTS_CONTAINER (g), GTS_CONTAINEE (ns->n));
1031   gts_allow_floating_gnodes = FALSE;
1032   gts_container_add (GTS_CONTAINER (g), GTS_CONTAINEE (n1));
1033   gts_container_add (GTS_CONTAINER (g), GTS_CONTAINEE (n2));
1034
1035   if (g_hash_table_lookup (bbg, ns->n)) {
1036     g_hash_table_remove (bbg, ns->n);
1037     if (gts_gnode_degree (n1, g1) > 0)
1038       g_hash_table_insert (bbg, n1, n1);
1039     if (gts_gnode_degree (n2, g1) > 0)
1040       g_hash_table_insert (bbg, n2, n2);
1041   }
1042 }
1043
1044 /**
1045  * gts_graph_bisection_new:
1046  * @wg: a #GtsWGraph.
1047  * @ntry: the number of tries for the graph growing algorithm.
1048  * @mmax: the number of unsucessful moves for the refinement algorithm.
1049  * @nmin: the minimum number of nodes of the coarsest graph.
1050  * @imbalance: the maximum relative imbalance allowed between the
1051  * weights of both halves of the partition.
1052  *
1053  * An implementation of a multilevel bisection algorithm as presented
1054  * in Karypis and Kumar (1997). A multilevel hierarchy of graphs is
1055  * created using the #GtsPGraph object. The bisection of the coarsest
1056  * graph is created using the gts_graph_ggg_bisection() function. The
1057  * graph is then uncoarsened using gts_pgraph_down() and at each level
1058  * the bisection is refined using gts_graph_bisection_bkl_refine().
1059  *
1060  * Returns: a new #GtsGraphBisection of @wg.
1061  */
1062 GtsGraphBisection * gts_graph_bisection_new (GtsWGraph * wg,
1063                                              guint ntry,
1064                                              guint mmax,
1065                                              guint nmin,
1066                                              gfloat imbalance)
1067 {
1068   GtsGraph * g;
1069   GtsPGraph * pg;
1070   GtsGraphBisection * bg;
1071   gdouble cost;
1072
1073   g_return_val_if_fail (wg != NULL, NULL);
1074
1075   g = GTS_GRAPH (wg);
1076   pg = gts_pgraph_new (gts_pgraph_class (), g,
1077                        gts_gnode_split_class (),
1078                        gts_wgnode_class (),
1079                        gts_wgedge_class (),
1080                        nmin);
1081
1082   bg = gts_graph_ggg_bisection (g, ntry);
1083 #ifdef DEBUG
1084   fprintf (stderr, "before size: %5d weight: %5g cuts: %5d cweight: %5g\n",
1085            gts_container_size (GTS_CONTAINER (bg->g1)),
1086            gts_graph_weight (bg->g1),
1087            gts_graph_edges_cut (bg->g1),
1088            gts_graph_edges_cut_weight (bg->g1));
1089   g_assert (gts_graph_bisection_check (bg));
1090 #endif
1091   while ((cost = gts_graph_bisection_bkl_refine (bg, mmax, imbalance))) {
1092 #ifdef DEBUG
1093     fprintf (stderr, "  cost: %g\n", cost);
1094     g_assert (gts_graph_bisection_check (bg));
1095 #endif
1096   }
1097 #ifdef DEBUG
1098   fprintf (stderr, "after  size: %5d weight: %5g cuts: %5d cweight: %5g\n",
1099            gts_container_size (GTS_CONTAINER (bg->g1)),
1100            gts_graph_weight (bg->g1),
1101            gts_graph_edges_cut (bg->g1),
1102            gts_graph_edges_cut_weight (bg->g1));
1103 #endif
1104   while (gts_pgraph_down (pg, (GtsFunc) bisection_children, bg)) {
1105 #ifdef DEBUG
1106     fprintf (stderr, "before size: %5d weight: %5g cuts: %5d cweight: %5g\n",
1107              gts_container_size (GTS_CONTAINER (bg->g1)),
1108              gts_graph_weight (bg->g1),
1109              gts_graph_edges_cut (bg->g1),
1110              gts_graph_edges_cut_weight (bg->g1));
1111 #endif
1112     while ((cost = gts_graph_bisection_bkl_refine (bg, mmax, imbalance))) {
1113 #ifdef DEBUG
1114       fprintf (stderr, "  cost: %g\n", cost);
1115       g_assert (gts_graph_bisection_check (bg));
1116 #endif
1117     }
1118 #ifdef DEBUG
1119     fprintf (stderr, "after  size: %5d weight: %5g cuts: %5d cweight: %5g\n",
1120              gts_container_size (GTS_CONTAINER (bg->g1)),
1121              gts_graph_weight (bg->g1),
1122              gts_graph_edges_cut (bg->g1),
1123              gts_graph_edges_cut_weight (bg->g1));
1124 #endif
1125   }
1126   gts_object_destroy (GTS_OBJECT (pg));
1127
1128   return bg;
1129 }
1130
1131 /**
1132  * gts_graph_bisection_destroy:
1133  * @bg: a #GtsGraphBisection.
1134  * @destroy_graphs: controls graph destruction.
1135  *
1136  * Frees all the memory allocated for @bg. If @destroy_graphs is %TRUE
1137  * the graphs created by @bg are destroyed.
1138  */
1139 void gts_graph_bisection_destroy (GtsGraphBisection * bg,
1140                                   gboolean destroy_graphs)
1141 {
1142   g_return_if_fail (bg != NULL);
1143
1144   g_hash_table_destroy (bg->bg1);
1145   g_hash_table_destroy (bg->bg2);
1146
1147   if (destroy_graphs) {
1148     gts_object_destroy (GTS_OBJECT (bg->g1));
1149     gts_object_destroy (GTS_OBJECT (bg->g2));
1150   }
1151
1152   g_free (bg);
1153 }
1154
1155 static void recursive_bisection (GtsWGraph * wg,
1156                                  guint np,
1157                                  guint ntry,
1158                                  guint mmax,
1159                                  guint nmin,
1160                                  gfloat imbalance,
1161                                  GSList ** list)
1162 {
1163   if (np == 0)
1164     *list = g_slist_prepend (*list, wg);
1165   else {
1166     GtsGraphBisection * bg =
1167       gts_graph_bisection_new (wg, ntry, mmax, nmin, imbalance);
1168     GtsGraph * g1 = bg->g1;
1169     GtsGraph * g2 = bg->g2;
1170
1171     gts_object_destroy (GTS_OBJECT (wg));
1172     gts_graph_bisection_destroy (bg, FALSE);
1173     recursive_bisection (GTS_WGRAPH (g1), np - 1, ntry, mmax, nmin, imbalance,
1174                          list);
1175     recursive_bisection (GTS_WGRAPH (g2), np - 1, ntry, mmax, nmin, imbalance,
1176                          list);
1177   }
1178 }
1179
1180 /**
1181  * gts_graph_recursive_bisection:
1182  * @wg: a #GtsWGraph.
1183  * @n: the number of bisection levels.
1184  * @ntry: the number of tries for the graph growing algorithm.
1185  * @mmax: the number of unsucessful moves for the refinement algorithm.
1186  * @nmin: the minimum number of nodes of the coarsest graph.
1187  * @imbalance: the maximum relative imbalance allowed between the
1188  * weights of both halves of the partition.
1189  *
1190  * Calls gts_graph_bisection_new() recursively in order to obtain a
1191  * 2^@n partition of @wg.
1192  *
1193  * Returns: a list of 2^@n new #GtsGraph representing the partition.
1194  */
1195 GSList * gts_graph_recursive_bisection (GtsWGraph * wg,
1196                                         guint n,
1197                                         guint ntry,
1198                                         guint mmax,
1199                                         guint nmin,
1200                                         gfloat imbalance)
1201 {
1202   GtsGraphBisection * bg;
1203   GtsGraph * g1, * g2;
1204   GSList * list = NULL;
1205
1206   g_return_val_if_fail (wg != NULL, NULL);
1207   g_return_val_if_fail (n > 0, NULL);
1208
1209   bg = gts_graph_bisection_new (wg, ntry, mmax, nmin, imbalance);
1210   g1 = bg->g1;
1211   g2 = bg->g2;
1212   gts_graph_bisection_destroy (bg, FALSE);
1213   recursive_bisection (GTS_WGRAPH (g1), n - 1, ntry, mmax, nmin, imbalance,
1214                        &list);
1215   recursive_bisection (GTS_WGRAPH (g2), n - 1, ntry, mmax, nmin, imbalance,
1216                        &list);
1217
1218   return list;
1219 }