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 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 <algorithm>
#include <common/types.h>
#include <common/globals.h>
#include <sparse/general.h>
#include <math.h>
#include <time.h>
#include <sparse/SparseMatrix.h>
#include <mingle/edge_bundling.h>
#include <mingle/ink.h>
#include <mingle/agglomerative_bundling.h>
#include <mingle/nearest_neighbor_graph.h>
#include <string.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) {
  int n = A->n, i;
 
  assert(SparseMatrix_is_symmetric(A, true));
 
  if (!A) return {};
  assert(A->m == 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 = NULL, *jp = NULL, 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 i, 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 == 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 (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 (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 (i = 0; i < R->m; i++){
      pick = &(ja[ia[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 (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);
}