partition.c

Go to the documentation of this file.

/* GTS - Library for the manipulation of triangulated surfaces
 * Copyright (C) 1999 Stéphane Popinet
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#include <math.h>

#include "gts.h"

/* #define DEBUG */

/* Graph partition */

guint gts_graph_partition_edges_cut (GSList * partition)
{
  guint cuts = 0;

  while (partition) {
    cuts += gts_graph_edges_cut (partition->data);
    partition = partition->next;
  }

  return cuts/2;
}

gfloat gts_graph_partition_edges_cut_weight (GSList * partition)
{
  gfloat weight = 0.;

  while (partition) {
    weight += gts_graph_edges_cut_weight (partition->data);
    partition = partition->next;
  }

  return weight/2.;
}

void gts_graph_partition_print_stats (GSList * partition,
                                      FILE * fp)
{
  GtsRange weight;
  GSList * i;

  g_return_if_fail (partition != NULL);
  g_return_if_fail (fp != NULL);

  gts_range_init (&weight);
  i = partition;
  while (i) {
    gts_range_add_value (&weight, gts_graph_weight (i->data));
    i = i->next;
  }
  gts_range_update (&weight);

  fprintf (fp, 
           "# parts: %d\n"
           "#   edge cuts: %5d edge cuts weight: %5g\n"
           "#   weight: ",
           g_slist_length (partition),
           gts_graph_partition_edges_cut (partition),
           gts_graph_partition_edges_cut_weight (partition));
  gts_range_print (&weight, fp);
  fputc ('\n', fp);
}

gfloat gts_graph_partition_balance (GSList * partition)
{
  gfloat wmin = G_MAXFLOAT;
  gfloat wmax = - G_MAXFLOAT;

  g_return_val_if_fail (partition != NULL, 0.);

  while (partition) {
    gfloat weight = gts_graph_weight (partition->data);
    if (weight < wmin)
      wmin = weight;
    if (weight > wmax)
      wmax = weight;
    partition = partition->next;
  }
  return wmax - wmin;
}

GSList * gts_graph_partition_clone (GSList * partition)
{
  GSList * cparts = NULL;

  while (partition) {
    cparts = 
      g_slist_prepend (cparts, 
                       gts_object_clone (GTS_OBJECT (partition->data)));
    partition = partition->next;
  }
  return cparts;
}

void gts_graph_partition_destroy (GSList * partition)
{
  GSList * i = partition;

  while (i) {
    gts_object_destroy (GTS_OBJECT (i->data));
    i = i->next;
  }
  g_slist_free (partition);
}

static void find_smallest_degree (GtsGNode * n, gpointer * data)
{
  GtsGNode ** nmin = data[0];
  GtsGraph * g = data[1];
  guint * min = data[2];
  guint degree = gts_gnode_degree (n, g);

  if (degree < *min) {
    *min = degree;
    *nmin = n;
  }
}

static gint graph_comp_weight (GtsGraph * g1, GtsGraph * g2)
{
  if (gts_graph_weight (g1) > gts_graph_weight (g2))
    return 1;
  return -1;
}

static void partition_update (GSList * list, GtsGraph * g)
{
  GSList * i;
  GtsGraph * g1;
  GtsHeap * size_heap;
  gboolean reinit = TRUE;

  /* initialize traversals */
  i = list;
  while (i) {
    GtsGNode * seed = GTS_OBJECT (i->data)->reserved;
    GTS_OBJECT (seed)->reserved = 
      gts_graph_traverse_new (g, seed, GTS_BREADTH_FIRST, reinit);
    reinit = FALSE;
    i = i->next;
  }
  
  size_heap = gts_heap_new ((GCompareFunc) graph_comp_weight);
  i = list;
  while (i) {
    gts_heap_insert (size_heap, i->data);
    i = i->next;
  }
  while ((g1 = gts_heap_remove_top (size_heap))) {
    GtsGraphTraverse * t = GTS_OBJECT (GTS_OBJECT (g1)->reserved)->reserved;
    GtsGNode * n = gts_graph_traverse_next (t);
    if (n) {
      gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
      gts_heap_insert (size_heap, g1);
    }
  }
  gts_heap_destroy (size_heap);

  /* destroy traversals */
  i = list;
  while (i) {
    GtsGNode * seed = GTS_OBJECT (i->data)->reserved;
    gts_graph_traverse_destroy (GTS_OBJECT (seed)->reserved);
    GTS_OBJECT (seed)->reserved = NULL;
    i = i->next;
  }
}

static void better_seed (GtsGNode * n, gpointer * data)
{
  guint * sum = data[0];
  GtsGNode ** seed = data[1];
  GtsGraph * g = data[2];
  guint sum1 = gts_graph_distance_sum (g, n);
  
  if (sum1 < *sum) {
    *sum = sum1;
    *seed = n;
  }
}

static GtsGNode * graph_new_seed (GtsGraph * g, GtsGNode * seed)
{
  guint sum = gts_graph_distance_sum (g, seed);
  gpointer data[3];
  GtsGNode * new_seed = seed;

  data[0] = &sum;
  data[1] = &new_seed;
  data[2] = g;
  gts_gnode_foreach_neighbor (seed, g, (GtsFunc) better_seed, data);

  return new_seed;
}

GSList * gts_graph_bubble_partition (GtsGraph * g, 
                                     guint np, 
                                     guint niter,
                                     GtsFunc step_info,
                                     gpointer data)
{
  GSList * list = NULL, * seeds = NULL;
  GtsGNode * seed = NULL;
  guint min = G_MAXINT/2 - 1;
  gpointer info[3];
  GtsGraph * g1;
  gboolean changed = TRUE;

  g_return_val_if_fail (g != NULL, NULL);
  g_return_val_if_fail (np > 0, NULL);

  info[0] = &seed;
  info[1] = g;
  info[2] = &min;
  gts_container_foreach (GTS_CONTAINER (g), 
                         (GtsFunc) find_smallest_degree,
                         info);
  if (seed == NULL)
    return NULL;

  g1 = GTS_GRAPH (gts_object_new (GTS_OBJECT (g)->klass));
  gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
  list = g_slist_prepend (list, g1);
  GTS_OBJECT (g1)->reserved = seed;
  seeds = g_slist_prepend (seeds, seed);

  while (--np && seed)
    if ((seed = gts_graph_farthest (g, seeds))) {
      g1 = GTS_GRAPH (gts_object_new (GTS_OBJECT (g)->klass));
      gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
      list = g_slist_prepend (list, g1);
      GTS_OBJECT (g1)->reserved = seed;
      seeds = g_slist_prepend (seeds, seed);
    }
  g_slist_free (seeds);
  
  partition_update (list, g);

  while (changed && niter--) {
    GSList * i;

    changed = FALSE;
    i = list;
    while (i) {
      GtsGraph * g1 = i->data;
      GtsGNode * seed = GTS_OBJECT (g1)->reserved;
      GtsGNode * new_seed = graph_new_seed (g1, seed);
      if (new_seed != seed) {
        changed = TRUE;
        GTS_OBJECT (g1)->reserved = new_seed;
      }
      i = i->next;
    }

    if (changed) {
      i = list;
      while (i) {
        GtsGraph * g1 = i->data;
        GtsGNode * seed = GTS_OBJECT (g1)->reserved;

        gts_object_destroy (GTS_OBJECT (g1));
        i->data = g1 = GTS_GRAPH (gts_object_new (GTS_OBJECT (g)->klass));
        gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
        GTS_OBJECT (g1)->reserved = seed;
        i = i->next;
      }
      partition_update (list, g);
      if (step_info)
        (* step_info) (list, data);
    }
  }
  g_slist_foreach (list, (GFunc) gts_object_reset_reserved, NULL);
  return list;
}

/* Graph bisection */

static gdouble node_cost (GtsGNode * n, gpointer * data)
{
  GtsGraph * g = data[0];
  GtsGraph * g1 = data[1];
  GSList * i = GTS_SLIST_CONTAINER (n)->items;
  gdouble cost = 0.;

  while (i) {
    GtsGEdge * e = i->data;
    GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, e);

    if (gts_containee_is_contained (GTS_CONTAINEE (n1), GTS_CONTAINER (g))) {
      if (gts_containee_is_contained (GTS_CONTAINEE (n1), GTS_CONTAINER (g1)))
        cost -= gts_gedge_weight (e);
      else 
        cost += gts_gedge_weight (e);
    }
    i = i->next;
  }

  return cost;
}

static void add_neighbor (GtsGNode * n, GtsEHeap * heap)
{
  if (GTS_OBJECT (n)->reserved == n)
    return;
  if (GTS_OBJECT (n)->reserved)
    gts_eheap_remove (heap, GTS_OBJECT (n)->reserved);
  GTS_OBJECT (n)->reserved = gts_eheap_insert (heap, n);
}

static void add_unused (GtsGNode * n, GtsGraph * g2)
{
  if (GTS_OBJECT (n)->reserved)
    GTS_OBJECT (n)->reserved = NULL;
  else
    gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
}

static gdouble degree_cost (GtsGNode * n, GtsGraph * g)
{
  return gts_gnode_degree (n, g); 
}

static void add_seed (GtsGNode * n, GtsEHeap * heap)
{
  gts_eheap_insert (heap, n);
}

static void boundary_node1 (GtsGNode * n, GtsGraphBisection * bg)
{
  GSList * i = GTS_SLIST_CONTAINER (n)->items;

  while (i) {
    GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
    if (gts_containee_is_contained (GTS_CONTAINEE (n1), 
                                    GTS_CONTAINER (bg->g2))) {
      g_hash_table_insert (bg->bg1, n, n1);
      return;
    }
    i = i->next;
  }
}

static void boundary_node2 (GtsGNode * n, GtsGraphBisection * bg)
{
  GSList * i = GTS_SLIST_CONTAINER (n)->items;

  while (i) {
    GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
    if (gts_containee_is_contained (GTS_CONTAINEE (n1), 
                                    GTS_CONTAINER (bg->g1))) {
      g_hash_table_insert (bg->bg2, n, n1);
      return;
    }
    i = i->next;
  }
}

static void check_bg (GtsGNode * n, gpointer * data)
{
  GHashTable * bg = data[0];
  GtsGraph * g = data[1];
  gboolean * ok = data[2];
  guint * nb = data[3];
  guint nn = gts_gnode_degree (n, g);

  if (nn > 0)
    (*nb)++;
  if ((nn > 0 && !g_hash_table_lookup (bg, n)) ||
      (nn == 0 && g_hash_table_lookup (bg, n))) {
    g_warning ("nn: %d lookup: %p\n",
               nn, g_hash_table_lookup (bg, n));
    *ok = FALSE;
  }
}

gboolean gts_graph_bisection_check (GtsGraphBisection * bg)
{
  gboolean ok = TRUE;
  guint nb;
  gpointer data[4];

  g_return_val_if_fail (bg != NULL, FALSE);

  nb = 0;
  data[0] = bg->bg1;
  data[1] = bg->g2;
  data[2] = &ok;
  data[3] = &nb;
  gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) check_bg, data);
  g_return_val_if_fail (g_hash_table_size (bg->bg1) == nb, FALSE);

  nb = 0;
  data[0] = bg->bg2;
  data[1] = bg->g1;
  gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) check_bg, data);
  g_return_val_if_fail (g_hash_table_size (bg->bg2) == nb, FALSE);

  return ok;
}

GtsGraphBisection * gts_graph_ggg_bisection (GtsGraph * g, guint ntry)
{
  gfloat size, bestcost = G_MAXFLOAT, smin;
  GtsGraph * bestg1 = NULL, * bestg2 = NULL;
  gboolean balanced = FALSE;
  GtsEHeap * degree_heap;
  GtsGNode * seed;
  GtsGraphBisection * bg;

  g_return_val_if_fail (g != NULL, NULL);

  bg = g_malloc (sizeof (GtsGraphBisection));
  bg->g = g;

  size = gts_graph_weight (g)/2.;
  smin = 0.9*size;

  degree_heap = gts_eheap_new ((GtsKeyFunc) degree_cost, g);
  gts_eheap_freeze (degree_heap);
  gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_seed, degree_heap);
  gts_eheap_thaw (degree_heap);

  while (ntry && ((seed = gts_eheap_remove_top (degree_heap, NULL)))) {
    GtsGraph * g1, * g2;
    GtsGNode * n;
    gdouble cost;
    gpointer data[2];
    GtsEHeap * heap;
  
    g1 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
                        g->node_class, g->edge_class);
    g2 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
                        g->node_class, g->edge_class);
    
    data[0] = g;
    data[1] = g1;
    heap = gts_eheap_new ((GtsKeyFunc) node_cost, data);

    gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (seed));
    GTS_OBJECT (seed)->reserved = seed;
    gts_gnode_foreach_neighbor (seed, g, (GtsFunc) add_neighbor, heap);

    while ((n = gts_eheap_remove_top (heap, &cost)))
      if (gts_graph_weight (g1) + gts_gnode_weight (n) <= size) {
        gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
        GTS_OBJECT (n)->reserved = n;
        gts_gnode_foreach_neighbor (n, g, (GtsFunc) add_neighbor, heap);
      }
      else
        GTS_OBJECT (n)->reserved = NULL;
    gts_eheap_destroy (heap);
    
    gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_unused, g2);

    cost = gts_graph_edges_cut_weight (g1);
    if (!bestg1 || 
        (!balanced && gts_graph_weight (g1) >= smin) ||
        (cost < bestcost && gts_graph_weight (g1) >= smin)) {
      if (bestg1)
        bestcost = cost;
      if (bestg1)
        gts_object_destroy (GTS_OBJECT (bestg1));
      if (bestg2)
        gts_object_destroy (GTS_OBJECT (bestg2));
      bestg1 = g1;
      bestg2 = g2;
      if (gts_graph_weight (g1) >= smin)
        balanced = TRUE;
    }
    else {
      gts_object_destroy (GTS_OBJECT (g1));
      gts_object_destroy (GTS_OBJECT (g2));
    }

    ntry--;
  }
  gts_eheap_destroy (degree_heap);

#ifdef DEBUG
  fprintf (stderr, "bestcost: %5g g1: %5g|%5d g2: %5g|%5d\n",
           bestcost, 
           gts_graph_weight (bestg1), 
           gts_container_size (GTS_CONTAINER (bestg1)),
           gts_graph_weight (bestg2), 
           gts_container_size (GTS_CONTAINER (bestg2)));
#endif

  g_assert (bestg1 != NULL);
  bg->g1 = bestg1;
  g_assert (bestg2 != NULL);
  bg->g2 = bestg2;
  
  /* boundary nodes */
  bg->bg1 = g_hash_table_new (NULL, NULL);
  gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) boundary_node1, bg);
  bg->bg2 = g_hash_table_new (NULL, NULL);
  gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) boundary_node2, bg);

  return bg;
}

GtsGraphBisection * gts_graph_bfgg_bisection (GtsGraph * g, guint ntry)
{
  gfloat size, bestcost = G_MAXFLOAT, smin;
  GtsGraph * bestg1 = NULL, * bestg2 = NULL;
  GtsEHeap * degree_heap;
  GtsGNode * seed;
  GtsGraphBisection * bg;

  g_return_val_if_fail (g != NULL, NULL);

  bg = g_malloc (sizeof (GtsGraphBisection));
  bg->g = g;

  size = gts_graph_weight (g)/2.;
  smin = 0.9*size;

  degree_heap = gts_eheap_new ((GtsKeyFunc) degree_cost, g);
  gts_eheap_freeze (degree_heap);
  gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_seed, degree_heap);
  gts_eheap_thaw (degree_heap);

  while (ntry && ((seed = gts_eheap_remove_top (degree_heap, NULL)))) {
    GtsGraph * g1, * g2;
    GtsGNode * n;
    gdouble cost;
    GtsGraphTraverse * t = gts_graph_traverse_new (g, seed, 
                                                   GTS_BREADTH_FIRST, TRUE);
    
    g1 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
                        g->node_class, g->edge_class);
    g2 = gts_graph_new (GTS_GRAPH_CLASS (GTS_OBJECT (g)->klass),
                        g->node_class, g->edge_class);

    while ((n = gts_graph_traverse_next (t)))
      if (gts_graph_weight (g1) + gts_gnode_weight (n) <= size) {
        gts_container_add (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
        GTS_OBJECT (n)->reserved = n;
      }
    gts_graph_traverse_destroy (t);
    
    gts_container_foreach (GTS_CONTAINER (g), (GtsFunc) add_unused, g2);

    cost = gts_graph_edges_cut_weight (g1);
    if (!bestg1 || (cost < bestcost && gts_graph_weight (g1) >= smin)) {
      if (bestg1)
        bestcost = cost;
      if (bestg1)
        gts_object_destroy (GTS_OBJECT (bestg1));
      if (bestg2)
        gts_object_destroy (GTS_OBJECT (bestg2));
      bestg1 = g1;
      bestg2 = g2;
    }
    else {
      gts_object_destroy (GTS_OBJECT (g1));
      gts_object_destroy (GTS_OBJECT (g2));
    }

    ntry--;
  }
  gts_eheap_destroy (degree_heap);

#ifdef DEBUG
  fprintf (stderr, "bestcost: %5g g1: %5g|%5d g2: %5g|%5d\n",
           bestcost, 
           gts_graph_weight (bestg1), 
           gts_container_size (GTS_CONTAINER (bestg1)),
           gts_graph_weight (bestg2), 
           gts_container_size (GTS_CONTAINER (bestg2)));
#endif

  bg->g1 = bestg1;
  bg->g2 = bestg2;
  
  /* boundary nodes */
  bg->bg1 = g_hash_table_new (NULL, NULL);
  gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) boundary_node1, bg);
  bg->bg2 = g_hash_table_new (NULL, NULL);
  gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) boundary_node2, bg);

  return bg;
}

static gdouble node_move_cost1 (GtsGNode * n, GtsGraphBisection * bg)
{
  return gts_gnode_move_cost (n, bg->g1, bg->g2);
}

static gdouble node_move_cost2 (GtsGNode * n, GtsGraphBisection * bg)
{
  return gts_gnode_move_cost (n, bg->g2, bg->g1);
}

static void build_heap (GtsGNode * n, GtsEHeap * heap)
{
  GTS_OBJECT (n)->reserved = gts_eheap_insert (heap, n);
}

gdouble gts_graph_bisection_kl_refine (GtsGraphBisection * bg,
                                       guint mmax)
{
  GtsEHeap * h1, * h2;
  GtsGNode * n;
  guint nm = 0, i;
  GtsGNode ** moves;
  gdouble bestcost = 0., totalcost = 0., best_balance;

  g_return_val_if_fail (bg != NULL, 0.);
  g_return_val_if_fail (mmax > 0, 0.);

  h1 = gts_eheap_new ((GtsKeyFunc) node_move_cost1, bg);
  gts_eheap_freeze (h1);
  gts_container_foreach (GTS_CONTAINER (bg->g1), (GtsFunc) build_heap, h1);
  gts_eheap_thaw (h1);

  h2 = gts_eheap_new ((GtsKeyFunc) node_move_cost2, bg);
  gts_eheap_freeze (h2);
  gts_container_foreach (GTS_CONTAINER (bg->g2), (GtsFunc) build_heap, h2);
  gts_eheap_thaw (h2);

  moves = g_malloc (sizeof (GtsGNode *)*mmax);
  best_balance = fabs (gts_graph_weight (bg->g1) - gts_graph_weight (bg->g2));

  do {
    GtsGraph * g1, * g2;
    gdouble cost;

    if (gts_graph_weight (bg->g1) > gts_graph_weight (bg->g2)) {
      n = gts_eheap_remove_top (h1, &cost);
      g1 = bg->g1;
      g2 = bg->g2;
    }
    else {
      n = gts_eheap_remove_top (h2, &cost);
      g1 = bg->g2;
      g2 = bg->g1;
    }
    if (n) {
      GSList * i;

      GTS_OBJECT (n)->reserved = NULL;
      gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
      gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));

      totalcost += cost;
      if (totalcost < bestcost) {
        bestcost = totalcost;
        nm = 0;
      }
      else if (totalcost == bestcost) {
        gdouble balance = fabs (gts_graph_weight (g1) - gts_graph_weight (g2));

        if (balance < best_balance) {
          best_balance = balance;
          nm = 0;
        }
      }        
      else
        moves[nm++] = n;

      i = GTS_SLIST_CONTAINER (n)->items;
      while (i) {
        GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
        if (GTS_OBJECT (n1)->reserved && 
            gts_containee_is_contained (GTS_CONTAINEE (n1), 
                                        GTS_CONTAINER (bg->g))) {
          GtsEHeap * h = 
            gts_containee_is_contained (GTS_CONTAINEE (n1), 
                                        GTS_CONTAINER (bg->g1)) ? h1 : h2;
          gts_eheap_remove (h, GTS_OBJECT (n1)->reserved);
          GTS_OBJECT (n1)->reserved = gts_eheap_insert (h, n1);
        }
        i = i->next;
      }
    }
  } while (n && nm < mmax);

  gts_eheap_foreach (h1, (GFunc) gts_object_reset_reserved, NULL);
  gts_eheap_foreach (h2, (GFunc) gts_object_reset_reserved, NULL);
  gts_eheap_destroy (h1);
  gts_eheap_destroy (h2);

  /* undo last nm moves */
  for (i = 0; i < nm; i++) {
    GtsGNode * n = moves[i];
    GtsGraph * g1 = 
      gts_containee_is_contained (GTS_CONTAINEE (n),
                                  GTS_CONTAINER (bg->g1)) ? bg->g1 : bg->g2;
    GtsGraph * g2 = g1 == bg->g1 ? bg->g2 : bg->g1;
    
    gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
    gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
  }
  g_free (moves);

  return bestcost;
}

static void build_bheap (GtsGNode * n, GtsGNode * n1, GtsEHeap * heap)
{
  GTS_OBJECT (n)->reserved = gts_eheap_insert (heap, n);
}

static void update_neighbors (GtsGNode * n, GtsGraphBisection * bg,
                              GtsEHeap * h1, GtsEHeap * h2)
{
  GSList * i;

  i = GTS_SLIST_CONTAINER (n)->items;
  while (i) {
    GtsGNode * n1 = GTS_GNODE_NEIGHBOR (n, i->data);
    if (gts_containee_is_contained (GTS_CONTAINEE (n1), 
                                    GTS_CONTAINER (bg->g))) {
      GtsEHeap * h;
      GtsGraph /* * g1,*/ * g2;
      GHashTable * bg1;

      if (gts_containee_is_contained (GTS_CONTAINEE (n1),
                                      GTS_CONTAINER (bg->g1))) {
        h = h1;
        //g1 = bg->g1;
        g2 = bg->g2;
        bg1 = bg->bg1;
      }
      else {
        h = h2;
        //g1 = bg->g2;
        g2 = bg->g1;
        bg1 = bg->bg2;
      }
      g_hash_table_remove (bg1, n1);
      if (h && GTS_OBJECT (n1)->reserved && GTS_OBJECT (n1)->reserved != n1) {
        gts_eheap_remove (h, GTS_OBJECT (n1)->reserved);
        GTS_OBJECT (n1)->reserved = NULL;
      }
      if (gts_gnode_degree (n1, g2)) {
        g_hash_table_insert (bg1, n1, n1);
        if (h && GTS_OBJECT (n1)->reserved != n1)
          GTS_OBJECT (n1)->reserved = gts_eheap_insert (h, n1);
      }
    }
    i = i->next;
  }  
}

gdouble gts_graph_bisection_bkl_refine (GtsGraphBisection * bg,
                                        guint mmax,
                                        gfloat imbalance)
{
  GtsEHeap * h1, * h2;
  GtsGNode * n;
  guint nm = 0, i;
  GtsGNode ** moves;
  gdouble bestcost = 0., totalcost = 0., best_balance;
  gboolean balanced = FALSE;

  g_return_val_if_fail (bg != NULL, 0.);
  g_return_val_if_fail (mmax > 0, 0.);
  g_return_val_if_fail (imbalance >= 0. && imbalance <= 1., 0.);

  h1 = gts_eheap_new ((GtsKeyFunc) node_move_cost1, bg);
  gts_eheap_freeze (h1);
  g_hash_table_foreach (bg->bg1, (GHFunc) build_bheap, h1);
  gts_eheap_thaw (h1);

  h2 = gts_eheap_new ((GtsKeyFunc) node_move_cost2, bg);
  gts_eheap_freeze (h2);
  g_hash_table_foreach (bg->bg2, (GHFunc) build_bheap, h2);
  gts_eheap_thaw (h2);

  moves = g_malloc (sizeof (GtsGNode *)*mmax);
  imbalance *= gts_graph_weight (bg->g);
  best_balance = fabs (gts_graph_weight (bg->g1) - gts_graph_weight (bg->g2));
  if (best_balance <= imbalance)
    balanced = TRUE;

  do {
    GtsGraph * g1, * g2;
    GHashTable * bg1, * bg2;
    gdouble cost;

    if (gts_graph_weight (bg->g1) > gts_graph_weight (bg->g2)) {
      n = gts_eheap_remove_top (h1, &cost);
      g1 = bg->g1;
      g2 = bg->g2;
      bg1 = bg->bg1;
      bg2 = bg->bg2;
    }
    else {
      n = gts_eheap_remove_top (h2, &cost);
      g1 = bg->g2;
      g2 = bg->g1;
      bg1 = bg->bg2;
      bg2 = bg->bg1;
    }
    if (n) {
      gdouble balance;
        
      GTS_OBJECT (n)->reserved = n;
      gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
      gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
      g_hash_table_remove (bg1, n);
      if (gts_gnode_degree (n, g1))
        g_hash_table_insert (bg2, n, n);

      update_neighbors (n, bg, h1, h2);

      totalcost += cost;
      balance = fabs (gts_graph_weight (g1) - gts_graph_weight (g2));
      
      if (!balanced && balance <= imbalance) {
        bestcost = totalcost;
        best_balance = balance;
        balanced = TRUE;
        nm = 0;
      }
      else if (totalcost < bestcost && 
               (balance < best_balance || balance <= imbalance)) {
        bestcost = totalcost;
        best_balance = balance;
        nm = 0;
      }
      else if (totalcost == bestcost && balance < best_balance) {
        best_balance = balance;
        nm = 0;
      }
      else
        moves[nm++] = n;
    }
  } while (n && nm < mmax);

  gts_container_foreach (GTS_CONTAINER (bg->g), 
                         (GtsFunc) gts_object_reset_reserved, NULL);
  gts_eheap_destroy (h1);
  gts_eheap_destroy (h2);

  /* undo last nm moves */
  for (i = 0; i < nm; i++) {
    GtsGNode * n = moves[i];
    GtsGraph * g1, * g2;
    GHashTable * bg1, * bg2;

    if (gts_containee_is_contained (GTS_CONTAINEE (n),
                                    GTS_CONTAINER (bg->g1))) {
      g1 = bg->g1;
      g2 = bg->g2;
      bg1 = bg->bg1;
      bg2 = bg->bg2;
    }
    else {
      g1 = bg->g2;
      g2 = bg->g1;
      bg1 = bg->bg2;
      bg2 = bg->bg1;
    }
    
    gts_container_add (GTS_CONTAINER (g2), GTS_CONTAINEE (n));
    gts_container_remove (GTS_CONTAINER (g1), GTS_CONTAINEE (n));
    g_hash_table_remove (bg1, n);
    if (gts_gnode_degree (n, g1))
      g_hash_table_insert (bg2, n, n);

    update_neighbors (n, bg, NULL, NULL);
  }
  g_free (moves);

  return bestcost;
}

/* Multilevel partitioning */

static void bisection_children (GtsGNodeSplit * ns, GtsGraphBisection * bg)
{
  GtsGraph * g, * g1;
  GHashTable * bbg;
  GtsGNode * n1 = GTS_GNODE_SPLIT_N1 (ns);
  GtsGNode * n2 = GTS_GNODE_SPLIT_N2 (ns);

  if (gts_containee_is_contained (GTS_CONTAINEE (ns->n),
                                  GTS_CONTAINER (bg->g1))) {
    g = bg->g1;
    g1 = bg->g2;
    bbg = bg->bg1;
  }
  else {
    g = bg->g2;
    g1 = bg->g1;
    bbg = bg->bg2;
  }

  gts_allow_floating_gnodes = TRUE;
  gts_container_remove (GTS_CONTAINER (g), GTS_CONTAINEE (ns->n));
  gts_allow_floating_gnodes = FALSE;
  gts_container_add (GTS_CONTAINER (g), GTS_CONTAINEE (n1));
  gts_container_add (GTS_CONTAINER (g), GTS_CONTAINEE (n2));

  if (g_hash_table_lookup (bbg, ns->n)) {
    g_hash_table_remove (bbg, ns->n);
    if (gts_gnode_degree (n1, g1) > 0)
      g_hash_table_insert (bbg, n1, n1);
    if (gts_gnode_degree (n2, g1) > 0)
      g_hash_table_insert (bbg, n2, n2);
  }
}

GtsGraphBisection * gts_graph_bisection_new (GtsWGraph * wg,
                                             guint ntry,
                                             guint mmax,
                                             guint nmin,
                                             gfloat imbalance)
{
  GtsGraph * g;
  GtsPGraph * pg;
  GtsGraphBisection * bg;
  gdouble cost;

  g_return_val_if_fail (wg != NULL, NULL);

  g = GTS_GRAPH (wg);
  pg = gts_pgraph_new (gts_pgraph_class (), g, 
                       gts_gnode_split_class (),
                       gts_wgnode_class (),
                       gts_wgedge_class (),
                       nmin);

  bg = gts_graph_ggg_bisection (g, ntry);
#ifdef DEBUG
  fprintf (stderr, "before size: %5d weight: %5g cuts: %5d cweight: %5g\n",
           gts_container_size (GTS_CONTAINER (bg->g1)),
           gts_graph_weight (bg->g1),
           gts_graph_edges_cut (bg->g1),
           gts_graph_edges_cut_weight (bg->g1));
  g_assert (gts_graph_bisection_check (bg));
#endif
  while ((cost = gts_graph_bisection_bkl_refine (bg, mmax, imbalance))) {
#ifdef DEBUG
    fprintf (stderr, "  cost: %g\n", cost);
    g_assert (gts_graph_bisection_check (bg));
#endif
  }
#ifdef DEBUG
  fprintf (stderr, "after  size: %5d weight: %5g cuts: %5d cweight: %5g\n",
           gts_container_size (GTS_CONTAINER (bg->g1)),
           gts_graph_weight (bg->g1),
           gts_graph_edges_cut (bg->g1),
           gts_graph_edges_cut_weight (bg->g1));
#endif
  while (gts_pgraph_down (pg, (GtsFunc) bisection_children, bg)) {
#ifdef DEBUG
    fprintf (stderr, "before size: %5d weight: %5g cuts: %5d cweight: %5g\n",
             gts_container_size (GTS_CONTAINER (bg->g1)),
             gts_graph_weight (bg->g1),
             gts_graph_edges_cut (bg->g1),
             gts_graph_edges_cut_weight (bg->g1));         
#endif
    while ((cost = gts_graph_bisection_bkl_refine (bg, mmax, imbalance))) {
#ifdef DEBUG
      fprintf (stderr, "  cost: %g\n", cost);
      g_assert (gts_graph_bisection_check (bg));
#endif
    }
#ifdef DEBUG
    fprintf (stderr, "after  size: %5d weight: %5g cuts: %5d cweight: %5g\n",
             gts_container_size (GTS_CONTAINER (bg->g1)),
             gts_graph_weight (bg->g1),
             gts_graph_edges_cut (bg->g1),
             gts_graph_edges_cut_weight (bg->g1));
#endif
  }
  gts_object_destroy (GTS_OBJECT (pg));

  return bg;
}

void gts_graph_bisection_destroy (GtsGraphBisection * bg,
                                  gboolean destroy_graphs)
{
  g_return_if_fail (bg != NULL);

  g_hash_table_destroy (bg->bg1);
  g_hash_table_destroy (bg->bg2);

  if (destroy_graphs) {
    gts_object_destroy (GTS_OBJECT (bg->g1));
    gts_object_destroy (GTS_OBJECT (bg->g2));
  }

  g_free (bg);
}

static void recursive_bisection (GtsWGraph * wg,
                                 guint np,
                                 guint ntry,
                                 guint mmax,
                                 guint nmin,
                                 gfloat imbalance,
                                 GSList ** list)
{
  if (np == 0)
    *list = g_slist_prepend (*list, wg);
  else {
    GtsGraphBisection * bg = 
      gts_graph_bisection_new (wg, ntry, mmax, nmin, imbalance);
    GtsGraph * g1 = bg->g1;
    GtsGraph * g2 = bg->g2;

    gts_object_destroy (GTS_OBJECT (wg));
    gts_graph_bisection_destroy (bg, FALSE);
    recursive_bisection (GTS_WGRAPH (g1), np - 1, ntry, mmax, nmin, imbalance,
                         list);
    recursive_bisection (GTS_WGRAPH (g2), np - 1, ntry, mmax, nmin, imbalance,
                         list);
  }
}

GSList * gts_graph_recursive_bisection (GtsWGraph * wg,
                                        guint n,
                                        guint ntry,
                                        guint mmax,
                                        guint nmin,
                                        gfloat imbalance)
{
  GtsGraphBisection * bg;
  GtsGraph * g1, * g2;
  GSList * list = NULL;

  g_return_val_if_fail (wg != NULL, NULL);
  g_return_val_if_fail (n > 0, NULL);
  
  bg = gts_graph_bisection_new (wg, ntry, mmax, nmin, imbalance);
  g1 = bg->g1;
  g2 = bg->g2;
  gts_graph_bisection_destroy (bg, FALSE);
  recursive_bisection (GTS_WGRAPH (g1), n - 1, ntry, mmax, nmin, imbalance, 
                       &list);
  recursive_bisection (GTS_WGRAPH (g2), n - 1, ntry, mmax, nmin, imbalance,
                       &list);

  return list;
}