agglomerative_bundling.cpp

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 v2.0
 * which accompanies this distribution, and is available at
 * https://www.eclipse.org/org/documents/epl-2.0/EPL-2.0.html
 *
 * Contributors: Details at https://graphviz.org
 *************************************************************************/
 
#include "config.h"
 
#include <algorithm>
#include <cassert>
#include <common/globals.h>
#include <common/types.h>
#include <math.h>
#include <mingle/agglomerative_bundling.h>
#include <mingle/edge_bundling.h>
#include <mingle/ink.h>
#include <mingle/nearest_neighbor_graph.h>
#include <sparse/SparseMatrix.h>
#include <sparse/general.h>
#include <string.h>
#include <time.h>
#include <vector>
 
enum { MINGLE_DEBUG = 0 };
 
namespace {
struct Agglomerative_Ink_Bundling {
  Agglomerative_Ink_Bundling(int level_, int n_, SparseMatrix A_,
                             const std::vector<pedge> &edges_)
      : level(level_), n(n_), A(A_), edges(edges_) {}
 
  int level; /* 0, 1, ... */
  int n;
  SparseMatrix A; /* n x n matrix, where n is the number of edges/bundles in
                     this level */
  SparseMatrix R0 = nullptr; /* this is basically R[level - 1].R[level -
                      2]...R[0], which gives the map of bundling i to the
                      original edges: first row of R0 gives the nodes on the
                      finest grid corresponding to the coarsest node 1, etc */
  SparseMatrix R = nullptr;  
  std::vector<double>
      inks; /* amount of ink needed to draw this edge/bundle. Dimension n. */
  double total_ink = -1; /* amount of ink needed to draw this edge/bundle.
                       Dimension n. */
  std::vector<pedge>
      edges; /* the original edge info. This does not vary level to level and
                 is of dimenion n0, where n0 is the number of original edges */
  bool delete_top_level_A = false; /*whether the top level matrix should be
                              deleted on garbage collecting the grid */
};
} // namespace
 
using aib_t = std::vector<Agglomerative_Ink_Bundling>;
 
static aib_t Agglomerative_Ink_Bundling_init(SparseMatrix A,
                                             const std::vector<pedge> &edges,
                                             int level) {
  assert(A != nullptr);
  int n = A->n, i;
 
  assert(SparseMatrix_is_symmetric(A, true));
 
  assert(A->m == (size_t)n);
  Agglomerative_Ink_Bundling grid(level, n, A, edges);
  if (level == 0) {
    double total_ink = 0;
    for (i = 0; i < n; i++) {
      grid.inks.push_back(ink1(edges[i]));
      total_ink += grid.inks[i];
    }
    grid.total_ink = total_ink;
  }
  return aib_t{grid};
}
 
static void Agglomerative_Ink_Bundling_delete(aib_t &grid) {
  for (Agglomerative_Ink_Bundling &a : grid) {
    if (a.A) {
      if (a.level == 0) {
        if (a.delete_top_level_A)
          SparseMatrix_delete(a.A);
      } else {
        SparseMatrix_delete(a.A);
      }
    }
    /* on level 0, R0 = NULL, on level 1, R0 = R */
    if (a.level > 1)
      SparseMatrix_delete(a.R0);
    SparseMatrix_delete(a.R);
  }
}
 
static void Agglomerative_Ink_Bundling_establish(aib_t &grid, int *pick,
                                                 double angle_param,
                                                 double angle) {
  /* pick is a work array of dimension n, with n the total number of original
   * edges */
  SparseMatrix A = grid.front().A;
  int n = grid.front().n, level = grid.front().level, nc = 0;
  int *ia = A->ia, *ja = A->ja;
  int i, j, k, jj, jc, jmax, ni, nj, npicks;
  const std::vector<pedge> &edges = grid.front().edges;
  const std::vector<double> &inks = grid.front().inks;
  double inki, inkj;
  double gain, maxgain, minink, total_gain = 0;
  int *ip = nullptr, *jp = nullptr, ie;
  std::vector<std::vector<int>>
      cedges; /* a table listing the content of bundled edges in the coarsen
grid. cedges[i] contain the list of origonal edges that make up the bundle i in
the next level */
  double ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
  point_t meet1, meet2;
 
  if (Verbose > 1)
    fprintf(stderr, "level ===================== %d, n = %d\n",
            grid.front().level, n);
  cedges.resize(n);
  std::vector<double> cinks(n, 0.0);
 
  if (grid.front().level > 0) {
    ip = grid.front().R0->ia;
    jp = grid.front().R0->ja;
  }
 
  assert(n == A->n);
  std::vector<int> matching(n, UNMATCHED);
 
  for (i = 0; i < n; i++) {
    if (matching[i] != UNMATCHED)
      continue;
 
    /* find the best matching in ink saving */
    maxgain = 0;
    minink = -1;
    jmax = -1;
    for (j = ia[i]; j < ia[i + 1]; j++) {
      jj = ja[j];
      if (jj == i)
        continue;
 
      /* ink saving of merging i and j */
      if ((jc = matching[jj]) == UNMATCHED) {
        /* neither i nor jj are matched */
        inki = inks[i];
        inkj = inks[jj];
        if (ip && jp) {             /* not the first level */
          ni = ip[i + 1] - ip[i];   // number of edges represented by i
          nj = ip[jj + 1] - ip[jj]; // number of edges represented by jj
          memcpy(pick, &jp[ip[i]], sizeof(int) * ni);
          memcpy(pick + ni, &jp[ip[jj]], sizeof(int) * nj);
        } else { /* first level */
          pick[0] = i;
          pick[1] = jj;
          ni = nj = 1;
        }
        if (MINGLE_DEBUG)
          if (Verbose)
            fprintf(stderr, "ink(%d)=%f, ink(%d)=%f", i, inki, jj, inkj);
      } else {
        /* j is already matched. Its content is on cedges[jc] */
        inki = inks[i];
        inkj = cinks[jc];
        if (MINGLE_DEBUG)
          if (Verbose)
            fprintf(stderr, "ink(%d)=%f, ink(%d->%d)=%f", i, inki, jj, jc,
                    inkj);
        if (ip) {
          ni = ip[i + 1] - ip[i]; // number of edges represented by i
          memcpy(pick, &jp[ip[i]], sizeof(int) * ni);
        } else {
          ni = 1;
          pick[0] = i;
        }
        nj = cedges[jc].size();
        npicks = ni;
        for (k = 0; k < nj; k++) {
          pick[npicks++] = cedges[jc][k];
        }
      }
 
      npicks = ni + nj;
      ink1 =
          ink(edges, npicks, pick, &ink0, &meet1, &meet2, angle_param, angle);
      if (MINGLE_DEBUG) {
        if (Verbose) {
          fprintf(stderr, ", if merging {");
          for (k = 0; k < npicks; k++)
            fprintf(stderr, "%d,", pick[k]);
          fprintf(stderr, "}, ");
          fprintf(stderr, " ink0=%f, ink1=%f", inki + inkj, ink1);
        }
      }
 
      gain = inki + inkj - ink1;
      if (MINGLE_DEBUG)
        if (Verbose)
          fprintf(stderr, " gain=%f", gain);
      if (gain > maxgain) {
        maxgain = gain;
        minink = ink1;
        jmax = jj;
        if (MINGLE_DEBUG)
          if (Verbose)
            fprintf(stderr, "maxgain=%f", maxgain);
      }
      if (MINGLE_DEBUG)
        if (Verbose)
          fprintf(stderr, "\n");
    }
 
    /* now merge i and jmax */
    if (maxgain > 0) {
      /* a good bundling of i and another edge jmax is found */
      total_gain += maxgain;
      jc = matching[jmax];
      if (jc ==
          UNMATCHED) { /* i and j both unmatched. Add j in the table first */
        if (MINGLE_DEBUG)
          if (Verbose)
            printf("maxgain=%f, merge %d with best edge: %d to form coarsen "
                   "edge %d. Ink=%f\n",
                   maxgain, i, jmax, nc, minink);
        matching[i] = matching[jmax] = nc;
        if (ip) {
          for (k = ip[jmax]; k < ip[jmax + 1]; k++) {
            ie = jp[k];
            cedges[nc].push_back(ie);
          }
        } else {
          cedges[nc].push_back(jmax);
        }
        jc = nc;
        nc++;
      } else { /*j is already matched */
        if (MINGLE_DEBUG)
          if (Verbose)
            printf("maxgain=%f, merge %d with existing cluster %d\n", maxgain,
                   i, jc);
        matching[i] = jc;
        grand_total_ink -= cinks[jc]; /* ink of cluster jc was already added,
                                         and will be added again as part of a
                                         larger cluster with i, so dicount */
      }
    } else { /*i can not match/bundle successfully */
      if (MINGLE_DEBUG)
        if (Verbose)
          fprintf(stderr, "no gain in bundling node %d\n", i);
      assert(maxgain <= 0);
      matching[i] = nc;
      jc = nc;
      minink = inks[i];
      nc++;
    }
 
    /* add i to the appropriate table */
    if (ip) {
      for (k = ip[i]; k < ip[i + 1]; k++) {
        ie = jp[k];
        cedges[jc].push_back(ie);
      }
    } else {
      cedges[jc].push_back(i);
    }
    cinks[jc] = minink;
    grand_total_ink += minink;
    grand_total_gain += maxgain;
 
    if (MINGLE_DEBUG) {
      if (Verbose) {
        fprintf(stderr, " coarse edge[%d]={", jc);
        for (const int &cedge : cedges[jc]) {
          fprintf(stderr, "%d,", cedge);
        }
        fprintf(stderr, "}, grand_total_gain=%f\n", grand_total_gain);
      }
    }
  }
 
  if (nc >= 1 && total_gain > 0) {
    /* now set up restriction and prolongation operator */
    SparseMatrix P, R, R1, R0, B, cA;
    double one = 1.;
 
    R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
    for (i = 0; i < n; i++) {
      jj = matching[i];
      SparseMatrix_coordinate_form_add_entry(R1, jj, i, &one);
    }
    R = SparseMatrix_from_coordinate_format(R1);
    SparseMatrix_delete(R1);
    P = SparseMatrix_transpose(R);
    B = SparseMatrix_multiply(R, A);
    if (!B)
      return;
    cA = SparseMatrix_multiply(B, P);
    if (!cA)
      return;
    SparseMatrix_delete(B);
    SparseMatrix_delete(P);
    grid.front().R = R;
 
    level++;
    aib_t cgrid = Agglomerative_Ink_Bundling_init(cA, edges, level);
 
    if (grid.front().R0) {
      R0 = SparseMatrix_multiply(R, grid.front().R0);
    } else {
      assert(grid.front().level == 0);
      R0 = R;
    }
    cgrid.front().R0 = R0;
    cgrid.front().inks = cinks;
    cgrid.front().total_ink = grand_total_ink;
 
    if (Verbose > 1)
      fprintf(stderr,
              "level %d->%d, edges %d -> %d, ink %f->%f , gain = %f, or %f\n",
              grid.front().level, cgrid.front().level, grid.front().n,
              cgrid.front().n, grid.front().total_ink, grand_total_ink,
              grid.front().total_ink - grand_total_ink, grand_total_gain);
    assert(fabs(grid.front().total_ink - cgrid.front().total_ink -
                grand_total_gain) <= 0.0001 * grid.front().total_ink);
 
    Agglomerative_Ink_Bundling_establish(cgrid, pick, angle_param, angle);
    grid.insert(grid.end(), cgrid.begin(), cgrid.end());
 
  } else {
    if (Verbose > 1)
      fprintf(stderr,
              "no more improvement, orig ink = %f, gain = %f, stop and final "
              "bundling found\n",
              grand_total_ink, grand_total_gain);
    /* no more improvement, stop and final bundling found */
  }
}
 
static aib_t Agglomerative_Ink_Bundling_new(SparseMatrix A0,
                                            const std::vector<pedge> &edges,
                                            double angle_param, double angle) {
  /* give a link of edges and their nearest neighbor graph, return a multilevel
   * of edge bundling based on ink saving */
  SparseMatrix A = A0;
 
  if (!SparseMatrix_is_symmetric(A, false) || A->type != MATRIX_TYPE_REAL) {
    A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
  }
  aib_t grid = Agglomerative_Ink_Bundling_init(A, edges, 0);
 
  std::vector<int> pick(A0->m);
 
  Agglomerative_Ink_Bundling_establish(grid, pick.data(), angle_param, angle);
 
  if (A != A0)
    grid.front().delete_top_level_A = true; // be sure to clean up later
 
  return grid;
}
 
static void agglomerative_ink_bundling_internal(
    int dim, SparseMatrix A, std::vector<pedge> &edges, int nneighbors,
    int *recurse_level, int MAX_RECURSE_LEVEL, double angle_param, double angle,
    double *current_ink, double *ink00) {
 
  int j, jj, k;
  int *ia, *ja;
  int *pick;
  SparseMatrix R;
  double ink0, ink1;
  point_t meet1, meet2;
  double wgt_all;
  [[maybe_unused]] const double TOL = 0.0001;
  clock_t start;
 
  (*recurse_level)++;
  if (Verbose > 1)
    fprintf(stderr,
            "agglomerative_ink_bundling_internal, recurse level ------- %d\n",
            *recurse_level);
 
  assert(A->m == (size_t)A->n);
 
  start = clock();
  aib_t grid = Agglomerative_Ink_Bundling_new(A, edges, angle_param, angle);
  if (Verbose > 1)
    fprintf(stderr, "CPU for agglomerative bundling %f\n",
            (double)(clock() - start) / CLOCKS_PER_SEC);
  ink0 = grid.front().total_ink;
 
  /* find coarsest */
  ink1 = grid.back().total_ink;
 
  if (*current_ink < 0) {
    *current_ink = *ink00 = ink0;
    if (Verbose > 1)
      fprintf(stderr, "initial total ink = %f\n", *current_ink);
  }
  if (ink1 < ink0) {
    *current_ink -= ink0 - ink1;
  }
 
  if (Verbose > 1)
    fprintf(stderr,
            "ink: %f->%f, edges: %d->%d, current ink = %f, percentage gain "
            "over original = %f\n",
            ink0, ink1, grid.front().n, grid.back().n, *current_ink,
            (ink0 - ink1) / *ink00);
 
  /* if no meaningful improvement (0.0001%), out, else rebundle the middle
   * section */
  if ((ink0 - ink1) / *ink00 < 0.000001 || *recurse_level > MAX_RECURSE_LEVEL) {
    /* project bundles up */
    R = grid.back().R0;
    if (R) {
      ia = R->ia;
      ja = R->ja;
      for (size_t i = 0; i < R->m; i++) {
        pick = &ja[ia[i]];
 
        if (MINGLE_DEBUG)
          if (Verbose)
            fprintf(stderr, "calling ink2...\n");
        ink1 = ink(edges, ia[i + 1] - ia[i], pick, &ink0, &meet1, &meet2,
                   angle_param, angle);
        if (MINGLE_DEBUG)
          if (Verbose)
            fprintf(stderr, "finish calling ink2...\n");
        assert(fabs(ink1 - grid.back().inks[i]) <= std::max(TOL, TOL * ink1) &&
               ink1 - ink0 <= TOL);
        assert(ink1 < 1000 * ink0); /* assert that points were found */
        wgt_all = 0.;
        if (ia[i + 1] - ia[i] > 1) {
          for (j = ia[i]; j < ia[i + 1]; j++) {
            /* make this edge 4 points, insert two meeting points at 1 and 2,
             * make 3 the last point */
            jj = ja[j];
            pedge_double(edges[jj]); /* has to call pedge_double twice: from 2
                                        points to 3 points to 5 points. The last
                                        point not used, may be improved later */
            pedge_double(edges[jj]);
            pedge &e = edges[jj];
 
            e.x[1 * dim] = meet1.x;
            e.x[1 * dim + 1] = meet1.y;
            e.x[2 * dim] = meet2.x;
            e.x[2 * dim + 1] = meet2.y;
            e.x[3 * dim] = e.x[4 * dim];
            e.x[3 * dim + 1] = e.x[4 * dim + 1];
            e.npoints = 4;
            e.wgts = std::vector<double>(4, e.wgt);
            wgt_all += e.wgt;
          }
          for (j = ia[i]; j < ia[i + 1]; j++) {
            pedge &e = edges[ja[j]];
            e.wgts[1] = wgt_all;
          }
        }
      }
    }
  } else {
    int ne, npp, l;
    SparseMatrix A_mid;
    double wgt;
 
    /* make new edges using meet1 and meet2.
 
        call Agglomerative_Ink_Bundling_new
 
        inherit new edges to old edges
     */
 
    R = grid.back().R0;
    assert(
        R &&
        grid.size() >
            1); /* if ink improved, we should have gone at leat 1 level down! */
    ia = R->ia;
    ja = R->ja;
    ne = R->m;
    std::vector<pedge> mid_edges(ne);
    std::vector<double> xx(4 * ne);
    for (size_t i = 0; i < R->m; i++) {
      pick = &ja[ia[i]];
      wgt = 0.;
      for (j = ia[i]; j < ia[i + 1]; j++)
        wgt += edges[j].wgt;
      if (MINGLE_DEBUG)
        if (Verbose)
          fprintf(stderr, "calling ink3...\n");
      ink1 = ink(edges, ia[i + 1] - ia[i], pick, &ink0, &meet1, &meet2,
                 angle_param, angle);
      if (MINGLE_DEBUG)
        if (Verbose)
          fprintf(stderr, "done calling ink3...\n");
      assert(fabs(ink1 - grid.back().inks[i]) <= std::max(TOL, TOL * ink1) &&
             ink1 - ink0 <= TOL);
      assert(ink1 < 1000 * ink0); /* assert that points were found */
      xx[i * 4 + 0] = meet1.x;
      xx[i * 4 + 1] = meet1.y;
      xx[i * 4 + 2] = meet2.x;
      xx[i * 4 + 3] = meet2.y;
      mid_edges[i] = pedge_wgt_new(2, dim, &xx.data()[i * 4], wgt);
    }
 
    A_mid = nearest_neighbor_graph(ne, std::min(nneighbors, ne), xx);
 
    agglomerative_ink_bundling_internal(dim, A_mid, mid_edges, nneighbors,
                                        recurse_level, MAX_RECURSE_LEVEL,
                                        angle_param, angle, current_ink, ink00);
    SparseMatrix_delete(A_mid);
 
    /* patching edges with the new mid-section */
    for (size_t i = 0; i < R->m; i++) {
      // middle section of edges that will be bundled again
      const pedge &midedge = mid_edges[i];
      npp = midedge.npoints + 2;
      for (j = ia[i]; j < ia[i + 1]; j++) {
        jj = ja[j];
        pedge_wgts_realloc(edges[jj], npp);
        pedge &e = edges[jj];
 
        assert(e.npoints == 2);
        for (l = 0; l < dim; l++) { /* move the second point to the last */
          e.x[(npp - 1) * dim + l] = e.x[1 * dim + l];
        }
 
        for (k = 0; k < midedge.npoints; k++) {
          for (l = 0; l < dim; l++) {
            e.x[(k + 1) * dim + l] = midedge.x[k * dim + l];
          }
          if (k < midedge.npoints - 1) {
            if (!midedge.wgts.empty()) {
              e.wgts[k + 1] = midedge.wgts[k];
            } else {
              e.wgts[k + 1] = midedge.wgt;
            }
          }
        }
        e.wgts[npp - 2] =
            e.wgts[0]; // the last interval take from the 1st interval
 
        e.npoints = npp;
      }
    }
 
    for (int i = 0; i < ne; i++)
      pedge_delete(mid_edges[i]);
  }
 
  Agglomerative_Ink_Bundling_delete(grid);
}
 
void agglomerative_ink_bundling(int dim, SparseMatrix A,
                                std::vector<pedge> &edges, int nneighbor,
                                int MAX_RECURSE_LEVEL, double angle_param,
                                double angle) {
  int recurse_level = 0;
  double current_ink = -1, ink0;
 
  ink_count = 0;
  agglomerative_ink_bundling_internal(dim, A, edges, nneighbor, &recurse_level,
                                      MAX_RECURSE_LEVEL, angle_param, angle,
                                      &current_ink, &ink0);
 
  if (Verbose > 1)
    fprintf(stderr,
            "initial total ink = %f, final total ink = %f, inksaving = %f "
            "percent, total ink_calc = %f, avg ink_calc per edge = %f\n",
            ink0, current_ink, (ink0 - current_ink) / ink0, ink_count,
            ink_count / (double)A->m);
}