closest.c

Go to the documentation of this file.

/*************************************************************************
 * Copyright (c) 2011 AT&T Intellectual Property 
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v1.0
 * which accompanies this distribution, and is available at
 * https://www.eclipse.org/legal/epl-v10.html
 *
 * Contributors: Details at https://graphviz.org
 *************************************************************************/
 
#include <assert.h>
#include <limits.h>
#include <neatogen/kkutils.h>
#include <neatogen/closest.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdlib.h>
#include <util/alloc.h>
#include <util/gv_math.h>
#include <util/list.h>
#include <util/sort.h>
 
/*****************************************
** This module contains functions that  **
** given a 1-D layout construct a graph **
** where close nodes in this layout are **
** adjacent                             **
*****************************************/
 
typedef struct {
    /* this structure represents two nodes in the 1-D layout */
    size_t left;  
    size_t right; 
    double dist;                /* distance between the nodes in the layout */
} Pair;
 
#define LT(p,q) ((p).dist < (q).dist)
#define EQ(p,q) ((p).dist == (q).dist)
 
typedef LIST(Pair *) pairs_t;
 
static void push(pairs_t *s, Pair x) {
 
  // copy the item to save it on the stack
  Pair *copy = gv_alloc(sizeof(x));
  *copy = x;
 
  LIST_PUSH_BACK(s, copy);
}
 
static bool pop(pairs_t *s, Pair *x) {
 
  // is the stack empty?
  if (LIST_IS_EMPTY(s)) {
    return false;
  }
 
  // remove the top and pass it back to the caller
  Pair *top = LIST_POP_BACK(s);
  *x = *top;
 
  // discard the previously saved copy
  free(top);
 
  return true;
}
 
#define sub(h,i) (*LIST_GET((h), (i)))
 
/* An auxulliary data structure (a heap) for 
 * finding the closest pair in the layout
 */
typedef struct {
    Pair *data;
    size_t heapSize;
    size_t maxSize;
} PairHeap;
 
#define left(i) (2*(i))
#define right(i) (2*(i)+1)
#define parent(i) ((i)/2)
#define insideHeap(h,i) ((i)<h->heapSize)
#define greaterPriority(h,i,j) \
  (LT(h->data[i],h->data[j]) || ((EQ(h->data[i],h->data[j])) && (rand()%2)))
 
static void heapify(PairHeap *h, size_t i) {
    size_t largest;
    while (1) {
        size_t l = left(i);
        size_t r = right(i);
        if (insideHeap(h, l) && greaterPriority(h, l, i))
            largest = l;
        else
            largest = i;
        if (insideHeap(h, r) && greaterPriority(h, r, largest))
            largest = r;
        if (largest == i)
            break;
 
        SWAP(&h->data[largest], &h->data[i]);
        i = largest;
    }
}
 
static void freeHeap(PairHeap * h)
{
    free(h->data);
}
 
static void initHeap(PairHeap *h, double *place, size_t *ordering, size_t n) {
    Pair edge;
 
    h->heapSize = n == 0 ? 0 : (n - 1);
    h->maxSize = h->heapSize;
    h->data = gv_calloc(h->maxSize, sizeof(Pair));
 
    for (size_t i = 0; n != 0 && i < n - 1; i++) {
        edge.left = ordering[i];
        edge.right = ordering[i + 1];
        edge.dist = place[ordering[i + 1]] - place[ordering[i]];
        h->data[i] = edge;
    }
    for (size_t j = (n - 1) / 2; n != 0 && j != SIZE_MAX; j--) {
        heapify(h, j);
    }
}
 
static bool extractMax(PairHeap * h, Pair * max)
{
    if (h->heapSize == 0)
        return false;
 
    *max = h->data[0];
    h->data[0] = h->data[h->heapSize - 1];
    h->heapSize--;
    heapify(h, 0);
    return true;
}
 
static void insert(PairHeap * h, Pair edge)
{
    size_t i = h->heapSize;
    if (h->heapSize == h->maxSize) {
        size_t newSize = h->maxSize * 2;
        h->data = gv_recalloc(h->data, h->maxSize, newSize, sizeof(Pair));
        h->maxSize = newSize;
    }
    h->heapSize++;
    h->data[i] = edge;
    while (i > 0 && greaterPriority(h, i, parent(i))) {
        SWAP(&h->data[i], &h->data[parent(i)]);
        i = parent(i);
    }
}
 
static int cmp(const void *a, const void *b, void *context) {
  const size_t *x = a;
  const size_t *y = b;
  const double *place = context;
 
  if (place[*x] < place[*y]) {
    return -1;
  }
  if (place[*x] > place[*y]) {
    return 1;
  }
  return 0;
}
 
static void find_closest_pairs(double *place, size_t n, int num_pairs,
                               pairs_t* pairs_stack) {
    /* Fill the stack 'pairs_stack' with 'num_pairs' closest pairs int the 1-D layout 'place' */
    PairHeap heap;
    size_t *left = gv_calloc(n, sizeof(size_t));
    size_t *right = gv_calloc(n, sizeof(size_t));
    Pair pair = {0}, new_pair;
 
    /* Order the nodes according to their place */
    size_t *ordering = gv_calloc(n, sizeof(size_t));
    size_t *inv_ordering = gv_calloc(n, sizeof(size_t));
 
    for (size_t i = 0; i < n; i++) {
        ordering[i] = i;
    }
    gv_sort(ordering, n, sizeof(ordering[0]), cmp, place);
    for (size_t i = 0; i < n; i++) {
        inv_ordering[ordering[i]] = i;
    }
 
    /* Initialize heap with all consecutive pairs */
    initHeap(&heap, place, ordering, n);
 
    /* store the leftmost and rightmost neighbors of each node that were entered into heap */
    for (size_t i = 1; i < n; i++) {
        left[ordering[i]] = ordering[i - 1];
    }
    for (size_t i = 0; n != 0 && i < n - 1; i++) {
        right[ordering[i]] = ordering[i + 1];
    }
 
    /* extract the 'num_pairs' closest pairs */
    for (int i = 0; i < num_pairs; i++) {
 
        if (!extractMax(&heap, &pair)) {
            break;              /* not enough pairs */
        }
        push(pairs_stack, pair);
        /* insert to heap "descendant" pairs */
        size_t left_index = inv_ordering[pair.left];
        size_t right_index = inv_ordering[pair.right];
        if (left_index > 0) {
            size_t neighbor = ordering[left_index - 1];
            if (inv_ordering[right[neighbor]] < right_index) {
                /* we have a new pair */
                new_pair.left = neighbor;
                new_pair.right = pair.right;
                new_pair.dist = place[pair.right] - place[neighbor];
                insert(&heap, new_pair);
                right[neighbor] = pair.right;
                left[pair.right] = neighbor;
            }
        }
        if (right_index < n - 1) {
            size_t neighbor = ordering[right_index + 1];
            if (inv_ordering[left[neighbor]] > left_index) {
                /* we have a new pair */
                new_pair.left = pair.left;
                new_pair.right = neighbor;
                new_pair.dist = place[neighbor] - place[pair.left];
                insert(&heap, new_pair);
                left[neighbor] = pair.left;
                right[pair.left] = neighbor;
            }
        }
    }
    free(left);
    free(right);
    free(ordering);
    free(inv_ordering);
    freeHeap(&heap);
}
 
static void add_edge(vtx_data * graph, int u, int v)
{
    for (size_t i = 0; i < graph[u].nedges; i++) {
        if (graph[u].edges[i] == v) {
            /* edge already exist */
            return;
        }
    }
    /* add the edge */
    graph[u].edges[graph[u].nedges++] = v;
    graph[v].edges[graph[v].nedges++] = u;
    if (graph[0].ewgts != NULL) {
        graph[u].ewgts[0]--;
        graph[v].ewgts[0]--;
    }
}
 
static void construct_graph(size_t n, pairs_t *edges_stack,
                            vtx_data **New_graph) {
    /* construct an unweighted graph using the edges 'edges_stack' */
    vtx_data *new_graph;
 
    /* first compute new degrees and nedges; */
    int *degrees = gv_calloc(n, sizeof(int));
    size_t top = LIST_SIZE(edges_stack);
    size_t new_nedges = 2 * top + n;
    Pair pair;
    int *edges = gv_calloc(new_nedges, sizeof(int));
    float *weights = gv_calloc(new_nedges, sizeof(float));
 
    for (size_t i = 0; i < n; i++) {
        degrees[i] = 1;         /* save place for the self loop */
    }
    for (size_t i = 0; i < top; i++) {
        pair = sub(edges_stack, i);
        degrees[pair.left]++;
        degrees[pair.right]++;
    }
 
    /* copy graph into new_graph: */
    for (size_t i = 0; i < new_nedges; i++) {
        weights[i] = 1.0;
    }
 
    *New_graph = new_graph = gv_calloc(n, sizeof(vtx_data));
    for (size_t i = 0; i < n; i++) {
        new_graph[i].nedges = 1;
        new_graph[i].ewgts = weights;
        new_graph[i].edges = edges;
        assert(i <= INT_MAX);
        *edges = (int)i; // self loop for Lap
        *weights = 0;           /* self loop weight for Lap */
        weights += degrees[i];
        edges += degrees[i];    /* reserve space for possible more edges */
    }
 
    free(degrees);
 
    /* add all edges from stack */
    while (pop(edges_stack, &pair)) {
        assert(pair.left <= INT_MAX);
        assert(pair.right <= INT_MAX);
        add_edge(new_graph, (int)pair.left, (int)pair.right);
    }
}
 
void
closest_pairs2graph(double *place, int n, int num_pairs, vtx_data ** graph)
{
    /* build a graph with with edges between the 'num_pairs' closest pairs in the 1-D space: 'place' */
    pairs_t pairs_stack = {0};
    assert(n >= 0);
    find_closest_pairs(place, (size_t)n, num_pairs, &pairs_stack);
    construct_graph((size_t)n, &pairs_stack, graph);
    LIST_FREE(&pairs_stack);
}